Cost-Effective Hydraulic Tomography Surveys for Predicting Flow and

Apr 21, 2009 - To make the comparison more specific, this study focuses on the ... c0 represents the initial concentration in the entire modeling doma...
0 downloads 0 Views 789KB Size
Download PDF
Environ. Sci. Technol. 2009, 43, 3720–3727

Cost-Effective Hydraulic Tomography Surveys for Predicting Flow and Transport in Heterogeneous Aquifers C H U E N - F A N I , * ,† T I A N - C H Y I J I M Y E H , ‡,§ AND JUI-SHENG CHEN† Graduate Institute of Applied Geology, National Central University, 300 Chungda Road, Chungli City, Taoyuan 32001, Taiwan, Department of Hydrology and Water Resources, The University of Arizona, Tucson, Arizona 85721, and Department of Resources Engineering, National Cheng Kung University, 1 University Road, Tainan City, 70101 Taiwan

Received August 28, 2008. Revised manuscript received March 27, 2009. Accepted March 30, 2009.

This study shows how a cost-effective hydraulic tomography survey (HTS) and the associated data estimator can be used to characterize flow and transport in heterogeneous aquifers. The HTS is an improved field hydraulic test that accounts for responses of hydraulic stresses caused by pumping or injection events at different locations of an aquifer. A sequential data assimilation method based on a cokriging algorithm is then used to map the aquifer hydraulic conductivity (K). This study uses a synthetic two-dimensional aquifer to assess the accuracy of predicted concentration breakthrough curves (BTCs) on the basis of the K fields estimated by geometric mean, kriging, and HTS. Such K fields represent different degrees of flow resolutions as compared with the synthetically generated one. Without intensive experiments to calibrate accurate dispersivities at sites, the flow field based on the HTS K field can yield accurate predictions of BTC peaks and phases. On the basis of calculating mean absolute and square errors for estimated K fields, numerical assessments on the HTS operation strategy show that more pumping events will generally lead to more accurate estimations of K fields, and the pump locations need to be installed in high K zones to maximize the delivery of head information from pumps to measurement points. Additionally, the appropriate distances of installed wells are suggested to be less than one-third of the ln(K) correlation length in x direction.

1. Introduction The behavior of nonreactive solute transport in aquifers has been recognized to be strongly influenced by spatial distributions of hydrogeological properties. Such properties, including hydraulic conductivity (K) and porosity (n), can create irregular flow paths, which considerably disperse plumes during plume migrations (1-4). Characterizing dispersivities at a particular site is essential to any effort in predicting the propagations of contaminated plumes at that location. Theoretical and experimental investigations have found that the dispersivities at fields are generally greater than those normally obtained in laboratory column tests. As * Corresponding author phone: 886-3-4227151, extension 65874; fax: 886-3-4263127; e-mail: [email protected]. † National Central University. ‡ The University of Arizona. § National Cheng Kung University. 3720

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 43, NO. 10, 2009

a result, laboratory measurements of dispersivities cannot be employed to predict the contaminant transport at field scales. Instead dispersivities evaluated from field scale tracer tests are usually suggested for predicting site-specific transport problems (2-7). Characterization of single solute transport at a particular site commonly employs classical advection dispersion equations (CADEs), associated with calibrated flow field and dispersion coefficients. Because of technical difficulties in capturing velocity variations in detail, many investigators focused on reproducing the transport process of a contaminated plume using the concept of macrodispersion (1, 5, 7). The macrodispersion coefficient (De) evaluated on the basis of such a concept is generally considered to be scale dependent. The determination of scale-dependent De has been the subject of research during the past three decades (see reviews in refs 2, 3, 5, and 8 -11). Using De to characterize contaminant transport is convenient for us because intensive investigations have developed a number of analytical and numerical models for CADEs (see reviews in refs 5 and 8-12). However, experimental investigations commonly found the behavior of anomalous transport in heterogeneous aquifers (13-16). Such behavior generally exhibits early arrival times and long tails for concentration breakthrough curves (BTCs) and cannot be adequately predicted by CADEs (13-15). Therefore, more general but complex approaches such as those that consider fractional order advection dispersion (FAD) and mobile immobile zones (MIZs) have been proposed recently to characterize anomalous transport in heterogeneous aquifers (12, 13, 16, 17). Note that the dispersion coefficients used for FAD or MIZ are still the Des and generally need to be calibrated with observations at sites of interest. Although intensive investigations to characterize contaminant transport in heterogeneous aquifers have developed a number of theories to better interpret plume migrations, our ability to accurately predict contaminant transport today is still limited. This is because of technical limitations to resolving aquifer heterogeneities for most practical problems. Traditional pumping tests associated with measurements of groundwater levels only provide hydraulic properties around wells, and such properties are considered to be locally averaged (18-20). As a result, the information of velocity variations is not taken into account when one uses the averaged flow field to characterize contaminant transport in a heterogeneous aquifer. Although the De at a particular site provides an alternative to compensate our capability to describe the velocity variations for that site, such an approach usually requires routine tracer experiments to determine a De value. Moreover, if the paths of contaminated plumes cover a large area, the determination of De will be costly when one attempts to obtain the representative De values for a variety of scales in the investigation site. Motivated by the technical difficulties of characterizing solute transport in heterogeneous aquifers, this study presents a cost-effective technology that enables us to characterize solute transport with high resolution. More specifically, a hydraulic tomography survey (HTS), which integrates field measurements from multiple cross-hole pumping or injection tests, will be used to improve the resolution of flow patterns at sites of interest. Such crosshole tests are cost effective because the number of installed wells can be significantly reduced. The HTS accounts for the responses of hydraulic stresses caused by pumping or injection tests at different locations (19-25). A sequential successive linear estimator (SSLE) based on a cokriging 10.1021/es8024098 CCC: $40.75

 2009 American Chemical Society

