A Hybrid Finite-Difference and Analytic Element Groundwater Model

Regional finite-difference models tend to have large cell sizes, often on the order of 1–2 km on a side. Although the regional flow patterns in deeper formations may be adequately represented by such a model, the intricate surface water and groundwater interactions in the shallower layers are not. Several stream reaches and nearby wells may occur in a single cell, precluding any meaningful modeling of the surface water and groundwater interactions between the individual features. We propose to replace the upper MODFLOW layer or layers, in which the surface water and groundwater interactions occur, by an analytic element model (GFLOW) that does not employ a model grid; instead, it represents wells and surface waters directly by the use of point-sinks and line-sinks. For many practical cases it suffices to provide GFLOW with the vertical leakage rates calculated in the original coarse MODFLOW model in order to obtain a good representation of surface water and groundwater interactions. However, when the combined transmissivities in the deeper (MODFLOW) layers dominate, the accuracy of the GFLOW solution diminishes. For those cases, an iterative coupling procedure, whereby the leakages between the GFLOW and MODFLOW model are updated, appreciably improves the overall solution, albeit at considerable computational cost. The coupled GFLOW–MODFLOW model is applicable to relatively large areas, in many cases to the entire model domain, thus forming an attractive alternative to local grid refinement or inset models.     

A Semi‐Structured MODFLOW‐USG Model to Evaluate Local Water Sources to Wells for Decision Support

In order to better represent the configuration of the stream network and simulate local groundwater‐surface water interactions, a version of MODFLOW with refined spacing in the topmost layer was applied to a Lake Michigan Basin (LMB) regional groundwater‐flow model developed by the U.S. Geological. Regional MODFLOW models commonly use coarse grids over large areas; this coarse spacing precludes model application to local management issues (e.g., surface‐water depletion by wells) without recourse to labor‐intensive inset models. Implementation of an unstructured formulation within the MODFLOW framework (MODFLOW‐USG) allows application of regional models to address local problems. A “semi‐structured” approach (uniform lateral spacing within layers, different lateral spacing among layers) was tested using the LMB regional model. The parent 20‐layer model with uniform 5000‐foot (1524‐m) lateral spacing was converted to 4 layers with 500‐foot (152‐m) spacing in the top glacial (Quaternary) layer, where surface water features are located, overlying coarser resolution layers representing deeper deposits. This semi‐structured version of the LMB model reproduces regional flow conditions, whereas the finer resolution in the top layer improves the accuracy of the simulated response of surface water to shallow wells. One application of the semi‐structured LMB model is to provide statistical measures of the correlation between modeled inputs and the simulated amount of water that wells derive from local surface water. The relations identified in this paper serve as the basis for metamodels to predict (with uncertainty) surface‐water depletion in response to shallow pumping within and potentially beyond the modeled area, see Fienen et al. (2015a).     

Generating MODFLOW grids from boundary representation solid models

Complex stratigraphy can be difficult to simulate in MODFLOW models. MODFLOW uses a structured grid that requires that each grid layer be continuous throughout the model domain. This makes it difficult to explicitly represent common features such as pinchouts and embedded seams in a MODFLOW model. In this paper, we describe a method for automatically generating MODFLOW-compatible grids from boundary-representation solid models. Solid models are data structures developed originally for computer-aided design applications that define the geometry of three-dimensional objects. Solid models can be used to represent arbitrarily complex stratigraphy. The elevations defined by the solids are then extracted from the solids in a manner that preserves the continuous-layer requirement imposed by MODFLOW. Two basic approaches are described: The first method adjusts the MODFLOW grid dimensions (layer elevations) to fit the solid model boundaries, and the second method creates a regular MODFLOW grid and adjusts the material properties to match the changes in stratigraphy. One of the main benefits of using solid models to define stratigraphy for MODFLOW models is that it provides a grid-independent definition of the layer elevations that can be used to immediately re-create the MODFLOW grid geometry after any change to the grid resolution.     

Simulating ground water-lake interactions: approaches and insights

Approaches for modeling lake-ground water interactions have evolved significantly from early simulations that used fixed lake stages specified as constant head to sophisticated LAK packages for MODFLOW. Although model input can be complex, the LAK package capabilities and output are superior to methods that rely on a fixed lake stage and compare well to other simple methods where lake stage can be calculated. Regardless of the approach, guidelines presented here for model grid size, location of three-dimensional flow, and extent of vertical capture can facilitate the construction of appropriately detailed models that simulate important lake-ground water interactions without adding unnecessary complexity. In addition to MODFLOW approaches, lake simulation has been formulated in terms of analytic elements. The analytic element lake package had acceptable agreement with a published LAKI problem, even though there were differences in the total lake conductance and number of layers used in the two models. The grid size used in the original LAKI problem, however, violated a grid size guideline presented in this paper. Grid sensitivity analyses demonstrated that an appreciable discrepancy in the distribution of stream and lake flux was related to the large grid size used in the original LAKI problem. This artifact is expected regardless of MODFLOW LAK package used. When the grid size was reduced, a finite-difference formulation approached the analytic element results. These insights and guidelines can help ensure that the proper lake simulation tool is being selected and applied.     

