MODSharp: Regional-scale numerical model for quantifying groundwater flux and contaminant discharge into the coastal zone

In this paper the development of a quasi-three-dimensional numerical model that can be used for quantifying groundwater inputs and associated contaminant discharged from coastal aquifers into the coastal zone at a regional scale is presented. The present model is called MODSharp. In order to handle problems at a regional scale, the sharp interface approach which is used for conceptualising seawater intrusion, is applied to this model. This model can be used for the simulation of groundwater flow and contaminant transport in layered coastal aquifers at a regional scale. The method of characteristics is used to solve the advection-dispersion equation, which governs contaminant transport in coastal aquifers. In this study, MODSharp is used to investigate the influence of Sea Water Intrusion Interface (SWII) on the temporal and spatial variations of contaminant flux from coastal aquifers into the sea. Thus, a large number of simulations for different scenarios are performed. For the case that the land-ward boundary condition is constant head, it is shown that simplification and neglecting the effects of SWII causes an erroneous estimate of the speed of the contaminant plume movement towards sea. The value of hydraulic conductivity is shown to have a significant effect on the amount of discharged contaminant. Finally, it is concluded that over-simplification of the sea-ward boundary condition in numerical simulations, causes an incorrect estimate of temporal and spatial variations of the discharged contaminant into coastal water.     

Vectorized simulation of groundwater flow and streamline transport

We describe a modeling suite of Matlab functions for simulating nonpoint source (NPS) pollution in groundwater aquifers. The NPS model simulates groundwater flow and contaminant transport from a large array (order of 10² – 10⁷) of spatially distributed sources with time-varying pollution strength to a similarly large array of spatially distributed production wells (receptors) using the streamline transport approach. The code solves three equations: steady-state groundwater flow, particle tracking, and transient advection dispersion contaminant transport. The code performs convolution integration in its predictive step. Written in highly efficient vectorized form to avoid time consuming “for/while” loops, the code is also suitable for other groundwater flow and transport problems. The code is verified against analytical solutions and finite element software Comsol. An application illustrates 200 years of transient nitrate transport in the 2000 km² Tule River aquifer sub-basin of the Central Valley, California, with 9000 individual nitrate sources and 1920 wells.     

A complex variable boundary element method for a class of boundary value problems in anisotropic thermoelasticity

A boundary element method based on the Cauchy's integral formulae, called the complex variable boundary element method (CVBEM), is proposed for the numerical solution of boundary value problems governing plane thermoelastic deformations of anisotropic elastic bodies. The method is applicable for a wide class of problems which do not involve inertia or coupling effects and can be easily and efficiently implemented on the computer. It is applied to solve specific test problems.     

Numerical treatments of the interface discontinuity in solid-water mass transfer problems

The material discontinuity profile along the separation surface between the liquid and solid phases may represent a difficulty for the numerical solutions in studies of groundwater flow and contaminant transport in porous media and fractured systems. The present paper focuses on the dispersion of an hydrosoluble pollutant where the numerical procedures are required to describe the discontinuity of the concentration profile and at the same time preserve the continuity of the mass flux. The analysis shows that, whenever the use of shifted computational grid offers the possibility of a particularly elegant approach, the use of nonshifted computational grid should be preferred.     

Numerical simulation of contaminant transport in groundwater using software tools of $$r^3t$$

When a realistic modelling of radioactive contaminant transport in flowing groundwater is required, very large systems of coupled partial and ordinary differential equations can arise that have to be solved numerically. For that purpose, the software package \(r^3t\) is developed in which several advanced numerical methods are implemented to solve such models efficiently and accurately. Using software tools of \(r^3t\) one can treat successfully nontrivial mathematical problems like advection-dominated system with different retardation of transport for each component and with nonlinear Freundlich sorption and/or precipitation. Additionally, long time simulations on complex 3D geological domains using unstructured grids can be realized. In this paper we introduce and summarize the most important and novel features of numerical simulation for radioactive contaminant transport in porous media when using \(r^3t\).     

Flux-based method of characteristics for contaminant transport in flowing groundwater

A numerical method for treating advection-dominated contaminant transport in flowing groundwater is described. This method combines advantages of numerical discretizations by finite volume methods (like local mass conservation and the positivity of solutions) and by methods of characteristics (like larger time steps and reduced artificial numerical dispersion). For one-dimensional problems the method can produce equivalent algebraic systems as the finite volume Eulerian-Lagrangian localized adjoint method [13] and the flux-based modified method of characteristics [23] (and some other methods). An extension of the "flux-based methods of characteristics" for complex transport problems on multidimensional unstructured computational grids is the main contribution of this paper. Numerical results are included for a well established test example using a flux-based method of characteristics with aligned finite volumes.     

