Theory of transient streaming potentials associated with axial-symmetric flow in unconfined aquifers

Summary

We present a semi-analytical solution for the transient streaming potential response of an unconfined aquifer to continuous constant rate pumping. We assume that flow occurs without leakage from the unit below a transverse anisotropic aquifer and neglect flow in the unsaturated zone by treating the water-table as a moving material boundary. In the development of the solution to the streaming potential problem, we impose insulating boundary conditions at land surface and the lower boundary of the lower confining unit. We solve the problem exactly in the double Laplace—Hankel transform space and obtain the inverse transforms numerically. The solution is used to analyse transient streaming potential data collected during dipole hydraulic tests conducted at the Boise Hydrogeophysical Research Site in 2007 June. This analysis yields estimates of aquifer hydraulic parameters. The estimated hydraulic parameters, namely, hydraulic conductivity, transverse hydraulic anisotropy, specific storage and specific yield, compare well to published estimates obtained by inverting drawdown data collected at the field site.

Hydrogeophysics, Hydrology, Permeability and porosity

1 Introduction

In natural porous media, pore water is characterized by an excess of electrical charges that counterbalance the deficiency of charge at the pore water/mineral interface (e.g. Leroy et al. 2007). When water is flowing through such a charged porous material, the excess of electrical charge present in the pore water is dragged along with the pore water. This creates a source current, which in turn is responsible for an electrical field in the conducting material (e.g. Chandler 1981; Sill 1983; Titov et al. 2002). The source current is known as the streaming current, and the associated electrical potential is referred to as streaming potential. This is an example of an electrokinetic effect. The self-potential method, of which the streaming potential method discussed herein is a special case, is a passive method aimed at mapping or recording the resulting distribution of the electrical potential at the ground surface (possibly in boreholes) with a network of electrodes connected to a multichannel, high impedance (10–100 MegaOhm) and high sensitivity (0.1 mV) voltmeter as, for example, in Rizzo et al. (2004).

The self-potential method, therefore, offers a non-intrusive geophysical method that is sensitive to the flow of the ground water. Bogoslovsky & Ogilvy (1973) were among the first to establish that streaming potential anomalies are generated during steady-state pumping test experiments. Recently, several works have shown that clear streaming potential signals are generated during pumping (Darnet et al. 2003; Rizzo et al. 2004; Suski et al. 2004; Titov et al. 2005; Maineult et al. 2008), infiltration (Suski et al. 2006) and dipole aquifer (Jardani et al. 2009) tests. However, general analytical and semi-analytical solutions for such problems have been missing in the near-surface geophysics and hydrogeology literature.

Recently, Malama et al. (2009) developed a semi-analytical solution for the confined aquifer problem. This solution was successfully tested on the streaming potential data presented by Rizzo et al. (2004). In this study, as an extension to the work of Malama et al. (2009), we develop a semi-analytical solution of the self-potential field associated with transient hydraulic response of an unconfined aquifer to continuous constant rate pumping through a fully penetrating well assumed to be a line sink. This analytical solution is tested using data recorded at the Boise Hydrogeophysical Research Site (BHRS) to determine aquifer hydraulic conductivity, specific storage and specific yield, and is shown to fit measured data well and yield parameter values that compare well to published results for the research site.

2 Mathematical Formulation

We consider the streaming potentials that arise due to fluid flow towards a fully penetrating line sink in an unconfined aquifer of finite vertical extent and infinite radial extent. The aquifer is assumed to be homogeneous, transverse anisotropic (radial and vertical hydraulic conductivities are not equal) and bounded from below by an impermeable layer. A schematic of the subsurface considered in this study is shown Fig. 1(a). The governing equation for the fluid flow problem, in dimensionless form, is

1

where s_D = s/H_c is dimensionless aquifer drawdown, s = h(r, z, 0) −h(r, z, t) is drawdown (m), h_D = h/H_c, h is hydraulic head (m), r_D = r/b₂, z_D = z/b₂, t_D = α_rt/b²₂, r is the radial distance from the pumping well (m), z is the vertical position from the base of the aquifer (m), t is elapsed time since onset of pumping test (s), κ = K_z/K_r is the anisotropy ratio, α_r = K_r/S_s is the radial hydraulic diffusivity (m² s⁻¹), K_r and K_z are radial and vertical hydraulic conductivities (m s⁻¹), respectively, S_s is the specific storage (m⁻¹), b₂ is the initial saturated thickness (m), H_c = Q/(4πb₂K_r) is the characteristic head of the system (m) and Q is the volumetric pumping rate (m³ s⁻¹). Eq. (1) is solved subject to the initial condition

