Semi-analytical solutions for solute transport in fractured porous media using a strip source of finite width

Transient and steady-state analytical solutions are derived to investigate solute transport in a fractured porous medium consisting of evenly spaced, parallel discrete fractures. The solutions incorporate a finite width strip source, longitudinal and transverse dispersion in the fractures, source decay, aqueous phase decay, one-dimensional diffusion into the matrix, sorption to fracture walls, and sorption within the matrix. The solutions are derived using Laplace and Fourier transforms, and inverted by interchanging the order of integration and utilizing a numerical Laplace inversion algorithm. The solutions are verified for simplified cases by comparison to solutions derived by Batu [Batu V. A generalized two-dimensional analytical solution for hydrodynamic dispersion in bounded media with the first-type condition at the source. Wat Resour Res 1989;25(6):1125] and Sudicky and Frind [Sudicky EA, Frind EO. Contaminant transport in fractured porous media: analytical solutions for a system of parallel fractures. Wat Resour Res 1982;18(6):1634]. The application of the solutions to a fractured sandstone demonstrates that narrower source widths and larger values of transverse dispersivity both lead to lower downstream concentrations in the fractures and shorter steady-state plumes. The incorporation of aqueous phase decay and source concentration decay both lead to lower concentrations and shorter plumes, with even moderate amounts of decay significantly shortening the persistence of contamination.     

Modeling plume behavior for nonlinearly sorbing solutes in saturated homogeneous porous media

Transport of a sorbing solute in a two-dimensional steady and uniform flow field is modeled using a particle tracking random walk method. The solute is initially introduced from an instantaneous point source. Cases of linear and nonlinear sorption isotherms are considered. Local pore velocity and mechanical dispersion are used to describe the solute transport mechanisms at the local scale. The numerical simulation of solute particle transport yields the large scale behavior of the solute plume. Behavior of the plume is quantified in terms of the center-of-mass displacement distance, relative velocity of the center-of-mass, mass breakthrough curves, spread variance, and longitudinal skewness. The nonlinear sorption isotherm affects the plume behavior in the following way relative to the linear isotherm: (1) the plume velocity decreases exponentially with time; (2) the longitudinal variance increases nonlinearly with time; (3) the solute front is steepened and tailing is enhanced     

An Improved Method for Estimation of Biodegradation Rate with Field Data

The method of Buscheck and Alcantar (1995) is shown to be valid for the concentration data observed in wells along the centerline and near a source of relatively large area as compared with the lateral dispersion of the plume. When such a condition is not satisfied, the method of Buscheck and Alcantar (1995) overestimates the biodegradation rate due to its usage of a solution to a one-dimensional advection-dispersion equation that does not consider lateral dispersion and finite source size. It is also shown that the concentration of a two- or three-dimensional plume does not always decrease exponentially with the distance from the source along a plume centerline. An improved method based on multidimensional solutions of Domenico (1987) is proposed. In using the improved method, observed concentrations are normalized by the error functions embedded in the solutions to two- or three-dimensional advection-dispersion equations for a continuous source of finite size. With the normalized concentrations and following the procedure of Buscheck and Alcantar (1995), one may estimate the biodegradation rate, which is usually smaller than that obtained with the original data. It is demonstrated with the field data presented by Buscheck and Alcantar (1995) at the Fairfax Terminal that the method of Buscheck and Alcantar (1995) indeed overestimates the biodegradation rate by 21% and 65% for a two-dimensional plume and for a three-dimensional plume, respectively.     

Variable-density groundwater flow and solute transport in irregular 2D fracture networks

