Bayesian approach for three-dimensional aquifer characterization at the Hanford 300 Area

This study presents a stochastic, threedimensional characterization of a heterogeneous hydraulic conductivity field within the Hanford 300 Area, Washington, USA, by assimilating large-scale, constant-rate injection test data with small-scale, three-dimensional electromagnetic borehole flowmeter (EBF) measurement data. We first inverted the injection test data to estimate the transmissivity field, using zeroth-order temporal moments of pressure buildup curves. We applied a newly developed Bayesian geostatistical inversion framework, the method of anchored distributions (MAD), to obtain a joint posterior distribution of geostatistical parameters and local log-transmissivities at multiple locations. The unique aspects of MAD that make it suitable for this purpose are its ability to integrate multiscale, multi-type data within a Bayesian framework and to compute a nonparametric posterior distribution. After we combined the distribution of transmissivities with depthdiscrete relative-conductivity profile from the EBF data, we inferred the three-dimensional geostatistical parameters of the log-conductivity field, using the Bayesian model-based geostatistics. Such consistent use of the Bayesian approach throughout the procedure enabled us to systematically incorporate data uncertainty into the final posterior distribution. The method was tested in a synthetic study and validated using the actual data that was not part of the estimation. Results showed broader and skewed posterior distributions of geostatistical parameters except for the mean, which suggests the importance of inferring the entire distribution to quantify the parameter uncertainty. Correspondence to: Y. Rubin (rubin@ce.berkeley.edu)


