Stochastic Inversion of a Tomographic Pumping Test: Identifying Conductivity Horizontal Correlation and Longitudinal Macrodispersivity

Two important properties for modeling flow and transport in heterogeneous aquifers are the horizontal integral scale (Ih) of hydraulic logconductivity (Y = ln K) and the longitudinal macrodispersivity (αL). Estimating the former generally requires K measurements in many wells, which are typically not available. In this work, we present a method for estimating Ih and αL from hydraulic tomography (HT) pumping tests. The approach is based on stochastic inversion, assuming head and K are stationary random space functions. Head variance and mean head gradient are calculated from measured head data at large distances from the pumping well and for large times from the start of pumping. These are then used in a theoretical formula for head variance in uniform mean flow to obtain Ih and subsequently αL, requiring only previous knowledge of the vertical integral scale (Iv) and the logconductivity variance σY2 $\left({\sigma }_{Y}^{2}\right)$ . The method is applied to data from an HT pumping test at the Boise Hydrogeophysical Research site and results for Ih and αL are in the range previously reported in the literature. Particularly promising results are found for αL, which is seen to be robust, with relatively little variation for different values of Iv and σY2 ${\sigma }_{Y}^{2}$ .


Introduction
The spatial variability of aquifer hydraulic conductivity K is generally associated with an enhanced rate of solute spreading in groundwater transport, coined as macrodispersion.This is an important topic related to a number of applications such as contaminant and colloidal transport in aquifers, kinetics of geochemical reactions, and groundwater bioremediation.In view of the interest in this topic, a considerable effort has been invested in developing various methods of aquifer heterogeneity characterization using field measurements.In the present paper we focus on the hydraulic tomography (HT) method.It consists of conducting pumping at varying intervals along a number of wells and measuring the pressure response by transducers in numerous discrete intervals of observation wells in and surrounding the volume of investigation.The identification of K spatial distribution from measurements of the pressure head H is known as a solution to the inverse problem.Two main inversion approaches were advanced in the past, as discussed in the following.
The first approach consists of casting the flow equations in a numerical form with K values at the nodes of the spatial grid regarded as unknowns, while the measured H are known at a set of nodes.This HT analysis entails solving an inverse problem with thousands of unknown variables.While theory on HT inverse problems dates back to early publications such as Carrera and Neuman (1986), Gottlieb and Dietrich (1995) and Yeh et al. (1996), it is still a very active field of study (e.g., Cardiff & Barrash, 2011;Hochstetler et al., 2016;Zha et al., 2018;Zhao & Illman, 2018); a thorough review of HT literature and field studies is presented in the work of Cardiff and Barrash (2011).Nevertheless, there are various problematic aspects of the HT inverse methods, the first and foremost of which is that the problem is ill posed (Bohling & Butler, 2010).In addition, to achieve accurate estimation many data measurements are required, which are often not available.The result of HT inversion is usually presented as a spatial map of hydraulic conductivity, and these maps are often misinterpreted as a deterministic output, despite that the method is in fact Bayesian and only offers a "most likely" result.The method is also extremely demanding in terms of computational effort.Finally, the HT results are also difficult to verify in field experiments since pointwise measurements, or independent tests for comparison, are usually scarce and most of the literature focuses on validation using synthetic examples.Thus, only a limited number of studies include application to high resolution 3D field experiments (Berg & Illman, 2011;Cardiff et al., 2012Cardiff et al., , 2013;;Hochstetler et al., 2016;Tiedeman & Barrash, 2020).
Abstract Two important properties for modeling flow and transport in heterogeneous aquifers are the horizontal integral scale (I h ) of hydraulic logconductivity (Y = ln K) and the longitudinal macrodispersivity (α L ).Estimating the former generally requires K measurements in many wells, which are typically not available.In this work, we present a method for estimating I h and α L from hydraulic tomography (HT) pumping tests.The approach is based on stochastic inversion, assuming head and K are stationary random space functions.Head variance and mean head gradient are calculated from measured head data at large distances from the pumping well and for large times from the start of pumping.These are then used in a theoretical formula for head variance in uniform mean flow to obtain I h and subsequently α L , requiring only previous knowledge of the vertical integral scale (I v ) and the logconductivity variance  (  2

𝑌𝑌
) . The method is applied to data from an HT pumping test at the Boise Hydrogeophysical Research site and results for I h and α L are in the range previously reported in the literature.Particularly promising results are found for α L , which is seen to be robust, with relatively little variation for different values of I v and   2  .
In addition to the drawbacks of HT discussed above, some important applications, such as prediction of solute plume spreading by groundwater flow, require global parameters characterizing the permeability spatial variability rather than a pointwise distribution of K.In the same vein, a typical procedure entails generating K values in a large zone of the aquifer, extending past the one investigated by the HT, for instance using geostatistical tools and this also requires a statistical characterization of the K distribution.Then the salient question is whether it is possible to identify these global parameters directly from the distribution of the head measurements, bypassing the need for a pointwise inversion.
To address these questions, a second approach, which we have adopted in the past and pursue here, is a statistical analysis.It entails modeling the aquifer logconductivity field Y = ln K as a stationary random space function which is then characterized by a few statistical parameters: the geometric mean 〈Y〉 = ln K G , the variance   2  and the horizontal I h and vertical I v integral scales, respectively, associated with the assumed axi-symmetric auto-correlation function.They are exhaustive if Y is modeled as multi-Gaussian.The stochastic inversion aims at deriving these parameters from the measured H, which is regarded as a realization of a random field as well.Subsequently, realizations of K values can be spatially generated using standard algorithms such as Sequential Gaussian Simulation.
This work deals with the second approach discussed above.While the general idea has been presented in geostatistical literature (e.g., Firmani et al., 2006;Zech et al., 2015), the current approach is based on discussions in Indelman et al. (1996) and Indelman (2001).It was further developed in recent years (Cheng et al., 2019(Cheng et al., , 2020) ) and applied to a synthetic case in Bellin et al. (2020).It was only very recently, in Cheng et al. (2022a) that the stochastic inversion was implemented in a field experiment of HT, as described in the following.The solution for the head H 0 (x, y, z, t) in a homogeneous unconfined aquifer depends on the three properties: K, the specific storativity S s and the specific yield S y .The approach of Cheng et al. (2022a) was to fit each measured H(t) curve at a given location with an H 0 (t) model, where the three parameters in H 0 are coined as equivalent properties.In simple words, the equivalent properties: K eq , S s,eq , and S y,eq , are the constant properties of an aquifer whose head response fits the measured head response in the actual heterogeneous aquifer.For heterogeneous aquifers the equivalent properties vary in space and the aim of Cheng et al. (2022a) was precisely to analyze the spatial variation of their statistical moments, that is, the mean and variance (a brief recapitulation of the method and results is given in Section 2.3).In principle, the equivalent properties represent a kind of spatial average over the volume extending from the pumping source to the observation point.At a sufficiently large distance R the mean K eq value levels off and tends to the aquifer constant effective conductivity in mean uniform flow, which is a key parameter for modeling flow under natural gradient conditions.The identification of the effective properties was one of the main results of Cheng et al. (2022a).However, there was no estimation of correlation (i.e., integral) scales or macrodispersivity in that work and results were given only for mean and variance of aquifer equivalent properties.
The aim of the present study is to develop a procedure for estimating logconductivity horizontal integral scale I h and longitudinal macrodispersivity α L from HT data and to implement the procedure on the Boise Hydrogeophysical Research Site (BHRS) using data from the 2011 Boise Hydraulic Tomography Test (BHTT).We seek the identification of the parameters which quantify the local K statistics, namely K G ,   2  , I h , I v , by using the same measurements from the BHTT which were used in the analysis of Cheng et al. (2022a).The approach suggested in previous literature (e.g., Bellin et al., 2020) is to derive a solution of the direct problem for the mean head 〈H〉, which satisfies the appropriate boundary and initial conditions and then to identify the parameters by a best fit with the measured mean head   .We attempted to implement this approach by using the solution derived in Indelman (2003) for a constant rate pumping source in an elastic aquifer.Due to the complexity of the problem, the theoretical solution for 〈H〉 was derived by a first order approximation in   2  .We attempted to apply the inverse procedure by matching the solution of Indelman (2003) to the mean head computed from measurements in BHTT for early pumping times in which flow is storativity dominated (i.e., elastic aquifer).This did not succeed due to the difficulty of separating in the measured signal the leading order term from the one related to heterogeneity.In other words, we found that for the BHRS, the heterogeneity is not significant enough to clearly quantify its impact on the mean head during early pumping times.
To deal with this issue we turn to the measurement based head variance   2  , which is much more promising than 〈H〉 since it is entirely related to heterogeneity.The theoretical first order approximation of   2  for flow to a source in an unbounded domain was derived by Severino et al. (2008) in terms of six quadratures, which is much too computationally demanding, rendering its use untenable.Severino (2011), Severino et al. (2019), and Severino 10.1029/2023WR036256 3 of 18 and Cuomo (2020) could reduce the number of quadratures for a line source which represents a fully penetrating well in a confined aquifer, but this is not the case we are interested in (i.e., unconfined aquifer and pumping in an interval).To obtain a theoretical expression related to   2  that can be calculated, we implement an additional approximation, namely assuming that the theoretical   2  pertains to a uniform mean flow driven by the local mean horizontal head gradient  ℎ = ∕ (Di Dato et al., 2019).We use the head measurements at the largest times and distances from the pumping source to calculate   2  , where we expect our assumption of uniform mean flow to be most applicable.Under these conditions, a simple analytical expression provides a relationship between the measurement based   One of the motivations for focusing on identification of I h is that it is typically challenging to estimate because of the scarcity in direct measurements of K in the horizontal direction.The encouraging result presented in this work is that the range of identified I h is in agreement with previous estimates in the literature.An even more promising finding is that the value of the longitudinal macrodispersivity in mean uniform flow   =  2  ℎ is robust, that is, practically independent of the selected   2  values.Moreover, it is in agreement with the values pertinent to a large number of transport experiments analyzed by Zech et al. (2023) for aquifers of a mild level of heterogeneity.It is emphasized that to the best of our knowledge, the present study is a first attempt to identify flow and transport parameters by stochastic inversion for a continuous pumping field test.In view of the encouraging results we conclude the paper by proposing types of future tomographic tests which are tailored to the purpose of identifying transport parameters.
The plan of the paper is as follows.In Section 2.1 we present the method developed in this work.Section 2.2 provides general information on the BHRS and details of the HT field experiment conducted in 2011.This is followed by a summary of previous relevant analysis carried out using the BHTT data in Section 2.3.Sections 3.1 and 3.2 present the implementation of the method presented in Section 2.1 on the BHTT data, with Section 3.1 focusing on the mean head and variance, while Section 3.2 are the results for I h and α L .In Section 3.3, we consider the possibility of reducing the amount of field data and evaluate its impact on the results.Finally, Section 4 presents a short summary, the main conclusions of this work and future recommendations.