2

the far-field boundary condition

3

the no-flow condition at the base of the aquifer

4

and the linearized kinematic condition at the water-table (Neuman 1972)

5

where α_D = κ/ϑ, ϑ = S_y/(b₂S_s) and S_y is specific yield or drainable porosity, which is a measure of the fraction of the bulk volume a saturated porous medium would yield when the water is allowed to drain out under the action of gravity. The complete non-linear form of the kinematic condition is given in Bear (1972). The linearized form given here was shown to simulate the classical delayed gravity response of unconfined aquifer drawdown to constant rate pumping by Neuman (1972). The complete non-linear form of the kinematic condition arises due to the fact that the water-table is a moving material boundary whose position is unknown a priori. The linearization is based on the assumption that the drawdowns are small relative to the initial saturated thickness of the aquifer (Neuman 1972). A major consequence of the linearization is that values of specific yield estimated with the model of Neuman (1972) are typically smaller (see Section 5.2) than what would be expected from porosity values determined by other methods (Moench 1994).

Figure 1.

Schematic of (a) the multilayered subsurface including an unconfined aquifer and (b) the three-layer conceptual model used to develop solution.

The Neumann boundary condition at the pumping well (line sink), in dimensional form, is

6

which, in non-dimensional form, becomes

7

The exact analytical solution to this flow problem in double Hankel—Laplace transform space was obtained by Neuman (1972), and is

8

where η² = (p+a²)/κ, p and a are the Laplace (in t_D) and Hankel (in r_D) transform parameters, respectively, and

9

For the streaming potential response, we consider a three-layer conceptual model consisting of an aquifer bounded below by a more electrically conductive but hydraulically impermeable unit over a highly resistive bedrock, and above by the unsaturated zone, as shown in Fig. 1(b). The bounding unit below the aquifer is typically a clay aquitard, hence the higher electrical conductivity and lower permeability assumed here. In this study, we limit the development to formations where hydraulic anisotropy is more important than electrical conductivity anisotropy, and assume that the latter can be neglected. The governing equation for the transient streaming potential response of the ith layer is (e.g. Revil et al. 2003)

10

where j_i is the electric current density (A m⁻²), and i = 1, 2, 3 with i = 1 being the unsaturated zone, i = 2 the unconfined aquifer and i = 3 the layer below the aquifer (see Fig. 1b). Revil et al. (2003) and references therein, have shown that

11

where σ_i is the electrical conductivity of the i th layer (S m⁻¹), E_i = −∇φ_i is the electric field (V m⁻¹), φ_i = ϕ_i−ϕ_0,i is the electric potential change (V) in i th layer due to pumping in one of the layers, ϕ_0,i is the potential at t = 0, j_s,i is the electric current density due to fluid flow in one of the layers and, in two-dimensions, is given by (e.g. Sailhac & Marquis 2001)

12

where γ is the specific weight of water (N m⁻³), ℓ_i = −σ_iC′_i/γ is the electrokinetic coupling term [m²(V s)⁻¹], C′_i = ∂φ_i/∂h is the streaming potential coupling coefficient (V m⁻¹), K_i is the hydraulic conductivity tensor and q_i = −K_i∇h_i = K_i∇s_i is the Darcy fluid flux (m s⁻¹). Hence, substituting eqs (11), (12) and Darcy's law into eq. (10), for cylindrical coordinates, yields the following governing equations for the streaming potential response in the aquifer (i = 2), given here in dimensionless form

13

and for i = 1, 3:

14

where, for i = 1, 2, 3, φ_D,i = φ_i/Φ_c and Φ_c = (γℓ₂/σ₂)H_c. Note that the term is square brackets in eq. (13) does not appear in eq. (14) because it is assumed in this study that there is no fluid flow in layers 1 and 3. We solve these equations subject to the initial condition

15

and the boundary conditions

16

at far-field in all layers,

17

(insulating condition) at ground surface,

18

