Evaluation of Bias Associated with Capture Maps Derived from Nonlinear Groundwater Flow Models

The impact of groundwater withdrawal on surface water is a concern of water users and water managers, particularly in the arid western United States. Capture maps are useful tools to spatially assess the impact of groundwater pumping on water sources (e.g., streamflow depletion) and are being used more frequently for conjunctive management of surface water and groundwater. Capture maps have been derived using linear groundwater flow models and rely on the principle of superposition to demonstrate the effects of pumping in various locations on resources of interest. However, nonlinear models are often necessary to simulate head‐dependent boundary conditions and unconfined aquifers. Capture maps developed using nonlinear models with the principle of superposition may over‐ or underestimate capture magnitude and spatial extent. This paper presents new methods for generating capture difference maps, which assess spatial effects of model nonlinearity on capture fraction sensitivity to pumping rate, and for calculating the bias associated with capture maps. The sensitivity of capture map bias to selected parameters related to model design and conceptualization for the arid western United States is explored. This study finds that the simulation of stream continuity, pumping rates, stream incision, well proximity to capture sources, aquifer hydraulic conductivity, and groundwater evapotranspiration extinction depth substantially affect capture map bias. Capture difference maps demonstrate that regions with large capture fraction differences are indicative of greater potential capture map bias. Understanding both spatial and temporal bias in capture maps derived from nonlinear groundwater flow models improves their utility and defensibility as conjunctive‐use management tools.     

Estimates of Horizontal Groundwater Flow Velocities in Boreholes

Currently, monitoring tools can be deployed in observation boreholes to better assess groundwater flow, flux of dissolved contaminants and their mass discharge in an aquifer. The relationship between horizontal water velocity in observation boreholes and Darcy fluxes in the surrounding aquifer has been studied for natural flow conditions (i.e., no pumping). Interpretation of measurements taken with dilution tests, the colloidal borescope, the Heat Pulse Flowmeter, and other techniques require the conversion of observed borehole velocity u to aquifer Darcy flux q_∞ . This conversion is typically done through a proportionality factor α = u/q_∞ . In experimental studies as well as in theoretical developments, reported values of α vary almost three orders of magnitude (from 0.5 to 10). This large variability in reported values of α could be explained by: (1) unclear distinction between Darcy flux and water seepage velocity, (2) unclear definition of water velocity in the borehole, (3) effects of well screen and the presence of the measurement device itself on the observable velocities, and (4) hydraulic conditions in the borehole annulus. We address (1), (2) from a conceptual/theoretical perspective, and (3) by means of numerical simulations. We show that issue (1) in low porosity aquifers can yield to order-of-magnitude discrepancies in estimates of q_∞ ; (2) may result in discrepancies of up to 50%, and (3) can cause differences up to 20% of water velocity in the borehole void space compared to the theoretical case of an open borehole.     

Effects of Improper Characterization of Aquifer Thickness on Estimates of Hydraulic Conductivity from Pumping Tests

Pumping test data for surficial aquifers are commonly analyzed under the assumption that the base of the aquifer corresponds to the bottom of the test wells (i.e., the aquifer is truncated). This practice can lead to inaccurate hydraulic conductivity estimates, resulting from the use of low saturated thickness values with transmissivity estimates, and not accounting for the effects of partially penetrating wells. Theoretical time-drawdown data were generated at an observation well in a hypothetical unconfined aquifer for various values of saturated thickness and were analyzed by standard curve-matching techniques. The base of the aquifer was assumed to be the bottom of the pumping and observation wells. The overestimation of horizontal hydraulic conductivity was found to be directly proportional to the error in assumed saturated thickness, and to the (actual) ratio of vertical to horizontal hydraulic conductivity (K_v/K_h). Inaccurately high estimates of hydraulic conductivity obtained by aquifer truncation can lead to overestimates of ground water velocity and contaminant plume spreading, narrow capture zone configuration estimates, and overestimates of available ground water resources.     

On Solving Groundwater Flow and Transport Models with Algebraic Multigrid Preconditioning

