Stochastic analysis of the hydraulic conductivity estimated for a heterogeneous aquifer via numerical modelling

The paper aims to evaluate the impacts of the average hydraulic conductivity of the heterogeneous aquifer on the estimated hydraulic conductivity using the observations from pumping tests. The results of aquifer tests conducted at a karst aquifer are first introduced. A MODFLOW groundwater flow model was developed to perform numerical pumping tests, and the heterogeneous hydraulic conductivity (K) field was generated using the Monte Carlo method. The K was estimated by the Theis solution for an unconfined aquifer. The effective hydraulic conductivity (Ke) was calculated to represent the hydraulic conductivity of a heterogeneous aquifer. The results of numerical simulations demonstrate that Ke increase with the mean of hydraulic conductivity (EK), and decrease with the coefficient of variation of the hydraulic conductivity (Cv). The impact of spatial variability of K on the estimated Ke at two observation wells with smaller EK is less significant compared to the cases with larger EK.


INTRODUCTION
Study of the heterogeneity aquifers is important to the analysis of groundwater flow and contaminant transport.The heterogeneity of an aquifer was shown as the spatial variation of hydraulic conductivity (K).Sudicky (1986) reported that K is a key factor in controlling groundwater flow and solute transport and it can vary significantly over short distances.Moench et al. (2001) noted that the differences (errors) between simulated and measured drawdowns in groundwater models are caused by local variations of hydraulic properties, primarily in K.In particular, karst aquifers are more spatially varied compared with porous medium; the K values of a karst aquifer are of higher variability than that of a porous medium (Illman, 2006).Sanchez-Vila et al. (2006) summarized the evolution of the concept of representative hydraulic conductivity in the literature.The need for relationships between the ensemble values of interpreted parameters and the actual statistics of the spatial random function was emphasized by Sanchez-Vila et al. (2006).Moench et al. (2001) first quantified the impacts of the heterogeneity of an aquifer on the estimates of hydraulic conductivity, and indicated that the effects of heterogeneity for a mildly heterogeneous aquifer on the results of the aquifer tests were small.To date, the differences between the estimated K and the real ones are still unknown for aquifers with high heterogeneity.In this study, therefore, the main focus is on assessing the impacts of the aquifer heterogeneity on the estimated K, through taking into account different degrees of heterogeneity.
In general, the understanding of heterogeneity is based on numerous field investigations, such as pumping tests, slug tests, ground tomography and other geophysical work; however, it is impossible to perform field investigations at each test site in a small scale to analyse the hydraulic properties of an aquifer.Numerical simulation is more convenient to conduct than other field methods (Zhang et al., 2007).
This paper aims to analyse the uncertainty of the estimated K values in a heterogeneous karst aquifer.First, the K values based on the field tests conducted in a karst aquifer are introduced and used for the following groundwater numerical simulation.Then, numerous hypothetic numerical tests were conducted to study the uncertainty of estimated K in a heterogeneous aquifer.In the numerical simulation, the hydraulic heads of an observation well were recorded during a pumping test to estimate a K value by an optimization algorithm.Finally, numerous numerical tests with different mean and variation of K values were performed to evaluate the impacts of the simulated hydraulic conductivity fields on the estimates of K.

Estimates of K from in situ pumping tests
Three pumping-recovery tests were conducted in the unconfined karst aquifer within the Houzhai basin from 13 to 21 August 2008.The Houzhai basin is near the town of Puding, southwest China.Drawdowns from the pumping and recovery periods were recorded every 10 minutes at the beginning of the tests.At the end of the pumping-recovery tests, the changes of drawdown were minimal so that the groundwater flow in the aquifer was regarded as steady-state.The Dupuit solution for an unconfined aquifer was used for steady-state.Transient drawdowns collected at each observation well in the pumping-recovery tests were used to estimate the hydraulic parameters using the Moench solution (Lu et al., 2011).In total, 13 K estimates were obtained based on the pumping-recovery tests.The probability distribution of these K estimates (Fig. 1) indicates the pumping K is normally distributed.

Generation of the random K fields in numerical simulation
Generally, the K estimated through a pumping test is affected by the mean and the variation of K, and the observation well locations.In numerical analyses, the mean of the K for each grid, labelled as EK, represents the average K of the aquifer and the coefficient of variation, Cv, represents the variability of the K value within the model domain.The K field was generated randomly according to the specific EK and Cv.In a karst area, the spatial correlation exists slightly because of the strong heterogeneity of aquifer.Therefore, the Monte Carlo simulation without spatial correlation is applicable for the karst aquifer in the study.