Introduction
Hydrogeological characterization plays a key role in various projects involving groundwater flow and contaminant transport.A detailed three-dimensional (3-D) description of spatial variability in subsurface hydraulic properties is imperative for predicting water and solute movement in the subsurface (Rubin, 2003).Recent focus on geochemical and microbiological reactions in field studies, for example, requires flow parameters to be fully characterized a priori for testing their research hypotheses (Scheibe et al., 2001;Scheibe and Chien, 2003;Fienen et al., 2004).
One of the main challenges in hydrogeological characterization is to integrate datasets of different types and scales.Typical field studies usually include two or more different complementary sources of information, which may include depth-discrete small-scale data such as core analysis, slug tests and electromagnetic borehole flowmeter (EBF) tests and large-scale data such as pumping tests and tracer tests.With stochastic modeling of flow and transport becoming increasingly common, it is important not only to combine best-fitted values from each dataset, but also to correctly quantify and weigh errors and uncertainty associated with different datasets, and to transfer the uncertainty to the final prediction (Maxwell et al., 1999;Hou and Rubin, 2005;De Barros et al., 2009).
To tackle this challenge, various researchers have applied Bayesian approaches to the problem of subsurface characterization (Copty et al., 1993;Woodbury and Rubin, 2000;Chen et al., 2001).Within a Bayesian framework, the probability density function of a parameter vector can be updated sequentially to include more datasets in a consistent manner.In addition, the resulting predictive distribution can H.Murakami et al.: Bayesian approach for three-dimensional aquifer characterization properly account for the parameter uncertainty inherent in estimating the parameter values from the data (Diggle and Ribeiro, 2002).Two recent developments in particular have increased the potential of the Bayesian approach for subsurface characterization: (1) Bayesian model-based geostatistics, and (2) the method of anchored distributions (MAD).
Bayesian model-based geostatistics, introduced by Kitanidis (1986) and Handcock and Stein (1993), assumes a parametric model for a spatial stochastic process and infers geostatistical structural parameters based on small-scale datasets or point measurements (Diggle and Ribeiro, 2006).While the traditional variogram approach determines bestfitted estimates of geostatistical structural parameters and their asymptotic confidence interval, the Bayesian modelbased approach yields a posterior distribution of the parameters.Diggle and Ribeiro (2006) showed that correlation parameters such as variance and scale follow non-Gaussian and skewed distributions, which suggests that the first two moments are not enough to characterize the distribution.
The method of anchored distributions (MAD) is a general Bayesian method for inverse modeling of spatial random fields that addresses complex patterns of spatial variability, multiple sources and scales of data available for characterizing the fields, and the complex relationships between observed and target variables (Zhang and Rubin, 2008a, b;Rubin et al., 2010).The central element of MAD is a new concept called "anchors".Anchors are devices for localizing large-scale data: they are used to convert large-scale, indirect data into local distributions of the target variables.The goal of the inversion is to determine the joint distribution of the anchors and structural parameters, conditioned on all of the measurements.The structural parameters describe large-scale trends of the target variable fields, whereas the anchors capture local heterogeneities.Following the inversion, the joint distribution of anchors and structural parameters can be directly used to generate random fields of the target variable(s).Different from most of the other inversion methods that determine a single best estimate of the field and asymptotic uncertainty bounds (Kitanidis, 1995;Zhu and Yeh, 2006;Ramarao et al., 1995), MAD yields a posterior distribution of the parameters.
In this paper, we assimilate EBF tests and constant-rate injection tests for characterizing a 3-D hydraulic conductivity K field at the Integrated Field Research Challenge (IFRC) site in the Hanford 300 Area (http://ifchanford.pnl.gov).Since the EBF tests yield only relative K values along each of the EBF test wells, we need a local transmissivity T value at each of the EBF test wells to convert the relative values to absolute K values (Javandel and Witherspoon, 1969;Molz et al., 1994;Young et al., 1998;Fienen et al., 2004).The local T values can be determined by inverting the largescale constant-rate injection tests.This assimilation requires us to quantify the uncertainty in T values based on the injection tests and to combine that uncertainty with the one from the EBF data.The particular difficulty in inverting injection-test or pumping-test data is the computational effort associated with a long time series.Li et al. (2005) and Zhu and Yeh (2006) have applied temporal moments of drawdowns to estimate T and storativity S fields.The sandbox experiment by Liu et al. (2007) has shown that the moment approach can successfully characterize a T field.The advantage of employing temporal moments is that we can compute them using steady-state flow equations, which can reduce the computational burden significantly.In addition, when the interest is limited to T , the zeroth-order temporal moment can eliminate the effects of the uncertainty in S, since it does not depend on S (Zhu and Yeh, 2006).
In the following sections, we first describe the site and the experimental procedure.We then present our approach, including the geostatistical inversion framework and the inference of the 3-D geostatistical parameters.After presenting the inversion results in a synthetic study to demonstrate and verify our approach, we discuss the results using the actual data at the site.

Site and experiment description
The Hanford 300 Area is located at the southern part of the Hanford Nuclear Reservation one mile north of Richland, Washington, USA.The IFRC site is located within the footprint of a former disposal facility for uranium-bearing liquid wastes known as the South Process Pond, approximately 250 m west from the Columbia River.As is shown in Fig. 1, the triangular well field consists of 25 wells fully screened Hydrol.Earth Syst.Sci., 14, 1989Sci., 14, -2001Sci., 14, , 2010 www.hydrol-earth-syst-sci.net/14/1989/2010/ through the saturated portion of the Hanford formation, ten wells partially screened at different depths, and one deep characterization well (Bjornstad et al., 2009).
In this study, we focus on the saturated portion of the highly permeable and coarse-grained Hanford formation, which is a shallow unconfined aquifer.The main lithology is a poorly sorted mixture, dominated by gravel up to boulder size, with lesser amounts of sand and silt (Bjornstad et al., 2009).It overlies the Ringold formation, the upper portion of which is a continuous low-permeability layer consisting of cohesive and compact, well-sorted fine sand to silty sand.The saturated thickness is variable over the site, ranging from about 5 m to 8 m, with daily and seasonal fluctuations of the water table in response to changes in the river stage.The prior estimates of hydraulic conductivity are 1000-100 000 m/day for the Hanford formation and 0.01-3.00m/day for the Ringold formation (Meyer et al., 2007).
The aquifer in the Hanford formation has been very difficult to characterize, since typical methods, such as permeameter tests and slug tests, do not provide reliable results due to the coarse-grained and highly permeable nature of the aquifer (Meyer et al., 2007;Newcomer, 2008;Vermeul et al., 2009).In addition, traditional analysis of pumping and injection tests yields only the averaged property over a large domain, since the zone-of-influence expands very rapidly.Combining EBF tests and injection tests through inverse modeling is one of the few feasible alternatives to characterize the 3-D heterogeneous structure of the aquifer.
Fourteen constant-rate injection tests were conducted, each of which had one injection well and 7 to 10 observation wells.All the wells used in the tests are fully screened over the saturated portion of the Hanford formation.The distance between the injection and observation wells ranges between 8 and 60 m.The injection duration and rate are approximately 20 min and 315-318 gpm (1.19-1.20 m 3 min −1 ), respectively.The preliminary analysis of the late-time curve data, using the Cooper-Jacob straightline method, has shown that most of the observation wells yield similar estimates for the T values in each test, which is considered to be the geometric mean of T , T G , over the entire well field, as is mathematically proved by Sánchez-Villa et al. (1999).It suggests that the zone-of-influence expands very rapidly and the conventional pumping test analysis yields only an effective property, smoothing out the local heterogeneity at the well field.
The EBF test data were obtained at 19 fully screened wells, which yielded 283 depth-discrete relative hydraulic conductivities with 0.3-0.6 m depth intervals.The pumping rate was 1.04-1.55gpm (3.94×10 −3 -3.94×10 −3 m 3 min −1 ), and kept constant during the test at each well.The vertical profiles indicated that the hydraulic conductivity over the central third of the Hanford formation was lower than its top and bottom thirds at many of the wells.Although the thickness and contact depths for this lower permeability material vary across the site, this general pattern was observed to some extent at most of the monitoring well locations.
The more detailed description of the site and data is available in Bjornstad et al. (2009) and Rockhold et al. (2010).

Methodology
For the 3-D characterization, we employ a two-step approach to combine the EBF and injection-test data.First, we invert the large-scale injection tests to characterize the 2-D T field.We apply MAD to invert the zeroth-order moments of pressure build-up curves at multiple observation wells.As a result, we obtain a joint posterior distribution of T at the EBF test wells.Second, we combine this distribution with the EBF data for determining the absolute K values.Instead of a single K value at each of the EBF data point, we obtain the distribution of K at each point.Based on the distribution of the absolute K, we infer the 3-D geostatistical parameters, using the Bayesian model-based geostatistics.
Compared to direct coupling of the EBF and pumping tests used in Li et al. (2008), this two-step approach has a significant computational advantage.This approach is possible because we can model the flow process during the injection tests as 2-D planar flow in the horizontal plane, due to the particular site conditions as the following.The coarsegrained and highly permeable nature of the aquifer caused the elastic response and drainage effect to occur very rapidly (less than 30 s after the injection started), so that the radial flow regime dominated the pressure buildup responses (Neuman, 1975).In addition, despite the large injection rate, the maximum pressure buildup at the nearest observation wells was less than several centimeters, which is much smaller than the aquifer thickness (5-8 m).Although the EBF tests suggested vertical heterogeneity in the saturated zone, Dagan et al. (2009) recently showed that Dupuit's assumption is still valid -when the aquifer thickness is not large compared to the vertical integral scale, and the ratio between the vertical and horizontal integral scale is large, which is the case at this site.

MAD framework
In this section, we summarize the Bayesian inversion framework, called (as indicated above) method of anchored distributions (MAD), which we use to invert the injection test data.This method was introduced by Zhang and Rubin (2008a, b) and Rubin et al. (2010) and is summarized here for completeness.
We denote a spatial random process by Y (x), where x is the space coordinate.We further denote an entire field of Y by a random vector Ỹ , and denote a realization of the field by ỹ.The dimension of Ỹ and ỹ is equal to the number of elements in the discretized field.The field Ỹ is defined through the vector of model parameters {θ,ϑ}.The θ part of this vector, called the structural parameter vector, includes a set of parameters designed to capture the global features of Ỹ , such as the mean of the field and correlation structures.The ϑ component of this vector consists of the anchored distributions, or anchors in short.Anchors are devices used to capture local features of Ỹ that cannot be captured by the global parameters θ.In their simplest form, an anchor would be error free measurements of Y .Other forms of anchors include measurements coupled with error distributions and/or anchors that are obtained by inversion.
The data z is a vector of multiple observations of a physical process.The data can be described by the following equation: where M is a known function, or a set of functions, numerical or analytical, of the spatial field, and ε is a vector of zero-mean errors.The goal of the inversion is, first, to derive a posterior distribution of the model parameters conditioned on the data z, p(θ,ϑ|z).This distribution then allows us to generate multiple realizations of the field Ỹ for prediction.Using Bayes' rule, we can define the posterior distribution of parameters as: ( where p(θ) is the prior distribution, p(ϑ|θ) is the anchor distribution given a structural parameter vector θ, and p(z|θ,ϑ) is the likelihood of data z given a parameter set {θ,ϑ}.
We estimate the likelihood p(z|θ,ϑ) using the Monte Carlo simulations.Since the model parameters {θ,ϑ} and the data z are connected through the field, we generate multiple conditional realizations of the field Ỹ for any given {θ,ϑ}; with each realization, the forward model provides a prediction of z in the form of z, according to Eq. (1).In other words, z is viewed as a measured outcome from a random process, whereas z is one of many possible realizations, given a particular parameter set of {θ,ϑ}.By generating random fields for a given parameter set {θ,ϑ} and simulating the physical process on each field, we obtain multiple realizations of z, i.e., an ensemble of z.The ensemble of z is then used for estimating the probability density function (pdf) of z.After determining the pdf, it is straightforward to calculate the density at a point z, p(z|θ,ϑ), which is the likelihood.
As the dimension of z increases, a larger number of realizations of z are necessary to estimate the pdf accurately, which increases the computational burden.To accommodate the large-dimensional data, we may divide the vector z into L segments as z = {z 1 ,z 2 ,...,z L }.We can decompose the likelihood into each segment as, In Eq. ( 3), we assume that the data segments z 1 , z 2 ,. . ., z L are conditionally independent for a given {θ,ϑ}, since we consider that {θ,ϑ} contains information equivalent to the data.This equality strictly holds when the data segments are independent of each other -for example, when the data locations are beyond the zone-of-influence or zone-ofcorrelation.It approximately holds when the data segments are only weakly correlated, such as with different types of data at the same site.As Hou and Rubin (2005) pointed out, assuming independence leads to higher entropy and makes the estimation less informative.

Specification of a 2-D geostatistical model
Here we specify the geostatistical model for the 2-D T field.Let Y (x) be natural-log transmissivity, lnT (x), at the location x = (x 1 , x 2 ) in the 2-D domain.We assume that a vector Y , containing Y at multiple locations x, follows a multivariate Gaussian distribution with exponential covariance.We define a structural parameter vector as θ = {µ,σ 2 ,φ}, including uniform mean µ, variance σ 2 , and integral scale φ, which are used at a geologically similar site (i.e.unconsolidated glacial materials) (Rubin, 2003;McLaughlin et al., 1993).We define a vector ϑ(x ϑ ) to represent a set of anchors.Since the anchors are a subset of the field, p(ϑ|θ) is a multivariate Gaussian distribution with mean µ and covariance σ 2 R (2-D) (x ϑ ,x ϑ ), where R (2-D) (x ϑ ,x ϑ ) is an autocorrelation matrix as a function of φ and the locations of ϑ, x ϑ .The distribution of Y conditioned on the structural parameters and anchors p(y|θ,ϑ) is a multivariate Gaussian distribution with conditional mean µ Y |ϑ and conditional covariance Y |ϑ , where the mean and covariance conditioned on the anchors are defined as where R (2-D) (x,x) is the auto-correlation matrix for Y , and R (2-D) (x,x ϑ ) is the cross-correlation matrix between Y and ϑ.

Specification of likelihood
We consider the data z consisting of L injection tests (l = 1, 2, . . ., L).We divide z into L segments as z = {z 1 ,z 2 ,...,z L }, where z l is the vector containing the zeroth-order moments at multiple observation wells in the l-th injection test.The governing equation and the temporal moment formulation are shown in Appendix A.
In order to determine the likelihood p(z|θ,ϑ), we first compute the likelihood in each injection test p(z l |θ,ϑ).
Using the ensemble of zl simulated on the multiple fields conditioned on each parameter set {θ,ϑ}, we calculate the mean and covariance to determine p(z l |θ,ϑ) (Robert and Casella, 2005).When we include the multiple injection tests in the inversion, we multiply the likelihoods of multiple tests, according to Eq. ( 3), to obtain the likelihood for the entire data p(z|θ,ϑ).We have observed that the zeroth-order moments from the different injection tests are not strongly correlated, so that we may use Eq.(3).

Placement of anchors
The success of MAD depends on placement of the anchors from two reasons.First, careful placement of anchors will maximize their ability to extract information from observations.Second, the computational burden is linked to the number of anchors, and a smaller number would improve computational efficiency.Hence, we need to place anchors (1) at sensitive locations to the data, (2) to capture local features of the field, and (3) according to the goal of the inversion.A detailed discussion is available in Rubin et al. (2010).Here we discuss the issues relevant to our inversion.
First, to find sensitive locations, we refer to the sensitivity analysis.Li et al. (2005) has formulated the sensitivity of zeroth-order moments to a lnT value at a specific location, using the adjoint-state method (Sun, 1994).For the current estimate of T field required in the sensitivity analysis, we may use a uniform T field, since our priors are only for the global parameters (i.e., mean, variance, and scale) and we do not have any local information before the injection test.In this case, sensitivity is high around observation well locations, which is consistent with findings by Castagna and Bellin (2009) and Vasco et al. (2000).
Second, to capture heterogeneity, we would ideally have more than one anchor per integral scale.Although the real integral scale is not known in advance, we may consider the minimum possible integral scale at the site.Anchors outside the well plot, far from any of the observation wells, are not effective in resolving spatial heterogeneity, so that we need fewer anchors outside the well plot.
Third, to achieve our goal, which is to obtain the lnT values at the EBF well locations, we need anchors at all EBF well locations.All the EBF wells are used as observation wells during the injection tests, so that we do not need additional anchors for this purpose.

3-D geostatistical model for hydraulic conductivity field
The 2-D inversion of the injection tests yielded a joint distribution of lnT at the EBF well locations.Since we placed anchors at all those locations, we can use the anchor distribution directly.Let us denote the lnT values at the EBF well location by a vector ϑ EBF , which is a subset of ϑ.Marginalizing the other parameters leads to the posterior distribution of ϑ EBF conditioned on the injection test data z as p(ϑ EBF |z).
Let K(x 1 , x 2 , x 3 ) and k(x 1 , x 2 , x 3 ) be the absolute and relative K values at the location x = (x 1 , x 2 , x 3 ) in the 3-D domain, respectively.Based on Javandel and Witherspoon (1969), we have the correlation between the absolute and relative K values as et al., 1994;Fienen et al., 2004), where b(x 1 , x 2 ) is the aquifer thickness at the horizontal location (x 1 , x 2 ).We can then determine the natural log-conductivity u = lnK at (x 1 , x 2 , x 3 ) by We use a N-vector k containing all the relative conductivity values from the EBF data at N locations x, and a Nvector u containing all the lnK values at the same locations as k.Equation ( 5) allows us to combine k and p(ϑ EBF |z) into p(u|k,z), which is the distribution of u conditioned on both the injection test data and the EBF data.
We construct a 3-D geostatistical model, assuming that u(x) follows a multivariate Gaussian distribution with mean β and covariance (η 2 R (3-D) (x,x)+ ν 2 I), where η 2 is the variance of variability in lnK, R (3-D) (x, x) is the auto-correlation matrix for u(x) as a function of the locations x, the horizontal integral scale λ h and the vertical integral scale λ v , I is the identity matrix of order N, and ν 2 is the nugget, which represents the EBF measurement errors.The structural parameter vector of the 3-D geostatistical model is {β,η 2 ,λ h ,λ v ,ν 2 }.Our goal here is to obtain a joint posterior distribution of the parameters conditioned on both data {k,z} through u: For the prior distribution of the parameters, we assume the Jeffreys prior for the mean and variance (Jeffreys, 1946), which is the least informative prior for those two parameters.The prior distribution of all the structural parameters is defined as For the rest of the prior distribution π(λ h ,λ v ,ν 2 ), we use an independent uniform distribution for each of {λ h ,λ v ,ν 2 } bounded by each set of the minimum and maximum possible values.Following Diggle and Ribeiro (Chapter 6, 2006), we obtain an analytical expression for p(β,η 2 ,λ h ,λ v ,ν 2 |u) (Appendix B).In addition, based on the joint distribution of the structural parameters and u, determined by p(β,η 2 ,λ h ,λ v ,ν 2 ,u|k,z) = p(β,η 2 ,λ h ,λ v ,ν 2 |u)p(u|k,z), we can sample multiple parameter sets {β,η 2 ,λ h ,λ v ,ν 2 ,u}, and generate multiple 3-D random fields conditioned on each of the parameter sets.