Coupling geological and numerical models to simulate groundwater flow and contaminant transport in fractured media

A new modeling approach is presented to improve numerical simulations of groundwater flow and contaminant transport in fractured geological media. The approach couples geological and numerical models through an intermediate mesh generation phase. As a first step, a platform for 3D geological modeling is used to represent fractures as 2D surfaces with arbitrary shape and orientation in 3D space. The advantage of the geological modeling platform is that 2D triangulated fracture surfaces are modeled and visualized before building a 3D mesh. The triangulated fractures are then transferred to the mesh generation software that discretizes the 3D simulation domain with tetrahedral elements. The 2D triangular fracture elements do not cut through the 3D tetrahedral elements, but they rather form interfaces with them. The tetrahedral mesh is then used for 3D groundwater flow and contaminant transport simulations in discretely fractured porous media. The resulting mesh for the 2D fractures and 3D rock matrix is checked to ensure that there are no negative transmissibilities in the discretized flow and transport equation, to avoid unrealistic results. To test the validity of the approach, flow and transport simulations for a tetrahedral mesh are compared to simulations using a block-based mesh and with results of an analytical solution. The fluid conductance matrix for the tetrahedral mesh is also analyzed and compared with known matrix values.     

A wireless sensor system for validation of real-time automatic calibration of groundwater transport models

In this paper, we present the use of a wireless sensor network in a lab for subsurface contaminant plume monitoring with the objective of automatic calibration of groundwater transport models. A tank configured to simulate an aquifer was used as a testbed, and a 2D model was created based on the setup. To simulate a contaminant plume, an ion tracer was injected into the tank. Sensor probes capable of detecting the plume were buried inside the tank, and wireless motes used to take readings from the sensors and relay data to a base station. More importantly, a run-time fault detection and diagnosis for abnormal sensor readings is designed and integrated into the data acquisition system. Further, an adaptive data collection technique is integrated that is able to provide evidence about the effectiveness of the groundwater transport model in use. Results from the tracer tests are presented, as well as lessons gained.     

A MATLAB method of lines template for transport equations

Many environmental problems involve diffusion and convection processes, which can be described by partial differential equations (PDEs). This paper will describe the development of a MATLAB template that generates a numerical solution to PDEs using the method of lines. The template will be applied to various unsaturated flow problems within soil physics to demonstrate the versatility of the method. In particular, the template will generate solutions for three cases (1) one-dimensional Richards’ equation for vertical infiltration; (2) coupled one-dimensional Richards’ equation and solute transport equation for horizontal water and contaminant flow; and (3) two-dimensional Richard’s equation for unsaturated flow over a complex geometry. Where possible, the results from the template will be compared against analytical solutions to determine the accuracy of the numerical solution. In addition, the paper will provide a discussion on possible extensions to the template and future directions.     

Application of the CVBEM to non-uniform St. Venant torsion

The complex-variable boundary-element method, or CVBEM, is used to approximate the stress distribution associated with non-uniform St. Venant torsion problems. By specifying either the normal- or tangential-force equilibrium equation in terms of the warping function or its conjugate, a Laplace equation results which is immediately tractable by using the CVBEM. A comparison of modeling results to known solutions indicates that the modeling technique is a useful approach for the estimation of interior stresses as well as for the evaluation of modeling error by means of an approximate boundary determined by the CVBEM approximation function.     

Modeling of 2D density-dependent flow and transport in porous media using finite volume method

A two-dimensional finite volume unstructured mesh method (FVUM) based on a triangular mesh is developed for modeling density-dependent flow and transport through saturated-unsaturated porous media. The combined flow and transport model can handle a wide range of real-word problems, including the simulation of flow alone, contaminant transport alone, and combined flow and transport. Saltwater intrusion problems and instability caused by denser water on the top were investigated in this paper. Because the fundamental mechanism causing saltwater intrusion most likely is caused by density-induced convection and dispersion, the developed model is used to assess the interplay between density-driven flow and subsurface media through which the saltwater intrusion occurs.The mathematical formulation of the model is comprised of fluid flow and solute transport equations, coupled by fluid density. In the specific case of saltwater intrusion and unstable brine transport problems, this set of governing equation is non-linear and requires iterative methods to solve them simultaneously. Three case studies, which include a wide range of physical conditions, are used for verification of presented model and comparison with previously published solutions from other researchers presents encouraging agreements.     