Forward groundwater flow modelling via MODFLOW
MODFLOW was used to simulate the hydraulic heads in response to the hypothetical groundwater pumping.The simulated aquifer was unconfined, and the dimension of the model was 2500 m in length, 2500 m in width and 50 m in depth.The model area was discretized into 100 columns and 100 rows in the x and y directions (see Fig. 2), respectively.Since the determination of the vertical K was not the focus of this study, the model had only one layer of 50 m depth, thereby inducing a two-dimensional groundwater flow system.A pumping well was placed at the centre of the model.The pumping test lasted for 24 hours with a constant pumping rate of 1728 m 3 /d, which is the actual rate of a pumping test conducted at Houzhai (Lu et al., 2013).Two observation wells, labelled Wmin and Wmax (see Fig. 2), were placed at different directions from the pumping well.Well Wmin was located in the short axis of the influence since pumping, and well Wmax was located in the long axis of the influence.

Estimation of the grid hydraulic conductivity and the effective hydraulic conductivity
Theis solution for the unconfined aquifer (Quasi-Theis method, Chen, 2001), was used in the inverse problem to calculate the hydraulic heads of the simulated aquifer, which is expressed as: where h0 is the hydraulic head prior to pumping, h(t) is the hydraulic head at time t after pumping, and W(u) is the well function.The gradient method reported by Chen and Ayers (1997) was regarded as the inverse method to estimate the K of the hypothetical aquifers, which was to minimize the square errors between the observed hydraulic heads from the forward modelling and the hydraulic heads using the Quasi-Theis solution.
The effective hydraulic conductivity, Ke, was generally regarded as the hydraulic conductivity of the heterogeneous aquifer (Renard and de Marsily, 1997), and Ke can be calculated from the discrete K values obtained from the previous parameter estimations.The equations proposed by Gutjahr et al. (1978) were adopted in this study because they are simple to conduct without considering the boundary conditions.Ke was calculated as follow, respectively,

Ke Kg σ
(2) where 2 Y σ is the variance of lnK, Kg is the geometric mean of K values, n is the dimensionality of media, N is the number of the K values.In particular, since the numerical modelling was twodimensional in the study, Ke is equal to Kg. EK and Cv for the K of model grids represent the average K and the variability of the hypothetical aquifer, respectively.In the numerical simulations, the EK values were set to be 10, 50 and 100 m/d, respectively, and the Cv values of given K at each mean value were set to be 0.1 (small variance), 0.5 (medium variance) and 1 (large variance).The K values of the grids in the model were assigned according to ten thousand K values from the random number generator at a time.Once the EK and Cv values were fixed, one hundred realizations were run to obtain an estimate of Ke for an observation well.The impacts of the aquifer heterogeneity were assessed by generation of the heterogeneous K fields, running the numerical pumping tests in these heterogeneous aquifers and performing inverse modelling for K estimation.Chen and Ayers (1997) noted that the gradient method has a good efficiency on computation and accuracy of the estimates of parameter.In their study, the gradient method was conducted using the datasets collected from the pumping test in field.The simulated hydraulic field was significantly affected by the setting of the numerical model.Therefore, the estimations of K in the homogenous aquifers were verified first.

Comparison of the given K with the estimated K for a homogenous aquifer
The parameter estimation based on the Quasi-Theis solution was initially conducted for homogeneous aquifers with K values ranged from 10 to 700 m/d.The estimates of K were compared with given values.To evaluate the impact of specific yield (Sy), numerous pumping tests were simulated while the Sy values were equal to 0.1 and 0.01.For the homogeneous aquifers, one observation well was located in the model to record the hydraulic heads for the parameter estimation (Fig. 3).
With Sy set to 0.01 and 0.1, respectively, the estimated K and given K values are plotted in Fig. 3, which shows a good match between the two groups of K values in the homogeneous aquifer.The estimated errors of K values are less than ±5% related to the given values (Fig. 3).

Evaluation of the inverse method in a heterogeneous aquifer
In addition, the accuracy of results using the inverse method in a heterogeneous aquifer was verified.Beck and Arnold (1977) reported a method which can give the confidence interval (CI) of the estimated parameter.For each set of the EK and Cv value, 100 numerical simulations were performed to evaluate the accuracy of the estimated K values.The estimated standard error (ESE) and CI of the estimated K represent the possible deviations (Beck and Arnold, 1977).Due to the large variability of the estimated K values, the length of confidence interval (LCI) is a better parameter to evaluate the accuracy of estimated K values.
Table 1 lists the results of ESE and LCI.The estimated K values differ significantly between each other.While the EK is fixed, ESE and LCI increase with the Cv value.The maximum ESE is less than 0.3 m/d and the LCI values for the estimated K are mostly less than 1 m/d.Therefore, in terms of solution algorithm, the estimated K of the heterogeneity aquifer is reliable in view of the small ESE and LCI values.