(insulating condition) at the base of layer 3. Here b_D,1 = b₁/b₂ and b_D,3 = b₃/b₂ are the normalized thicknesses of the layers 1 and 3. Note that the insulating condition in eq. (18) is applied at the contact between layer 3 and an underlying highly resistive bedrock (see Fig. 1a). It has been shown in Malama et al. (2009) that the line sink condition at the pumping well can be written as

19

for the streaming potential problem. Upon non-dimensionalization eq. (19) becomes

20

In dimensionless form, continuity conditions at the upper and lower boundaries of the aquifer become

21

22

23

24

where σ_D,i = σ_i/σ₂, the electrical conductivity of the i th layer normalized by that of layer i = 2. In the following section, the analytical solution to this problem is developed.

3 Analytical Solution

Taking the Laplace and Hankel transforms (see Appendix A for definition of the latter) of eqs (13) and (14) then and solving subject to the conditions given in eqs (15)–(24) leads to the following solutions for the Laplace—Hankel transforms of dimensionless electric potential

25

where

26

27

28

29

30

31

32

33

34

The inverse transforms were obtained numerically as discussed in Appendix B. This solution describes the transient streaming potential field generated by pumping water from an unconfined aquifer. It can be used to estimate aquifer hydraulic parameters from measurements of streaming potentials at ground surface as we demonstrate in Section 5 below.

4 Model Predicted Response

4.1 Response in a radially infinite domain

The response predicted by the above solution is plotted in Fig. 2, where (a) φ_D,1 and (b) φ_D,2 are plotted against t_D/r²_D for different vertical positions. The results were computed with dimensionless parameters κ = 1.0, ϑ = 31.3, σ_D,1 = 0.05, σ_D,3 = 40.0, b_D,1 = 0.25 and b_D,3 = 0.5. These parameter values are representative of an alluvial aquifer (e.g. the aquifer at the BHRS) bounded above by an unsaturated zone and below by an impermeable and more electrically conductive clay layer and a highly resistive Basalt bedrock. The results shown in Fig. 2 indicate that the maximum streaming potential response is at the water-table (z_D = 1.0). This is due to the desaturation of the pores at the water-table. It should be noted that streaming potentials arise due to drag of excess charges in the diffuse layer. This drag of charge is maximized when the pores are drained of water completely, as is the case at the water-table. Below the water-table the pore space remains saturated.

Figure 2.

Plot of dimensionless streaming potentials φ_D,i for layers (a) i = 1 and (b) i = 2, against t_D/r²_D, for different values of z_D.

For the zone above the water-table (see Fig. 2a) the streaming potential signal decreases upward from the water-table (increasing z_D), attaining a minimum at ground surface (z_D = 1 +b_D,1). This is due to the increasing distance from the source of the streaming potentials, namely flow in the saturated zone. For the saturated zone (see Fig. 2b), the predicted response decreases downward from the water-table and is minimum at the interface with the more electrically conductive unit below the aquifer. This is attributable to the Neumann boundary condition of no-flow imposed at this interface for the flow problem. The vertical fluid flow velocities vanish at this lower boundary, and thus, the resulting streaming currents are minimized.

Fig. 3 shows the model predicted dimensionless streaming potentials generated in layers i = 1, 2 for different values of the dimensionless storage parameter ϑ = S_y/(b₂S_s), using the parameters κ = 1.0, σ_D,1 = 0.05, σ_D,3 = 40.0, b_D,1 = 0.25 and b_D,3 = 0.5. The parameter value ϑ = 0 corresponds to the confined aquifer solution presented in Malama et al. (2009). The results presented here indicate that streaming potentials decrease in magnitude as one decreases the ratio of the volume of water that can be released by water-table decline (measured by S_y) to the volume of water released from elastic storage (b₂S_s). This implies that low streaming potential signals are to be expected when pumping tests are conducted in media with small specific yield values.

Figure 3.

Plot of dimensionless streaming potentials, φ_D,i, versus t_D/r²_D, for different values of the dimensionless storage parameter ϑ = S_y/(b₂S_s). (a) i = 1, z_D = 1 +b_D,1 and (b) i = 2, z_D = 0.5.

