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.     

Dissolution of nonaqueous phase liquid pools in anisotropic aquifers

A two-dimensional numerical transport model is developed to determine the effect of aquifer anisotropy and heterogeneity on mass transfer from a dense nonaqueous phase liquid (DNAPL) pool. The appropriate steady state groundwater flow equation is solved implicitly whereas the equation describing the transport of a sorbing contaminant in a confined aquifer is solved by the alternating direction implicit method. Statistical anisotropy in the aquifer is introduced by two-dimensional, random log-normal hydraulic conductivity field realizations with different directional correlation lengths. Model simulations indicate that DNAPL pool dissolution is enhanced by increasing the mean log-transformed hydraulic conductivity, groundwater flow velocity, and/or anisotropy ratio. The variance of the log-transformed hydraulic conductivity distribution is shown to be inversely proportional to the average mass transfer coefficient.     

Scale effect and calibration of contaminant transport models

The influence on solute transport of the small-scale spatial variation of aquifer hydraulic conductivity (K) was analyzed by comparing results from fine-grid (2 m by 2 m) simulations of a synthetic heterogeneous aquifer to those from coarse-grid (8 m by 4 m) simulations of an equivalent homogeneous aquifer. Realizations of the K field of the heterogeneous aquifer were generated, using the Monte Carlo approach, from a lognormal distribution with mean log K of 2 (K in m/d) and three levels of log K variance of 0.1, 0.5, and 1.0. Numerical simulation results show that the average standard deviation of point concentrations increased from 1.21 to 5.78 when the value of log K variance was increased from 0.1 to 1.0. The average discrepancy between modeled concentrations (obtained from a coarse-grid deterministic numerical simulation) and the actual mean point concentrations (obtained from fine-grid Monte Carlo numerical simulations) increased from 0.91 to 4.23 with the increase in log K variance. The results from this study illustrate the uncertainty in predictions from contaminant transport models due to their inability to simulate the effects of heterogeneities at scales smaller than the model grid.     

A quasi three-dimensional model for flow and transport in unsaturated and saturated zones: 1. Implementation of the quasi two-dimensional case

A quasi three-dimensional (QUASI 3-D) model is presented for simulating the subsurface water flow and solute transport in the unsaturated and in the saturated zones of soil. The model is based on the assumptions of vertical flow in the unsaturated zone and essentially horizontal groundwater flow. The 1-D Richards equation for the unsaturated zone is coupled at the phreatic surface with the 2-D flow equation for the saturated zone. The latter was obtained by averaging 3-D flow equation in the saturated zone over the aquifer thickness. Unlike the Boussinesq equation for a leaky-phreatic aquifer, the developed model does not contain a storage term with specific yield and a source term for natural replenishment. Instead it includes a water flux term at the phreatic surface through which the Richards equation is linked with the groundwater flow equation. The vertical water flux in the saturated zone is evaluated on the basis of the fluid mass balance equation while the horizontal fluxes, in that equation, are prescribed by Darcy law. A 3-D transport equation is used to simulate the solute migration. A numerical algorithm to solve the problem for the general quasi 3-D case was developed. The developed methodology was exemplified for the quasi 2-D cross-sectional case (QUASI2D). Simulations for three synthetic problems demonstrate good agreement between the results obtained by QUASI2D and two fully 2-D flow and transport codes (SUTRA and 2DSOIL). Yet, simulations with the QUASI2D code were several times faster than those by the SUTRA and the 2DSOIL codes.     

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.     

Non-Fickian mass transport in fractured porous media

