Probabilistic risk assessment of contaminant transport in groundwater and vapour intrusion following remediation of a contaminant source

Community-scale simulations were performed to investigate the risk to groundwater and indoor air receptors downgradient of a contaminated site following the remediation of a long-term source. Six suites of Monte Carlo simulations were performed using a numerical model that accounted for groundwater flow, reactive solute transport, soil gas flow, and vapour intrusion in buildings. The model was applied to a three-dimensional, community-scale (250 m × 1000 m × 14 m) domain containing heterogeneous, spatially correlated distributions of the hydraulic conductivity, fraction of organic carbon, and biodegradation rate constant, which were varied between realizations. Analysis considered results from both individual realizations as well as the suite of Monte Carlo simulations expressed through several novel, integrated parameters, such as the probability of exceeding a regulatory standard in either groundwater or indoor air. Results showed that exceedance probabilities varied considerably with the consideration of biodegradation in the saturated zone, and were less sensitive to changes in the variance of hydraulic conductivity or the incorporation of heterogeneous distributions of organic carbon at this spatial scale. A sharp gradient in exceedance probability existed at the lateral edges of the plumes due to variability in lateral dispersion, which defined a narrow region of exceedance uncertainty. Differences in exceedance probability between realizations (i.e., due to heterogeneity uncertainty) were similar to differences attributed to changes in the variance of hydraulic conductivity or fraction of organic carbon. Simulated clean-up times, defined by reaching an acceptable exceedance probability, were found to be on the order of decades to centuries in these community-scale domains. Results also showed that the choice of the acceptable exceedance probability level (e.g., 1 vs. 5 %) would likely affect clean up times on the order of decades. Moreover, in the scenarios examined here, the risk of exceeding indoor air standards was greater than that of exceeding groundwater standards at all times and places. Overall, simulations of coupled transport processes combined with novel spatial and temporal quantification metrics for Monte Carlo analyses, provide practical tools for assessing risk in wider communities when considering site remediation.     

Efficient uncertainty quantification techniques in inverse problems for Richards’ equation using coarse-scale simulation models

This paper concerns efficient uncertainty quantification techniques in inverse problems for Richards’ equation which use coarse-scale simulation models. We consider the problem of determining saturated hydraulic conductivity fields conditioned to some integrated response. We use a stochastic parameterization of the saturated hydraulic conductivity and sample using Markov chain Monte Carlo methods (MCMC). The main advantage of the method presented in this paper is the use of multiscale methods within an MCMC method based on Langevin diffusion. Additionally, we discuss techniques to combine multiscale methods with stochastic solution techniques, specifically sparse grid collocation methods. We show that the proposed algorithms dramatically reduce the computational cost associated with traditional Langevin MCMC methods while providing similar sampling performance.     

Polynomial chaos expansion for global sensitivity analysis applied to a model of radionuclide migration in a randomly heterogeneous aquifer

We perform global sensitivity analysis (GSA) through polynomial chaos expansion (PCE) on a contaminant transport model for the assessment of radionuclide concentration at a given control location in a heterogeneous aquifer, following a release from a near surface repository of radioactive waste. The aquifer hydraulic conductivity is modeled as a stationary stochastic process in space. We examine the uncertainty in the first two (ensemble) moments of the peak concentration, as a consequence of incomplete knowledge of (a) the parameters characterizing the variogram of hydraulic conductivity, (b) the partition coefficient associated with the migrating radionuclide, and (c) dispersivity parameters at the scale of interest. These quantities are treated as random variables and a variance-based GSA is performed in a numerical Monte Carlo framework. This entails solving groundwater flow and transport processes within an ensemble of hydraulic conductivity realizations generated upon sampling the space of the considered random variables. The Sobol indices are adopted as sensitivity measures to provide an estimate of the role of uncertain parameters on the (ensemble) target moments. Calculation of the indices is performed by employing PCE as a surrogate model of the migration process to reduce the computational burden. We show that the proposed methodology (a) allows identifying the influence of uncertain parameters on key statistical moments of the peak concentration (b) enables extending the number of Monte Carlo iterations to attain convergence of the (ensemble) target moments, and (c) leads to considerable saving of computational time while keeping acceptable accuracy.     

