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

Bwalya Malama, Kristopher L. Kuhlman, André Revil

► To cite this version:

Bwalya Malama, Kristopher L. Kuhlman, André Revil. Theory of transient streaming potentials associated with axial-symmetric flow in unconfined aquifers. Geophysical Journal International, 2009, 179 (2), pp.990-1003. 10.1111/j.1365-246X.2009.04336.x. insu-00498835

HAL Id: insu-00498835

https://insu.hal.science/insu-00498835

Submitted on 23 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Theory of transient streaming potentials associated with axial-symmetric flow in unconfined aquifers

Bwalya Malama,<fnref n="1" /> Kristopher L. Kuhlman<fnref n="2" /> and André Revil<fnref n="3" />,<fnref n="4" />

Accepted 2009 July 11. Received 2009 July 7; in original form 2009 April 22

## 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.

**Key words:** 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.

<footnote n="1">¹CGISS & Department of Geosciences, Boise State University, Boise, ID 83725, USA. E-mail: bmalama@cgiss.boisestate.edu</footnote>

<footnote n="2">²Repository Performance Department, Sandia National Laboratories, 4100 National Parks Highway, Carlsbad, NM 88220, USA</footnote>

<footnote n="3">³Colorado School of Mines, Department of Geophysics, Green Center, 1500 Illinois street, Golden, CO 80401, USA</footnote>

<footnote n="4">⁴Laboratoire de Geophysique Interne et Tectonophysique (LGIT), CNRS, Université de Savoie, Campus Scientifique, 73376 Le Bourget du Lac cedex, France</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1. Schematic of (a) the multilayered subsurface including an unconfined aquifer and (b) the three-layer conceptual model used to develop solution.</figcaption></figure>

## 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

$$ \frac{\partial s_D}{\partial t_D} = \frac{1}{r_D} \frac{\partial}{\partial r_D} \left( r_D \frac{\partial s_D}{\partial r_D} \right) + \kappa \frac{\partial^2 s_D}{\partial z_D^2}, \quad (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_2$, $z_D = z/b_2$, $t_D = \alpha_r t/b_2^2$, 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), $\kappa = K_z/K_r$ is the anisotropy ratio, $\alpha_r = K_r/S_s$ is the radial hydraulic diffusivity (m$^2$ s$^{-1}$), $K_r$ and $K_z$ are radial and vertical hydraulic conductivities (m s$^{-1}$), respectively, $S_s$ is the specific storage (m$^{-1}$), $b_2$ is the initial saturated thickness (m), $H_c = Q/(4\pi b_2 K_r)$ is the characteristic head of the system (m) and Q is the volumetric pumping rate (m$^3$ s$^{-1}$). Eq. (1) is solved subject to the initial condition

$$ s_D(r_D, z_D, t = 0) = 0, \quad (2) $$

the far-field boundary condition

$$ \lim_{r_D \to \infty} s_D(r_D, z_D, t_D) = 0, \quad (3) $$

the no-flow condition at the base of the aquifer

$$ \left. \frac{\partial s_D}{\partial z_D} \right|_{z_D=0} = 0, \quad (4) $$

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

$$ - \left. \frac{\partial s_D}{\partial z_D} \right|_{z_D=1} = \left. \frac{1}{\alpha_D} \frac{\partial s_D}{\partial t_D} \right|_{z_D=1}, \quad (5) $$

where $\alpha_D = \kappa/\vartheta$, $\vartheta = S_y/(b_2 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).

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

$$ \lim_{r \to 0} r \frac{\partial s}{\partial r} = - \frac{Q}{2\pi b_2 K_r}, \quad (6) $$

which, in non-dimensional form, becomes

$$ \lim_{r_D \to 0} r_D \frac{\partial s_D}{\partial r_D} = -2. \quad (7) $$

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