The paper provides an introduction to fundamental concepts of mathematical modeling of mass transport in fractured porous heterogeneous rocks. Keeping aside many important factors that can affect mass transport in subsurface, our main concern is the multi-scale character of the rock formation, which is constituted by porous domains dissected by the network of fractures. Taking into account the well-documented fact that porous rocks can be considered as a fractal medium and assuming that sizes of pores vary significantly (i.e. have different characteristic scales), the fractional-order differential equations that model the anomalous diffusive mass transport in such type of domains are derived and justified analytically. Analytical solutions of some particular problems of anomalous diffusion in the fractal media of various geometries are obtained. Extending this approach to more complex situation when diffusion is accompanied by advection, solute transport in a fractured porous medium is modeled by the advection-dispersion equation with fractional time derivative. In the case of confined fractured porous aquifer, accounting for anomalous non-Fickian diffusion in the surrounding rock mass, the adopted approach leads to introduction of an additional fractional time derivative in the equation for solute transport. The closed-form solutions for concentrations in the aquifer and surrounding rocks are obtained for the arbitrary time-dependent source of contamination located in the inlet of the aquifer. Based on these solutions, different regimes of contamination of the aquifers with different physical properties can be readily modeled and analyzed.     

Fitting a stochastic partial differential equation to aquifer data

The steady state two dimensional groundwater flow equation with constant transmissivities was studied by Whittle in 1954 as a stochastic Laplace equation. He showed that the correlation function consisted of a modified Bessel function of the second kind, order 1, multiplied by its argument. This paper uses this pioneering work of Whittle to fit an aquifer head field to unequally spaced observations by maximum likelihood. Observational error is also included in the model. Both the isotropic and anisotropic cases are considered. The fitted field is then calculated on a two dimensional grid together with its standard deviation. The method is closely related to the use of two-dimensional splines for fitting surfaces to irregularly spaced observations.     

Modeling Transport in Transient Ground-Water Flow: An Unacknowledged Approximation

Abstract. During unsteady or transient ground-water flow, the fluid mass per unit volume of aquifer changes as the potentiometric head changes, and solute transport is affected by this change in fluid storage. Three widely applied numerical models of two-dimensional transport partially account for the effects of transient flow by removing terms corresponding to the fluid continuity equation from the transport equation, resulting in a simpler governing equation. However, fluid-storage terms remaining in the transport equation that change during transient flow are, in certain cases, held constant in time in these models. For the case of increasing heads, this approximation, which is unacknowledged in these models'documentation, leads to transport velocities that are too high, and increased concentration at fluid and solute sources. If heads are dropping in time, computed transport velocities are too low. Using parameters that somewhat exaggerate the effects of this approximation, an example numerical simulation indicates solute travel time error of about 14 percent but only minor errors due to incorrect dilution volume. For horizontal flow and transport models that assume fluid density is constant, the product of porosity and aquifer thickness changes in time: initial porosity times initial thickness plus the change in head times the storage coefficient. This formula reduces to the saturated thickness in unconfined aquifers if porosity is assumed to be constant and equal to specific yield. The computational cost of this more accurate representation is insignificant and is easily incorporated in numerical models of solute transport.     

Effects of rate-limited mass transfer on water sampling with partially penetrating wells

Nonequilibrium concentration type curves are numerically developed and sensitivity analyses are performed to examine the relationships between effluent concentrations in partially penetrating monitoring/extraction wells, the vertical plume shape, and the mass transfer characteristics of the aquifer. The governing two-dimensional, axisymmetric nonequilibrium solute transport equation is solved in three stages using an operator-splitting approach. In the first two stages, the advection and dispersion terms are solved with the Eulerian-Lagrangian method, based on the backward method of characteristics for advection and the standard implicit Galerkin finite element method for dispersion. In the third step, the first-order, immobile-mobile domain mass transfer term is computed analytically for both two-site and lognormally distributed, multirate models. Effluent concentration variations with time and contour plots of the pore water concentration distribution in the aquifer are compared for a wide range of field- and laboratory-measured mass transfer rates, various plume shapes, and relevant physical/chemical parameter values, including pumping rate, vertical anisotropy ratio, retardation factor, and porosity. The simulation results show that rate-limited mass transfer can have a significant impact on sample and aquifer pore water concentrations during three-dimensional transport to a partially penetrating well. An alternative dimensionless form of the nonequilibrium solute transport equation is derived to illustrate the key parameter groupings that quantify rate-limited sorption effects and show the relative importance of individual parameters. A hypothetical field application example demonstrates the fitting of dimensional type curves to discrete-interval sampling data in order to evaluate the mass transfer characteristics of an aquifer and shows how type curve superposition can be used to model complex plume shapes.     