Organization of constant-rate injection test data
To demonstrate our approach, we used seven out of the fourteen constant-rate injection tests at the site for the synthetic study, and four for the real data analysis (injection at Well 2-09, 2-18, 2-24, and 3-24).The locations of injection and observation wells are well balanced within the IFRC site.Figure 2 shows the configuration of the injection and observation wells for each of the seven tests.
For each test, we calculated the zeroth-order moments at multiple observation wells by integrating the pressure buildup curves.Since the well field is located near the Columbia River, the water table changes according to the river stage fluctuation.Since the change was mostly linear within 20 min after the injection started, we removed the ambient head contribution by linear interpolation.For quantifying the measurement errors, we followed Nowak et al. (2006) and Li et al. (2008), who determined the errors based on fluctuation or noise in the pressure measurements.Since the noise in our case is contained within the range of the instrument resolution, we determined the standard deviation of measurement error based on the resolution of the instrument, 0.003 ft (0.09 cm) by integrating it over the injection duration.

Prior distribution for MAD inversion
For the prior distribution of the 2-D structural parameters θ, we used three independent uniform distributions bounded by the minimum and maximum values, as are shown in Table 1.The prior distributions of each parameter sufficiently cover possible values from the historical data at the site (Meyer et al., 2007) or literature values for similar geological formations (Rubin, 2003).The uniform distributions are considered to be less informative than Gaussian distributions, which have been commonly used in the Bayesian geostatistical inversion (Li et al., 2005).Three thousand sets of θ are generated from p(θ) using a quasi Monte-Carlo sampling method (Krommer and Ueberhuber, 1998).
As is shown in Fig. 3, we placed 44 anchors at all the well locations inside the well plot and at sparse locations outside the well plot, following the discussion in Sect.3.1.4.For each set of θ, we generated 12 sets of anchors ϑ from p(ϑ|θ), so that the number of prior parameter sets {θ,ϑ} is 36 000.We gradually added more parameter sets until we observed the posterior distribution converged -not changing along with increasing numbers of sets.

