A semi-analytical solution for transient streaming potentials associated with confined aquifer pumping tests

Summary

We consider the transient streaming potential response due to pumping from a confined aquifer through a fully penetrating line sink. Confined aquifer flow is assumed to occur without fluid leakage from the confining units. However, since confining units are typically clayey, and hence more electrically conductive than the aquifer, they are treated as non-insulating in our three-layer conceptual model. We develop a semi-analytical solution for the transient streaming potentials response of the aquifer and the confining units to pumping of the aquifer. The solution is fitted to field measurements of streaming potentials associated with an aquifer test performed at a site located near Montalto Uffugo, in the region of Calabria in Southern Italy. This yields an average hydraulic conductivity that compares well to the estimate obtained using only hydraulic head data. Specific storage is estimated with greater estimation uncertainty than hydraulic conductivity and is significantly smaller than that estimated from hydraulic head data. This indicates that specific storage may be a more difficult parameter to estimate from streaming potential data. The mismatch may also be due to the fact that only recovery streaming potential data were used here whereas head data for both production and recovery were used. The estimate from head data may also constitute an upper bound since head data were not corrected for pumping and observation wellbore storage. Estimated values of the electrical conductivities of the confining units compare well to those estimated using electrical resistivity tomography. Our work indicates that, where observation wells are unavailable to provide more direct estimates, streaming potential data collected at land surface may, in principle, be used to provide preliminary estimates of aquifer hydraulic conductivity and specific storage, where the latter is estimated with greater uncertainty than the former.

Electrical properties, Hydrogeophysics, Hydrology

1 Introduction

Streaming potentials (also commonly referred to as self-potentials) are observed when a fluid flows through a capillary tube or porous medium; they arise due to the existence of an electric double layer at the solid–fluid interface. When fluid flow occurs, current arises due to the drag of the excess of charge present in the Gouy-Chapman layer. The divergence this source current yields streaming potentials (Sill 1983). Because of the coupling between fluid flow and streaming potential, several workers (Titov et al. 2002; Rizzo et al. 2004; Suski et al. 2004) have attempted to use streaming potentials measured in the neighbourhood of a pumping well to estimate the hydraulic properties of the porous medium. For example, Revil et al. (2003) and Darnet et al. (2007) analysed data obtained by Bogoslovsky & Ogilvy (1973) in an attempt to determine the hydraulic head distribution associated with the subsurface flow problem as well as to estimate the subsurface hydraulic conductivity. Most such attempts have been made only for the case of steady-state flow and steady-state streaming potentials (Sailhac & Marquis 2001; Rizzo et al. 2004; Suski et al. 2004). Particularly, these workers have only considered steady-state flow conditions in developing solutions for analysing experimental data.

Rizzo et al. (2004) used a first-order analysis of a steady-state solution to obtain an approximate linear relation between streaming potential and drawdown in a confined aquifer during the recovery period (after pump shutdown). Using this relation they obtained an approximate equation for transient streaming potential that is valid only for small variations in the piezometric surface and at late time during recovery. The approximate linearized solution can only be used to analyse streaming potential data associated with the recovery period of a pumping test experiment and would require one to pump for long periods (without recording self-potential data), until a steady-state is attained. Their approach only yields estimates of hydraulic conductivity but not specific storage. Titov et al. (2005) used numerical modelling to analyse self-potential signals associated with a pumping test and to estimate aquifer hydraulic properties. Recently, Straface et al. (2007) used a method referred to as the successive linear estimator (SLE), which is an iterative geostatistical inversion scheme developed by Yeh et al. (1996) and Zhang & Yeh (1997), to estimate aquifer hydraulic properties from hydraulic head and self-potential signals associated with a pumping test, using the model of Rizzo et al. (2004).

