Mario Putti

Surface-subsurface flow modeling with path-based runoff routing, boundary condition-based coupling, and assimilation of multisource observation data

Received 21 October 2008; revised 2 September 2009; accepted 16 September 2009; published 13 February 2010. [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.

Surface‐subsurface model intercomparison: A first set of benchmark results to diagnose integrated hydrology and feedbacks

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.

MIXED FINITE ELEMENTS AND NEWTON-TYPE LINEARIZATIONS FOR THE SOLUTION OF RICHARDS' EQUATION

We present the development of a two-dimensional Mixed-Hybrid Finite Element (MHFE) model for the solution of the non-linear equation of variably saturated flow in groundwater on unstructured triangular meshes. By this approach the Darcy velocity is approximated using lowest-order Raviart–Thomas (RT0) elements and is ‘exactly’ mass conserving. Hybridization is used to overcome the ill-conditioning of the mixed system. The scheme is globally first-order in space. Nevertheless, numerical results employing non-uniform meshes show second-order accuracy of the pressure head and normal fluxes on specific grid points. The non-linear systems of algebraic equations resulting from the MHFE discretization are solved using Picard or Newton iterations. Realistic sample tests show that the MHFE-Newton approach achieves fast convergence in many situations, in particular, when a good initial guess is provided by either the Picard scheme or relaxation techniques. Copyright

Newtonian nudging for a Richards equation-based distributed hydrological model

Abstract The objective of data assimilation is to provide physically consistent estimates of spatially distributed environmental variables. In this study a relatively simple data assimilation method has been implemented in a relatively complex hydrological model. The data assimilation technique is Newtonian relaxation or nudging, in which model variables are driven towards observations by a forcing term added to the model equations. The forcing term is proportional to the difference between simulation and observation (relaxation component) and contains four-dimensional weighting functions that can incorporate prior knowledge about the spatial and temporal variability and characteristic scales of the state variable(s) being assimilated. The numerical model couples a three-dimensional finite element Richards equation solver for variably saturated porous media and a finite difference diffusion wave approximation based on digital elevation data for surface water dynamics. We describe the implementation of the data assimilation algorithm for the coupled model and report on the numerical and hydrological performance of the resulting assimilation scheme. Nudging is shown to be successful in improving the hydrological simulation results, and it introduces little computational cost, in terms of CPU and other numerical aspects of the model’s behavior, in some cases even improving numerical performance compared to model runs without nudging. We also examine the sensitivity of the model to nudging term parameters including the spatio-temporal influence coefficients in the weighting functions. Overall the nudging algorithm is quite flexible, for instance in dealing with concurrent observation datasets, gridded or scattered data, and different state variables, and the implementation presented here can be readily extended to any of these features not already incorporated. Moreover the nudging code and tests can serve as a basis for implementation of more sophisticated data assimilation techniques in a Richards equation-based hydrological model.

Picard and Newton linearization for the coupled model for saltwater intrusion in aquifers

Abstract Difficulties in the numerical solution of the partial differential equations governing seawater intrusion in aquifers arise from the coupling between the flow and transport equations and from the nonlinear aspects of this coupling. Several linearization approaches are discussed for the solution of the nonlinear system which results from a finite element discretization of the coupled equations. It is first shown that the most commonly used solution method can be viewed as a Picard linearization applied to the transport equation, with the coupling resolved by iteration over the two governing equations. The full Newton scheme for solving the coupled problem produces a Jacobian of size 2 N × 2 N , where N is the number of nodes in the discretization of both the flow and transport equations. To reduce the size and complexity of the full Newton scheme, a partial Newton method is proposed, which, like the Picard approach, produces matrix systems of size N × N . This scheme applies Newton linearization to the transport equation, and conventional iteration to resolve the coupling. Results from two- and three-dimensional test simulations show that the partial Newton scheme gives improved convergence and robustness compared to Picard linearization, especially for highly advective problems or large density ratios. Both approaches involve the solution of a symmetric (flow) and a nonsymmetric (transport) system of equations, and it is shown that the per iteration CPU cost for the partial Newton method is not significantly greater than that of the Picard scheme.

Physically based modeling in catchment hydrology at 50: Survey and outlook

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.

Finite Element Approximation of the Diffusion Operator on Tetrahedra

Linear Galerkin finite element discretizations of the Laplace operator produce nonpositive stiffness coefficients for internal element edges of two-dimensional Delaunay triangulations. This property, also called the positive transmissibility (PT) condition, is a prerequisite for the existence of an M-matrix and ensures that nonphysical local extrema are not present in the solution. For tetrahedral elements, it has already been shown that the linear Galerkin approach does not in general satisfy the PT condition. We propose a modification of the three-dimensional Galerkin scheme that, if a Delaunay triangulation is used, satisfies the PT condition for internal edges and, if further conditions on the boundary are specified, yields an

Observation and modeling of catchment-scale solute transport in the hydrologic response: A tracer study

M

A tracer test in a shallow heterogeneous aquifer monitored via time-lapse surface electrical resistivity tomography

-matrix. The proposed approach can also be extended to the general diffusion operator with nonconstant or anisotropic coefficients.

Mesh locking effects in the finite volume solution of 2-D anisotropic diffusion equations

The coherent description of water flow and solute transport within heterogeneous hydrologic media (e.g., hillslopes or entire catchments) in response to external rainfall forcings represents a challenge in hydrological modeling. In this paper the mechanisms determining the mobilization and transport of solutes in soils through the paths of runoff formation are investigated by means of a tracer experiment conducted within an instrumented hillslope draining into a tributary of the Dese river basin (northeastern Italy). The response of the test catchment during a rainfall event occurred at the beginning of December 2006 has been analyzed by employing two different chemical tracers: nitrates from diffuse agricultural sources (NO 3-) and lithium from a point injection (Li+). Rainfall depths, streamflows, and pressure heads within different soil horizons were also collected. The observed hydro-chemical response of the Piovega Tre Comuni has been reproduced by a Mass Response Function solute transport model. Our results evidence the origins of the chromatographic effects seen in the systems response. In particular, the deep subsurface component of runoff is found to dominate the long-term behavior of the hydrograph and proves responsible for exporting relatively large amounts of solutes from soil. Experimental evidences and modeling results also suggest that the behavior of the breakthrough curve in stream waters may be strongly affected by the persistence of the rainfall forcing and by the relative magnitude of the rainfall-driven hydrologic contribution with respect to the background noise unrelated to soil moisture dynamics.