Published on Web 04/21/2009

algorithm is then employed to integrate information from head and direct K measurements and to reconstruct the distribution of K at a particular site (19, 21, 26). Consequently, the resolution of K distribution is improved, and the irregular flow pattern for a specific site is then closely captured. Unlike previous approaches that focus on reproducing the transport process with effective dispersivities, Des for different scales, the presented technology aims to improve the resolution of flow fields. Such flow fields associated with rough estimations of dispersivities and porosity, e.g., the geometric means of dispersivities and porosity, can significantly improve the prediction accuracy of nonreactive contaminant transport in heterogeneous aquifers.

2. Numerical Methods To make the comparison more specific, this study focuses on the steady flow and nonreactive solute transport problem in a vertical, two-dimensional, heterogeneous, and confined aquifer. The following sections will simply present the fundamental concept and estimation procedures associated with corresponding governing equations and assumptions for predicting flow and transport in aquifers. Detailed derivations of equations and numerical methods have been reported by many investigators (8, 11, 12, 19, 21, 26) and will not be repeated here. 2.1. Flow and Transport Model. Assuming steady state flow in a heterogeneous confined aquifer, the groundwater flow equation can be formally expressed as ∇ · [K(x) ∇ h(x)] + Q(x) ) 0

(1)

subject to boundary conditions h(x)|ΓD ) hD

and

[K(x) ∇ h(x)] · n|ΓN ) qN

where h(x) is the hydraulic head, ∇ is the del operator, K(x) is the hydraulic conductivity, and Q(x) represents the sources or sinks that are applied to aquifers. Cartesian coordinate x (x ) 1 and 3) represents x and z directions in the modeling area. hD is the prescribed head at the Dirichlet boundary ΓD, and qN is the specific flux at the Neumann boundary ΓN. Notation n is a unit vector normal to the boundary ΓN. In this study, eq 1 is simulated with the Galerkin finite element algorithm developed by Srivastava and Yeh (27). Likewise, the conservative solute transport equation for saturated porous media is locally governed by the CADE and is written as (8, 9, 11) ∂c(x, t) ) -v(x) ∇ c(x, t) + ∇ · [De(x) ∇ c(x, t)] + Qc(x, t) (2) ∂t subject to initial and boundary conditions c(x, 0)|Ω ) c0, c(x, t)|ΓD ) cD, and [De(x) ∇ c(x, t)] · n|ΓN ) qc where v(x) ) -K(x)∇h(x)/n(x) is the seepage velocity evaluated from flow fields, and n(x) is the effective porosity. The c(x,t) is the volumetric solute concentration measured in the liquid phase, and Qc(x,t) represents the rates where the volumetric solute concentration are injected (source) or extracted (sink) from the aquifer system. Notation c0 represents the initial concentration in the entire modeling domain Ω, cD is the specified concentration at Dirichlet boundary, and qc is the dispersive flux at the Neumann boundary. The De(x) in eq 2 is the macrodispersion coefficient evaluated on the basis of the seepage velocity and is given by (27) De(x) ) [RL(x) - RT(x)]

vivj + D0(x), vj

i, j ) 1, 3

(3)