Method for Identifying I h and α L
In this section we detail the method for estimating logconductivity longitudinal integral scale and the macrodispersivity from HT data.Overall, the method consists of calculating the head variance from the pumping test field measurements, denoted by   2  , and comparing with a theoretical expression for head variance, denoted by   2  .We note that the method is described here only in general terms, while many additional details are provided in the following sections during its implementation on the BHTT data.
We consider that a tomographic constant rate pumping test has been carried out in an unconfined aquifer of thickness D. Figure 1 depicts the layout of the model corresponding to the field test, with pumping source and observation points appearing alongside the coordinate system.Notice that the pumping location is assumed to be at a point in space, though in reality this is the center of the pumping well segment.Considering a Cartesian coordinate system, the pumping source is assumed to be located at x = y = 0 and z = Z (note that Z < 0), where z = 0 is the top of the aquifer (e.g., the initial water table).We also use a cylindrical coordinate system appropriate for the axisymmetric nature of the problem, with  = (  2 +  2 ) 1∕2 , φ = tan −1 (y/x), z and the pumping location is therefore at R = 0, z = Z.The tomographic test data consists of many head time series H(t) at the various observation points (R, φ, z) and for various tests conducted at different pumping locations Z.The measured pressure head is therefore denoted as H Z (R, φ, z, t).We further assume that previous analysis has been conducted on the HT data and estimates of   2  and I v have been made.To model the pumping test at late times, a theoretical model has been previously developed and used extensively, particularly in Cheng et al. (2022a).The model (H 0 ) assumes a homogeneous aquifer of constant hydraulic conductivity K, specific yield S y and a point source discharge Q located at z = Z with a lower boundary condition of no flow: ∂H 0 /∂z = 0 at z = −D (D is aquifer thickness), and a linearized water table boundary condition of S y ∂H 0 /∂t + K∂H 0 /∂z = 0 at z = 0.The late time assumption allows to neglect the impact of the elastic storage.
The details of the problem and the derivation of the solution have been presented in Dagan (1966) and the final expression is given by were J 0 is the zero'th order Bessel function of the first kind.
An analytical expression of   2  pertaining to mean uniform flow driven by a given constant head gradient J = (J h , J v ) in an unbounded domain, has been derived by Dagan (1989, Equation 3.7.16 therein).Thus, with H = 〈H〉 + h and 〈H〉 = H 0 = −J ⋅ x, the first-order approximation of   2  = ⟨ℎ 2 ⟩ (where 〈〉 denotes averaging) is given by where (3) Equation 2 was obtained for an axi-symmetric exponential logconductivity autocorrelation and additional details of the derivation can be found in Appendix A.
We propose to apply Equation 2 to field data from HT tests in an approximate manner by taking advantage of the asymptotic properties (in space and time) of the mean solution H 0 (Equation 1) for a source pumping from an unconfined aquifer.It was shown by Dagan and Lessoff (2011) that for a fixed R and large t the solution H 0 (Equation 1) tends to the steady one for a source in a confined aquifer.Thus, for sufficiently large R and t, the head becomes vertically uniform and head gradient varies only horizontally with distance from the well, that is, In other words, for a sufficiently large R and t, when I h ≪ R, the flow is slowly varying and can be approximated locally by a steady uniform flow.This approximation will be thoroughly examined and discussed in Sections 3.1 and 3.2.Substituting J v = 0 in Equation 2, the key formula for identification of the statistical parameters becomes where f and F h are defined in Equation 3.
We now detail the steps of the procedure for estimating I h and α L .The first step of the method consists of smoothing the head data signals H(t) by fitting them with a smooth function and filtering out any H(t) curves that appear nonphysical, as will be further discussed in Section 2.2.The second step is to isolate the filtered head data H(t) which pertains to observation locations with approximately the same depth as the pumping location, that is, z ≈ Z.These observation locations are then grouped together based on their distance from the pumping source and their depth, that is, H(t) are grouped with similar R and z ≈ Z values.This is further discussed in Section 2.3.Now the grouped head data have similar R and z, only differing by their φ coordinate.These can be averaged to obtain mean head  (  ) .We note that we remove the subscript Z from the head notation from here on, since we consider only z ≈ Z.Similarly, the grouped measured head data can be averaged to obtain the head variance as follows where the averaging is over the data with varying φ.
The third step of the method is to calculate the value of   2  from the measurements which are approximately in the regime of uniform mean flow with negligible vertical gradient (J v ≈ 0).This is done by taking only the data from the largest times t = t max and largest distance from the pumping source R = R max .These head measurements are used in Equation 5  max) .
The fourth step is to obtain the horizontal head gradient J h from the theoretical solution of H 0 given by Equation 1. H 0 is used to represent the field data average head  (  ) by the following approach.We first calculate the equivalent properties K eq and S y,eq for each head data H(R max , φ, z, t) by a best fit of H 0 to H(t), for the late time portion of the curves, at each location.This is further discussed in Sections 2.3 and 3.1.After averaging over φ we arrive at K eq and S y,eq for each measurement depth z.Substituting these K eq and S y,eq values back in Equation 1provides us with H 0 (R max , z, t) which represents  (max, , ) .Now we can calculate the derivative of the average head at the largest available time, that is, After averaging J h (z) over z we obtain one final value of J h (R max , t max ) for the measurements with the largest time and distance from the pumping well.
The fifth and last step of the procedure is to substitute the value of max) in the left hand side of Equation 4 and values of J h (R max , t max ),   2  and I v (the latter two are assumed known from previous analysis) in the right hand side of Equation 4. This allows to calculate I h .We note that Equations 1 and 4 are central to the method described above and each pertain to different flow conditions; the former describes well flow in a homogeneous medium, while the latter applies to uniform mean flow in a heterogeneous aquifer.However, each are carefully used for a different objective, that is, the former for the horizontal head gradient dominated by the mean flow and the latter for the head variance resulting from heterogeneity.
While the parameter I h is of general interest, its main application is to solute transport and more precisely to its impact on the longitudinal macrodispersivity α L .The latter quantifies the rate of spreading of solute plumes due to the conductivity spatial variability and it has been derived by a first-order approximation in   2  (e.g., Dagan, 1989).Thus, after a transient period, α L stabilizes at the asymptotic constant value given by Despite the first order assumption, this relationship was shown to apply also to large values of   2  (Fiori et al., 2017) based on accurate numerical simulations.Furthermore, it was also shown to be approximately obeyed in a few transport field experiments conducted in highly monitored aquifers (Zech et al., 2023).

