Stochastic inversion of sequential hydraulic tests for transient and highly permeable unconfined aquifer systems

Introduction Conclusions References


Conclusions References
Tables Figures

Back Close
Full

Abstract
A hydraulic tomography survey (HTS) is a conceptually improved technique that has been recognized to be efficient for estimating high-resolution aquifer parameters.
Based on the concept of HTS, this study presents a modified stochastic inverse model for estimating hydraulic conductivity (K ) and specific yield (Sy) in shallow and highly permeable unconfined aquifers.A well field with 15 fully screened wells was developed for the purpose of model implementations.In this study a synthetic example was first employed to assess the accuracy of the inverse model.We then implemented the model to field-scale, cross-hole injection tests in a shallow and highly permeable unconfined aquifer near the middle reach of the Wu River in central Taiwan.To assess the effect of constant head boundary conditions on the estimation results, two additional modeling domains were evaluated based on the same field data from the injection tests.Results for the synthetic example show that the modified inverse model can reproduce well the predefined geologic features of the unconfined aquifer.The inverse model can estimate accurately the ln K patterns and magnitudes.However, slightly fewer details of the ln Sy field are obtained due to the insensitivity of transient hydraulic stresses for specified sampling times.Model implementations of field-scale injection tests show that the model can estimate ln K and ln Sy fields with high spatial resolution.The estimated K and Sy values for the test site vary by one order of magnitude, indicating a relatively homogeneous aquifer for the tested well field.Results based on three different modeling domains show similar patterns and magnitudes of ln K and ln Sy near the well locations.This result suggests that the case with domain 40 m × 20 m should be sufficient for the injection tests at the well field.

Introduction
Traditional approaches to obtain hydraulic conductivity (K ) and the storage coefficient (S) for aquifers rely on slug and/or pumping tests and the associated well-defined analytical or semianalytical solutions.These approaches can efficiently estimate either averaged aquifer parameters (pumping/injection tests) between wells or point measurements (slug tests) near wells.These approaches are recognized to be insensitive to aquifer strata connectivity because the sizes of representative control volumes are relatively small (Bohling et al., 2007;Bohling and Butler, 2010).Previous investigations have recognized that the great challenge to characterize groundwater flow and transport remains the limited measurements.The small number of measurements for most practical problems have thus motivated the development of sequential cross-hole tests and associated inverse models in order to maximize the usefulness of measurements (e.g., Bohling et al., 2007;Cardiff et al., 2009;Zhu andYeh, 2005, 2006;Ni et al., 2009).Tomographic approaches such as hydraulic and pneumatic tomography surveys have been proven to be the feasible techniques for estimating aquifer properties with high spatial resolution on different scales (e.g., Gottlieb and Dietrich, 1995;Yeh and Zhang, 1996;Yeh and Liu, 2000;Vesselinov et al., 2001;Liu et al., 2002;Tsai et al., 2005;Zhu andYeh, 2005, 2006;Hernandez et al., 2006;Bohling et al., 2007;Illman et al., 2008;Ni and Yeh, 2008;Cardiff et al., 2009;Ni et al., 2009;Bohling, 2009;Yeh et al., 2009;Bohling and Butler, 2010;Brauchler et al., 2010;Cardiff and Barrash, 2011).Over the years, many models have been developed for analyzing hydraulic tests on the scales of laboratory or field experiments (e.g., Vesselinov et al., 2001;Liu et al., 2002Liu et al., , 2007;;Zhu and Yeh, 2005;Straface et al., 2007Straface et al., , 2011;;Cardiff et al., 2009Cardiff et al., , 2012;;Illman et al., 2010Illman et al., , 2012;;Schoniger et al., 2012).With the improvement of measuring and computing technologies, previous investigations have focused on the developments of efficient inverse models for problems with realistic scales and complexities.Gottlieb and Dietrich (1995) were the first to propose the concept of hydraulic tomography.They employed a least-squares optimization algorithm to estimate Introduction

Conclusions References
Tables Figures