where RL(x) is the longitudinal dispersivity, RT(x) is the transverse dispersivity, vi and vj are seepage velocities in different directions, and vj is the magnitude of the seepage velocity. Here D0(x) is the effective molecular diffusion. In this study D0(x) in all simulations is assumed to be negligible, this implies that the dispersion is dominated by the advective transport and mechanical dispersion (12). To solve eq 2, we employ the modified method of characteristic (MMOC) developed by Srivastava and Yeh (27). Such a method was also developed on the basis of the Galerkin finite element scheme. To make the numerical assessment reliable, in this study all the forward and inverse simulation cases we used had the same numerical schemes associated with the same spatial and temporal discretizations. 2.2. Stochastic Inversion of Hydraulic Conductivity. The key process in this study is to map a high-resolution K field. Here, a sequential successive linear estimator (SSLE) is employed to conduct the inversion of hydraulic conductivity. The approach of SSLE is an extension of the SLE (successive linear estimator) approach (19, 21, 26). The major advantage of the SSLE is that the sequential processes involved in the SLE can significantly relax the computational burden, especially for a large number of measurements taken from different stress events. The SLE approach is essentially a Bayes approach that uses a cokriging algorithm to evaluate mean parameter fields conditioned on available point data as well as geologic and hydrologic structures (i.e., spatial covariance functions of parameters and hydraulic heads and their crosscovariance functions) (26, 28-31). Different from traditional cokriging algorithms, SSLE uses a linear estimator based on differences between observed and simulated heads successively to update both conditional means and covariances of the estimates such that the nonlinear relation between information and parameters is considered. The general concept of the SSLE can be simply expressed as (4) ξ(1) ) λTζobs + µTεobs * (5) ξ(k+1) ) ξ(k) + ω(k)(εobs - ε(k)) * * (k) where ξ is the kth estimation of parameter (i.e., the hydraulic* conductivity K), ζobs and εobs represent the differences of the parameters (i.e., hydraulic conductivity) and hydraulic heads at observation locations, and λ, µ, and ω(k) are cokriging weighting vectors evaluated by stochastic simulations. On the basis of the observations from Kitanidis and Vomvoris (28) and subsequent work (29, 31), the cokriging can be performed when a general covariance model is assumed and the parameters for that model are found through cross-validation techniques. This process can reduce the dependence on known parameter values at certain points and allows incorporation of geologic intuition in a Bayesian sense. The estimation procedure starts with a weighted linear combination of direct measurements (or primary variable) of the parameters (i.e., the direct K measurements) and hydraulic head (or secondary variable) at different locations to obtain the first estimation of the parameters (i.e., eq 4). The weight µ is calculated on the basis of statistical moments (namely, means and covariances) of parameters, i.e., the covariances of head in space and the cross-covariances between head and hydraulic conductivity. The first estimation of hydraulic conductivity is then used in the mean equation to calculate the head values at observation locations (i.e., forward simulation). The differences between the observed and simulated heads are determined subsequently (i.e., ζobs and εobs). A weighted linear combination of these differences is then used in eq 5 to improve the previous estimations. Iterations between the forward simulation and estimation continue until the improvement of the estimations converged to a specified criterion. Note that the computations of VOL. 43, NO. 10, 2009 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

3721

FIGURE 2. Schematic representation of the synthetic two-dimensional aquifer. Red O are pumping locations. Blue 0 are head observation points. Orange 4 are concentration monitoring points.

FIGURE 1. Generated hydraulic and transport properties for the two-dimensional synthetic aquifer. covariance and cross-covariance of heads and hydraulic conductivity are not straightforward. In this study the firstorder approximation associated with the adjoint state method are employed to evaluate those covariance and crosscovariance (the derivation of adjoint state method can be found in detail in the Supporting Information). The computational cost depends highly on the numbers of nodes in computational domain and the number of head observations. More detailed procedures and the associated computational issues can be found in Zhang and Yeh (21) and Zhu and Yeh (19).

3. Illustrative Examples The objective of HTS is to make the optimal use of wells in an aquifer. The presented cross-hole pumping tests can create more head responses by switching pumping locations at different depths of wells. Because the pumping and monitoring locations are isolated by well packers, the vertical flow patterns will generally exist during pumping tests. Therefore, a vertical profile model is considered here to illustrate the concept and procedures of HTS and the associated simulations. Hereafter, HTS is used to represent the procedures of HTS and the associated SSLE model. 3.1. Generation of the Synthetic Aquifer. To fully control flow and transport conditions in modeling areas, we create our illustrative example synthetically. This study considers a confined, two-dimensional profile model with the area of 60 m × 20 m. The left and right sides of the aquifer are bounded by constant hydraulic heads with the values of 100.6 m and 100 m, respectively. The distributions of hydraulic conductivity, porosity, and longitudinal and transverse dispersivities are generated using a spectral algorithm associated with known means, variances, and correlation lengths in an exponential covariance model (32). Figure 1 shows the generated distributions of such hydrogeological parameters. In Figure 1, the geometric mean of K is 1.0 m/day, 3722

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 43, NO. 10, 2009

and the geometric mean of porosity (n) is 0.3 m3/m3. The variance of ln(K) is 1.0, and the variance of ln(n) is 0.01. Here the geometric mean value of longitudinal dispersivity (RL) is 1.0 m, while the geometric mean value of transverse dispersivity (RT) is 0.1 m. The variances of ln(RL) and ln(RT) are assumed to be 0.1 to make ranges of the generated dispersivities within same orders. In addition, all the generated random fields have same correlation length (λ) of 40 m for the x direction and 5 m for the z direction. Note that all the random fields are generated independently, representing the worst scenario for predicting flow and transport in such an aquifer. The distributions of K, porosity, and longitudinal and transverse dispersivities are considered to be the true fields in the modeling area. Our purpose here is to predict the flow and solute transport using limited head observations from wells. 3.2. Generation of Synthetic Observations. To prepare the synthetic measurements from the modeling area, we employ a Galerkin finite element scheme to obtain a steady state flow field (i.e., solution to eq 1). The modeling area is discretized into a total of 1200 elements (60 × 20). Each element is 1 m × 1 m in length. Figure 2 is the schematic diagram of the experimental setup for the locations of pumping, monitoring, and the concentration source zone. Three wells are installed at x ) 10 m, 30 m, and 50 m, respectively. Each well is considered to be a fully penetrated well, i.e., the well screen is opened through the entire aquifer. Moreover, the direct K measurements along each well are assumed to be known on the basis of the sampled well logs. A total of 27 direct K measurements (nine for each well) will be used for estimating kriging K and HTS K fields. Different from the kriging estimation, the HTS couples direct K measurements with head values observed from different pumping stresses (i.e., cross-hole pumping tests). To make the optimal use of information from well hydraulic tests, there are 10 pneumatic packers installed in each well to isolate the flows from vertical directions. Consequently, during each pumping test there are nine head changes in each well (marked by squares in Figure 2). In the illustrative example, eight pumping events are specified and are marked by circles in Figure 2. Except for the head values at pumping locations, eight pumping events can create a total of 216 head changes in the entire modeling area. The eight pumping events associated with the head and K measurements at well locations are assumed to be the known information used later for predicting aquifer hydraulic conductivity. Same as the finite element mesh for predicting the K field, Figure 2 also shows the locations for source zone and concentration monitoring ports. Three monitoring locations are located in the middle part of each well (i.e., z ) 10m). On the basis of the same spatial and temporal spaces, boundary conditions, and the initial location of the released plume, we conducted MMOC simulations based on the flow fields estimated with geometric mean, kriging, HTS, and true K fields. Here the BTCs obtained on the basis of the generated