Forward simulation in MAD
Figure 3 shows the 2-D computational domain used for the forward simulations.For simulating the zeroth-order temporal moments on multiple random fields, we followed the approach by Firmani et al. (2006), since the mathematical expression for the zeroth-order moments is the same as the one for steady-state flow toward a well.
The computational grid size is uniform equal to 0.2φ min in both x 1 and x 2 directions during the field generation.During the flow simulations, the grid blocks near the injection well are divided into non-uniform grid from 0.04φ min at the injection well location to 0.2φ min at a distance of 0.8φ min , satisfying the condition that the ratio between the neighboring block size should not exceed 1.5 (Firmani et al., 2006).We determined the domain size such that the observation wells were 2φ max away from the boundaries to reduce the boundary effect.During the flow simulation, we added another buffer zone with width φ max and uniform T equal to T G between the field and the boundaries for further reducing the boundary effect.We intended to satisfy 3φ max between any observation wells and the boundaries, based on the theory developed by Rubin and Dagan (1988).We used the SGSIM program from GSLIB (Deutsch and Journel, 1998) to generate spatially correlated Gaussian random fields conditioned on eachf set of {θ,ϑ}.We then simulated zeroth-order moments on each field, using a finiteelement method with linear elements.We used 250 realizations of random fields and moments for each {θ,ϑ} to estimate the likelihood p(z|θ, ϑ).We gradually added more realizations, and found that the likelihood values converged with 250 realizations.The 9 000 000 forward simulations took 60 000 computational hours in total.It was divided into several batches, and used up to 9000 cores on the Franklin supercomputer at the National Energy Research Scientific Computing Center (Berkeley, USA), each core of which is a 2.3 GHz single AMD Opteron processor.