The effect of the hydraulic conductivity anisotropy ratio, κ = K_z/K_r, is shown in Fig. 4(a). The model results shown in the figure indicate that κ has significant impact of the temporal evolution of streaming potentials. This is because the parameter has significant effect on the contribution of vertical flow in the aquifer to the streaming potential response at ground surface. Fig. 4(b) shows the drawdown response of the aquifer for different values of κ. The plots in the figure show that high potentials and drawdowns are associated with small values of κ at early and intermediate times. This is due to the fact that for values of κ < 1.0, a given pumping rate induces large vertical hydraulic gradients to compensate for small values of K_z, which in turn lead to large values of streaming potentials at early and intermediate times. For κ > 1.0 the converse is true as low vertical hydraulic gradients develop at early time due to large K_z values. At late time, when flow in the aquifer is predominantly radial, both drawdown and streaming potential responses merge to the response controlled by the radial hydraulic conductivity.

Figure 4.

Plot of dimensionless (a) streaming potentials at land surface, φ_D,1(r_D, z_D = 1 +b_D,1, t_D) and (b) aquifer drawdown, s_D(r_D, z_D, t_D), against t_D/r²_D, for different values of the hydraulic conductivity anisotropy ratio κ = K_z/K_r.

Surface measurements of streaming potentials along a transect, as in Bogoslovsky & Ogilvy (1973), show a cone of potentials about the pumping well. Fig. 5(a) is a plot of model predicted dimensionless potentials at ground surface against the radial distance from the pumping well. The model results are plotted for different values of the dimensionless time, t_D, showing the temporal evolution (outward propagation) of the cone of potentials around the pumping well. The figure shows that streaming potentials decay rapidly with radial distance from the pumping well. The corresponding temporal evolution of the cone of vertically averaged aquifer drawdown is shown in Fig. 5(b). Whereas the cones of drawdown show a sharp peak at r_D = 0, those of the streaming potentials predicted at land surface show a slope that approaches zero as r_D→ 0. This is due to the symmetry boundary condition in eq. (19) imposed at r_D = 0 for layer i = 1.

Figure 5.

Plot of (a) dimensionless streaming potentials at ground surface, φ_D,1(r_D, z_D = 1 +b_D,1, t_D) and (b) the depth averaged unconfined aquifer drawdown, ⟨s_D(r_D, t_D)⟩, against r_D, at different values of the dimensionless time, t_D.

The streaming potential response at ground surface for different values of the normalized thickness of the unsaturated zone (or depth to water-table) is plotted in Fig. 6, using the parameters ϑ = 31.3, κ = 1.0, σ_D,1 = 0.05, σ_D,3 = 40.0 and b_D,3 = 0.5. According to these results, there is about an order of magnitude decrease in the strength of the streaming potential signal at ground surface for two orders of magnitude increase in the depth to the water-table. Hence, for sites where the unsaturated zone is very thick improvements in electrode sensitivity and by the use of advanced noise filtering would be crucial to separate the streaming potentials associated with the pumping test of interest from background noise.

Figure 6.

Plot of dimensionless streaming potentials, φ_D,i, versus t_D/r²_D, for different values of the dimensionless thickness, b_D,1, of the zone above the water-table.

4.2 Response in an aquifer bounded by a constant head or no-flow boundary

It is not uncommon for a shallow fluvial aquifer to be bounded by a river, which acts as a constant head boundary, or by an impermeable fault that acts as a no-flow boundary. Such a flow domain can be conceptualized as semi-infinite. The solution for the radially infinite flow domain can be used to predict behaviour in the semi-infinite domain by invoking the principle of linear superposition in conjunction with the method of images (Carslaw & Jaeger 1959; Kresic 2007). To simulate the constant head (zero drawdown) boundary, an image source (injection well) is placed at a distance equal to that of the sink (pumping well) from the river boundary as shown in Fig. 7 (with the x-axis oriented perpendicular to the boundary). To simulate a no-flow boundary one uses an image sink. This approach can be extended to the streaming potential response since the governing equation of the streaming potential problem is also linear. Hence, one can use eq. (25) to compute the response at the point x_obs = (x_obs, y_obs, z_obs) according to

35

where and r and r* are the respective radial distances from the real and image wells. Fig. 8 shows the predicted behaviour for different values of d_D = d/b₂, where d is the distance to the boundary from the pumping well. The model results show that the effect of the boundary is to reduce the SP signal if the boundary condition is constant head, and to increase the signal if it is no-flow. As d_D increases, the solution approaches that for the infinite domain.

Figure 7.

Sketch of real and image well to simulate flow in semi-infinite domain bounded by a constant head or impermeable boundary, where (x_obs, y_obs) are coordinates of a measurement location.

Figure 8.