To quantitatively compare the estimation errors based on the limited measurements, we further draw the scatterplots of the estimated versus true K fields (Figure 5) and velocity fields (Figures S1 and S2 of the Supporting Information for x and z direction velocities, Vx and Vz), along with the mean absolute error normal L1 and mean square error normal L2, which are defined as n

L1 )

FIGURE 3. (a) Estimated K field using kriging with 27 direct K measurements along the three boreholes. (b) Estimated K field using HTS. K, RL, RT, and n fields (Figure 1a-d) are again considered to be our true BTCs for comparison purposes. Other flow fields obtained on the basis of geometric mean, kriging, and HTS K fields use geometric mean values of RL, RT, and n to be the input transport parameters.

4. Results and Discussion 4.1. Predictions of Hydraulic Conductivity K. On the basis of the known head measurements with eight pumping events, we can perform the inverse modeling using SSLE. In this study, the SSLE has been incorporated with an interactive interface VSAFT (variably saturated flow and transport in 2D) developed by using the Microsoft Visual C++ integrated development environment (IDE). The convergence criteria for the estimation processes of SSLE are presented in detail in the Supporting Information. Kriging is a typical algorithm used for many site investigations. To prepare the statistical parameters for kriging estimation, this study employs the GAMV program (33), associated with a nonlinear least-squares algorithm to obtain the best fit of variances and ranges for the exponential semivariogram model. On the basis of the 27 K observations from three wells, we obtain the ranges of 20.1 m and 5.1 m for x and z directions. The ln(K) variance values for the two directions are, respectively, 0.53 and 0.77. Figure 3 compares the hydraulic conductivity field estimated by the kriging estimation (Figure 3a) and that estimated by HTS (Figure 3b). Comparing the K fields in Figure 3 with the one in Figure 1a, we found that the HTS can closely capture the true K distribution. However, the kriging estimation only gives a relatively smooth K distribution pattern. Many investigations have found similar results that the kriging estimation provides smooth K distributions (19, 21). On the basis of different K fields, the head distributions are shown in Figure 4. In Figure 4, the head field is very accurate to the HTS compared with the true case. However, the flow field obtained by the kriged K field show a slight difference.

FIGURE 4. Steady state flow fields estimated with true, kriging, and HTS K fields.

∑

1 |χ - χ*i | n i)1 i

n

L2 )

and

∑

1 (χ - χ*i )2 n i)1 i

(6)

where χi and χ*i represent the true and estimated parameter (i.e., ln(K), Vx, or Vz). In eq 6, i indicates the element number, and n is the total number of elements for the modeling area. In addition to L1 and L2, the correlations based on the estimated and true parameters are also calculated. The correlation has the following mathematical formation n

∑ (χ - χj)(χ* - χ*) i

Correl(χ, χ*) )

i

i)1

∑ n

(7) (χi - χj)2(χ*i - χ*)2

i)1

where χj and χ* represent the geometric means of parameters. Comparing kriging and HTS estimations with the true conductivity values at each element, the results indicate that our algorithm, based on HTS and the associated SSLE model, produces accurate estimations of K and the velocities in x and z directions (Figure 5b and Figures S1b and S2b of the Supporting Information). However, the kriging algorithm provides bias estimations of K (Figure 5a) and velocity fields (Figures S1a and S2a of the Supporting Information). Such bias estimations for kriging can be caused by the source of 27 K direct measurements. Here most K measurements used for the kriging estimation are relatively high K. Therefore, kriging estimation gives estimation biased to high K values. To further evaluate the spatial structure of estimated K fields, we calculated the x direction variograms for three different K fields using a GAMV program in the GSLIB geostatistical software (Figure S3 of the Supporting Information) (33). We found that the K field estimated by HTS shows good agreement with the synthetically generated one. Although the minimal x direction distance of available K measurements at well locations is 20 m, the HTS results show good predictions on the range and variance for the K statistical structure. Because the kriging estimation fixes the variogram model and statistical structure, the statistical structure of the kriged K field may reproduce well the x direction variance and range. However, a great difference of variogram values are calculated when the distances are greater than 20 m (i.e., kriging range). The results in Figures 4 and 5 and Figures S1-S3 of the Supporting Information reveal that good fittings on head measurements may not lead to good calibrations on hydraulic conductivity and velocity fields. In order to obtain a better estimation of K and velocity fields, the correlation between head and conductivity must be considered simultaneously. This is HTS involvment in the data assimilation process. 4.2. Predictions of Transport with Different Degrees of K Resolutions. Many previous investigations have recognized that the spreading of transport behavior is strongly influenced by spatial variation in hydraulic conductivity K (5, 8, 11). Therefore, this study will directly focus on the assessment of different degrees of K resolutions in solute transport predicting. We have estimated a K field with HTS and a kriging estimation in Predictions of Hydraulic Conductivity K. Here we further employ these predicted K fields to assess the VOL. 43, NO. 10, 2009 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

