Sensitivity of advective travel time of contaminants to correlated formation porosities Jianting Zhu Abstract An approach is presented to quantify sensitivity of advective contaminant travel time to porosities of hydrogeologic units (HGUs) along the flowpaths when the porosities of different HGUs are correlated. The approach is an extension of the previous sensitivity analysis technique for independent input porosity cases. It is applicable in situations where a calibrated groundwater flow model exists, but a full contaminant transport model is not available. Three sensitivity indices are introduced based on the decomposition of covariance between the advective contaminant travel time and individual input porosities of HGUs. When the input HGU porosities are correlated, the three sensitivity indices quantify the total, intrinsic and correlated contributions from each individual HGU porosity, which should be considered in order to determine the relative importance to the uncertainty in advective contaminant travel time of the input HGU porosities contributing either independently or correlatedly. Two simple one-dimensional flow examples are presented to illustrate the applications of the approach to scenarios when each HGU has distinct porosity and situations of spatially variable porosity field. The approach is applicable to more complex multi-dimensional cases where advective contaminant travel time can be calculated based on simulated flow results from groundwater flow models. Keywords Contamination . Solute transport . Porosity correlation . Covariance decomposition . Sensitivity index

Introduction Porosity of hydrogeologic formation is an important parameter in estimating contaminant travel time and its Received: 10 April 2011 / Accepted: 18 September 2011 Published online: 11 October 2011 * Springer-Verlag 2011 J. Zhu ()) Desert Research Institute, Nevada System of Higher Education, 755 East Flamingo Road, Las Vegas, NV 89119, USA e-mail: [email protected] Tel.: +1-702-862-5416 Fax: +1-702-862-5427 Hydrogeology Journal (2012) 20: 135–141

uncertainty in groundwater contaminant migration scenarios (e.g., Zhu et al. 2006). The hydrogeologic formation in a region is often categorized as a composition of multiple hydrogeologic units (HGUs). An HGU is defined to have reasonably distinct hydrologic properties because of its geological and structural characteristics. The porosities of HGUs are particularly significant in applications when estimations of contaminant travel times are needed, but groundwater contaminant transport models are not available and only groundwater flow models are available. Groundwater flow models provide flux results which can be used in conjunction with available porosity information to determine groundwater seepage velocity. The seepage velocity then allows advective contaminant travel times to be estimated. The advective contaminant travel times are often required in many applications such as designing monitoring networks and remediation strategies of contaminant sites and determining the associated risks of groundwater contamination. Large uncertainty in many important parameters, including the formation porosity, may exist due to lack of knowledge and insufficient field measurements. Therefore, the HGU porosity value within any HGU is usually considered to be a random variable following a probability density function. As a result, the contaminant travel time is also a random variable. Sensitivity analysis is concerned with understanding how the input variations impact the changes of the output. In sensitivity analysis, the inputs are formally treated as random variables with specified distributions, and consequently, the output is also a random variable with a probability distribution. The characterization of the output distribution, given the input probability distribution, is the goal of uncertainty analysis. The assessment of the relative importance of the inputs in the output is the objective of sensitivity analysis. While the input parameters often exhibit correlations in typical field conditions, and can significantly affect the predictions and the associated uncertainties of flow and contaminant transport (Pohlmann et al. 2002; Lemke et al. 2004; Manache and Melching 2008), existing sensitivity analysis methods of flow and transport typically adopt the assumption of independent parameters (e.g., Boateng and Cawlfield 1999; Boateng 2007; Zhu et al. 2010). Understanding the contribution of each parameter and joint contribuDOI 10.1007/s10040-011-0788-0

136

tions of correlated parameters to output uncertainties is critical to uncertainty reduction. Therefore, it is desirable to incorporate the parameter correlations into the sensitivity analysis and to investigate the effects of parameter correlations on the results of uncertainty and sensitivity analysis (Fang et al. 2004; Jacques et al. 2006). Zhu et al. (2010) introduced an approach to quantify uncertainty and sensitivity of advective travel time to the porosities of HGUs along groundwater flow paths, assuming that the input porosities are independent. This study seeks to extend the approach in Zhu et al. (2010) to quantify the sensitivity of advective contaminant travel time to the porosities of HGUs when the porosities of different HGUs are correlated. The focus is on the impact of porosity correlations on the sensitivity of advective contaminant transport time. Specifically, the importance of an HGU porosity is quantitatively assessed through the covariance between the overall advective travel time of contaminants and the portion that is contributed from an individual HGU porosity. The covariance is then decomposed into two sensitivity indices, respectively reflecting the contribution of the HGU porosity due to its own variance, which is referred to as the intrinsic sensitivity index of the HGU porosity, and the contribution due to the correlation between the HGU porosity and other HGU porosities, which is called the correlated sensitivity index of the HGU porosity.