In this work, we develop a semi-analytical solution for the transient streaming potential response of a three-layered system, consisting of an aquifer and two confining units, due to pumping of the aquifer. In the solution developed here confined aquifer flow is assumed to occur without fluid leakage from the adjacent confining units. However, given that confined aquifers are typically bounded by more electrically conductive clay or clay-rich units, we develop the solution using a three-layer conceptual model where the confining units are treated as non-insulating. The three-layered conceptual model is a realistic simplification of complex layered hydrogeologic systems if the hydrostratigraphic units above and below a confined aquifer can be lumped into a single layer with averaged hydraulic and electrical parameters. Whereas others have endeavoured to solve this problem numerically (e.g. Titov et al. 2005, using the finite difference method), analytical and semi-analytical approaches offer significant advantages, enumerated by Li & Neuman (2007), namely: the solution being representable in dimensionless form, rendering it general rather than site specific; revealing dimensionless parameters and space–time coordinates that control system behaviour, which may otherwise remain unidentified; obviating the need to construct computational grids and compute results across the entire grid at all times of interest and, generally, rendering parameter estimation easier, more stable and computationally efficient. Additionally, such solutions can be used to provide a benchmark for numerical models.

The solution was applied to field measurements of streaming potentials associated with the recovery period of an aquifer test reported in Rizzo et al. (2004), yielding average hydraulic conductivity values that compare well to those obtained by Rizzo et al. (2004). Estimates of specific storage were also obtained but with greater estimation uncertainty than estimates of hydraulic conductivity. In fact, the estimates of specific storage differed from those obtained by Rizzo et al. (2004) from direct head measurements by three orders of magnitude. This difference, coupled with the larger relative estimation variances indicates that specific storage may be a more difficult parameter to estimate using transient recovery streaming potential data. It should also be noted that the larger specific storage estimated from hydraulic head data may be due, in part, to the fact that hydraulic head data were not corrected for pumping and observation wellbore storage effects, and due to the fact that only recovery streaming potential data were used whereas head data for both the pumping and recovery phases were used in Rizzo et al. (2004). Using pumping phase self-potential data may improve the correspondence between estimates of specific storage from self-potential and hydraulic head data.

In addition to obtaining estimates of hydraulic conductivity and specific storage, we obtained estimates of the electrical conductivities of the upper and lower confining units that compare well to the values estimated by Rizzo et al. (2004) using electrical resistivity tomography. Since transient self-potential data used in parameter estimation are usually obtained at land surface and instrumentation is only minimally invasive, the solution developed in this work has the potential for rapidly yielding preliminary estimates of aquifer hydraulic properties where hydraulic head data from observation wells are unavailable.

2 Mathematical Formulation

We consider the streaming potentials that arise due to fluid flow toward a fully penetrating line sink in a confined aquifer of infinite radial extent. The governing equation for the fluid flow problem is

1

where s₁ = h₁(r, 0) −h₁(r, t) is drawdown (m), h₁ is hydraulic head (m), α = K₁/S_s,1 is hydraulic diffusivity of the porous medium (m² s⁻¹), K₁ is hydraulic conductivity (m s⁻¹), S_s,1 is specific storage (m⁻¹) and (r, t) are space–time coordinates. Eq. (1) is solved subject to the initial condition

2

the far-field boundary condition

3

and the pumping well (line sink) condition

4

where b₁ is the thickness of the confined aquifer (m) and Q is the pumping rate (m³ s⁻¹). The solution to this flow problem, due to Theis (1935), is

5

where s_D,1(x) = E₁(x) is the exponential integral (Abramowitz & Stegun 1972) and x = r²/(4αt).

For the streaming potential response, we consider a three-layer conceptual model consisting of an aquifer with more electrically conductive, but hydraulically impermeable confining units above and below, as shown in Fig. 1(b). The three-layered conceptual model is a realistic simplification of the more general hydrogeologic system depicted in Fig. 1(a), if the hydrostratigraphic units above and below a confined aquifer can be lumped into a single layer with averaged hydraulic and electrical parameters. The governing equation for the transient streaming potential response of the ith layer is (Revil et al. 2003)

6

where j_i is the electric current density (Am⁻²) and i = 1, 2, 3. It has been shown by many workers (e.g. Revil et al. 2003 and references therein) that

7