Can Contaminant Transport Models Predict Breakthrough?

A solute breakthrough curve measured during a two-well tracer test was successfully predicted in 1986 using specialized contaminant transport models. Water was injected into a confined, unconsolidated sand aquifer and pumped out 125 feet (38.3 m) away at the same steady rate. The injected water was spiked with bromide for over three days; the outflow concentration was monitored for a month. Based on previous tests, the horizontal hydraulic conductivity of the thick aquifer varied by a factor of seven among 12 layers. Assuming stratified flow with small dispersivities, two research groups accurately predicted breakthrough with three-dimensional (12-layer) models using curvilinear elements following the arc-shaped flowlines in this test.
Can contaminant transport models commonly used in industry, that use rectangular blocks, also reproduce this breakthrough curve? The two-well test was simulated with four MODFLOW-based models, MT3D (FD and HMOC options), MODFLOWT, MOC3D, and MODFLOW-SURFACT.
Using the same 12 layers and small dispersivity used in the successful 1986 simulations, these models fit almost as accurately as the models using curvilinear blocks. Subtle variations in the curves illustrate differences among the codes. Sensitivities of the results to number and size of grid blocks, number of layers, boundary conditions, and values of dispersivity and porosity are briefly presented. The fit between calculated and measured breakthrough curves degenerated as the number of layers and/or grid blocks decreased, reflecting a loss of model predictive power as the level of characterization lessened. Therefore, the breakthrough curve for most field sites can be predicted only qualitatively due to limited characterization of the hydrogeology and contaminant source strength.     

On the aquifer's integrated balance equations

A general methodology is presented for describing transport phenomena in porous media at a macroscopic level. Then, these macroscopic balance equations are integrated (or averaged) along the vertical for confined, leaky and phreatic aquifers.The results are employed to derive (averaged) aquifer equations for the flow of water and of a solute (hydrodynamic dispersion). It is shown that in all cases, the resulting equation is identical to that derived on the basis of an assumption of horizontal flow (the Dupuit assumption).Macrodispersion, occurring at the aquifer level, is discussed and appropriate coefficients are proposed.     

Simulation of long-term transient density-dependent transport in groundwater

A numerical model for the economical simulation of long-term transient response in density-dependent transport problems is introduced. Although a classical Galerkin finite element approach is used, emphasis on optimum efficiency throughout the development results in a scheme that is found to be significantly less costly than comparable existing schemes. This advantage in efficiency increases the scope of simulation problems that can be handled within the constraints of a limited research budget. Some distinctive aspects are the elimination of static quantities in the fluid continuity equation, achieved by the introduction of equivalent freshwater head, and the elimination of numerical integration, achieved by the deliberate choice of linear elements. As a result of this choice, fluid velocities are discontinuous across the element boundaries. It is shown, however, that the solution obtained with discontinuous velocities approaches that obtained with continuous velocities as the grid is refined, and that the two types of solutions give essentially the same results when the elements are in the same size range. The model is applied to simulate the complete transient response for a well-known problem of seawater intrusion in a confined aquifer. The simulation is performed with both the constant dispersion coefficient used by previous researchers, and a more physically realistic velocity-dependent dispersion coefficient. Responses are found to be substantially different for the two types of coefficients, with the velocity-dependent dispersion coefficient producing much slower convergence to a state of dynamic equilibrium, and a much more pointed saltwater toe, which at the bottom of the aquifer tends to a sharp interface at equilibrium. Finally, it is shown by means of large-scale applications that the model is capable of efficiently simulating the long-term transient response in systems of practical significance.     

Modified Theis equation by considering the bending effect of the confining unit

