A conceptual–numerical model to simulate hydraulic head in aquifers that are hydraulically connected to surface water bodies

In this paper, we present a conceptual‐numerical model that can be deduced from a calibrated finite difference groundwater‐flow model, which provides a parsimonious approach to simulate and analyze hydraulic heads and surface water body–aquifer interaction for linear aquifers (linear response of head to stresses). The solution of linear groundwater‐flow problems using eigenvalue techniques can be formulated with a simple explicit state equation whose structure shows that the surface water body–aquifer interaction phenomenon can be approached as the drainage of a number of independent linear reservoirs. The hydraulic head field could be also approached by the summation of the head fields, estimated for those reservoirs, defined over the same domain set by the aquifer limits, where the hydraulic head field in each reservoir is proportional to a specific surface (an eigenfunction of an eigenproblem, or an eigenvector in discrete cases). All the parameters and initial conditions of each linear reservoir can be mathematically defined in a univocal way from the calibrated finite difference model, preserving its characteristics (geometry, boundary conditions, hydrodynamic parameters (heterogeneity), and spatial distribution of the stresses). We also demonstrated that, in practical cases, an accurate solution can be obtained with a reduced number of linear reservoirs. The reduced computational cost of these solutions can help to integrate the groundwater component within conjunctive use management models. Conceptual approximation also facilitates understanding of the physical phenomenon and analysis of the factors that influence it. A simple synthetic aquifer has been employed to show how the conceptual model can be built for different spatial discretizations, the parameters required, and their influence on the simulation of hydraulic head fields and stream–aquifer flow exchange variables. A real‐world case was also solved to test the accuracy of the proposed approaches, by comparing its solution with that obtained using finite‐difference MODFLOW code. Copyright © 2011 John Wiley & Sons, Ltd.     

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.     

Hydraulic-Head Formulation for Density-Dependent Flow and Transport

Density-dependent flow and transport solutions for coastal saltwater intrusion investigations, analyses of fluid injection into deep brines, and studies of convective fingering and instabilities of denser fluids moving through less dense fluids typically formulate the groundwater flow equation in terms of pressure or equivalent freshwater head. A formulation of the flow equation in terms of hydraulic head is presented here as an alternative. The hydraulic-head formulation can facilitate adaptation of existing constant-density groundwater flow codes to include density-driven flow by avoiding the need to convert between freshwater head and hydraulic head within the code and by incorporating density-dependent terms as a compartmentalized “correction” to constant-density calculations already performed by the code. The hydraulic-head formulation also accommodates complexities such as unconfined groundwater flow and Newton-Raphson solution schemes more readily than the freshwater-head formulation. Simulation results are presented for four example problems solved using an implementation of the hydraulic-head formulation in MODFLOW.     

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.     

Transient Solutions to Groundwater Mounding in Bounded and Unbounded Aquifers

In this study, the well‐known Hantush solution procedure for groundwater mounding under infinitely long infiltration strips is extended to finite and semi‐infinite aquifer cases. Initially, the solution for infinite aquifers is presented and compared to those available in literature and to the numerical results of MODFLOW. For the finite aquifer case, the method of images, which is commonly used in well hydraulics, is used to be able to represent the constant‐head boundaries at both sides. It is shown that a finite number of images is enough to obtain the results and sustain the steady state. The effect of parameters on the growth of the mound and on the time required to reach the steady state is investigated. The semi‐infinite aquifer case is emphasized because the growth of the mound is not symmetric. As the constant‐head boundary limits the growth, the unbounded side grows continuously. For this reason, the groundwater divide shifts toward the unbounded side. An iterative solution procedure is proposed. To perform the necessary computations a code was written in Visual Basic of which the algorithm is presented. The proposed methodology has a wide range of applicability and this is demonstrated using two practical examples. The first one is mounding under a stormwater dispersion trench in an infinite aquifer and the other is infiltration from a flood control channel into a semi‐infinite aquifer. Results fit very well with those of MODFLOW.     

Modifications to the Conduit Flow Process Mode 2 for MODFLOW-2005

As a result of rock dissolution processes, karst aquifers exhibit highly conductive features such as caves and conduits. Within these structures, groundwater flow can become turbulent and therefore be described by nonlinear gradient functions. Some numerical groundwater flow models explicitly account for pipe hydraulics by coupling the continuum model with a pipe network that represents the conduit system. In contrast, the Conduit Flow Process Mode 2 (CFPM2) for MODFLOW-2005 approximates turbulent flow by reducing the hydraulic conductivity within the existing linear head gradient of the MODFLOW continuum model. This approach reduces the practical as well as numerical efforts for simulating turbulence. The original formulation was for large pore aquifers where the onset of turbulence is at low Reynolds numbers (1 to 100) and not for conduits or pipes. In addition, the existing code requires multiple time steps for convergence due to iterative adjustment of the hydraulic conductivity. Modifications to the existing CFPM2 were made by implementing a generalized power function with a user-defined exponent. This allows for matching turbulence in porous media or pipes and eliminates the time steps required for iterative adjustment of hydraulic conductivity. The modified CFPM2 successfully replicated simple benchmark test problems.     

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.     

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.     