Description of BHRS Hydraulic Tomography Pumping Test
The BHRS is a heavily investigated site located in proximity to the Boise River, about 15-km upstream from downtown Boise, ID, USA.It mainly consists of a fluvial aquifer composed of coarse cobble, gravel and sand sediments, confined below by a clay layer and unconfined from above.The thickness of the aquifer is approximately 16 m, with slight variations.The site is described in detail and characterized in numerous previous publications, for example, Barrash and Clemo (2002), Barrash and Reboulet (2004), Moret et al. (2004), Barrash et al. (2006), Moret et al. (2006), Clement and Knoll (2006), Clement and Barrash (2006), Clement et al. (2006) ln  , is reported to have values of 0.36-0.92with an overall value reported to be 0.49 (Barrash & Cardiff, 2013).
We consider data obtained in a field campaign in the summer of 2011 and published previously in Cardiff et al. (2013), Cheng et al. (2022a).Partially penetrating pumping tests were carried out using various 10 cm diameter wells located onsite.The wells onsite are depicted in Figure 2a and consist of 13 wells distributed in a pattern of two rings; an outer ring (wells C1-C6) and an inner ring (wells B1-B6), surrounding the central well A1.Constant rate pumping was carried out, in turn, at various depths with 1 m intervals in three pumping wells: A1, B1, and C1.Furthermore, pressure was recorded in different combinations of seven surrounding wells: B3, C1, C2, C3, C4, C5, and C6.The test setup for an example case in which pumping is carried out at well B1 and observations in wells B3, C3, C4, C5, and C6 is illustrated in Figure 2b.Each observation well was equipped with a packer-and-port string consisting of seven 1 m open intervals separated by a 1 m inflatable packer above and below.This allowed for independent measurements of head as a function of time H(t) at successive 1 m intervals in each observation well.Altogether, 2,472 measurements of H(t), corresponding to different observation and pumping locations, were available for data analysis.Pumping flow rate and duration varied between the different tests, with rates in the range 25.35-60.68L/min and duration between 13 and 15 min.
The data processing procedure is described in detail in Cheng et al. (2022a).The H(t) curves are first smoothed by fitting them with a MATLAB Lowess Smoothing function.Then, the data is filtered by removing curves which are suspected to be nonphysical (see Cheng et al., 2022a for more details).Finally, we consider only head measurements which are approximately at the same depth as the pumping interval location.This was found to help simplify the analysis, reducing the impact of boundaries and allowing to focus on horizontal flow.Thus, if we assume a Cartesian grid in which the pumping location is at x = 0, y = 0, and z = Z, we can describe a head measurement as H Z (x, y, z, t), where z is the elevation of the measurement point and Z is the elevation of the pumping location (see Figure 1).Only H Z (x, y, z, t) data in which |z − Z| < 1 m are used in the analysis of this work and therefore the subscript Z will not be necessary from here on, as we remember that z ≈ Z.We note that the pumping and measurement vertical locations are described as points in space (z and Z) although they represent a meter long interval.This has been found to be a sufficient representation for modeling purposes.After this data selection we are left with only 235 measured head curves, as detailed in Table 1 (see next section for additional information).

