[1] 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 diffusive drainage terms of the HSB model become less important for the fast draining divergent hillslopes. These results VOL. 39, NO. 11, 1316, doi:10.1029/2002WR001728, 2003 have important implications for the use of simplified versions of the H...

[1] A distributed physically based model incorporating novel approaches for the representation of surface-subsurface processes and interactions is presented. A path-based description of surface flow across the drainage basin is used, with several options for identifying flow directions, for separating channel cells from hillslope cells, and for representing stream channel hydraulic geometry. Lakes and other topographic depressions are identified and specially treated as part of the preprocessing procedures applied to the digital elevation data for the catchment. Threshold-based boundary condition switching is used to partition potential (atmospheric) fluxes into actual fluxes across the land surface and changes in surface storage, thus resolving the exchange fluxes, or coupling, between the surface and subsurface modules. Nested time stepping allows smaller steps to be taken for typically faster and explicitly solved surface runoff routing, while a mesh coarsening option allows larger grid elements to be used for typically slower and more compute-intensive subsurface flow. Sequential data assimilation schemes allow the model predictions to be updated with spatiotemporal observation data of surface and subsurface variables. These approaches are discussed in detail, and the physical and numerical behavior of the model is illustrated over catchment scales ranging from 0.0027 to 356 km 2 , addressing different hydrological processes and highlighting the importance of describing coupled surfacesubsurface flow.Citation: Camporese, M., C. Paniconi, M. Putti, and S. Orlandini (2010), Surface-subsurface flow modeling with path-based runoff routing, boundary condition-based coupling, and assimilation of multisource observation data, Water Resour.

There are a growing number of large-scale, complex hydrologic models that are capable of simulating integrated surface and subsurface flow. Many are coupled to land-surface energy balance models, biogeochemical and ecological process models, and atmospheric models. Although they are being increasingly applied for hydrologic prediction and environmental understanding, very little formal verification and/or benchmarking of these models has been performed. Here we present the results of an intercomparison study of seven coupled surface-subsurface models based on a series of benchmark problems. All the models simultaneously solve adapted forms of the Richards and shallow water equations, based on fully 3-D or mixed (1-D vadose zone and 2-D groundwater) formulations for subsurface flow and 1-D (rill flow) or 2-D (sheet flow) conceptualizations for surface routing. A range of approaches is used for the solution of the coupled equations, including global implicit, sequential iterative, and asynchronous linking, and various strategies are used to enforce flux and pressure continuity at the surface-subsurface interface. The simulation results show good agreement for the simpler test cases, while the more complicated test cases bring out some of the differences in physical process representations and numerical solution approaches between the models. Benchmarks with more traditional runoff generating mechanisms, such as excess infiltration and saturation, demonstrate more agreement between models, while benchmarks with heterogeneity and complex water table dynamics highlight differences in model formulation. In general, all the models demonstrate the same qualitative behavior, thus building confidence in their use for hydrologic applications.

Picard iteration is a widely used procedure for solving the nonlinear equation governing flow in variably saturated porous media. The method is simple to code and computationally cheap, but has been known to fail or converge slowly. The Newton method is more complex and expensive (on a per-iteration basis) than Picard, and as such has not received very much attention. Its robustness and higher rate of convergence, however, make it an attractive alternative to the Picard method, particularly for strongly nonlinear problems. In this paper the Picard and Newton schemes are implemented and compared in one-, two-, and three-dimensional finite element simulations involving both steady state and transient flow. The eight test cases presented highlight different aspects of the performance of the two iterative methods and the different factors that can affect their convergence and efficiency, including problem size, spatial and temporal discretization, initial solution estimates, convergence error norm, mass lureping, time weighting, conductivity and moisture content characteristics, boundary conditions, seepage faces, and the extent of fully satarated zones in the soft. Previous strategies for enhancing the performance of the heard and Newton schemes are revisited, and new ones are suggested. The strategies include chord slope approximations for the derivatives of the characteristic equations, relaxing convergence requirements along seepage faces, dynamic time step control, nonlinear relaxation, and a mixed Picard-Newton approach. The tests show that the Picard or relaxed Picard schemes are often adequate for solving Richards' equation, but that in cases where these fail to converge or converge slowly, the Newton method should be used. The mixed Picard-Newton approach can effectively overcome the Newton scheme's sensitivity to initial solution estimates, while comparatively poor performance is reported for the various chord slope approximations. Finally, given the reliability and efficiency of current conjugate gradient-like methods for solving linear nonsymmetric systems, the only real drawback of using Newton rather than Picard iteration is the algebraic complexity and computational cost of assembling the derivative terms of the JacobJan matrix, and it is suggested that both methods can be effectively implemented and used in numerical models of Richards' equation. great popularity because it is the most intuitive linearization of Richards' equation, is computationally inexpensive on a per-iteration basis, and preserves symmetry of the discrete system of equations. However, the method may diverge under certain conditions, as has been observed empirically [e.g., Huyakorn et al., 1984; Celia et al., 1990] and verified theoretically [Aldama and Paniconi, 1992]. The Newton scheme, also known as Newton-Raphson iteration, yields nonsymmetric system matrices and is more complex and expensive than Picard linearization, though it achieves a higher rate of convergence and can be more robust than Picard for certain types of proble...

Integrated, process‐based numerical models in hydrology are rapidly evolving, spurred by novel theories in mathematical physics, advances in computational methods, insights from laboratory and field experiments, and the need to better understand and predict the potential impacts of population, land use, and climate change on our water resources. At the catchment scale, these simulation models are commonly based on conservation principles for surface and subsurface water flow and solute transport (e.g., the Richards, shallow water, and advection‐dispersion equations), and they require robust numerical techniques for their resolution. Traditional (and still open) challenges in developing reliable and efficient models are associated with heterogeneity and variability in parameters and state variables; nonlinearities and scale effects in process dynamics; and complex or poorly known boundary conditions and initial system states. As catchment modeling enters a highly interdisciplinary era, new challenges arise from the need to maintain physical and numerical consistency in the description of multiple processes that interact over a range of scales and across different compartments of an overall system. This paper first gives an historical overview (past 50 years) of some of the key developments in physically based hydrological modeling, emphasizing how the interplay between theory, experiments, and modeling has contributed to advancing the state of the art. The second part of the paper examines some outstanding problems in integrated catchment modeling from the perspective of recent developments in mathematical and computational science.