Predicted SP response for a no flow and constant head boundary at different values of d_D, the normalized distance from pumping well to boundary.

5 Application to Field Data

5.1 Background

The data analysed here were obtained during a two week period in 2007 June at the BHRS where a series of pumping tests monitored with observation wells and two electrodes arrays for streaming potentials and electrical resistivity tomography, were conducted (Cardiff et al. 2009; Jardani et al. 2009). The multilayered shallow subsurface system at the BHRS consists of an unsaturated zone with an average thickness of about 3 m, and an unconfined aquifer separated from a deeper confined aquifer by a clay layer (Barrash & Reboulet 2004). Geophysical surveys at the site indicate that the deeper aquifer is underlain with a highly resistive basalt bedrock, consistent with the regional geological setting (Squires et al. 1992; Othberg 1994; Barrash et al. 1997). The shallow unsaturated zone and unconfined aquifer at the site, which has a vertical extent not exceeding 20 m, consists of unconsolidated cobble and sand fluvial deposits (Barrash & Reboulet 2004). The clay unit, which underlies the unconfined aquifer, is continuous at the BHRS and has a vertical extent of about 10 m, estimated from results of geophysical investigations using seismic and electrical methods. Data and inversion results from previous geophysical experiments conducted at the site indicate that the confined aquifer below the clay unit has a thickness of about 20–25 m and consists of unconsolidated silty sands. The site is bounded to the west by the Boise River. Presently, there are 18 wells in place at the BHRS, all of which are screened over the entire saturated thickness of the unconfined aquifer. Of the 18 wells at the BHRS, 13 are located within a 20 m diameter central area of the site (Fig. 9). These centrally located wells are arranged in two near-concentric rings around a central well, and provide a high density of subsurface geological, geophysical and hydrologic information for detailed understanding of the shallow unconfined aquifer at several spatial scales (Clement et al. 1999; Barrash et al. 2006).

Figure 9.

Aerial photo of the BHRS and existing wells (inset) at the site (after Barrash & Reboulet 2004).

Ten dipole tests, in which water was pumped from one well and injected into another, were conducted in 2007 June at the BHRS. In this study, we analyse data from the test conducted on June 27, where water was pumped from well C1 and injected into well C4. Streaming potentials were measured with a Keithley 2701 multichannel voltmeter, and 80 non-polarizing Pb/PbCl2 (Petiau) electrodes (Petiau 2000) were used as discussed in Jardani et al. (2009). These were scanned every 30 s during the course of each dipole experiment. A reference electrode was placed at a distance of 90 m from the central well area of the BHRS where the experiments were conducted. Seven additional electrodes were placed in boreholes. Fig. 10 is a map of the positions of the electrodes used for streaming potential measurement, stainless steel electrodes for the resistivity tomography, and wells in pressure transducers were placed for direct head measurement. The initial streaming potential distribution, ϕ_0,i, recorded prior to dipole pumping/injection tests is also shown in the figure. Ordinary kriging was used to interpolate the electrode data yielding the smooth field shown in the figure.

Figure 10.

Map of the positions of Petiau electrodes (*) used for streaming potential measurement, stainless steel electrodes (○) for the resistivity tomography and wells (•) in the central area of the BHRS. The streaming potential distribution prior to dipole pumping/injection tests is also shown. The dashed line separates two areas of relatively higher and lower transmissivities (after Jardani et al. 2009).

5.2 Parameter estimation results

The data are analysed using eq. (35) for flow in a semi-infinite domain with a constant head boundary. Jardani et al. (2009) report laboratory results of material properties in which they determined, using the methodology used by Suski et al. (2006) and Bolève et al. (2007), that the average electrolyte conductivity at the BHRS is 2.21 × 10⁻² S m⁻¹ (at 25°) and the streaming potential coupling coefficient is C = −γℓ₂/σ₂ = − 13.4 ± 3.4 mV m⁻¹. These values were used in the estimation of aquifer parameters discussed here.

Eq. (35) was used with the parameter estimation software PEST (Doherty 2002) to estimate aquifer hydraulic parameter K_r, K_z, S_s and S_y. To simulate the dipole setup, the effect of the injection well and its image (due to river boundary) were linearly superposed with the effect of the pumping well and its image as discussed in Section 4.2 above. The data analysed here were collected on June 27 during the dipole test where water was pumped from well C1 and injected into C4. The flow rate during the tests was maintained at a constant rate of 4.1 × 10⁻¹ m³ s⁻¹ (65 gpm). The saturated thickness, which corresponds to the initial head, of the aquifer during the tests was 16 m. The maximum drawdowns and build-ups, recorded in the pumping and injection wells, respectively, averaged about 75 cm.