where σ_i is the electrical conductivity of the ith 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 = (γℓ_i/K_i) q_i is the electric current density due to fluid flow in one of the layers, γ is the specific weight of water (N m⁻³), ℓ_i is the streaming current coupling coefficient (m² V⁻¹ s⁻¹), and q_i = −K_i∇h_i = K_i∇s_i is the Darcy fluid flux (m s⁻¹).

Figure 1.

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

Substituting eq. (7) into eq. (6), in light of the radial flow assumption adopted above and the non-insulating nature of the upper and lower confining units, leads to

8

It should also be noted that for flow in an aquifer satisfying the solution of Theis (1935), a basic assumption adopted is that there is no fluid leakage from the confining units into the aquifer. This assumption further implies no fluid movement within the confining units, in which case the last term on the left-hand side of eq. (8) vanishes identically for i = 2, 3.

For the aquifer (i = 1)eq. (8) is solved subject to the initial condition

9

the far-field boundary condition

10

and the line sink condition

11

The condition given by eq. (11) defines the electrical sink/source due to groundwater extraction/injection through the pumping well. For the upper and lower confining units (i = 2, 3),eq. (8) is solved subject to the initial condition

12

the far-field boundary condition

13

the conditions for no pumping well at the centre of the confining units (since the pumping well is assumed to be screened only in the aquifer),

14

and

15

and the insulation boundary conditions

16

for the upper confining unit, and

17

for the lower confining unit, where b₂ is the z-coordinate of the upper boundary of the upper confining unit, and −b₃ is the z-coordinate of the lower boundary of the lower confining unit, see Fig. 1. Eqs (16) and (17) imply that the half-spaces above the upper unit and below the lower unit are insulating. This is based on the assumption that the half-space above the upper confining unit is the infinitely resistive atmosphere, and that below the lower confining unit is highly resistive unweathered bedrock.

Given that the confining units are non-insulating at their respective common boundaries with the aquifer, the following electrical potential and normal flux continuity conditions are imposed at these two boundaries:

18

19

20

21

3 Analytical Solution in Laplace–Hankel Transform Space

To solve the self-potential response problem described above, we first rewrite eq. (8) in dimensionless form as

22

where r_D = r/b₁, z_D = z/b₁, t_D = αt/b²₁, φ_D,i = φ_i/Φ_c, Φ_c = H_c(γℓ₁/σ₁), β_i = (ℓ_i/ℓ₁)/σ_D,i, with β₁ ≡ 1.0, and σ_D,i = σ_i/σ₁. For aquifer flow toward a pumping well, it is convenient to set H_c = Q/(4πb₁K₁). Under the condition of no fluid leakage from the confining units into the aquifer, the parameters β₂ and β₃ do not play a role in the self-potential response of the system to pumping, since, as discussed above, the third term on the left-hand side of eq. (22) vanishes identically. These parameters would only influence the self-potential response when fluid flow within the confining units cannot be neglected.

In dimensionless form, the initial and boundary conditions become

23

24

25

26

27

where b_D,i = b_i/b₁. In dimensionless form, continuity conditions at aquifer-confining layer boundaries become

28

29

30

31

Taking the Laplace and Hankel transforms (see Appendix A for definition of the latter) of eq. (22) and solving subject to the conditions given in eqs (23)–(31) leads to the following solutions for the Laplace–Hankel transforms of dimensionless electric potential in layers 1, 2 and 3:

32

where

33

34

35

36

37

38

39

40

41

42

43

where p and a are the Laplace and Hankel transform parameters, respectively. The inverse double Laplace–Hankel transform of the change in potential in both the aquifer and confining units due to pumping is

44

45

and

46

Eqs (45) and (46) are evaluated numerically. The computer programme, written in C++, is available from the authors upon request.

4 Model-Predicted Response

The predicted response in the aquifer for different values of z_D is shown in Fig. 2. The dimensionless parameter values used are σ_D,2 = σ_D,3 = 1 × 10³. The figure shows significant vertical variation in self-potential in the aquifer (an order of magnitude difference between z_D = 0 and z_D = 0.5) despite the fact that flow is entirely radial. This variation with z_D is attributable to charge inflow from the confining units. Charge inflow from the confining units also leads to steady-state late-time response of confined aquifer electric potential.

