Modeling Steady Sea Water Intrusion with Single‐Density Groundwater Codes

Steady interface flow in heterogeneous aquifer systems is simulated with single‐density groundwater codes by using transformed values for the hydraulic conductivity and thickness of the aquifers and aquitards. For example, unconfined interface flow may be simulated with a transformed model by setting the base of the aquifer to sea level and by multiplying the hydraulic conductivity with 41 (for sea water density of 1025 kg/m³). Similar transformations are derived for unconfined interface flow with a finite aquifer base and for confined multi‐aquifer interface flow. The head and flow distribution are identical in the transformed and original model domains. The location of the interface is obtained through application of the Ghyben‐Herzberg formula. The transformed problem may be solved with a single‐density code that is able to simulate unconfined flow where the saturated thickness is a linear function of the head and, depending on the boundary conditions, the code needs to be able to simulate dry cells where the saturated thickness is zero. For multi‐aquifer interface flow, an additional requirement is that the code must be able to handle vertical leakage in situations where flow in an aquifer is unconfined while there is also flow in the aquifer directly above it. Specific examples and limitations are discussed for the application of the approach with MODFLOW. Comparisons between exact interface flow solutions and MODFLOW solutions of the transformed model domain show good agreement. The presented approach is an efficient alternative to running transient sea water intrusion models until steady state is reached.     

Discontinuous Steady‐State Analytical Solutions of the Boussinesq Equation and Their Numerical Representation by Modflow

Known analytical solutions of groundwater flow equations are routinely used for verification of computer codes. However, these analytical solutions (e.g., the Dupuit solution for the steady‐state unconfined unidirectional flow in a uniform aquifer with a flat bottom) represent smooth and continuous water table configurations, simulating which does not pose any significant problems for the numerical groundwater flow models, like MODFLOW. One of the most challenging numerical cases for MODFLOW arises from drying‐rewetting problems often associated with abrupt changes in the elevations of impervious base of a thin unconfined aquifer. Numerical solutions of groundwater flow equations cannot be rigorously verified for such cases due to the lack of corresponding exact analytical solutions. Analytical solutions of the steady‐state Boussinesq equation, associated with the discontinuous water table configurations over a stairway impervious base, are presented in this article. Conditions resulting in such configurations are analyzed and discussed. These solutions appear to be well suited for testing and verification of computer codes. Numerical solutions, obtained by the latest versions of MODFLOW (MODFLOW‐2005 and MODFLOW‐NWT), are compared with the presented discontinuous analytical solutions. It is shown that standard MODFLOW‐2005 code (as well as MODFLOW‐2000 and older versions) has significant convergence problems simulating such cases. The problems manifest themselves either in a total convergence failure or erroneous results. Alternatively, MODFLOW‐NWT, providing a good match to the presented discontinuous analytical solutions, appears to be a more reliable and appropriate code for simulating abrupt changes in water table elevations.     

Application of an integrated surface water‐groundwater model to multi‐aquifers modeling in Choushui River alluvial fan,Taiwan

This research proposes a combination of SWAT and MODFLOW, MD‐SWAT‐MODFLOW, to address the multi‐aquifers condition in Choushui River alluvial fan, Taiwan. The natural recharge and unidentified pumping/recharge are separately estimated. The model identifies the monthly pumping/recharge rates in multi‐aquifers so that the daily streamflow can be simulated correctly. A multi‐aquifers condition means a subsurface formation composed of at least the unconfined aquifer, the confined aquifer, and an in‐between aquitard. In such a case, the variation of groundwater level is related to pumping/recharge activities in vertically adjacent aquifer and the river‐aquifer interaction. Both factors in turn affect the streamflow performance. Results show that MD‐SWAT‐MODFLOW performs better than SWAT alone in terms of simulated streamflow, especially during low flow period, when pumping/recharge rates are properly estimated. A sensitivity analysis of individual parameter suggests that the vertical leakance may be the most sensitive among all investigated MODFLOW parameters in terms of the estimated pumping/recharge among aquifers, and the Latin‐Hypercube‐One‐factor‐At‐a‐Time sensitivity analysis indicates that the hydraulic conductivity of channel is the most sensitive to the model performance. It also points out the necessity to simultaneously estimate pumping/recharge rates in multi‐aquifers. The estimated net pumping rate can be treated as a lower bound of the actual local pumping rate. As a whole, the model provides the spatio‐temporal groundwater use, which gives the authorities insights to manage groundwater resources. Copyright © 2012 John Wiley & Sons, Ltd.     