The Theis equation has been widely used to study the transient movement of groundwater as a result of pumping in a confined aquifer. It is well known that the observed drawdown at early times has an obvious departure from the theoretical drawdown based on the Theis equation. The Theis equation was derived under the assumption that total stress in the aquifer was constant and the mechanical behavior of the confining unit was neglected. However, most geological formations, especially those which are well consolidated, have rigidity and therefore may bend like a plate to a certain extent. The increase in the effective stress in the aquifer due to pumping may not contribute entirely to the compression of the aquifer, but may be partially cancelled out by bending of the overlying aquitard. This means only a part of the total stress is used to compact the aquifer, or the aquifer cannot produce as much water as estimated from the Theis equation. This paper investigated the impact of the bending effect of the confining unit on drawdown. An analytical model which couples flow in the aquifer and bending of the confining unit was presented. The theory is based on elastic plates and solutions were given to the drawdown of groundwater level and deflection of the overlying formation. The drawdown estimated from the new equation was compared with that from the Theis equation. It can be concluded that drawdown from the Theis equation is less than the drawdown predicted by including the bending effect of the confining unit. Both a hypothetical example and a field pumping test in Shandong Province, China, were used to demonstrate the bending effect of the confining unit in the analysis of pumping test data. This paper demonstrated that the initial disagreement between observed drawdown and the Theis solution could be caused by the bending effect of the confining unit, a phenomenon not well addressed in traditional pumping test analysis. A quantitative understanding of this phenomenon can provide improved guidelines for analyzing drawdown data in a confined aquifer.     

Probabilistic modeling of aquifer heterogeneity using reliability methods

A probabilistic model of groundwater contaminant transport is presented. The model is based on coupling first- and second-order reliability methods (FORM and SORM) with a two-dimensional finite element solution of groundwater transport equations. Uncertainty in aquifer media is considered by modeling hydraulic conductivity as a spatial random field with a prescribed correlation structure. FORM and SORM provide the probability that a contaminant exceeds a target level at a well, termed the probability of failure. Sensitivity of the probability of failure to basic variabilities in grid block conductivity is also obtained, at no additional computational effort. The effect of the choice of the predetermined target level at the observation well is provided, along with its effect on the relevant sensitivity information. Considerable saving in computational time was achieved by superimposing a coarse random variables mesh on a finer numerical mesh. The presence of regions of lower conductivity on the probabilistic event is analyzed, and the regions in which conductivity most affects the results are identified.     

An investigation of numerical grid effects in parameter estimation

Modern ground water characterization and remediation projects routinely require calibration and inverse analysis of large three-dimensional numerical models of complex hydrogeological systems. Hydrogeologic complexity can be prompted by various aquifer characteristics including complicated spatial hydrostratigraphy and aquifer recharge from infiltration through an unsaturated zone. To keep the numerical models computationally efficient, compromises are frequently made in the model development, particularly, about resolution of the computational grid and numerical representation of the governing flow equation. The compromise is required so that the model can be used in calibration, parameter estimation, performance assessment, and analysis of sensitivity and uncertainty in model predictions. However, grid properties and resolution as well as applied computational schemes can have large effects on forward-model predictions and on inverse parameter estimates. We investigate these effects for a series of one- and two-dimensional synthetic cases representing saturated and variably saturated flow problems. We show that "conformable" grids, despite neglecting terms in the numerical formulation, can lead to accurate solutions of problems with complex hydrostratigraphy. Our analysis also demonstrates that, despite slower computer run times and higher memory requirements for a given problem size, the control volume finite-element method showed an advantage over finite-difference techniques in accuracy of parameter estimation for a given grid resolution for most of the test problems.     

Groundwater,salt and contaminant transport in a coastal aquifer system with a low-permeability layer in the intertidal zone