$$ \bar{s}_D^*(a, z_D, p) = \frac{2}{p(p + a^2)} \left[ 1 - \frac{\cosh(\eta z_D)}{\Delta_1} \right], \quad (8) $$
where $\eta^2 = (p + a^2)/\kappa$, p and a are the Laplace (in $t_D$) and Hankel (in $r_D$) transform parameters, respectively, and

$$ \Delta_1 = \cosh(\eta) + \frac{\eta\alpha_D}{p} \sinh(\eta). \qquad (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)

$$ \nabla \cdot \mathbf{j}_i = 0, \qquad (10) $$

where $\mathbf{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

$$ \mathbf{j}_i = \sigma_i \mathbf{E}_i + \mathbf{j}_{s,i}, \qquad (11) $$

where $\sigma_i$ is the electrical conductivity of the i th layer (S m⁻¹), $\mathbf{E}_i = -\nabla\phi_i$ is the electric field (V m⁻¹), $\phi_i = \phi_i - \phi_{0,i}$ is the electric potential change (V) in i th layer due to pumping in one of the layers, $\phi_{0,i}$ is the potential at $t = 0$, $\mathbf{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)

$$ \mathbf{j}_{s,i} = \gamma \ell_i \mathbf{K}_i^{-1} \mathbf{q}_i, \qquad (12) $$

where $\gamma$ is the specific weight of water (N m⁻³), $\ell_i = -\sigma_i C'_i / \gamma$ is the electrokinetic coupling term [m²(V s⁻¹)], $C'_i = \partial\phi_i / \partial h$ is the streaming potential coupling coefficient (V m⁻¹), $\mathbf{K}_i$ is the hydraulic conductivity tensor and $\mathbf{q}_i = -\mathbf{K}_i \nabla h_i = \mathbf{K}_i \nabla 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

$$ \frac{1}{r_D} \frac{\partial}{\partial r_D} \left( r_D \frac{\partial \phi_{D,2}}{\partial r_D} \right) + \frac{\partial^2 \phi_{D,2}}{\partial z_D^2} - \left[ \frac{1}{r_D} \frac{\partial}{\partial r_D} \left( r_D \frac{\partial s_D}{\partial r_D} \right) + \frac{\partial^2 s_D}{\partial z_D^2} \right] = 0, \qquad (13) $$

and for $i = 1, 3$:

$$ \frac{1}{r_D} \frac{\partial}{\partial r_D} \left( r_D \frac{\partial \phi_{D,i}}{\partial r_D} \right) + \frac{\partial^2 \phi_{D,i}}{\partial z_D^2} = 0, \qquad (14) $$

where, for $i = 1, 2, 3$, $\phi_{D,i} = \phi_i / \Phi_c$ and $\Phi_c = (\gamma \ell_2 / \sigma_2) 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

$$ \phi_{D,i}(r_D, z_D, t_D = 0) = 0, \qquad (15) $$

and the boundary conditions

$$ \lim_{r_D \to \infty} \phi_{D,i}(r_D, z_D, t_D) = 0, \qquad (16) $$

at far-field in all layers,

$$ \left. \frac{\partial \phi_{D,1}}{\partial z_D} \right|_{z_D=1+b_{D,1}} = 0, \qquad (17) $$

(insulating condition) at ground surface,

$$ \left. \frac{\partial \phi_{D,3}}{\partial z_D} \right|_{z_D=-b_{D,3}} = 0, \qquad (18) $$

(insulating condition) at the base of layer 3. Here $b_{D,1} = b_1 / b_2$ and $b_{D,3} = b_3 / b_2$ 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

$$ \lim_{r \to 0} r \frac{\partial \phi_i}{\partial r} = \begin{cases} 0 & i = 1, 3 \\ -\frac{Q}{2\pi b_2 K_r} \frac{\gamma \ell_2}{\sigma_2} & i = 2 \end{cases}, \qquad (19) $$

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

$$ \lim_{r_D \to 0} r_D \frac{\partial \phi_{D,i}}{\partial r_D} = \begin{cases} 0 & i = 1, 3 \\ -2 & i = 2 \end{cases}. \qquad (20) $$

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

$$ \phi_{D,1}(r_D, z_D = 1, t_D) = \phi_{D,2}(r_D, z_D = 1, t_D), \qquad (21) $$

$$ \phi_{D,3}(r_D, z_D = 0, t_D) = \phi_{D,2}(r_D, z_D = 0, t_D), \qquad (22) $$
$$ \sigma_{D,1} \left. \frac{\partial \phi_{D,1}}{\partial z_D} \right|_{z_D=1} = \left. \frac{\partial \phi_{D,2}}{\partial z_D} \right|_{z_D=1}, \qquad (23) $$

$$ \sigma_{D,3} \left. \frac{\partial \phi_{D,3}}{\partial z_D} \right|_{z_D=0} = \left. \frac{\partial \phi_{D,2}}{\partial z_D} \right|_{z_D=0}. \qquad (24) $$

where $\sigma_{D,i} = \sigma_i / \sigma_2$, 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

$$ \bar{\phi}_{D,i}^*(a, z_D, p) = \begin{cases} \left( \frac{1}{\Delta_2} [\xi_1 + \xi_2 f_1(a)] - \frac{\bar{s}_D^*(a, 1, p)}{\cosh(ab_{D,1})} \right) \cosh[a(1 + b_{D,1} - z_D)] & i = 1 \\ \frac{1}{\Delta_2} [\xi_1 g_1(a, z_D) + \xi_2 g_2(a, z_D)] - \bar{s}_D^*(a, z_D, p) & i = 2, \\ \frac{1}{\Delta_2} [\xi_2 - \bar{s}_D^*(a, 0, p) f_2(a)] \cosh[a(b_{D,3} + z_D)] & i = 3 \end{cases} \qquad (25) $$

where

$$ \xi_1 = \sigma_{D,3} \bar{s}_D^*(a, 0, p) \sinh(ab_{D,3}) \qquad (26) $$

$$ \xi_2 = \bar{s}_D^*(a, 1, p)[\sigma_{D,1} \sinh(ab_{D,1}) - (p/a) \cosh(ab_{D,1})] \qquad (27) $$

$$ f_1(a) = [\cosh(ab_{D,3}) \cosh(a) + \sigma_{D,3} \sinh(ab_{D,3}) \sinh(a)] / \cosh(ab_{D,1}) \qquad (28) $$

$$ f_2(a) = \sigma_{D,1} \cosh(a) \sinh(ab_{D,1}) + \sinh(a) \cosh(ab_{D,1}) \qquad (29) $$

$$ g_1(a, z_D) = \cosh(ab_{D,1}) \cosh[a(1 - z_D)] + \sigma_{D,1} \sinh(ab_{D,1}) \sinh[a(1 - z_D)] \qquad (30) $$

$$ g_2(a, z_D) = \cosh(ab_{D,3}) \cosh(az_D) + \sigma_{D,3} \sinh(ab_{D,3}) \sinh(az_D) \qquad (31) $$

$$ \Delta_2 = |A_1 \sinh(a) + A_2 \cosh(a)| \qquad (32) $$

$$ A_1 = \cosh(ab_{D,1}) \cosh(ab_{D,3}) + \sigma_{D,1} \sigma_{D,3} \sinh(ab_{D,1}) \sinh(ab_{D,3}) \qquad (33) $$

$$ A_2 = \sigma_{D,1} \sinh(ab_{D,1}) \cosh(ab_{D,3}) + \sigma_{D,3} \cosh(ab_{D,1}) \sinh(ab_{D,3}). \qquad (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) $\phi_{D,1}$ and (b) $\phi_{D,2}$ are plotted against $t_D/r_D^2$ for different vertical positions. The results were computed with dimensionless parameters $\kappa = 1.0$, $\vartheta = 31.3$, $\sigma_{D,1} = 0.05$, $\sigma_{D,3} = 40.0$, $b_{D,1} = 0.25$ and

<figure><img src="image_3.png" /><figcaption>Figure 2. Plot of dimensionless streaming potentials $\phi_{D,i}$ for layers (a) $i = 1$ and (b) $i = 2$, against $t_D/r_D^2$, for different values of $z_D$.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Plot of dimensionless streaming potentials, $\phi_{D,i}$, versus $t_D/r_D^2$, for different values of the dimensionless storage parameter $\vartheta = S_y/(b_2S_s)$. (a) $i = 1, z_D = 1 + b_{D,1}$ and (b) $i = 2, z_D = 0.5$.</figcaption></figure>

$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.

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 $\vartheta = S_y/(b_2S_s)$, using the parameters $\kappa = 1.0, \sigma_{D,1} = 0.05, \sigma_{D,3} = 40.0, b_{D,1} = 0.25$ and $b_{D,3} = 0.5$. The parameter value $\vartheta = 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_2S_s$). This implies that low streaming potential signals are to be expected when pumping tests are conducted in media with small specific yield values.