Figure 2.

Semi-log plot of the dimensionless streaming potential response of the aquifer, φ_D,1, against t_D/r²_D for different values of z_D with σ_D,2 = σ_D,3 = 10³.

Fig. 3 shows the predicted response in the upper confining unit, at three values of z_D, in both log–log and semi-log space. The semi-log plot shows that at late time, the slope of φ_D,2 is equal to the slope of the function u_D/(σ_D,2 + σ_D,3), where . It should be noted that, at late time, u_D ≈ −[ε + ln(x)]. Hence, in dimensional form one obtains

47

which would be useful for estimating the hydraulic conductivity of the aquifer. It should also be noted by comparing Figs 2 and 3 that, whereas the change in electric potential attains a steady-state in the confined aquifer, the same is not the case for that in the confining units. This is due to the fact that the charge flux into the aquifer from the confining units is balanced by the charge flux out of the aquifer through the pumping well. In contrast, there is no source of charge flux into the confining units to balance the outward flux into the aquifer.

Figure 3.

(a) Log–log and (b) semi-log plot of the dimensionless streaming potential response of the upper confining unit, φ_D,2, against t_D/r²_D for different values of z_D with σ_D,2 = σ_D,3 = 10³.

Fig. 4(a) is a plot of the dimensionless streaming potential response of the upper confining unit, φ_D,2, against dimensionless radial distance, r_D, at different values of dimensionless time, t_D. The figure shows the temporal evolution of the cone of potentials around the pumping well. The cones of the potentials closely mimic those of drawdown in the confined aquifer around the pumping well, as can be seen by comparing plots (a) and (b) of Fig. 4. This is to be expected since the extraction of water at the pumping well is the forcing function for the self-potential response.

Figure 4.

Plot of (a) the dimensionless streaming potential response of the upper confining unit, φ_D,2, and (b) the dimensionless drawdown in the aquifer, s_D,1, against dimensionless radial distance, r_D, at different values of dimensionless time, t_D.

For the case where a pump is operated from t_D = 0 to t_D = τ_D, the streaming potential response of the upper confining unit for both the pumping and recovery periods in the upper confining unit is given by

48

where φ_D,2(r_D, z_D, t_D − τ_D) ≡ 0 for t_D < τ_D. For large values of t_D − τ_D, eqs (48) and (47) lead to the following result

49

A solution of this form was used by Rizzo et al. (2004) in their analysis of recovery data; it is a special case of the more general analytical solution developed above. Fig. 5 shows the response predicted by eq. (48) for different values of τ_D.

Figure 5.

Log–log plot of the dimensionless streaming potential response of the upper confining unit, φ_D,2, against t_D/r²_D for different values of τ_D with σ_D,2 = σ_D,3 = 10³.

5 Application to Field Data

The model developed here was fitted to field data reported in Rizzo et al. (2004), which was obtained at a test site located near Montalto Uffugo, in the region of Calabria in Southern Italy. The aquifer at the site, which we treat as confined, is a silty sand layer extending from a mean depth of 11–55 m. It is bounded above by a shale formation that is overlain with heterogeneous gravels in a silty sand matrix. A shale substratum lies below the aquifer formation. A schematic of the subsurface at the test site showing the major hydrostratigraphic units is shown in Fig. 6. Electrical resistivity tomography results obtained by Rizzo et al. (2004) and reproduced here in Fig. 7, show that the different hydrostratigraphic units at the field site are not of uniform thicknesses. In the conceptual model used to develop the solution in this work, we assume that such units are of uniform thicknesses. Our solution should thus be understood to be an approximation of actual system behaviour. Additional details of the geology of the site, and on monitoring of the hydraulic and streaming potential responses, may be found in Rizzo et al. (2004).

Figure 6.

A schematic of the subsurface at test site near Montalto Uffugo, in the region of Calabria in Southern Italy, showing the major hydrostratigraphic units (after Rizzo et al. 2004).

Figure 7.

Electrical resistivity tomography at the test site (after Rizzo et al. 2004).