3-D geostatistical model
After calculating p(u|k,z) from the relative K values and lnT values at the EBF well locations, we generated a thousand sets of u from p(u|k,z).For each set of u, we computed a posterior distribution p(β,η 2 ,λ h ,λ v ,ν 2 |u), based on the uniform prior distribution of λ h , λ v and ν 2 bounded by the values shown in Table 2.We then integrate the distribution numerically to determine p(β,η 2 ,λ h ,λ v ,ν 2 |k,z).

Results and discussion
We first tested the MAD and numerical setting in a synthetic study for inverting the injection test data, and then we applied it to the actual data at the Hanford site.

Synthetic study for 2-D transmissivity field
We generated a synthetic reference 2-D lnT field with a 2-D structural parameter set {µ,σ 2 ,φ} = {−1.8,1.5,20},shown in Fig. 4. We obtained maximum likelihood estimates of the true parameters as {−1.76,1.46,20.0},with near-exhaustive sampling (one out of every five points) (GeoR package by Ribeiro and Diggle, 2001).We then calculated the zerothorder moments on the reference field and superposed a zero-mean independent Gaussian measurement error, which has the same variance as the actual data from the study site.
Our inversion process is based on the same sets of injection and observation wells as the actual experiments conducted at the IFRC site (Fig. 2).We also combined the different www.hydrol-earth-syst-sci.net/14/1989/2010/ Hydrol.Earth Syst.Sci., 14, 1989-2001, 2010  Figure 5 shows the marginal posterior distributions of the 2-D geostatistical structural parameters {µ,σ 2 ,φ} based on the various number of tests, with their corresponding true values.While the mean has a symmetric Gaussian-like distribution, the variance and scale has broad and skewed distributions.The results are improved with increasing number of tests up to three tests, i.e. the posterior distributions become narrower and biases are reduced.The improvement by ad-ditional tests is more significant for the variance and scale, which suggests that the estimation of variance and scale requires more observations.The improvement, however, is not significant for more than three tests, and the distributions based on four to seven tests are very close to each other, which would suggest that the effect of increasing number of tests could be saturated due to the measurement errors and redundancy of information in the data.Although we may expect tighter distributions with more information, some of the uncertainty cannot be eliminated due to measurement errors and the limited number of observation wells.In addition, the same observation wells were used repeatedly for several tests.Yeh and Li (2000) and Zhu and Yeh (2005) also reported that increasing the number of pumping tests did not improve the estimation above a certain number (three to four tests in their cases).
To examine the effect of anchors and evaluate the random fields generated based on the posterior distributions, we generated 200 000 fields from the posterior distribution of parameters (5000 posterior sample parameter sets with 40 fields per parameter set), and compared the ensemble with the true field.Two cases are studied: one based on a single test (injection at Well 2-18) and the other based on three tests (injection at Well 2-09, 2-24, 3-24).
Figure 6 shows the mean and 98% confidence interval of the predicted lnT fields along the centerline of the well field as shown in Fig. 4 (Line A-B).The centerline also corresponds to the line passing through large variability, from high near the top to low in the middle of the well plot.The figure also includes five realizations that depict the level of heterogeneity in the randomly-generated fields to be used for stochastic simulations.The mean field and random fields along the line all follows a general trend of the true field, especially so with more tests assimilated.If there were no anchors, the mean field would be a flat line at the global mean, and the random fields would be distributed around the flat line.Therefore the deviation from such a straight line is attributed to anchors that capture local heterogeneity.The uncertainty bounds are found to be tighter near the center, where there are more observation wells and more anchors.Increasing the number of tests not only reduces the uncertainty, but also reduces the bias by moving the mean field closer to the true field.