We analyse the data from electrodes situated in the neighbourhood of the injection well (C4), which is where the streaming potential response was the most pronounced during the experiments. By ‘neighbourhood of injection well’ is meant points within a radial distance (relative to the injection well) equal to half the separation distance between the injection and pumping wells of a particular dipole set-up. For the data analysed here, where the separation distance between wells C1 and C4 at the BHRS is 16.8 m, this implies electrodes with a radius of 8.4 m from the injection well (C4). The reasons for the asymmetry in the streaming potential response (i.e. larger potential responses in the neighbourhood of injection well than in the neighbourhood of pumping well) during the dipole tests have been elucidated by Jardani et al. (2009). The asymmetry has been attributed to the Reynolds number of Re > 1.0 associated with the inertial laminar flow regime in the vicinity of the pumping wells in these experiments (Jardani et al. 2009).

Fig. 11 shows the fit of the model to the transient streaming potential data at the conclusion of the parameter estimation procedure. The dashed lines indicate bounds of two standard deviations of the residuals on model fit at optimality. The R²-value, which is measure of goodness-of-fit of the model to the data, is included in the plots shown in the figure. The hydraulic parameter values obtained from the streaming potential data are summarized in Table 1 and those determined directly from drawdown data, using the model of Neuman (1972), are shown in Table 2. The drawdown data used to estimate the parameters in Table 2 were obtained from separate classical pumping tests performed at the BHRS to avoid, without loss of generality, the extra complexity in the analysis introduced by dipole tests, and also to obtain independent (by a different set of tests) estimates of hydraulic parameters.

Figure 11.

Inversion results showing model fits to data collected on 2007 June 27 and the estimates of aquifer parameters using data from selected electrodes. The dashed lines are bounds on model fit of two standard deviations of the residuals.

Table 1.

Hydraulic parameter values estimated from streaming potential data for the indicated electrodes (r′ is radial distance from injection well (C4) to respective electrode).

Table 2.

Hydraulic parameter values estimated with the model of Neuman (1972).

To provide a measure of parameter estimation uncertainty, we summarize the ‘composite sensitivities’ of simulated responses to the parameters, at the completion of the optimization process, in Table 3 and the parameter estimation variances normalized by the estimated parameter values in Table 4. The parameter estimation variances are normalized to allow for meaningful comparison of the estimation uncertainties of parameters whose values differ by orders of magnitude. The sensitivities in Table 3 are not the elements of the Jacobian matrix. Rather, they are the so-called ‘composite sensitivities’ of simulated responses to the parameters. They are the normalized magnitudes of the column vectors of the Jacobian matrix, with each element of the vector multiplied by the corresponding weight of the observation data with which it is associated (Doherty 2002).

Table 3.

Composite sensitivities of simulated responses to estimated parameters at optimality.

Table 4.

Normalized estimation variances for parameters estimated from streaming potential data.

5.3 Discussion of results

The results in Fig. 11 show that the model, at parameter estimation optimality, fits the field data well as indicated by R²-values of about 0.99 for all data analysed. Most of the measured streaming potentials fall with the bounds of two standard deviations at optimality. The hydraulic conductivities estimated from transient streaming potential data are of the same order of magnitude as the published values for the aquifer at the BHRS and those reported Table 2, which were obtained directly from head data collected in observation wells B2, B4, B5 and B6 (see data and model fits in Fig. 12) during a pumping test conducted at the research site. For instance, Barrash et al. (2006) report aquifer hydraulic conductivities, estimated from pumping test drawdown data, averaging about K_r = 7 × 10⁻⁴ m s⁻¹ and K_z = 3 × 10⁻⁴ m s⁻¹, and Table 2 gives average values of K_r = 4.6 × 10⁻⁴ m s⁻¹ and K_z = 6.5 × 10⁻⁴ m s⁻¹. The values estimated here from streaming potential data average K_r = 2 × 10⁻⁴ m s⁻¹ and K_z = 3 × 10⁻⁴ m s⁻¹. The values of the specific storage estimated from streaming potential data are on average larger than those estimated from drawdown data. Barrash et al. (2006) report a specific storage value of the order of 1 × 10⁻⁴ m⁻¹, and results in Table 2 have an average value of S_s = 2.6 × 10⁻⁴ m⁻¹. The average value estimated here from streaming potential data is of the order of 1 × 10⁻³ m⁻¹.