various HGUs are correlated. From Eq. (1), the mean advective travel time can be expressed as, N X

Y ¼

aj ’j

ð2Þ

j¼1

The variance of the advective travel time is, N  2 X 2Y ¼ Y 2  Y ¼ a2j 2’j þ 2 j¼1



N X N X

aj ak jk ’j ’k

ð3Þ

j¼1 k¼jþ1

where σφj is the standard deviation of φj, ρik is the correlation coefficient of φi and φk. In general, the total contribution from HGU i porosity φi can be quantitatively assessed through the covariance between the output (i.e., the advective travel time) and the portion of that travel time that is contributed from the porosity of HGU i (i.e., aiφi) as follows,

Covðai ’i ; Y Þ ¼ Cov ai ’i ;

N X

! aj ’j

j¼1

¼ ai

Sensitivity indices of porosities

N X

aj rij s ’i s ’j

ð4Þ

j¼1

The advective contaminant travel time is related to the seepage velocity of groundwater flows. The seepage velocity along the flowpath X, Va(X), is related to the Darcy flux Vd(X) via porosity (φ) by Va(X) = Vd(X)/φ. The advective travel time through a small segment ΔX within an HGU j along the flowpath is then equal to [ΔX/ Vd(X)]φj. Therefore, the contaminant travel time within any HGU is linearly related to the porosity of this HGU. The total travel time of the contaminant is then simply a linear combination of porosities of all HGUs that it has travelled through,

Y ¼

N X

where Cov denotes the covariance operator. Based on this concept, three sensitivity indices can then be introduced for the HGU porosity φi as follows,

Si ¼

ai ¼

aj ’j

Covðai ’i ; Y Þ s 2Y

ð1Þ

N P j¼1

N P j¼1

a2j s 2’j

þ2

aj rij s ’i s ’j

N P N P j¼1 k¼jþ1

ð5Þ aj ak rjk s ’j s ’k

j¼1

where Y is the advective travel time of contaminant; φj is the porosity of HGU j; and N is the number of HGUs; aj (j = 1, …, N) are the coefficients which are determined based on the Darcy-flux results from the groundwater flow model. The coefficient aj for HGU j can be calculated from the travel time by setting φj to 1 for HGU j and to a very small number (e.g. 10–10) for the other HGUs (Zhu et al. 2010). In this study, the sensitivity analysis approach in Zhu et al. (2010) is extended to the case when input porosities of Hydrogeology Journal (2012) 20: 135–141

ð1Þ

Si

¼ ¼

Varðai ’i Þ s 2Y N P j¼1

a2j s 2’j þ 2

a2i s 2’i N P N P j¼1 k¼jþ1

ð6Þ aj ak rjk s ’j s ’k

DOI 10.1007/s10040-011-0788-0

137 ð2Þ

Si

ð1Þ

¼ Si  Si

ai ¼

N P j¼1

a2j s 2’j

N P j¼1:j6¼i

þ2

aj rij s ’i s ’j

N P N P j¼1 k¼jþ1

ð7Þ aj ak rjk s ’j s ’k

In the preceding expressions, Var denotes variance operator. Si represents the portion of variance of advective travel time Y that comes from the HGU i porosity quantitatively through the term aiφi, which is called the total sensitivity index of φi. Sð1Þ i is the contribution of φi due to the variance of the term aiφi, which is referred to as the intrinsic sensitivity index of φi. The remaining part of the total sensitivity index Si, i.e., Sð2Þ i as expressed in Eq. (7), is due to the correlation between HGU i porosity φi and other HGU porosities, which is called the correlated sensitivity index of φi. More specifically, the covariance Cov(aiφi, Y) is the total contribution of aiφi, consisting of its intrinsic portion a2i s 2ai reflecting aiφi’s contribution to the variations of the output advective travel time of contaminant Y, and its correlated N P portion ai j¼1:j6¼i ajij ’ ’ reflecting the influence of the interaction between aiφi and other HGUs through the correlated random porosities. In the case when all HGU porosities are independent (i.e., ρij = 0 for i ≠ j), Sð2Þ i is equal to zero, and Sð1Þ in Eq. (6) reduces to the square of the i standardized coefficients presented in Zhu et al. (2010), which apportioned the advective travel time variances to that of the individual input HGU i porosity. i

j