A New Package in MODFLOW to Simulate Unconfined Groundwater Flow in Sloping Aquifers

The nonhorizontal‐model‐layer (NHML) grid system is more accurate than the horizontal‐model‐layer grid system to describe groundwater flow in an unconfined sloping aquifer on the basis of MODFLOW‐2000. However, the finite‐difference scheme of NHML was based on the Dupuit‐Forchheimer assumption that the streamlines were horizontal, which was acceptable for slope less than 0.10. In this study, we presented a new finite‐difference scheme of NHML based on the Boussinesq assumption and developed a new package SLOPE which was incorporated into MODFLOW‐2000 to become the MODFLOW‐SP model. The accuracy of MODFLOW‐SP was tested against solution of Mac Cormack (1969). The differences between the solutions of MODFLOW‐2000 and MODFLOW‐SP were nearly negligible when the slope was less than 0.27, and they were noticeable during the transient flow stage and vanished in steady state when the slope increased above 0.27. We established a model considering the vertical flow using COMSOL Multiphysics to test the robustness of constrains used in MODFLOW‐SP. The results showed that streamlines quickly became parallel with the aquifer base except in the narrow regions near the boundaries when the initial flow was not parallel to the aquifer base. MODFLOW‐SP can be used to predict the hydraulic head of an unconfined aquifer along the profile perpendicular to the aquifer base when the slope was smaller than 0.50. The errors associated with constrains used in MODFLOW‐SP were small but noticeable when the slope increased to 0.75, and became significant for the slope of 1.0.     

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.     

Effect of elevation resolution on evapotranspiration simulations using MODFLOW

Surface elevations represented in MODFLOW head-dependent packages are usually derived from digital elevation models (DEMs) that are available at much high resolution. Conventional grid refinement techniques to simulate the model at DEM resolution increases computational time, input file size, and in many cases are not feasible for regional applications. This research aims at utilizing the increasingly available high resolution DEMs for effective simulation of evapotranspiration (ET) in MODFLOW as an alternative to grid refinement techniques. The source code of the evapotranspiration package is modified by considering for a fixed MODFLOW grid resolution and for different DEM resolutions, the effect of variability in elevation data on ET estimates. Piezometric head at each DEM cell location is corrected by considering the gradient along row and column directions. Applicability of the research is tested for the lower Rio Grande (LRG) Basin in southern New Mexico. The DEM at 10 m resolution is aggregated to resampled DEM grid resolutions which are integer multiples of MODFLOW grid resolution. Cumulative outflows and ET rates are compared at different coarse resolution grids. Results of the analysis conclude that variability in depth-to-groundwater within the MODFLOW cell is a major contributing parameter to ET outflows in shallow groundwater regions. DEM aggregation methods for the LRG Basin have resulted in decreased volumetric outflow due to the formation of a smoothing error, which lowered the position of water table to a level below the extinction depth.     

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.     

Studying the Flow Dynamics of a Karst Aquifer System with an Equivalent Porous Medium Model

The modeling of groundwater flow in karst aquifers is a challenge due to the extreme heterogeneity of its hydraulic parameters and the duality in their discharge behavior, that is, rapid response of highly conductive karst conduits and delayed drainage of the low‐permeability fractured matrix after recharge events. There are a number of different modeling approaches for the simulation of the karst groundwater dynamics, applicable to different aquifer as well as modeling problem types, ranging from continuum models to double continuum models to discrete and hybrid models. This study presents the application of an equivalent porous model approach (EPM, single continuum model) to construct a steady‐state numerical flow model for an important karst aquifer, that is, the Western Mountain Aquifer Basin (WMAB), shared by Israel and the West‐Bank, using MODFLOW2000. The WMAB was used as a catchment since it is a well‐constrained catchment with well‐defined recharge and discharge components and therefore allows a control on the modeling approach, a very rare opportunity for karst aquifer modeling. The model demonstrates the applicability of equivalent porous medium models for the simulation of karst systems, despite their large contrast in hydraulic conductivities. As long as the simulated saturated volume is large enough to average out the local influence of karst conduits and as long as transport velocities are not an issue, EPM models excellently simulate the observed head distribution. The model serves as a starting basis that will be used as a reference for developing a long‐term dynamic model for the WMAB, starting from the pre‐development period (i.e., 1940s) up to date.     