Figure 12.

Drawdown data and fits of the Neuman (1972) model to the data collected during a pumping test conducted at the BHRS. For these data, the aquifer was pumped from well A1 at a constant rate of 2.84 × 10⁻³ m³ s⁻¹.

The variability in parameter estimates reported in Table 1 is probably due to aquifer and shallow subsurface heterogeneity. It should be noted that though the parameter values reported here were estimated from point streaming potential data, they represent spatially averaged values on support volumes that encompass the pumping well and the measurement point. Such support volumes typically cover the heterogeneities in the space between the observation point and the pumping well (Sanchez-Vila & Tartakovsky 2007). That the aquifer at the BHRS is heterogeneous, has been discussed in great detail in the work of Barrash & Clemo (2002).

The estimates of specific yield shown in Table 1 are similar to those obtained by Barrash et al. (2006) and those given in Table 2, and fall in the range characteristic of unconfined aquifer models that use the linearized kinematic condition at the water-table in the manner of Neuman (1972). It has been observed in the hydrogeology literature (Nwankwor et al. 1984, 1992; Moench 1994) that treating the water-table in the manner of Neuman (1972) leads to estimates of specific yield that are smaller than what would be expected from porosity values determined by other methods. Alternative models have been proposed that modify the boundary condition at the water-table (Moench 1994) or account for flow in the unsaturated zone above the water-table (Tartakovsky & Neuman 2007). Such approaches as that of Moench (1994) are empirical and introduce additional parameters that have no clear physical meaning. The approach of Tartakovsky & Neuman (2007), though based on sound theory, adopts a linearization of Richards equation that can only be solved analytically assuming the unsaturated zone is vertically infinite. Hence, in this study, we opt for the much simpler approach of linearizing the kinematic condition at the water-table according to Neuman (1972).

It should also be noted that parameter values estimated from the streaming potential data (Table 1) show more variability from electrode to electrode than those estimated from drawdown data (Table 2). The streaming potential data are point data whereas the drawdown data are depth-averaged in each well, since all the wells at the research site are screened across the entire saturated thickness of the aquifer. The depth-averaged data have the effect of smoothing out the effect of heterogeneities. In addition to greater variability of parameter estimates, there is also greater variation among the streaming potential responses (Fig. 11) recorded at the electrodes than in the hydraulic responses (Fig. 12). For instance, streaming potential data at electrodes e-6, e-12 and e-13 show monotonic increase with time whereas those recorded at electrodes e-3–e-5 show a rapid increase initially to some maximum potential, followed by a slow decay with time. The data that show monotonic growth with time also yield larger values of the hydraulic conductivity anisotropy ratio (κ > 1.0). These results agree well with model simulated behaviour plotted in Fig. 4, showing the effect of κ on streaming potential response. The fact that the field data yield such varied anisotropy ratios over such small distances about the injection well is an indication of the significant subsurface heterogeneity at the research site.

The relative sensitivities summarized in Table 3 show that radial hydraulic conductivity is on average the most sensitive of the four hydraulic parameters to streaming potentials. Owing to its higher sensitivity, it is estimated with the lowest uncertainty as can be seen from the values of the estimation variances in Table 4. This is to be expected as the pumping/injection wells used in the dipole experiment are all fully penetrating and induce greater hydraulic gradients in the radial direction. Additionally, the estimation variances of hydraulic conductivity associated with the electrodes at which the streaming potentials increased monotonically with time, are on average three orders of magnitude smaller than those associated with the other electrodes. These are also the data that yielded hydraulic anisotropy ratios that are larger than 1.0. The estimation variance for specific storage was smallest at electrode e-12, which was closest (r′ = 1.24 m) to the injection well (C4). The variability in estimation variances may also be attributed to subsurface heterogeneity as well as to differing noise content of the data collected with the different electrodes.

6 Conclusion