Identification of heterogeneous aquifer transmissivity using an AE-based method

Determination of hydraulic head, H, as a function of spatial coordinates and time, in ground water flow is the basis for aquifer management and for prediction of contaminant transport. Several computer codes are available for this purpose. Spatial distribution of the transmissivity, T(x,y), is a required input to these codes. In most aquifers, T varies in an erratic manner, and it can be characterized statistically in terms of a few moments: the expected value, the variance, and the variogram. Knowledge of these moments, combined with a few measurements, permits one to estimate T at any point using geostatistical methods. In a review of transmissivity data from 19 unconsolidated aquifers, Hoeksema and Kitanidis (1985) identified two types of the logtransmissivity Y= ln(T) variations: correlated variations with variance sigma2Yc and correlation scale, I(Y), on the order of kilometers, and uncorrelated variations with variance sigma2Yn. Direct identification of the logtransmissivity variogram, Gamma(Y), from measurements is difficult because T data are generally scarce. However, many head measurements are commonly available. The aim of the paper is to introduce a methodology to identify the transmissivity variogram parameters (sigma2Yc, I(Y), and sigma2Yn) using head data in formations characterized by large logtransmissivity variance. The identification methodology uses a combination of precise numerical simulations (carried out using analytic element method) and a theoretical model. The main objective is to demonstrate the application of the methodology to a regional ground water flow in Eagle Valley basin in west-central Nevada for which abundant transmissivity and head measurements are available.     

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.     

Tide-induced head fluctuations in coastal aquifers of variable thickness

In this work, a new analytical solution to describe tide-induced head fluctuations in aquifers of variable thickness is presented. The proposed model assumes a finite and confined aquifer with a thickness that increases or decreases quadratically with the distance to the coast. A closed-form analytical solution is obtained by solving a boundary-value problem with both a separation of variables method and a change of variables method. This solution is a generalization of the solution obtained by Cuello et al., Hydrogeological Journal, 2017, 25, 1509–1515. The analytical solution is expressed in terms of the wedging parameter, a parameter that depends on the length and thicknesses at the coast and at the inland edge of the aquifer. Positive values of the wedging parameter describe aquifers with increasing thickness towards land and negative values describe aquifers with a decreasing thickness in the inland direction. The comparison of the new solution and the solution for a finite aquifer with constant thickness indicates that the sign of the wedging parameter enhances or decreases the amplitude of the tide-induced signal. However, the differences in time-lag between both solutions are negligible near the coast. The slope factor, which quantifies the inconsistencies between aquifer diffusivities estimated from attenuation and time-lag data, is computed and analysed. Near the coast, slope factor values greater than one are obtained for negative wedging parameters while slope factor values less than one are obtained for positive wedging parameters. The analysis of the new solution also indicates that more reliable estimates of the hydraulic diffusivity can be obtained from time-lag data.     

Uncertainty estimation of pathlines in ground water models

A method is proposed to estimate the uncertainty of the location of pathlines in two-dimensional, steady-state confined or unconfined flow in aquifers due to the uncertainty of the spatially variable unconditional hydraulic conductivity or transmissivity field. The method is based on concepts of the semianalytical first-order theory given in Stauffer et al. (2002, 2004), which allows estimates of the lateral second moment (variance) of the location of a moving particle. However, this method is reformulated in order to account for nonuniform recharge and nonuniform aquifer thickness. One prominent application is the uncertainty estimation of the catchment of a pumping well by considering the boundary pathlines starting at a stagnation point. In this method, the advective transport of particles is considered, based on the velocity field. In the case of a well catchment, backtracking is applied by using the reversed velocity field. Spatial variability of hydraulic conductivity or transmissivity is considered by taking into account an isotropic exponential covariance function of log-transformed values with parameters describing the variance and correlation length. The method allows postprocessing of results from ground water models with respect to uncertainty estimation. The code PPPath, which was developed for this purpose, provides a postprocessing of pathline computations under PMWIN, which is based on MODFLOW. In order to test the methodology, it was applied to results from Monte Carlo simulations for catchments of pumping wells. The results correspond well. Practical applications illustrate the use of the method in aquifers.     

A computer program for flow-log analysis of single holes (FLASH)

A new computer program, FLASH (Flow-Log Analysis of Single Holes), is presented for the analysis of borehole vertical flow logs. The code is based on an analytical solution for steady-state multilayer radial flow to a borehole. The code includes options for (1) discrete fractures and (2) multilayer aquifers. Given vertical flow profiles collected under both ambient and stressed (pumping or injection) conditions, the user can estimate fracture (or layer) transmissivities and far-field hydraulic heads. FLASH is coded in Microsoft Excel(5) with Visual Basic for Applications routines. The code supports manual and automated model calibration.     