Illustrative examples Two examples are presented to demonstrate the applications of the sensitivity analysis approach using a simple one-dimensional steady-state flow domain shown in Fig. 1. For this simple configuration, the coefficients in Eq. (1) can be easily calculated as aj = Lj/q = fjL/q, where Li is the length of HGU j, q is the Darcy flux in the onedimensional domain, fj (= Li/L) is the length fraction of HGU j. In this case, since aj is simply proportional to fj,

the variance of the advective travel time is generally quadratic functions of the length fractions and individual HGU porosity variances, and a linear function of correlation coefficients among individual HGU porosities, as seen from Eq. (3). As q and L cancel out in Eqs. (5)– (7), it can be concluded that the sensitivity indices are only related to the length fraction fj and the statistics of porosities, including standard deviation and correlation coefficients. In other words, the coefficients aj in Eqs. (5)– (7) can be simply replaced by the length fraction fj for this special domain configuration. In the first example, the one-dimensional domain has only four separate HGUs, which is used to illustrate the applicability of the proposed approach to scenarios when each HGU has distinct porosity. Table 1 list parameters used in the example including values of saturated hydraulic conductivities, HGU length fractions and porosity statistics, which represent the base case of this example. For the base case, porosities for all HGUs are assumed to be equally correlated with the correlation coefficient of 0.5. Note that the saturated hydraulic conductivities are required to calculate Darcy flux. In this simple illustrative example, it is hoped to demonstrate the influence of various variables on the sensitivity results by deviating them from the base case as shown in Table 1. For this simple example, the Darcy flux and the advective travel time, which are the essential quantities in the approach, can be determined analytically. The steady-state Darcy flux is equal to the product of prescribed hydraulic gradient and the HGU-length weighted harmonic mean of hydraulic conductivities of all individual HGUs. For more complicated field applications, the sensitivity analysis procedure will be the same, but the Darcy flux and the corresponding advective travel time may need to be determined using more sophisticated groundwater flow models and particle tracking techniques. For field applications, the hydraulic conductivities of HGUs are likely to be either measured directly or derived from hydraulic head measurements using groundwater flow models. Figure 2 shows the sensitivity indices for all the four HGUs for the base case. It is apparent that HGU 1 dominates travel-time uncertainty as it has the largest standard deviation (or variance) of its porosity. HGU 1 has not only the largest total sensitivity index, but also the

Fig. 1 One-dimensional example of flow domain and hydrogeologic units. For steady-state flow, Darcy flux q is the same for all hydrogeologic units. Each hydrogeologic unit could have its own mean porosity and variance. The porosities of hydrogeologic units can be correlated

Hydrogeology Journal (2012) 20: 135–141

DOI 10.1007/s10040-011-0788-0

138 Table 1 Saturated hydraulic conductivity, length fraction and porosity statisticsa for the base case of one-dimensional example HGU No.

Ks (m/day)

Length fraction

Mean porosity

SD

1 2 3 4

12 1.2 0.12 0.012

0.25 0.25 0.25 0.25

0.4 0.2 0.1 0.05

0.4 0.2 0.1 0.05

a

In addition to mean and standard deviation (SD) given in the table, the correlation coefficients are all 0.5 for the base case

largest portion of the total index from the intrinsic sensitivity index. For other HGUs, the total sensitivity indices are mainly from the correlated portion since their intrinsic sensitivity indices are small. It is shown that the porosity of HGU 1 is the most dominant contributor to the advective travel-time uncertainty. Figure 3 shows the results of sensitivity indices as functions of correlation coefficients between HGU 1 porosity (φ1) and HGU 2 porosity (φ2) when the correlation coefficients among other HGU porosities remain the same as the base case at 0.5. As expected, as the correlation between φ1 and φ2 increases, their respective intrinsic indices decrease while correlated indices for both HGUs increases. However, it is interesting to note that the net results indicate the total sensitivity index of φ1 decreases while that of φ2 increases. The results imply that the less important HGU benefits from its correlation with the more dominant HGU, while the uncertainty contribution of more dominant HGU decreases by the correlation with the less important HGU. As a result, these two total indices move close to each other as their correlation increases. For the more dominant HGU porosity, the increase of its correlated index is slower than the decrease of its intrinsic index. The net result is the decrease of its total sensitivity index for the more dominant HGU. In contrast, the opposite is true for the less important HGU and the net result is the increase of total sensitivity index for the less dominant HGU. Figure 4 plots the results of the sensitivity indices as functions of coefficient of variation (CV) of HGU 1 porosity when the standard deviations of other HGU