Numerical Modeling of Degenerate Equations in Porous Media Flow

In this paper is introduced a new numerical formulation for solving degenerate nonlinear coupled convection dominated parabolic systems in problems of flow and transport in porous media by means of a mixed finite element and an operator splitting technique, which, in turn, is capable of simulating the flow of a distinct number of fluid phases in different porous media regions. This situation naturally occurs in practical applications, such as those in petroleum reservoir engineering and groundwater transport. To illustrate the modelling problem at hand, we consider a nonlinear three-phase porous media flow model in one- and two-space dimensions, which may lead to the existence of a simultaneous one-, two- and three-phase flow regions and therefore to a degenerate convection dominated parabolic system. Our numerical formulation can also be extended for the case of three space dimensions. As a consequence of the standard mixed finite element approach for this flow problem the resulting linear algebraic system is singular. By using an operator splitting combined with mixed finite element, and a decomposition of the domain into different flow regions, compatibility conditions are obtained to bypass the degeneracy in order to the degenerate convection dominated parabolic system of equations be numerically tractable without any mathematical trick to remove the singularity, i.e., no use of a parabolic regularization. Thus, by using this procedure, we were able to write the full nonlinear system in an appropriate way in order to obtain a nonsingular system for its numerical solution. The robustness of the proposed method is verified through a large set of high-resolution numerical experiments of nonlinear transport flow problems with degenerating diffusion conditions and by means of a numerical convergence study.     

Application of MATLAB and Python optimizers to two case studies involving groundwater flow and contaminant transport modeling

One approach for utilizing geoscience models for management or policy analysis is via a simulation-based optimization framework—where an underlying model is linked with an optimization search algorithm. In this regard, MATLAB and Python are high-level programming languages that implement numerous optimization routines, including gradient-based, heuristic, and direct-search optimizers. The ever-expanding number of available algorithms makes it challenging for practitioners to identify optimizers that deliver good performance when applied to problems of interest. Thus, the primary contribution of this paper is to present a series of numerical experiments that investigated the performance of various MATLAB and Python optimizers. The experiments considered two simulation-based optimization case studies involving groundwater flow and contaminant transport. One case study examined the design of a pump-and-treat system for groundwater remediation, while the other considered least-squares calibration of a model of strontium (Sr) transport. Using these case studies, the performance of 12 different MATLAB and Python optimizers was compared. Overall, the Hooke-Jeeves direct search algorithm yielded the best performance in terms of identifying least-cost and best-fit solutions to the design and calibration problems, respectively. The IFFCO (implicit filtering for constrained optimization) direct search algorithm and the dynamically dimensioned search (DDS) heuristic algorithm also consistently yielded good performance and were up to 80% more efficient than Hooke-Jeeves when applied to the pump-and-treat problem. These results provide empirical evidence that, relative to gradient- and population-based alternatives, direct search algorithms and heuristic variants, such as DDS, are good choices for application to simulation-based optimization problems involving groundwater management.     

Density-driven flow and solute transport problems. A 2-D numerical model based on the network simulation method

The network simulation method, based on the formal equivalence between physical systems and electrical networks, solves numerical problems of relatively mathematical complexity in a versatile, efficient and computationally fast way. In this paper, the method is applied for the first time to the design of a general purpose model for simulating two-dimensional transient density-driven flow and solute transport through porous media, a mathematical model made up by coupled, nonlinear differential equations. Using the Boussinesq approximation and the stream function formulation, the model is used to solve two typical problems related with groundwater flows. Isochlor concentration and stream function curves are presented and successfully compared with those of other authors. Simulation is carried out using the digital computer program Pspice with relatively low computing times.     

Superconvergence of a Mixed Covolume Method for Elliptic Problems

We consider a mixed covolume method for a system of first order partial differential equations resulting from the mixed formulation of a general self-adjoint elliptic problem with a variable full diffusion tensor. The system can be used to model the transport of a contaminant carried by a flow. We use the lowest order Raviart-Thomas mixed finite element space. We show the first order convergence in L ² norm and the superconvergence in certain discrete norms both for the pressure and velocity. Finally some numerical examples illustrating the error behavior of the scheme are provided. Supported by the National Natural Science Foundation of China under grant No. 10071044 and the Research Fund of Doctoral Program of High Education by State Education Ministry of China.     

lr2dinv: A finite-difference model for inverse analysis of two-dimensional linear or radial groundwater flow

