Edward A. Sudicky

Migration of contaminants in groundwater at a landfill: A case study: 1. Groundwater flow and plume delineation

Abstract A landfill-derived contaminant plume with a maximum width of ∼600 m, a length of ∼700 m and a maximum depth of 20 m in an unconfined sand aquifer was delineated by means of a monitoring network that includes standpipe piezometers, multilevel point-samplers and bundle-piezometers. The extent of detectable contamination caused by the landfill, which began operation in 1940 and which became inactive in 1976, was determined from the distributions of chloride, sulfate and electrical conductance in the sand aquifer, all of which have levels in the leachate that are greatly above those in uncontaminated groundwater. The maximum temperature of groundwater in the zone of contamination beneath the landfill is 12°C, which is 4–5°C above background. The thermal plume in the aquifer extends ∼150 m downgradient from the centre of the landfill. A slight transient water-table mound exists beneath the landfill in the late spring and summer in response to snowmelt and heavy rainfall. Beneath the landfill, the zone of leachate contamination extends to the bottom of the aquifer, apparently because of transient downward components of hydraulic gradient caused by the water-table mound and possibly because of the higher density and lower viscosity of the contaminated water. Values of hydraulic conductivity, which show variations due to local heterogeneity, were obtained from slug tests of piezometers, from pumping tests and from laboratory tests. Because of the inherent uncertainty in the aquifer parameter values, the 38-yr. frontal position of the plume calculated using the Darcy equation with the assumption of plug flow can differ from the observed frontal position by many hundreds of metres, although the use of mean parameter values produces a close agreement. The width of the plume is large relative to the width of the landfill and can be accounted for primarily by variable periods of lateral east- and westward flow caused by changes in water-table configuration due to the variable nature of recharge. Northward from the landfill, the vertical thickness of the plume decreases and the top of the plume is farther below the water table. The thickness of the zone of uncontaminated groundwater above the plume increases northward as the area of recharge of uncontaminated water downflow from the landfill increases. Because dispersion in the vertical direction is weak, there is very little mixing between the overlying zone of recharge water and the contaminant plume. Concentration profiles are irregular beneath and near the landfill and become smooth downgradient where the maximum concentrations are much less than those beneath landfill. These features are attributed to a strong influence of longitudinal dispersion. The plume passes beneath a small shallow stream near the landfill without significant influence on the stream.

Three-dimensional analysis of variably-saturated flow and solute transport in discretely-fractured porous media

A discrete fracture, saturated-unsaturated numerical model is developed where the porous matrix is represented in three dimensions and fractures are represented by two-dimensional planes. This allows a fully three-dimensional description of the fracture network connectivity. Solute advection and diffusion in the porous matrix are also directly accounted for. The variably-saturated flow equation is discretized in space using a control volume finite-element technique which ensures fluid conservation both locally and globally. Because the relative permeability and saturation curves for fractures may be highly nonlinear, and in strong contrast to those of the matrix, the robust Newton-Raphson iteration method is implemented according to the efficient procedure of Kropinski (1990) and Forsyth and Simpson (1991) to solve the variably-saturated flow equation. Upstream weighting of the relative permeabilities is used to yield a monotone solution that lies in the physical range and adaptive time stepping further enhances the efficiency of the solution process. A time-marching Galerkin finite-element technique is used to discretize the solute transport equation. Although the methodology is developed in a finite-element framework, a finite-difference discretization for both groundwater flow and solute transport can be mimicked through a manipulation of the influence coefficient technique. The use of an ILU-preconditioned ORTHOMIN solver permits the fast solution of matrix equations having tens to hundreds of thousands of unknowns. Verification examples are presented along with illustrative problems that demonstrate the complexity of variably-saturated flow and solute transport in fractured systems.

Migration of contaminants in groundwater at a landfill: A case study: 4. A natural-gradient dispersion test