Numerical simulations of variable-density flow and solute transport have been conducted to investigate dense plume migration for various configurations of 2D fracture networks. For orthogonal fractures, simulations demonstrate that dispersive mixing in fractures with small aperture does not stabilize vertical plume migration in fractures with large aperture. Simulations in non-orthogonal 2D fracture networks indicate that convection cells form and that they overlap both the porous matrix and fractures. Thus, transport rates in convection cells depend on matrix and fracture flow properties. A series of simulations in statistically equivalent networks of fractures with irregular orientation show that the migration of a dense plume is highly sensitive to the geometry of the network. If fractures in a random network are connected equidistantly to the solute source, few equidistantly distributed fractures favor density-driven transport. On the other hand, numerous fractures have a stabilizing effect, especially if diffusive transport rates are high. A sensitivity analysis for a network with few equidistantly distributed fractures shows that low fracture aperture, low matrix permeability and high matrix porosity impede density-driven transport because these parameters reduce groundwater flow velocities in both the matrix and the fractures. Enhanced molecular diffusion slows down density-driven transport because it favors solute diffusion from the fractures into the low-permeability porous matrix where groundwater velocities are smaller. For the configurations tested, variable-density flow and solute transport are most sensitive to the permeability and porosity of the matrix, which are properties that can be determined more accurately than the geometry and hydraulic properties of the fracture network, which have a smaller impact on density-driven transport.     

Effective dispersion in conditioned transmissivity fields

We analyze the impact of conditioning to measurements of hydraulic transmissivity on the transport of a conservative solute. The effects of conditioning on solute transport are widely discussed in the literature, but most of the published works focuses on the reduction of the uncertainty in the prediction of the plume dispersion. In this study both ensemble and effective plume moments are considered for an instantaneous release of a solute through a linear source normal to the mean flow direction, by taking into account different sizes of the source. The analysis, involving a steady and spatially inhomogeneous velocity field, is developed by using the stochastic finite element method. Results show that conditioning reduces the ensemble moment in comparison with the unconditioned case, whereas the effective dispersion may increase because of the contribution of the spatial moments related to the lack of stationarity in the flow field. As the number of conditioning points increases, this effect increases and it is significant in both the longitudinal and transverse directions. Furthermore, we conclude that the moment derived from data collected in the field can be assessed by the conditioned second-order spatial moment only with a dense grid of measured data, and it is manifest for larger initial lengths of the plume. Nevertheless, it seems very likely that the actual dispersion of the plume may be underestimated in practical applications.     

The significance of flow in the matrix of the Chalk unsaturated zone

The significance of flow in the matrix of the Chalk unsaturated zone, in comparison with flow in fractures, has been the subject of much debate. In this article, important elements of the literature are discussed in detail and several simple modelling analyses based on steady-state flow are presented. A study of the sensitivity of solute spreading to fracture spacing in models that ignore matrix flow shows that this latter assumption is generally incompatible with observed solute profiles, unless unrealistically small fracture spacings are assumed. The effect of air phase continuities (e.g. bedding planes) on matrix flow has also been examined. These discontinuities are frequently interrupted by points of connectivity between matrix blocks. An issue therefore is the relationship between connectivity and its effect on inter-block conductance. A simple analysis of the Laplace equation shows that just 1% connectivity represents an effective pathway equivalent to 18% of the local rock hydraulic conductivity. Obviously, when there is no fracture flow, solute spreading is significantly reduced. However, dual permeability model simulations show that matrix flow reduces solute spreading in the presence of persistent fracture flow. All of the above studies suggest that flow in the matrix of the Chalk unsaturated zone is significant and that ignoring it may result in a serious misunderstanding of the system.     

Analytical solutions for transport of decaying solutes in rivers with transient storage

Analytical solutions are presented for solute transport in rivers including the effects of transient storage and first order decay. The solute transport model considers an advection–dispersion equation for transport in the main channel linked to a first order mass exchange between the main channel and the transient storage zones. In case of a conservative tracer, it is shown that different analytical solutions presented in the literature are mathematically identical. For non-conservative solutes, first order decay reactions are considered with different reaction rate coefficients in the main river channel and in the dead zones. New analytical solutions are presented for different boundary conditions, i.e. instantaneous injection in an infinite river reach, and variable concentration time series input in a semi-infinite river reach. The correctness and accuracy of the analytical solutions is verified by comparison with the OTIS numerical model. The results of analytical and numerical approaches compare favourably and small differences can be attributed to the influence of boundary conditions. It is concluded that the presented analytical solutions for solute transport in rivers with transient storage and solute decay are accurate and correct, and can be usefully applied for analyses of tracer experiments and transport characteristics in rivers with mass exchange in dead zones.     

Evaluation of volatilization as a natural attenuation pathway for MTBE