Screening‐Level Sensitivity Analysis for the Design of Pump‐and‐Treat Systems

Pump‐and‐treat systems can prevent the migration of groundwater contaminants and candidate systems are typically evaluated with groundwater models. Such models should be rigorously assessed to determine predictive capabilities and numerous tools and techniques for model assessment are available. While various assessment methodologies (e.g., model calibration, uncertainty analysis, and Bayesian inference) are well‐established for groundwater modeling, this paper calls attention to an alternative assessment technique known as screening‐level sensitivity analysis (SLSA). SLSA can quickly quantify first‐order (i.e., main effects) measures of parameter influence in connection with various model outputs. Subsequent comparisons of parameter influence with respect to calibration vs. prediction outputs can suggest gaps in model structure and/or data. Thus, while SLSA has received little attention in the context of groundwater modeling and remedial system design, it can nonetheless serve as a useful and computationally efficient tool for preliminary model assessment. To illustrate the use of SLSA in the context of designing groundwater remediation systems, four SLSA techniques were applied to a hypothetical, yet realistic, pump‐and‐treat case study to determine the relative influence of six hydraulic conductivity parameters. Considered methods were: Taguchi design‐of‐experiments (TDOE); Monte Carlo statistical independence (MCSI) tests; average composite scaled sensitivities (ACSS); and elementary effects sensitivity analysis (EESA). In terms of performance, the various methods identified the same parameters as being the most influential for a given simulation output. Furthermore, results indicate that the background hydraulic conductivity is important for predicting system performance, but calibration outputs are insensitive to this parameter (K_BK). The observed insensitivity is attributed to a nonphysical specified‐head boundary condition used in the model formulation which effectively “staples” head values located within the conductivity zone. Thus, potential strategies for improving model predictive capabilities include additional data collection targeting the K_BK parameter and/or revision of model structure to reduce the influence of the specified head boundary.     

Field Study of Subsurface Heterogeneity with Steady‐State Hydraulic Tomography

Remediation of subsurface contamination requires an understanding of the contaminant (history, source location, plume extent and concentration, etc.), and, knowledge of the spatial distribution of hydraulic conductivity (K) that governs groundwater flow and solute transport. Many methods exist for characterizing K heterogeneity, but most if not all methods require the collection of a large number of small‐scale data and its interpolation. In this study, we conduct a hydraulic tomography survey at a highly heterogeneous glaciofluvial deposit at the North Campus Research Site (NCRS) located at the University of Waterloo, Waterloo, Ontario, Canada to sequentially interpret four pumping tests using the steady‐state form of the Sequential Successive Linear Estimator (SSLE) ( Yeh and Liu 2000 ). The resulting three‐dimensional (3D) K distribution (or K‐tomogram) is compared against: ( 1 ) K distributions obtained through the inverse modeling of individual pumping tests using SSLE, and ( 2 ) effective hydraulic conductivity (K_eff) estimates obtained by automatically calibrating a groundwater flow model while treating the medium to be homogeneous. Such a K_eff is often used for designing remediation operations, and thus is used as the basis for comparison with the K‐tomogram. Our results clearly show that hydraulic tomography is superior to the inversions of single pumping tests or K_eff estimates. This is particularly significant for contaminated sites where an accurate representation of the flow field is critical for simulating contaminant transport and injection of chemical and biological agents used for active remediation of contaminant source zones and plumes.     

Bayesian approach to contaminant source characterization in water distribution systems: adaptive sampling framework

Bayesian analysis can yield a probabilistic contaminant source characterization conditioned on available sensor data and accounting for system stochastic processes. This paper is based on a previously proposed Markov chain Monte Carlo (MCMC) approach tailored for water distribution systems and incorporating stochastic water demands. The observations can include those from fixed sensors and, the focus of this paper, mobile sensors. Decision makers, such as utility managers, need not wait until new observations are available from an existing sparse network of fixed sensors. This paper addresses a key research question: where is the best location in the network to gather additional measurements so as to maximize the reduction in the source uncertainty? Although this has been done in groundwater management, it has not been well addressed in water distribution networks. In this study, an adaptive framework is proposed to guide the strategic placement of mobile sensors to complement the fixed sensor network. MCMC is the core component of the proposed adaptive framework, while several other pieces are indispensable: Bayesian preposterior analysis, value of information criterion and the search strategy for identifying an optimal location. Such a framework is demonstrated with an illustrative example, where four candidate sampling locations in the small water distribution network are investigated. Use of different value-of-information criteria reveals that while each may lead to different outcomes, they share some common characteristics. The results demonstrate the potential of Bayesian analysis and the MCMC method for contaminant event management.     