The effect of the hydraulic conductivity anisotropy ratio, $\kappa = K_z/K_r$, is shown in Fig. 4(a). The model results shown in the figure indicate that $\kappa$ 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 $\kappa$. The plots in the figure show that high potentials and drawdowns are associated with small values of $\kappa$ at early and intermediate times. This is due to the fact that for values of $\kappa < 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 $\kappa > 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

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Plot of dimensionless (a) streaming potentials at land surface, $\phi_{D,i}(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^2$, for different values of the hydraulic conductivity anisotropy ratio $\kappa = K_z/K_r$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Plot of (a) dimensionless streaming potentials at ground surface, $\phi_{D,1}(r_D, z_D = 1 + b_{D,1}, t_D)$ and (b) the depth averaged unconfined aquifer drawdown, $\langle s_D(r_D, t_D) \rangle$, against $r_D$, at different values of the dimensionless time, $t_D$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Plot of dimensionless streaming potentials, $\phi_{D,i}$, versus $t_D/r_D^2$, for different values of the dimensionless thickness, $b_{D,1}$, of the zone above the water-table.</figcaption></figure>

aquifer is predominantly radial, both drawdown and streaming potential responses merge to the response controlled by the radial hydraulic conductivity.

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 \to 0$. This is due to the symmetry boundary condition in eq. (19) imposed at $r_D = 0$ for layer $i = 1$.

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 $\vartheta = 31.3$, $\kappa = 1.0$, $\sigma_{D,1} = 0.05$, $\sigma_{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.

## 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
<figure><img src="image_8.png" /><figcaption>**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_{\text{obs}}$, $y_{\text{obs}}$) are coordinates of a measurement location.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**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.</figcaption></figure>

compute the response at the point $\mathbf{x}_{\text{obs}} = (x_{\text{obs}}, y_{\text{obs}}, z_{\text{obs}})$ according to

$$ \phi_{D,\text{obs}}^{(i)}(\mathbf{x}_{\text{obs}}, t) = \phi_{D,i}(r, z_{\text{obs}}, t) \pm \phi_{D,i}(r^*, z_{\text{obs}}, t), \quad (35) $$

where $\phi_{D,i} = \mathcal{H}_0^{-1}\mathcal{L}^{-1}\{\phi_{D,i}^*\}$ 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_2$, 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.

# 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><img src="image_10.png" /><figcaption>**Figure 9.** Aerial photo of the BHRS and existing wells (inset) at the site (after Barrash & Reboulet 2004).</figcaption></figure>

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, φ₀ᵢ, 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.

## 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ᵣ, Kₛ, Sₛ and Sᵣ. 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
<figure><img src="image_11.png" /><figcaption>**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).</figcaption></figure>

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.

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).