Analytical methods for transient flow to a well in a confined-unconfined aquifer

Concurrent existence of confined and unconfined zones of an aquifer can arise owing to ground water withdrawal by pumping. Using Girinskii's potential function, Chen (1974, 1983) developed an approximate analytical solution to analyze transient ground water flow to a pumping well in an aquifer that changes from an initially confined system to a system with both unconfined and confined regimes. This article presents the details of the Chen model and then compares it with the analytical model developed by Moench and Prickett (1972) for the same problem. Hypothetical pumping test examples in which the aquifer undergoes conversion from confined to water table conditions are solved by the two analytical models and also a numerical model based on MODFLOW. Comparison of the results suggests that the solutions of the Chen model give better results than the Moench and Prickett model except when the radial distance is very large or aquifer thickness is large compared with drawdown.     

Incorporating Super‐Diffusion due to Sub‐Grid Heterogeneity to Capture Non‐Fickian Transport

Numerical transport models based on the advection‐dispersion equation (ADE) are built on the assumption that sub‐grid cell transport is Fickian such that dispersive spreading around the average velocity is symmetric and without significant tailing on the front edge of a solute plume. However, anomalous diffusion in the form of super‐diffusion due to preferential pathways in an aquifer has been observed in field data, challenging the assumption of Fickian dispersion at the local scale. This study develops a fully Lagrangian method to simulate sub‐grid super‐diffusion in a multidimensional regional‐scale transport model by using a recent mathematical model allowing super‐diffusion along the flow direction given by the regional model. Here, the time randomizing procedure known as subordination is applied to flow field output from MODFLOW simulations. Numerical tests check the applicability of the novel method in mapping regional‐scale super‐diffusive transport conditioned on local properties of multidimensional heterogeneous media.     

A Stable and Efficient Numerical Algorithm for Unconfined Aquifer Analysis

The nonlinearity of equations governing flow in unconfined aquifers poses challenges for numerical models, particularly in field-scale applications. Existing methods are often unstable, do not converge, or require extremely fine grids and small time steps. Standard modeling procedures such as automated model calibration and Monte Carlo uncertainty analysis typically require thousands of model runs. Stable and efficient model performance is essential to these analyses. We propose a new method that offers improvements in stability and efficiency and is relatively tolerant of coarse grids. It applies a strategy similar to that in the MODFLOW code to the solution of Richard's equation with a grid-dependent pressure/saturation relationship. The method imposes a contrast between horizontal and vertical permeability in gridblocks containing the water table, does not require "dry" cells to convert to inactive cells, and allows recharge to flow through relatively dry cells to the water table. We establish the accuracy of the method by comparison to an analytical solution for radial flow to a well in an unconfined aquifer with delayed yield. Using a suite of test problems, we demonstrate the efficiencies gained in speed and accuracy over two-phase simulations, and improved stability when compared to MODFLOW. The advantages for applications to transient unconfined aquifer analysis are clearly demonstrated by our examples. We also demonstrate applicability to mixed vadose zone/saturated zone applications, including transport, and find that the method shows great promise for these types of problem as well.     

Developing Potentiometric Surfaces and Flow Fields with a Head‐Specified MODFLOW Model