Stochastic analysis of estuarine hydraulics

A new methodology is presented for the solution of the stochastic hydraulic equations characterizing steady, one-dimensional estuarine flow. The methodology is predicated on quasi-linearization, perturbation methods, and the finite difference approximation of the stochastic differential operators. Assuming Manning's roughness coefficient is the principal source of uncertainty in the model, stochastic equations are presented for the water depths and flow rates in the estuarine system. Moment equations are developed for the mean and variance of the water depths. The moment equations are compared with the results of Monte Carlo simulation experiments. The results confirm that for any spatial location in the estuary that (1) as the uncertainty in the channel roughness increases, the uncertainty in mean depth increases, and (2) the predicted mean depth will decrease with increasing uncertainty in Manning'sn. The quasi-analytical approach requires significantly less computer time than Monte Carlo simulations and provides explicit     

Delineating a recharge area for a spring using numerical modeling, Monte Carlo techniques, and geochemical investigation

Recharge areas of spring systems can be hard to identify, but they can be critically important for protection of a spring resource. A recharge area for a spring complex in southern Wisconsin was delineated using a variety of complementary techniques. A telescopic mesh refinement (TMR) model was constructed from an existing regional-scale ground water flow model. This TMR model was formally optimized using parameter estimation techniques; the optimized "best fit" to measured heads and fluxes was obtained by using a horizontal hydraulic conductivity 200% larger than the original regional model for the upper bedrock aquifer and 80% smaller for the lower bedrock aquifer. The uncertainty in hydraulic conductivity was formally considered using a stochastic Monte Carlo approach. Two-hundred model runs used uniformly distributed, randomly sampled, horizontal hydraulic conductivity values within the range given by the TMR optimized values and the previously constructed regional model. A probability distribution of particles captured by the spring, or a "probabilistic capture zone," was calculated from the realistic Monte Carlo results (136 runs of 200). In addition to portions of the local surface watershed, the capture zone encompassed areas outside of the watershed--demonstrating that the ground watershed and surface watershed do not coincide. Analysis of water collected from the site identified relatively large contrasts in chemistry, even for springs within 15 m of one another. The differences showed a distinct gradation from Ordovician-carbonate-dominated water in western spring vents to Cambrian-sandstone-influenced water in eastern spring vents. The difference in chemistry was attributed to distinctive bedrock geology as demonstrated by overlaying the capture zone derived from numerical modeling over a bedrock geology map for the area. This finding gives additional confidence to the capture zone calculated by modeling.     

A Model-Averaging Method for Assessing Groundwater Conceptual Model Uncertainty

This study evaluates alternative groundwater models with different recharge and geologic components at the northern Yucca Flat area of the Death Valley Regional Flow System (DVRFS), USA. Recharge over the DVRFS has been estimated using five methods, and five geological interpretations are available at the northern Yucca Flat area. Combining the recharge and geological components together with additional modeling components that represent other hydrogeological conditions yields a total of 25 groundwater flow models. As all the models are plausible given available data and information, evaluating model uncertainty becomes inevitable. On the other hand, hydraulic parameters (e.g., hydraulic conductivity) are uncertain in each model, giving rise to parametric uncertainty. Propagation of the uncertainty in the models and model parameters through groundwater modeling causes predictive uncertainty in model predictions (e.g., hydraulic head and flow). Parametric uncertainty within each model is assessed using Monte Carlo simulation, and model uncertainty is evaluated using the model averaging method. Two model-averaging techniques (on the basis of information criteria and GLUE) are discussed. This study shows that contribution of model uncertainty to predictive uncertainty is significantly larger than that of parametric uncertainty. For the recharge and geological components, uncertainty in the geological interpretations has more significant effect on model predictions than uncertainty in the recharge estimates. In addition, weighted residuals vary more for the different geological models than for different recharge models. Most of the calibrated observations are not important for discriminating between the alternative models, because their weighted residuals vary only slightly from one model to another.     