## 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^2$-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 \times 10^{-4} \text{ m s}^{-1}$ and $K_z = 3 \times 10^{-4} \text{ m s}^{-1}$, and Table 2 gives average values of $K_r = 4.6 \times 10^{-4} \text{ m s}^{-1}$ and $K_z = 6.5 \times 10^{-4} \text{ m s}^{-1}$. The values estimated here from streaming potential data average $K_r = 2 \times 10^{-4} \text{ m s}^{-1}$ and $K_z = 3 \times 10^{-4} \text{ m s}^{-1}$. 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 \times 10^{-4} \text{ m}^{-1}$, and results in Table 2 have an average value of $S_s = 2.6 \times 10^{-4} \text{ m}^{-1}$. The average value estimated here from streaming potential data is of the order of $1 \times 10^{-3} \text{ m}^{-1}$.

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
<figure><img src="image_12.png" /><figcaption>**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.</figcaption></figure>

<table><caption>**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).</caption><thead><tr><th>Electrode</th><th>r' (m)</th><th>K<sub>r</sub> (×10<sup>-4</sup> m s<sup>-1</sup>)</th><th>S<sub>s</sub> (×10<sup>-3</sup> m<sup>-1</sup>)</th><th>S<sub>y</sub></th><th>κ</th></tr></thead><tbody><tr><td>3</td><td>1.90</td><td>2.14</td><td>38.2</td><td>0.103</td><td>0.13</td></tr><tr><td>4</td><td>2.92</td><td>1.30</td><td>1.66</td><td>0.226</td><td>0.21</td></tr><tr><td>5</td><td>5.26</td><td>1.06</td><td>6.70</td><td>0.028</td><td>0.08</td></tr><tr><td>6</td><td>7.60</td><td>1.02</td><td>1.89</td><td>0.007</td><td>3.74</td></tr><tr><td>12</td><td>1.24</td><td>2.74</td><td>14.2</td><td>0.007</td><td>2.08</td></tr><tr><td>13</td><td>2.43</td><td>1.17</td><td>2.69</td><td>0.013</td><td>6.81</td></tr></tbody></table>

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
<table><caption>**Table 2.** Hydraulic parameter values estimated with the model of Neuman (1972).</caption><thead><tr><th>Well</th><th>r<br>(m)</th><th>K<sub>r</sub><br>(10<sup>-4</sup> m s<sup>-1</sup>)</th><th>S<sub>y</sub><br>(10<sup>-4</sup> m<sup>-1</sup>)</th><th>κ</th><th>S<sub>y</sub></th></tr></thead><tbody><tr><td>B2</td><td>3.51</td><td>4.64</td><td>2.22</td><td>1.31</td><td>0.052</td></tr><tr><td>B4</td><td>2.62</td><td>4.61</td><td>3.48</td><td>2.02</td><td>0.090</td></tr><tr><td>B5</td><td>3.65</td><td>4.64</td><td>2.16</td><td>1.25</td><td>0.062</td></tr><tr><td>B6</td><td>2.69</td><td>4.57</td><td>2.52</td><td>1.08</td><td>0.080</td></tr></tbody></table>

