E. A. Sudicky

Nonreactive and reactive solute transport in three‐dimensional heterogeneous porous media: Mean displacement, plume spreading, and uncertainty

The field-scale transport of reactive and nonreactive solutes by groundwater in a statistically anisotropic aquifer is examined by means of high-resolution, three-dimensional numerical solutions of the steady state flow and transient advection-dispersion equations. The presence of physical and chemical heterogeneities in the aquifer media is modeled with the use of a geostatistical description of the hydraulic conductivity and chemical distribution coefficient. The geostatistical parameters describing the spatial variations of the log-transformed hydraulic conductivity fields, ln [K(x)], are chosen to resemble those of the Borden aquifer. For sorbing solutes the spatial variation in the log-transformed distribution coefficient fields, ln [Kd(x)], is generated so as to have a geometric mean and variance similar to that estimated by Durant (1986) for the organic chemical tetrachlorethane with Borden sand. It is further assumed that the ln [K(x)] and ln [Kd(x)] fields are inversely correlated to each other and that they possess the same spatial correlation structure. Five realizations of media which incorporate these characteristics are generated by means of the Fourier spectral technique of Robin et al. (1993). It is shown that joint K(x) and Kd(x) variability can impart a large-scale “pseudokinetic” behavior, in that the ensemble mean bulk retardation factor can increase with time and plume displacement distance, even though sorption is modeled as being linear and instantaneous at the scale of any single heterogeneity and the flow field is assumed to be steady state. From realization to realization, however, the apparent bulk retardation factor can either increase dramatically at early time in a manner similar to that observed by Roberts et al. (1986) during the Borden tracer test or it can decrease. At large time the ensemble mean velocity of the centroid of the reactive plume is close to a value given by the mean fluid velocity divided by the arithmetic mean of the locally variable retardation factors. Local-scale transport nonidealities, such as intraparticle diffusion and relatively rapid kinetic sorption, are shown to have minimal influence on the plume centroid velocity and macrodispersivity, relative to the effect of aquifer heterogeneity. The results of the numerical simulations further demonstrate that the first-order stochastic analyses of Gelhar and Axness (1983), Dagan (1988), and Naff (1990) tend to overestimate the actual field-scale longitudinal spreading of a nonreactive solute. This result is believed to occur because these analyses include the artificial effect of plume centroid dispersion about the ensemble mean position, in addition to the actual spreading of each plume realization about its center of mass. When the effect of plume centroid dispersion is added to the numerical simulation results, a reasonable agreement with the solution of Naff (1990), which takes into consideration local-scale dispersion, is obtained for the longitudinal macrodispersivity. In the transverse direction, the ensemble mean spreading agrees reasonably well with the solution of Dagan (1988). It is also shown that the longitudinal macrodispersivity of a reactive solute can be enhanced relative to that of a nonreactive one. Its actual value, however, is overpredicted by the first-order solution of Garabedian (1987) for the reasons given above. Also explored, in a preliminary fashion, are some issues related to prediction uncertainty for both reactive and nonreactive solutes migrating through heterogeneous aquifers.

Cross‐correlated random field generation with the direct Fourier Transform Method

This paper presents a computer algorithm that is capable of cogenerating pairs of three-dimensional, cross-correlated random fields. The algorithm produces random fields of real variables by the inverse Fourier transform of a randomized, discrete three-dimensional spectral representations of the variables. The randomization is done in the spectral domain in a way that preserves the direct power and cross-spectral density structure. Two types of cross spectra were examined. One type specifies a linear relationship between the two fields, which produces the same correlation scales for both variables but different variances. The second cross spectrum is obtained from a specified transfer function and the two power spectra, and it produces fields with different correlation scales. For both models the degree of correlation is specified by the coherency. A delay vector can also be specified to produce an out-of-phase correlation between the two fields. The algorithm is very efficient computationally, is relatively easy to use, and does not produce the lineation problems that can be encountered with the turning bands method. Perhaps most important, this random field generator is capable of co-generating cross-correlated random fields.

Spatial Variability of Strontium Distribution Coefficients and Their Correlation With Hydraulic Conductivity in the Canadian Forces Base Borden Aquifer