3723

FIGURE 5. Scatterplots of estimated K fields versus the true one. (a) K field estimated by kringing algorithm. (b) K field estimated by HTS. improvement in predicting solute transport. In addition to the true, kriged, and HTS K fields, the geometric mean K field is added in transport simulations for illustration purposes. Because of the limited information of aquifer properties, such K fields represent three different degrees of accuracy (or resolution). The geometric mean K field has the lowest resolution because the single mean value of ln(K) is applied to the entire aquifer. Kriging estimation incorporates the direct K measurements and the associated geostatistical structure for an aquifer. Therefore, the resolution of the kriged K field is higher than that of the geometric mean K field. The HTS integrates information from K measurements and head observations from pumping tests and provides the highest resolution of the K field. Although many investigations have focused on the evaluations of dispersivity RL and RT on the basis of the observations of BTCs, plume size, travel times, and heterogeneity of ln(K) (2, 3, 10, 34), the estimated dispersivities using these proposed approaches show a variety of ranges. For example, if we solely rely on the true BTCs from the middle point of each well (i.e., z ) 10 m) and adjust the CADE model with a geometric mean flow field to fit the true BTCs, we obtain the representative RL values of 0.8, 1.5, and 2.1 m for locations at x ) 10, 30, and 50 m, respectively. This result is consistent with many previous studies indicating that the dispersivity values increase with travel distances. Furthermore, the macrodispersivity formulas proposed by Gelhar et al. (35) will lead to a longitudinal dispersivity value of 6.3 m and a transverse dispersivity value of 0.4 m. Such results require that the correlation length (λ) of ln(K) fluctuations must be significantly larger than the values of dispersivities RL and RT (3). In addition, the formulas are only applicable for late times (i.e., long travel distances). Other formulas such as Neuman (2), Xu and Eckstein (34), and Schulze-Makuch (10) give the longitudinal dispersivity values of 7.9, 1.86, and 1.065 m, respectively. All of these formulas were based on the experimental data organized by Gelhar et al. (3). In this study, the dispersivities evaluated by different formulas (i.e., formulas proposed by Gelhar et al. (35), Neuman (2), Xu and Eckstein (34), and Schulze-Makuch (10)) are consistently of the same orders (i.e., 1 m for RL and 0.1 m for RT). We then use the geometric mean values of RL and RT to be our dispersivities applied to cases with geometric mean, kriging, and HTS K fields. Except for the true K field case, the distributions of porosity (n) for cases with geometric mean, kriging, and HTS K fields are also assumed to include the geometric mean values for all of the modeling areas. For practical applications, such mean values can be obtained approximately by analyzing direct K measurements from well logs or well hydraulic tests. Although the analyses of well 3724

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 43, NO. 10, 2009

FIGURE 6. Comparison distances.

of

BTCs

at

different

monitoring

logs and hydraulic tests may introduce uncertainties on estimating the geometric mean values, the effects of such uncertainty on the final predictions of K fields are not considered in this study. To compare the results based on different K fields, the monitoring location at the middle portion of each well (i.e., z ) 10 m) was selected. Figure 6 shows the BTCs obtained at the middle locations of three wells. Simulation results show that the K field estimated by HTS obtained the best prediction of BTCs at all three locations when compared with the true case. The flow field based on the HTS K field can yield accurate BTC peaks and phases at different locations. Here, the scale-

TABLE 1. Simulation Scenarios for Examples to Assess the Effects of Different Degrees of K Heterogeneity on Final Predictions of HTS K Fieldsa test geostatistical parameter test parameter values fixed parameters

ln(K) variance

ln(K) x correlation length (m)

0.1, 0.5, 1.0, and 2.0 ln(K) x correlation length ) 40 m ln(K) z correlation length ) 5 m

20, 40, 60, and 80 ln(K) variance ) 1.0 ln(K) z correlation length ) 5 m

a

The second column shows the test cases with varied variances and fixed values of ln(K) x and z correlation lengths, while the third column shows the test cases with varied ln(K) x correlation lengths and fixed ln(K) variance and z correlation lengths.