We have developed a program for inverse analysis of two-dimensional linear or radial groundwater flow problems. The program, lr2dinv, uses standard finite difference techniques to solve the groundwater flow equation for a horizontal or vertical plane with heterogeneous properties. In radial mode, the program simulates flow to a well in a vertical plane, transforming the radial flow equation into an equivalent problem in Cartesian coordinates. The physical parameters in the model are horizontal or x-direction hydraulic conductivity, anisotropy ratio (vertical to horizontal conductivity in a vertical model, y-direction to x-direction in a horizontal model), and specific storage. The program allows the user to specify arbitrary and independent zonations of these three parameters and also to specify which zonal parameter values are known and which are unknown. The Levenberg–Marquardt algorithm is used to estimate parameters from observed head values. Particularly powerful features of the program are the ability to perform simultaneous analysis of heads from different tests and the inclusion of the wellbore in the radial mode. These capabilities allow the program to be used for analysis of suites of well tests, such as multilevel slug tests or pumping tests in a tomographic format. The combination of information from tests stressing different vertical levels in an aquifer provides the means for accurately estimating vertical variations in conductivity, a factor profoundly influencing contaminant transport in the subsurface.     

Solution of contaminant transport with equilibrium and non-equilibrium adsorption

In this paper we present a numerical approximation scheme for the solution of contaminant transport problems with diffusion and adsorption in equilibrium and non-equilibrium mode. The method is based on time stepping and operator splitting. The non-linear transport is solved semi-analytically via the multiple Riemann problem, the non-linear diffusion by a finite volume method and by Newton’s type of linearisation, and finally the reaction part, incorporating the non-equilibrium adsorption, is transformed to an integral equation which is solved numerically using time discretization. Various results of numerical experiments are shown, and the method is applied to the practical dual-well problem.     

A combined finite and infinite element approach for modeling spherically symmetric transient subsurface flow

This paper presents a finite element-infinite element coupling approach for modeling a spherically symmetric transient flow problem in a porous medium of infinite extent. A finite element model is used to examine the flow potential distribution in a truncated bounded region close to the spherical cavity. In order to give an appropriate artificial boundary condition at the truncated boundary, a transient infinite element, that is developed to describe transient flow in the exterior unbounded domain, is coupled with the finite element model. The coupling procedure of the finite and infinite elements at their interface is described by means of the boundary integro-differential equation rather than through a matrix approach. Consequently, a Neumann boundary condition can be applied at the truncated boundary to ensure the C¹-continuity of the solution at the truncated boundary. Numerical analyses indicate that the proposed finite element-infinite element coupling approach can generate a correct artificial truncated boundary condition to the finite element model for the unbounded flow transport problem.     

A new package for simulating periodic boundary conditions in MODFLOW and SEAWAT

Modeling of coastal groundwater systems is a challenging problem due to their highly dynamic boundary conditions and the coupling between the equations for groundwater flow and solute transport. A growing number of publications on aquifers subject to tides have demonstrated various modeling approaches, ranging from analytical solutions to comprehensive numerical models. The United States Geological Survey code SEAWAT has been a popular choice in studies of this type. Although SEAWAT allows the incorporation of time-variant boundary conditions, the implementation of tidal boundaries is not straightforward, especially when a seepage face develops during falling tide. Here, a new package is presented, called the periodic boundary condition (PBC) package, that can be incorporated into MODFLOW and SEAWAT to overcome the difficulties encountered with tidal boundaries. It dynamically updates the boundary conditions for head and concentration during the simulation depending on a user-defined tidal signal and allows for the development of a seepage face. The package has been verified by comparing it to four different published models of tidally influenced groundwater systems of varying complexity. Excellent agreement was obtained in all cases. The new package is an important extension to the existing capabilities of MODFLOW and SEAWAT with respect to simulating periodic boundary conditions.     

Flux-based methods of characteristics for coupled transport equations in porous media

Mathematical models of transport of radioactive contaminants in flowing groundwater involve large systems of coupled advection dominated transport equations. High-resolution explicit finite volume methods, if applied to advective part of model and combined with appropriate numerical methods for diffusion-dispersion-reaction part, can offer precise and monotone numerical solutions, but they require small time steps. This paper describes Flux-Based Methods Of Characteristics that are extension of explicit finite volume methods, that have no restriction on time steps and that produce numerical solutions with valid discrete minimum and maximum principle. Such particular method was implemented in software package R³T (Retardation, Reaction, Radionuclides and Transport) and it was used successfully to solve large systems of coupled transport equations with different retardation factors for transport.