Numerical modelling of mass transport in hydrogeologic environments: performance comparison of the Laplace Transform Galerkin and Arnoldi modal reduction schemes

The Laplace Transform Galerkin (LTG) method and the Arnoldi modal reduction method (AMRM) have been implemented in finite element schemes designed to solve mass transport problems in porous media by Sudicky [Sudicky, E.A., Water Resour. Res., 25(8) (1989) 1833-46] and Woodbury et al [Woodbury, A.D., Dunbar, W.S., & Nour-Omid, B., Water Resour. Res., 26(10) (1990) 2579-90]. In this work, a comparative analysis of the two methods is performed with attention focused on efficiency and accuracy. The analysis is performed over one- and two-dimensional domains composed of homogeneous and heterogeneous material properties. The results obtained using homogeneous material properties indicate that for a given mesh design the LTG method maintains a higher degree of accuracy than does the AMRM. However, in terms of efficiency, the Arnoldi attains a pre-defined level of accuracy faster than does the LTG method. It is also shown that for problems involving homogeneous material properties the solution obtained using the LTG method on a coarse mesh is comparable in terms of solution time and accuracy to that obtained using the AMRM on a fine mesh. Comparisons similar to those performed using homogeneous material properties are also performed for the case where the hydraulic conductivity field is heterogeneous. For this case, the level of accuracy achieved by the AMRM and the LTG method are similar. However, as with the analysis involving homogeneous material properties, the AMRM is found to be more efficient than the LTG method. It is also shown that for heterogeneous material properties, use of the LTG method under high grid Peclet conditions can be potentially problematic.     

Using data assimilation method to calibrate a heterogeneous conductivity field conditioning on transient flow test data

A data assimilation method is developed to calibrate a heterogeneous hydraulic conductivity field conditioning on transient pumping test data. The ensemble Kalman filter (EnKF) approach is used to update model parameters such as hydraulic conductivity and model variables such as hydraulic head using available data. A synthetical two-dimensional flow case is used to assess the capability of the EnKF method to calibrate a heterogeneous conductivity field by assimilating transient flow data from observation wells under different hydraulic boundary conditions. The study results indicate that the EnKF method will significantly improve the estimation of the hydraulic conductivity field by assimilating continuous hydraulic head measurements and the hydraulic boundary condition will significantly affect the simulation results. For our cases, after a few data assimilation steps, the assimilated conductivity field with four Neumann boundaries matches the real field well while the assimilated conductivity field with mixed Dirichlet and Neumann boundaries does not. We found in our cases that the ensemble size should be 300 or larger for the numerical simulation. The number and the locations of the observation wells will significantly affect the hydraulic conductivity field calibration.     

Groundwater response to tidal fluctuations in a leaky confined coastal aquifer with a finite length

The study on the hydraulic properties of coastal aquifers has significant implications both in hydrological sciences and environmental engineering. Although many analytical solutions are available, most of them are based on the same basic assumption that assumes aquifers extend landward semi‐infinitely, which does not necessarily reflect the reality. In this study, the general solutions for a leaky confined coastal aquifer have been developed that consider both finitely landward constant‐head and no‐flow boundaries. The newly developed solutions were then used to examine theoretically the joint effects of leakage and aquifer length on hydraulic head fluctuations within the leaky confined aquifer, and the validity of using the simplified solution, which assumes the aquifer is semi‐infinite. The results illustrated that the use of the simplified solution may cause significant errors, depending on joint effects of leakage and aquifer length. A dimensionless characteristic parameter was then proposed as an index for judging the applicability of the simplified solution. In addition, practical application of the general solution for the constant‐head inland boundary was used to characterize the hydraulic properties of a leaky confined aquifer using the data collected from a field site at the Seine River estuary, France, and the versatility of the general solution was further justified.     

Two-dimensional stochastic analysis of flow in leaky confined aquifers subject to spatial and periodic leakage

This paper addresses the issue of flow in heterogeneous leaky confined aquifers subject to leakage. The leakage into the confined aquifer is driven by spatial and periodic fluctuations of water table in an overlying phreatic aquifer. The introduction of leakage leads to non-uniformity in the mean head gradient and results in nonstationarity in hydraulic head and velocity fields. Therefore, a nonstationary spectral approach based on Fourier–Stieltjes representations for the perturbed quantities is adopted to account for the spatial variability of nonstationary head fields. Closed-form expressions for the variances of hydraulic head and specific discharge are developed in terms of statistical properties of hydraulic parameters. The results indicate that the spatiotemporal variations in leakage leads to enhanced variability of the hydraulic head and of the specific discharge, which increase with distance from any arbitrary reference point. The coefficient of leakage and the spatial structure of log transmissivity field and of the amplitude of water table fluctuation are critical in quantifying the variability of the hydraulic head and of the specific discharge.     

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.     

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.