Distribution coefficients (Kd), defined as the ratio of the concentration of solute associated with the solids to the concentration in solution, are widely used in the prediction of reactive solute transport. With the advent of stochastic approaches to describe solute transport, there is a need to examine the spatial distribution of Kd, and its correlation with the hydraulic conductivity (K). Distribution coefficients were measured in triplicates for strontium on 1279 subsamples of cores from Canadian Forces Base Borden for which K measurements were available. The Kd values ranged from 4.4 to 29.8 mL/g, with a mean of 9.9 and standard deviation of 2.89 mL/g. The standard error on the triplicate means was 0.95 mL/g or approximately 10% of the mean. The spatial behavior of Kd and K (expressed as In (Kd) and ln (K)) was examined in three directions: horizontally along two orthogonal transects and vertically. The two variables each behaved nearly identically in the two horizontal directions, suggesting horizontal isotropy. Horizontally, ln (Kd) appeared as “white noise” suggesting that the horizontal spacing between cores (1 m) was too large to detect any self-correlation. The distribution coefficient displayed increasing power spectral density with increasing scale in the vertical direction, while In (K) showed these trends in all directions. Depending on the model used, the, correlation lengths obtained by least squares fits of the power spectra varied from 1 to 7.5 m horizontally and from 10 to 30 cm vertically for ln (K); and from 30 cm to 2 m horizontally and from 30 to 70 cm vertically for ln (Kd). The ln (Kd) values showed a significant but very weak negative overall correlation with ln (K) at the 99.95% confidence level. The cross-spectral and coherency analysis showed that the sign and degree of correlation between ln (Kd) and ln (K) depended on the scale and direction considered. The correlations in all directions and at all scales were weak, and could not always be declared significantly different than zero at the 95% confidence level.

Tritium and Helium 3 Isotope Ratios for Direct Estimation of Spatial Variations in Groundwater Recharge

In the absence of strong dispersion the ratio of atmospherically derived tritium (3H) and its stable daughter, helium 3 (3He), can be used to accurately date shallow groundwater with ages ranging from 0 to about 50 years. The effects of dispersion on 3H/3He ages are examined by numerical simulation in simple one- and two-dimensional flow systems. For a 3H input function that is relatively constant in time the difference between 3H/3He ages and purely advective groundwater travel times is controlled by the relative magnitude of the groundwater velocity, the 3H decay constant, and the dispersivity of the porous media. In many shallow sandy aquifers the groundwater velocity is sufficiently large and the dispersivity sufficiently small that the difference between the true groundwater travel time and the 3H/3He age near the water table can be neglected (<10% difference). This suggests that measurements of the 3H/3He age gradient near the water table can be used to estimate recharge rates and vertical velocity values with a high degree of spatial resolution. In contrast to 3H activities, modeled 3H/3He ages are relatively insensitive to uncertainties with respect to the 3H input function. Thus combined measurements of 3H and 3He offer considerable advantage over measurements of 3H alone. The ability to obtain reliable estimates of spatial variations in recharge rates combined with accurate groundwater ages in shallow unconfined aquifers will undoubtedly provide valuable constraints when calibrating models of groundwater flow and solute transport in such systems.

Application of a fully‐integrated surface‐subsurface flow model at the watershed‐scale: A case study

[1]xa0Results are presented in which a physically-based, three-dimensional model that fully integrates surface and variably-saturated subsurface flow processes is applied to the 75 km2 Laurel Creek Watershed within the Grand River basin in Southern Ontario, Canada. The primary objective of this study is to gauge the models ability to reproduce surface and subsurface hydrodynamic processes at the watershed scale. Our objective was first accomplished by calibrating the steady-state subsurface portion of the system to 50 observation wells where hydraulic head data were available, while simultaneously matching the stream baseflow discharge. The level of agreement between the observed and computed subsurface hydraulic head values, baseflow discharge and the spatial pattern of the surface drainage network indicates that the model captures the essence of the surface-subsurface hydraulic characteristics of the watershed. The calibrated model is then subjected to two time series of input rainfall data and the calculated discharge hydrographs are compared to the observed rainfall-runoff responses. The calculated and observed rainfall-runoff responses were shown to agree moderately well for both rainfall data series that were utilized. Additionally, the spatial and temporal responses of the watershed with respect to the overland flow areas contributing to streamflow and the surface-subsurface exchange fluxes across the land surface during rainfall inundation and subsequent drainage phases demonstrate the dynamic nature of the interaction occurring between the surface and subsurface hydrologic regimes. Overall, it is concluded that it is feasible to apply a fully-integrated, surface/variably-saturated subsurface flow model at the watershed scale and possibly larger scales.