Abstract A natural-gradient tracer test using a chloride solution with an initial injection volume of 0.7 m 3 was performed in the sandy aquifer at the Borden site. The solution was injected into five well points set ∼1 m below the water table in an uncontaminated zone situated above the contaminant plume at a location ∼450 m downflow from the landfill. Under conditions of natural groundwater flow, the tracer slug was monitored for a period of 4 months by withdrawing small-volume samples from points in a three-dimensional array of bundle-type multilevel samplers. Measurements of hydraulic head were obtained from a network of miniature piezometers. Soon after injection, the tracer slug gradually split into two halves, one half moving horizontally at an average velocity of 2.9 · 10 −6 ms −1 and the other horizontally at 8.2 · 10 −7 ms −1 . Although the split has been attributed to local lateral heterogeneity, the nature of the heterogeneity and its influence on the hydraulic-head distribution were not clearly distinguishable in the field data obtained before, during or after the test. The chloride patterns for each of the two halves of the tracer slug evolved into Gaussian forms although the patterns demonstrated some irregularity at early time. The relatively smooth Gaussian forms were unexpected because the aquifer has numerous small-scale heterogeneities observed in vertical cores obtained from the tracer zone and because the glaciofluvial origin of the aquifer suggests that heterogeneities are not laterally continuous. Simulated chloride distributions from a three-dimensional analytical solution to the advection-dispersion equation were fitted to the field data to obtain best-fit estimates of the values of longitudinal, transverse-horizontal and transverse-vertical dispersivity at various travel distances for each of the two halves of the tracer zone. This is the first known field test that has permitted the estimation of three principal dispersion coefficients in layered media. The longitudinal dispersivity was found to increase from 0.01 m at a distance of 0.75 m from the source to 0.08 m at 11.0 m. The transverse-horizontal dispersivity increased also to a value of 0.03 m at 11.0 m. Transverse-vertical dispersion was very weak and was accounted for by molecular diffusion. The relative lack of vertical dispersion is consistent with the shape of the plume of leachate contamination from the landfill. It was concluded that the observed increase in dispersivity along the path of migration is likely caused by heterogeneities. Information on the dispersion-controlling heterogeneities is not yet available as practical field methodologies for their identification and description have not yet been developed. Until such information is incorporated into mass-transport models, a realistic solution of the dispersion problem in heterogeneous media will remain elusive.

Surface‐subsurface model intercomparison: A first set of benchmark results to diagnose integrated hydrology and feedbacks

There are a growing number of large-scale, complex hydrologic models that are capable of simulating integrated surface and subsurface flow. Many are coupled to land-surface energy balance models, biogeochemical and ecological process models, and atmospheric models. Although they are being increasingly applied for hydrologic prediction and environmental understanding, very little formal verification and/or benchmarking of these models has been performed. Here we present the results of an intercomparison study of seven coupled surface-subsurface models based on a series of benchmark problems. All the models simultaneously solve adapted forms of the Richards and shallow water equations, based on fully 3-D or mixed (1-D vadose zone and 2-D groundwater) formulations for subsurface flow and 1-D (rill flow) or 2-D (sheet flow) conceptualizations for surface routing. A range of approaches is used for the solution of the coupled equations, including global implicit, sequential iterative, and asynchronous linking, and various strategies are used to enforce flux and pressure continuity at the surface-subsurface interface. The simulation results show good agreement for the simpler test cases, while the more complicated test cases bring out some of the differences in physical process representations and numerical solution approaches between the models. Benchmarks with more traditional runoff generating mechanisms, such as excess infiltration and saturation, demonstrate more agreement between models, while benchmarks with heterogeneity and complex water table dynamics highlight differences in model formulation. In general, all the models demonstrate the same qualitative behavior, thus building confidence in their use for hydrologic applications.

Comparative error analysis in finite element formulations of the advection-dispersion equation

Abstract The behaviour of numerical solutions of the one-dimensional advection-dispersion equation is investigated and comparisons between the consistent and the lumped formulations of Galerkin finite element schemes are made. Well-known criteria for the control of accuracy in the lumped (finite difference) formulation are reviewed. It is found that, because the numerical error produced by the consistent formulation is generally less than that produced by the lumped formulation, these criteria can also be used for the control of numerical dispersion in the consistent formulation. However, because the error in both types of solutions decreases in time when the discretization is invariant, the criteria can be relaxed with advancing simulation time. For the consistent formulation it is found that beyond some initial time period, the numerical error depends only on the temporal discretization. This suggests that constant accuracy can be maintained throughout the simulation period while allowing the time step length to grow.

Dynamics of groundwater recharge and seepage over the Canadian landscape during the Wisconsinian glaciation