<table><caption>**Table 3.** Composite sensitivities of simulated responses to estimated parameters at optimality.</caption><thead><tr><th>Electrode</th><th>sen (K<sub>r</sub>)<br>(×10<sup>-4</sup>)</th><th>sen (S<sub>s</sub>)<br>(×10<sup>-4</sup>)</th><th>sen (S<sub>y</sub>)<br>(×10<sup>-4</sup>)</th><th>sen (κ)<br>(×10<sup>-4</sup>)</th></tr></thead><tbody><tr><td>3</td><td>2.30</td><td>0.92</td><td>0.77</td><td>1.24</td></tr><tr><td>4</td><td>9.74</td><td>1.50</td><td>0.80</td><td>4.73</td></tr><tr><td>5</td><td>6.39</td><td>2.02</td><td>1.93</td><td>3.66</td></tr><tr><td>6</td><td>2.86</td><td>0.94</td><td>0.19</td><td>0.17</td></tr><tr><td>12</td><td>1.78</td><td>0.53</td><td>0.02</td><td>0.05</td></tr><tr><td>13</td><td>3.86</td><td>0.80</td><td>0.24</td><td>0.09</td></tr></tbody></table>

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

<table><caption>**Table 4.** Normalized estimation variances for parameters estimated from streaming potential data.</caption><thead><tr><th>Electrode</th><th>σ̂<sup>2</sup><sub>K<sub>r</sub></sub></th><th>σ̂<sup>2</sup><sub>S<sub>s</sub></sub></th><th>σ̂<sup>2</sup><sub>S<sub>y</sub></sub></th><th>σ̂<sup>2</sup><sub>κ</sub></th></tr></thead><tbody><tr><td>3</td><td>1.32 × 10<sup>-1</sup></td><td>1.01 × 10<sup>1</sup></td><td>1.07 × 10<sup>1</sup></td><td>1.66 × 10<sup>1</sup></td></tr><tr><td>4</td><td>4.47 × 10<sup>-1</sup></td><td>2.79 × 10<sup>-1</sup></td><td>8.65 × 10<sup>-1</sup></td><td>1.60 × 10<sup>0</sup></td></tr><tr><td>5</td><td>2.47 × 10<sup>-1</sup></td><td>8.09 × 10<sup>-2</sup></td><td>1.18 × 10<sup>0</sup></td><td>2.80 × 10<sup>-1</sup></td></tr><tr><td>6</td><td>1.09 × 10<sup>-4</sup></td><td>1.54 × 10<sup>-1</sup></td><td>3.00 × 10<sup>0</sup></td><td>2.56 × 10<sup>0</sup></td></tr><tr><td>12</td><td>6.85 × 10<sup>-4</sup></td><td>1.45 × 10<sup>-3</sup></td><td>5.94 × 10<sup>1</sup></td><td>1.62 × 10<sup>1</sup></td></tr><tr><td>13</td><td>3.99 × 10<sup>-4</sup></td><td>6.98 × 10<sup>0</sup></td><td>6.79 × 10<sup>1</sup></td><td>3.04 × 10<sup>1</sup></td></tr></tbody></table>
<figure><img src="image_13.png" /><figcaption>**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⁻¹.</figcaption></figure>

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.