Back Close
Full the spatial distribution of hydraulic conductivity in a synthetic two-dimensional saturated aquifer assuming steady-state flow.A similar tomography concept but using a stochastic approach was used in the study of Kitanidis (1995).Yeh and Liu (2000) developed a sequential stochastic inverse model called the sequential successive linear estimator (SSLE).Their model can be applied to hydraulic tests for aquifers under confined and steady-state conditions.The optimization algorithm in the SSLE model relies on the solution of sensitivity equations for flow and iteratively updating the cokriging weightings to obtain a conditional K field.The advantage of the sequential inclusion of pumping tests is the computational efficiency because the sequential pumping or injection events are calculated separately and are integrated by inner iteration loops.The sizes of the matrices for stochastic simulations are significantly reduced.Zhu and Yeh (2005) improved the SSLE to solve for transient flow problems.The SSLE inverse model has been tested successfully for laboratory experiments (Liu et al., 2002(Liu et al., , 2007;;Illman et al., 2008Illman et al., , 2010Illman et al., , 2012) ) and field studies (Straface et al., 2007(Straface et al., , 2011;;Illman et al., 2009).Detailed distributions of hydraulic conductivity (K ) and specific yield (Sy) in unconfined aquifers are critical in predicting near-surface contaminant transport and quantifying surface water-ground water interactions.It is recognized that the estimations of aquifer properties in unconfined aquifers are difficult due to the complexity of drawdown behavior to be fitted and the nonlinearity of the flow equation to be solved (Mao et al., 2011).A recent study by Cardiff et al. (2009)  unconfined aquifer systems.Their numerical experiments showed that the estimation of K is slightly degraded if joint estimations of storage parameters are required.However, such an effect is small compared with the degradation of estimates related to adopting a wrong geostatistical structure.The developed numerical model was implemented to a field-scale inverse problem at BHRS (Cardiff et al., 2012).
With the development of measurement technologies and associated analyzing algorithms, we expect that data in laboratories or at sites can be obtained with high quality and quantity.The sequential hydraulic test is a practically feasible technique for high-resolution aquifer characterizations.The model used for such aquifer test needs to meet the operation procedures and high-resolution requirements (Ni and Yeh, 2008;Ni et al., 2009).Our objectives of this study are three-fold: (1) to develop a cross-hole inverse model for characterizing groundwater flow in transient and unconfined aquifers, (2) to develop a well field and conduct a field-scale experiment following the characterization procedures used in the inverse model, and (3) to implement the developed inverse model to characterize the hydraulic tests from the well field.More specifically, in the study we modified the SSLE model developed by Zhu and Yeh (2005) and Ni and Yeh (2008) to estimate two-dimensional depth-averaged distributions of K and Sy in transient unconfined aquifer systems.This paper begins with the derivation of the SSLE for transient flow in an unconfined aquifer system.We then illustrate the developed inverse model with a synthetic example.Observations from field-scale experiments are employed to test the developed model in estimating ln K and ln Sy distributions for the well field in the middle reach of the Wu River in central Taiwan.

Flow equations
Consider transient groundwater flow in a depth-averaged, two-dimensional, and heterogeneous unconfined aquifer system.The governing equation for the flow in such an unconfined aquifer system can be represented by Introduction

Conclusions References
Tables Figures

Back Close
Full subject to the initial condition and the boundary conditions where h = h(x, t) is the hydraulic head [ ], x is the vector of Cartesian coordinates (x, y) T , and t is the time.The aquifer parameter T (x) = K (x)h(x) is the transmissivity where K (x) is the hydraulic conductivity [ /t], which is assumed homogeneous throughout the depth of the aquifer.The Sy = Sy(x) is the specific yield.q s = q s (x, t) is the depth-averaged source or sink term for the aquifer system.In Eqs.
(3) and (4), h f (x, t) is the specified head on the boundary segment Γ D , while q b = q b (x, t) is the flux across the Neumann boundary Γ N and n is an outward unit vector normal to Γ N .Equation ( 1) is a nonlinear equation for transient two-dimensional unconfined aquifer systems.To solve this equation, an iteration algorithm involving inner iterative loops at each time step is usually employed to obtain the solution of h.Such nonlinearity should be considered for the inversion procedures and will be discussed later in the following sections.

Optimization algorithm
Over the years many solutions to the inverse problem based on pumping or slug test data have been developed.These methods use observed hydraulic head h(x, t) to estimate the distribution of hydraulic conductivity K (x) based on the minimization of an Introduction

Conclusions References
Tables Figures

Back Close
Full objective function, [h(x, t) − ĥ(x, t)] 2 , where ĥ(x, t) is the simulated hydraulic head (e.g., Mclaughlin and Townley, 1996;Zimmerman et al., 1998;Vargas-Guzmân and Yeh, 1999).The approach of the SSLE follows the same concept but focuses on minimizing the residual head and sequentially incorporates the pumping or injection data from different stresses.Iterations are used to ensure that the estimations best fit the observations from different stresses.The detailed description for the sequential process can be found in Yeh and Liu (2000) and also in Zhu and Yeh (2005).
To estimate flow properties such as K and Sy in an unconfined aquifer, we assume that the spatial distributions of natural logarithm of K and Sy (i.e., ln K and ln Sy) in unconfined aquifers are stochastic processes.The K and Sy values are composed of mean values and the associated small perturbations.Mathematical formulas can be written as ln K = f +f and ln Sy = g+g , where f and g are the mean values and f and g are the perturbations.The transient head responses to an injection or pumping test are also treated as stochastic processes and can be represented by h = h + h , where h is the mean head and h is the perturbation.Substituting these stochastic variables into Eq.( 1), taking the conditional expectation, and conditioning on observations of hydraulic head and parameters such as K and Sy results in the following conditional mean equation for groundwater flow in unconfined aquifer systems: where K c , h c , and Sy c are conditional effective hydraulic conductivity, hydraulic head, and specific yield, respectively.Please note that the source/sink term q s (i.e., the pumping/injection in the aquifer) in this study is considered to be a deterministic process.
To obtain the K c and Sy c in Eq. ( 5), a cokriging algorithm is employed to integrate the available point data (i.e., heads, K , and Sy) with geostatistical structure.More specifically, in this study the sequential cokriging approach updates the conditional means and covariances based on differences between observed and simulated heads.The Introduction

Conclusions References
Tables Figures

Back Close
Full following steps present the estimation procedure for an unconfined aquifer system.
Note that the conceptual model and the numerical parameters are assumed to be ready for the estimation procedures.
1. Select a stress event and the associated head observations.Apply initial guesses of K and Sy values to the simulation domain.Typical values are the geometric means of K and Sy.Then the forward simulation will give the first attempt at simulated hydraulic heads for the simulation domain.
2. Collect the simulated heads at locations where the observations (i.e., the hydraulic heads) are available.Calculate the head differences between simulated and observed heads.These head differences are called residual heads in this study.
Incorporate the head differences and the observed K and Sy in a cokriging system into updated ln K and ln Sy fields.Note that the geostatistical parameters for the cokriging system are the initial guess values at this step.
3. Run a forward simulation based on the updated ln K and ln Sy fields.Then the differences between simulated and observed heads are again recalculated.This process can be considered as the first iteration for the first stress event.We treat this flow field as the conditional mean flow field (i.e., the solution of Eq. 5) and it will be used for stochastic simulation in the next step.
4. Use the adjoint state method to solve for sensitivity equations.The sensitivity matrices are then employed to calculate the cokriging covariance matrices for residual heads and the cross covariances for residual head, ln K , and ln Sy.
5. Update the ln K and ln Sy fields by using the cokriging interpolation and the associated covariance and cross-covariance matrices from the previous step.meets the predefined convergence criterion, a different stress event and the associated head observations are selected.Then the estimation process goes to step 3 and repeats steps 3 to 6 until all the stresses are involved in the estimation.Usually the change in cokriging error variances is the rigid bound to stop the iterations for each stress event.

Sensitivity estimations for covariance matrices
To calculate the covariance of heads and cross covariance of heads and ln K and ln Sy requires the determination of the sensitivity matrices.This study employs the adjoint state method to conduct sensitivity analyses (Sykes et al., 1985;Sun, 1994).
The detailed derivation of the adjoint state approach can be found in the study of Sykes et al. (1985) for a steady-state confined aquifer, in the study of Zhu and Yeh (2005) for a transient confined aquifer, and in the study of Ni and Yeh (2008) for gas-phase flows in unsaturated porous media.For a flow in a transient unconfined aquifer system, we derive the state sensitivities based on the governing Eqs.(1) to (4).The mathematical formula of adjoint state equation is the following formula: subject to the boundary and initial conditions In Eqs. ( 6) to ( 9), ψ is an arbitrary function that will be used in estimating the state sensitivity for the observation at location x k or at time t l , where the subscripts k and l represent the indices of locations and times for head observations.Notation δ( ) is the delta function.The solutions to Eqs. ( 6) to (9) require the conditional mean head equation (i.e., Eq. 5) to be solved.Additionally, the number of simulations for Eq. ( 6) depends on the number of spatial and temporal head observations.In this study the K and Sy are considered to be uncorrelated.Under this condition, the sensitivity of the head at location x k and time t l to ln K at location x n is given by where Ω and Γ represent the spatial and temporal domains.One can apply a similar procedure to obtain the sensitivity of the head at location x k and time t l to l nSy at location x n , which is In the adjoint state method Eqs. ( 10) and ( 11) are incorporated with covariances of ln K and ln Sy to obtain the cross-covariance matrices.For each pumping or injection test, evaluations of the sensitivities and statistical structure are the basis for estimating our conditional parameters K and Sy.

Numerical example
To verify the developed model for estimating K (or T ) and Sy in unconfined aquifers, we first tested the developed model with a synthetic example.The scale and values of aquifer parameters for the synthetic example are similar to those of our well field.Introduction

Conclusions References
Tables Figures

Back Close
Full The following sections present how we create the synthetic model and estimate the unconfined aquifer parameters ln K and ln Sy based on the generated synthetic head observations.

Model description
In this numerical example the size of the aquifer is 40 m × 20 m and has a constant aquifer thickness of 10 m. Figure 1 shows the generated ln K and ln Sy fields.In this study the fast Fourier transform algorithm (Deutsch and Journel, 1997) was employed to generate the ln K and ln Sy fields based on an exponential covariance model.The variances of ln K and ln Sy for the exponential covariance model were 1.0 and 0.1, respectively.Isotropic correlation lengths for the ln K and ln Sy fields were fixed at 10 m in the exponential covariance model.Under the predefined conditions for random field generation, the range of generated K values was from 5 to 200 m day −1 .However, the values of Sy varied between 0.05 and 0.16.
Figure 2 shows the conceptual model of the numerical example.The constant head conditions were specified at the north and south boundaries.There were 5 pumping wells and 10 monitoring wells in the simulation area.The uniform finite element mesh for this example was 1 m × 1 m.Such an element size required 800 elements and 1722 nodes to cover the entire simulation area.To generate the synthetic head observations for testing our model, this study conducted forward simulations based on a constant pumping rate applied at the locations of the pumping wells (as marked by circles).
For each pumping event, the transient head observations were then collected from the observation wells (as marked by triangles).The contour lines shown in Fig. 2 represent the solution for one of the pumping events.All these five pumping data sets were then integrated into the inverse model to estimate the ln K and ln Sy fields shown in Fig. 1.Introduction

Conclusions References
Tables Figures

Back Close
Full

Results and discussion of the numerical example
Figure 3 shows the results that are estimated by the developed inverse model.With a total of 25 iterations (5 iterations for each pumping stress), the developed inverse model can capture the patterns and magnitudes of the generated ln K and ln Sy fields (see Fig. 1).To quantitatively evaluate the accuracy of the inverse model, the scatterplots in Fig. 4 show the element-by-element comparison of simulated and true values.
The correlation of the generated and estimated parameters is based on the following formula: where θ and θ * represent the geometric means of parameters.Based on the definition of the correlation value, this value represents the overall normalized difference between simulated and generated ln K or ln Sy fields.However, the mean absolute error normal (L1) and mean square error normal (L2) values reflect the overall relative errors between the estimated and the generated ln K or ln Sy field.Note that the value of ln K variance for generating the ln K field was 1.0, while the variance for generating the ln Sy field was only 0.1.There is an order of magnitude difference between the ln K and ln Sy variances.The small differences of correlation values between ln K and ln Sy estimations imply that the transient head data might be insensitive for the estimation of the ln Sy field.
To test the effect of transient head observations on the ln K estimation, we treated the last data point as the steady-state head value for each pumping stress and reran the inverse model.The result of the ln K estimation showed an identical ln K distribution pattern.The correlation value was slightly higher than the one that was obtained based Introduction

Conclusions References
Tables Figures

Back Close
Full Figure 5 shows the comparison of observed and simulated head data at wells for five pumping events.Note that the observed data were the generated head observations based on the synthetic ln K and ln Sy fields.Most observation heads fall in the range of 7.8 m-8.0 m because the boundary conditions in the synthetic example had been fixed in this range.The values below 7.8 m show the observation data near or at the pumping wells.Figure 6 shows the distribution of error variances for ln K and ln Sy from our inverse model.In the synthetic example, ln K and ln Sy values at the locations of monitoring and pumping wells were fixed.These ln K and ln Sy values are considered to be the conditioning points and are incorporated in the cokriging iterations.The results clearly reflect relatively low error variance values at these locations.Note that the error variance values for ln Sy were generally smaller than those for ln K because we used a relatively small variance to generate the ln Sy field.We found that the specified head conditions along the north and south boundaries can lead to small ln K estimation errors along these two boundaries.Such small error variance values along specified head boundaries significantly control the developments of estimation errors for ln K field inversions.However, in the situation of ln Sy estimations, the relatively small error variances are only obtained near the conditioning points (i.e., at the pumping and monitoring wells).

Field-scale implementation
The previous section has shown that the developed inverse model can estimate synthetic ln K and ln Sy fields with accurate patterns and magnitudes.In most situations, the operation of pumping tests can be replaced by injection tests, depending on the conditions at the site of interest.In this study we developed a well field in a highly per- in the synthetic example.However, due to the limited saturated thickness at the test site and too small well diameters for large-capacity pumps, cross-hole injection tests were conducted to obtain head observations.

Site description
The rivers in Taiwan flow mainly in the east-west direction because of high mountain areas in the central portion of the island.Important characteristics of the river flows in Taiwan are the short concentration times and large river slopes.This generally leads to large grain sizes for river deposits, and the river discharges are strongly variable over time.The study area is located in the middle reach of the Wu River in central Taiwan.
Figure 7 shows the location of our well field at the test site.The unconfined aquifer materials for the test site are mainly alluvium deposits.The Holocene alluvial deposits are mainly composed of gravel, sand, and silt, with thicknesses varying from 7 to 10 m at the test site.Based on site investigations and analyses of soil samples from well logs, these gravel deposits can be categorized as GW-GM (well-graded gravel with sand and silt) according to the Unified Soil Classification System (USCS).The Pliocene Cholan Formation was identified below the unconfined aquifer.The Cholan Formation in this area is composed of interbedded sandstone and shale.Due to the relatively small hydraulic conductivity of the Cholan Formation, the interface of the Cholan Formation and the alluvium deposits is considered to be the bottom boundary for our numerical model.
The average groundwater table in the test site is 2.4 m below ground level.Figure 8 shows the distribution of developed wells at the experimental site.The well names starting with "DP" represent the injection wells, while the names with "DW" indicate the monitoring wells.The injection wells were installed with 4-inch PVC tubes and opened through the unconfined aquifer, while the monitoring well were installed with 2-inch PVC tubes.In this study a total of 15 wells was installed in the well field for our field-scale injection tests.Introduction

Conclusions References
Tables Figures

Back Close
Full

Cross-hole injection test
Our well field is located in the middle reach of the Wu River.The aquifer is highly permeable and the injection test at the research site requires a high-capacity pump (approximately 400 L min −1 ) to produce sufficient head changes for observations.The water used for injection was sourced from the irrigation system near the well field.
A constant flow rate of 400 L min −1 irrigation water was pumped from the irrigation channel and injected into the injection wells during the tests.
Figure 9 shows the time series of selected head observations from our injection tests.We measured the head values both in monitoring and injection wells during the injection tests.With sequentially switching the injection wells there were 14 head variation time series for each injection event.The subfigures on the left column (Fig. 9a, c, and e) show the head observations at injection wells (i.e., the DP wells), and the subfigures in the right column (Fig. 9b, d, and f) show the selected head time series at monitoring wells (i.e., the DW wells).Note that the heads at injection wells were not included in our simulations.During the injection test, we found that the head drawdowns at the injection wells increase dramatically at the beginning (less than 60 s) of the injection.Only the DP9801 well shows a gradual increase of heads at the injection well (Fig. 9a).However, the monitoring wells always show gradually increasing heads.During the injection tests the head observations at injection wells were monitored in real time and the injections were stopped if the head changes at injection wells showed relatively small changes of heads.In a unique case for injection at the DP9801 well, the power generator unexpectedly shut off after 500 s during the injection experiment.The observations for this particular situation were also included in the parameter estimation process.

Parameter estimations
In this study the boundary and initial conditions of the simulation area were defined based on the interpolation and extrapolation of water level measurements in the well Introduction

Conclusions References
Tables Figures

Back Close
Full field (see Fig. 8).The values of constant head boundaries on the day (30 October 2009) that we conducted the injection tests were slightly smaller than those that we used for our synthetic example.In the model we assumed that no K and Sy observation was predefined in the parameter estimation.The sampling time interval for head observations in this study was fixed to one second.Based on the data quantity of our cross-hole experiments, it is computationally difficult to include all the head time series in the inverse model because of the storage and calculation of covariance matrices.Here we follow the suggestions by Ni and Yeh (2009) and Cardiff and Barrash (2011) indicating that the key observations in early, middle, and the late time intervals are sufficient to reproduce the high-resolution aquifer parameters well.From each monitoring well, we then selected two head observations from the early time interval (less than 60 s), two head observations from the middle time interval (between 60 and 200 s), and two observations from the late time interval (between 500 to 1000 s).The artificial selection of head observations was based on our experience of parameter estimations.Different selections of head observations may lead to slightly different parameter estimations.
However, the patterns of the parameter distributions should be similar.
Figure 10 shows the distributions of the estimated ln K and ln Sy for the test site.The hydraulic conductivity (K ) values vary from 10 to 100 (m day −1 ), while the specific yield (Sy) values are in the range of 0.02-0.16.The patterns of the ln K are complex, and the high ln K zones are located along the northern boundary and near the well locations.
The pattern of ln Sy distribution is relatively simple.A high ln Sy zone was found near wells DP9804 and DW9806.It is reasonable to suspect that these high ln K or ln Sy zones may be relevant to the developments of the well field.In this study we are not focusing on the validation of well installation processes at this site.
Figure 11 shows cokriging error variance after 25 iterations (5 iterations for each injection event) for ln K (Fig. 11a) and ln Sy (Fig. 11b).Because there was no conditioning point assigned at the well location, the results can only reflect the estimation error variances based on the head observations.The results of Fig. 11a show high error variance areas in the central portion of the modeling domain.Although the low-ln K Introduction

Conclusions References
Tables Figures

Back Close
Full areas show slightly large error variances, the differences between the error variances are very small (see Figs. 10a and 11a).Similar results hold for the ln Sy estimations shown in Fig. 10b and 11b.The constant head boundary condition plays an important role in controlling the development of the ln K error variances.To evaluate the effect of the constant head boundaries on the estimation results, two additional models with expanded domain sizes were conducted, and the results are compared in detail in the next section.

Boundary effect on parameter estimations
To evaluate the effect of the constant head boundaries on the estimation results, we artificially expanded our simulation domains and applied the same data sets.For this purpose, two additional models were developed.Figure 12 shows the sizes for three different models.The two additional models were named Model-1 for the largest 70 m × 50 m model and Model-2 for the 50 m × 30 m model.The original model with the size of 40 m × 20 m was named Model-3.Because head observations were not available for assigning the constant head boundary conditions, we then used the boundary values that were extrapolated based on head observations from the well field.Detailed information for all the boundary conditions is shown in Fig. 12.We used the same size of the computational mesh for the two large-sized models.Based on the mesh size of 1 m × 1 m, the number of elements is 3500 for Model-1 and 1500 for Model-2.The number of elements in the original model (Model-1) is 800.It is worth mentioning here the computational issue for the developed model and the expanded simulation domains.In general about 80 % of the computational time is used for estimating sensitivity matrices by employing the adjoint state method.Because the head changes are not significant, two to three iterations are generally required to obtain convergent solutions for nonlinear mean head and adjoint state equations (i.e., Eqs. 5 and 6).Based on our workstation with an Intel Xeon 2.0 GHz CPU and 12 GB memory the computational times for three different models are about 30 min for Model-3, 50 min for Model-2, and 80 min for Model-1.This result clearly shows that the size of simulation domain dramatically increases the computational effort required to solve for the mean head and adjoint state equations.However, in this study the improvement of estimated ln K and ln Sy fields for expanded models is not significant.

Conclusions
This study developed a two-dimensional inverse model for transient and unconfined aquifers based on the concept of sequential cokriging interpolation proposed by Zhu and Yeh (2005).Here the governing equations for stochastic simulations of flow in unconfined aquifer systems were reformulated and the associated programs were modified to solve for the nonlinear flow and adjoint state equations.The numerical example in this study showed that the developed model can detect detailed spatial variations of ln K and ln Sy based on the generated aquifer parameters and the associated head observations.With selected head observations in early, middle, and late times from monitoring wells, the sequential inclusion of head measurements from different Figures pumping/injection events can significantly reduce the sizes of covariance matrices for ln K and ln Sy to head measurements.The distribution of estimation error (or cokriging error) variances for the numerical example show significantly different mechanisms in ln K and ln Sy estimations.With the exception of the predefined conditioning points, the ln K estimation error variances are strongly constrained by the constant head boundaries.Unlike the results of ln K estimation error variances, the ln Sy error variances are not affected by the boundaries, implying that the estimations of ln Sy distributions could be source dominated and are highly relevant to transient processes.However, estimation of ln K can be boundary dominated because the constant head boundary conditions are deterministically defined such that the error variances are forced to zero at constant head boundaries.This situation is similar to the concept of forward flow and transport simulations.The distances of the boundaries should be far enough to eliminate the boundary effect if the natural hydrogeological boundaries are not available.Typical approaches to obtain reliable simulation results are to conduct additional models with different simulation domains but using the same aquifer parameters and the associated hydrological stresses.
In this study we have developed a well field for field-scale injection tests.The study area is located in the middle reach of the Wu River in central Taiwan.The unconfined aquifer materials for the test site are mainly alluvium deposits, and such alluvium deposits are composed of gravel, sand, and silt, with thicknesses varying from 7 to 10 m.The Pliocene Cholan Formation was identified under the unconfined aquifer.Due to the relatively small hydraulic conductivity of the Cholan Formation, the interface of the Cholan Formation and the alluvium deposits was considered to be the bottom boundary of the inverse model.A total of 15 wells -5 injection wells and 10 monitoring wells -was installed in the test site.The depth-averaged field injection tests were conducted on the day of 30 October 2009.
Based on the 5 injection tests and the associated head observations from the well field, we implemented the inverse model to estimate the two-dimensional (depthaveraged) ln K and ln Sy distributions for the well field.Simulation results show that Introduction

Conclusions References
Tables Figures

Back Close
Full Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | employed the Dupuit-Forchheimer assumption to analyze pumping test data from fully screened wells in an unconfined aquifer in Boise Hydrogeophysical Research Site (BHRS).The quasi-linear Bayesian geostatistical inverse method developed by Kitanidis (1995) was used to optimize the distribution of the depth-averaged hydraulic conductivity.The potential-based approach allows the unconfined nonlinear flow equation to be transformed to a linear one and be solved by a typical optimization algorithm.Cardiff and Barrash (2011) summarized the investigations of hydraulic tomography surveys (HTS) for two-and threedimensional problems and used numerical experiments to assess the effect of the storage coefficient and the geostatistical parameters on estimating K values in transient Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 6. Check convergence criteria for the working stress event.This study uses the summation of residual heads and the cokriging error variance as the convergence criteria.If the summation of the residual heads or the change in cokriging variance Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | on transient head data.Such a small difference may reflect the effect of the redundant transient head data.
meable aquifer in central Taiwan to test the developed model for field-scale problems.The scale and the boundary conditions for the experiments were similar to those shown Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Figures
Figures 13 and 14 show the results based on the same head observations at well locations.We are especially interested in the areas within the well fields (i.e., the areas surrounded by the wells) because these areas have the head observations available for ln K and ln Sy estimations.High ln K values are located in the northern portions of the well fields.With the changes of modeling domains the patterns of estimated ln K fields near the well fields remain similar.Away from the well fields the ln K values rely on the updated statistical structures.The unexpectedly high and low ln K values outside the well field might be caused by the numerical accuracy.A similar conclusion 14965 Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 3 .
Fig. 3.The estimated (a) ln K and (b) ln Sy fields for the numerical example.

Fig. 4 .Fig. 6 .
Fig. 4. Scatterplots for direct comparison of generated and estimated parameter values.Left: the generated vs. estimated ln K values and the associated error analysis.Right: the generated vs. estimated ln Sy values and the associated error analysis.L1 and L2 represent the mean absolute error normal and mean square error normal, respectively.

Fig. 11 .
Fig. 11.The distributions of the cokriging error variance for (a) ln K and (b) ln Sy after 25 iterations (5 iterations for each injection).