IFRC data analysis
After we gained confidence from the synthetic study, we applied the same scheme to the data from the IFRC site.Since the synthetic study indicated that four tests would be enough, we used up to four tests (the same sets): one injection test (injection at Well 2-18), two tests (injection at Well 2-09 and 2-24), three tests (injection at Well 2-09, 2-24, and 3-24) and four tests (injection at Well 2-09, 2-18, 2-24, and 3-24).Since the true values are unknown in this case, we validated the posterior distribution by comparing predictions with the testing set data, i.e., the observations not included in the inversion.It is a common procedure in statistics to divide the dataset into a training set (i.e., data used in inversion) and testing set (i.e., data used for testing or validating the inversion result).
Figure 7 shows the marginal posterior distributions of the three 2-D structural parameters for the 2-D lnT field at the IFRC site.These three plots show a similar feature to the synthetic study in Fig. 5 such as a Gaussian-like distribution for the mean, broad and skewed distributions for the variance and scale, and the effect of increasing the number of injection tests in the inversion.For the first validation of the 2-D lnT inversion, we generated 200 000 fields (5000 posterior sample parameter sets with 40 fields per parameter set) based on one test (injection at Well 2-09), two tests (injection at Well 2-09 and 3-24) and three tests (injection at Well 2-09, 2-24 and 3-24).We then predicted the zeroth-order temporal moments in the injection test at Well 2-18, which was not the part of the estimation.Figure 8 shows the marginal predictive distributions of the observed moments at two observation wells, based on one, two and three tests, compared with the actual data.The true value is contained within the range of values defining the predictive distributions, and the combination of the three tests improves the prediction by narrowing the distributions.
As another validation, we obtained the maximum a posteriori (MAP) estimate of the geometric mean of T , T G = exp(µ), in Fig. 7, which is 0.52 m 2 s −1 .We compared this value with T G estimated from the Cooper-Jacob analysis (Sánchez-Villa et al., 1999), in which we fitted the late-time pressure build-up curves at multiple observation wells in the injection test at Well 2-18.The 95% confidence bound of T G was 0.38-0.57m 2 s −1 .As we expected, our estimate of T G corresponded to the estimates based on conventional analysis for a large-scale injection test.
Figure 9 shows the marginal distribution for the 3-D geostatistical structural parameters conditioned on the EBF data and injection test data.For the horizontal scale, vertical scale, and nugget, the upper and lower bounds of the xaxis correspond to the bounds of the prior distributions.We can see that the marginal posterior distributions of the structural parameters are skewed except for the mean, which suggests that the entire distribution is necessary to quantify the parameter uncertainty.The horizontal scale has a particularly broad distribution, which is not zero at the bounds.This is because available data is insufficient for narrowing down the distribution (Diggle and Ribeiro, 2002;Hou and Rubin, 2005).Inferring the scale parameters requires different lags (i.e.distances) among the data points.Although we have many different lags for the vertical scale along the boreholes, spacing between the wells restricts variation in horizontal lags for the horizontal scale.It suggests the importance of setting reasonable bounds for the prior distribution based on the information from geologically similar sites.
We also compared the distributions based on the different number of tests included in the injection test inversion.As it turned out, increasing the number of injection tests did not reduce the parameter uncertainty in the 3-D structural parameters as significantly as in the 2-D parameters.This is because the uncertainty and sparseness of the EBF data obscures additional information of the increasing number of injection tests in the 3-D spatial inference.These findings are consistent with Li et al. (2008), who also found that the 3-D characterization of the aquifer was limited by the EBF data density.
Figure 10 shows the 3-D mean field, based on the 5000 parameter sets generated from the distribution p(β,η 2 ,λ h ,λ v ,ν 2 ,u|k,z).We can see the high-low-high layers in lnK along the centerline of the field, which is consistent with the observations in the tracer tests later conducted at the site (Rockhold et al., 2010;Zachara, 2010).