Volatilization and diffusion through the unsaturated zone can be an important pathway for natural attenuation remediation of methyl tert-butyl ether (MTBE) at gasoline spill sites. The significance of this pathway depends primarily on the distribution of immiscible product within the unsaturated zone and the relative magnitude of aqueous-phase advection (ground water recharge) to gaseous-phase diffusion. At a gasoline spill site in Laurel Bay, South Carolina, rates of MTBE volatilization from ground water downgradient from the source are estimated by analyzing the distribution of MTBE in the unsaturated zone above a solute plume. Volatilization rates of MTBE from ground water determined by transport modeling ranged from 0.0020 to 0.0042 g m(-2)/year, depending on the assumed rate of ground water recharge. Although diffusive conditions at the Laurel Bay site are favorable for volatilization, mass loss of MTBE is insignificant over the length (230 m) of the solute plume. Based on this analysis, significant volatilization of MTBE from ground water downgradient from source areas at other sites is not likely. In contrast, model results indicate that volatilization coupled with diffusion to the atmosphere could be a significant mass loss pathway for MTBE in source areas where residual product resides above the capillary zone. Although not documented, mass loss of MTBE at the Laurel Bay site due to volatilization and diffusion to the atmosphere are predicted to be two to three times greater than mass loading of MTBE to ground water due to dissolution and recharge. This result would imply that volatilization in the source zone may be the critical natural attenuation pathway for MTBE at gasoline spill sites, especially when considering capillary zone limitations on volatilization of MTBE from ground water and the relative recalcitrance of MTBE to biodegradation.     

Comparison of alternative upscaled model formulations for simulating DNAPL source dissolution and biodecay

This paper compares the performance of analytical and numerical approaches for modeling DNAPL dissolution with biodecay. A solution derived from a 1-D advective transport formulation (“Parker” model) is shown to agree very closely with high resolution numerical solutions. A simple lumped source mass balance solution in which with decay is assumed proportional to DNAPL mass (“Falta1” model) over- or underpredicts aqueous phase biodecay depending on the magnitude of the exponential factor governing the relationship between dissolution rate and DNAPL mass. A modification of the Falta model that assumes decay proportional to the source exit concentration is capable of accurately simulating source behavior with strong aqueous phase biodecay if model parameters are appropriately selected or calibrated (“Falta2” model). However, parameters in the lumped models exhibit complex interdependencies that cannot be quantified without consideration of transport processes within the source zone. Combining the Falta2 solution with relationships derived from the Parker model was found to resolve these limitations and track the numerical model results. A method is presented to generalize the analytical solutions to enable simulation of partial mass removal with changes in source parameters over time due to various remedial actions. The algorithm is verified by comparison with numerical simulation results. An example application is presented that demonstrates the interactions of partial mass removal, enhanced biodecay, enhanced mass transfer and source zone flow reduction applied at various time periods on contaminant flux reduction. Increasing errors that arise in numerical solutions with coarse discretization and high decay rates are shown to be controlled by using an adjusted decay coefficient derived from the Parker analytical solution.     

A Practical Modeling Framework for Non‐Fickian Transport and Multi‐Species Sequential First‐Order Reaction

Many studies indicate that small‐scale heterogeneity and/or mobile–immobile mass exchange produce transient non‐Fickian plume behavior that is not well captured by the use of the standard, deterministic advection‐dispersion equation (ADE). An extended ADE modeling framework is presented here that is based on continuous time random walk theory. It can be used to characterize non‐Fickian transport coupled with simultaneous sequential first‐order reactions (e.g., biodegradation or radioactive decay) for multiple degrading contaminants such as chlorinated solvents, royal demolition explosive, pesticides, and radionuclides. To demonstrate this modeling framework, new transient analytical solutions are derived and are inverted in Laplace space. Closed‐form, steady‐state, multi‐species analytical solutions are also derived for non‐Fickian transport in highly heterogeneous aquifers with linear sorption–desorption and matrix diffusion for use in spreadsheets. The solutions are general enough to allow different degradation rates for the mobile and immobile zones. The transient solutions for multi‐species transport are applied to examine the effects of source remediation on the natural attenuation of downgradient plumes of both parent and degradation products in highly heterogeneous aquifers. Results for representative settings show that the use of the standard, deterministic ADE can over‐estimate cleanup rates and under‐predict the cleanup timeframe in comparison to the extended ADE analytical model. The modeling framework and calculations introduced here are also applied for a 30 year groundwater cleanup program at a site in Palm Bay, Florida. The simulated plume concentrations using the extended ADE exhibited agreement with observed long concentration tails of trichloroethene, cis 1,2 DCE, and VC that remained above cleanup goals.     