Investigating changes in an aquifer system often involves comparison of observed heads from different synoptic measurements, generally with potentiometric surfaces developed by hand or a statistical approach. Alternatively, head‐specified MODFLOW models, in which constant head cells simulate observed heads, generate gridded potentiometric surfaces that explicitly account for Darcy's Law and mass balance. We developed a transient head‐specified MODFLOW model for the stratified Cambrian‐Ordovician sandstone aquifer system of northeastern Illinois to analyze flow within its 275 m deep cone of depression. Potentiometric surfaces were developed using static heads from production wells regardless of open interval; hence assuming no vertical head difference. This assumption was tested against steady‐state, head‐specified models of each sandstone strata for 1980 and 2014. The results indicate that the original conceptual model was appropriate in 1980 but not 2014, where a vertical head difference had developed at the center of the cone of depression. For earlier years, when the head difference was minimal, the transient head‐specified model compared well with a traditional, flow‐specified model. In later years, the transient head‐specified model overestimated removal of water from storage. MODFLOW facilitates the development of a time‐series of potentiometric surfaces and can easily be modified to test the impacts of different conceptual models, such as assumptions on vertical head differences. For this study of a deep confined aquifer, MODFLOW also offers advantages in generating potentiometric surfaces and flow fields over statistical interpolation techniques, although future research is needed to assess its performance in other settings.     

Simulating Seepage into Mine Shafts and Tunnels with MODFLOW

In cases when an equivalent porous medium assumption is suitable for simulating groundwater flow in bedrock aquifers, estimation of seepage into underground mine workings (UMWs) can be achieved by specifying MODFLOW drain nodes at the contact between water bearing rock and dewatered mine openings. However, this approach results in significant numerical problems when applied to simulate seepage into an extensive network of UMWs, which often exist at the mine sites. Numerical simulations conducted for individual UMWs, such as a vertical shaft or a horizontal drift, showed that accurate prediction of seepage rates can be achieved by either applying grid spacing that is much finer than the diameter/width of the simulated openings (explicit modeling) or using coarser grid with cell sizes exceeding the characteristic width of shafts or drifts by a factor of 3. Theoretical insight into this phenomenon is presented, based on the so-called well-index theory. It is demonstrated that applying this theory allows to minimize numerical errors associated with MODFLOW simulation of seepage into UMWs on a relatively coarse Cartesian grid. Presented examples include simulated steady-state groundwater flow from homogeneous, heterogeneous, and/or anisotropic rock into a vertical shaft, a horizontal drift/cross-cut, a ramp, two parallel drifts, and a combined system of a vertical shaft connected to a horizontal drift.     

Use of General Purpose Graphics Processing Units with MODFLOW

To evaluate the use of general‐purpose graphics processing units (GPGPUs) to improve the performance of MODFLOW, an unstructured preconditioned conjugate gradient (UPCG) solver has been developed. The UPCG solver uses a compressed sparse row storage scheme and includes Jacobi, zero fill‐in incomplete, and modified‐incomplete lower‐upper (LU) factorization, and generalized least‐squares polynomial preconditioners. The UPCG solver also includes options for sequential and parallel solution on the central processing unit (CPU) using OpenMP. For simulations utilizing the GPGPU, all basic linear algebra operations are performed on the GPGPU; memory copies between the central processing unit CPU and GPCPU occur prior to the first iteration of the UPCG solver and after satisfying head and flow criteria or exceeding a maximum number of iterations. The efficiency of the UPCG solver for GPGPU and CPU solutions is benchmarked using simulations of a synthetic, heterogeneous unconfined aquifer with tens of thousands to millions of active grid cells. Testing indicates GPGPU speedups on the order of 2 to 8, relative to the standard MODFLOW preconditioned conjugate gradient (PCG) solver, can be achieved when (1) memory copies between the CPU and GPGPU are optimized, (2) the percentage of time performing memory copies between the CPU and GPGPU is small relative to the calculation time, (3) high‐performance GPGPU cards are utilized, and (4) CPU‐GPGPU combinations are used to execute sequential operations that are difficult to parallelize. Furthermore, UPCG solver testing indicates GPGPU speedups exceed parallel CPU speedups achieved using OpenMP on multicore CPUs for preconditioners that can be easily parallelized.     