Colloid‐facilitated contaminant transport in discretely fractured porous media: 1. Numerical formulation and sensitivity analysis

A two-dimensional numerical model is developed that incorporates the mechanism of colloid-facilitated transport in discretely fractured porous media. The numerical model accounts for aqueous phase contaminant transport in the fractures and the porous matrix, colloid transport in the fractures, and sorption of the solute. Deep-bed filtration of the colloids is accounted for, and the solute is allowed to sorb on both the mobile and filtered colloids. The numerical formulation allows for either equilibrium or kinetic sorption reactions onto the fracture walls, the matrix solids, and the mobile and filtered colloids according to either a Langmuir or a Freundlich isotherm. The results of a series of simulations involving a system of parallel fractures explore the importance of mobile colloids on contaminant migration and indicate that if sorption onto the colloids is a slow kinetic process, then the mobile colloids may lead to significantly enhanced contaminant migration.

Mechanisms Controlling Vacuum Extraction Coupled With Air Sparging for Remediation of Heterogeneous Formations Contaminated by Dense Nonaqueous Phase Liquids

The numerical model CompFlow is used to study the mechanisms controlling vacuum extraction, coupled with air sparging, as a means for remediation of heterogeneous formations contaminated with dense nonaqueous phase liquids (DNAPLs). Two dominant mechanisms are demonstrated to control this remediation technology. First, at early times, the gas phase directly contacts the DNAPL, particularly in the unsaturated zone, causing relatively rapid transfer of contaminant from the nonaqueous phase to the gas phase and subsequent removal by the vacuum extractor. Second, at later times, remediation is controlled by the transfer of contaminant from the nonaqueous phase to the aqueous phase below the water table. During this time the vacuum extractor pumps both liquid and vaporized water in the aqueous and gas phases. This causes the contaminant that is dissolved in the aqueous phase to migrate vertically upward across the permeability layers toward the vacuum extractor where it is removed. This intermediate to late time removal mechanism is shown to be controlled by contaminant dissolution, which is a slower transfer process than the direct DNAPL vaporization that occurs at early time. Our analysis indicates that as long as both air and water are actively flushed through the DNAPL zone, both early-time vaporization and intermediate- to late-time dissolution are effective mechanisms leading to the removal of the DNAPL. We show that it may be possible to design the remedial system so as to reduce its performance sensitivity to geologic heterogeneity. A lack of sensitivity of a remedial design to heterogeneity is highly desirable because a robust design implies that the degree of site characterization required for reasonable success will be less than that needed for a less robust scheme.

Influence of Leaky Boreholes on Cross‐Formational Groundwater Flow and Contaminant Transport

Abandoned and improperly sealed boreholes, monitoring wells, and water supply wells are common features at many contaminated sites. These features can act as conduits that transmit contaminants between aquifers separated by otherwise continuous aquitards. In this work the leaky boreholes are represented as highly conductive one-dimensional line elements superimposed onto a mesh of three-dimensional finite elements representing the porous medium. Simulation results are presented for a series of scenarios involving a simple hydrogeologic setting composed of an upper confined aquifer, a middle aquitard, and a lower confined aquifer. The simulations examine the effect of varying the borehole properties, and vertical hydraulic gradient across the aquitard, and the borehole location. The results show that a contaminant can rapidly migrate downward along a leaky borehole and create an extensive plume in the lower aquifer, even if the borehole is filled with aquifer sediments. If the borehole is an open feature across the aquitard, the entire plume, or a significant portion of it, that is migrating into the upper aquifer can be diverted into the lower one if the vertical hydraulic gradient across the aquitard can be diverted into the lower one if the vertical hydraulic gradient across the aquitard is sufficiently strong.morexa0» The peak concentration arriving at a pumping well located in the lower aquifer and the time of arrival are functions of the proximity of the leaky borehole to the pumping well and its angular offset from the central flow line passing through the surficial source. Overpressurization of the lower aquifer due to injection can overcome downward preexisting hydraulic gradients across the aquitard such that contaminants can rapidly migrate upward along the leaky borehole and cause contamination of the otherwise protected upper aquifer. 26 refs., 7 figs., 1 tab.«xa0less