Summary of Previous Data Analysis
Herein is a brief description of the data analysis and results of Cheng et al. (2022a), which will be the starting point for this work.The goal of that work was to estimate aquifer equivalent property statistical moments from the measured data.A procedure was presented and carried out for the BHTT data leading to estimates of mean and variance of the equivalent properties: K eq , S s,eq , and S y,eq .Equivalent properties at a given point are defined as those pertaining to a homogeneous medium, under the same boundary conditions and pumping source, which lead to approximately the same head response H(t) as the one measured in the field experiment at that point.Results showed that average equivalent properties decrease with horizontal distance from the pumping well and appear to stabilize at sufficiently large distances, in line with existing theory for K eq .The squared coefficient of variation showed similar behavior, with values indicating a weakly heterogeneous aquifer.Estimated values are presented in Table 2 of Cheng et al. (2022a) and are seen to be in agreement with literature values for BHRS.
A main component of the method proposed in Cheng et al. (2022a) was a separation of each measured head signal H(t) into three time periods.First, is the early time period (0 < t < t E , where t E is the endpoint of the early time period), in which water is drawn from the elastic storage and the impact of the water table is negligible.Second, is the late time period (t > t L , where t L is the initial time of the late time period), in which water is drawn mainly from water table drainage and the elastic storage is negligible.Third, is the intermediate time period (t E < t < t L ) in which both elastic storage and water table drainage are significant.It was shown that separate analysis of the early and late time periods is beneficial because it reduces the number of parameters which need to be identified in an inverse procedure, making the process more efficient and accurate, avoiding the intermediate time period in which models are more error prone.Furthermore, it was shown that results for mean and variance of K eq are consistent between the early time and late time analysis, that is, the approach is robust and reliable.
In this work, our goal is to obtain additional aquifer characterization which was not achieved in Cheng et al. (2022a), namely the horizontal correlation scale of conductivity and the macrodispersivity.This will be carried out by analysis of the late time period data.In Table 1 we present the distribution and quantity of H(t) curves from the BHTT data, which were used for our analysis.The data is clustered into groups within a range of 2 m depths and approximately 0.8 m for distances from the pumping source (R), for example, for Z = −6 m and R = 9.6 m, the range is actually −7 < z, Z < −5 and roughly 9.1 < R < 9.9 m with 10 available H(t) curves from all the pumping tests.This is the remaining data after the processing and selection carried out for a late time period analysis, that is, this was the data previously used in Cheng et al. (2022a) for estimating the moments of K eq and S y,eq .It should be noted that for the data of Table 1, it was found that t L ≈ 200 s.

Background to Implementation of New Method
The stochastic solution to the inverse problem, which we are interested in, consists of identifying the statistical structural parameters characterizing Y = ln K, namely K G ,   2  , I h , I v (or the associated f = I v /I h ), with the aid of measured H statistical parameters.The general approach forwarded in the literature is to derive the solution of the direct problem for the random H (or the associated K eq ), to obtain its statistical moments and subsequently to identify the structural parameters by a best fit of measured and theoretical H moments (see, e.g., Dagan, 1989).We attempted to implement this approach on the BHTT data without success, as explained in the next paragraph.
In the case of a three-dimensional pumping source representing a well, one of the approaches was to derive 〈K eq 〉 as function of R in steady flow and then to identify the structural parameters by a best fit with the measurement based   (e.g., Bellin et al., 2020;Indelman et al., 1996).A similar approach can be taken by using the solution of Indelman (2003) for time dependent 〈H〉 considering a pumping source with constant discharge in an elastic confined aquifer.The two approaches discussed above are based on solutions obtained by a first-order approximation in   2  and have not been tested against field experiments.We have attempted to implement both approaches using the BHTT data and this, as far as we are aware, is a first such endeavor.Our first attempt was to use   calculated by Cheng et al. (2022a, Figure 7 therein), however, it was not conducive since there was no clear trend of   with R. In the second attempt, we used the Indelman (2003) solution and tried to fit the theoretical 〈H〉 with the measurement based   as functions of time in the early period.This attempt was also not successful since the measured   did not allow for separation of the zero order approximation and the heterogeneity induced term of order   2  which make up the theoretical solution.A more promising approach is to fit the theoretical head variance   2  and the measurement based one   2  , since they are directly related to heterogeneity.Deriving   2  for steady source flow by Severino et al. (2008) involved performing 6 quadratures numerically, making it not applicable to our procedure due to computational cost.The number of quadratures was reduced to 4 by Severino et al. (2019) for a line source representing a fully penetrating well in a confined aquifer, which is still numerically demanding and does not apply to the configuration of the BHRS.We therefore turn to the mean uniform flow approximation and the method described in Section 2.1.