The experiment was conducted in 2003 July and involved pumping continuously at a constant rate of Q = 2.7 × 10⁻³ m³ s⁻¹ for a period of τ = 5855 min (∼4 d) from the confined aquifer. Streaming potential data were collected continuously, beginning 21 min before pump shut-off and continuing for several hours of the hydraulic recovery period. The layout of the electrodes used to collect the self-potential data is shown in Fig. 8. The pumping well (P5) and observation wells used to collect hydraulic head data are also shown. In fitting the model to observed data, we use the parameter values, reported in Rizzo et al. (2004), of 0.0915 S m⁻¹, 10 and 0.8 mV m⁻¹, respectively, for the electrical conductivity of water (σ_w), the formation factor (F) and the parameter γℓ₁/σ₁.

Figure 8.

Layout of the electrodes used to collect SP data at the test site. The pumping well (P5) and observation wells used to collect hydraulic head data are also shown (after Rizzo et al. 2004).

The non-linear parameter estimation software PEST (Doherty 2001) was used to jointly estimate the parameters K₁, S_s,1, σ₂ and σ₃ using the dimensional form of eq. (48) by minimizing the sum of squared residuals between observed and model-predicted self-potentials at each electrode. The noise in the data, as discussed in Rizzo et al. (2004), may be attributable to telluric currents and induction effects from a powerline crossing the field. Only multiples of the 50 Hz component of the noise were filtered out during data acquisition using a Fourier transform and low-pass filter. Despite the noise, Rizzo et al. (2004) showed that the decrease in the observed self-potential signal coincided with the drop in observed drawdown after cessation of pumping. This is an unambiguous indication that the measured self-potential signals are due to recovery of the aquifer piezometric surface. The parameter estimation results presented below were obtained by individually fitting the model to data collected at each electrode. Results for only six electrodes (three from each of the two lines shown in Fig. 8) are presented here for brevity. The six electrodes were selected to cover the domain of interest in an average sense.

5.1 Parameter estimation results

The estimated parameter values for data from electrodes 11, 13, 20, 35, 40 and 47 are given in Table 1. The mean values of K₁ and S_s,1, estimated from hydraulic head data by Rizzo et al. (2004), were 2.8 × 10⁻⁶ m s⁻¹ and 1.1 × 10⁻⁴ m⁻¹, respectively. The corresponding values estimated here using streaming potential data obtained with electrodes 11, 13, 20, 35, 40 and 47 are K = 2.2 × 10⁻⁶ m s⁻¹ and S_s = 4.7 × 10⁻⁷ m⁻¹. The estimated values of K given in the table are comparable to those obtained by Rizzo et al. (2004). Additionally, estimates of the electrical conductivities of the upper and lower confining units, σ₂ and σ₃, were obtained. They are also listed in Table 1. The average values of σ₂ = 5 × 10⁻² and σ₃ = 4 × 10⁻² S m⁻¹, compare well to the values reported in Rizzo et al. (2004, see Fig. 7) that were obtained using electrical resistivity tomography, where electrical resistivity is the reciprocal of electrical conductivity. Their results indicate that the confining units, which comprise shale and heterogeneous gravels in a silty sand matrix, have electrical conductivity values in the range 0.02–0.1 S m⁻¹.

Table 1.

Estimated parameter values at the indicated electrodes. r_e is the radial distance from the electrode to the pumping well.

Fig. 9 shows the fit of the solution given by eq. (48) to the measured potential change during recovery for electrodes listed in Table 1. The model fits to the data shown in the figure were obtained by individually fitting the model to data collected at each electrode. As can be seen in the figure, the solution fits the data well with relatively large coefficients of correlation (R² > 0.8). Table 2 gives the normalized estimation variances of the parameters listed in Table 1. The estimation variances are normalized by the estimated parameter values to allow for meaningful comparison of the estimation uncertainties of parameters whose values differ by orders of magnitude.

Figure 9.

Fit of model-predicted response to field data. The dashed lines are one standard deviation bounds on model fit (data after Rizzo et al. 2004).

Table 2.

Normalized estimation variances associated with parameter estimates reported in Table 1.