REFERENCES

Bailey, D.H., Jeyabalan, K., Li, X.S., 2000. A comparison of three high-precision quadrature schemes, *Exp. Math.*, **14**(3), 317–329.

Barrash, W. & Clemo, T., 2002. Hierarchical geostatistics and multifacies systems: Boise hydrogeophysical research site, Boise, Idaho *Water Resour. Res.*, **38**(10), 1196, doi:10.1029/2002WR001436.

Barrash, W. & Reboulet, E., 2004. Significance of porosity for stratigraphy and textural composition in subsurface, coarse fluvial deposits: Boise hydrogeophysical research site, *Geol. soc. Am. Bull.*, **116**, 1059–1073, doi:10.1130/B25370.1.

Barrash, W., Morin, R.H. & Gallegos, D.M., 1997. Lithologic, hydrologic, and petrophysical characterization of an unconsolidated cobble-and-sand aquifer, capital station site, Boise, Idaho, in *Proceedings of the 32nd Symposium on Engineering Geology and Geotechnical Engineering*. pp. 307–323, Idaho State University, Pocatello, ID.

Barrash, W., Clemo, T., Fox, J.J. & Johnson, T.C., 2006. Field, laboratory, and modeling investigation of the skin effect at wells with slotted casing, *J. Hydrol.*, **326**(1–4), 181–198, doi:10.1016/j.jhydrol.2005.10.029.

Bear, J., 1972. *Dynamics of Fluids in Porous Media*, Dover, New York.

Bogoslovsky, V.A. & Ogilvy, A.A., 1973. Deformations of natural electric fields near drainage structures. *Geophys. Prospect.*, **21**, 716–723.

Bolève, A., Crespy, A., Revil, A., Janod, J. & Mattiuzzo, J.L., 2007. Streaming potentials of granular media: influence of the Dukhin and Reynolds numbers, *J. geophys. Res.*, **112**, B08204, doi:10.1029/2006JB004673.