Analysis of Head Mean and Variance Based on BHTT Measurements
In this section we implement steps 2 and 3 of the method described in Section 2.1.The uniform mean flow approximation adopted in our method requires considering data from large times and distances from the pumping well.Therefore, we selected from Table 1 the two distances R 1 = 9.6 m and R 2 = 12.4 m.These are the largest distances which still provide a significant number of data sets.Furthermore, we concentrate on head measurements at the large times t 1 = 450 s and t 2 = 600 s; t 2 is the total duration of the test while t 1 was selected for the purpose of comparison.It is recalled that for both t j (j = 1, 2) the regime is that of pumping mainly from the water table specific yield, that is, the late time period (see Section 2.3).
To derive the mean head   (, , ) , we observe that in Table 1 for each R and z the data bank consists of a number of head signals H(R, z, φ, t) which differ in the value of the azimuthal angle φ (remembering that Z ≃ z).Due to the assumed axi-symmetric stationarity, the mean head is obtained by averaging over φ.For example, taking R 1 = 9.6 m and z 3 = −6 m, there are 10 available time dependent measured signals which are averaged to arrive at  (1, 3, ) .Subsequently, the resulting four profiles over the vertical  (, , ) , for i, j = 1, 2 and z k = −2k (k = 1, …, 8), are depicted in Figure 3.It is emphasized that the values for z 1 = −2 m, z 7 = −14 m, and z 8 = −16 m are based on a small number of measurements (see Table 1); furthermore, these measurement locations are close to the aquifer boundaries and influenced by the boundary conditions.For this reason we use only the data in the central layer between the depths −12 m ≤ z k ≤ −6 m, for which the approximation of locally uniform flow is rather justified (see discussion in Section 3.2).
The homogeneous aquifer solution presented in Equation 1 can be used to represent the average head in a heterogeneous aquifer by substituting K and S y with the equivalent properties: K eq and S y,eq .Modeling the BHTT data with Equation 1 implies a point source is used to represent the actual interval, that is, replacing the line source with a point source at the center of the interval.In Cheng et al. (2022a), this was found to be an accurate approximation for R > 4.2 m in Table 1.Following the procedure of Cheng et al. (2022a), the values of the equivalent properties K eq and S y,eq for each measured head signal H(R, φ, z, t) indicated in Table 1 was obtained by a best fit of the measured H(t) with the solution H 0 (R, z, t, Z), where Q was taken as the field test pumping rate value for an interval with center at Z = z.Furthermore, we averaged the values of K eq and S y,eq over φ and over −12 m ≤ z k ≤ −6 m for all the signals of sources located at R 1 or R 2 , with the resulting values given in the legend of Figure 3.These values are consistent with the mean ones presented in Cheng et al. (2022a, Figures 7a and 10a).The resulting four plots of H 0 (R, z, t, Z) as function of z for R 1,2 and t 1,2 are presented in Figure 3 (this is denoted H 0 (R i , z, t j ), where z ≈ Z).It is seen that the agreement between the measurement based   (solid lines) and the H 0 profiles (dashed lines), in the interval of interest −12 m ≤ z k ≤ −6 m, is quite good.Furthermore, for both measured and computed profiles there is practically no difference between those pertaining to t 1 = 450 s and t 2 = 600 s.This result strengthens our assumption that a quasi-steady regime is achieved for t > 450 s at the selected R = R 1,2 .
We examine now the main parameter of interest, namely the head variance, which quantifies the impact of heterogeneity.The calculation procedure is described in Section 2.1.For each of the 4 values of R i , t j and the different 2  = 2 × 10 −6 m 2 for R 2 = 12.4 m, respectively, valid for t > 450 s.These will serve us in the following in order to identify I h and α L .

Identification of I h and α L at BHRS
In Section 2.1, we discussed the approximation which allows us to obtain a simple theoretical expression for   2  (see Equation 4).It involves assuming sufficiently R and t so that the flow is slowly varying and can be approximated locally by a uniform flow.In the previous section, we found H 0 (R, z, t) representing the mean measured head   in the late period for the largest R 1,2 and the largest t 1,2 in the central zone −12m < z < −6 m (see Figure 3).Thus, the idea is to use this H 0 function to calculate J h via Equation 6.In words, we assume that locally the mean flow behaves like a steady uniform one driven by the head gradient prevailing in the late period for the largest R and t attained in the BHTT.The validity of this simplifying assumption is examined next.

Mean Head Gradient (J v and J h ) for Large t and R
We turn to evaluate whether the mean flow pertaining to the BHTT at the values of R 1,2 and t 1,2 in the central zone, can by regarded as locally uniform with negligible vertical head gradient, as far as the computation of   2  is concerned.We will use H 0 to represent   , as depicted in Figure 3.While a full theoretical solution for   2  in source flow may provide a definitive answer, at present we shall only examine the behavior of J h = ∂H 0 /∂R and J v = ∂H 0 /∂z.These two expressions were derived analytically by differentiation of Equation 1 and then the plots were drawn by plugging in values of Q, K eq , and S y,eq pertaining to the BHTT (K eq = 2.4 × 10 −4 m/s and S y,eq = 0.06 are taken from Cheng et al. (2022a)) and examining their behavior for the values of R 1,2 , t 1,2 of interest.Thus, in Figure 5a we depict the dependence of J v upon R at the boundaries of the central layer z = z 3 = −6 m and z = z 6 = −12 m at times t 1 = 450 s and t 2 = 600 s.First, it is clear that the flow is quasi-steady, as the results for t 1 and t 2 practically coincide, particularly for the lower layer (z = −12 m).Furthermore, |J v | decreases toward zero with increasing R and for large distances (R = R 1,2 ) a fairly constant value of J v is seen with varying R and z.This behavior is in line with uniform mean flow in which J v is approximately zero and constant with z and R.
In a similar manner we display the variation of J h with R in Figure 5b.Again there is little difference between J h at t 1,2 , and J h drops with R, varying slowly for R > R 1 = 9.6 m.However, J h slightly varies between z 3 and z 6 and also for R > R 1 , as the assumption of uniform flow is only an approximation.In the computations which follow we adopt as representative values the average J h across z which is close to J h at z = (z 3 + z 6 )/2 = −9 m.The pertinent values are J v = −10 −3 , J h = −2.9× 10 −3 for R 1 = 9.6 m and J v = −8.8× 10 −4 , J h = −2.1 × 10 −3 for R 2 = 12.4 m.It is emphasized that the contributions of the vertical J v and horizontal J h mean head components to   2  of Equation 2 are proportional to   2  and   2 ℎ , respectively.Thus, for the aforementioned values, the relative contribution  (∕ℎ) 2 is equal to 0.12 and 0.17 for R = 9.6 m and R = 12.4 m, respectively.In view of these relatively small values and the various approximations discussed, we may neglect the contribution of J v altogether and regard the flow as horizontal with J h given by the above values.