River Seepage Conductance in Large‐Scale Regional Studies

Flow exchange between surface and groundwater is of great importance be it for beneficial allocation and use of water resources or for the proper exercise of water rights. In large‐scale regional studies, most numerical models use coarse grid sizes, which make it difficult to provide an accurate depiction of the phenomenon. In particular, a somewhat arbitrary leakance coefficient in a third type (i.e., Cauchy, General Head) boundary condition is used to calculate the seepage discharge as a function of the difference of head in the river and in the aquifer, whose value is often found by calibration. A different approach is presented to analytically estimate that leakance coefficient. It is shown that a simple equivalence can be deduced from the analytical solution for the empirical coefficient, so that it provides the accuracy of the analytical solution while the model maintains a very coarse grid, treating the water‐table aquifer as a single calculation layer. Relating the empirical leakance coefficient to the exact conductance, derived from physical principles, provides a physical basis for the leakance coefficient. Factors such as normalized wetted perimeter, degree of penetration of the river, presence of a clogging layer, and anisotropy can be included with little computational demand. In addition the river coefficient in models such as MODFLOW, for example, can be easily modified when grid size is changed without need for recalibration.     

A Spatially Periodic Solute Boundary for MT3DMS and PHT3D

The assumption of spatial repetition is commonly made when producing bedform scale models of the hyporheic zone. Two popular solute transport codes, MT3DMS and PHT3D, do not currently provide the necessary boundary condition required to simulate spatial periodicity in hyporheic zone transport problems. In this study, we develop a spatially periodic boundary (SPB) for solutes that is compatible with a SPB that was previously developed for MODFLOW to simulate the flow component of spatially periodic problems. The approach is ideal for simulating groundwater flow and transport patterns under repeating surface features, such as ripples or dunes on the bottom of a lake or stream. The appropriate block‐centered finite‐difference approach to implement the boundary is presented and the necessary source code modifications are discussed. The performance of the solute SPB, operating in conjunction with the groundwater flow SPB, is explored through comparison of a multi‐bedform hyporheic‐zone model with a single bedform variant. The new boundary conditions perform well in situations where both dispersive effects and lateral seepage flux in the underflow regime beneath the hyporheic zone are minimal.     

Unstructured Gridding for MODFLOW from Prior Groundwater Flow Models: A New Paradigm

This work introduces a new unstructured gridding approach relying on feedback from a previous groundwater flow model. All cells in a relatively coarse model using a rectilinear grid are recursively subdivided following a cell wise specific discharge-based indicator to generate quadtree, octree or Voronoï grids. This technique leverages the full potential of the latest MODFLOW engines. The suitability of this approach is demonstrated on challenging single and multilayered heterogeneous formations. The proposed method is straightforward to implement in existing software packages. It supports iterative updating of groundwater flow models from the legacy rectilinear to unstructured grids.     

Application of a perfectly matched layer in seismic wavefield simulation with an irregular free surface

Recently, an effective and powerful approach for simulating seismic wave propagation in elastic media with an irregular free surface was proposed. However, in previous studies, researchers used the periodic condition and/or sponge boundary condition to attenuate artificial reflections at boundaries of a computational domain. As demonstrated in many literatures, either the periodic condition or sponge boundary condition is simple but much less effective than the well‐known perfectly matched layer boundary condition. In view of this, we intend to introduce a perfectly matched layer to simulate seismic wavefields in unbounded models with an irregular free surface. We first incorporate a perfectly matched layer into wave equations formulated in a frequency domain in Cartesian coordinates. We then transform them back into a time domain through inverse Fourier transformation. Afterwards, we use a boundary‐conforming grid and map a rectangular grid onto a curved one, which allows us to transform the equations and free surface boundary conditions from Cartesian coordinates to curvilinear coordinates. As numerical examples show, if free surface boundary conditions are imposed at the top border of a model, then it should also be incorporated into the perfectly matched layer imposed at the top‐left and top‐ right corners of a 2D model where the free surface boundary conditions and perfectly matched layer encounter; otherwise, reflections will occur at the intersections of the free surface and the perfectly matched layer, which is confirmed in this paper. So, by replacing normal second derivatives in wave equations in curvilinear coordinates with free surface boundary conditions, we successfully implement the free surface boundary conditions into the perfectly matched layer at the top‐left and top‐right corners of a 2D model at the surface. A number of numerical examples show that the perfectly matched layer constructed in this study is effective in simulating wave propagation in unbounded media and the algorithm for implementation of the perfectly matched layer and free surface boundary conditions is stable for long‐time wavefield simulation on models with an irregular free surface.     

Simulation of Lake-Groundwater Interaction under Steady-State Flow