Received 30 May 2007; accepted 12 October 2007; published 15 February 2008. [1] Pleistocene glaciations and their associated dramatic climatic conditions are suspected to have had a large impact on the groundwater flow system over the entire North American continent. Because of the myriad of complex flow-related processes involved during a glaciation period, numerical models have become powerful tools for examining groundwater flow system evolution in this context. In this study, a series of key processes pertaining to coupled groundwater flow and glaciation modeling, such as densitydependent (i.e., brine) flow, hydromechanical loading, subglacial infiltration, isostasy, and permafrost development, are included in the numerical model HydroGeoSphere to simulate groundwater flow over the Canadian landscape during the Wisconsinian glaciation (�� 120 ka to present). The primary objective is to demonstrate the immense impact of glacial advances and retreats during the Wisconsinian glaciation on the dynamical evolution of groundwater flow systems over the Canadian landscape, including surface-subsurface water exchanges (i.e., recharge and discharge fluxes) in both the subglacial and the periglacial environments. It is shown that much of the infiltration of subglacial meltwater occurs during ice sheet progression and that during ice sheet regression, groundwater mainly exfiltrates on the surface, in both the subglacial and periglacial environments. The average infiltration/exfiltration fluxes range between 0 and 12 mm/a. Using mixed, ice sheet thickness–dependent boundary conditions for the subglacial environment, it was estimated that 15–70% of the meltwater infiltrated into the subsurface as recharge, with an average of 43%. Considering the volume of meltwater that was generated subsequent to the last glacial maximum, these recharge rates, which are related to the bedrock type and elastic properties, are historically significant and therefore played an immense role in the evolution of groundwater flow system evolution over the Canadian landmass over the last 120 ka. Finally, it is shown that the permafrost extent plays a key role in the distribution of surface-subsurface interaction because the presence of permafrost acts as a barrier for groundwater flow.

Three-dimensional plume source reconstruction using minimum relative entropy inversion

In this paper we extend the minimum relative entropy (MRE) method to recover the source-release history of a three dimensional plume. This extension is carried out in an analytic framework, and in order to qualify as a linear inverse problem the various transport parameters such as dispersivity and the like are considered to be known. In addition, the groundwater flow system is assumed to be steady and uniform. The contributions of this paper include an explanation how MRE can be used as a measure of resolution in linear inversion, a reporting of a three dimensional analytic solution for mass transport in a steady one dimensional velocity field for a variable in-time source loading, an estimation the source-release history for synthetically generated data sets, and an application of the methodology to a case-study problem at the Gloucester Landfill in Ontario, Canada. We found that the relative entropy measure is useful in indicating the reduction in uncertainty between the posterior and prior pdfs as a result of the new information provided by the physical constraints and data. Using the individual model-parameter relative entropies as a measure of resolution, one can make quantitative judgments about which part of the history is likely to be well resolved. Comparing inversion results for synthetic aquifers with one well and two sample points with only one sample point indicates that temporal data at several wells allows for a superior reconstruction of the release history. We investigate the potential benefits of locating sample points on an spatial rather than temporal basis. Results show that early part of the release history is poorly recovered. Comparing these results with one well and two sample points indicates that temporal data at a few wells allows for a better reconstruction of the release history. An incomplete time record is also investigated. Results show that early part of the release history prior to the commencement of measurements is poorly recovered. It is essential that as much as possible of the time histories of plumes be monitored if the entire release history is to be determined. The MRE approach is used to reconstruct the release history of a 1,4-dioxane plume measured at the Gloucester Landfill in Ontario, Canada. The recovered release history is fairly narrow and generally flat in shape, although two peaks are evident. Neither peak is particularly well defined, judging from the resolution curve and analysis of the confidence ranges. One of the peaks coincides with year 1979, close to the year 1978 in which a large spill was noted. However, a considerable variation of possible source releases is possible.

Multicomponent simulation of wastewater-derived nitrogen and carbon in shallow unconfined aquifers. I. Model formulation and performance