porosities remain the same as the base case shown in Table 1. As expected, the total index and intrinsic index of φ1 increase with increasing the CV of φ1. However, it is interesting to note that the extents of their respective increases are different. The intrinsic index increases largely as a quadratic function of standard deviation of φ1, as shown in Eq. (6), which implies that the increase of the intrinsic index is slow when the standard deviation of φ1 is small. The

(a)

(b)

0.6

(c)

Sensitivity Index

0.5 Total

0.4

Intrinsic Correlated

0.3 0.2 0.1 0.0

1

2

3

4

HGU No. Fig. 2 Total, intrinsic and correlated sensitivity indices for each hydrogeologic unit (HGU) for the base case Hydrogeology Journal (2012) 20: 135–141

Fig. 3 Sensitivity indices as functions of correlation coefficients between hydrogeologic unit (HGU) 1 porosity and HGU 2 porosity: a total sensitivity index, b intrinsic sensitivity index, and c correlated sensitivity index DOI 10.1007/s10040-011-0788-0

139

(a)

(a)

(b)

(b)

(c)

Fig. 4 Sensitivity indices as functions of coefficient of variation (CV) of hydrogeologic unit (HGU) 1 porosity: a total sensitivity index, b intrinsic sensitivity index, and c correlated sensitivity index

increase of total sensitivity index is mainly from the correlated index when the CV of φ1 is small. Upon further increase in the CV of φ1, the extent of increase in the intrinsic sensitivity index of φ1 increases rapidly due to its largely quadratic function of the standard deviation of φ1. As a result, the correlated sensitivity index shows a non-monotonic behavior as seen in Fig. 4c. Figure 5 shows the results of the sensitivity indices as functions of length fraction of HGU 1. Note that the results shown in Fig. 5 are produced by increasing the length fraction of HGU 1 and decreasing the fractions of the other three HGUs equally. Similar to the impact of standard deviation (or CV) of φ1 as discussed earlier, the interactions of intrinsic and correlated sensitivity indices determine how the total sensitivity index of φ1 increases as HGU 1 becomes a more and more significant portion in the domain and eventually completely dominant as all uncertainty of the advective contaminant travel time is from the variance of HGU 1 porosity φ1. Hydrogeology Journal (2012) 20: 135–141

(c)

Fig. 5 Sensitivity indices as functions of length fraction of hydrogeologic unit (HGU) 1: a total sensitivity index, b intrinsic sensitivity index, and c correlated sensitivity index

In the second example, the one-dimensional domain is discretized into 20 HGUs. Note that the discretization can be done into any numbers for spatially variable porosity fields. This example is presented to illustrate the applicability of the proposed approach to situations of spatially variable porosity field. In this example, the domain is assumed to be stationary with spatially variable porosity. The correlation coefficient is generated based on the exponential correlation function (Gutjahr and Gelhar 1981; Zhu and Satish 2001),   rij ¼ exp xij =l

ð8Þ

where ξij is the separation distance between HGU i and HGU j, 1 is a parameter representing a correlation range. Note that the mean and variance do not affect sensitivity indices of the travel time; they only influence mean and variance of the advective contaminant travel time. DOI 10.1007/s10040-011-0788-0

140

The variations of the sensitivity indices in relation to the locations in the domain are shown in Fig. 6. The sensitivity indices are symmetric about the middle point of the domain and the total sensitivity index reaches the largest in the middle of the domain (Fig. 6a). The variations of the total sensitivity index with locations are due to the correlated portion of the sensitivity index, as the intrinsic sensitivity index is constant throughout the domain shown in Fig. 6b. Note that the intrinsic sensitivity index is only related to the variance of porosity which is constant throughout the domain for a stationary stochastic porosity field. The middle point is the most correlated location in the domain as it correlates in both directions, which results in the highest correlated sensitivity index as shown in Fig. 6c. The variations of the

(a)

(b)

sensitivity indices with the correlation range are illustrated in Fig. 7 at a few selected locations. Since the intrinsic sensitivity index is constant throughout the domain, the three red curves of intrinsic sensitivity index results in Fig. 7 are identical for the three locations. As expected, the intrinsic sensitivity index decreases rapidly as the correlation range increases. When the correlation range decreases, the stochastic porosity field approaches white noise and both the total sensitivity index and the intrinsic sensitivity index reach 0.05 (i.e., = 1/N where N is the total number of HGUs in the domain), while the correlated sensitivity index approaches zero at all locations. Two simple examples have been presented for illustrative purposes—the first example involves porous formation consisting of a few segments with distinct porosities, while the second example illustrates how the approach can be used in spatially variable continuous field where spatial correlation structure of porosity is defined by a correlation function. It should be pointed out that the approach is applicable in any situations where spatial correlation and variance of formation porosity are prescribed no matter whether the porosity field is stationary or non-stationary. Since there are no restrictions on the spatial structure of porosity in applying the approach, the approach can also be similarly applied for more general conditions where the domain may be divided into a few large HGUs and within each HGU the porosity is spatially variable and correlated.