Summary
In this paper, we presented a Bayesian approach for characterizing a 3-D K field by assimilating the EBF and constantrate injection tests.We employed a two-step approach -first inverting the constant-rate injection test data for obtaining the joint distribution of local T values at the EBF well locations, and then converting the EBF data to local K values for the 3-D characterization.For the injection test inversion, we used MAD, which is a newly developed Bayesian geostatistical inversion framework.We inverted zeroth-order moments of pressure buildups at multiple observation wells, which can eliminate uncertainty in a storage coefficient, as well as significantly reduce computational cost.
In a synthetic study, we first showed that MAD could successfully infer the geostatistical parameters and predict the 2-D lnT field.As we included more tests, we could further reduce the uncertainty, and better capture the local heterogeneity.The improvement, however, was saturated at three to four tests due to the measurement errors and redundancy of information in the data, which is consistent with the previous studies (Yeh and Li, 2000;Zhu and Yeh, 2005).
By applying the method to the actual data, we obtained the posterior distribution of geostatistical structural parameters and the anchor values of the 2-D lnT field for the Hanford 300 Area IFRC site.We validated the result using the predictive distribution of the zeroth-moments in the injection test that were not part of the inversion.In addition, the MAP estimate of the mean lnT coincided with the T G value obtained from the Cooper-Jacob analysis, which confirmed our method's consistency with conventional pumping test analysis.
We then combined the relative K values from the EBF data with the distribution of lnT at EBF wells so that we obtain the distribution of the depth-discrete absolute K in the 3-D domain.The uncertainty in T (2-D) is consistently carried on into the lnK values (3-D) as a probability distribution.We thus constructed a 3-D geostatistical model for the lnK field using the model-based Bayesian geostatistical approach.
We demonstrated the advantages of MAD such that MAD was directly connected to the stochastic forward simulations, and it directly inferred the joint distribution of the parameters to be used as an input of the simulations.Compared to the other inverse modeling methods that yield only bestestimates and confidence bounds, MAD fully quantifies the parametric uncertainty as statistical distributions, which is necessary to capture skewed distributions found for the variance and scale parameters.
For the field application, we showed that combining EBF and injection tests is promising for characterizing a heterogeneous lnK field in a coarse-grained and highly permeable aquifer, where the other conventional techniques fail to provide information of local heterogeneity.We found that the 3-D characterization is restricted by the EBF measurement density, in the sense that increasing the amount of depthaveraged information from the injection tests did not contribute significantly to narrowing down the posterior distribution of the 3-D geostatistical parameters.