Unconfined linear flow to a horizontal well

The validity of a previously proposed but untested modification to equations for flow to a horizontal well is assessed using a specially developed finite-difference model. This modification extends confined flow equations to allow the head in the well and the saturated depth at the well to be estimated in unconfined conditions. The study is limited to the case of two-dimensional flow with no flow in the direction parallel to the line of the well. The results show that the modified equations for both a finite unconfined aquifer and, by inference, an infinite unconfined aquifer are adequately accurate for practical application.     

Improved water table dynamics in MODFLOW

The standard formulation of a block-centered finite-difference model, such as MODFLOW, uses the center of the cell as the location of a cell node. Simulations of a dynamic water table can be improved if the node of a cell containing the water table is located at the water table rather than at the center of the cell. The LPF package of MOD-FLOW-2000 was changed to position a cell's node at the water table in convertible cells with a water table. Improved accuracy in the upper regions of an unconfined aquifer is demonstrated for pumping from a partially penetrating well. The change introduces a nonlinearity into the solution of the flow equations that results in slightly slower convergence of the flow solution, 7% slower in the presented demonstration. Accuracy of simulations is improved where vertical flow is dominated by a moving water table, but not when a large water table gradient dominates over the water table movement.     

Grid cell distortion and MODFLOW's integrated finite-difference numerical solution

The ground water flow model MODFLOW inherently implements a nongeneralized integrated finite-difference (IFD) numerical scheme. The IFD numerical scheme allows for construction of finite-difference model grids with curvilinear (piecewise linear) rows. The resulting grid comprises model cells in the shape of trapezoids and is distorted in comparison to a traditional MODFLOW finite-difference grid. A version of MODFLOW-88 (herein referred to as MODFLOW IFD) with the code adapted to make the one-dimensional DELR and DELC arrays two dimensional, so that equivalent conductance between distorted grid cells can be calculated, is described. MODFLOW IFD is used to inspect the sensitivity of the numerical head and velocity solutions to the level of distortion in trapezoidal grid cells within a converging radial flow domain. A test problem designed for the analysis implements a grid oriented such that flow is parallel to columns with converging widths. The sensitivity analysis demonstrates MODFLOW IFD's capacity to numerically derive a head solution and resulting intercell volumetric flow when the internal calculation of equivalent conductance accounts for the distortion of the grid cells. The sensitivity of the velocity solution to grid cell distortion indicates criteria for distorted grid design. In the radial flow test problem described, the numerical head solution is not sensitive to grid cell distortion. The accuracy of the velocity solution is sensitive to cell distortion with error <1% if the angle between the nonparallel sides of trapezoidal cells is <12.5 degrees. The error of the velocity solution is related to the degree to which the spatial discretization of a curve is approximated with piecewise linear segments. Curvilinear finite-difference grid construction adds versatility to spatial discretization of the flow domain. MODFLOW-88's inherent IFD numerical scheme and the test problem results imply that more recent versions of MODFLOW 2000, with minor modifications, have the potential to make use of a curvilinear grid.     

Flow‐Through Stream Modeling with MODFLOW and MT3D: Certainties and Limitations