Low-permeability layer (LPL), formed by natural deposit or artificial reclamation and commonly found below the intertidal zone of coastal groundwater system, can retard the ingress of seawater and contaminants, and shorten the travel time of the land-sourced contaminant to the marine environment compared with a homogenous sandy coastal aquifer. However, there is limited understanding on how an intertidal LPL, a condition occurred in a coastal aquifer at Moreton Bay, Australia, influences the groundwater and contaminant transport across the shallow beach aquifer system. We characterized the aquifer hydrological parameters, monitored the in situ groundwater heads, and constructed a 2-D numerical model to analyses the cross-shore hydrological processes in this stratified system. The calibrated model suggests that in the lower aquifer, the inland-source fresh groundwater flowed horizontally towards the sea, upwelled along the freshwater–saltwater interface, and exited the aquifer at the shore below the LPL. Whereas in the upper aquifer, the tidally driven seawater circulation formed a barrier that prevented fresh groundwater from horizontal transport and discharge to the beach above the LPL, thereby directing its leakage to the lower aquifer. A contaminant represented by a conservative tracer was ‘released’ the upper aquifer in the model and results showed that the spreading extent of the contaminant plume, the maximum rate of contaminant discharge to the ocean, and its plume length decreased compared with a simulation case in a homogenous sandy aquifer. Sensitivity analysis was also conducted to investigate the characteristics of the LPL, including its continuity and hydraulic conductivity, which were found to vary along the beach at Moreton Bay. The result shows that with a lower hydraulic conductivity and continuous layer of LPL reduced the groundwater exchange and contaminant transport between upper and lower aquifer. The findings from the combined field and modelling investigations on the impact of an intertidal LPL on coastal aquifer systems highlight its significant implications to alter the groundwater and mass transport across the land–ocean interface.     

Beitrag zur analytischen Beschreibung von Transportprozessen im Grundwasserleiter bei diffusem Stoffeintrag

A model is described for the relation between water quantity and water quality in the aquifer with the diffuse substance input by the gravitation water being taken into account. The model is based on a partial differential equation of the second order for a flow pipe. For constant transport parameters, a sink of the kinetics of the first order and a time-dependent source term an exact analytical solution is presented and explained by a model. At a diffuse substance input, the dispersion has only a slight influence on the transport of substances contained in water. Potential applications of these solutions to different problems are mentioned.     

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.     

An isotopic investigation of a near-surface groundwater system

ABSTRACT

The distribution of environmental tritium, deuterium and oxygen-18 in the unsaturated zone and the underlying sandy phreatic aquifer was studied throughout 1981 in an area of high pine forests in the Rhine valley near Heidelberg. The observed vertical distribution of isotopes in the unsaturated zone can be satisfactorily explained by the combined use of a multi-cell model for moisture transport and an evapotranspiration model. The distribution in the underlying aquifer of the tracer input at the water table obtained using this method is found by replacing the total vertical diffusion coefficient in the diffusion equation with the dispersion coefficient. In this way observed tritium profiles are satisfactorily simulated for the period 1966–1981. The stable isotope profile in the unsaturated zone however remains largely unexplained due to inadequate data on the stable isotope content of precipitation over the investigated area.     

Approximate solution for solute transport during spherical-flow push-pull tests

An approximate analytical solution to the advection-dispersion equation was derived to describe solute transport during spherical-flow conditions in single-well push-pull tests. The spherical-flow case may be applicable to aquifer tests conducted in packed intervals or partially penetrating wells. Using results of two-dimensional numerical simulations, we briefly illustrate the applicability of the derived spherical-flow solution and provide a comparison with its cylindrical-flow counterpart. Good agreement between simulated extraction-phase breakthrough curves and the spherical-flow solution was found when the length of the injection/extraction region was small compared to both aquifer thickness and maximum solute frontal position at the end of the injection phase. On the other hand, discrepancies between simulated breakthrough curves and the spherical-flow solution increased with increasing anisotropy in hydraulic conductivities. Several inherent limitations embedded in its derivation such as assumptions of isotropy and homogeneity warrant the cautious use of the spherical-flow solution.