Probabilistic data integration for characterization of spatial distribution of residual LNAPL

The spatial distribution of residual light non-aqueous phase liquid (LNAPL) is an important factor in reactive solute transport modeling studies. There is great uncertainty associated with both the areal limits of LNAPL source zones and smaller scale variability within the areal limits. A statistical approach is proposed to construct a probabilistic model for the spatial distribution of residual NAPL and it is applied to a site characterized by ultra-violet-induced-cone-penetration testing (CPT–UVIF). The uncertainty in areal limits is explicitly addressed by a novel distance function (DF) approach. In modeling the small-scale variability within the areal limits, the CPT–UVIF data are used as primary source of information, while soil texture and distance to water table are treated as secondary data. Two widely used geostatistical techniques are applied for the data integration, namely sequential indicator simulation with locally varying means (SIS–LVM) and Bayesian updating (BU). A close match between the calibrated uncertainty band (UB) and the target probabilities shows the performance of the proposed DF technique in characterization of uncertainty in the areal limits. A cross-validation study also shows that the integration of the secondary data sources substantially improves the prediction of contaminated and uncontaminated locations and that the SIS–LVM algorithm gives a more accurate prediction of residual NAPL contamination. The proposed DF approach is useful in modeling the areal limits of the non-stationary continuous or categorical random variables, and in providing a prior probability map for source zone sizes to be used in Monte Carlo simulations of contaminant transport or Monte Carlo type inverse modeling studies.     

Groundwater Modeling with Nonlinear Uncertainty Analyses to Enhance Remediation Design Confidence

Soil and groundwater contamination are often managed by establishing on‐site cleanup targets within the context of risk assessment or risk management measures. Decision‐makers rely on modeling tools to provide insight; however, it is recognized that all models are subject to uncertainty. This case study compares suggested remediation requirements using a site‐specific numerical model and a standardized analytical tool to evaluate risk to a downgradient wetland receptor posed by on‐site chloride impacts. The base case model, calibrated to observed non‐pumping and pumping conditions, predicts a peak concentration well above regulatory criteria. Remediation scenarios are iteratively evaluated to determine a remediation design that adheres to practical site constraints, while minimizing the potential for risk to the downgradient receptor. A nonlinear uncertainty analysis is applied to each remediation scenario to stochastically evaluate the risk and find a solution that meets the site‐owner risk tolerance, which in this case required a risk‐averse solution. This approach, which couples nonlinear uncertainty analysis with a site‐specific numerical model provides an enhanced level of knowledge to foster informed decision‐making (i.e., risk‐of‐success) and also increases stakeholder confidence in the remediation design.     

Estimating the role of a ditch network in groundwater recharge in a Mediterranean catchment using a water balance approach

Water balance variables were monitored in a farmed Mediterranean catchment characterized by a dense ditch network to allow for the separate estimation of the diffuse and concentrated recharge terms during flood events. The 27 ha central part of the catchment was equipped with (i) rain gauges, (ii) ditch gauge stations, (iii) piezometers, (iv) neutron probes, and (v) an eddy covariance mast including a 3D sonic anemometer and a fast hygrometer. The water balance was calculated for two autumnal rain and flood events. We also estimated the uncertainty of this approach with Monte Carlo simulations. Results show, that although ditch area represents only 6% of the total study area, concentrated recharge appeared to be the main source of groundwater recharge. Indeed, it was 40–50% of the total groundwater recharge for autumnal events, which are the major annual recharge events. This indicate that both, concentrated and diffuse recharge should be taken into account in any hydrological modeling approach for Mediterranean catchments. This also means that, since they collect overland flow that is often largely contaminated by chemicals, ditches may be a place where groundwater contamination is likely to occur. The uncertainty analysis indicates that recharge estimates based on water balance exhibit large uncertainty ranges. Nevertheless, Monte Carlo simulations showed that concentrated recharge was higher than expected based on their area.     

Groundwater–surface water interaction in Lake Nasser,Southern Egypt