Concluding remarks The uncertainty of hydrogeologic units’ (HGUs) porosity is an important factor on advective contaminant travel time and its uncertainty. This study presented an approach to apportion the total sensitivity index of individual HGU porosities into the intrinsic sensitivity index and the correlated sensitivity index when HGU porosities are correlated. The approach is applicable for any porous media as long as the variances and correlation coefficient

(c)

Fig. 6 Sensitivity indices as functions of dimensionless distance (x/L): a total sensitivity index, b intrinsic sensitivity index, and c correlated sensitivity index Hydrogeology Journal (2012) 20: 135–141

Fig. 7 Sensitivity indices as functions of dimensionless correlation range (1 /L) at three selected dimensionless distances of 0.025, 0.275 and 0.525. These three dimensionless distances are shown inside parenthesis in the legends. The three red curves of intrinsic sensitivity index results are identical for the three locations DOI 10.1007/s10040-011-0788-0

141

matrix of formation porosity can be prescribed. The values of all these three sensitivity indices need to be considered in order to quantitatively determine the relative importance of the input HGU porosities either independently or correlatedly. For the case of independent input HGU porosities, these indices reduce to a single sensitivity index presented in Zhu et al. (2010). Two simple examples have been presented to illustrate the application of the proposed approach and the results demonstrate that the partition of the advective contaminant travel time sensitivity is a combined outcome of interactions among multiple factors even for the simple illustrative examples shown in this study. The proposed approach is able to quantitatively capture these interactions and apportion the overall sensitivity of contaminant travel time to the individual HGU porosities and their possible correlations. Acknowledgements Funding support from the Maki Chair program at Desert Research Institute, Nevada, USA is greatly acknowledged.

References Boateng S (2007) Probabilistic unsaturated flow along the textural interface in three capillary barrier models. J Environ Eng 133 (11):1024–1031 Boateng S, Cawlfield JD (1999) Two-dimensional sensitivity analysis of contaminant transport in the unsaturated zone. Ground Water 37(2):185–193

Hydrogeology Journal (2012) 20: 135–141

Fang S, Gertner GZ, Anderson AA (2004) Estimation of sensitivity coefficients of nonlinear model input parameters which have a multinormal distribution. Comput Phys Commun 157:9–16 Gutjahr AL, Gelhar LW (1981) Stochastic models of subsurface flow: infinite versus finite domains and stationarity. Water Resour Res 17(2):337–350 Jacques J, Lavergne C, Devictor N (2006) Sensitivity analysis in presence of model uncertainty and correlated inputs. Reliab Eng Syst Safe 91:1126–1134 Lemke LD, Abriola LM, Lang JR (2004) Influence of hydraulic property correlation on predicted dense nonaquesous phase liquid source zone architecture, mass recovery and contaminant flux. Water Resour Res 40:W12417. doi:10.1029/ 2004WR003061 Manache G, Melching CS (2008) Identification of reliable regression- and correlation-based sensitivity measures for importance ranking of water-quality model parameters. Environ Modell Softw 23:549–562 Pohlmann KF, Hassan AE, Chapman JB (2002) Modeling densitydriven flow and radionuclide transport at an underground nuclear test: uncertainty analysis and effect of parameter correlation. Water Resour Res 38(5):1059. doi:10.1029/ 2001WR001047 Zhu J, Satish MG (2001) Stochastic analysis of groundwater flow in aquifers with leakage. Stoch Environ Res Risk Assess 15:228– 243 Zhu J, Pohlmann KF, Ye M, Shafer DS, Carroll RWH (2006) Groundwater pathways from regions upgradient of potential Yucca Mountain repository. International High-Level Radioactive Waste Management Conference, 30 April–4 May 2006, Las Vegas, NV, pp 623–629 Zhu J, Pohlmann KF, Chapman JB, Russell CE, Carroll RWH, Shafer DS (2010) Sensitivity of solute advective travel time to porosities of hydrogeologic units. Ground Water 48(3):442– 447. doi:10.1111/j.1745-6584.2009.00664.x

DOI 10.1007/s10040-011-0788-0