Numerical Modeling of Coupled Ground and Surface Water Flow and Transport

Funded by the National Science Foundation, Division of Mathematical Sciences

Principal Investigators

Clint Dawson, Professor of Aerospace Engineering and Engineering Mechanics, The University of Texas at Austin clint@ices.utexas.edu

Mary F. Wheeler, Professor of Aerospace Engineering and Engineering Mechanics and Professor of Petroleum and Geosystems Engineering, The University of Texas at Austin mfw@ices.utexas.edu

Objectives

Sustainable management and protection of water resources is one of the key problems of the 21st century. Historically, ground and surface water have been modeled in isolation, without properly accounting for the interactions between the two. Recent field studies, however, in the United States and abroad, point to the necessity for properly accounting for the coupling between the two regimes.

The objective of this proposal is to develop advanced numerical methods for coupling surface water flow with complex subsurface hydrosystems to accurately simulate flow, transport and reaction processes over large space and time scales. As the relevant processes strongly differ in the various subdomains of a hydrosystem, different model concepts must be chosen for the subdomains, and special coupling methods must be developed to account for the interaction processes between these subdomains as well as for moving subdomain boundaries. The resulting numerical methods will involve efficient discretization techniques well adapted to the dominant processes in each subdomain, as well as fast and robust solvers. Representative benchmark problems based on laboratory and field experiments will be formulated and will validate the numerical schemes.

The proposed work includes:

• Development of coupling conditions for mass and momentum in ground water and surface water regions. This includes mathematical formulation and algorithmic and software tool development between different domains and models.
• Formulation of numerical techniques and analysis of algorithms.
• Design of relevant benchmarks based on laboratory and field data for verficiation and evaluation of numerical algorithms.
• Educational and professional outreach and cross-disciplinary training of future scientists and engineers through workshops, video and web-based technologies.

Recent papers resulting from this project

• C. Dawson, Analysis of discontinuous finite element methods for ground water/surface water coupling, SIAM Journal on Numerical Analysis, 52, pp. 63-88, 2006.
• H. Li, M. Farthing, C. Dawson, and C.T. Miller, Local discontinuous Galerkin approximations to Richards' equation, Advances in Water Resources, in press.
• E.J. Kubatko, J.J. Westerink and C. Dawson, hp Discontinous Galerkin methods for advection-dominated problems in shallow water, Computer Methods in Applied Mechanics and Engineering, 196, pp. 437-451, 2006.
• V. Aizinger and C. Dawson, The local discontinuous Galerkin method for three dimensional shallow water flow, Computer Methods in Applied Mechanics and Engineering, 196, pp. 734-746, 2007.
• E.J. Kubatko, J.J. Westerink and C. Dawson, Semi-discrete discontinuous Galerkin methods and stage exceeding order, strong stability preserving Runge-Kutta time discretizations, Journal of Computational Physics, in press.

Major presentations related to this work

• Westerink, J., Kubatko, E. and Dawson, C., ``An unstructured grid morphodynamic model with a discontinuous Galerkin method for bed evolution,'' The 3rd International Workshop on Unstructured Grid Numerical Modelling of Coastal Shelf and Ocean Flows, Toulouse, France, Sept. 2004.
• Dawson, C., ``Discontinuous Galerkin Finite Element Methods for Shallow Water Flows,'' Workshop on Modelling and Computation in Environmental Science, Hohenwart, Germany, October 2004.
• Dawson, C., ``Discretization techniques for coupled flow and transport,'' Center for Applied Scientific Computing, Lawrence Livermore National Laboratory, Livermore, CA, November 2004.
• Li, H., Farthing, M., Dawson, C. and Miller, C., ``Local discontinuous Galerkin and variable step size, variable order time integration for Richards equation,'' American Geophysical Union, San Francisco, CA, December 2004.
• Dawson, C., ``Recent advances in surface water modeling,'' Plenary lecture, 7th SIAM Conference on Mathematical and Computational Issues in the Geosciences, Avignon, France, June 2005.
• Dawson, C., ``Numerical methods for ground water/surface water coupling,'' 7th SIAM Conference on Mathematical and Computational Issues in the Geosciences, Avignon, France, June 2005.
• Dawson, C., ``Discontinuous Galerkin methods for 2-D and 3-D Shallow Water Equations,'' Keynote lecture, 8th U.S. National Congress on Computational Mechanics, Austin, July 2005.
• Kubatko, E., Westerink J., and Dawson, C., ``hp discontinuous Galerkin methods for shallow water flow and transport,'' 8th U.S. National Congress on Computational Mechanics, Austin, July 2005.
• Wheeler, M.F., ``Mathematical issues in multiphysics couplings,'' Modeling and Computation in Environmental Science, Hohenwart Forum, Hohenwart, Germany, October, 2004.
• Wheeler, M.F., ``Physics-driven algorithms for modeling subsurface phenomena,'' Thirteenth Conference on Finite Elements for Flow Problems, University of Wales Swansea, Swansea, Wales, April, 2005.
• Wheeler, M.F., ``Physics-driven algorithms for modeling subsurface phenomena,'' SIAM Conference on Mathematical and Computational Issues in the Geosciences, Avignon, France, June, 2005.
• Dawson, C., ``Numerical simulation of coupled ground water/surface water flow and transport,'' AMS/SIAM/MAA Joint Mathematics Meeting, San Antonio, January 2006.
• Dawson, C., ``DG methods for shallow water modeling,'' Department of Defense PET workshop, Engineering Research and Development Center, Vicksburg, MS, Feb. 2006.
• Dawson, C., ``Modeling coastal hydrodynamics and hurricane storm surges,''NSF-CBMS Conference, University of Nevada-Las Vegas, May 2006.
• Dawson, C., ``Discontinuous and continuous Galerkin methods for shallow water and hurricane storm surges,'' MAFELAP, Brunel University, United Kingdom, June 2006.
• Dawson, C., ``Discontinuous Galerkin methods for coupled ground water/surface water flow and transport,'' World Congress on Computational Mechanics 7, Los Angeles, July 2006.
• Bunya, S., Westerink, J., Kubatko, E. and Dawson, C., `` A new wetting and drying algorithm for discontinuous Galerkin solutions to the shallow water equations,'' World Congress on Computational Mechanics 7, Los Angeles, July 2006.

