Modeling multiaquifer wells with MODFLOW

Multiaquifer wells, i.e., wells that are open across more than one aquifer, can have a profound effect on the hydraulics of a ground water system. These wells change the physical system by establishing direct hydraulic links between otherwise isolated strata. Several methods are available to simulate multiaquifer wells in the context of comprehensive ground water flow simulators. In this paper, we review four methods to represent multiaquifer wells with the widely used code MODFLOW. These methods include a specialized code developed, but never formally released, by the U.S. Geological Survey (USGS), the Multi-Aquifer Well (MAW1) Package. An expanded implementation of the techniques in the MAW1 Package has been incorporated in the Multi-Node Well Package released recently by the USGS (Halford and Hanson 2002). We examine the performance of the methods in the context of a benchmarking study against the analytical solutions of Papadopulos (1966) and Sokol (1963). Our results demonstrate that results obtained with the MAW1 Package closely match exact solutions for pumping and nonpumping conditions, using both coarse and refined grids.     

Modeling Surface Water-Groundwater Interaction with MODFLOW: Some Considerations

The accuracy with which MODFLOW simulates surface water-groundwater interaction is examined for connected and disconnected losing streams. We compare the effect of different vertical and horizontal discretization within MODFLOW and also compare MODFLOW simulations with those produced by HydroGeoSphere. HydroGeoSphere is able to simulate both saturated and unsaturated flow, as well as surface water, groundwater and the full coupling between them in a physical way, and so is used as a reference code to quantify the influence of some of the simplifying assumptions of MODFLOW. In particular, we show that (1) the inability to simulate negative pressures beneath disconnected streams in MODFLOW results in an underestimation of the infiltration flux; (2) a river in MODFLOW is either fully connected or fully disconnected, while in reality transitional stages between the two flow regimes exist; (3) limitations in the horizontal discretization of the river can cause a mismatch between river width and cell width, resulting in an error in the water table position under the river; and (4) because coarse vertical discretization of the aquifer is often used to avoid the drying out of cells, this may result in an error in simulating the height of the groundwater mound. Conditions under which these errors are significant are investigated.     

Numerical Modeling of Slug Tests with MODFLOW Using Equivalent Well Blocks

Slug tests are a widely used technique to estimate aquifer hydraulic parameters and the test data are generally interpreted with analytical solutions under various assumptions. However, these solutions are not convenient when slug tests are required to be analyzed in a three‐dimensional model for complex aquifer‐aquitard systems. In this study, equivalent well blocks (EWB) are proposed in numerical modeling of slug test data with MODFLOW. Multi‐well slug tests in partially penetrating wells with skin zones can be simulated. Accuracy of the numerical method is demonstrated by benchmarking with analytical solutions. The EWB method is applied in a case study on slug tests in aquitards in the Pearl River Delta, China.     

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 surface and ground water mixing in the hyporheic zone using MODFLOW and MT3D

We used a three-dimensional MODFLOW model, paired with MT3D, to simulate hyporheic zones around debris dams and meanders along a semi-arid stream. MT3D simulates both advective transport and sink/source mixing of solutes, in contrast to particle tracking (e.g. MODPATH), which only considers advection. We delineated the hydrochemically active hyporheic zone based on a new definition, specifically as near-stream subsurface zones receiving a minimum of 10% surface water within a 10-day travel time. Modeling results indicate that movement of surface water into the hyporheic zone is predominantly an advective process. We show that debris dams are a key driver of surface water into the subsurface along the experimental reach, causing the largest flux rates of water across the streambed and creating hyporheic zones with up to twice the cross-sectional area of other hyporheic zones. Hyporheic exchange was also found in highly sinuous segments of the experimental reach, but flux rates are lower and the cross-sectional areas of these zones are generally smaller. Our modeling approach simulated surface and ground water mixing in the hyporheic zone, and thus provides numerical approximations that are more comparable to field-based observations of surface–groundwater exchange than standard particle-tracking simulations.     

Spatial Periodic Boundary Condition for MODFLOW