This paper aims to assess MODFLOW and MT3D capabilities for simulating the spread of contaminants from a river exhibiting an unusual relationship with an alluvial aquifer, with the groundwater head higher than the river head on one side and lower on the other (flow‐through stream). A series of simulation tests is conducted using a simple hypothetical model so as to characterize and quantify these limitations. Simulation results show that the expected contaminant spread could be achieved with a specific configuration composed of two sets of parameters: (1) modeled object parameters (hydraulic groundwater gradient, hydraulic conductivity values of aquifer and streambed), and (2) modeling parameters (vertical discretization of aquifer, horizontal refinement of stream modeled with River [RIV] package). The influence of these various parameters on simulation results is investigated, and potential complications and errors are identified. Contaminant spread from stream to aquifer is not always reproduced by MT3D due to the RIV package's inability to simulate lateral exchange fluxes between stream and aquifer. This paper identifies the need for a MODFLOW streamflow package allowing lateral stream‐aquifer interactions and streamflow routine calculations. Such developments could be of particular interest for modeling contaminated flow‐through streams.     

Modeling axisymmetric flow and transport

Unmodified versions of common computer programs such as MODFLOW, MT3DMS, and SEAWAT that use Cartesian geometry can accurately simulate axially symmetric ground water flow and solute transport. Axisymmetric flow and transport are simulated by adjusting several input parameters to account for the increase in flow area with radial distance from the injection or extraction well. Logarithmic weighting of interblock transmissivity, a standard option in MODFLOW, can be used for axisymmetric models to represent the linear change in hydraulic conductance within a single finite-difference cell. Results from three test problems (ground water extraction, an aquifer push-pull test, and upconing of saline water into an extraction well) show good agreement with analytical solutions or with results from other numerical models designed specifically to simulate the axisymmetric geometry. Axisymmetric models are not commonly used but can offer an efficient alternative to full three-dimensional models, provided the assumption of axial symmetry can be justified. For the upconing problem, the axisymmetric model was more than 1000 times faster than an equivalent three-dimensional model. Computational gains with the axisymmetric models may be useful for quickly determining appropriate levels of grid resolution for three-dimensional models and for estimating aquifer parameters from field tests.     

Accelerating Groundwater Flow Simulation in MODFLOW Using JASMIN‐Based Parallel Computing

To accelerate the groundwater flow simulation process, this paper reports our work on developing an efficient parallel simulator through rebuilding the well‐known software MODFLOW on JASMIN (J Adaptive Structured Meshes applications Infrastructure). The rebuilding process is achieved by designing patch‐based data structure and parallel algorithms as well as adding slight modifications to the compute flow and subroutines in MODFLOW. Both the memory requirements and computing efforts are distributed among all processors; and to reduce communication cost, data transfers are batched and conveniently handled by adding ghost nodes to each patch. To further improve performance, constant‐head/inactive cells are tagged and neglected during the linear solving process and an efficient load balancing strategy is presented. The accuracy and efficiency are demonstrated through modeling three scenarios: The first application is a field flow problem located at Yanming Lake in China to help design reasonable quantity of groundwater exploitation. Desirable numerical accuracy and significant performance enhancement are obtained. Typically, the tagged program with load balancing strategy running on 40 cores is six times faster than the fastest MICCG‐based MODFLOW program. The second test is simulating flow in a highly heterogeneous aquifer. The AMG‐based JASMIN program running on 40 cores is nine times faster than the GMG‐based MODFLOW program. The third test is a simplified transient flow problem with the order of tens of millions of cells to examine the scalability. Compared to 32 cores, parallel efficiency of 77 and 68% are obtained on 512 and 1024 cores, respectively, which indicates impressive scalability.     

Air and water flows induced by pumping tests in unconfined aquifers with low‐permeability zones