Impacts of the level of aquifer heterogeneity on the inversely estimated K
When EK = 10, 50 and 100 m/d and Cv = 0.1, 0.5 and 1, respectively, the Ke variation with EK for the two observation wells are shown in Fig. 4. When Cv is 0.1, the Ke values are the closer to EK rather than other two Cv values.In terms of the well Wmax, the Ke values are 9.94, 47.03 and 82.82 m/d, respectively, corresponding to Cv = 0.1 and EK = 10, 50 and 100 m/d, respectively.In the same cases, the Ke values are 9.85, 37.65 and 92.28 m/d for the well Wmin.The larger Cv value represents the larger variability of the given hydraulic conductivity field.Because a smaller Cv indicates that the field of K has a small discrepancy compared to a homogenous field, the estimated Ke values are closer to the EK value.Moench et al. (2001) also reported that the calculated hydraulic conductivity values from the aquifer tests were very similar to the geometric mean values (equal to Ke, here) calculated from many small-scale measurements made in the same aquifer.It should be noted that the variability of the estimates at small scale for that aquifer was relatively small (Cv of ln K = 0.24).Therefore, this affirmed that the estimated Ke values are close to the average hydraulic conductivity value for the quasi homogeneous aquifer value with a small Cv.When EK is assigned to be 10 m/d constantly, the difference of the estimations of Ke between Cv = 1 and Cv = 0.1 is 4.81 m/d for well Wmin, and the corresponding difference for well Wmax is 2.76 m/d.Therefore, the relative differences to the EK value (10 m/d) are 48.1% for well Wmin, and 27.6% for well Wmax.However, as the EK value changes to 100 m/d, the relative difference of Ke between Cv = 1 and Cv = 0.1 rises to 65.97% and 49.79% for well Wmin and well Wmax, respectively.The estimated Ke values have low difference between any two of the three Cv values while the EK value is small (see Fig. 4).When the EK value is small, the impact of variability of K on the estimated K is less significant compared to the cases with larger EK values.These results indicate that Ke is affected by the spatial variability of K at each grid and the influence of the variability will be stronger as the average K increases.
Comparing Fig. 4(a) and (b), the estimated Ke values at well Wmin are different to the estimated ones at well Wmax.That is because the local heterogeneity of the aquifer around the two observation wells differs.The differences between the estimations at well Wmin and well Wmax indicate that the groundwater flows at different locations can behave differently.

SUMMARY AND CONCLUSIONS
With the increasing variation of the given K field, Ke departs far from the given mean values.When EK remains constant, Ke drops as Cv increases.Moreover if the aquifer has a larger EK value, Ke has a larger decrease with the same Cv value.If the hydraulic conductivity of the karst aquifer increases and the heterogeneity of the aquifer is fixed, the estimated hydraulic conductivity will increase with the average hydraulic conductivity of the aquifer.However, the increasing heterogeneity of the aquifer leads to the decrease of estimated hydraulic conductivity.The effective hydraulic conductivity (Ke) will be affected significantly by the aquifer heterogeneity.Therefore, the estimated Ke has big uncertainty while the aquifer has strong heterogeneity.
Consequently, the comprehensive characterization via some aquifer tests at specific locations cannot be regarded as the effective hydraulic conductivity of the highly heterogeneous aquifer.A large number of field investigations and exploration of the aquifer need be conducted to obtain the accurate spatial patterns of the field heterogeneity.

Fig. 1
Fig.1The fitted curve of K estimates from pumping tests with a normal cumulative probability distribution.

Fig. 2
Fig.2The planar map of the model domain with a random K field. 2

Fig. 3
Fig. 3 Verification of the determination of hydraulic conductivities with Sy = 0.01 and Sy = 0.1 in a homogenous aquifer.Top and bottom dash lines represent the 1.05 times and the 0.95 times the given value, respectively.

Fig. 4
Fig.4The effective hydraulic conductivity varied by the mean of given hydraulic conductivity at well Wmin and well Wmax with respect to Cv = 0.1, 0.5 and 1, respectively.

Figure 4
Figure 4 shows that the estimated Ke values at well Wmin and well Wmax increase with EK, but decrease with Cv.For instance, the Ke value derived for well Wmin is 22.56 m/d for the K field of EK = 50 m/d and Cv = 1.The Ke value increases to 28.68 m/d while the Cv value decreases to 0.5, and Ke increases to 26.31 m/d while the EK value increases to 100 m/d.Ke of an aquifer is generally affected by the spatial patterns of low-and high-permeability zones.Thus the Ke value of a heterogeneous aquifer with a high average hydraulic conductivity and strong variability is probably less than the Ke value of a weaker heterogeneous aquifer with a low average hydraulic conductivity.That brings significant uncertainty to the estimated K for a heterogeneous aquifer.When EK is assigned to be 10 m/d constantly, the difference of the estimations of Ke between Cv = 1 and Cv = 0.1 is 4.81 m/d for well Wmin, and the corresponding difference for well Wmax

Table 1
The estimated standard error (ESE) and the length of confidence interval (LCI) (t0.95) for two observation wells.Mean and SD represent the mean and standard deviation of results through 100 realizations at each EK and Cv value, respectively.