Assessing human health risk of groundwater DNAPL contamination by quantifying the model structure uncertainty

Journal of Hydrology 584 (2020) 124690

Contents lists available at ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

Research papers

Assessing human health risk of groundwater DNAPL contamination by quantifying the model structure uncertainty

T

Yue Pan, Xiankui Zeng , Hongxia Xu, Yuanyuan Sun, Dong Wang, Jichun Wu ⁎

Key Laboratory of Surficial Geochemistry, Ministry of Education, School of Earth Sciences and Engineering, Nanjing University, Nanjing, China

ARTICLE INFO

ABSTRACT

This manuscript was handled by C. Corradini, Editor-in-Chief, with the assistance of Felipe de Barros, Associate Editor

The contamination of groundwater with dense non-aqueous phase liquids (DNAPL) can pose long-term risks to human health. Due to the complexity of geological environments and the properties of DNAPL, a key step in human health risk (HHR) assessment is the establishment of a reliable DNAPL transport model. In this study, we evaluate the impact of the model structure uncertainty on DNAPL transport modeling for HHR assessment. First, Gaussian process regression (GPR) was used to learn the structure error of the DNAPL transport model through statistical simulations, and then, Markov chain Monte Carlo (MCMC) simulation was used to calibrate the parameters of the physical model and GPR. On the basis of two case studies, namely, a sandbox experiment and a synthetic DNAPL transport model investigation, the results demonstrated that ignoring the model structure error will lead to deviations during parameter calibrations and model predictions. This is because the model parameters could be overfitted to compensate for the model structure error during the calibration process. Results of variance decomposition verified that the model structure error was the main uncertainty source in DNAPL transport modeling for both of the case studies. The DNAPL transport modeling and HHR assessment conducted by quantifying the uncertainties of the model parameters and structure simultaneously was found to be more reliable than that by calibrating only the model parameters. In addition, the HHR assessment of groundwater DNAPL contamination without considering the model structure error could lead to misleading depictions of the key risk zone. We propose that GPR is an effective technique to quantitatively describe the structure error of DNAPL transport models.

Keywords: Model structure uncertainty Gaussian process regression Human health risk DNAPL transport Markov chain Monte Carlo simulation

1. Introduction In many regions around the world, groundwater is an important source of water for humans and the ecological environment. Presently, along with rapid industrialization, the contamination of groundwater is threatening the water security of many communities (Burri et al., 2019; Rodríguez-Lado et al., 2013). The groundwater contaminants are diverse, such as heavy metal ions, pesticide residues, and industrial organic compounds (Lapworth et al., 2012; Teh et al., 2016; Teijon et al., 2010). Some of these contaminants are soluble in groundwater, while some are insoluble. The insoluble liquid contaminants are also called non-aqueous phase liquids (NAPLs). If the contamination is heavier than water, it is referred to as dense non-aqueous phase liquids (DNAPLs) (Kueper et al., 1993). Some DNAPLs, e.g., perchloroethylene (PCE), are toxic, and these can pose serious hazards to human health and the ecological environment (Henri et al., 2016). A key step in assessing human health risk (HHR) is to calculate the contaminant concentration that an individual is exposed to (de Barros

⁎

and Fiori, 2014; Maxwell and Kastenberg, 1999; Teh et al., 2016). It is difficult to monitor DNAPL transport in groundwater because it can pass through the unsaturated zone and downward by gravity in the saturated zone until it reaches the bottom of an aquifer (Kueper et al., 1993; Moreno-Barbero, 2005). Moreover, some DNAPLs can be transferred into the aqueous phase as a dissolution component, and these chemicals can be transported in groundwater as solutes. The DNAPL transport process involves the transport of DNAPL components and associated dissolved solutes in the aqueous phase, as well as DNAPL vapors in the gas phase (Falta et al., 1995; Lehame et al., 2003). Presently, numerical models are often used to quantitatively delineate the process of DNAPL transport in groundwater (Giraud et al., 2018; Parker and Park, 2004; Tatti et al., 2019). The reliability of DNAPL transport numerical models, especially in regard to the DNAPL component and dissolved solutes in the aqueous phase, has attracted widespread concerns among the groundwater/hydrology community. Quantifying the modeling uncertainty of a DNAPL transport model is necessary for assessing the distribution of DNAPL and HHR. Over the

Corresponding author. E-mail address: [email protected] (X. Zeng).

https://doi.org/10.1016/j.jhydrol.2020.124690 Received 25 December 2019; Received in revised form 5 February 2020; Accepted 13 February 2020 Available online 14 February 2020 0022-1694/ © 2020 Elsevier B.V. All rights reserved.

Journal of Hydrology 584 (2020) 124690

Y. Pan, et al.