Numerical analysis of solute migration through fractured clayey deposits into underlying aquifers

This study examines the degree to which vertically fractured clayey aquitards protect an underlying aquifer from near-surface sources of contamination. Evidence from field studies indicates that many clayey aquitards previously assumed to be unfractured at depth have vertical fractures that actively transmit groundwater. Calculated groundwater velocities within these fractures can exceed several meters per day for apertures in the range of those reported to occur in clays. If this is the case, significant quantities of dissolved contaminants originating in near-surface zones may move rapidly downward through the fracture network into underlying aquifers. A two-dimensional Laplace transform Galerkin finite element model (Sudicky and McLaren, this issue) was used to assess the importance of idealized planar vertical fractures on flow and solute transport through a 15-m-thick aquitard and subsequent plume evolution within an underlying aquifer in which groundwater flows horizontally. Sensitivity analyses indicate that solute transport through the aquitard is strongly affected by downward flow in widely spaced deep fractures with apertures as small as 10 μm and that a plume of large lateral extent and significant concentration forms within the underlying aquifer in a few tens of years for apertures of the order of 20 μm or much sooner if the apertures are 50 μm. This finding suggests that small hydraulically active fractures that are very difficult to detect or identify within an otherwise protective clayey barrier may cause rapid and large-scale contamination of groundwater beneath the barrier.

Heterogeneity in hydraulic conductivity and its role on the macroscale transport of a solute plume: From measurements to a practical application of stochastic flow and transport theory

[1]xa0The spatial variability of hydraulic conductivity in a shallow unconfined aquifer located at North Bay, Ontario, composed of glacial-lacustrine and glacial-fluvial sands, is examined in exceptional detail and characterized geostatistically. A total of 1878 permeameter measurements were performed at 0.05 m vertical intervals along cores taken from 20 boreholes along two intersecting transect lines. Simultaneous three-dimensional (3-D) fitting of Ln(K) variogram data to an exponential model yielded geostatistical parameters for the estimation of bulk hydraulic conductivity and solute dispersion parameters. The analysis revealed a Ln(K) variance equal to about 2.0 and 3-D anisotropy of the correlation structure of the heterogeneity (λ1, λ2, and λ3 equal to 17.19, 7.39, and 1.0 m, respectively). Effective values of the hydraulic conductivity tensor and the value of the longitudinal macrodispersivity were calculated using the theoretical expressions of Gelhar and Axness (1983). The magnitude of the longitudinal macrodispersivity is reasonably consistent with the observed degree of longitudinal dispersion of the landfill plume along the principal path of migration. Variably saturated 3-D flow modeling using the statistically derived effective hydraulic conductivity tensor allowed a reasonably close prediction of the measured water table and the observed heads at various depths in an array of piezometers. Concomitant transport modeling using the calculated longitudinal macrodispersivity reasonably predicted the extent and migration rates of the observed contaminant plume that was monitored using a network of multilevel samplers over a period of about 5 years. It was further demonstrated that the length of the plume is relatively insensitive to the value of the longitudinal macrodispersivity under the conditions of a steady flow in 3-D and constant source strength. This study demonstrates that the use of statistically derived parameters based on stochastic theories results in reliable large-scale 3-D flow and transport models for complex hydrogeological systems. This is in agreement with the conclusions reached by Sudicky (1986) at the site of an elaborate tracer test conducted in the aquifer at the Canadian Forces Base Borden.