Cardiff, M., Barrash, W., Kitanidis, P.K., Malama, B., Revil, A., Straface, S. & Rizzo, E., 2009. A potential-based inversion of unconfined steady-state hydraulic tomography, *Ground Water*, **47**(2), 259–270.

Carslaw, H.S. & Jaeger J.C., 1959. *Conduction of Heat in Solids*, 2nd edn, Oxford University Press, New York, pp. 273.

Chandler, R., 1981. Transient streaming potential measurements on fluid-saturated porous structures: an experimental verification of Biot's slow wave in the quasi-static limit, *J. acoust. Soc. Am.*, **70**(1), 116–121.

Clement, W.P., Knoll, M.D., Liberty, L.M., Donaldson, P.R., Michaels, P.M., Barrash, W. & Pelton, J.R., 1999. Geophysical surveys across the Boise Hydrogeophysical Research Site to determine geophysical parameters of a shallow alluvial aquifer, in *Proceedings of SAGEEP99: Application of Geophysics to Engineering and Environmental Problems*, Oakland, CA, March, pp. 399–408.

Darnet, M., Marquis, G. & Sailhac, P., 2003. Estimating aquifer hydraulic properties from the inversion of surface Streaming Potential (SP) anomalies, *Geophys. Res. Lett.* **30**(13), 1679–1682, doi:10.1029/2003GL017631.

Doherty, J., 2002. *Manual for PEST*, 5th edn, Watermark Numerical Computing, Australia.

Jardani, A. et al., 2009. Reconstruction of the water table from self-potential data: a Bayesian approach, *Ground Water*, **47**(2), 213–227.

Leroy, P., Revil, A., Altmann, S. & Tournassat, C., 2007. Modeling the composition of the pore water in a clay-rock geological formation (Callovo-Oxfordian, France), *Geochim. Cosmochim. Acta*, **71**(5), 1087–1097, doi:10.1016/j.gca.2006.11.009.

Kresic, N., 2007. *Hydrogeology and Groundwater Modeling*, 2nd edn, CRC Press, New York, pp. 192.

Maineult, A., Strobach, E. & Renner, J., 2008. Self-potential signals induced by periodic pumping tests, *J. geophys. Res.*, **113**, B01203, doi:10.1029/2007JB005193.

Malama, B., Revil, A. & Kuhlman, K.L., 2009. A semi-analytical solution for transient streaming potentials associated with confined aquifer pumping tests, *Geophys. J. Int.*, **176**, 1007–1016, doi:10.1111/j.1365-246X.2008.04014.x.

Moench, A.F., 1994. Specific yield as determined by type-curve analysis of aquifer-test data, *Ground Water*, **32**(6), 949–947.

Neuman, S.P., 1972. Theory of flow in unconfined aquifers considering delayed response of the water table, *Water Resour. Res.*, **8**(4), 1031–1045.

Neuman, S.P. & Witherspoon, P.A., 1968. Theory of flow in aquicludes adjacent to slightly leaky aquifers, *Water Resour. Res.*, **4**(1), 103–112.

Nwankwor, G.I., Cherry, J.A. & Gillham, R.W., 1984. A comparative study of specific yield determinations for a shallow sand aquifer, *Ground Water*, **22**(6), 764–772.

Nwankwor, G.I., Gillham, R.W., van der Kamp, G. & Akindunni, F.F., 1992. Unsaturated and saturated flow in response to pumping of an unconfined aquifer: field evidence of delayed drainage, *Ground Water*, **30**(5), 690–700.

Othberg, K.L., 1994. Geology and geomorphology of the Boise valley and adjoining areas, western snake river plain, Idaho, *Idaho Geol. Surv. Bull.* **29**, 1–54.

Petiau, G., 2000. Second generation of lead–lead chloride electrodes for geophysical applications, *Pageoph*, **157**, 357–382.

Press, W.H., Teukolsky, S.A., Vetterling, W.T. & Flannery, B.P., 2007. *Numerical Recipies*, 3rd edn, Cambridge University Press.