Matrix diffusion-derived plume attenuation in fractured bedrock

Matrix diffusion can attenuate the rate of plume migration in fractured bedrock relative to the rate of ground water flow for both conservative and nonconservative solutes of interest. In a system of parallel, equally spaced constant aperture fractures subject to steady-state ground water flow and an infinite source width, the degree of plume attenuation increases with time and travel distance, eventually reaching an asymptotic level. The asymptotic degree of plume attenuation in the absence of degradation can be predicted by a plume attenuation factor, beta, which is readily estimated as R' (phi(m)/phi(f)), where R' is the retardation factor in the matrix, phi(m) is the matrix porosity, and phi(f) is the fracture porosity. This dual-porosity relationship can also be thought of as the ratio of primary to secondary porosity. Beta represents the rate of ground water flow in fractures relative to the rate of plume advance. For the conditions examined in this study, beta increases with greater matrix porosity, greater matrix fraction organic carbon, larger fracture spacing, and smaller fracture aperture. These concepts are illustrated using a case study where dense nonaqueous phase liquid in fractured sandstone produced a dissolved-phase trichloroethylene (TCE) plume approximately 300 m in length. Transport parameters such as matrix porosity, fracture porosity, hydraulic gradient, and the matrix retardation factor were characterized at the site through field investigations. In the fractured sandstone bedrock examined in this study, the asymptotic plume attenuation factors (beta values) for conservative and nonconservative solutes (i.e., chloride and TCE) were predicted to be approximately 800 and 12,210, respectively. Quantitative analyses demonstrate that a porous media (single-porosity) solute transport model is not appropriate for simulating contaminant transport in fractured sandstone where matrix diffusion occurs. Rather, simulations need to be conducted with either a discrete fracture model that explicitly incorporates matrix diffusion, or a dual-continuum model that accounts for mass transfer between mobile and immobile zones. Simulations also demonstrate that back diffusion from the matrix to fractures will likely be the time-limiting factor in reaching ground water cleanup goals in some fractured bedrock environments.     

Multidimensional analytical models for isotope ratios in groundwater pollutant plumes of organic contaminants undergoing different biodegradation kinetics

This work presents analytical models which are able to predict contours of concentrations and isotope ratios of organic pollutants in homogeneous aquifers. Four analytical solutions of the advective–dispersive transport equation for reactive transport from the literature differing in assumptions regarding biodegradation kinetics were used. Stable isotope ratios are computed after modelling the individual reactive transport of isotopic species in the aquifer, which respond differently to fractionation by biodegradation or sorption. The main finding of this study is that the isotope ratios in the plumes are very sensitive to the assumptions underlying the biodegradation kinetics in the models. When biodegradation occurs throughout the core of the plume as first-order reaction, the transversal gradients in isotope ratios are smooth. When biodegradation occurs in a bi-molecular reaction with an electron acceptor (modelled by double-Monod kinetics), steep transversal isotope gradients are predicted. When the reaction rates approach instantaneous reaction along the plume fringes, isotope shifts in the core of the plume disappear. A model incorporating plume and fringe degradation produces the most plausible predictions of isotope ratios in this study. It is shown furthermore that isotope fractionation by sorption causes an even different pattern of isotope ratios, with positive shifts restricted to near the forerunning front of an expanding plume. The models developed in this work can serve for the validation of numerical models and may be incorporated in natural attenuation support systems such as e.g. BIOSCREEN.     

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.     

Analytical solutions for solute transport in saturated porous media with semi-infinite or finite thickness

Three-dimensional analytical solutions for solute transport in saturated, homogeneous porous media are developed. The models account for three-dimensional dispersion in a uniform flow field, first-order decay of aqueous phase and sorbed solutes with different decay rates, and nonequilibrium solute sorption onto the solid matrix of the porous formation. The governing solute transport equations are solved analytically by employing Laplace, Fourier and finite Fourier cosine transform techniques. Porous media with either semi-infinite or finite thickness are considered. Furthermore, continuous as well as periodic source loadings from either a point or an elliptic source geometry are examined. The effect of aquifer boundary conditions as well as the source geometry on solute transport in subsurface porous formations is investigated.     