Iterative solvers preconditioned with algebraic multigrid have been devised as an optimal technology to speed up the response of large sparse linear systems. In this work, this technique was implemented in the framework of the dual delineation approach. This involves a single groundwater flow linear solution and a pure advective transport solution with different right-hand sides. The new solver was compared with other preconditioned iterative methods, the MODFLOW's GMG solver, and direct sparse solvers. Test problems include two- and three-dimensional benchmarks spanning homogeneous and highly heterogeneous and anisotropic formations. For the groundwater flow problems, using the algebraic multigrid preconditioning speeds up the numerical solution by one to two orders of magnitude. The algebraic multigrid preconditioner efficiency was preserved for the three dimensional heterogeneous and anisotropic problem unlike for the MODFLOW's GMG solver. Contrarily, a sparse direct solver was the most efficient for the pure advective transport processes such as the forward travel time simulations. Hence, the best sparse solver for the more general advection-dispersion transport equation is likely to be Péclet number dependent. When equipped with the best solvers, processing multimillion grid blocks by the dual delineation approach is a matter of seconds. This paves the way for its routine application to large geological models. The paper gives practical hints on the strategies and conditions under which algebraic multigrid preconditioning would remain competitive for the class of nonlinear and/or transient problems.     

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.     

Including Vertical Fault Structures in Layered Groundwater Flow Models

Accurate representation of groundwater flow and solute transport requires a sound representation of the underlying geometry of aquifers. Faults can have a significant influence on the structure and connectivity of aquifers, which may allow permeable units to connect, and aquifers to seal when juxtaposed against lower permeability units. Robust representation of groundwater flow around faults remains challenging despite the significance of faults for flow and transport. We present a methodology for the inclusion of faults utilizing the unstructured grid features of MODFLOW-USG and MODFLOW 6. The method focuses on the representation of fault geometries using non-neighbor connections between juxtaposed layers. We present an illustration of the method for a synthetic fluvial aquifer. The combined impact of the heterogeneous aquifer and fault offset is clearly visible where channel features at different depths in the aquifer were connected at the fault. These results highlight the importance of representing fault features in groundwater flow models.     

A Correction on Coastal Heads for Groundwater Flow Models

We introduce a simple correction to coastal heads for constant‐density groundwater flow models that contain a coastal boundary, based on previous analytical solutions for interface flow. The results demonstrate that accurate discharge to the sea in confined aquifers can be obtained by direct application of Darcy's law (for constant‐density flow) if the coastal heads are corrected to ((α + 1)/α)h_s ? B/2α, in which h_s is the mean sea level above the aquifer base, B is the aquifer thickness, and α is the density factor. For unconfined aquifers, the coastal head should be assigned the value . The accuracy of using these corrections is demonstrated by consistency between constant‐density Darcy's solution and variable‐density flow numerical simulations. The errors introduced by adopting two previous approaches (i.e., no correction and using the equivalent fresh water head at the middle position of the aquifer to represent the hydraulic head at the coastal boundary) are evaluated. Sensitivity analysis shows that errors in discharge to the sea could be larger than 100% for typical coastal aquifer parameter ranges. The location of observation wells relative to the toe is a key factor controlling the estimation error, as it determines the relative aquifer length of constant‐density flow relative to variable‐density flow. The coastal head correction method introduced in this study facilitates the rapid and accurate estimation of the fresh water flux from a given hydraulic head measurement and allows for an improved representation of the coastal boundary condition in regional constant‐density groundwater flow models.     

Review of Modeling Approaches to Groundwater Flow in Deformed Carbonate Aquifers

We discuss techniques to represent groundwater flow in carbonate aquifers using the three existing modeling approaches: equivalent porous medium, conduit network, and discrete fracture network. Fractures in faulted stratigraphic successions are characterized by dominant sets of sub-vertical joints. Grid rotation is recommended using the equivalent porous medium to match higher hydraulic conductivity with the dominant orientation of the joints. Modeling carbonate faults with throws greater than approximately 100 m is more challenging. Such faults are characterized by combined conduit-barrier behavior. The barrier behavior can be modeled using the Horizontal Flow Barrier Package with a low-permeability vertical barrier inserted to represent the impediment of horizontal flow in faults characterized by sharp drops of the piezometric surface. Cavities can occur parallel to the strike of normal faults generating channels for the groundwater. In this case, flow models need to account for turbulence using a conduit network approach. Channels need to be embedded in an equivalent porous medium due to cavities a few centimeters large, which are present in carbonate aquifers even in areas characterized by low hydraulic gradients. Discrete fracture network modeling enables representation of individual rock discontinuities in three dimensions. This approach is used in non-heavily karstified aquifers at industrial sites and was recently combined with the equivalent porous medium to simulate diffusivity in the matrix. Following this review, we recommend that the future research combines three practiced modeling approaches: equivalent porous medium, discrete fracture network, and conduit network, in order to capture structural and flow aspects in the modeling of groundwater in carbonate rocks.     