Results

Development and analysis of a coupled shallow water/vadose zone model

We have formulated and implemented discontinuous Galerkin (DG) based approaches for coupling the depth-averaged shallow water equations with ground water flow equations, assumed to be modeled by Richards' equation. The coupling of the two models is realized through the assumption of continuity of water pressure and water flux through the interface between the surface and subsurface. These conditions are imposed weakly in the formulation.

Under the assumption of fully saturated flow, which removes the nonlinearities in Richards' equation, Dawson has developed a priori error estimates for the coupled approach in a continuous-time setting. In this approach, a local discontinuous Galerkin (LDG) method is used to model ground water flow. These estimates are suboptimal in ground water pressure, but optimal in velocity. The loss of order is due to the nonlinearities in the shallow water model and the approximation of pressure boundary conditions in the ground water model. Note that these estimates make very weak assumptions on the relationship between the mesh used in the ground water domain and that used in the surface water domain; in particular, these meshes need not align at the coupling interface. Furthermore, the methods extend to the case where a mixed finite element method is used in the ground water region.

A code which simulates coupled surface/subsurface flow based on the model described above has been developed. The code uses a DG discretization in the surface water domain, and an LDG discretization of Richards' equation in the subsurface. As part of this effort, we have extensively studied features of the LDG method for Richards' equation. In collaboration, with C.T. Miller and his group at the University of North Carolina-Chapel Hill, we have investigated implementation issues related to the LDG method, including the choice of numerical fluxes and the use of adaptive, higher-order BDF time-stepping approaches. Extensive one-dimensional studies were performed in 2006. We are continuing our collaboration with Miller and his group for similar studies in higher spatial dimensions.

In collaboration with student Marco Iglesias-Hernandez, we have performed extensive verification and validation studies of the LDG method for Richards' equation in both one and two space dimensions. These include comparisons to well-known analytical solutions and spatial/time-step refinement studies. Furthermore, we have applied the code to infiltration problems in two-dimensional heterogeneous media to study the effects of entry pressure. For example, we have considered a two-dimensional porous medium consisting of a coarse sand with high conductivity and low entry pressure, and a fine sand with low conductivity and high entry pressure. The domain is shown here, with the hatched region illustrating the fine sand block. In the following scenarios, the domain is initially unsaturated and water infiltrates through the top boundary of the domain. The question is whether or not the water infiltrates into the find sand region. In the following figures, we show results for two different choices of relative permeabilities and entry pressures in the fine and coarse sand regions. In the first case, the fine sand region has low relative permeability, and as can been seen in the following movie of water saturation profiles up to time .031 days (click here), the water flows around the fine sand region. In the second scenario, the fine sand has lower entry pressure and thus water is able to infiltrate into the fine sand region. The following movie illustrates this effect.

These studies were very computationally intensive, and as Richards' equation involves the solution of large nonlinear systems of equations, which in turn leads to the repeated solution of large linear systems of equations, we have worked extensively in improving the linear solvers in the code. This work has included incorporating more efficient, sparse-storage data structures and ILU preconditioners.

With these improvements to the efficiency of the ground water model, we are currently performing more in-depth studies of surface water/ground water coupling under different scenarios. For example, see the following movie. In this model we have a one-dimensional shallow water channel overlying a two-dimensional vadose zone. The dimensions of the domain are approximately 100 cm by 90-100 cm, with the top of the domain varying in height. The initial conditions are unsaturated in the ground water domain with a very thin layer (approximately 1 cm) of still water in the shallow water domain at the top. After the initial time, a surge of water moves into the shallow water channel and water begins to infiltrate into the ground water region. The movie shows the change in hydraulic head from initial time to time 200 minutes.

New algorithms for modeling flow in the subsurface