Characterization and Source History Modeling Using Low‐k Zone Profiles at Two Source Areas

Like tree rings, high‐resolution soil sampling of low‐permeability (low‐k) zones can be used to evaluate the style of source history at contaminated sites (i.e., historical pattern of concentration and composition vs. time since releases occurred at the interface with the low‐k zone). This is valuable for the development of conceptual site model (CSM) and can serve as an important line of evidence supporting monitored natural attenuation (MNA) as a long‐term remedy. Source histories were successfully reconstructed at two sites at Naval Air Station Jacksonville using a simple one‐dimensional (1D) model. The plume arrival time and historical composition were reconstructed from the time initial releases that were suspected to occur decades earlier. At the first site (Building 106), the source reconstructions showed relatively constant source concentrations, but significant attenuation over time in the downgradient plume in the transmissive zone, suggesting MNA may not be an appropriate remedy if source control is a requirement, but attenuation processes are clearly helping to maintain plume stability and reduce risk. At the second site (Building 780), source concentrations in the transmissive zone showed an approximately a one order of magnitude over time, but apparently less attenuation in the downgradient plume. The source reconstruction method appeared to reflect site remediation efforts (excavation, soil vapor extraction) implemented in the 1990s. Finally, a detailed analysis using molecular biological tools, carbon isotopes, and by‐products suggests that most degradation activity is associated with high‐k zones but not with low‐k zones at these source areas. Overall, the source reconstruction methodology provided insight into historical concentration trends not obtainable otherwise given the limited long‐term monitoring data.     

Irreversible phosphorus sorption in septic system plumes?

The mobility of phosphorus (P) in septic system plumes remains a topic of debate because of the considerable reactivity of this constituent. In this study, a septic system plume in Ontario was monitored over a 16-year period with detail that clearly shows the advancing frontal portion of the P plume. This monitoring record provides insight into the extent of secondary P attenuation in the ground water zone beyond that available from previous studies. A P plume 16 m in length developed over the monitoring period with PO(4)-P concentrations (3 to 6 mg/L) that approached the concentrations present under the tile bed. Simulations using an analytical model showed that when first-order solute decay was considered to account for the possibility of secondary P attenuation in the ground water zone, field values could only be matched when decay was absent or occurred at an exceedingly slow rate (half-life greater than 30 years). Thus, hypothesized secondary P attenuation mechanisms such as slow recystallization of sorbed P into insoluble metal phosphate minerals, diffusion into microsites, or kinetically slow direct precipitation of P minerals such as hydroxyapatite were inactive in the ground water zone at this site or occurred at rates that were too slow to be observed in the context of the current 16-year study. Desorption tests on sediment samples from below the tile bed indicated a PO(4) distribution coefficient (K(d)) of 4.8, which implies a P retardation factor of 25, similar to the field apparent value of 37 determined from model calibrations. This example of inactive secondary P attenuation in the ground water zone shows that phosphorus in some ground water plumes can remain mobile and conservative for decades. This has important implications for septic systems located in lakeshore environments when long-term usage scenarios are considered.     

Eulerian spatial moments for solute transport in three-dimensional heterogeneous,dual-permeability media

A Eulerian analytical method is developed for nonreactive solute transport in heterogeneous, dual-permeability media where the hydraulic conductivities in fracture and matrix domains are both assumed to be stochastic processes. The analytical solution for the mean concentration is given explicitly in Fourier and Laplace transforms. Instead of using the fast fourier transform method to numerically invert the solution to real space (Hu et al., 2002), we apply the general relationship between spatial moments and concentration (Naff, 1990; Hu et al., 1997) to obtain the analytical solutions for the spatial moments up to the second for a pulse input of the solute. Owing to its accuracy and efficiency, the analytical method can be used to check the semi-analytical and Monte Carlo numerical methods before they are applied to more complicated studies. The analytical method can be also used during screening studies to identify the most significant transport parameters for further analysis. In this study, the analytical results have been compared with those obtained from the semi-analytical method (Hu et al., 2002) and the comparison shows that the semi-analytical method is robust. It is clearly shown from the analytical solution that the three factors, local dispersion, conductivity variation in each domain and velocity convection flow difference in the two domains, play different roles on the solute plume spreading in longitudinal and transverse directions. The calculation results also indicate that when the log-conductivity variance in matrix is 10 times less than its counterpart in fractures, it will hardly influence the solute transport, whether the conductivity field is matrix is treated as a homogeneous or random field.     