Bootstrap Calibration and Uncertainty Estimation of Downhole NMR Hydraulic Conductivity Estimates in an Unconsolidated Aquifer

Characterization of hydraulic conductivity (K) in aquifers is critical for evaluation, management, and remediation of groundwater resources. While estimates of K have been traditionally obtained using hydraulic tests over discrete intervals in wells, geophysical measurements are emerging as an alternative way to estimate this parameter. Nuclear magnetic resonance (NMR) logging, a technology once largely applied to characterization of deep consolidated rock petroleum reservoirs, is beginning to see use in near‐surface unconsolidated aquifers. Using a well‐known rock physics relationship—the Schlumberger Doll Research (SDR) equation—K and porosity can be estimated from NMR water content and relaxation time. Calibration of SDR parameters is necessary for this transformation because NMR relaxation properties are, in part, a function of magnetic mineralization and pore space geometry, which are locally variable quantities. Here, we present a statistically based method for calibrating SDR parameters that establishes a range for the estimated parameters and simultaneously estimates the uncertainty of the resulting K values. We used co‐located logging NMR and direct K measurements in an unconsolidated fluvial aquifer in Lawrence, Kansas, USA to demonstrate that K can be estimated using logging NMR to a similar level of uncertainty as with traditional direct hydraulic measurements in unconsolidated sediments under field conditions. Results of this study provide a benchmark for future calibrations of NMR to obtain K in unconsolidated sediments and suggest a method for evaluating uncertainty in both K and SDR parameter values.     

Re-Purposing Groundwater Flow Models for Age Assessments: Important Characteristics

Groundwater flow model construction is often time-consuming and costly, with development ideally focused on a specific purpose, such as quantifying well capture from water bodies or providing flow fields for simulating advective transport. As environmental challenges evolve, the incentive to re-purpose existing groundwater flow models may increase. However, few studies have evaluated which characteristics of groundwater flow models deserve greatest consideration when re-purposing models for groundwater age and advective transport simulations. In this paper, we compare simulated age metrics produced by three MODFLOW-MODPATH models of the same area but with differing levels of complexity (layering and heterogeneity). Comparisons are made at three watershed scales (HUC 8 to HUC 12). Groundwater age metrics, specifically the young fraction and median age of the young and old fractions, are used for evaluation because they relate to intrinsic susceptibility of aquifers and are simpler to interpret than full age distributions used for advective transport. Results indicate that: (1) the young fraction is less sensitive to model layering than the median age of young and old fractions, suggesting that simple models may suffice for basic intrinsic susceptibility assessments; (2) water table mounding and associated discharge into partially penetrating boundaries, such as head-water streams, is important for simulating both the young fraction and the median age of the young fraction; and (3) the influence of partially penetrating head-water streams is maintained regardless of the porosity distribution. Results of this work should aid modelers with evaluating the appropriateness of re-purposing existing groundwater flow models for age simulations.     

Multimodel Validity Assessment of Groundwater Flow Simulation Models Using Area Metric Approach

We demonstrate the application of the Area Metric developed by Ferson et al. (2008) for multimodel validity assessment. The Area Metric quantified the degree of models' replicative validity: the degree of agreement between the observed data and the corresponding simulated outputs represented as their empirical cumulative distribution functions (ECDFs). This approach was used to rank multiple representations of a case study groundwater flow model of a landfill by their Area Metric scores. A multimodel approach allows to account for uncertainties that may either be epistemic (from lack of knowledge) or aleatory (from variability inherent in the system). The Area Metric approach enables explicit incorporation of model uncertainties, epistemic as well as aleatory, into validation assessment. The proposed approach informs understanding of the collected data and that of the model domain. It avoids model overfitting to a particular system state, and in fact is a blind assessment of the models' validity: models are not adjusted, or updated, to achieve a better numerical fit. This approach assesses the degree of models' validity, in place of the typical binary model validation/invalidation process. Collectively, this increases confidence in the model's representativeness that, in turn, reduces risk to model users.