Hillslope response to rainfall remains one of the central problems of catchment hydrology. Flow processes in a one-dimensional sloping aquifer can be described by Boussinesq's hydraulic groundwater theory. Most hillslopes, however, have complex three-dimensional shapes that are characterized by their plan shape, profile curvature of surface and bedrock, and the soil depth. Field studies and numerical simulation have shown that these attributes are the most significant topographic controls on subsurface flow and saturation along hillslopes. In this paper the Boussinesq equation is reformulated in terms of soil water storage rather than water table height. The continuity and Darcy equations formulated in terms of storage along the hillslope lead to the hillslope-storage Boussinesq (HSB) equation for subsurface flow. Solutions of the HSB equation account explicitly for plan shape of the hillslope by introducing the hillslope width function and for profile curvature through the bedrock slope angle and the hillslope soil depth function. We investigate the behavior of the HSB model for different hillslope types (uniform, convergent, and divergent) and different slope angles under free drainage conditions after partial initial saturation (drainage scenario) and under constant rainfall recharge conditions (recharge scenario). The HSB equation is solved by means of numerical integration of the partial differential equation. We find that convergent hillslopes drain much more slowly compared to divergent hillslopes. The accumulation of moisture storage near the outlet of convergent hillslopes results in bell-shaped hydrographs. In contrast, the fast draining divergent hillslopes produce highly peaked hydrographs. In order to investigate the relative importance of the different terms in the HSB equation, several simplified nonlinear and linearized versions are derived, for instance, by recognizing that the width function of a hillslope generally shows smooth transition along the flow direction or by introducing a fitting parameter to account for average storage along the hillslope. The dynamic response of these reduced versions of the HSB equation under free drainage conditions depend strongly on hillslope shape and bedrock slope angle. For flat slopes (of the order of 5%), only the simplified nonlinear HSB equation is able to capture the dynamics of subsurface flow along complex hillslopes. In contrast, for steep slopes (of the order of 30%), we see that all the reduced versions show very similar results compared to the full version. It can be concluded that the complex derivative terms of width with respect to flow distance play a less dominant role with increasing slope angle. Comparison with the hillslope-storage kinematic wave model of Troch et al. <2002> shows that the diffusivedrainage terms of the HSB model become less important for the fast draining divergent hillslopes. These results have important implications for the use of simplified versions of the HSB equation in landscapes and for the development of appropriate analytical solutions for subsurface flow along complex hillslopes.