MODFLOW is one of the most popular groundwater simulation tools available; however, the development of lake modules that can be coupled with MODFLOW is lacking apart from the LAK3 package. This study proposes a new approach for simulating lake - groundwater interaction under steady-state flow, referred to as the sloping lakebed method (SLM). In this new approach, discretization of the lakebed in the vertical direction is independent of the spatial discretization of the aquifer system, which can potentially solve the problem that the lake and groundwater are usually simulated at different scales. The lakebed is generalized by a slant at the bottom of each lake grid cell, which can be classified as fully submerged, dry, and partly submerged. The SLM method accounts for all lake sources and sinks, establishing a governing equation that can be solved using Newton's method. A benchmarking case study was conducted using a modified model setup in the LAK3 user manual. It was found that when there is a sufficient number of layers at the top of the groundwater model, SLM simulates an almost identical groundwater head as the LAK3-based model; when the number of layers decreases, SLM is unaffected while LAK3 may be at a risk of giving unrealistic results. Additionally, the SLM can reflect the relationship between the simulated lake surface area and lake water depth more accurately. Therefore, the SLM method is a promising alternative to the LAK3 package when simulating lake - groundwater interaction.     

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.     

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.     

Vertical Discretization Impact in Numerical Modeling of Matrix Diffusion in Contaminated Groundwater

Understanding the effects of contaminants that can diffuse into low-permeability (“low-k”) zones is crucial for effective groundwater remedial decision-making. Because low-k zones can serve as low-level sources of contamination to more transmissive zones over time, an accurate evaluation of the impacts of matrix diffusion at contaminated sites is vital. This study compared numerical groundwater flow and transport simulations using MODFLOW/RT3D at a hypothetical site using three cases, each with increasing discretization of the vertical 10-m thick domain: (1) a coarse multilayer heterogeneous grid based on one layer for each of four different hydrogeological units, (2) a “low-resolution” discretization approach where the low-k units were divided into several sublayers giving the model 10 layers, and (3) a “high-resolution” numerical model with 199 layers that are a few centimeters thick. When comparing the results of each case, significant differences were observed between the discretizations used, even though all other model input data were identical. The conventional grid models (Cases 1 and 2) appeared to underestimate groundwater plume concentrations by a factor ranging from 1.1 to 36 when compared to the high-resolution grid model (Case 3), and underestimated predicted cleanup times by more than a factor of 10 for some of the hypothetical sampling points in the modeling domain. These results validate the implication of Chapman et al. (2012), that conventional vertical discretization of numerical groundwater flow and transport models at contaminated sites (with layers that are greater than 1 m thick) can lead to significant errors when compared to more accurate high-resolution vertical discretization schemes (layers that are centimeters thick).     

Analytic Element Modeling of Steady Interface Flow in Multilayer Aquifers Using AnAqSim

This paper presents the analytic element modeling approach implemented in the software AnAqSim for simulating steady groundwater flow with a sharp fresh‐salt interface in multilayer (three‐dimensional) aquifer systems. Compared with numerical methods for variable‐density interface modeling, this approach allows quick model construction and can yield useful guidance about the three‐dimensional configuration of an interface even at a large scale. The approach employs subdomains and multiple layers as outlined by Fitts (2010) with the addition of discharge potentials for shallow interface flow (Strack 1989). The following simplifying assumptions are made: steady flow, a sharp interface between fresh‐ and salt water, static salt water, and no resistance to vertical flow and hydrostatic heads within each fresh water layer. A key component of this approach is a transition to a thin fixed minimum fresh water thickness mode when the fresh water thickness approaches zero. This allows the solution to converge and determine the steady interface position without a long transient simulation. The approach is checked against the widely used numerical codes SEAWAT and SWI/MODFLOW and a hypothetical application of the method to a coastal wellfield is presented.     

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.     

Nonlinear Flow Process: A New Package to Compute Nonlinear Flow in MODFLOW

A new MODFLOW package (Nonlinear Flow Process; NLFP) simulating nonlinear flow following the Forchheimer equation was developed and implemented in MODLFOW‐2005. The method is based on an iterative modification of the conductance calculated and used by MODFLOW to obtain an effective Forchheimer conductance. The package is compatible with the different layer types, boundary conditions, and solvers as well as the wetting capability of MODFLOW. The correct implementation is demonstrated using four different benchmark scenarios for which analytical solutions are available. A scenario considering transient flow in a more realistic setting and a larger model domain with a higher number of cells demonstrates that NLFP performs well under more complex conditions, although it converges moderately slower than the standard MODFLOW depending on the nonlinearity of flow. Thus, this new tool opens a field of opportunities to groundwater flow simulation with MODFLOW, especially for core sample simulation or vuggy karstified aquifers as well as for nonlinear flow in the vicinity of pumping wells.