Fig. 1 .
Fig. 1.Site map of the IFRC site (The coordinate system follows the convention used at the Hanford site).

Fig. 2 .
Fig. 2. Configuration of injection and observation wells in each test used in this paper.The reference point of local coordinates is at (594 164 m, 115 976 m) in the Hanford coordinates.

Fig. 3 .
Fig. 3. Anchor locations in the domain for the constant-rate injection test inversion.The reference point of local coordinates is at (594 164 m, 115 976 m) in the Hanford coordinates.

Fig. 4 .
Fig. 4. Reference field for the synthetic study.Line A-B is used for the transect.

Fig. 5 .
Fig. 5. Marginal posterior distributions of the structural parameters (mean, variance and scale) in the synthetic study, with their corresponding true values.The ones based on the different number of tests are compared.

Fig. 6 .
Fig. 6.Comparison among the reference field, the mean field and the 98% confidence interval of the generated fields, along the center line of the IFRC well field (Line A-B in Fig. 4), for the inversion based on (a) one test (injection at Well 2-18) and (b) three tests (injection at Wells 2-09, 2-24 and 3-24).The gray lines are the realizations of the fields.

Fig. 7 .
Fig. 7. Marginal posterior distributions of the structural parameters (mean, variance and scale) for the Hanford IFRC site data.

Fig. 8 .
Fig. 8.Comparison between the zeroth-order moments observed at Well 2-09 and 3-24 in the injection test at Well 2-18, and predictive posterior distributions from the inversion, including the different numbers of injection tests.

Fig. 9 .
Fig. 9. Marginal posterior distributions of 3-D geostatistical structural parameters of lnK values at the Hanford IFRC site, based on the different number of injection test.

Fig. 10 .
Fig. 10.3-D mean lnK field in the saturated portion of the Hanford formation.The black dots and lines represent the well locations.The reference point of local coordinates is at (594 164 m, 115 976 m) in the Hanford coordinates.

Table 1 .
The lower and upper bounds of the prior distribution for the structural parameters of the 2-D transmissivity field.Table2.1 and Table 2.2 * The upper bound and lower bounds of K multiplied by the average aquifer thickness 7.62 m.

Table 2 .
The lower and upper bounds of the prior distribution for the horizontal scale, vertical scale and nugget variance of the 3-D hydraulic conductivity field.