past several decades, a large number of researchers have focused their attention on uncertainty quantifications of water resources and environmental models (Dai et al., 2019; Hamraz et al., 2015; Højberg and Refsgaard, 2005; Mustafa et al., 2018). Generally, this research has involved quantifying the modeling uncertainty caused by two types of sources, namely, the model parameters and structure. The model parameter uncertainty is concerning because model outputs are sensitive to some model parameters. Currently, several methods are commonly used to deal with parameter uncertainty, e.g., the least square regression (LSR) method (Doherty et al., 2010; Hill and Tiedeman, 2007), generalized likelihood uncertainty estimation (GLUE) (Beven and Binley, 1992; Hassan et al., 2008), and Markov chain Monte Carlo (MCMC) simulation (Hassan et al., 2009; Zeng et al., 2016a). In practice, the groundwater model structure, e.g., aquifer structure, boundary conditions, and chemical reactions of the contaminants, is always complex and unknown to the model users (Neuman and Wierenga, 2003; Troldborg et al., 2007). Because of the limited observational information, these model conditions are often simplified when establishing a DNAPL transport model (Refsgaard et al., 2012). The simplification of the model structure can cause a discrepancy between the model results and the actual hydrogeological conditions, which will lead to unreliable model outputs, e.g., DNAPL concentrations. Moreover, if the model structure uncertainty is not considered during model calibration, the model parameters could be overfitted to compensate for the model structure error (Gaganis, 2009). The study of model structure uncertainty has received attention in the past decade (Enemark et al., 2019). Currently, two types of methods are usually used to treat and quantify the model structure uncertainty, i.e., Bayesian model averaging (BMA) (Hoeting et al., 1999) and datadriven error modeling (Kennedy and O'Hagan, 2001). In BMA, the modeling uncertainty is quantified through assembling the outputs of all plausible model structures that are weighted based on their performance on reproducing observations (Raftery et al., 2005; Zeng et al., 2018). The predictions of BMA are more reliable than that of a single model because BMA incorporates the information on both model parameters and the structure. Troldborg et al. (2010) used four groundwater conceptual models to account for the uncertainty caused by the source zone and geological structure, and then, the mass discharge of a contaminated site was estimated using BMA. However, the implementation of BMA is limited to several aspects, e.g., the assignment of prior weights, building sufficient and independent model structures, and the accurate estimation of model weights (Cao et al., 2018; Lu et al., 2015). The data-driven error modeling method aims to build a statistical model to represent the model structure error through machine learning algorithms, e.g., Gaussian process regression (GPR). Kennedy and O'Hagan (2001) integrated Bayesian calibration and GPR to calibrate the parameter uncertainty and model inadequacy, and the prediction uncertainty was quantified by the inferred posterior mean and variance. Xu and Valocchi (2015) developed a Bayesian approach to calibrate the parameter uncertainty and Gaussian process error model simultaneously for a groundwater flow model, and they demonstrated that integrating the Gaussian process error model could effectively reduce the predictive bias. Because of the low demand for prior information on the model structure, the data-driven error modeling method shows great flexibility in the uncertainty quantification of the model structure (Xu et al., 2017a). The results of HHR assessment depend on the estimation of the contamination’s distribution. Thus, it is necessary to quantify the groundwater model structure error to obtain a reliable concentration distribution and scientifically sound HHR assessment. In this study, we integrate the model structure uncertainty quantification into the HHR assessment of groundwater DNAPL contamination, and the model structure error is described by a data-driven error modeling method (i.e., GPR). In addition, the impact of the model structure uncertainty on DNAPL transport modeling and the HHR assessment was also

evaluated through two case studies, a two-dimension DNAPL sandbox experiment and a synthetic three-dimension DNAPL transport model. The two case studies verified the validity of the GPR method. Moreover, the Bayesian theory based uncertainty quantification often requires tens of thousands of model executions to calibrate the model parameters and structure, and the resulting heavy computation time becomes an obstacle to HHR assessment. In order to overcome this burden, the sparse grid stochastic collocation technique (Bungartz and Dirnstorfer, 2003) was used to build cheap-to-run surrogates for DNAPL transport models in this study. Furthermore, the HHR assessment framework proposed in this study is independent of the specific groundwater DNAPL transport model, and thus, it is applicable to any water environmental model. The remainder of this paper is organized as follows. The framework of the HHR assessment, DNAPL transport model, MCMC simulation, GPR, variance decomposition, and sparse grid surrogate technique are briefly introduced in Section 2. Section 3 describes the settings of the two case studies, namely, a sandbox experiment and a synthetic model. Section 4 describes the implementation of the methods. The results and discussion of the two case studies are presented in Section 5, and the conclusions are given in Section 6. 2. Methodology 2.1. Human health risk assessment We applied the U.S. Environmental Protection Agency (EPA) Risk Assessment Guidance for Superfund (RAGS) to assess the HHR (EPA, 2009). The contaminant PCE is carcinogenic, and long-term intake can lead to human endocrine system and immune system dysfunction, liver and kidney function damage, and even death (Guyton et al., 2014). In the RAGS model, carcinogenic risk is quantified by estimating the increased probability of an individual suffering from cancer in a lifetime (Kumar et al., 2015; Saleh et al., 2019). The risk is calculated as follows:

Risk (x ) = 1

(1)

eCPF × ADD (x )

where CPF represents the cancer potency factor [kg d/mg], and ADD represents the average daily dose [mg/(kg d)]. x represents the control plane at the flow direction. As shown in Fig. 1, control planes represent the environmentally sensitive locations that are densely populated or contain human inhabitants (Henri et al., 2016). The EPA defines acceptable risk as values between 10−4 and 10−6 (EPA, 2009), and we set the acceptable risk value as 10−6 in this work for clarity. That is to say, if the risk is higher than 10−6, an individual would be subjected to unacceptable threats from the contaminant. In the RAGS model, the ADD is an important parameter that is related to many factors such as the contaminant concentration and individual exposure level. The exposure pathways are various, e.g., dermal sorption from showering, ingestion of tap water, inhalation of air (de Barros and Rubin, 2008; Kumar et al., 2015; Maxwell and Kastenberg, 1999). This work focuses on the impact of transport modeling of DNAPL in groundwater on the HHR assessment. Thus, only the ingestion of tap water was considered in this study, and the ADD was calculated as follows:

ADD (x ) = c¯ (x )

IR ED × EF BW AT

(2)

where IR represents the ingestion rate of water [L/d], BW represents the body weight [kg], ED represents the exposure duration [y], EF represents the exposure frequency [d/y], AT represents the averaging time [d], and c¯ represents the maximum average control plane concentrations of PCE (Maxwell et al., 1999; Siirila and Maxwell, 2012). The term c¯ is calculated as follows:

c¯ (x ) = max t>0

{ ED1

t + ED t

c (t , x ) d

}

(3)

where c(t, x) represents the flux-averaged concentration at the control 2

Journal of Hydrology 584 (2020) 124690

Y. Pan, et al.

Fig. 1. Overview of the DNAPL contamination and HHR targets.

plane position x. In this work, the parameters in Eqs. ((1) and (2)) (except c¯ ) were set to experience values.

krw =

Sw 1

Swr Swr

Sg

Sgr

1

Swr

m

(7) m

2.2. DNAPL transport simulation

krg =

T2VOC (TOUGH2 Volatile Organic Compound) is a module of the TOUGH2 family code. It is based on the integral finite-difference method and has been widely used to solve multiphase contaminant transport problems (Mccray and Falta, 1997; Niessner and Helmig, 2007). Generally, T2VOC assumes a multiphase system with the following three mass components (ω): air (a), water (w), and organic chemicals (n). These chemicals are distributed in one or more of the following three phases (γ): gas (g), aqueous (w), and NAPL (n). The balance equation for DANPL transport is written as follows (Falta et al., 1995):

where Sγ represents the saturation occupied by the γ phase, and Sγr represents the residual saturation occupied by the γ phase. The relationship between the capillary pressure and saturation is written as follows (Parker and Lenhard, 1987):

d dt

Pcgw =

M dVn = Vn

F ·nd

n

+

q dVn

(4)

Vn

n

Pcgn =

F

=

k

kr

( P

g)

1

Sg

Sw

Snr

1

Sg

Swr

Snr

(1

Sg

Swr 1

1 1

Snr )(1 Swr

Swr Sw Sw )

wg nw

Sw 1

Swr Swr

Snr Snr

n n 1

1 n

n n 1

1

1 n

1

(9) wg gn

1 n

n

n 1

1

(10)

where Pcgn represents the gas-DNAPL capillary pressure [Pa], Pcgw represents the gas-water capillary pressure [Pa], αgn represents the gasPCE inlet pressure reciprocal, αnw represents the PCE-water inlet pressure reciprocal, and n represents the fitting parameter. In T2VOC, hydrodynamic dispersion is neglected and molecular diffusion is simulated for the gas phase only (Basu et al., 2008; Falta et al., 1995). Mass transfer is depicted by advection in all phases (Bennett et al., 1998; James and Oldenburg, 1997). Fluid phases are assumed to be in local chemical equilibrium for each element, so that interphase mass transfer among all phases can be described by equilibrium partitioning relations (Gudbjerg et al., 2005; Neriah and Paster, 2019). The water in the NAPL phase is neglected. The model accounts for (a) interphase mass transfer, (b) adsorption of NAPLs to the solid phase, and (c) decay of NAPLs by biodegradation (Falta et al., 1995). The relevant mechanisms of interphase mass transfer for NAPL include NAPL evaporation, dissolution, and gaseous-aqueous mass transfer (Dafny, 2017). In this work, we only considered DNAPL transport in a saturated aquifer, so we were only concerned with the DNAPL component and DNAPL dissolution transport in the aqueous phase. For the DNAPL dissolution transport in the aqueous phase, T2VOC produces a velocitydependent numerical dispersion. It has a comparable effect to hydrodynamic dispersion (Dafny, 2017). T2VOC has been successfully implemented in the modeling of DNAPL transport and dissolution processes in many research projects (Bennett et al., 1998; James and Oldenburg, 1997).

(5)

where k represents the absolute permeability [m2], krγ represents the relative permeability of phase γ [m2], ργ represents the γ phase density [kg/m3], ηγ represents the γ phase dynamic viscosity [cP], Pγ is the fluid pressure in phase γ [Pa], and g is the gravitational acceleration vector [m/s2]. Moreover, the DNAPL transport process is also affected by the relative permeability and capillary pressure between the phases. In this work, we apply the Stone function to describe the relationship between the relative permeability and saturation. The Parker and Lenhard function can be used to illustrate the relationship between the capillary pressure and saturation. The relationship between the relative permeability and saturation is written as follows (Stone, 1970):

krn =

gn

Sw + Sn Swr Snr 1 Swr Snr

Sw + Sn Swr Snr 1 Swr Snr

where Mω represents the mass of component ω in the unit pore volume [kg], Vn represents the arbitrary flow region, Γn represents the surface area [m2], Fω represents the flux of component ω into volume Vn, n represents the inward unit normal vector, and qω represents the rate of heat generation in the unit pore volume. The motion of each phase of the fluid is under the influence of pressure, gravity, effective permeability, and phase viscosity following the multiphase extension of Darcy’s law:

F =

wg

(8)

Snr Snr m

(6) 3

Journal of Hydrology 584 (2020) 124690

Y. Pan, et al.

2.3. Markov chain Monte Carlo simulation

Ci, j = q (mi , mj )

The relationship between the model f and the observations D can be written as:

=

where ε represents the residuals errors. Based on Bayes’ theorem, the posterior distribution p(θ|D) of model parameters θ is described as:

p ( |D ) =

L ( |D ) p ( ) L ( |D ) p ( ) d

(12)

where L(θ|D) represents the joint likelihood function of observation data D, and parameters θ, θ[θ1, θ2,…, θd] represent d-dimensional model parameters. p(θ) represents the prior distribution of parameters θ, and it is usually specified by the user according to expert experience or pilot experiments. Because the likelihood function of nonlinear model f is very complicated, it is intractable to derive the analytical expression of the posterior distribution p(θ|D). Yet the MCMC method provides an effective way to construct the posterior distribution via Markov chain samples (Haario et al., 2001). It builds a single or parallel Markov chains to search for the probability space of parameters. Markov chain samples evolve from the prior parameter distribution to posterior distribution gradually by incorporating all information. In this work, the DiffeRential Evolution Adaptive Metropolis (DREAMZS) algorithm was used to perform MCMC simulation, and this algorithm has been widely used to calibrate model parameters in high-dimensional and complex hydrology problems (Wöhling et al., 2012; Zeng et al., 2016a). The details of this algorithm can be found in earlier studies (Laloy and Vrugt, 2012; Vrugt et al., 2009). The Gaussian likelihood function is often chosen as the expression of the joint likelihood:

L ( |D) = (2 )

N /2 |

|

1/2 exp

1 (D 2

p( ,

1 (D

f ( ))

p (D | ,

log(p (D | , 1 = (D 2

(i, j = 1, 2, ......,N )

(15)

(16)

)p ( )p ( )

)) µ)T C 1 (D

f( )

f( )

1 log |C| 2

µ)

N log(2 ) 2 (17)

The same as described in Section 2.3, we applied MCMC simulation to estimate the posterior distributions of θ and Φ jointly. After the posterior samples of θ and Φ are generated, the predictions of f(θ) and b (m,Φ) for a new input x* (denoted as f*(θ) and b*(m,Φ), respectively) can be obtained. The probability distribution of b*(m, Φ) obeys the normal distribution N (b¯ , Cb (b )) , and

b¯ = E (b |D

f ( ),

C

) = µ + C T C 1 (D

f( )

µ)

(18) (19)

T C 1C

Here, the asterisk superscript represents the variable associated with x*, μ* means μ(x*,Φ), C** is the prior covariance matrix, C* is calculated as Ci, j = k (x i , xj ) , and C** is calculated as Ci, j = k (x i , x j ) . As a result, the predicted output at x* is expressed as follows:

(13)

(20)

D =f +b +

In summary, the implementation of the GPR based model structure uncertainty analysis can be summarized as follows: (i) define the prior distribution of the unknown hyper-parameters Φ and model parameters θ; (ii) use MCMC simulation to obtain the posterior distribution of Φ and θ, and then, calculate f*, C, C*, and C** based on the obtained posterior samples and Eqs. ((15) and (17)); (iii) use multivariate normal sampling to attain b*i and δ*i, and then, calculate the predicted D*i by using Eqs. (18)–(20).

2.4. Gaussian process regression As mentioned in Section 2.3, the residuals error ε represents the deviation between model output and observation. Generally, this deviation is derived from three sources including observation error, biased model parameters θ, and defected model structure. The uncertainty derived from observation and model parameters is well calibrated through assuming the error structure of residuals and Bayesian inference respectively in MCMC. In order to quantify the model structure error explicitly, the relationship between observation D and model f is expressed as:

)+

2

+

where p(D|θ,Φ) represents the likelihood function; here, we chose the logarithmic form to describe:

where N represents the number of observations, Σ represents the residual’s covariance matrix, and f(θ) represents the output of a nonlinear model f (e.g., DNAPL transport model).

D = f ( ) + b (m ,

|D )

Cb (b ) = C f ( ))T

mj )