The objective of this study was to develop a semi-analytical solution to the problem of transient streaming potentials associated with pumping water from an unconfined aquifer. We adopted a three-layer conceptual model, consisting of a homogeneous unconfined aquifer and an impermeable layer below the aquifer, all having finite vertical extent and infinite radial extent. The aquifer was taken to be transverse anisotropic in terms of hydraulic conductivity but isotropic in terms of electrical conductivity. Flow in the unsaturated zone was neglected and the water-table was treated as a moving material boundary in the manner of Neuman (1972). The method of images was used to extend the solution to semi-infinite (in the horizontal direction) domains bounded by a constant head or constant flux boundary condition. Model computed results show that the streaming potential response at ground surface is sensitive to aquifer hydraulic anisotropy as well as to the depth from ground surface to the aquifer.

The solution was applied to field data collected during dipole (pumping-injection) aquifer tests conducted in 2007 June at the BHRS. The variability in temporal evolution of streaming potentials observed in the data is attributable to heterogeneity in the hydraulic conductivity field. It is strongly dependent on the hydraulic anisotropy ratio κ. Model simulations as well as model application to field data indicate that for κ≥ 1.0, streaming potentials increase monotonically with time. On the other hand, for κ < 1.0, when large vertical hydraulic gradients develop at early-time to compensate for low values of vertical hydraulic conductivity, the streaming potentials increase rapidly initially, before showing slow decay to values controlled by late-time radial flow.

Application of the solution to field data yielded estimates of hydraulic conductivity, hydraulic anisotropy, specific storage and specific yield. Estimates of hydraulic conductivity and specific yield were found to compare well with published values for the research site determined from pumping tests. Estimates of specific storage were also obtained but these were, on average, an order of magnitude larger than published results from pumping tests conducted at the research site. These large values may be due to the fact that flow in the unsaturated zone above the water-table is neglected. Since electrodes are placed at ground surface, they are sensitive to the capacity of the unsaturated zone to store water. Nevertheless, the study presented here demonstrates that one can estimate the hydraulic conductivity, specific storage and specific yield of unconfined aquifers from transient self-potential data. Since such measurements are usually conducted on the surface and instrumentation is only minimally invasive, the solution has the potential for rapidly yielding preliminary estimates of aquifer hydraulic properties where observation wells are unavailable for direct drawdown response measurement.

Acknowledgments

The study presented here was supported in part by EPA grant X-960041-01-0. We thank the two anonymous reviewers whose insightful comments helped improve the paper. We also thank Dr. Warren Barrash for organizing the experimental aspects of the study presented here. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy's National Nuclear Security Administration under contract DE-AC04-94AL85000.

Appendices

Appendix A: The Hankel Transform

The zero-order Hankel transform, f*(a), of a function, f(r_D), which we refer to in this study simply as the Hankel transform, is given by

(A1)

where a is the real-valued Hankel parameter and J₀ is the zero-order Bessel function of the first kind. The inverse Hankel transform of f*(a) is defined as

(A2)

A particularly useful relation, adopted from Neuman & Witherspoon (1968) and used in this study, is

(A3)

Appendix B: Numerical Inversion of Laplace and Hankel Transforms

The numerical approach used to invert the double Laplace-Hankel inversion is an improvement on the method in Malama et al. (2009). In this study the Laplace transform is now inverted using the Stehfest algorithm (Stehfest 1970) and a portion of the Hankel integral is integrated using an accelerated form of tanh-sinh quadrature (Takahasi & Mori 1974). The infinite integral in eq. (A2) for inverting the Hankel transform is evalutated in two steps (Wieder 1999). One finite (0 ≤ a ≤ j_0,1), the other approaching infinite in length (j_0,1 ≤ a ≤ j_0,m), where j_0,n is the nth zero of J₀(ar_D) and m is typically 10–12. The finite portion is evaluated using tanh-sinh quadrature Takahasi & Mori (1974) where each level k of abcissa-weight pairs uses a step size h = 4/2^k and a number of abcissa N = 2^k, similar to the approach of Bailey et al. (2000). Evaluating the function at abcissa related to level k, the function evalutations at lower levels k−j are simply the 2^k−jth terms in the series (only requiring the weights to be computed for each level, reusing the function evaluations from the largest k). Richardson's deferred approach to the limit (Press et al. 2007) is used to extrapolate the approximate solutions to the limit as h→ 0 using a polynomial extrapolation. The infinite portion of the Hankel integral is integrated between successive j_0,n using Gauss—Lobatto quadrature and these alternating sums are accelerated using Wynn's ε-algorithm.