dependent dispersion coefficient becomes insignificant because the dispersion mechanism has been mainly involved in the more accurate flow field. The BTC results based on the kriged K field show wrong predictions of BTC peaks and underestimations of phase speeds. Moreover, the phase differences increase with the distances between monitoring points and the original plume location. Note that the variogram model and associated variances and ranges used in the kriging estimation are based on the 27 K measurements from well logs. For most practical problems, information such as the variogram model and parameters are often limited. Therefore, in a realistic situation, the prediction of BTC based on the kriged K field may obtain even lower accuracy than that shown in this study. The BTCs simulated based on geometric mean K fields present acceptable phase speeds but show overestimations of the BTC peak values. The overestimations of peak values increase with increasing monitoring distances from the original plume location (around 50% and 80% for x ) 30 and 50 m, respectively). This result is consistent with the conclusion from previous investigations in that the dispersivity values need to be increased with increasing travel distances to adequately reflect the plume spreading caused by flow variations in heterogeneous aquifers (2, 3). Although in most practical problems the plume patterns are generally unavailable, it is of interest to explore more on the concentration distributions rather than just the BTC point values. To clearly compare the plume patterns, we plotted the concentration distributions at time ) 500 and 1000 days in Figure S4 of the Supporting Information. The results based on true and HTS K fields show similar concentration patterns. However, the patterns obtained based on the kriged and geometric mean K fields are significantly different from those obtained based on the true K field. The results from Figures 6 and S4 of the Supporting Information indicate that the measurement and analysis procedures can significantly affect the final predictions of solute transport. With higher resolutions of K fields, the solute transport processes can be closely captured by using rough parameter values such as the mean values of porosity and dispersivities. In addition, investigators who use numerical models to characterize flow and transport can reduce the effort on conducting tracer tests at fields for evaluating transport parameters. Our comparison results imply that the K distribution in our numerical example is relatively more important than other parameters such as porosity and dispersivities to perform a good prediction of solute transport. 4.3. HTS Sampling Strategy. Predictions of Transport with Different Degrees of K Resolutions has illustrated the accuracy of BTCs and plume patterns obtained from HTS K field. The K field used for the transport prediction is based on eight pumping events associated with 216 head measurements. For the most practical applications, it is of interest to understand the locations and the number of pumping events that may affect the final estimations of the K fields. Additionally, the available direct K measurements and the degrees of K heterogeneity are site specific and can significantly affect the procedures and accuracy to characterize

flow and contaminant transport. Here, we aim to provide a general insight into the sampling strategy of the HTS that can make characterization of flow and transport more efficient. A variety of cases are used to illustrate the accuracy of K field estimations based on different sets of pumping strategies, numbers of direct K measurements, and degrees of K heterogeneity. 4.3.1. HTS Pumping Strategy. On the basis of the same K field in previous sections, all pumping scenarios are assumed to have the same head measurement ports for each pumping event (i.e., total of 27 head measurements for each pumping event) and have the same pumping rate -10 m3/ day for each pump. Results and analysis for different pumping scenarios are presented in detail in Figure S5 of the Supporting Information. Here, we summarize the results of the HTS pumping strategy and have the following observations: (1) The pumping locations where more stress (pressure) produces changes for head measurements will lead to a more accurate estimation of K fields, (2) more pumps associated with more pumping events will generally lead to a more accurate estimation of the K field; however, the appropriate number of pumps depends on the pumping rate and the K distribution in an aquifer, and (3) Following the observation in (2), the pump locations should be in high K zones (or layers) of an aquifer to maximize the delivery of head information from pumps to measurement points. On the basis of the observations of the numerical experiment for the HTS operation strategy, we found that the appropriate number of pumps for the presented aquifer is two or three placed in the high K zone, and these pumps have to be operated individually, i.e., cases (f) and (j) in Figure S5 of the Supporting Information. This conclusion agrees well with those proposed by Snodgrass and Kitanidis (29) and Yeh and Liu (26). 4.3.2. HTS Sampling Strategy for Different Degrees of K Heterogeneity. To provide systematical assessment for the HTS in estimating K fields with different degrees of K heterogeneity, we conducted numerical experiments based on the same modeling domain (60 m × 20 m), boundary conditions, well locations (x )10, 30, and 50 m), and head observation locations (nine observations for each well). However, the K fields are generated independently using a spectral random field generation algorithm associated with an exponential covariance model. On the basis of the analysis of results from the pumping strategy, the assessments of the sampling strategy only use three pumping locations (at z ) 2.5 m for each well) associated with three pumping events. The variances and correlation lengths to generate K fields are allowed to vary in a variety of ranges. Table 1 lists the scenarios of K heterogeneity for the test examples. For each variance or correlation length value, we generate 25 equally likely unconditional K fields. Similar to the previous sections in generating synthetic head observations, each K field is then estimated by the HTS on the basis of the head observations from 27 ports. The L1 and L2 values for 25 estimations are then averaged to show the mean performance of HTS for a specified variance or correlation length value. In this study a total of 200 K fields are generated and then VOL. 43, NO. 10, 2009 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

3725