Methodology for comparing source and plume remediation alternatives

It is often difficult at contaminated sites to decide whether remediation effort should be focused on the contaminant source, the dissolved plume, or on both zones. The decision process at these sites is hampered by a lack of quantitative tools for comparing remediation alternatives. A new screening-level mass balance approach is developed for simulating the transient effects of simultaneous ground water source and plume remediation. The contaminant source model is based on a power function relationship between source mass and source discharge, and it can consider partial source remediation at any time after the initial release. The source model serves as a time-dependent mass flux boundary condition to a new analytical plume model, where flow is assumed to be one dimensional, with three-dimensional dispersion. The plume model simulates first-order sequential decay and production of several species, and the decay rates and parent/daughter yield coefficients are variable functions of time and distance. This new method allows for flexible simulation of natural attenuation or remediation efforts that enhance plume degradation. The plume remediation effort may be temporary or delayed in time, limited in space, and it may have different chemical effects on different contaminant species in the decay chain.     

Analytical solutions to the transient,unsaturated transport of water and contaminants through horizontal porous media

We present a range of analytical solutions to the combined transient water and solute transport for horizontal flow. We adopt the concept of a scale and time dependent dispersivity used for contaminant transport in aquifers and apply it to transient, unsaturated horizontal flow to develop similarity solutions for both constant solute concentration and solute flux boundary conditions. Through the use of a specific form of the water profile as used by Brutsaert [Water Resour Res 1968:4;785], the solute profiles can be reduced to a simple quadrature. We also derive a solution for the instantaneous injection of water and solute into a horizontal media for an arbitrary dispersivity. It is found that the solute concentration remains constant in both space and time as the water redistributes, suggesting that the solute does not disperse relative to the water.     

Effects of pore-scale dispersion,degree of heterogeneity,sampling size,and source volume on the concentration moments of conservative solutes in heterogeneous formations

Pore-scale dispersion (PSD), aquifer heterogeneity, sampling volume, and source size influence solute concentrations of conservative tracers transported in heterogeneous porous formations. In this work, we developed a new set of analytical solutions for the concentration ensemble mean, variance, and coefficient of variation (CV), which consider the effects of all these factors. We developed these models as generalizations of the first-order solutions in the log-conductivity variance of point concentration proposed by [Fiori A, Dagan G. Concentration fluctuations in aquifer transport: a rigorous first-order solution and applications. J Contam Hydrol 2000;45(1–2):139–163]. Our first-order solutions compare well with numerical simulations for small and moderate formation heterogeneity and from small to large sampling and source volumes. However, their performance deteriorates for highly heterogeneous formations. Successively, we used our models to study the interplay among sampler size, source volume, and PSD. Our analysis shows a complex and important interaction among these factors. Additionally, we show that the relative importance of these factors is also a function of plume age, of aquifer heterogeneity, and of the measurement location with respect to the mean plume center of gravity. We found that the concentration moments are chiefly controlled by the sampling volume with pore-scale dispersion playing a minor role at short times and for small source volumes. However, the effect of the source volume cannot be neglected when it is larger than the sampling volume. A different behavior occurs for long periods, which may be relevant for old contaminations, or for small injection volumes. In these cases, PSD causes a significant dilution, which is reflected in the concentration statistics. Additionally, at the center of the mean plume, where high concentrations are most likely to occur, we found that sampling volume and PSD are attenuating mechanisms for both concentration ensemble mean and coefficient of variation, except at very large source and sampler sizes, where the coefficient of variation increases with sampler size and PSD. Formation heterogeneity causes a faster reduction of the ensemble mean concentrations and a larger uncertainty at the center of the mean plume. Therefore, our results highlight the importance of considering the combined effect of formation heterogeneity, exposure volume, PSD, source size, and measurement location in performing risk assessment.