One of the most common methods to dispose of domestic wastewater involves the release of septic effluent from drains located in the unsaturated zone. Nitrogen from such systems is currently of concern because of nitrate contamination of drinking water supplies and eutrophication of coastal waters. The objectives of this study are to develop and assess the performance of a mechanistic flow and reactive transport model which couples the most relevant physical, geochemical and biochemical processes involved in wastewater plume evolution in sandy aquifers. The numerical model solves for variably saturated groundwater flow and reactive transport of multiple carbon- and nitrogen-containing species in a three-dimensional porous medium. The reactive transport equations are solved using the Strang splitting method which is shown to be accurate for Monod and first- and second-order kinetic reactions, and two to four times more efficient than sequential iterative splitting. The reaction system is formulated as a fully kinetic chemistry problem, which allows for the use of several special-purpose ordinary differential equation (ODE) solvers. For reaction systems containing both fast and slow kinetic reactions, such as the combined nitrogen-carbon system, it is found that a specialized stiff explicit solver fails to obtain a solution. An implicit solver is more robust and its computational performance is improved by scaling of the fastest reaction rates. The model is used to simulate wastewater migration in a 1-m-long unsaturated column and the results show significant oxidation of dissolved organic carbon (DOC), the generation of nitrate by nitrification, and a slight decrease in pH.

Numerical simulation of multiphase flow and phase partitioning in discretely fractured geologic media

Abstract The three-dimensional compositional model CompFlow has been extended to allow the simulation of the multiphase advective, dispersive and diffusive flux of non-aqueous phase liquid (NAPL) contaminants in a discrete-fracture network, allowing for phase partitioning and dynamic interactions between the fracture network and the surrounding low-permeability rock matrix. The approach used to couple fluxes between the fractures and matrix allows representation of capillary pressure differences within the fractures and matrix and makes no assumption of equilibrium hydraulic conditions between the two. The model is verified for the case of aqueous-phase solute transport by comparison with an analytical solution. An example problem is presented involving the migration of a dense non-aqueous phase liquid (DNAPL) consisting of trichlorethylene (TCE) in a single vertical fracture within a low-permeability material with significant matrix porosity. The simulation results demonstrate that matrix diffusion acts to transfer significant amounts of contaminant to the matrix in the aqueous phase. After the DNAPL source is removed, the NAPL ultimately disappears from the fracture due to partitioning of contaminant into the aqueous phase with concomitant matrix diffusion. It is shown that as the porosity of the matrix increases, the rate of migration of the TCE DNAPL front within fractures is retarded, due to dissolution and matrix diffusion. The sensitivity of DNAPL migration within the fracture to the form of the relative permeability relationship is also discussed. The model is then used to highlight the potential for deep DNAPL penetration through a vertical cross-section consisting of a shallow unconfined sand aquifer and a deeper sand aquifer separated by a layer of fractured clay. Vertical fractures through the clay that hydraulically connects the shallow and deep aquifers are shown to be capable of transmitting both dissolved contaminant and DNAPL to the underlying aquifer, with DNAPL travel times through the 5-m thick clay unit being on the order of days. For the scenarios examined, the DNAPL entering the lower aquifer via the fractured clay unit may completely dissolve at the aquifer–aquitard interface due to lateral groundwater flow in the lower aquifer, provided the aperture of the vertical fractures in the clay layer are less than about 30 μm; however, some DNAPLs may penetrate into the lower aquifer and exist as a non-aqueous phase if the fracture apertures in the clay layer are 50 μm in size or larger. In this latter case, the presence of the DNAPL in the lower aquifer acts as a persistent source of groundwater contamination and can produce an extensive plume in the direction of groundwater flow.

Density-dependent solute transport in discretely-fractured geologic media : is prediction possible?

Abstract The development of a dense solute plume in a fractured geologic medium can be highly irregular due to both the complexity of the fracture network as well as the presence of convection cells that may arise as a result of the density contrast between the invading solute and the ambient groundwater. A two-dimensional numerical model has been developed here to investigate density-dependent groundwater flow and solute transport in geologic materials that contain discrete fractures in order to examine some of the complex forms into which plumes can evolve, particularly with regard to fracture–matrix interactions. Results from simulations which involve parallel vertical fractures show that the evolution of the solute plume is affected by the development of convection cells in the porous matrix blocks between the vertical fractures. In a geologic medium containing a network of regularly spaced horizontal and vertical fractures, complex migration pathways can develop that are unexpected even though the geometry and interconnectivity of the fractures are known a priori. Downward solute migration can occur in some vertical fractures, while upward migration of less dense fluid can occur in others with transient circulation patterns developing in the intervening porous matrix. Because of the inherent uncertainty associated with fracture delineation, and because of the irregular nature of unstable dense plumes, deterministic prediction of dense-plume migration pathways and travel times in fractured geologic media will be subject to considerable uncertainty.