Mary Wheeler and her students, postdocs and several collaborators have recently investigated new approaches for modeling subsurface flows using multipoint flux approximation (MPFA) methods and DG methods. With Ivan Yotov from University of Pittsburgh, Wheeler has analyzed MPFA methods for general quadrilateral and hexahedral elements by showing their equivalence to a BDM mixed finite element method. This approach is useful for this project because to accurately couple subsurface/surface water flows requires modeling geometrically complex interfaces between the two regimes. These are the first convergence proofs to be obtained for the MPFA method.

Wheeler and several collaborators have been investigating multiscale DG methods for coupling multiple physics and/or multiple domains (which may have their own grid and time-step) through interfaces. Here one must develop appropriate transmission or matching conditions on the interface. One approach is the use of mortar finite element methods. In this approach, the simulation domain is first divided into a series of non-overlapping subdomains (blocks) depending on efficient representation of the problem such as geometric irregularities and variations of physical, chemical and biological properties. Each block may have specific physical, chemical and biological processes in a more regular subregion, compared with the original problem. A local physical model is solved in each block using the most efficient algorithm. We fill the interfaces between blocks with mortars elements of a finite element space called the mortar space. All blocks are then coupled by weakly forcing continuity of fluxes (or other relevant quantities) on the block interfaces using the mortar spaces. Mortar spaces provide a flexible way to couple different discretizations on non-matching multi-block grids as well as different physical models on different parts of the simulation domain. They also lend themselves to multiscale resolution, as one can couple highly refined regions where one wants to capture fine scale phenomena, with more coarsely refined regions through the use of a mortar space. Recently the couplings over subdomains of different discretizations of DG with DG using higher order mortars have been shown to be multiscale methods yielding optimal global rates of convergence. In addition the DG couplings provide a promising approach as domain decomposition solvers.

Wheeler and former student Owen Eslinger have developed new DG methods for modeling air-water flow. This model is based on a full two phase formulation with compressibility and capillary pressure included, and is manipulated to form a pressure equation and a saturation equation. The Kirchoff transformation is then applied to the saturation equation. An interior penalty Galerkin method is applied to the pressure equation, and the LDG method is applied to the transformed saturation equation. This work was the basis for Eslinger's Ph.D. thesis, which was completed in 2005.

Improvements to DG methods for 2D and 3D shallow water flow

Dawson has continued his collaboration with Joannes Westerink and Shintaro Bunya at the University of Notre Dame, and Ethan Kubatko, a recent Ph.D. graduate at Notre Dame and current ICES post-doc at UT Austin, on the development of DG methods for shallow water flow. Recent developments include the implementation of improved time-stepping methods based on strong-stability-preserving Runge-Kutta time discretizations (SSP-RK), and the development of mass-conserving methods for modeling wetting and drying phenomena. These algorithms are being incorporated into a newly developed, h-p adaptive DG code for depth-integrated shallow water models, which will be used in future ground water/surface water coupling studies. Through the use of stage-exceeding-order SSP-RK methods, substantial improvements in the efficiency of SSP-RK time discretizations for DG-based shallow water models have been observed. Normally a k-th order RK method requires computing k stages. By increasing the number of stages to s > k, methods have been derived which enlarge the region of stability. The cost required to compute additional stages can be offset by increases in the time-step. Our findings show for example that for k=2 and k=3, the optimal number of stages in terms of efficiency is s=3 and s=5, respectively, while for k=4, using s > 4 stages can improve efficiency up to 25% over s=4. These methods have been implemented and are currently being used in our 2D shallow water code.

Dawson is also continuing collaboration with former student Vadym Aizinger on three-dimensional shallow water modeling. A stability analysis for the LDG method applied to the three-dimensional shallow water flow equations has been dervied, assuming hydrostatic pressure, including the effects of modeling the free surface. We are currently collaborating with Aizinger in extending this model to handle turbulent and density-driven (baroclinic) flows. We plan to use this model in future ground water/surface water couplings where the surface water circulation may be driven by vertical variations in salinity and/or temperature. In these cases, traditional depth-averaged models are not appropriate.

Analysis and numerical methods for wetland flow models

With student Mauricio Santillana, Dawson has been investigating models for shallow water circulation in wetlands. The movitation for this work is that wetlands have unique features, and thus traditional shallow water models may not be appropriate in thes regions. Our plan is to study models for wetland flow analytically and numerically, and to further study the coupling of these models to depth-integrated shallow water and ground water models. Wetland flow models assume the acceleration terms in the horizontal momentum equations are neglible, that flow is dominated by gravitational forces and friction, and thus horizontal water fluxes and water height can be related using Manning's equation or similar alternative expressions. The resulting flow model is a nonlinear diffusion equation. Santillana is studying this model both analytically and numerically. Preliminary existence results for solutions have been proved, and stability of continuous and discontinuous Galerkin approximations have been demonstrated. We are currently investigating a priori error estimates for this model in the context of either continuous or discontinuous approximations; but these estimates are complicated by the fact that the diffusion coefficient is not bounded below in general, and possibly nonlinear in both the solution and its gradient. A numerical model has also been developed, and preliminary verification studies have shown that the model reproduces results reported in the literature for a standard one-dimensional test problem. We expect that this work will form a large part of Santinilla's Ph.D. thesis.

<< Back to Current Funded Projects