Table 3 gives a summary of the statistics of the residuals associated with the electrodes listed in Table 1. The statistics listed in the table are the sum of squared residuals , mean residual (μ_ε), variance of the residuals (σ²_ε) and the maximum residual [max(ε)]. The means of the residuals are on the order of 1 μV, which is about three orders of magnitude smaller than the measured electric potentials. In Fig. 9, the dashed lines represent bounds of one standard deviation on the fitted solution. Most of the self-potential measurements fall within these bounds. From Fig. 9 and the summary statistics in Table 3, it is clear that eq. (48) fits the data well.

Table 3.

Summary statistics of the residuals for the indicated electrodes.

6 Conclusion

The focus of this work was to present a semi-analytical solution to the problem of transient streaming potentials associated with pumping water from a confined aquifer. We adopted a three-layer conceptual model, consisting of a homogeneous aquifer and homogeneous impermeable confining units. In reality, the aquifer may be heterogeneous and the confining units may be multilayered and heterogeneous. For homogeneous multilayered confining units, the electrical properties of the units may be averaged to obtain the three-layered conceptual model used here. The solution indicates that the constant slope of the late-time surface self-potential data may be used to provide estimates of aquifer hydraulic conductivity, if estimates of the electrical properties of the aquifer and the confining units are available from other geophysical methods.

The solution was applied to field measurements obtained by Rizzo et al. (2004), yielding average values of 2.2 × 10⁻⁶ m s⁻¹ and 4.7 × 10⁻⁷ m⁻¹ for hydraulic conductivity and specific storage, respectively. The estimation standard deviations of these parameters are given in Table 1. Using hydraulic head data, Rizzo et al. (2004) estimated these parameter values to be 2.8 × 10⁻⁶ m s⁻¹ and 1.1 × 10⁻⁴ m⁻¹, respectively. Whereas using eq. (49), one can only estimate hydraulic conductivity from streaming potential data, as in Rizzo et al. (2004), we have demonstrated here that one can also obtain estimates of specific storage using the model developed in this work.

It should be noted, however, that the values of specific storage estimated from the transient recovery data of streaming potentials are significantly smaller than those obtained by Rizzo et al. (2004) from direct measurements of head using observation wells. The normalized estimation variances reported in Table 2 indicate that the estimation uncertainty associated with specific storage is larger than that associated with hydraulic conductivity. This may indicate that specific storage is a more difficult parameter to estimates with the analytical approach presented here. It is also worth noting that the larger values of specific storage estimated from hydraulic head data may partly be due to the fact that hydraulic head data were not corrected for pumping well and observation well storage. The observation and pumping wells were 20 cm in diameter, and hence, had significant wellbore storage that would lead one to overestimate aquifer specific storage. Mucha & Paulikova (1986) demonstrated that if the effect of wellbore storage is not accounted for in the forward model used in parameter estimation, it can lead to significant (one or more orders of magnitude depending on wellbore radius) overestimation of aquifer specific storage. This is due to the fact that wellbore storage leads to a lag in drawdown response, which for the classical Theis (1935) solution translates into high aquifer storage. Another reason for the mismatch may be that only recovery self-potential data were used whereas head data for both the pumping and recovery phases were used in Rizzo et al. (2004). Using pumping phase self-potential data may improve the correspondence between self-potential and hydraulic head based estimates of specific storage.

In addition to yielding estimates of hydraulic conductivity and specific storage, the solution developed in this work yielded estimates of the electrical conductivities of the upper and lower confining units that compare well to the values estimated by Rizzo et al. (2004) using electrical resistivity tomography. This demonstrates that one can in principle, estimate, not only the hydraulic conductivity, but also the specific storage of the aquifer, albeit with greater estimation uncertainty, and the electrical conductivities of the upper and lower confining units using only transient self-potential measurements. 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 hydraulic head data from observation wells are unavailable.

Acknowledgments

The work 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 Salvatore Straface and Enzo Rizzo for providing us with the data analysed in this work.

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 work 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 particular relation, adopted from Neuman & Witherspoon (1968), used in this work, is

(A3)

Appendix B: Nomenclature