A cross‐sectional model, based on the two dimensional groundwater flow equation of Edelman, was applied at seven transects distributed over four geological cross sections to estimate groundwater heads and recharge from/or groundwater discharge to Lake Nasser. The lake with a length of 500 km and an average width of 12 km was created over the period 1964–1970, the time for constructing the Aswan High Dam (AHD). The model, constrained by regional‐scale groundwater flow and groundwater head data in the vicinity of the lake, was successfully calibrated to timeseries of piezometeric heads collected at the cross sections in the period 1965–2004. Inverse modeling yielded high values for the horizontal hydraulic conductivity in the range of 6.0 to 31.1 m day^?1 and storage coefficient between 0.01 and 0.40. The results showed the existence of a strong vertical anisotropy of the aquifer. The calibrated horizontal permeability is systematically higher than the vertical permeability (≈1000:1). The calibrated model was used to explore the recharge from/or groundwater discharge to Lake Nasser at the seven transects for a 40‐year period, i.e. from 1965 to 2004. The analysis for the last 20‐year period, 1985–2004, revealed that recharge from Lake Nasser reduced by 37% compared to the estimates for the first 20‐year period, 1965–1984. In the period 1965–2004, seepage of Lake Nasser to the surrounding was estimated at 1.15 × 10⁹ m³ year^?1. This led to a significant rise of the groundwater table. Variance‐based sensitivity and uncertainty analysis on the Edelman results were conducted applying quasi‐Monte Carlo sequences (Latin Hypercube sampling). The maximum standard deviation of the total uncertainty on the groundwater table was 0.88 m at Toshka (west of the lake). The distance from the lake, followed by the storage coefficient and hydraulic conductivity, were identified as the most sensitive parameters. Copyright © 2012 John Wiley & Sons, Ltd.     

Using data assimilation method to calibrate a heterogeneous conductivity field and improve solute transport prediction with an unknown contamination source

Hydraulic conductivity distribution and plume initial source condition are two important factors affecting solute transport in heterogeneous media. Since hydraulic conductivity can only be measured at limited locations in a field, its spatial distribution in a complex heterogeneous medium is generally uncertain. In many groundwater contamination sites, transport initial conditions are generally unknown, as plume distributions are available only after the contaminations occurred. In this study, a data assimilation method is developed for calibrating a hydraulic conductivity field and improving solute transport prediction with unknown initial solute source condition. Ensemble Kalman filter (EnKF) is used to update the model parameter (i.e., hydraulic conductivity) and state variables (hydraulic head and solute concentration), when data are available. Two-dimensional numerical experiments are designed to assess the performance of the EnKF method on data assimilation for solute transport prediction. The study results indicate that the EnKF method can significantly improve the estimation of the hydraulic conductivity distribution and solute transport prediction by assimilating hydraulic head measurements with a known solute initial condition. When solute source is unknown, solute prediction by assimilating continuous measurements of solute concentration at a few points in the plume well captures the plume evolution downstream of the measurement points.     

Global sensitivity analysis for a numerical model of radionuclide migration from the RRC “Kurchatov Institute” radwaste disposal site

Today, in different countries, there exist sites with contaminated groundwater formed as a result of inappropriate handling or disposal of hazardous materials or wastes. Numerical modeling of such sites is an important tool for a correct prediction of contamination plume spreading and an assessment of environmental risks associated with the site. Many uncertainties are associated with a part of the parameters and the initial conditions of such environmental numerical models. Statistical techniques are useful to deal with these uncertainties. This paper describes the methods of uncertainty propagation and global sensitivity analysis that are applied to a numerical model of radionuclide migration in a sandy aquifer in the area of the RRC “Kurchatov Institute” radwaste disposal site in Moscow, Russia. We consider 20 uncertain input parameters of the model and 20 output variables (contaminant concentration in the observation wells predicted by the model for the end of 2010). Monte Carlo simulations allow calculating uncertainty in the output values and analyzing the linearity and the monotony of the relations between input and output variables. For the non monotonic relations, sensitivity analyses are classically done with the Sobol sensitivity indices. The originality of this study is the use of modern surrogate models (called response surfaces), the boosting regression trees, constructed for each output variable, to calculate the Sobol indices by the Monte Carlo method. It is thus shown that the most influential parameters of the model are distribution coefficients and infiltration rate in the zone of strong pipe leaks on the site. Improvement of these parameters would considerably reduce the model prediction uncertainty.     

