A plume-scale mass balance is developed to assess the natural attenuation (NA) of dissolved organic contaminants in fractured, dual-porosity aquifers. This methodology can be used to evaluate contaminant distribution within the aquifer, plume source term, contaminant biodegradation and plume status. The approach is illustrated for a site on the UK Upper Chalk aquifer impacted by petroleum fuel containing methyl tert-butyl ether (MTBE) and tert-amyl methyl ether (TAME). Variability in site investigation data and uncertainty in the mass balance were assessed using probabilistic analysis. The analysis shows that benzene, toluene, m/p-xylene and o-xylene (BTEX) compounds are biodegraded primarily by denitrification and sulfate reduction in the aquifer, with an equivalent plume-scale first-order biodegradation rate of 0.49 a−1. Other biodegradation processes are less important. Sorption contributes to hydrocarbon attenuation in the aquifer but is less important for MTBE and TAME. Uncertainty in the plume source term and site hydrogeological parameters had the greatest effect on the mass balance. The probabilistic analysis enabled the most likely long-term composition of the plume source term to be deduced and provided a site-specific estimate of contaminant mass flux for the prediction of plume development. The mass balance methodology provides a novel approach to improve NA assessments for petroleum hydrocarbons and other organic contaminants in these aquifer settings.

Supplementary material: Plume organic chemistry, aquifer hydraulic conductivity and probability distribution functions used for mass balance inputs are available at https://doi.org/10.6084/m9.figshare.c.7016429

Thematic collection: This article is part of the Monitoring the aquifers collection available at: https://www.lyellcollection.org/topic/collections/monitoring-the-aquifers