2

where Ψ represents the indicator function; when i = j, Ψ equals one, and otherwise, it equals zero. λ, σ2, and σδ2 are three hyper-parameters. λ represents the characteristic length scale, σ2 controls the marginal variance of b(m), and σδ2 describes the errors except for the model structure error, such as the measurement error. Generally, the prior mean μ is assumed as zero, and the uncorrelated errors obey independent and identical normal distributions. The model parameters θ and Gaussian process model hyper-parameters Φ(λ, σ2, σδ2) are jointly calibrated through Bayesian inference. In addition, the prior probability is specified as p(θ, Φ) = p(θ)p(Φ) and the prior distributions of θ and Φ are assumed to be independent. According to Bayes’ theorem, the posterior probability is formed as:

(11)

D=f( )+

mj )T (mi

(m i

2 exp

2.5. Variance decomposition The variance decomposition method can be used to quantify the contributions of uncertainty sources (e.g., model parameter and structure error, measurement error) to the total predictive uncertainty (Kumar et al., 2019; Xu et al., 2017b). According to the law of total variance, the variance of the prediction D* can be calculated as:

(14)

where b(m, Φ) represents the model structure error, m represents the inputs of the model, Φ represents tuning parameters, and δ represents the observation error. The data-driven machine learning provides an effective way to characterize model structure errors, and the GPR method is used in this study (Kennedy and O'Hagan, 2001; Xu et al., 2014). The prior distribution of the Gaussian process function b(m,Φ) obeys a multivariate Gaussian distribution N(μ,C). The Gaussian process is specified by a mean function μ(m) = E[b(m)] and a covariance function q(m, m’) = E [(b(m) - μ(m))(b(m’) - μ(m’))]. In this work, the squared exponential covariance function is applied:

V[D ] = E + E [Vb

|

[E

,b ,

, | ,b

[V |

,b ,

[D | , b ,

[D | , b ]]] + V [Eb ,

]] , |

[D | ]]

(21)

where V represents the variance, and E represents the expectation. The other notations are the same as described in Sections 2.3 and 2.4. The left side of Eq. (21) is the variance of D*i, i = 1, … , M, which means the total predictive uncertainty, and the posterior predictive samples are obtained by Eq. (20). The first term of the right side of Eq. (21) represents the measurement uncertainty. The second and third terms 4

Journal of Hydrology 584 (2020) 124690

Y. Pan, et al.

represent the uncertainty of the model structure and parameters, respectively. More details about variance decomposition and the calculation process can be found in Xu et al. (2017a), Dai and Ye (2015), and Hoeting et al. (1999). 2.6. Sparse grid surrogate technique Because of the massive number of computations required during the model calibration process, we implemented sparse grid stochastic collocation to construct a surrogate model, thereby replacing the original DNAPL transport numerical model. Sparse grid (SG) surrogate models are built from a series of Lagrange interpolants. For a nonlinear function, e.g., the DNAPL transport model f(θ), the one-dimensional hierarchical Lagrange interpolant (Lr) is defined as: L

Lr L (f ( )) =

L

Lr i (f ( )) = i=0

mi

i=0 j=0

coij

i j(

), i = 0, …, L;

(22)

where L represents the resolution level, ΔLr represents the incremental interpolation operator, mi represents the number of SG nodes for level i, coji represents the corresponding interpolation coefficient, and φii(θ) represents the interpolation basis function. For a multi-dimension model f(θ), θ[θ1, θ2,…, θd], the multi-dimensional hierarchical interpolation (mLr) is performed similarly to that of a one-dimensional function, and the SG nodes are generated in the manner of the Smolyak rule (Smolyak, 1963); the required SG nodes are significantly less than those of a full-tensor product grid. In order to further reduce the number of SG nodes, i.e., the costs of building the surrogate, some adaptive SG node refinement strategies have been developed, e.g., local adaptive SG (Ma and Zabaras, 2009b; Zeng et al., 2016b), dimensional adaptive SG (Gerstner and Griebel, 2003; Zeng et al., 2012), local-dimensional adaptive SG (Jakeman and Roberts, 2011; Jakeman and Roberts, 2013). In this study, local-dimensional adaptive refinement was used. More details about the construction of the adaptive SG surrogate, e.g., calculating SG coefficients, defining the SG basis function, can be found in Ma and Zabaras (2009a), Pflüger et al. (2010), and Zhang et al. (2013). 3. Descriptions of the two case studies Two case studies were used in this study, namely, a sandbox DNAPL transport experiment and a synthetic groundwater DNAPL transport model investigation. For each case study, the true model structure and the defective model structure caused by oversimplification or limited information, respectively, are both described. 3.1. DNAPL transport sandbox experiment Fig. 2. Sketch map of the sandbox experiment (a is the structure of the sandbox aquifer, b is the LTM saturation observations, and c is the structure of the conceptual model used to simulate DNAPL transport).

In this section, a two-dimensional sandbox experiment model is used to simulate DNAPL transport (Fig. 2(a)). The sandbox was 60.0 cm in length, 45.0 cm in height, and 1.6 cm in thickness. The left and right boundaries were set to a Neumann boundary, and the flow velocity was kept at 1.5 m/d by use of a peristaltic pump. The sandbox was filled with quartz with a particle size of 0.60–0.85 mm (sand A). Moreover, five rectangular lenses consisting of low permeability quartz with a particle size of 0.18–0.22 mm (sand B) were placed in the sandbox. The injection point was located above the uppermost lenses at x = 30.0 cm, y = 40.0 cm, and z = 0.8 cm. The PCE injection rate was 0.5 mL/min, and the total amount of PCE in the sandbox was 40 mL. We monitored PCE saturation in the sandbox for the observations, and these data were measured by the light transmission method (LTM). Based on the LTM image (Fig. 2(b)), we can see that the PCE bypassed the four lenses and was transported to the bottom of the box. The experiment observation period was 80 min, and the sampling interval was 5 min. We used the observations before 60th min to calibrate the numerical transport model and used the latter

observations (from 61th min to 80th) to validate this model. In practice, it is difficult to accurately characterize the structure of a groundwater aquifer. This study assumes that the established groundwater conceptual model deviates from the true aquifer structure, i.e., the left lense was missed by the conceptual model due to inadequate geological drilling (Fig. 2(c)), and this defective model structure was used to simulate the transport of DNAPL. 3.2. Synthetic DNAPL transport model The synthetic three-dimensional groundwater model had three aquifers, as shown in Fig. 3(a). The model domain was 1000 m long, 100 m wide, and 20 m thick. We set three sands for the three aquifers, respectively, e.g., sand 1, sand 2, and sand 3. The permeability of three 5

Journal of Hydrology 584 (2020) 124690

Y. Pan, et al.

Fig. 3. (a) The true model structure, and (b) the simplified model structure.

sands was ranked as follows: sand 1 > sand 2 > sand 3. As shown in Fig. 3, the upper (sand 1) and lower (sand 2) aquifers were separated by an uneven and thin aquitard (sand 3), e.g., a mudstone aquitard. The distribution of the middle model layer was discontinuous, and the thickness was not uniform. For each layer, the porous medium was regarded as homogenous and isotropic. The left and right boundaries were constant water heads with hydraulic gradients of 0.01. The other model boundaries were all set to impermeable. When assessing the HHR, we assumed a scenario that a total of 85 kg PCE had been poured into the aquifers 50 days earlier. The pour site was around x = 50 m, y = 50 m at the top of the first layer. Therefore, the aquifers had already been contaminated by PCE before the HHR assessment. Fig. 4 exhibits the initial PCE concentration distribution in aquifer. In actual situations, we are unable to attain sufficient information on the distribution of groundwater aquifers, especially for thin and discontinuous layers. As a result, a defective model structure is often used to represent such groundwater systems. In this case study, we assumed that the all aquifers were simplified to uniform layers because of the limited borehole information on the thin aquitard. This simplified model structure had the same model domain and boundary conditions with that of the true model structure, as presented in Fig. 3(b). In addition, the thicknesses of the three even model layers were set to 9.0 m, 2.0 m, and 9.0 m, respectively, from top to bottom. In this case study, the observation data were generated by running the “true” model. Two observation wells were placed in the aquifer to monitor the concentration of PCE, and these were equipped with screens at the first layer (obs1) and the third layer (obs2), as shown in Fig. 4(a) and (c). The concentrations in the wells were recorded once a month over five years. Thus, for each well, there were 60 concentration records. As terminology, we use obs1-1 to obs1-60 to present the 60 observations monitored in obs1, and similar codes for the data for obs2. The total number of observations was 120. These observations were

eroded with Gaussian random error with a mean of 0 and variance of 3% in relation to the observed data, and then, the results were used as conditional data. For each observation well, the conditional data were divided into two parts; the former 50 observations were used for model calibration, and the latter 10 observations were used for model validation. 4. Implementation of the methodology The DNAPL transport modeling based on the defective model structure was conducted for each groundwater DNAPL transport case (the sandbox experiment and the synthetic model). For each case study, the DNAPL transport model’s parameter uncertainty was quantified by MCMC, and the GPR method was used for quantifying the model structure uncertainty. The variance decomposition method was used to quantify different uncertainty sources’ contributions to the total predictive uncertainty. This uncertainty quantification method was verified by the two groundwater examples, and then, the synthetic model was implemented to conduct the HHR assessment. 4.1. Uncertainty quantification of the sandbox DNAPL transport model In this case study, the DNAPL’s saturation distribution of the sandbox experiment was predicted in two ways; the first was to quantify only the parameter uncertainty of the DNAPL transport model, and the second was to quantify both the model parameter and structure uncertainty. According to Section 3.1, we applied T2VOC to simulate the PCE transport in the sandbox. The simulation was generalized into a twodimensional x-z section. The x and z directions were uniformly discretized into 30 and 45 grids, respectively. For the uncertainty quantification of model parameters, five important model parameters were regarded as random variables, including the permeability of sand A (ka)

Fig. 4. The initial PCE concentration distributions in the aquifers; (a), (b), and (c) represent the x-y sections at z = 12 m, z = 10 m, and z = 4 m, respectively. 6

Journal of Hydrology 584 (2020) 124690

Y. Pan, et al.

and sand B (kb), the residual NAPL phase saturations of sand A (Snra) and sand B (Snrb), and the viscosity of PCE (ηn); the meaning of these parameters can be found in Section 2.2. The prior distributions of these parameters were defined as uniform distributions, and the distribution ranges were defined as [2.0 × 10−10, 5.0 × 10−10] for ka, [8.0 × 10−12, 3.5 × 10−11] for kb, [0.05, 0.12] for Snra, [0.05, 0.12] for Snrb, and [0.6, 1.2] for ηn. The MCMC simulations were used to infer posterior parameters’ distributions through Bayes’ theory, and DREAMZS was used as the searching algorithm. Three parallel chains were applied to conduct the MCMC simulations, and the iteration length of each chain was set to 12,000 with a burn-in period of 3000. After the MCMC simulations were completed, the generated posterior distributions of the model parameters were used for predicting the DNAPL saturation distribution. For the uncertainty quantification of the model structure, in addition to quantifying the model parameters’ uncertainty, the GPR method was used to describe the model structure error. The random model parameters and their prior distributions were the same as those described above. In addition, the squared exponential covariance function was used to build the GPR error model in this study (Eq. (15)), including three hyper-parameters (i.e., the characteristic length scale λ, the standard deviation σ, and the standard deviation of measurement error σδ). The prior distribution was defined as a gamma distribution G(5.0, 5.0) for λ, and the range was set to [0.2, 2.0]. An exponential distribution with a mean of 0.25 was defined as the prior distribution of σ, and the range was set to [0.0002, 2.0]. The σδ follows a uniform prior distribution on [0.002, 0.5]. The model parameters and GPR model’s hyper-parameters were jointly calibrated through MCMC simulations. Moreover, the MCMC simulations were implemented with the same setting to that used for the parameter uncertainty quantification above. After the posterior predictive samples were obtained by the MCMC method, we conduct variance decomposition to quantify the contributions of model parameters, model structure, and measurement error to the predictive uncertainty. Generally, one forward run of this DNAPL transport model (by T2VOC code) takes more than 100 s on a personal computer (PC, 8-core i7-6700 3.40 GHz processors). In order to overcome the problem of heavy time consumption caused by the large number of model executions in MCMC simulations, i.e., 72,000 forward model runs would have been required for this case study, the SG technique was used to build surrogate models. Specifically, we built an SG surrogate model to emulate the response relationship between the model output (PCE saturation) and model input (five random model parameters). The original DNAPL transport model was executed for 4801 times to build this surrogate model, i.e., the cost of building this SG surrogate model was 4801 forward model runs. After the surrogate model was constructed, we tested the surrogate accuracy by comparing the outputs (PCE saturation) of the surrogate model and the original DNAPL transport model with 150 random model parameter sets. Fig. 5 illustrates the comparison between them; the RMSE (root mean square error) between the outputs of these two models was approximately 1.0 × 10−3. In addition, the R2 (coefficient of determination) between the outputs of these two models was more than 0.997. Obviously, the SG surrogate model showed sufficient accuracy in emulating the response relationship between PCE saturation and model parameters. The performance of each DNAPL transport modeling was quantitatively measured by the log-score metric (Hoeting et al., 1999), which was defined as

model parameters) in this study. A smaller log-score value indicates a better predictive performance. 4.2. Uncertainty quantification of the synthetic DNAPL transport model In this case study, the HHR of synthetic groundwater DNAPL contamination was assessed with two approaches, i.e., the DNAPL transport was modeled by only considering the model parameter uncertainty and by considering both the model parameter and structure uncertainty, respectively. The metric of log-score was used to measure the performance of each DNAPL transport modeling effort. According to Section 3.2, we built a three-dimensional PCE transport model by using T2VOC. The x, y, and z directions were uniformly discretized into 30, 8, and 11 grids, respectively. The procedures for applying the MCMC and GPR methods were similar to those in Section 4.1. The unknown model parameters in the synthetic model included the permeability and PCE-water inlet pressure reciprocal of the two model layers (sand 1 and sand 2); the meanings of these parameters are described in Section 2.2. The prior ranges of these four model parameters were specified as [1.0 × 10−10, 5.0 × 10−10] for the permeability of sand 1 (k1), [40, 60] for the PCE-water inlet pressure reciprocal of sand 1 (αnw1), [8.0 × 10−11, 2.0 × 10−10] for the permeability of sand 1 (k2), and [20, 40] for the PCE-water inlet pressure reciprocal of sand 2 (αnw2). In addition, the prior uniform distributions were assumed for these parameters. Three parallel Markov chains were set in the MCMC simulations, and the iteration length and burn-in length were set to 9000 and 3000, respectively, for each chain. The three hyper-parameters of the GPR model had the same prior probability functions to those described in Section 4.1, and the distribution ranges were defined as [0.2, 4.0] for λ, [0.0002, 0.2] for σ, and [0.002, 0.5] for σδ. Moreover, the setting of MCMC simulations was the same as that used for the parameter uncertainty analysis above. After that, the DNAPL concentration distribution was generated based on the obtained parameters’ posterior distributions and GPR, and then, these data were used for the HHR assessment. The variance decomposition method was used to quantify the contributions of uncertainty sources (e.g., model parameters, model structure, and measurement error) to the predictive uncertainty. For the two HHR assessments (i.e., with and without GPR, respectively) of this synthetic case study, some parameters of the HHR model were set to constants according to earlier research (de Barros et al., 2009; Henri et al., 2016), e.g., the cancer potency factor (CPF) was 0.0002 [kg d/mg], the ingestion rate (IR) was 1.6 [L/d], the body weight (BW) was 70 [kg], the exposure duration (ED) was 30 [year], the exposure frequency (EF) was 350 [d/year], the average time of the expected lifetime (AT) was 25,550 [d]. Moreover, 27 control planes were equally spaced along the × direction. The 1st control plane was located at x = 67 m, which was close to the contaminant source zone, and the interval of each control plane was 33 m. The HHRs were assessed at these planes based on Eqs. ((1) and (2)). A complete MCMC simulation would have demanded 27,000 iterations to generate posterior parameter samples for each HHR assessment. As a result, a total of 54,000 DNAPL transport model executions would have been required for this case study. One forward run of this DNAPL transport model (by T2VOC code) would take more than 2 min on a PC (8-core i7-6700 3.40 GHz processors). Similar to Section 4.1, we constructed a cheap-to-run surrogate model based on the SG technique for the original time-consuming DNAPL transport model to overcome the computational burden. We built the response relationship between the PCE concentration and the four model parameters. The original DNAPL transport model was executed a total of 2206 times at the SG nodes. The results of the surrogate accuracy tests are displayed in Fig. 6. The RMSE between the outputs of these two models (PCE concentration) was always smaller than 1.0 × 10−4. In addition, the R2 between the outputs of these two models was always higher than 0.9999. As a result, the SG surrogate model had sufficient accuracy to

n

log

score =

log p ( i |D) i=1

(23)

where Δi means a predictive quantity (e.g., DNAPL concentration), p (Δi|D) is the posterior predictive probability of Δi (i = 1, … n) under observation D, and n means the number of validation observations (or 7

Journal of Hydrology 584 (2020) 124690

Y. Pan, et al.

Fig. 5. Comparisons between the surrogate model and original sandbox DNAPL transport model at some observation points.

Fig. 6. Comparisons between the surrogate model and original synthetic DNAPL transport model at some observation points.

replace the original model.

error (i.e., the simplification of left lense) was considered or not. Then, the inferred DNAPL saturation predictions were evaluated and compared at some observation points (those not used for model calibration).

5. Results and discussion In this section, the DNAPL transport modeling and/or HHR assessments are evaluated in two ways, i.e., according to whether the model structure uncertainty was quantified or not for each case study. Moreover, the contributions of DNAPL transport modeling uncertainty sources to model predictions (e.g., DNAPL concentrations) were quantified based on variance decomposition.

5.1.1. Posterior distributions of model parameters The blue lines in Fig. 7 present the posterior distributions of model parameters inferred by MCMC simulations without considering the model structure error. The horizontal axis denotes the parameter’s prior range (hereafter the same). Clearly, the high probability density region was distributed in the middle of the prior range for each parameter. In addition, due to the sandbox experiment being a real case, the true values of the model parameters were not available. Thus, it was infeasible to evaluate the accuracy of posterior parameter distributions. The red lines in Fig. 7 present the posterior distributions of five model parameters and three hyper-parameters by incorporating the

5.1. DNAPL transport modeling of the sandbox experiment The DNAPL transport modeling of the sandbox experiment was conducted in two ways, i.e., according to whether the model structure 8

Journal of Hydrology 584 (2020) 124690

Y. Pan, et al.

Fig. 7. Model parameters’ posterior distributions for the sandbox experiment case study; the red and blue lines represent the DNAPL transport that was simulated with and without the uncertainty quantification of the model structure, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

GPR error model into MCMC simulations. We found that these five model parameters’ posterior distributions were apparently different from those inferred by only quantifying the model parameter uncertainty (the blue line). In addition, some parameters’ posterior distributions (e.g., ka, kb, Snra) displayed the characteristics of a multimodal distribution, which indicates that the parameter space was very complex. Therefore, model structure error had a significant impact on the parameter calibration. Hence, ignoring the model structure error may lead to deviations in parameter calibrations; this is because the model parameters could be overfitted to compensate for the model structure error. Moreover, the plots of Fig. 7(f)–(h) illustrate that the posterior distributions of three hyper-parameters converged to narrower ranges compared to their prior distributions.

uncertainty. The red lines of Fig. 8 present the DNAPL saturation predictions inferred by incorporating the GPR results. We found that all of the validation observations were predicted with high accuracy (i.e., high predictive probability). Compared with the blue lines in Fig. 8, it was obvious that the performance of DNAPL saturation predictions was improved significantly by quantifying the DNAPL transport model’s structural uncertainty, except for one observation (x = 37 cm, y = 21 cm, t = 80 min). Moreover, the log-score was equal to −9.98 for all validation observations, which was much smaller than that without considering the model structure error (log-score = −0.79). Therefore, it was effective to improve the predictive performance of the DNAPL transport model by quantifying the model structure error.

5.1.2. DNAPL transport modeling In order to evaluate the accuracy of DNAPL saturation predictions, the inferred posterior distributions of DNAPL saturation at the validation observations were compared with their true values (the vertical black dash dot line). The blue lines of Fig. 8 exhibit the predictions of DNAPL saturation at six validation observations without considering the model structure error. We found that only one observation (x = 37 cm, y = 21 cm, t = 80 min) was predicted with high accuracy, i.e., the predicted distribution had a high probability density at the observed saturation. The other five observations were predicted with apparent deviations, i.e., the predicted distribution had a low probability density at these observed saturations. In addition, the prediction of DNAPL saturation showed similar performance at other validation observations, and these data are not shown here because of the limited space; the log-score was equal to −0.79 for all of the validation observations. Therefore, the DNAPL saturation distribution cannot be reliably predicted by only accounting for the model parameter

5.1.3. Variance decomposition for the sandbox experiment Fig. 9 exhibits the contribution fractions of uncertainty sources (model parameter, model structure, and measurement error) to the DNAPL transport predictive uncertainty at 20 observations. Clearly, the model structure error was the major uncertainty source in the sandbox experiment model, and it accounted for more than 80% of the total uncertainty. The measurement error contributed generally less than 1% to the predictive uncertainty. The model parameters contributed less than 5% to the predictive uncertainty for most observations, except at O-5 and O-12 where the fractions were larger than 10%. Noteworthy, according to Fig. 8, we can see that the prediction distributions at O-5 and O-12 showed better performance than the other observations when ignoring the model structure error. This was because model parameters made more contributions to the predictive uncertainty for these two observations, while ignoring the model structure error had less of a negative impact on model predictions. Moreover, the model structure error could be partly compensated by calibrating the 9

Journal of Hydrology 584 (2020) 124690

Y. Pan, et al.

Fig. 8. Predictive distributions of DNAPL saturation; the red and blue lines represent the DNAPL transport that was simulated with and without the uncertainty quantification of the model structure, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

model parameters.

densities at the corresponding true values. The log-score was equal to 4.17 for all of the model parameters. Thus, these model parameters cannot be calibrated reliably by using the defective model structure and only considering the model parameter’s uncertainty. As GPR was used to describe the model structure error, the seven parameters, which included four model parameters and three hyperparameters, were calibrated by MCMC simulations. The posterior distributions of these seven parameters are described in Fig. 11. Compared with Fig. 10, we found that the parameters’ posterior distributions were closer to their true values, and these posterior distributions showed relatively higher probability densities at the corresponding true values. The log-score was equal to −17.12 for all of the model parameters, which was smaller than that without considering the model structure error (log-score = 4.17). Therefore, ignoring the model structure error could lead to deviations in the model parameters’ calibrations, and these parameters can be better identified by quantifying the model structure uncertainty. In addition, the posterior distributions of the three hyper-parameters converged to much narrower ranges than their

5.2. HHR assessments of the synthetic case study The DNAPL transport modeling and HHR assessment for the synthetic case were conducted in two ways, i.e., according to whether the model structure error (i.e., the simplification of the middle aquifer) was quantified or not. Then, the inferred DNAPL predictions and HHRs were evaluated and compared. 5.2.1. Posterior distributions of model parameters Fig. 10 presents the posterior distributions of the model parameters by MCMC simulations without incorporating the GPR results. The vertical black dash line marks the true parameter value (hereafter the same). The horizontal axis denotes the parameter’s prior range. Clearly, all of the model parameters (except for αnw1) were calibrated with apparent deviations by the corresponding posterior distributions, and these parameters’ posterior distributions showed very low probability

Fig. 9. Contribution fractions of uncertainty sources (model parameter, model structure, and measurement error) to the predictive uncertainty at 20 observations for the sandbox experiment. 10

Journal of Hydrology 584 (2020) 124690

Y. Pan, et al.

simulations, and these predictive distributions always deviated significantly from the true values. In addition, the predictions of DNAPL concentrations at other observations displayed similar features to Fig. 12. The log-score was equal to 0 for all of the validation observations (0*log0 = 0), which means all observations were predicted with a zero probability density at their true values. Therefore, it is inadvisable to describe the DNAPL concentration based on a DNAPL transport model with structure error, although parameter uncertainty quantification is performed. The parameter calibration cannot fully compensate for the uncertainty caused by the model structure error. Fig. 13 presents the obtained predictive distributions of PCE concentrations when considering the model structure uncertainty. The results clearly showed that all of these predictive distributions were concentrated at the true values of observations. Compared with Fig. 12, the accuracy of these predictions was improved significantly by incorporating the GPR error model, and these predictive distributions showed apparently higher probability densities at the corresponding true values. The log-score was equal to −14.67 for all of the validation observations, which was smaller than that without quantifying the model structure uncertainty (log-score = 0). Thus, it is more reliable to simulate the transport of DNAPL contamination by quantifying the uncertainties of model parameters and structure simultaneously than that by calibrating only model parameters. In addition, the GPR model was found to be an effective technique to describe the model structure error. Generally, the predictive distributions with the uncertainty quantification of the model structure had wider ranges than those with only quantifying the model parameters’ uncertainty, and some distributions even had a minimum value less than zero. This was because the addition of the model structure error term brought about an increase in the prediction variance.

Fig. 10. Posterior distributions of the model parameters for the synthetic case without considering the model structure uncertainty.

corresponding prior distributions. 5.2.2. DNAPL transport modeling Fig. 12 shows the obtained predictive distributions of PCE concentrations without considering the model structure error. The performance of DNAPL transport modeling was evaluated at validation observations that were not used for model calibration. Obs1-1 denotes the 1st observation well at the 1st observation moment (hereafter the same). Because of limited space, only six typical observations and their predictions are displayed in this case study. The results clearly showed that all of the observations could not be effectively predicted by the

Fig. 11. Posterior distributions of model parameters and hyper-parameters for the synthetic case study considering the model structure uncertainty. 11

Journal of Hydrology 584 (2020) 124690

Y. Pan, et al.

Fig. 12. DNAPL concentrations’ distributions predicted without considering the model structure uncertainty for the synthetic case.

Fig. 13. DNAPL concentrations’ distributions predicted with the uncertainty quantification of the model structure for the synthetic case.

5.2.3. Variance decomposition for the synthetic model Fig. 14 illustrates the contribution fractions of uncertainty sources to the predictive uncertainty of the synthetic model at 20 observations (Obs1-51–60, Obs2-51–60). We found that the model structure error accounted for more than 90% of the total uncertainty for all observations of both of the two wells. As a result, the model structure error dominated the predictive uncertainty of the synthetic model, and the model parameters’ calibration had a limited effect on model predictions. This is why the model parameters and predictions were calibrated with significant deviations by ignoring the model structure error, as shown by Fig. 12.

5.2.4. HHR assessment The predictive distributions of groundwater contaminant concentrations obtained were used for the HHR assessment. The increased probability of an individual suffering from cancer in a lifetime was calculated as shown in Eq. (1). According to the RAGS (EPA, 2009), the acceptable risk value was defined as 10−6 in this work. Fig. 15 presents the probability densities of HHR assessed by using the PCE concentration distributions without addressing the model structure uncertainty. The vertical black line represents the acceptable risk value (10 −6). The control planes were numbered from 1 to 27 in the order of distance from the contaminant source; plane 1 was closest to the source, and plane 27 was farthest from the source. Because of 12

Journal of Hydrology 584 (2020) 124690

Y. Pan, et al.

Fig. 14. Contribution fractions of uncertainty sources (model parameter, model structure, and measurement error) to the predictive uncertainty at 20 observations for the synthetic model.

space limitations, here we only display six of the control planes, e.g., from 11th to 16th. Obviously, the HHR decreased as the distance of the control plane from the contaminant source increased. In addition, for the control plane that was closer to the contaminant source than the 13th control plane, the HHR value with the highest probability density was always higher than 10−6. This indicates that necessary actions must be taken to prevent this contamination. For the control plane that was farther away from the source than the 13th control plane, the HHR value decreased gradually and was less than 10−6. Therefore, the 13th control plane was called the “critical control plane” based on this HHR assessment. Fig. 16 exhibits the probability densities of HHRs assessed with the GPR results. Similar to the tendency shown in Fig. 15, the HHR decreased as the distance of the control plane from the contaminant source increased. However, the HHRs assessed by incorporating GPR results were usually slightly higher than those obtained by neglecting the model structure error. Moreover, the ranges of distributions were

much wider, which indicates that HHR assessed with considering the model structure error had a larger variance. Based on the assessment of HHR, the critical control plane was the 15th control plane. Compared with the results of Fig. 15, the critical control plane had changed from the 13th to the 15th control plane when quantifying the results by considering the model structure error. Although it is not feasible to perform a straightforward evaluation of which way of HHR assessment (e.g., with or without the uncertainty quantification of the model structure) is more reasonable, we believe that more reliable contaminant concentration predictions will results in more reliable HHR assessments. Thus, the HHR assessed by quantifying the model structure uncertainty is more reliable than that obtained by neglecting the model structure error, and in this study, the difference between the two ways of HHR assessments was not negligible. Specifically, in the synthetic case study, if the HHR is assessed by neglecting the model structure error, this would result in underestimation of the risk and misleading depictions of the key risk zone, e.g., the 13th

Fig. 15. Probability densities of HHR at some control planes; the DNAPL concentrations were predicted without considering the model structure uncertainty. 13

Journal of Hydrology 584 (2020) 124690

Y. Pan, et al.

Fig. 16. Probability densities of HHR at some control planes; the DNAPL concentrations were predicted considering the model structure uncertainty.

without considering the model structure error may lead to misleading depictions of the key risk zone. It is more reliable to assess HHR by considering the model structure uncertainty than by neglecting the model structure error. (3) The GPR method is an effective technique to quantitatively describe the structure error of the DNAPL transport model.

and 14th control planes would be regarded as safe but these areas were actually at high risk to human health. According to the RAGS model (Eqs. (1) and (2)), HHR is affected by physiological, behavioral, and hydrogeological factors. The objective of this work was to investigate the effect of the model structure error on the HHR assessment for groundwater DNAPL contamination, and only hydrogeological factors were evaluated in this research. Currently, some studies have evaluated physiological and behavioral factors in HHR assessment, e.g., de Barros and Rubin (2008) found that at lower risk values and longer travel times, HHR is less sensitive to physiological parameters. Physiological and behavioral factors can be incorporated in the model parameters’ uncertainty in a future study. Moreover, the degradation of PCE was not considered in this study. The PCE would degrade into trichloroethylene (TCE), dichloroethylene (DCE), and vinyl chloride (VC), and finally, the contaminants would transform into non-toxic ethane in effect (Henri et al., 2015). All of these sub-products except ethane could pose HHR, and HHR assessment would be improved if the degradation of PCE was taken into account in groundwater models.

The proposed HHR assessment framework in this study is independent of any specific case studies, e.g., the GPR based model structure error quantification, the stochastic collocation SG surrogate method and thus, the proposed approach can be used for any complex water environmental models. The GPR method was used to learn the model structure error by means of data-driven statistical simulations. Certainly, other machine learning methods, e.g., support vector regression (Xu et al., 2014), artificial neural networks (Lu and Ricciuto, 2019), can also be used for this purpose. The squared exponential (SE) covariance function was applied to perform the GPR, and the learning performance of the GPR in regard to the structure error of the DNAPL transport model potentially could be improved by using other forms of covariance functions, e.g., the combination of SE and rational quadratic covariance functions (Williams and Rasmussen, 2006). In addition, the mass balance of the GPR based DNAPL transport modeling was not enforced during MCMC simulation in this study, and a recalibration strategy (Xu and Valocchi, 2015) would be an effective way to address this problem. This will be the focus of our future work.

6. Conclusion The HHR caused by DNAPL contamination in groundwater has aroused widespread public attention. We evaluated the impact of the model structure uncertainty on DNAPL transport modeling and an HHR assessment by using a data-driven error model to describe the model structure uncertainty. Based on two case studies, namely, a sandbox experiment and synthetic DNAPL transport model, the key findings are summarized as follows:

Declaration of Competing Interest 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.

(1) Ignoring the model structure error will lead to deviations in the model parameters’ posterior distributions, and the model parameter will be overfitted to compensate for the model structure error. The predictions of DNAPL transport model through quantifying the uncertainties of model parameters and structure simultaneously are more reliable than those obtained by calibrating only the model parameters. (2) The HHR assessment of groundwater DNAPL contamination

Acknowledgements This study was supported by the National Natural Science Foundation of China (41730856), and the Fundamental Research Funds for the Central Universities (0206-14380085). 14

Journal of Hydrology 584 (2020) 124690

Y. Pan, et al.

References

Hill, M., Tiedeman, C., 2007. Effective Calibration of Groundwater Models, With Analysis of Data, Sensitivities, Predictions, and Uncertainty, New York. Hoeting, J.A., Madigan, D., Raftery, A.E., Volinsky, C.T., 1999. Bayesian model averaging: a tutorial. Stat. Sci. 14 (4), 382–401. Højberg, A.L., Refsgaard, J.C., 2005. Model uncertainty – parameter uncertainty versus conceptual models. Water Sci. Technol. 52 (6), 177–186. https://doi.org/10.2166/ wst.2005.0166. Jakeman, J., Roberts, S., 2011. Local and Dimension Adaptive Sparse Grid Interpolation and Quadrature. Jakeman, J., Roberts, S., 2013. Local and Dimension Adaptive Stochastic Collocation for Uncertainty Quantification, 88, 181–203 pp. James, A.L., Oldenburg, C.M., 1997. Linear and Monte Carlo uncertainty analysis for subsurface contaminant transport simulation. Water Resour. Res. 33 (11), 2495–2508. Kennedy, M.C., O'Hagan, A., 2001. Bayesian calibration of computer models. J. Roy. Stat. Soc. B 63, 425–450. https://doi.org/10.1111/1467-9868.00294. Kueper, B.H., Redman, D., Starr, R.C., Reitsma, S., Mah, M., 1993. A field experiment to study the behavior of tetrachloroethylene below the water table: spatial distribution of residual and pooled DNAPL. Groundwater 31 (5), 756–766. https://doi.org/10. 1111/j.1745-6584.1993.tb00848.x. Kumar, B., Verma, V.K., Sharma, C.S., Akolkar, A.B., 2015. Estimation of toxicity equivalency and probabilistic health risk on lifetime daily intake of polycyclic aromatic hydrocarbons from urban residential soils. Hum. Ecol. Risk Assess. 21 (2), 434–444. https://doi.org/10.1080/10807039.2014.921530. Kumar, D., Singh, A., Jha, R.K., Sahoo, S.K., Jha, V., 2019. A variance decomposition approach for risk assessment of groundwater quality. Expo. Health 11 (2), 139–151. https://doi.org/10.1007/s12403-018-00293-6. Laloy, E., Vrugt, J.A., 2012. High-dimensional posterior exploration of hydrologic models using multiple-try DREAMzs and high-performance computing. Water Resour. Res. 48.Artn, W01526. https://doi.org/10.1029/2011wr010608. Lapworth, D.J., Baran, N., Stuart, M.E., Ward, R.S., 2012. Emerging organic contaminants in groundwater: a review of sources, fate and occurrence. Environ. Pollut. 163, 287–303. https://doi.org/10.1016/j.envpol.2011.12.034. Lehame, S.A., Lerner, D.N., Kueper, B.H., Smith, J.W.N., Wealthall, G.P., 2003. An Illustrated Handbook of DNAPL Transport and Fate in the Subsurface. Environment Agency, Almondsbury, Bristol. Lu, D., Ricciuto, D., 2019. Efficient surrogate modeling methods for large-scale Earth system models based on machine-learning techniques. Geosci. Model Dev. 12 (5), 1791–1807. https://doi.org/10.5194/gmd-12-1791-2019. Lu, D., Ye, M., Curtis, G.P., 2015. Maximum likelihood Bayesian model averaging and its predictive analysis for groundwater reactive transport models. J. Hydrol. 529, 1859–1873. Ma, X., Zabaras, N., 2009a. An adaptive hierarchical sparse grid collocation algorithm for the solution of stochastic differential equations. J. Comput. Phys. 228 (8), 3084–3113. https://doi.org/10.1016/j.jcp.2009.01.006. Ma, X., Zabaras, N., 2009b. An efficient Bayesian inference approach to inverse problems based on an adaptive sparse grid collocation method. Inverse Probl 25 (3). https:// doi.org/10.1088/0266-5611/25/3/035013. Artn035013. Maxwell, R.M., Kastenberg, W.E., 1999. Stochastic environmental risk analysis: an integrated methodology for predicting cancer risk from contaminated groundwater. Stoch. Environ. Res. Risk Assess. 13 (1–2), 27–47. https://doi.org/10.1007/ s004770050030. Maxwell, R.M., Kastenberg, W.E., Rubin, Y., 1999. A methodology to integrate site characterization information into groundwater-driven health risk assessment. Water Resour. Res. 35 (35), 2841–2856. https://doi.org/10.1029/1999WR900103. Mccray, J.E., Falta, R.W., 1997. Numerical simulation of air sparging for remediation of NAPL contamination. Groundwater 35 (1), 99–110. Moreno-Barbero, E., 2005. Simulation and Performance Assessment of Partitioning Tracer Tests in Heterogeneous Aquifers, 11, 395–404 pp. Mustafa, S.M.T., Nossent, J., Ghysels, G., Huysmans, M., 2018. Estimation and impact assessment of input and parameter uncertainty in predicting groundwater flow with a fully distributed model. Water Resour. Res. 54 (9), 6585–6608. https://doi.org/10. 1029/2017WR021857. Neriah, A.B., Paster, A., 2019. Enhancing groundwater remediation in air sparging by changing the pulse duration. Ground Water Monit. Remediat. 39 (1), 43–53. https:// doi.org/10.1111/gwmr.12308. Neuman, S.P., Wierenga, P.J., 2003. A Comprehensive Strategy of Hydrogeologic Modeling and Uncertainty Analysis for Nuclear Facilities and Sites (NUREG/CR6805) In: D.-. Division of Systems Analysis and Regulatory Effectiveness Office of Nuclear Regulatory Research U.S. Nuclear Regulatory Commission Washington (Editor), Washington, DC. Niessner, J., Helmig, R., 2007. Multi-scale modeling of three-phase–three-component processes in heterogeneous porous media. Adv. Water Resour. 30 (11), 2309–2325. https://doi.org/10.1016/j.advwatres.2007.05.008. Parker, J.C., Lenhard, R.J., 1987. A model for hysteretic constitutive relations governing multiphase flow: 1. Saturation-pressure relations. Water Resour. Res. 23 (12), 2187–2196. https://doi.org/10.1029/WR023i012p02187. Parker, J.C., Park, E., 2004. Modeling field-scale dense nonaqueous phase liquid dissolution kinetics in heterogeneous aquifers. Water Resour. Res. 40 (5). https://doi. org/10.1029/2003wr002807. Pflüger, D., Peherstorfer, B., Bungartz, H.-J., 2010. Spatially adaptive sparse grids for high-dimensional data-driven problems. J. Complex 26 (5), 508–522. https://doi. org/10.1016/j.jco.2010.04.001. Raftery, A.E., Gneiting, T., Balabdaoui, F., Polakowski, M., 2005. Using Bayesian model averaging to calibrate forecast ensembles. Mon. Weather Rev. 133 (5), 1155–1174. https://doi.org/10.1175/mwr2906.1.

Basu, N.B., Rao, P.S.C., Falta, R.W., Annable, M.D., Jawitz, J.W., Hatfield, K., 2008. Temporal evolution of DNAPL source and contaminant flux distribution: impacts of source mass depletion. J. Contam. Hydrol. 95 (3–4), 93–109. https://doi.org/10. 1016/j.jconhyd.2007.08.001. Bennett, D.H., James, A.L., McKone, T.E., Oldenburg, C.M., 1998. On uncertainty in remediation analysis: variance propagation from subsurface transport to exposure modeling. Reliab. Eng. Syst. Saf. 62 (1–2), 117–129. https://doi.org/10.1016/S09518320(97)00160-9. Beven, K., Binley, A., 1992. The future of distributed models: model calibration and uncertainty prediction. Hydrol. Process. 6 (3), 279–298. https://doi.org/10.1002/ hyp.3360060305. Bungartz, H.J., Dirnstorfer, S., 2003. Multivariate quadrature on adaptive sparse grids. Computing 71 (1), 89–114. https://doi.org/10.1007/s00607-003-0016-4. Burri, N.M., Weatherl, R., Moeck, C., Schirmer, M., 2019. A review of threats to groundwater quality in the anthropocene. Sci. Total Environ. 684, 136–154. https:// doi.org/10.1016/j.scitotenv.2019.05.236. Cao, T.T., Zeng, X.K., Wu, J.C., Wang, D., Sun, Y.Y., Zhu, X.B., Lin, J., Long, Y.Q., 2018. Integrating MT-DREAMzs and nested sampling algorithms to estimate marginal likelihood and comparison with several other methods. J. Hydrol. 563, 750–765. https://doi.org/10.1016/j.jhydrol.2018.06.055. Dafny, E., 2017. TCE longevity in the vadose zone and loading to the groundwater—The case of episodic NAPL releases from near-surface source. Environ. Technol. Innov. 7, 128–140. https://doi.org/10.1016/j.eti.2016.12.007. Dai, H., Chen, X., Ye, M., Song, X., Hammond, G., Hu, B., Zachara, J.M., 2019. Using Bayesian networks for sensitivity analysis of complex biogeochemical models. Water Resour. Res. 55 (4), 3541–3555. https://doi.org/10.1029/2018WR023589. Dai, H., Ye, M., 2015. Variance-based global sensitivity analysis for multiple scenarios and models with implementation using sparse grid collocation. J. Hydrol. 528, 286–300. https://doi.org/10.1016/j.jhydrol.2015.06.034. de Barros, F.P.J., Fiori, A., 2014. First-order based cumulative distribution function for solute concentration in heterogeneous aquifers: theoretical analysis and implications for human health risk assessment. Water Resour. Res. 50 (5), 4018–4037. https://doi. org/10.1002/2013WR015024. de Barros, F.P.J., Rubin, Y., 2008. A risk-driven approach for subsurface site characterization. Water Resour. Res. 44 (1), 568–569. https://doi.org/10.1029/ 2007WR006081. de Barros, F.P.J., Rubin, Y., Maxwell, R.M., 2009. The concept of comparative information yield curves and its application to risk-based site characterization. Water Resour. Res. 45 (45), 192–200. Doherty, J., Brebber, L. and Whyte, P., 2010. PEST: Model-independent parameter estimation user manual, Watemark Comput, Corinda, Australia. Enemark, T., Peeters, L.J.M., Mallants, D., Batelaan, O., 2019. Hydrogeological conceptual model building and testing: a review. J. Hydrol. 569, 310–329. https://doi. org/10.1016/j.jhydrol.2018.12.007. EPA, U.S., 2009. Risk Assessment Guidance for Superfund Volume I: Human Health Evaluation Manual (Part F, Supplemental Guidance for Inhalation Risk Assessment). In: U.S. EPA (Editor), Washington DC. Falta, R.W., Pruess, K., Finsterle, S., Battistelli, A., 1995. T2VOC user‘s guide. Gaganis, P., 2009. Model calibration/parameter estimation techniques and conceptual model error. In: Uncertainties in Environmental Modelling and Consequences for Policy Making. Springer Netherlands, Dordrecht, pp. 129–154. Gerstner, T., Griebel, M., 2003. Dimension-adaptive tensor–product quadrature. Computing 71 (1), 23. Giraud, Q., Gonçalvès, J., Paris, B., Joubert, A., Colombano, S., Cazaux, D., 2018. 3D numerical modelling of a pulsed pumping process of a large DNAPL pool: In situ pilotscale case study of hexachlorobutadiene in a keyed enclosure. J. Contam. Hydrol. 214, 24–38. https://doi.org/10.1016/j.jconhyd.2018.05.005. Gudbjerg, J., Heron, T., Sonnenborg, T.O., Jensen, K.H., 2005. Three dimensional numerical modeling of steam override observed at a fulllogcale remediation of an unconfined aquifer. Ground Water Monit. R 25 (3), 116–127. Guyton, K.Z., Hogan, K.A., Scott, C.S., Cooper, G.S., Bale, A.S., Kopylev, L., Barone, S., Makris, S.L., Glenn, B., Subramaniam, R.P., Gwinn, M.R., Dzubow, R.C., Chiu, W.A., 2014. Human health effects of tetrachloroethylene: key findings and scientific issues. Environ. Health Perspect. 122 (4), 325–334. https://doi.org/10.1289/ehp.1307359. Haario, H., Saksman, E., Tamminen, J., 2001. An adaptive Metropolis algorithm. Bernoulli 7 (2), 223–242. Hamraz, B.S., Akbarpour, A., Bilondi, M.P., Tabas, S.S., 2015. On the assessment of ground water parameter uncertainty over an arid aquifer. Arab. J. Geosci. 8 (12), 10759–10773. https://doi.org/10.1007/s12517-015-1935-z. Hassan, A.E., Bekhit, H.M., Chapman, J.B., 2008. Uncertainty assessment of a stochastic groundwater flow model using GLUE analysis. J. Hydrol. 362 (1–2), 89–109. https:// doi.org/10.1016/j.jhydrol.2008.08.017. Hassan, A.E., Bekhit, H.M., Chapman, J.B., 2009. Using Markov Chain Monte Carlo to quantify parameter uncertainty and its effect on predictions of a groundwater flow model. Environ. Model. Soft. 24 (6), 749–763. https://doi.org/10.1016/j.envsoft. 2008.11.002. Henri, C., Fernandezgarcia, D., De Barros, F., 2015. Probabilistic Human Health Risk Assessment of Chemical Mixtures: Hydro-Toxicological Interactions and Controlling Factors. Water Resour. Res. 4086–4108. https://doi.org/10.1002/2014WR016717. Henri, C.V., Fernàndezgarcia, D., De Barros, F.P.J., 2016. Assessing the joint impact of DNAPL source-zone behavior and degradation products on the probabilistic characterization of human health risk. Adv. Water Resour. 88, 124–138. https://doi.org/ 10.1016/j.advwatres.2015.12.012.

15

Journal of Hydrology 584 (2020) 124690

Y. Pan, et al. Refsgaard, J.C., Christensen, S., Sonnenborg, T.O., Seifert, D., Højberg, A.L., Troldborg, L., 2012. Review of strategies for handling geological uncertainty in groundwater flow and transport modeling. Adv. Water Resour. 36, 36–50. https://doi.org/10. 1016/j.advwatres.2011.04.006. Rodríguez-Lado, L., Sun, G., Berg, M., Zhang, Q., Xue, H., Zheng, Q., Johnson, C.A., 2013. Groundwater arsenic contamination throughout China. Science 341 (6148), 866–868. https://doi.org/10.1126/science.1237484. Saleh, H.N., Panahande, M., Yousefi, M., Asghari, F.B., Oliveri Conti, G., Talaee, E., Mohammadi, A.A., 2019. Carcinogenic and non-carcinogenic risk assessment of heavy metals in groundwater wells in neyshabur plain. Iran. Biol. Trace Elem. Res. 190 (1), 251–261. https://doi.org/10.1007/s12011-018-1516-6. Siirila, E.R., Maxwell, R.M., 2012. A new perspective on human health risk assessment: development of a time dependent methodology and the effect of varying exposure durations. Sci. Total Environ. 431 (5), 221–232. https://doi.org/10.1016/j.scitotenv. 2012.05.030. Smolyak, S.A., 1963. Quadrature and interpolation formulas for tensor products of certain classes of functions, Dokl. Akad. Nauk SSSR, pp. 123. Stone, H.L., 1970. Probability model for estimating three-phase relative permeability. J. Pet. Technol. 22 (02), 214–218. https://doi.org/10.2118/2116-PA. Tatti, F., Petrangeli Papini, M., Torretta, V., Mancini, G., Boni, M.R., Viotti, P., 2019. Experimental and numerical evaluation of Groundwater Circulation Wells as a remediation technology for persistent, low permeability contaminant source zones. J. Contam. Hydrol. 222, 89–100. https://doi.org/10.1016/j.jconhyd.2019.03.001. Teh, T., Nik Norulaini, N.A.R., Shahadat, M., Wong, Y., Mohd Omar, A.K., 2016. Risk assessment of metal contamination in soil and groundwater in Asia: a review of recent trends as well as existing environmental laws and regulations. Pedosphere 26 (4), 431–450. https://doi.org/10.1016/S1002-0160(15)60055-8. Teijon, G., Candela, L., Tamoh, K., Molina-Díaz, A., Fernández-Alba, A.R., 2010. Occurrence of emerging contaminants, priority substances (2008/105/CE) and heavy metals in treated wastewater and groundwater at Depurbaix facility (Barcelona, Spain). Sci. Total Environ. 408 (17), 3584–3595. https://doi.org/10.1016/j. scitotenv.2010.04.041. Troldborg, L., Refsgaard, J.C., Jensen, K.H., Engesgaard, P., 2007. The importance of alternative conceptual models for simulation of concentrations in a multi-aquifer system. Hydrogeol. J. 15 (5), 843–860. Troldborg, M., Nowak, W., Tuxen, N., Bjerg, P.L., Helmig, R., Binning, P.J., 2010. Uncertainty evaluation of mass discharge estimates from a contaminated site using a fully Bayesian framework. Water Resour. Res. 46 (12). https://doi.org/10.1029/ 2010wr009227.

Vrugt, J.A., ter Braak, C.J.F., Diks, C.G.H., Robinson, B.A., Hyman, J.M., Higdon, D., 2009. Accelerating Markov chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling. Int. J. Nonlin. Sci. Num. 10 (3), 273–290. Williams, C.K., Rasmussen, C.E., 2006. Gaussian Processes for Machine Learning. MIT Press Cambridge, MA. Wöhling, T., Bidwell, V.J., Barkle, G.F., 2012. Dual-tracer, non-equilibrium mixing cell modelling and uncertainty analysis for unsaturated bromide and chloride transport. J. Contam. Hydrol. https://doi.org/10.1016/j.jconhyd.2012.08.001. 140–141, 150–163. Xu, T., Valocchi, A.J., 2015. A Bayesian approach to improved calibration and prediction of groundwater models with structural error. Water Resour. Res. 51 (11), 9290–9311. https://doi.org/10.1002/2015wr017912. Xu, T., Valocchi, A.J., Choi, J., Amir, E., 2014. Use of machine learning methods to reduce predictive error of groundwater models. Groundwater 52 (3), 448–460. https://doi. org/10.1111/gwat.12061. Xu, T., Valocchi, A.J., Ye, M., Liang, F., 2017a. Quantifying model structural error: efficient Bayesian calibration of a regional groundwater flow model using surrogates and a data-driven error model. Water Resour. Res. 53 (5). Xu, T., Valocchi, A.J., Ye, M., Liang, F., Lin, Y.-F., 2017b. Bayesian calibration of groundwater models with input data uncertainty. Water Resour. Res. 53 (4), 3224–3245. https://doi.org/10.1002/2016wr019512. Zeng, L., Shi, L., Zhang, D., Wu, L., 2012. A sparse grid based Bayesian method for contaminant source identification. Adv. Water Resour. 37, 1–9. https://doi.org/10. 1016/j.advwatres.2011.09.011. Zeng, X.K., Wu, J.C., Wang, D., Zhu, X.B., 2016a. Assessing the pollution risk of a groundwater source field at western Laizhou Bay under seawater intrusion. Environ. Res. 148, 586–594. https://doi.org/10.1016/j.envres.2015.11.022. Zeng, X.K., Ye, M., Burkardt, J., Wu, J.C., Wang, D., Zhu, X.B., 2016b. Evaluating two sparse grid surrogates and two adaptation criteria for groundwater Bayesian uncertainty quantification. J. Hydrol. 535, 120–134. https://doi.org/10.1016/j.jhydrol. 2016.01.058. Zeng, X.K., Ye, M., Wu, J.C., Wang, D., Zhu, X.B., 2018. Improved nested sampling and surrogate-enabled comparison with other marginal likelihood estimators. Water Resour. Res. 54 (2), 797–826. https://doi.org/10.1002/2017wr020782. Zhang, G.N., Lu, D., Ye, M., Gunzburger, M., Webster, C., 2013. An adaptive sparse-grid high-order stochastic collocation method for Bayesian inference in groundwater reactive transport modeling. Water Resour. Res. 49 (10), 6871–6892. https://doi.org/ 10.1002/Wrcr.20467.

16