Revil, A., Naudet, V., Nouzaret, J. & Pessel, M., 2003. Principles of electrography applied to self-potential electrokinetic sources and hydrogeological applications, *Water Resour. Res.*, **39**(5), 1114, doi:10.1029/2001WR000916.

Rizzo, E., Suski, B. & Revil, A., 2004. Self-potential signals associated with pumping tests experiments, *J. geophys. Res.*, **109**, B10203, doi:10.1029/2004/JB003049.

Sailhac, P. & Marquis, G., 2001. Analytic potential for the forward and inverse modeling of SP anomalies caused by subsurface fluid flow, *Geophys. Res. Lett.*, **28**(9), 1851–1854.

Sanchez-Vila, X., Tartakovsky, D.M., 2007. Ergodicity of pumping tests, *Water Resour. Res.*, **43**, W03414, doi:10.1029/2006WR005241.

Sill, W.R., 1983. Self-potential modeling from primary flows, *Geophysics*, **48**(1), 76–86.

Squires, E., Wood, S.H. & Osiensky, J.L., 1992. Hydrogeologic framework of the Boise aquifer system, Ada county, Idaho, Technical Report, Idaho Water Resources Research Institute, Moscow, Idaho.

Stehfest, H., 1970. Algorithm 368: numerical inversion of Laplace transforms, *Commun. ACM*, **13**(1), 47–49.

Suski, B., Rizzo, E. & Revil, A., 2004. A sandbox experiment of self-potential signals associated with a pumping test, *Vadose Zone J.*, **3**, 1193–1199.

Suski, B., Revil, A., Titov, Konosavsky, Voltz, Dangès & Huttel 2006. Monitoring of an infiltration experiment using the self-potential method, *Water Resour. Res.*, **42**, W08418, doi:10.1029/2005WR004840.

Takahasi, H. & Mori, M., 1974. Double exponential formulas for numerical integration, *Publ. Res. Inst. Math. Sci.*, **9**, 721–741.

Tartakovsky, G.D. & Neuman, S.P., 2007. Three-dimensional saturated-unsaturated flow with axial symmetry to a partially penetrating well in a compressible unconfined aquifer, *Water Resour. Res.*, **43**, 1–17.

Titov, K., Ilyin, Y., Konosavsky, P. & Levitski, A., 2002. Electrokinetic spontaneous polarization in porous media: petrophysics and numerical modeling, *J. Hydrol.*, **267**, 207–216.

Titov, K., Revil, A., Konosavsky, P., Straface, S. & Troisi, S., 2005. Numerical modeling of self-potential signals associated with a pumping test experiment, *Geophys. J. Int.*, **162**, 641–650.

Wieder, T., 1999. Algorithm 794: numerical Hankel Transform by the Fortran Program HANKEL, *ACM Trans. Math. Software*, **25**(2), 240–250.
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

$$ \mathcal{H}_0\{f(r_D)\} = f^*(a) = \int_0^\infty r_D J_0(ar_D) f(r_D) dr_D, \qquad (A1) $$

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

$$ \mathcal{H}_0^{-1}\{f^*(a)\} = f(r_D) = \int_0^\infty a J_0(ar_D) f^*(a) da. \qquad (A2) $$

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

$$ \mathcal{H}_0\left\{\frac{1}{r_D} \frac{\partial}{\partial r_D} \left(r_D \frac{\partial f}{\partial r_D}\right)\right\} = -a^2 f^* - \lim_{r_D \to 0} r_D \frac{\partial f}{\partial r_D}. \qquad (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 evaluated in two steps (Wieder 1999). One finite ($0 \le a \le j_{0,1}$), the other approaching infinite in length ($j_{0,1} \le a \le j_{0,m}$), where $j_{0,n}$ is the nth zero of $J_0(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 evaluations at lower levels $k-j$ are simply the $2^{k-j}$th 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 \to 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 $\epsilon$-algorithm.