Monitored natural attenuation (MNA) is an established risk-based management approach for groundwater contaminated with petroleum hydrocarbons (Wiedemeier et al. 1999; Rivett and Thornton 2008). Technical guidance is available to assist with the implementation of MNA (e.g. Wiedemeier et al. 1995, 1999; Environment Agency 2000; ASTM 2016) but this is difficult to apply at sites with complex source histories and typically provides an indirect estimate of contaminant mass loss using rate constants obtained from limited groundwater monitoring (Wilson et al. 2004; Hosseini 2009; Thornton et al. 2013). Rate constants derived in this way may be in error or neglect the variability in contaminant concentration (Ravi et al. 1997; Environment Agency 2002; Hosseini 2009). However, mass balances can provide significant improvements in process understanding and assessment of contaminant fate in plumes (e.g. King et al. 1999; Thornton et al. 2001). Furthermore, the performance assessment of natural attenuation (NA) commonly uses deterministic analysis to predict plume behaviour and risk to receptors (Ravi et al. 1997; Wiedemeier et al. 1999; Hosseini 2009; Zhang 2010). This seldom considers the variability in environmental data, and in many cases site investigation for MNA assessments does not yield the extensive dataset required to evaluate uncertainty in predictions (McNab and Dooher 1998; Environment Agency 2001a, b; H. Liang et al. 2011; S.H. Liang et al. 2011; Beck and Mann 2012). Sources of uncertainty include knowledge of representative parameter values, the underlying conceptual model and observation data (US EPA 1999; Environment Agency 2001a, b; Anderson et al. 2007; Winter and Tartakovsky 2008; Zhang 2010; H. Liang et al. 2011; S.H. Liang et al. 2011; Beck and Mann 2012; Middlemis and Peeters 2018; La Cecilia et al. 2020). Deterministic analysis may also result in overly conservative predictions with higher costs to achieve remediation objectives (Anderson et al. 2007; H. Liang et al. 2011; S.H. Liang et al. 2011; Rivera-Velasquez et al. 2013). Instead, performance assessment and evaluations of contaminant fate and transport should include probabilistic analysis to account for incomplete conceptual understanding of the system interactions and uncertainty in parameter values describing processes (Woodbury et al. 1995; Environment Agency 2001a; Bredehoeft 2003, 2005; Bolster et al. 2009; Beck and Mann 2012). This approach is often incorporated within environmental risk assessments (e.g. Anderson et al. 2007; Ho et al. 2007; Winter and Tartakovsky 2008; Liang 2009; Zhang 2010; D. Zhang et al. 2010; P. Zhang et al. 2010; Rodak and Silliman 2012; Yang et al. 2018; Gonzalez et al. 2020; La Cecilia et al. 2020; Troldborg et al. 2022) but has seldom been used to interpret the NA of plumes at field scale (Thornton et al. 2001; Hosseini 2009; Liang 2009; O'Malley and Vesselinov 2014; Xie et al. 2016).

This paper presents a mass balance methodology integrated with probabilistic analysis for the field-scale performance assessment of petroleum hydrocarbon NA in fractured, dual-porosity aquifers. The methodology is based on a conceptual model that describes the distribution and attenuation of potentially biodegradable contaminants in a dual-porosity aquifer system, and enables the calculation of mass balances for contaminant plumes in groundwater using data obtained from a typical site investigation and groundwater monitoring study for MNA. It is illustrated for a petroleum fuel-release site on the UK Upper Chalk aquifer. The Chalk is the most important aquifer in the UK, accounting for 60% of groundwater used in England and Wales (Allen et al. 1997; MacDonald and Allen 2001). It is a soft, white limestone (c. 99% CaCO3) that forms a deep consolidated, fractured dual-porosity aquifer, with minor marl and flint layers but very low mineral oxide (MnO2, Fe(OH)3) content (Bottrell et al. 2010). The aquifer matrix has a high primary porosity, typically 30–45% in the Upper Chalk (Jones and Robins 1999), but low effective permeability (Allen et al. 1997). The matrix porewater has a long residence time compared with the fracture water, and within the saturated zone can be considered essentially as immobile storage relative to the timescale of plume development investigated here (Younger and Elliot 1995; Bloomfield 1996; Allen et al. 1997). The matrix is dissected by fractures that contribute minor secondary porosity (typically only 1%) but which have very high permeability and may contribute 99% of the flow (Price 1996; Allen et al. 1997; MacDonald and Allen 2001). Groundwater flow and solute transport in the system occurs primarily in the fractures, which have high transmissivity, whereas the low transmissivity matrix provides storage for the bulk of the groundwater (Bloomfield 1996; Allen et al. 1997; MacDonald and Allen 2001). Solute exchange between the matrix and fractures occurs by molecular diffusion (Younger and Elliot 1995).

The purpose of this paper is to evaluate the mass balance methodology in this contaminated aquifer setting, considering dual-porosity effects on contaminant distribution and attenuation, and the contribution of parameter uncertainty in predictions of NA performance. The specific objectives were to:

  • undertake a plume-scale mass balance to evaluate the contaminant distribution and NA of dissolved organic contaminants in groundwater at field scale for the site under study; and

  • evaluate the effect of variability in parameter inputs on the mass balance for the site and identify key parameters responsible for uncertainty in the performance prediction for NA, using probabilistic methods.

The site is a former retail filling station located on the Upper Chalk aquifer in southern England (Fig. 1). Petroleum fuel was accidently released from an underground storage tank (UST) in February 1999. This has contaminated the saturated zone of the aquifer beneath the site with light non-aqueous phase liquid (LNAPL), which has migrated through vertical fractures up to 15 m below the water table, in response to the pressure head developed by the connected vertical height of LNAPL in fractures (Wealthall et al. 2002). The theoretical basis for this behaviour is well understood (Mercer and Cohen 1990; Newell et al. 1995; Hardisty et al. 1998) and relatively small volumes of LNAPL within vertical or subvertical fractures can produce significant LNAPL pressure heads that result in penetration into the saturated zone (Hardisty et al. 1998, 2003; Thornton and Wealthall 2008; CL:AIRE 2014). Dissolution of this residual LNAPL has contaminated the groundwater with petroleum hydrocarbons, including benzene, toluene, m/p-xylene, o-xylene (BTEX), the fuel additives methyl tert-butyl ether (MTBE), tert-amyl methyl ether (TAME) and tert-butyl alcohol (TBA), and other aromatic compounds (Wealthall et al. 2002; Spence et al. 2005a).

The aquifer is unconfined in the vicinity of the site and overlain by c. 7 m of drift deposits. The mean water table is c. 20–22 m below ground level. The groundwater flow direction is southeasterly, approximately parallel to the main road adjacent to the filling station (Spence et al. 2005a) (Fig. 1). A network of single-screen (c. 9 m screen length) monitoring wells (MW7, MW9, MW10, MW11, MW13, MW22 and MW24) and multilevel groundwater sampling (MLS) wells (MW14, MW15 MW16, MW17, MW18, MW19, MW20 and MW23), which allow level-discrete sampling over the vertical interval of the plume (Einarson and Cherry 2002; Thornton and Wealthall 2008), were installed at the site (Fig. 1). MW14 was installed upgradient of the site to characterize the background groundwater quality. Similarly, MW23 was installed in the saturated zone underneath the former leaking UST. All monitoring wells are referenced relative to this datum (as metres below forecourt level (m bfl)) using ordnance survey data.

Groundwater quality surveys of dissolved contaminant distribution were completed at intervals of 1–3 months after the fuel release. These show that the locations of MW10 and MW11 are uncontaminated and not on the plume flow path. Three comprehensive surveys of groundwater chemistry were completed for this study on all monitoring wells at annual intervals between 2001 and 2003.

Site investigation

The site investigation undertaken at the site is described in Wealthall et al. (2002). The methods used for groundwater sampling and analysis are explained in Spence et al. (2005a). Cored boreholes were drilled in the saturated zone at the locations of MW15, MW16, MW17, MW20 and MW23 using rotary water flush methods, whereas rock core from the unsaturated zone at MW23 was recovered using percussion drilling techniques (Spence et al. 2005a, b).

Vertical profiles of solute distribution were obtained from the MLS completions after installation of a continuous multichannel tubing (CMT) system, containing seven monitoring ports fitted with 10 cm screens (Einarson and Cherry 2002; Thornton and Wealthall 2008). The MLS completions were installed to a depth of 55 m, with monitoring intervals of between 1 and 2 m in the plume wells (MW15, MW16, MW17, MW19 and MW23) and up to 7.6 m in the upgradient well (MW14).

Groundwater samples for chemical analysis were collected according to Spence et al. (2005a). Dissolved hydrocarbon, MTBE, TAME and TBA concentrations were determined using the method of Dewsbury et al. (2003) and other chemical parameters were determined using standard methods (Spence et al. 2005a). Rock cores were sampled anaerobically at intervals of c. 1 m for porewater chemistry, as described by Spence et al. (2005b). This analysis also provided an estimate of the saturated porosity of samples over the same interval as the porewater chemistry. The chemical composition of the porewater samples was determined using the methods previously indicated. Rock core from MW20 (uncontaminated location) was sampled at 1 m intervals and analysed for its solid-phase organic carbon content using standard methods (van Reeuwijk 1987).

Conceptual model and assumptions

A conceptual model to estimate the contaminant mass balance in the saturated zone was developed to assess the fate of dissolved-phase contaminants released to groundwater from the dissolution of residual LNAPL trapped in fractures below the water table. The model is presented graphically in Figure 2 and the calculation of the different terms in the model is described in detail below. In essence, the model considers saturated zone processes only, so the net release is the total spill mass (Mspill) less any product recovered from the unsaturated zone (e.g. by soil vapour extraction (SVE)).

In the saturated zone contaminant is initially present as a separate non-aqueous phase that dissolves into flowing water present in fractures (MSZ fract dissolved). This input over the history of the plume is estimated from groundwater monitoring well data, within the term MSZ dissolved released. The current mass of dissolved contaminants in the saturated zone (MSZ dissolved) is also estimated using the monitoring well data. Dissolved contaminants in the fracture water may biodegrade (included in the term MSZ degraded) or diffuse into porewater within the Chalk matrix (Mresidual ED). Similarly, electron acceptors (EA) within the aquifer matrix bounded by the plume that have been consumed by biodegradation in the fractures is also included (Mresidual EA). Biodegradation within fracture waters is assumed to be driven by electron acceptors advected into the source zone from upgradient (Madvected EA) and electron acceptors stored in the matrix porewater within the plume volume that diffuse into fracture waters as biodegradation lowers EA concentrations in the fracture waters (Mdiffusion EA). Biodegradation is considered to only occur within the fractures, as pore throat sizes in the Chalk are typically too small for microbes to access (e.g. Fredrickson et al. 1997; Bottrell et al. 2000; Rebata-Landa and Santamarina 2006; Thornton et al. 2006; Bartlett et al. 2010). Contaminant that has diffused into the Chalk matrix porewater is therefore treated as an immobile reservoir, part of which may partition by sorption to the small amount of organic matter present in the Chalk matrix (MSZ sorbed). The overall mass of contaminant biodegraded in the saturated zone is represented by MSZ biodegraded.

The relevant input terms governing the distribution of dissolved phase fractions are:
where MSZ dissolved released is the mass dissolved from the trapped LNAPL to form the plume, and MSZ dissolved, MSZ fract dissolved, MSZ sorbed and MSZ biodegraded represent the dissolved mass present in the aquifer matrix, fractures, sorbed to the aquifer matrix and biodegraded in the groundwater, respectively.
MSZ biodegraded includes contributions from the following sources:
where Madvected EA represents biodegradation due to the consumption of aqueous electron acceptors (O2, NO3, SO4) supplied by the advection of uncontaminated groundwater into the plume, Mresidual EA reflects the biodegradation due to Mn(IV)/Fe(III) oxide reduction, SO4 reduction and methanogenesis within the plume, and Mdiffusion EA reflects the biodegradation due to the consumption of aqueous electron acceptors originally present in the aquifer matrix bounded by the plume. Biodegradation reactions are assumed to be instantaneous for petroleum hydrocarbons (Wiedemeier et al. 1999).

The plume is represented as a box model, in which reactions describing contaminant attenuation are coupled with transport processes to determine the spatial and temporal distribution of mass in the system. Following Anderson et al. (2007), the plume width, depth and length are assumed to be defined, respectively, by the distance between contaminated monitoring wells transverse to the flow path, the depth of dissolved contaminants in MLS boreholes close to the source and the maximum downgradient migration of dissolved contaminants, determined from the MLS installations. Monitoring data show that the groundwater flow direction is consistent over the reference period of the plume (December 1999–July 2003).

Inputs in the mass balance are estimated from measurements of the groundwater hydrochemistry, aquifer hydrogeology, geochemical properties and other site data. The groundwater chemistry of MW7 and MW13 is assumed to represent the potential plume source composition. This is justified as the MLS profiles show that contaminants are not present in the saturated zone below the former leaking tank (at the location of MW23), and MW7, 16 m down hydraulic gradient of MW23, contains the highest concentrations of dissolved hydrocarbons in groundwater at the site.

Plume source term input

This input was evaluated using monthly time-series data for contaminant concentrations in MW7 (Fig. 3) and MW13 (Fig. 4). These monitoring wells have the highest dissolved contaminant concentrations and are immediately down hydraulic gradient of the UST location (MW23). The total mass of dissolved contaminants released into the plume over this period, MSZ dissolved released, was obtained from the following relationship:
where Csource (M/L3) is the concentration of each contaminant in MW7 or MW13, W is the plume width (L), I is the hydraulic gradient (L/L), K is the aquifer hydraulic conductivity (L/T), thicknessK interval is the packer test interval (aquifer thickness) used in the hydraulic tests to determine the aquifer hydraulic conductivity at successive depths (L), and T is the time period of the analysis (December 1999 – July 2003).

To estimate the contaminant mass flux from MW7 or MW13, contaminant concentrations were assumed to be consistent over the vertical interval of aquifer monitored by each well, to the full depth of the plume. This is reasonable as these monitoring wells have 9 m screens and vertical mixing of groundwater is expected within the casing, leading to homogenization of vertical contaminant concentration gradients across the screen (Puls and Paul 1997; Schreiber and Bahr 1999). The temporal variation in dissolved contaminant concentrations in MW7 and MW13 is likely to result from groundwater intermittently accessing fresh LNAPL trapped in fractures in the unsaturated zone and an increased flux of dissolved contaminants through a zone of more intensive fracturing in the shallow saturated zone (deduced from rock core fracture logs and tracer tests), due to seasonal variations in recharge and a fluctuating water table (Hardisty et al. 1998; Bottrell et al. 2010). These are characteristic features of the Chalk aquifer properties and hydrogeology (Allen et al. 1997). There is a significant variation in the aquifer vertical hydraulic conductivity, which will determine the vertical variation in the contaminant mass flux over the depth of the plume. The hydraulic conductivity profile for MW15 (see Supplementary Fig. 1), the MLS borehole closest to MW7 and MW13 (Fig. 1), was used in this calculation. The plume width, W, is assumed to be the distance between the contaminated monitoring wells, MW13 and MW19, transverse to the plume flow path.

Dissolved contaminant mass

Inputs of MSZ dissolved and MSZ fract dissolved were estimated using xz contour plots of contaminant concentrations obtained from the June 2002 groundwater monitoring survey (Fig. 4). The contaminant concentrations in 2002 were greater than in the 2001 and 2003 surveys, representing the maximum input for the mass balance. Due to different porosity terms, MSZ fract dissolved was estimated using:
where εf is the fracture porosity (dimensionless), Cresidual ED is the mean residual electron donor (organic contaminant) concentration for a given contour interval (M/L3) and V is the volume of aquifer contained within the contour interval (L3). A fracture porosity value of 0.01 was used in the calculations for MSZ fract dissolved, as typical for the Upper Chalk aquifer (Price 1996). Values of MSZ dissolved were estimated using:
where εm is the matrix porosity (dimensionless) and all other terms are as for MSZ fract dissolved. Values of εm were obtained from analyses (n = 17) of rock core collected below the location of the former UST at MW23 (Spence et al. 2005b). Only εm values for the monitored depth of the plume (20–40 m) were used.

Sorbed contaminant mass

Sorption of hydrophobic organic compounds such as petroleum hydrocarbons to the aquifer matrix can be approximated by the following expression (Schwarzenbach et al. 1993):
where S is sorbed solute concentration (Msolute/Msorbent), C is dissolved solute concentration (Msolute/L3solution) and Kd is the equilibrium distribution coefficient (L3solution/Msorbent). Kd can be used to estimate a retardation factor, Rf, describing the effect of sorption on solute transport in the porous medium (Fetter 1999):
where ρd is the bulk density of the porous medium (M/L3) and ε is the porosity (dimensionless). The contribution of fracture porosity was considered by replacing ε with (εmεf) in equation (7). The sorbed mass of a solute can then be estimated according to (King et al. 1999):
Equation (8) was used to estimate MSZ sorbed for each contaminant in the plume, using estimates of MSZ dissolved (equation 5) and Rf (equation 7). It was assumed that sorption occurs after diffusion into the matrix and is negligible for fracture surfaces. This is justified since contaminants diffuse from the fractures into the aquifer matrix, which has a significantly greater surface area. There is excellent agreement between vertical profiles of organic contaminants in matrix porewater from rock cores and fracture water chemistry in MLS wells installed in the same borehole (Spence et al. 2005b). This suggests that there is chemical equilibrium between contaminant concentrations in the fracture water and the matrix porewater, supporting the use of a linear sorption model in the calculations. Values of Kd for contaminants in the plume were estimated indirectly by two methods using empirical relationships from the literature, as described below.
For low solute concentrations, Kd, considered in terms of the organic carbon-normalized distribution coefficient, Koc, can be expressed as a linear function of the fraction of particulate organic carbon (foc) on the sorbent (aquifer matrix) using the relationship (Schwarzenbach et al. 1993):
where Koc has units of L3/M and foc is dimensionless. Published values of Koc were used for organic compounds (Fetter 1999). However, for 1,2,4-trimethylbenzene (1,2,4-TMB) and 1,3,5-trimethylbenzene (1,3,5-TMB) Koc was estimated using the general relationship (Schwarzenbach and Westall 1981):
where Kow is the octanol–water partition coefficient for the organic compound, and a and b are regression coefficients. Values of 1.01 and −0.72 were used, respectively, for the coefficients a and b in equation (10), based on the empirical expression derived for aromatic hydrocarbons in Schwarzenbach et al. (1993, p. 274, equations 11–22a), to estimate Koc for 1,2,4-TMB and 1,3,5-TMB. The Koc estimates were then used to calculate Kd values for these two compounds using equation (10). The foc values used to estimate organic contaminant Kd were obtained from chemical analysis of uncontaminated rock core collected from MW20 (Fig. 1) (Bottrell et al. 2010).

Biodegraded contaminant mass

Estimates of Madvected EA were obtained by comparing profiles of dissolved O2, NO3 and SO4 in the upstream (uncontaminated) MLS (MW14) with MW15, the MLS immediately downgradient of the plume source (Fig. 1). Concentrations of dissolved O2 and NO3 were found to be below detection limits in MW15 and dissolved SO4 was depleted, relative to MW14, on each annual survey (Spence et al. 2005a). This indicates complete consumption of O2 and NO3 supplied to the plume by advection, with partial consumption of SO4 in the plume. The contribution of each electron acceptor in the budget for Madvected EA was quantified using the following expression:
where Cbackground EA is the concentration (M/L3) of O2, NO3 and SO4 in uncontaminated groundwater (MW14), and the other terms are as in equation (3). Dissolved electron acceptor concentrations obtained from two monthly and three annual groundwater quality surveys of MW14 between 2000 and 2003 were used in the calculation of this input. It is assumed that concentrations of O2, NO3 and SO4 in groundwater entering the plume at the location of MW15 are identical to those in the profile at MW14 over the plume history. There is good agreement between profiles of these electron acceptors obtained from the annual surveys of MW14, suggesting that this is reasonable (Spence et al. 2005a). The hydraulic conductivity profile for MW15 (Supplementary Fig. 1) was used to calculate the mass flux of each electron acceptor entering the rear of the plume; this input is then a function of the aquifer K profile in MW15 and the electron acceptor profile in MW14. The electron acceptor mass flux profile obtained for the monitoring intervals in MW14 was then integrated using the trapezoid rule over the full depth (c. 15 m) of the plume profile at MW15.
Estimates of Mresidual EA are based on the residual concentration of Mn(II), Fe(II) and CH4 in the plume, which are higher than those in the uncontaminated groundwater at MW14. This input was determined by integration of the contoured xz concentration plots for the dissolved Mn(II), Fe(II) and CH4 distribution in 2002, using the volume of aquifer contained within a contoured interval and a mean solute concentration for that contour interval, according to:
where Cresidual EA is the mean concentration of the residual species for a given contour interval (M/L3). Although SO4 reduction occurs in the plume, SO4 was never depleted (Spence et al. 2005a) and this residual SO4 must be considered in the estimation of MSZ biodegraded. The depletion of SO4 in groundwater is strongly correlated with the distribution of BTEX (Spence et al. 2005a), and Mresidual EA for SO4 was therefore estimated with equation (12) by integration of the SO4 contour plot over the area occupied by the benzene plume in 2002. This year corresponds to the maximum spatial distribution of dissolved BTEX in the aquifer.

Mdiffusion EA is the contribution to MSZ biodegraded from the consumption of aqueous electron acceptors (O2, NO3 and SO4) formerly present in the uncontaminated aquifer matrix but now depleted by biodegradation in the plume. These electron acceptors diffuse from the matrix to support biodegradation in the fractures. This is supported by the excellent agreement between the matrix porewater and fracture water profiles of these species (Spence et al. 2005b). It indicates that there is rapid chemical equilibrium between the fracture water and the matrix porewater, and that the original reservoir of O2 and NO3 in the matrix porewater has diffused into the fractures and been consumed by biodegradation. The exception is SO4, which remains as a residual (incompletely used) fraction, with lower concentrations in the plume than in the uncontaminated aquifer.

In estimating Mdiffusion EA it was assumed that the aquifer matrix now occupied by the plume originally had dissolved O2, NO3 and SO4 concentrations (Cbackground EA) equivalent to those in uncontaminated groundwater at MW14. Also, the relationships noted above between the contaminant and aqueous electron acceptor distribution in the matrix porewater, fracture water and plume contour plots suggest that aqueous electron acceptors in the matrix will be depleted over a volume of aquifer corresponding to the BTEX plume. Values of Mdiffusion EA were therefore estimated for O2, NO3 and SO4 using the following expression:
where Vbenzene plume is the volume of the aquifer (L3) occupied by the benzene plume in 2002 (representing the maximum spatial distribution of dissolved BTEX contaminants).

Probabilistic analysis

Probabilistic analyses were performed using the Oracle Crystal Ball® software suite. Different scenarios were evaluated to assess the variation in mass balance predictions due to variabilities in parameter values. This was achieved using Monte Carlo methods and probability distribution functions (pdfs) assigned to each data range in Table 1. The significance of parameter uncertainty on the mass balance was also evaluated using a sensitivity analysis for each input. Uncertainty is assessed in the sensitivity analysis using pdfs assigned to parameters and is propagated through the various embedded forecasts in the simulation. In this way the probabilistic analysis assessed uncertainty in the conceptual model and input data used in the mass balance.

Monte Carlo simulations were run for 10 000 realizations and the significance of forecasts was evaluated at the 95% confidence interval (Environment Agency 2001b; Anderson et al. 2007). Specific pdfs describing parameter inputs for the Monte Carlo simulations were obtained by statistical analysis of the site data using the Crystal Ball® software (Supplementary Table 1) to account for the site-specific variability in parameter values (McNab and Dooher 1998; Soutter and Musy 1998; Wang et al. 2004; Anderson et al. 2007; H. Liang et al. 2011; S.H. Liang et al. 2011). However, in some cases an assumed pdf was assigned to a parameter input where data were limited according to guidance and previous studies (Environment Agency 2001a, b; Liang 2009; H. Liang et al. 2011; S.H. Liang et al. 2011).

Plume source term

Estimates of MSZ dissolved released were evaluated for three scenarios, considering time-variant contaminant concentrations for MW7 only, MW13 only and MW7 combined with MW13. This evaluated uncertainty in the conceptual model for the plume source term. Four values (0.0021, 0.0022, 0.0024 and 0.0031) of the hydraulic gradient, Ix, obtained from separate groundwater elevation surveys were used in this analysis. A uniform pdf was assigned to this range (Environment Agency 2001a, b), which was also used to calculate other input terms requiring this parameter.

Plume residuals

The mass of plume residuals (dissolved organic contaminants and electron acceptors/biodegradation products) was estimated by assigning a triangular pdf to each set of concentration (Cresidual ED or Cresidual EA) contours, using the lower, arithmetic mean and upper solute concentrations as the minimum, most likely and maximum values, respectively, in the distribution (Thornton et al. 2001). An extreme value probability distribution was obtained by fitting Crystal Ball® to the data range for matrix porosity, εm. This distribution was truncated at the minimum and maximum values of εm, as these defined the range used.

Sorbed contaminant mass

Three parameter values, foc, εm and ρb, were evaluated in the calculation of MSZ sorbed. Crystal Ball® was used to fit a Weibull pdf, truncated at the minimum and maximum values that define the range, to 29 values of foc, obtained by analysis of rock core from MW20. The values and pdfs for εm used in the calculation of plume residuals (see above) were included to estimate this input. Values of ρb were obtained using an expression that includes εm, and so the same probability distribution (Extreme value) was used for this parameter in the Monte Carlo simulations.

Biodegraded contaminant mass

In the calculation of Madvected EA all data from two monthly and three annual groundwater quality surveys of MW14 between 2000 and 2003 were used to determine the electron acceptor concentration, Cbackground EA, used in equation (12). This accounts for temporal variability in the background electron acceptor mass flux entering the plume via advection of uncontaminated groundwater. In this case Cbackground EA was obtained by fitting a pdf to the time-series data for each monitoring point on MW14 (seven in total) using Crystal Ball®. This provides different values of Cbackground EA as a function of depth and time over a significant part of the plume history. Estimates of Mdiffusion EA also included the pdfs fitted by Crystal Ball® to Cbackground EA.

Contaminant mass balance

The mass balance for the saturated zone of the aquifer is shown in Table 2. This includes estimates of the total contaminant mass and the mass of BTEX only, where relevant. Mean values of predictions and the output range between the 95% confidence level are provided for each forecast. All masses are in kilograms.

Values of MSZ dissolved released predicted using the MW7 data are higher than those using the MW13 data. These estimates vary by orders of magnitude at the 95% confidence interval. However, mean values of MSZ dissolved released predicted using a combination of MW7 and MW13 data vary much less. The residual contaminants in the saturated zone include the dissolved mass in the matrix (MSZ dissolved), fractures (MSZ fract dissolved) and sorbed to the matrix (MSZ sorbed). Residual dissolved contaminants in the matrix represent the largest fraction, with approximately 50% (110 kg) comprising BTEX alone. The dissolved contaminant mass in fractures is very small, due to the low fracture porosity. Most (94% based on mean values) of the contaminant fraction sorbed to the aquifer matrix is accounted for by BTEX. This group also accounts for c. 62% of the total residual contaminants (205 kg) in the saturated zone (Table 2).

Plume balance

The consumption of electron acceptors in the saturated zone (MSZ biodegraded) is equivalent to the biodegradation of 1595 kg of BTEX, with a relatively restricted range (Table 2). This is based on the reaction stoichiometry for individual biodegradation processes (Wiedemeier et al. 1999) and assumes that oxidant consumption results primarily from BTEX biodegradation, which is supported by stable isotope studies of contaminant biodegradation (Spence et al. 2005a).

The sum of the total residual contaminant and biodegraded hydrocarbon fractions can be compared with the predicted total dissolved contaminant mass released in the saturated zone, according to equation (1). The latter implies that the dissolved contaminant mass released into the groundwater from LNAPL dissolution should be balanced by the residual dissolved contaminant mass distributed between the different aquifer fractions and the contaminant mass biodegraded within the plume. The respective input terms in equation (14) are summarized in the ‘plume balance’ (Table 2). There is good agreement in the results when using the combined MW7 and MW13 data as the plume source term, particularly for the ‘all contaminant’ fraction. The residual and biodegraded contaminant mass is 1927 kg (range of 1647–2169 kg) and this brackets the mean value of MSZ dissolved released (1796 kg) predicted using the combined MW7 and MW13 data (Table 2).

Oxidant mass balance

The consumption of electron acceptors contributed by advection, matrix diffusion and plume residuals during biodegradation in the saturated zone (MSZ biodegraded) is presented in Table 3. An estimated 7681 kg of electron acceptors have been consumed, with aqueous oxidants dominating the budget.

The consumption of dissolved NO3 and SO4 accounts for 96% of the total electron acceptor budget in the plume, whereas Mn/Fe mineral oxides contribute very little (<1% in total). The contribution of methanogenesis in biodegradation is negligible. The mass balance suggests that the supply of dissolved O2 and NO3 by advection into the plume is more important for biodegradation than diffusion from the matrix, over the timescale of the analysis. The true contribution of advection and diffusion in the SO4 budget cannot be estimated because the residual SO4 in the plume may include contributions from both sources. Instead, the estimates of SO4 from advection and diffusion are maximum values over the timescale of the analysis, and the SO4 consumption is obtained by subtracting the residual plume SO4 estimate from the sum of these two components (Table 3).

Sensitivity analyses

The value of MSZ dissolved released was primarily sensitive to contaminant concentrations rather than the aquifer properties (hydraulic gradient). The aquifer matrix porosity was most important for estimates of MSZ dissolved, whereas estimates of MSZ sorbed are primarily affected by uncertainty in the foc of the aquifer matrix, which accounts for 93% of the variability. This is consistent with similar studies assessing the organic contaminant fate (e.g. Soutter and Musy 1998; Dubus et al. 2003) and is expected, as hydrocarbon sorption is strongly controlled by the organic matter content of the aquifer matrix.

Electron acceptor inputs in the mass balance are controlled primarily by the aquifer properties. Estimates of Madvected EA and Mdiffusion EA are influenced by, respectively, the hydraulic gradient and the matrix porosity, more than variability in the background groundwater electron acceptor concentration. Similarly, MSZ biodegraded also depended primarily on the variation in the aquifer hydraulic gradient through the influence on Madvected EA.

Regarding the plume balance, variability in the aquifer foc significantly (90%) controls the predicted total residual contaminants in the saturated zone (inputs 9 and 10 in Table 2). However, the sum of residual contaminants and equivalent BTEX biodegraded in the saturated zone (inputs 12 and 13 in Table 2) were most sensitive to variations in the aquifer hydraulic conductivity (79–81% of total variance). The variations in other aquifer and hydrochemical properties (e.g. matrix porosity, foc and background SO4 concentration) contributed much less (<15% in total) to the uncertainty in this prediction.

Sensitivity of mass balance to the uncertainty of input parameters

By identifying critical inputs that introduce most uncertainty in the mass balance, the sensitivity analysis indicates system parameters or aquifer properties that need greater characterization during a site investigation for performance assessment of NA using this methodology (Ho et al. 2007; D. Zhang et al. 2010; P. Zhang et al. 2010; Beck and Mann 2012). Additional data that account for a greater range in magnitude of a sensitive parameter or property may in turn significantly reduce the uncertainty in the mass balance input terms, providing a greater confidence in predictions (Hosseini 2009).

Most pdfs assigned to parameters were obtained by fitting the Crystal Ball® library to site-specific datasets. This means that the sensitivity analyses probably reflect the true effect of data variability on predictions rather than the assigned pdf (Anderson et al. 2007). A uniform distribution was assumed for the hydraulic gradient (Environment Agency 2001a, b) but a previous study by the authors shows that this property is relatively insensitive to the pdf assigned to it in a plume-scale mass balance (Thornton et al. 2001).

The mass balance is most sensitive to the aquifer hydraulic gradient, and less so to other aquifer properties such as foc, matrix porosity or the concentration of electron acceptors in background groundwater. This probably occurs because the hydraulic gradient typically exhibits significant variability in the Chalk aquifer (MacDonald and Allen 2001) and controls the groundwater flux through the plume. In turn, this controls the advective mass flux of aqueous electron acceptors into the plume and quantity of contaminants biodegraded.

Estimated plume source term

The plume source term was evaluated using contaminant concentrations in MW7 and MW13. These are the highest in the saturated zone but with considerable spatial–temporal variation. Both monitoring wells are transverse to the plume flow path (Fig. 1) and could therefore contribute dissolved contaminants to the plume. For this reason, the mass balance considered MW7 and MW13 in isolation and in combination.

Contaminant concentrations in MW13 are generally much lower than those in MW7 (cf. Figs 2 and 3) and the MSZ dissolved released predicted for the MW13 data is less than the sum of residual and biodegraded contaminants in the plume (inputs 12 and 13 in Table 2). This suggests that the plume is not generated by contaminant concentrations found exclusively in MW13. The large range in MSZ dissolved released predicted for the MW7 data results from the Monte Carlo sampling of the data used in this forecast. However, it is unrealistic that the highest value will define most of the mass released from this location because contaminant concentrations in MW7 peak in October 2000 and are much lower over most of the simulation period. This also suggests that the plume source is not described by the MW7 data exclusively.

Intuitively, the plume should result from dissolved contaminants contributed by both MW7 and MW13. This is because there are contaminated MLS boreholes downgradient of MW7 (e.g. MW15 and MW19) and MW13 (e.g. MW17 and MW16) along the plume flow path (Fig. 1). The range in forecast values of MSZ dissolved released produced from a combination of MW7 and MW13 data is less than that for the MW7 data (Table 2). Again, the range predicted by the probabilistic analysis is based on the Monte Carlo sampling of the combined dataset and these extreme values may not be realistic for the plume. This is because the possibility of all parameters taking their extreme values simultaneously is less likely than the possibility of any one of them taking its extreme value individually (Environment Agency 2001a, b). It is more appropriate to compare the mean values of predicted estimates, although the simulation will include statistical analysis of the distribution. A statistical distribution can be fitted to the MSZ dissolved released forecast for MW7 and MW13 using Crystal Ball® (Fig. 5). The forecast has a log-normal distribution, which allows the plume source term in the saturated zone to be predicted (with respect to the dissolved contaminant load) from the time-series contaminant concentration data for these two monitoring wells. This can also be undertaken to forecast future concentrations or mass flux estimates of individual or grouped contaminants. The uncertainty in this forecast can be further reduced by including more time-series data from these monitoring wells in the simulation.

Residual plume contaminant mass

The residual contaminant mass in the saturated zone, represented by dissolved compounds in the matrix (MSZ dissolved), fractures (MSZ fract dissolved) and sorbed to the matrix (MSZ sorbed), comprises 332 kg of total contaminants (Input 9 in Table 2). The respective estimates for the BTEX group are much lower (Input 10 in Table 2). This difference occurs because there is a much higher mass of dissolved MTBE and TAME in the matrix, rather than sorbed to the aquifer compared with BTEX. This is shown in Table 2, where estimates of MSZ dissolved for BTEX are approximately half those of the total contaminants, whereas estimates of MSZ sorbed for BTEX are only slightly less than those for the total contaminants. Due to higher solubility, MTBE and TAME will sorb less to the aquifer matrix compared with the BTEX compounds. This is illustrated by the respective distribution coefficients (Kd) and corresponding retardation factors (Rf) estimated for these compounds (see Supplementary Table 2), which show that significant attenuation of the hydrocarbon plumes by sorption can be expected. However, there will be negligible retardation of the MTBE and TAME plumes. This is confirmed by comparing groundwater velocities from tracer tests with the size of the MTBE plume, which indicates minimal attenuation by sorption (Bottrell et al. 2010).

Electron acceptor consumption

Advection is predicted to supply most of the dissolved O2, NO3 and SO4 consumed by biodegradation in the plume (Table 3). However, consumption of these electron acceptors originally present in the matrix now contained within the plume is also important. These respective masses (Madvected EA and Mdiffusion EA) must be considered maximum possible inputs in each case because the residual SO4 present in the plume cannot be separated into components contributed by advection or matrix diffusion. However, if all the residual SO4 in the plume was contributed by diffusion, the estimated (adjusted) value of Mdiffusion EA (2482 kg) would still not exceed the value of Madvected EA (3570 kg).

The contribution of electron acceptors in the matrix is finite and limited to the volume of aquifer occupied by the plume. It is reasonable to assume that aqueous O2, NO3 and SO4 originally present in the uncontaminated aquifer matrix has diffused into the fractures to support contaminant biodegradation in the plume (Spence et al. 2005a). The absence or reduced concentration of O2, NO3 and SO4 in the plume indicates that the matrix has been depleted in these oxidants. When this reservoir is exhausted, biodegradation will be primarily sustained by the supply of these electron acceptors from advection. However, as the plume expands downgradient into the aquifer, these oxidants will be supplied by diffusion from the uncontaminated matrix to support contaminant biodegradation in the fractures (Spence et al. 2005a). The rate of plume expansion will then be a function of the type and quantity of electron acceptors in the matrix, the composition of the contaminant plume, the matrix–fracture mass transfer rate (controlled by diffusion), and biodegradation kinetics. More contaminant mass can be biodegraded by NO3 and SO4 than O2, and the relatively high concentrations of NO3 and SO4 in the matrix has supported significant attenuation of the BTEX plume by biodegradation.

The contribution of Mn(IV) reduction, Fe(III) reduction and methanogenesis in contaminant biodegradation is minor compared with aerobic respiration, denitrification and sulfate reduction. The consumption of mineral oxidants in the plume is less than 1% of the total electron acceptor budget. The estimate for Fe is a minimum value because dissolved Fe can precipitate with sulfide to form iron sulfide phases (FeS) under the SO4-reducing conditions in the plume (Spence et al. 2005a). However, metal oxide coatings were observed on some core samples, suggesting that the contribution of mineral oxidants to biodegradation in the saturated zone is small. Methanogenesis is not an important pathway for contaminant degradation in the plume (Table 2), most probably due to the greater availability of NO3 and SO4 in the groundwater.

Plume balance

The sum of all dissolved contaminants released into the groundwater (Input 1 in Table 2) is in good agreement with the residual and biodegraded mass of all contaminants in the plume (Input 12 in Table 2). The difference in the mass balance is only 7%.

As the estimate of MSZ dissolved released (Input 1 in Table 2) and plume balance (inputs 12 and 13 in Table 2) for the combined MW7 and MW13 scenario have been independently determined, but agree well in the mass balance, it is reasonable to assume that this scenario adequately represents the conceptual model for the plume source term. This is further evaluated in Figure 6, which shows overlays of MSZ dissolved released and the plume balance forecast for different scenarios. The distributions of each forecast are mutually exclusive when using the MW13 data to estimate MSZ dissolved released (Fig. 6a) but overlap for both the total and BTEX-only plume balance when the combined MW7 and MW13 data are used (Fig. 6b, c). This suggests that the plume balance reflects the release of dissolved contaminants from MW7 and MW13 in combination, and that the amount of BTEX biodegradation attributable to electron acceptor consumption in Table 2, MSZ biodegraded, is a reasonable estimate of BTEX mass turnover for the period of study.

A plume-scale biodegradation rate for BTEX can be estimated using the equivalent mass of BTEX biodegraded in the plume, MSZ biodegraded, and the plume balance (inputs 11 and 12 in Table 2). The plume balance gives a total contaminant load of 1927 kg in the saturated zone, of which 1595 kg was biodegraded between December 1999 and July 2003. Assuming a first-order kinetic model, this mass loss corresponds to a first-order biodegradation rate, λ, of 0.0013/day or 0.49 a−1 for BTEX. This corresponds to a half-life (T0.5) of 1.41 years for the BTEX plume. It should be noted that the estimated BTEX biodegradation rate is a spatially and temporally averaged value, scaled over the volume and history of the BTEX plume, for the time period assessed. However, the biodegradation rate compares well with estimates obtained for BTEX plumes in other studies using different methods (Suarez and Rifai 1999; Wiedemeier et al. 1999).

The mass balance methodology developed in this study supports field-scale performance assessment of NA for dissolved petroleum fuel components in fractured, dual-porosity aquifers. It provides a robust quantitative framework to estimate site-specific mass balances for groundwater plumes using site investigation data to evaluate the distribution of dissolved contaminants in the saturated zone, estimate plume source terms, contaminant biodegradation and plume status (shrinking, stable or expanding), and identify potential limits on NA performance (e.g. oxidant supply). It is a useful screening tool to evaluate sites where MNA is being considered with other management options and also prior to undertaking more sophisticated numerical fate and transport modelling for site assessment. The methodology is transferable to, and can be further developed for, other dual-porosity aquifer settings when supported by appropriate data collection and site characterization. The approach can be also used to assess the NA of other organic contaminants in such aquifers and, when combined with probabilistic analysis, can address the inherent uncertainty in site investigation data.

For the case study, probabilistic analysis was particularly helpful in deducing the likely plume source term from different scenarios and predicting the long-term plume source composition as a probability distribution. This is advantageous in providing a site-specific estimate of the contaminant mass flux for forward modelling of plume development and NA performance. It also showed that errors or uncertainty in the conceptual model for the plume source term had a significant influence on the overall plume mass balance and interpretation of NA. However, the sensitivity analysis enabled the parameters that have the greatest control on predictions to be identified, in order to prioritize the site investigation to reduce the uncertainty in data inputs.

The biodegradation rate estimated for the BTEX plume using the probabilistic analysis is consistent with that reported for other BTEX plumes in the literature. Biodegradation using NO3 and SO4 supplied by advection and matrix diffusion appears to be a significant process limiting BTEX migration in the aquifer. Other biodegradation processes are limited by oxidant availability. Sorption to the aquifer matrix is also important for the attenuation of BTEX but less so for MTBE and TAME, which are more mobile in the groundwater.

The authors gratefully acknowledge the assistance of the site owner and Entec UK Ltd in the completion of this work. The authors also thank the two reviewers for their thoughtful comments, which helped improve the manuscript.

SFT: conceptualization (lead), data curation (equal), formal analysis (lead), funding acquisition (lead), investigation (equal), methodology (lead), project administration (lead), writing – original draft (lead), writing – review & editing (equal); MJS: conceptualization (supporting), data curation (equal), formal analysis (supporting), investigation (equal), methodology (equal), writing – original draft (supporting), writing – review & editing (supporting); SHB: conceptualization (equal), data curation (equal), formal analysis (supporting), funding acquisition (equal), investigation (equal), methodology (equal), project administration (equal), writing – original draft (supporting), writing – review & editing (equal); KHS: conceptualization (supporting), data curation (supporting), formal analysis (supporting), investigation (equal), methodology (equal), writing – original draft (supporting), writing – review & editing (supporting).

This work was funded by the Engineering and Physical Sciences Research Council and UK Environment Agency in the project ‘Processes controlling the natural attenuation of fuel hydrocarbons and MTBE in the Chalk aquifer’ (GR/R02849/01).

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

All data generated or analysed during this study are included in this published article (and, if present, its Supplementary information files).

This is an Open Access article distributed under the terms of the Creative Commons Attribution 4.0 License (http://creativecommons.org/licenses/by/4.0/)