Small‐scale hyporheic zone (HZ) models often use a spatial periodic boundary (SPB) pair to simulate an infinite repetition of bedforms. SPB's are common features of commercially available multiphysics modeling packages. MODFLOW's lack of this boundary type has precluded it from being effectively utilized in this area of HZ research. We present a method to implement the SPB in MODFLOW by development of the appropriate block‐centered finite‐difference expressions. The implementation is analogous to MODFLOW's general head boundary package. The difference is that the terms on the right hand side of the solution equations must be updated with each iteration. Consequently, models that implement the SPB converge best with solvers that perform both inner and outer iterations. The correct functioning of the SPB condition in MODFLOW is verified by two examples. This boundary condition allows users to build HZ‐bedform models in MODFLOW, facilitating further research using related codes such as MT3DMS and PHT3D.     

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.     

Efficient Storage of Large MODFLOW Models

We present a methodology for storing the bulkier portions of a set of MODFLOW input and output files in a compressed binary format using the HDF5 library. This approach results in compression ratios of up to 99% with no significant time penalty. The highly compressed format is particularly beneficial when dealing with large regional models or Monte Carlo simulations. The strategy is focused on the list‐ and array‐based portions of the input files including the cell property and recharge arrays, and is compatible with models containing parameters, including pilot points. The utilities are based on a modified version of the MODFLOW code and are, therefore, compatible with any standard MODFLOW simulation. We present used cases and instructions on how to use the utilities.     

A Parallel PCG Solver for MODFLOW

In order to simulate large-scale ground water flow problems more efficiently with MODFLOW, the OpenMP programming paradigm was used to parallelize the preconditioned conjugate-gradient (PCG) solver with in this study. Incremental parallelization, the significant advantage supported by OpenMP on a shared-memory computer, made the solver transit to a parallel program smoothly one block of code at a time. The parallel PCG solver, suitable for both MODFLOW-2000 and MODFLOW-2005, is verified using an 8-processor computer. Both the impact of compilers and different model domain sizes were considered in the numerical experiments. Based on the timing results, execution times using the parallel PCG solver are typically about 1.40 to 5.31 times faster than those using the serial one. In addition, the simulation results are the exact same as the original PCG solver, because the majority of serial codes were not changed. It is worth noting that this parallelizing approach reduces cost in terms of software maintenance because only a single source PCG solver code needs to be maintained in the MODFLOW source tree.     

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.     

MODFLOW as a Configurable Multi-Model Hydrologic Simulator

MODFLOW 6 is the latest in a line of six “core” versions of MODFLOW released by the U.S. Geological Survey. The MODFLOW 6 architecture supports incorporation of additional hydrologic processes, in addition to groundwater flow, and allows interaction between processes. The architecture supports multiple model instances and multiple types of models within a single simulation, a flexible approach to formulating and solving the equations that represent hydrologic processes, and recent advances in interoperability, which allow MODFLOW to be accessed and controlled by external programs. The present version of MODFLOW 6 consolidates popular capabilities available in MODFLOW variants, such as the unstructured grid support in MODFLOW-USG, the Newton-Raphson formulation in MODFLOW-NWT, and the support for partitioned stress boundaries in MODFLOW-CDSS. The flexible multi-model capability allows users to configure MODFLOW 6 simulations to represent the local-grid refinement (LGR) capabilities available in MODFLOW-LGR, the multi-species transport capabilities in MT3DMS, and the coupled variable-density capabilities available in SEAWAT. This paper provides a new, holistic and integrated overview of simulation capabilities made possible by the MODFLOW 6 architecture, and describes how ongoing and future development can take advantage of the program architecture to integrate new capabilities in a way that is minimally invasive and automatically compatible with the existing MODFLOW 6 code.     

Parallelization of MODFLOW Using a GPU Library

A new method based on a graphics processing unit (GPU) library is proposed in the paper to parallelize MODFLOW. Two programs, GetAb_CG and CG_GPU, have been developed to reorganize the equations in MODFLOW and solve them with the GPU library. Experimental tests using the NVIDIA Tesla C1060 show that a 1.6‐ to 10.6‐fold speedup can be achieved for models with more than 10⁵ cells. The efficiency can be further improved by using up‐to‐date GPU devices.     

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.