Air flows from the atmosphere into an unconfined aquifer when the water table falls during pumping tests. Pumping test results in unconfined aquifers may be significantly affected by low‐permeability zones (LPZs) near the initial water table position, because they restrict the downward movement of air. A transient, three‐dimensional air–water two‐phase flow model is employed to investigate numerically the effects of local heterogeneity on pumping test results in unconfined aquifers. Two cases of local heterogeneities are considered herein: a LPZ around the pumping well and on one side of the pumping well. Results show that the drawdown with the LPZ is significantly greater than that of the homogeneous aquifer. The differences in drawdown are the most significant at intermediate times and gradually diminish at later times. The LPZ significantly reduces air flow from the atmosphere to the aquifer. The pore air velocity in the LPZ is very low. The air pressure at the observation point under the LPZ when air begins to enter is significantly lower than the air pressure of the homogeneous aquifer at the same point. After that, the air pressure increases quickly and then increases slowly. The time for the air pressure to reach the atmospheric pressure is significantly longer. Copyright © 2013 John Wiley & Sons, Ltd.     

Rainfall infiltration‐derived submarine groundwater discharge from multi‐layered aquifer system terminating at the coastline

This article investigates the quantity of submarine groundwater discharge (SGD) from a coastal multi‐layered aquifer system in response to constant rainfall infiltration. The system comprises an unconfined aquifer, a leaky confined aquifer and an aquitard between them and terminates at the coastline. An approximate analytical solution is derived based on the following assumptions: (i) flow is horizontal in the aquifers and vertical in the aquitard, and (ii) flow in the unconfined aquifer is described by nonlinear Boussinesq equation. The analytical solution is compared with numerical solutions of the strictly two‐dimensional nonlinear model to validate the model assumptions used for the analytical solution. The SGD from the leaky confined aquifer increases with the inland rainfall infiltration recharge and the specific leakage of aquitard. The maximum SGD ranges from 1·87 to 10·37 m³ per day per meter of shoreline when rainfall infiltration ranges from 18·2 to 182 mm/year and the specific leakage of aquitard varies from 10^?9 to 10^?1 l/day. Copyright © 2011 John Wiley & Sons, Ltd.     

A general analytical solution for flow to a single horizontal well by Fourier and Laplace transforms

The objective of this paper is to present an analytical solution for describing the head distribution in an unconfined aquifer with a single pumping horizontal well parallel to a fully penetrating stream. The Laplace-domain solution is developed by applying Fourier sine, Fourier and Laplace transforms to the governing equation as well as the associated initial and boundary conditions. The time-domain solution is obtained after taking the inverse Laplace transform along with the Bromwich integral method and inverse Fourier and Fourier sine transforms. The upper boundary condition of the aquifer is represented by the free surface equation in which the second-order slope terms are neglected. Based on the solution and Darcy’s law, the equation representing the stream depletion rate is then derived. The solution can simulate head distributions in an aquifer infinitely extending in horizontal direction if the well is located far away from the stream. In addition, the solution can also simulate head distributions in confined aquifers if specific yield is set zero. It is shown that the solution can be applied practically to evaluate flow to a horizontal well.     

Predicting Collector Well Yields with MODFLOW

Groundwater flow models are commonly used to design new wells and wellfields. As the spatial scale of the problem is large and much local‐scale detail is not needed, modelers often utilize two‐dimensional (2D) or quasi three‐dimensional models based on the Dupuit‐Forchheimer assumption. Dupuit models offer a robust set of tools for simulating regional groundwater flow including interactions with surface waters, the potential for well interference, and varying aquifer properties and recharge rates. However, given an assumed operating water level or drawdown at a well screen, Dupuit models systematically overpredict well yields. For design purposes, this discrepancy is unacceptable, and a method for predicting accurate well yields is needed. While published methods exist for vertical wells, little guidance is available for predicting yields in horizontal screens or collector wells. In plan view, a horizontal screen has a linear geometry, and will likely extend over several neighboring cells that may not align with rows or columns in a numerical model. Furthermore, the model must account for the effects of converging three‐dimensional (3D) flow to the well screens and hydraulic interference among the well screens; these all depend on the design of a specific well. This paper presents a new method for simulating the yield of angled or horizontal well screens in numerical groundwater flow models, specifically using the USGS code MODFLOW. The new method is compared to a detailed, 3D analytic element model of a collector well in a field of uniform flow.