estimated by the HTS. Note that we consider only the correlation lengths in the x direction because the vertical observations can always be increased by installing more packers in a well along the vertical direction. The results for assessments of different degrees of K heterogeneity are summarized in Figure S6 of the Supporting Information. As is expected, the estimation accuracy decreases with the increase in K variance values (Figure S6a of the Supporting Information). Similar to the limitation of other first-order methods, the adjoint state method used here to obtain head covariance and head and K cross correlations may not be able to represent well the correlations in highly heterogeneous K fields. For the changes of ln(K) correlation lengths in the x direction, we have conclusions similar to the work of Illman et al. (23). The estimation accuracy will be improved with the increased K correlation lengths in the x direction (Figure S6b of the Supporting Information). In Figure S6b of the Supporting Information, the simulation results reveal that the simulation errors are not significant when the ln(K) x correlation lengths are greater than 60 m. This distance covers three wells in the modeling area. On the basis of the analysis of different correlation lengths of K fields, this study suggests that an appropriate distance of installed wells is less than one-third of the ln(K) x correlation length. 4.3.3. Effect of Direct K Measurements on HTS Estimation Results. On the basis of the pumping strategy, here we consider cases (f) and (j) to assess the degradation of K fields due to the decrease in direct K measurement numbers. Case (f) presents two pumping locations with two pumping events, while case (j) has three pumping locations and conducts three pumping events individually. Here the number of direct K measurements is changed from 27 to 15, 9, 6, 3, 1, and 0. For a scenario that involves 15 direct K measurements, there are five K measurements in each well, and the measurement locations are distributed along the z direction (i.e., z ) 2, 6, 10, 14, and 18 m). Similar situations are applied to scenarios with nine, six, and three direct K measurements. The direct K measurements are three (i.e., z ) 2, 10, and 18 m), two (i.e., z ) 6 and 14 m), and one (i.e., z ) 10 m) for each well. For a special scenario such as a one K measurement, the K measurement location is in the middle portion (i.e., z ) 10 m) of the second well. The estimation results for cases (f) and (j) for different numbers of direct K measurements are summarized in Figures S7 and S8 of the Supporting Information. The results in Figure S7 of the Supporting Information show that the estimation errors increase gradually with decreased K measurement numbers. However, in each case the differences are not significant for scenarios with small numbers of total K measurements such as three, one, and zero. The worst scenario is a zero K measurement for each case. Under this scenario, the L1 and L2 for case (f) are 0.385 and 0.230, while the L1 and L2 for case (j) are 0.352 and 0.192. We found that more pumps associated with corresponding pumping events will give better results, although no direct K measurement is used for K estimations. Figure S8 of the Supporting Information shows the results with a zero K measurement for cases (f) and (j). The results present that an additional pumping location and pumping event applied to the second well in case (j) will improve the estimation results. Comparing the results in Figure S8 of the Supporting Information with those in panels f and j of Figure S5 of the Supporting Information, we found that the improvement with more pumping locations and events is significant when no direct K measurement is available. However, with sufficient direct K measurements (Figure S5f,j of the Supporting Information), more pumping locations and events have relatively little impact on the final K estimation. 3726

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 43, NO. 10, 2009

Acknowledgments This research was supported in part by the National Science Council of the Republic of China under Contract NSC 972116-M-008-002. T.-C.J.Y. acknowledges support from the National Cheng-Kung University during his visit in 2006, 2007, and 2008. Support from the Strategic Environmental Research and Development Program (SERDP) as well as funding from the National Science Foundation (NSF) through Grants EAR0229713, EAR-0229717, IIS-0431069, IIS-0431079, EAR0450336, and EAR-0450388 are acknowledged. We thank three anonymous reviewers for their valuable comments.

Supporting Information Available Figures S1-S8 and detailed descriptions of HTS convergence criteria, sensitivity analysis, and pumping strategies. This material is available free of charge via the Internet at http:// pubs.acs.org.

Literature Cited (1) Dagan, G. Time-dependent macrodispersion for solute transport in anisotropic heterogeneous aquifers. Water Resour. Res. 1988, 24, 1491–1500. (2) Neuman, S. P. Universal scaling of hydraulic conductivity and dispersivities in geological media. Water Resour. Res. 1990, 26 (8), 1749–1758. (3) Gelhar, L. W.; Welty, C.; Rehfeldt, K. R. A critical-review of data on field-scale dispersion aquifers. Water Resour. Res. 1992, 28 (7), 1955–1974. (4) Ptak, T.; Piepenbrink, M.; Martac, E. Tracer tests for the investigation of heterogeneous porous media and stochastic modelling of flow and transport - a review of some recent developments. J. Hydrol. 2004, 294 (1-3), 122–163. (5) Dagan, G.; Neuman, S. P. Subsurface Flow and Transport: A Stochastic Approach; Cambridge University Press: New York, 1997. (6) Mallants, D.; Espino, A.; Van Hoorick, M.; Feyen, J.; Vandenberghe, N.; Loy, W. Dispersivity estimates from a tracer experiment in a sandy aquifer. Ground Water 2000, 38 (2), 304– 310. (7) Sudicky, E. A. A natural gradient experiment on solute transport in a sand aquifer: Spatial variability of hydraulic conductivity and its role in the dispersion process. Water Resour. Res. 1986, 22 (13), 2069–2082. (8) Dagan, G. Flow and Transport in Porous Formations; SpringerVerlag: New York, 1989. (9) Rubin, Y. Applied Stochastic Hydrogeology; Oxford University Press: New York, 2003. (10) Schulze-Makuch, D. Longitudinal dispersivity data and implications for scaling behavior. Ground Water 2005, 43 (3), 443– 456. (11) Gelhar, L. W. Stochastic Subsurface Hydrology; Prentice-Hall: Upper Saddle River, NJ, 1993. (12) Zheng, C.; Bennett, G. D. Applied Contaminant Transport Modeling. 2nd ed.; Wiley-Interscience: New York, 2002. (13) Benson, D. A.; Wheatcraft, S. W.; Meerschaert, M. M. The fractional-order governing equation of Levy motion. Water Resour. Res. 2000, 36 (6), 1413–1423. (14) Berkowitz, B.; Scher, H.; Silliman, S. E. Anomalous transport in laboratory-scale, heterogeneous porous media. Water Resour. Res. 2000, 36 (1), 149–158. (15) Dentz, M.; Cortis, A; Scher, H; Berkowitz, B. Time behavior of solute in heterogeneous media: transition from anomalous to normal transport. Adv. Water Res. 2004, 27 (2), 155–173. (16) Cortis, A.; Berkowitz, B. Anomalous transport in “classical” soil and sand columns. Soil Sci. Soc. Am. J. 2004, 68 (5), 1539–1548. (17) Feehley, C. E.; Zheng, C. M.; Molz, F. J. A dual-domain mass transfer approach for modeling solute transport in heterogeneous aquifers: Application to the macrodispersion experiment (MADE) site. Water Resour. Res. 2000, 36 (9), 2501–2515. (18) Sellwood, S. M.; Healey, J. M.; Birk, S.; Butler, J. J. Direct-push hydrostratigraphic profiling: Coupling electrical logging and slug tests. Ground Water 2005, 43 (1), 19–29. (19) Zhu, J. F.; Yeh, T. C. J. Characterization of aquifer heterogeneity using transient hydraulic tomography. Water Resour. Res. 2005, 41 (7), W07028, doi:10.1029/2004WR003790. (20) Illman, W. A.; Craig, A. J.; Liu, X. Practical issues in imaging hydraulic conductivity through hydraulic tomography. Ground Water 2008, 46 (1), 120–132.