Bayesian model averaging assessment on groundwater management under model structure uncertainty

This study introduces Bayesian model averaging (BMA) to deal with model structure uncertainty in groundwater management decisions. A robust optimized policy should take into account model parameter uncertainty as well as uncertainty in imprecise model structure. Due to a limited amount of groundwater head data and hydraulic conductivity data, multiple simulation models are developed based on different head boundary condition values and semivariogram models of hydraulic conductivity. Instead of selecting the best simulation model, a variance-window-based BMA method is introduced to the management model to utilize all simulation models to predict chloride concentration. Given different semivariogram models, the spatially correlated hydraulic conductivity distributions are estimated by the generalized parameterization (GP) method that combines the Voronoi zones and the ordinary kriging (OK) estimates. The model weights of BMA are estimated by the Bayesian information criterion (BIC) and the variance window in the maximum likelihood estimation. The simulation models are then weighted to predict chloride concentrations within the constraints of the management model. The methodology is implemented to manage saltwater intrusion in the “1,500-foot” sand aquifer in the Baton Rouge area, Louisiana. The management model aims to obtain optimal joint operations of the hydraulic barrier system and the saltwater extraction system to mitigate saltwater intrusion. A genetic algorithm (GA) is used to obtain the optimal injection and extraction policies. Using the BMA predictions, higher injection rates and pumping rates are needed to cover more constraint violations, which do not occur if a single best model is used.     

A stochastic optimization framework for the restoration of an over-exploited aquifer

ABSTRACT

This study investigates the impact of hydraulic conductivity uncertainty on the sustainable management of the aquifer of Lake Karla, Greece, using the stochastic optimization approach. The lack of surface water resources in combination with the sharp increase in irrigation needs in the basin over the last 30 years have led to an unprecedented degradation of the aquifer. In addition, the lack of data regarding hydraulic conductivity in a heterogeneous aquifer leads to hydrogeologic uncertainty. This uncertainty has to be taken into consideration when developing the optimization procedure in order to achieve the aquifer’s sustainable management. Multiple Monte Carlo realizations of this spatially-distributed parameter are generated and groundwater flow is simulated for each one of them. The main goal of the sustainable management of the ‘depleted’ aquifer of Lake Karla is two-fold: to determine the optimum volume of renewable groundwater that can be extracted, while, at the same time, restoring its water table to a historic high level. A stochastic optimization problem is therefore formulated, based on the application of the optimization method for each of the aquifer’s multiple stochastic realizations in a future period. In order to carry out this stochastic optimization procedure, a modelling system consisting of a series of interlinked models was developed. The results show that the proposed stochastic optimization framework can be a very useful tool for estimating the impact of hydraulic conductivity uncertainty on the management strategies of a depleted aquifer restoration. They also prove that the optimization process is affected more by hydraulic conductivity uncertainty than the simulation process.

Editor Z.W. Kundzewicz; Guest editor S. Weijs     

Assessment of uncertainty associated with the estimation of well catchments by moment equations

Non-local stochastic moment equations are used successfully to analyze groundwater flow in randomly heterogeneous media. Here we present a moment equations-based approach to quantify the uncertainty associated with the estimation of well catchments. Our approach is based on the development of a complete second order formalism which allows obtaining the first statistical moments of the trajectories of conservative solute particles advected in a generally non-uniform groundwater flow. Approximate equations of moments of particles’ trajectories are then derived on the basis of a second order expansion in terms of the standard deviation of the aquifer log hydraulic conductivity. Analytical expressions are then obtained for the predictors of locations of mean stagnation points, together with their associated uncertainties. We implement our approach on heterogeneous media in bounded two-dimensional domains, with and without including the effect of conditioning on hydraulic conductivity information. The impact of domain size, boundary conditions, heterogeneity and non-stationarity of hydraulic conductivity on the prediction of a well catchment is explored. The results are compared against Monte Carlo simulations and semi-analytical solutions available in the literature. The methodology is applicable to both infinite and bounded domains and is free of distributional assumptions (and so applies to both Gaussian and non-Gaussian log hydraulic conductivity fields) and formally includes the effect of conditioning on available information.