Estimation of I h and α L
We now substitute the calculated values from Sections 3.1 and 3.2.1 in Equation 4.   1), for t and z indicated in the legend. 10.1029/2023WR036256 11 of 18 R 2 = 12.4 m.Then, we replace the corresponding values of J h = −2.9× 10 −3 and J h = −2.1 × 10 −3 for R 1 = 9.6 m and R 2 = 12.4 m, respectively, in the right hand side of Equation 4. The result is an equation linking the three parameters   2  , I h and I v .We selected I h as the one to be identified while adopting for   2  and I v values from the literature for the BHRS.This choice is motivated by practical considerations regarding aquifer characterization.Indeed, the common approach for estimating the properties consists of collecting K data in a few boreholes, where K is measured along the vertical.This data is typically used for estimates of   2  and I v based on the vertical variogram of Y.In contrast, estimation of I h requires information from many boreholes distributed in the horizontal plane to allow deriving the Y horizontal variogram, which is quite prohibitive in many applications.Furthermore, I h is the key parameter for identification of a L .For this reason we focus on searching I h for a few estimated values of the other parameters. which was estimated in that 3D HT inversion of field data from the BHRS.Finally, the highest value we considered is   2  = 0.5 , which is the common value attributed to the BHRS, obtained from slug test analysis carried out in Barrash and Cardiff (2013).Literature values of I v = 2 m have been presented in Cardiff et al. (2013) and I v = 1.2, 1.5 m in Cardiff et al. (2011).Here, we will consider the range of I v = 1 − 2 m and use a representative overall value of I v = 1.2 m which appears in Barrash and Cardiff (2013).
We begin the analysis with Figure 6, showing a general representative plot of Equation 4, that is, I h /I v as function of   2   2  for the considered R 1,2 and t 1,2 .It is seen that the value of t 1,2 is immaterial, as the curves coincide; however, the results are somewhat different for the R 1,2 values.It is not clear whether this represents a trend related to the approximation of uniform mean flow (which generally applies for R → ∞) or to the discrepancy between the spread of measured   2 ℎ values for the two R values, as observed in Figure 4.At any rate, the differences are not large and in order to assess the I h values pertinent to the BHRS we have averaged them over the two values of R and present the results as a function of I v for each of the selected   2  separately in Figure 7.It is seen that for the literature value of I v = 1.2 m the magnitude of I h varies between 2.7 and 7.2 m for   2  = 0.5 and   2  = 0.18 , respectively.These are in good agreement with values of I h previously stated in the literature for BHRS, which were found to be in the range of 4-10 m (Barrash & Clemo, 2002;Cardiff et al., 2011Cardiff et al., , 2013)).These results are also within the range of values identified in previous literature for other intensively studied aquifers (Table B1 in Zech et al., 2023).Even considering the entire range of 1 m < I v < 2 m in Figure 7, values of I h remain within the interval of roughly 2-8.5 m, showing robustness of the results.Next, we turn to estimating the value of α L from the BHTT, which is one of the main objectives of the present study.We achieved this by using Equation 7, that is, multiplying I h by   2  , and results are presented in Figure 8 as function of I v for the selected three values of   2  .The somewhat surprising and encouraging result is the robustness of α L as reflected by the closeness of its values, irrespective of   2  values.Thus, it is seen that α L ≈ 1.3 m for I v = 1.2 m while it varies between 0.8 and 1.6 m for the entire range of explored I v values.
Based on a thorough analysis of tens of transport field experiments Zech et al. (2023) have suggested to divide aquifers into three classes of heterogeneity level (see Figure 1 . For each level a range of α L values was recommended as a preliminary choice in applications, in absence of information on   2  and I h .Thus, for weak heterogeneity the mean value, based on data of 13 aquifers, is α L = 1.1 m (standard deviation 1.1 m).It is gratifying to find out that the value α L identified in the present study for the BHRS, α L ≈ 1.3 m, is indeed close to the one characterizing aquifers of weak heterogeneity, in agreement with the   2  < 1 criterion.

Impact of Reduced Field Data
The field pumping test configuration underlying the previous analysis includes three pumping (A1, B1, C1) and 8 observation wells (see Figure 2a), with a total number of 66 signals used for identification at R 1 = 9.6 m and R 2 = 12.4 m in the central z zone (see Table 1).The investment in the large number of wells and of measurements is justified for the experimental setting of the BHTT, for high-resolution investigations of contaminated sites and for in-situ remediation of contaminated source zones considering toxic plumes.However, for the application of the procedure to common plume-scale contamination issues (see next section) such an investment is not practical.It is therefore of interest to employ the extensive data basis in order to examine the impact of reducing the number of wells and measurements on the outcome of the analysis.
For the purpose of evaluating the impact of reducing the available field data, we have considered pumping along only a single well and measurements along only two observation wells.The pumping wells taken into account are A1 and B1, separately (see Figure 2a), with measured pressure signals originating from wells C4 and C5 (for pumping at A1) or C3 and C4 (for pumping at B1).This implies using only data pertaining to R = 9.6 m and results in a total number of measured H(t) in the central zone (−12 m < z < −6 m) of 14 for A1 and 28 for B1, respectively.The difference is due to the measurements filtered out, as explained in Section 2.2.It is emphasized that the combination of A1 and B1 signals (42 in total) exhaust their contribution to the total number used in the analysis of the preceding section, since the pumping of C1 did not take any part in the analysis for R = 9.6 m.This is observed by summing the contributing signals in the R = 9.6 m column in Table 1 to find the same number (44).However, the following analysis accounts for the impact of reducing the number of pumping wells and of measurements as we will consider each case separately; pumping from A1 with only 14 signals and pumping from B1 with only 24 signals.The following analysis follows the same steps of the preceding section which dealt with the ensemble of all three pumping wells.
The first step was the derivation of K eq and S y,eq for each H(t) by best fit of H 0 , toward calculation of H 0 representing the average measured head   .Averaging over z in the central zone led to the following values: K eq = 2.18, 2.28, 2.25 × 10 −4 m/s and S y,eq = 0.078, 0.051, 0.057 for A1, B1 and combined A1 and B1, respectively.It is seen that the new values for A1 and B1 are quite robust and close to those of the combined set, which have already been presented in the legend of Figure 3.The same is true for the plots of H 0 as function of z for t = 600 s (Figure 3) which rely on K eq and S eq values, although these are not plotted for A1 and B1 separately here for brevity.We plot the measurement averaged head  () in Figure 9 for R = 9.6 m and t = 600 s alongside the result using all of the data (A1 and B1), as well as the fitted H 0 using K eq and S y,eq from all the data (replicated from Figure 3).It is seen that in the central zone of −12 m < z < −6 m all the curves practically coincide.Although   based on the A1 data solely is slightly less regular, it is still quite close.
We can conclude at this point that the computed H 0 , as well as the measurement based mean head   for the total set, are close to   based on a single pumping well and two observation wells with 28 measured signals and in fairly good agreement even for 14 measured signals.Therefore, we can proceed with the calculations of   2 ℎ for these cases in the following.

The head variance 𝐴𝐴 𝜎𝜎
2  as function of z for the A1 and B1 head measurements is calculated and the results are presented in Figure 10.It is seen that as expected the variance is smaller for A1, with the average value over z of   2  = 3.4 × 10 −6 m 2 , while for B1 the result is   2  = 5.6 × 10 −6 m 2 ; the latter is close to the variance obtained from all measurements of   2  = 6 × 10 −6 m 2 , presented in Figure 4. Thus, the smaller number of A1 signals underestimates the total variance, whereas B1 and the two corresponding observation wells capture it.However, the salient question is what is the impact upon the identified parameters I h and α L .This is achieved by the same procedure that was carried out previously for all the data, using Equation 4 and substituting the above   2  values with the previous test parameters values   2  = 0.18 , 0.36, 0.50, and I v = 1.2 m.The final results are I h = 4.9, 2.5, 1.9 m for A1 and I h = 8, 4.1, 3 m for B1, as compared with I h = 7.2, 3.7, 2.7 m using all the data (Figure 7).Thus, A1 underestimates I h by about 30% the values for the complete set of data while B1 overestimates by 10%.The agreement can be viewed as quite satisfactory in view of the uncertainty affecting I h .
Finally, the key parameter   =  2  ℎ is practically independent of   2  , precisely like in Figure 8 for the total set of data.Comparison for I v = 1.2 m, yields α L ≈ 0.9 m for A1, 1.5 m for B1, and 1.3 m for the total data set (Figure 8).Again, the agreement can be considered as quite good in view of the imprecise determination of α L in field studies (e.g., Zech et al., 2023).

Summary, Conclusions, and Future Application
This work presents a method for estimating logconductivity horizontal integral scale I h and longitudinal macrodispersivity α L of a phreatic aquifer from measured head data obtained in a HT test consisting of pumping from intervals along a few wells and measurement of the head response H(t) by ports along observation wells.The method adopts a stochastic inversion approach: the hydraulic properties are regarded as stationary random space functions characterized by a few statistical parameters and the aim is to identify them with the aid of those derived from the measurements of the random head field.The focus is on Y = ln K, of assumed axi-symmetric auto-correlation, completely characterized by the statistical parameters K G (the geometric mean),   2  (logconductivity variance), and horizontal and vertical integral scales I h and I v , respectively.The main idea of the new method is to calculate the variance of the measured head   2  and the horizontal of the mean head J h from the HT test data and then substitute this in a theoretical relationship pertaining to mean uniform flow driven by a constant head gradient.The method is applied to data from an HT test carried out in 2011 at the BHRS, which included three pumping and eight observation wells.
We began our analysis by averaging the observed head on intervals of R (the horizontal radial distance from the pumping wells) and of z (the vertical distance from the water table at z = 0 extending to z = −D, the bottom), that is, over rings surrounding the pumping wells.This was done for the largest time t = 600 s (and t = 450 s for the sake of comparison) in which the flow regime is dominated by water table drainage.The random head is assumed to be stationary over the rings due to the axi-symmetric nature of the Y field and of the well flow.We focused on the two largest R and on the aquifer central interval −12 m < z < −6 m for which there were many available measurements.We could approximate quite well the mean head by the solution pertaining to a pumping source in a homogeneous phreatic aquifer provided we replaced K and S y by their equivalent values, in the spirit of Cheng et al. (2022a).Attempts to identify the Y statistical parameters by using the expected value of K eq or of H and a first-order solution in   2  of the direct problem were not conducive.Instead, we used the head variance    2  derived from measurements in the central zone.In absence of a simple first-order solution for the well flow problem, we made the additional assumption that at the relatively large R and t considered here   2  is related to the local mean head gradient and to the statistical parameters by the simple relationship pertaining to mean uniform flow in an unbounded domain.Exchanging the spatial   2  and the ensemble   2  , we could obtain a simple relationship between the parameters   2  , I h and I v .Since I h is the most difficult to derive by direct measurements of K due to few wells which are typically available, we concentrated on the estimation of I h and adopted a few values of   2  and I v from previous literature on the BHRS.The resulting range of identified I h values is in agreement with previous estimates, suggesting that the new method presented here is promising.
We then used the same data in order to estimate the longitudinal macrodispersivity by the formula   =  2  ℎ .It turned out to be quite robust, practically independent of the selected   2  values and weakly dependent on I v .The identified α L ≈ 1.3 m for the literature value of I v = 1.2 m is in agreement with the values for aquifers of a similar level of heterogeneity as quantified by   2  .This is an important and encouraging result toward application to solute transport prediction.The assumptions and requirements for implementing the new approach are as follows.First, we assume that H and K are stationary random space functions.Then we use data obtained at large horizontal distances from the pumping well and at long times, assuming that under these conditions the flow is approximately uniform in the mean with a negligible vertical head gradient.Furthermore, we approximate the local head gradients using an analytical solution with parameters  ⟨⟩ and  ⟨⟩ .Finally, we require previous estimations of  ⟨⟩ ,  ⟨⟩ ,   2  and I v , which can be calculated in various methods, for example, the first three are calculated as explained in Cheng et al. (2022a).
Since the present study is a first attempt to apply stochastic inversion to a continuous pumping field test, it is of interest to discuss the future use of the methodology for aquifer characterization.First, we recommend additional theoretical and numerical investigations in order to validate the results based on the various simplifying approximations adopted in the study.This may include developing a simplified first-order solution as well as numerical simulations for deriving the head variance in 3D radial flow, considering confined or phreatic aquifers.Nevertheless, assuming the present results are validated, it is possible to draw a few tentative recommendations following this study.
An envisioned possible characterization configuration to guide plume-scale investigations and treatments, following the results presented in Section 3.3, is a much less extensive and costly pumping test than the full tomographic test.Our results show that the pumping test could consist of only a pumping well and two observations wells at the same distances R (must be sufficiently large R) and at different horizontal angles.We found that for the BHTT this reduced setup, with a total number of pressure signals of around 20, could lead to estimates close to those obtained by the complete configuration.An overall suggested characterization procedure could consist of first measuring K values along the three wells by existing technologies (e.g., Dietrich & Leven, 2009) in order to derive the logconductivity mean and variance, as well as its vertical integral scale.Subsequently, pumping can be carried out along a few small intervals in one of the wells and pressure signals measured at same elevations in the two other observation wells.Based on the mean head and its variance, derived from the signals in the observation wells, the equivalent conductivity, and the resulting horizontal integral scale and longitudinal macrodispersivity shall be derived along the lines of the present study.

Appendix A: Derivation of the Head Variance in Mean Uniform Flow
Two key equations of this work are Equations 2 and 3, which relate the head variance   2  to the various parameters of the problem.It was derived in the book by Dagan (1989, Section 3.7.2).However, the reader interested in the details of the procedure will encounter difficulties since the developments rely on different preceding sections of the book.For this reason we recapitulate here the derivation in a complete and concise manner.

A1. The Log-Conductivity Variance
The random log-conductivity Y = ln K is written as Y = 〈Y〉 + Y′(x) where 〈〉 stands for average, 〈Y〉 = ln K G , K G is the geometric mean assumed to be constant and the random Y′( x  For the exponential ρ = exp(−r′) integration can be carried out analytically and

A2. The Head H(x) First Order Solution
The equation of steady flow satisfied by the random head ∇ ⋅ (K∇H) = 0 can be rewritten as ∇ 2 H + ∇Y′ ⋅ ∇H = 0.The head field is split into H = H 0 + h, where H 0 is the solution for a homogeneous medium (Y′ = 0) satisfying the boundary conditions and h is a perturbation associated with spatial variability.For a mean uniform flow stemming for instance from constant heads applied on two planar boundaries, we have H 0 = −J.x+ const where J is the constant head gradient.Thus, 〈H〉 = H 0 while ∇ 2 h + ∇Y′ ⋅ ∇h = J ⋅ ∇Y′, with h satisfying homogeneous boundary conditions.At first-order in Y′, h satisfies ∇ 2 h = J ⋅ ∇Y′ with neglected terms of order Y′ 2 .Application of the Fourier Transform to the last equation yields the solution  ĥ() = − (  ⋅ ∕ 2 ) Ŷ ′ () and by inversion  () = −(2) −3∕2 ∫ (  ⋅ ∕ 2 ) Ŷ ′ ()exp(− ⋅ ).

𝑯𝑯
From the expression of h above the variance at order   2  is given by The covariance of the Fourier Transforms of the stationary Y′ is given by where δ is the Dirac operator (see for instance Dagan (1989, Section 1.6)).Substitution in the preceding equation yields   2  =  (2) −3∕2  2  ∫ (  ⋅ ∕ 2 ) 2 ρ() , which is constant.We consider a head gradient J(J h , J h , J v ) that is, of horizontal and vertical components.The two components, for flow parallel to the bedding or normal to it, respectively, can be written as Integration can be carried out analytically by switching to spherical coordinates for k′:   ′ 1 =  ′ sin  cos  ,  ′ 2 =  ′ sin  sin  ,   ′ 3 =  ′ cos  , dk′ = k′ 2 sin θ dθdϕdk′ within the limits θ(0, π), ϕ(0, 2π), k′(0, ∞).Substitution of the above expression   ( ′ ) =  (  2 ℎ  ) (8∕) 1∕2 ( 1 +  ′2 ) −2 for an exponential covariance and the triple integration leads to the final expressions of Equations 2 and 3.

Figure 1 .
Figure 1.Schematic illustration of the pumping test model (taken from Cheng et al., 2022a) with pumping source located at (0, 0, Z) and measurement point at (R, φ, z).

Figure 2 .
Figure 2. The Boise Hydrogeophysical Research Site-well layout (taken from Cardiff et al., 2013).(a) The well distribution map from an aerial view.(b) The experimental setup for the case of pumping at B1 and five observation wells.Measurement locations (blue and green circles) and pumping locations (red circles) are indicated.

Figure 3 .
Figure 3. Mean head  as a function of depth z calculated from Boise Hydraulic Tomography Test data for two R and t values indicated in the legend.The fitted H 0 curves (dashed lines) with equivalent properties (K eq = 2.3 × 10 −4 m/s, S y,eq = 0.057 for R = 9.6 m; K eq = 1.9 × 10 −4 m/s, S y,eq = 0.048 for R = 12.4 m) are plotted for reference.The shaded area represents the interval of interest away from the boundaries −12 m ≤ z k ≤ −6 m.
z k we averaged  [ (, , , ) − (, , ) ] 2 over φ, as prescribed by Equation 5, to obtain   2  (, , ) .The results are summarized in Figure 4, which displays the values of   2  along the vertical.It is seen that as expected   2  varies in a more irregular manner with z k than   , with no clear trend and the scatter of values is larger for R 1 = 9.6 m than for R 2 = 12.4 m.Again, the results for t 1 = 450 s and t 2 = 600 s practically coincide, particularly for the latest time period of t 2 = 600.Averaging   2  over the interval of assumed uniformity −12 m ≤ z k ≤ −6 m we arrive at the 2 key values of   2  = × 10 −6 m 2 for R 1 = 9.6 m and

Figure 4 .
Figure 4. Variance of head from measured data for different depths z and for two R and t values indicated in the legend.The overall head variance average in the range of −12 m ≤ z k ≤ −6 m is presented with dashed lines.Circle size is proportional to the number of data points available.
2 in the left hand side of Equation 4 is replaced by its measured values of

Figure 5 .
Figure 5. Vertical (J v , plot (a)) and horizontal (J h , plot (b)) mean head gradient as a function of R, calculated by derivation of H 0 (Equation1), for t and z indicated in the legend.

Figure 6 .
Figure 6.Integral scale ratio I h /I v = 1/f as a function of     2  calculated via Equation 4 for R 1,2 and t 1,2 .Diamonds indicate results for values of   2   2  representative of the Boise Hydrogeophysical Research Site.

Figure 8 .
Figure 8. Macrodispersivity as a function of vertical integral scale for three values of   2  and for t = 600 s.

Figure 7 .
Figure 7. Horizontal integral scale (average results for R 1,2 ) as a function of vertical integral scale for three values of   2  and t = 600 s.

Figure 9 .
Figure 9. Mean head   as a function of depth z calculated from the reduced Boise Hydraulic Tomography Test data via pumping well A1 or B1 only.The   from all data and the fitted H 0 (dashed lines) are plotted for reference.

Figure 10 .
Figure 10.Variance of head for different depths z from the reduced Boise Hydraulic Tomography Test data via pumping well A1 or B1 only.The overall head variance average in the range of −12 m ≤ z k ≤ −6 m is presented with dashed lines.Circle size is proportional to the number of data points available.

𝑣𝑣]
) is the stationary fluctuation.The two-point covariance of Y′ is given by   = ⟨ ′ () ′ ( + )⟩ =  2  () .The auto-correlation ρ is modeled as axi-symmetric, with the common simplified representation ρ(r′), where   ′ 1∕2 , with I h and I v the horizontal and vertical integral scales, respectively, and f = I v /I h < 1, the anisotropy coefficient.A common form of ρ is the exponential ρ = exp(−r′), which will be adopted here without loss of generality (other types are given in Dagan (1989, Section 3.2.3)).We shall use frequently the Fourier Transform  Ĉ =  2  ρ()  , where as usual   () = (2) −3∕2 ∫ ()exp(.) , where dr = dr x dr y dr z , i is the imaginary unit and ∫ stands for the triple integral  ∫ yields ρ(r) = (2π) −3/2 ∫ρ(k) exp(−ik ⋅ r)dk.For the particular case of ρ function of r′ it is useful to adopt the transformation   ′  = ∕ℎ ,   ′  = ∕ℎ, to calculate   2  for each depth z and then results are averaged over different depths.The final result is a single variance value represented by 2 (max,

Table 1
Number of Data Points (After Filtering) for Late Time Period Analysis, Considering Each Pumping Elevation Z and Observation Location R (Within 0.8 m of Specified Values) and for |z − Z| < 1 m