(21) Zhang, J. Q.; Yeh, T. C. J. An iterative geostatistical inverse method for steady flow in the vadose zone. Water Resour. Res. 1997, 33 (1), 63–71. (22) Hao, Y.; Yeh, T.-C. J.; Illman, W. A.; Ando, K.; Hsu, K.-C. Hydraulic tomography for detecting fracture connectivity. Ground Water 2008, 46 (2), 183–192. (23) Illman, W. A.; Liu, X.; Takeuchi, S.; Yeh, T.-C. J.; Ando, K.; Saegusa, H. Hydraulic tomography in fractured granite: Mizunami underground research site, Japan. Water Resour. Res. 2009, 45, W01406, doi:10.1029/2007WR006715. (24) Liu, X.; Illman, W. A.; Craig, A. J.; Zhu, J.; Yeh, T.-C. J. Laboratory sandbox validation of transient hydraulic tomography. Water Resour. Res. 2007, 43, W05404, doi:10.1029/2006WR005144. (25) Straface, S.; Yeh, T.-C. J.; Zhu, J.; Troisi, S.; Lee, C. H. Sequential aquifer tests at a well field, Montalto Uffugo Scalo, Italy. Water Resour. Res. 2007, 43, W07432, doi:10.1029/2006WR005287. (26) Yeh, T. C. J.; Liu, S. Y. Hydraulic tomography: Development of a new aquifer test method. Water Resour. Res. 2000, 36 (8), 2095–2105. (27) Srivastava, R.; Yeh, T. C. J. A 3-dimensional numerical model for water flow and transport of chemically reactive solute through porous-media under variably saturated conditions. Adv. Water Res. 1992, 15 (5), 275–287. (28) Kitanidis, P. K.; Vomvoris, E. G. A geostatistical approach to the inverse problem in groundwater modeling (steady state) and

(29) (30) (31)

(32)

(33) (34) (35)

one-dimensional simulations. Water Resour. Res. 1983, 19 (3), 677–690. Snodgrass, M.; Kitanidis, P. K. Transmissivity identification through multidirectional aquifer stimulation. Stochastic Hydrol. Hydraul. 1998, 12 (5), 299–316. Yeh, T.C. J.; Lee, C. H. Time to change the way we collect and analyze data for aquifer characterization. Ground Water 2007, 45 (2), 116–118. Fienen, M. N.; Clemo, T. M.; Kitanidis, P. K. An interactive Bayesian geostatistical inverse protocol for hydraulic tomography. Water Resour. Res. 2008, W00B01, doi:10.1029/ 2007WR006730. Gutjahr, A. Fast Fourier Transform for Random Field Generation; New Mexico Project Report for Los Alamos Grant Contract 4-R58-2690R; New Mexico Institute of Mining and Technology: Socorro, NM, 1989. Deutsch, C. V.; Journel, A. G. GSLIB Geostatistical Software Library and User’s Guide., 2nd ed.; Oxford University Press: New York, 1997. Xu, M. J.; Eckstein, Y. Use of weighted least-squares method in evaluation of the relationship between dispersivity and fieldscale. Ground Water 1995, 33 (6), 905–908. Gelhar, L. W.; Axness, C. L. Three-dimensional stochastic analysis of macrodispersion. Water Resour. Res. 1983, 19 (1), 161–180.

ES8024098

VOL. 43, NO. 10, 2009 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

3727