A physically-based analytical model to describe effective excess charge for streaming potential generation in water saturated porous media

TL;DR

This paper introduces a physically-based analytical model for predicting effective excess charge in saturated porous media, improving understanding of streaming potential generation for water flow estimation in hydrogeology.

Contribution

It presents a novel flux-averaging model based on a bundle of capillary tubes with fractal pore-size distribution, linking electrokinetic effects to hydrogeological properties.

Findings

Successfully tested against experimental data

Explains empirical relation between excess charge and permeability

Reproduces charge evolution with electrolyte concentration

Abstract

Among the different contributions generating self-potential, the streaming potential is of particular interest in hydrogeology for its sensitivity to water flow. Estimating water flux in porous media using streaming potential data relies on our capacity to understand, model, and upscale the electrokinetic coupling at the mineral-solution interface. Different approaches have been proposed to predict streaming potential generation in porous media. One of these approaches is the flux averaging which is based on determining the excess charge which is effectively dragged in the medium by water flow. In this study, we develop a physically-based analytical model to predict the effective excess charge in saturated porous media using a flux-averaging approach in a bundle of capillary tubes with a fractal pore-size distribution. The proposed model allows the determination of the effective excess charge as a function of pore water ionic concentration and hydrogeological parameters like porosity, permeability and tortuosity. The new model has been successfully tested against different set of experimental data from the literature. One of the main findings of this study is the mechanistic explanation to the empirical dependence between the effective excess charge and the permeability that has been found by several researchers. The proposed model also highlights the link to other lithological properties and it is able to reproduce the evolution of effective excess charge with electrolyte concentrations.

Figures

Tables

Table 1. Table 1: Comparison between fitted and predicted hydraulic tortuosities

Material	$\tau$ from fit	${\tau }_e$ from Eq. (30)
Glass beads packs from Glover and Déry [2010]
Mean value over 12 packs	1.20	1.24
Samples from Pengra et al. [1999]
Fontainebleau A	1.95	1.59
Fontainebleau B	1.83	1.82
Fontainebleau C	3.24	3.03
Glass beads (200 $\mu$m)	1.90	1.64

Equations

$$\nabla \cdot \left(\sigma \nabla \phi \right)=\nabla \cdot \left({\hat{Q}}_v{\text{v}}_D\right)$$

$$N(R)={\left(\frac{R_{REV}}{R}\right)}^D,$$

$$dN=D{R}_{REV}^D{R}^{-D-1}dR.$$

$$\varphi =\frac{\text{Volume of pores}}{\text{Volume of REV}}=\frac{{\int }_{R_{min}}^{R_{max}}{V}_p(R)dN}{\pi R_{REV}^{2}L},$$

$$\varphi =\frac{\tau D}{R_{REV}^{2-D}(2-D)}\left(R_{max}^{2-D}-R_{min}^{2-D}\right),$$

$$v(R,r)=\frac{{\rho }_{w}g}{4\eta \tau }\left[R^2-(R-r{)}^2\right]\frac{\Delta h}{L},$$

$$\overline{v}(R)=\frac{{\rho }_{w}g}{8\eta \tau }{R}^2\frac{\Delta h}{L}.$$

$$V_Q={\int }_{R_{min}}^{R_{max}}\overline{v}(R)\pi R^{2}dN=\frac{{\rho }_{w}g}{8\eta \tau }\frac{\pi D{R}_{REV}^D}{4-D}\left(R_{max}^{4-D}-R_{min}^{4-D}\right)\frac{\Delta h}{L}.$$

$$V_Q=\pi R_{REV}^2\frac{{\rho }_{w}g}{\eta }k\frac{\Delta h}{L},$$

$$k=\frac{D}{8\tau (4-D)R_{REV}^{2-D}}\left(R_{max}^{4-D}-R_{min}^{4-D}\right).$$

$$\varphi =\frac{\tau D}{(2-D)R_{REV}^{2-D}}{R}_{max}^{2-D},$$

$$k=\frac{D}{8\tau (4-D)R_{REV}^{2-D}}{R}_{max}^{4-D}.$$

$$k=\gamma {\varphi }^{\frac{4-D}{2-D}},$$

$${\bar{Q}}_v(r)=N_A{e}_0{C}^0\left[e^{-\frac{e_0\psi (r)}{k_{B}T}}-e^{\frac{e_0\psi (r)}{k_{B}T}}\right],$$

$$\psi (r)=\zeta e^{-\frac{r}{l_D}},$$

$$l_D=\sqrt{\frac{ϵk_{B}T}{2N_A{C}^0{e}_0^2}},$$

$${\hat{Q}}_v^R=\frac{1}{\overline{v}(R)\pi R^2}{\int }_A{\overline{Q}}_v(r)v(R,r)dA,$$

$${\hat{Q}}_v^R=\frac{2}{\overline{v}(R)R^2}{\int }_0^R{\overline{Q}}_v(r)v(R,r)(R-r)dr.$$

$$e^{\pm \frac{e_0\psi (r)}{k_{B}T}}=1\pm \frac{e_0\psi (r)}{k_{B}T}+\frac{1}{2}{\left(\frac{e_0\psi (r)}{k_{B}T}\right)}^2\pm \frac{1}{6}{\left(\frac{e_0\psi (r)}{k_{B}T}\right)}^3.$$

$$\begin{array}[]{ll}\hat{Q}^{R}_{v}=&-\frac{8N_{A}e_{0}^{2}C^{0}\zeta}{k_{B}T(R/l_{D})^{4}}\left\{6-e^{-\frac{R}{l_{D}}}\left[\left(\frac{R}{l_{D}}\right)^{3}+3\left(\frac{R}{l_{D}}\right)^{2}+6\left(\frac{R}{l_{D}}\right)+6\right]\right\}\\ &+\frac{24N_{A}e_{0}^{2}C^{0}\zeta}{k_{B}T(R/l_{D})^{3}}\left\{2-e^{-\frac{R}{l_{D}}}\left[\left(\frac{R}{l_{D}}\right)^{2}+2\left(\frac{R}{l_{D}}\right)+2\right]\right\}\\ &-\frac{16N_{A}e_{0}^{2}C^{0}\zeta}{k_{B}T(R/l_{D})^{2}}\left\{1-e^{-\frac{R}{l_{D}}}\left[\left(\frac{R}{l_{D}}\right)+1\right]\right\}\\ &-\frac{4N_{A}e_{0}^{4}C^{0}\zeta^{3}}{3(k_{B}T)^{3}(3R/l_{D})^{4}}\left\{6-e^{-\frac{3R}{l_{D}}}\left[\left(\frac{3R}{l_{D}}\right)^{3}+3\left(\frac{3R}{l_{D}}\right)^{2}+6\left(\frac{3R}{l_{D}}\right)+6\right]\right\}\\ &+\frac{4N_{A}e_{0}^{4}C^{0}\zeta^{3}}{(k_{B}T)^{3}(3R/l_{D})^{3}}\left\{2-e^{-\frac{3R}{l_{D}}}\left[\left(\frac{3R}{l_{D}}\right)^{2}+2\left(\frac{3R}{l_{D}}\right)+2\right]\right\}\\ &-\frac{8N_{A}e_{0}^{4}C^{0}\zeta^{3}}{3(k_{B}T)^{3}(3R/l_{D})^{2}}\left\{1-e^{-\frac{3R}{l_{D}}}\left[\left(\frac{3R}{l_{D}}\right)+1\right]\right\}.\end{array}$$

$${\hat{Q}}_v^R=\frac{8N_A{e}_0{C}^0}{(R/l_D{)}^2}\left[-2\frac{e_0\zeta }{k_{B}T}-{\left(\frac{e_0\zeta }{3k_{B}T}\right)}^3\right].$$

$${\hat{Q}}_v^{REV}=\frac{1}{v_D\pi R_{REV}^2}{\int }_{R_{min}}^{R_{max}}{\hat{Q}}_v^R(R)\overline{v}(R)\pi R^{2}dN,$$

$${\hat{Q}}_v^{REV}=8N_A{e}_0{C}^0\left[-2\frac{e_0\zeta }{k_{B}T}-{\left(\frac{e_0\zeta }{3k_{B}T}\right)}^3\right]\frac{4-D}{2-D}{\left(\frac{l_D}{R_{max}}\right)}^2.$$

$${\hat{Q}}_v^{REV}=N_A{e}_0{C}^0{l}_D^2\left[-2\frac{e_0\zeta }{k_{B}T}-{\left(\frac{e_0\zeta }{3k_{B}T}\right)}^3\right]\frac{1}{{\tau }^2}\frac{\varphi }{k}.$$

$$\zeta (C^0)=a+b{\text{log}}_{10}(C^0),$$

$${\text{log}}_{10}({\hat{Q}}_v^{REV})=A_1+A_2{\text{log}}_{10}(k),$$

$$\begin{array}[]{ll}A_{1}=\text{log}_{10}\left\{\frac{N_{A}e_{0}C^{0}}{\gamma^{\frac{2-D}{4-D}}}\left[-2\frac{e_{0}\zeta}{k_{B}T}-\left(\frac{e_{0}\zeta}{3k_{B}T}\right)^{3}\right]\left(\frac{l_{D}}{\tau}\right)^{2}\right\},\\ A_{2}=-\frac{2}{4-D}.\end{array}$$

$${\hat{Q}}_v^{REV}=-\frac{C^{EK}\sigma \eta }{k},$$

$$\sigma =\frac{1}{F}\left({\sigma }_w+\frac{4mF{\Sigma }_S}{d}\right),$$

$$({\tau }_e{)}^2=F\varphi ,$$

Full text

(1) Consejo Nacional de Investigaciones Científicas y Técnicas, Facultad de Ciencias Astronómicas y Geofísicas, Universidad Nacional de La Plata, La Plata, Argentina

(2) Facultad de Ciencias Naturales y Museo, Universidad Nacional de La Plata, La Plata, Argentina

(3) Sorbonne Université, CNRS, EPHE, UMR 7619 Metis, F-75005, Paris, France

(+)Both authors contributed equally to this publication

(*)Corresponding author: luisg@fcaglp.unlp.edu.ar

This paper has been published in Journal of Geophysical Research – Solid Earth, please cite as:

L. Guarracino, D. Jougnot (2018) A physically-based analytical model to describe effective excess charge for streaming potential generation in water saturated porous media, Journal of Geophysical Research - Solid Earth, 123(1), 52-65, doi:10.1002/2017JB014873.

Keypoints

• •

Derivation of a physically-based analytical model to determine the effective excess charge

• •

Mechanistic explanation for the empirical dependence between effective excess charge and permeability

• •

The new model reproduces experimental data for different media and a broad range of ionic concentrations

1 Introduction

The self-potential (SP) method is a passive geophysical method based on the measurements of the electrical field which is naturally generated in the subsurface. The SP method was first proposed by Robert Fox in 1830 [Fox, 1830] and it is considered as one of the oldest geophysical methods. Although SP data are relatively easy to measure, the extraction of useful information is a non-trivial task since the recorded signals are a superposition of different SP components. In natural porous media, signals are mainly generated by electrokinetic (water flux) and electrochemical (ionic fluxes or redox reactions) phenomena. For more details of this method and for an overview of all possible SP sources we refer to Revil and Jardani [2013].

In the present study, we focus on the electrokinetic (EK) contribution to the SP, that is the part of the signal generated from the water flow in porous media (often referred to as streaming potential). The surface of the minerals that constitute porous media is generally electrically charged, creating an electrical double layer (EDL) containing an excess of charge that counterbalances the charge deficiency of the mineral surface [see Hunter, 1981, Leroy and Revil, 2004]. Figure 1a shows how the EDL is composed: a Stern layer that contains only counterions coating the mineral with a very limited thickness, and a diffuse layer that contains both counterions and coions but with a net excess charge. The shear plane, that can be approximated as the limit between the Stern layer and diffuse layer [e.g. Leroy and Revil, 2004], separates the mobile and immobile part of the water molecules when subjected to a pressure gradient. This plane is characterized by an electrical potential called $\zeta$-potential [see Hunter, 1981]. When the water flows through the pore it drags a fraction of the excess charge that give rise to a streaming current and an electrical potential field.

The EK phenomenon has been studied experimentally and theoretically for more than a century due its relevance in many practical applications in various fields (e.g., microfluidic, chemical engineering, fluid mechanics). In hydrogeophysics, its interest relies on the fact that SP data is sensitive to groundwater flow. Indeed, groundwater flow monitoring can otherwise be performed locally using rather intrusive methods (i.e. hydraulic head measurements between boreholes or in situ sensors, heat or chemical tracer evolution from a single well). However, given that SP only provides indirect measurements (electrical potential) of water flux, it is necessary to rely on a theoretical framework and petrophysical relationships that can predict the EK phenomenon at macroscopic scales.

Two main approaches to simulate streaming current generation in fully saturated porous media can be found in the literature. The classical approach relies on the use of the coupling coefficient, which is a rock-dependent property that relates the difference of hydraulic pressure to the difference in electrical potential. This property has first been described experimentally [e.g. Quincke, 1859, Dorn, 1880] and later quantified by Helmholtz [1879] and von Smoluchowski [1903]. The so-called Helmholtz-Smoluchowski coupling coefficient has been developed from a capillary-tube model and its final expression does not depend on the geometrical properties of the porous medium. Therefore, it has been used for any kind of medium under the assumption that the electrical conductivity of the mineral surface could be neglected. When it is not the case, alternative formula have been proposed by several researchers [e.g., Revil et al., 1999, Glover and Déry, 2010]. Various model using capillaries to predict the coupling coefficient under saturated and partially saturated conditions can be found in the literature [e.g. Rice and Whitehead, 1965, Ishido and Mizutani, 1981, Jackson, 2010, Jackson and Leinov, 2012, Thanh et al., 2017]. The second approach to simulate streaming current generation is more recent and focuses on the excess charge that is dragged by the water flow. The first reference to this approach in the English literature can be found in Kormiltsev et al. [1998] and more detailed theoretical frameworks can be found in Revil and Leroy [2004], Revil et al. [2007] and Revil [2016a, b]. In this approach, the coupling parameter is the excess charge that is effectively dragged by water in the porous media. Then, the streaming current and streaming potential distribution can be generated by multiplying the effective excess charge by the water velocity. Note that both approaches describe the same physics; the main difference between them relies on the parameter (coupling coefficient or excess charge density) used to describe the electrokinetic coupling between fluid flow and streaming potential generation.

In this work, we focus on the excess charge approach in the framework proposed by Sill [1983] and modified by Kormiltsev et al. [1998] and Revil et al. [2007], where SP signals can be directly related to the water velocity:

[TABLE]

where, $\sigma$ is the electrical conductivity (S m*-1*), $\varphi$ the electrical potential (V) and $\hat{Q}_{v}$ the excess charge density (C m*-3*) effectively dragged by the water flow at a given Darcy velocity $\textbf{v}_{D}$ (m s*-1*). From Eq. (1), it becomes clear that determining the effective excess charge is crucial to predict streaming potential signal in hydrosystems and to infer water flux from SP measurements.

Experimental evidences has shown that the effective excess charge depends on the porous medium permeability [Titov et al., 2002, Jardani et al., 2007, Bolève et al., 2012]. The first relationship to estimate the effective excess charge from permeability can be found in Titov et al. [2002] (Fig. 2 of their paper). Some years later, Jardani et al. [2007] proposed a more precise empirical relationship that has been proven to be very useful as it decreases the number of variables to be estimated and provides good estimates of efective excess charge for different types of porous media. This relationship has been used in many studies involving electrokinetic phenomena, such as dam leakage [e.g. Bolève et al., 2009, Ikard et al., 2014], surface-groundwater interaction [e.g. Linde et al., 2011], seismoelectric studies [e.g. Mahardika et al., 2012, Jougnot et al., 2013, Revil et al., 2015, Monachesi et al., 2015], hydraulically-active fracture identification [e.g. Roubinet et al., 2016] and permeability field characterization [e.g. Jardani and Revil, 2009, Soueid Ahmed et al., 2014]. However, Eq. (26) has been obtained regardless of pore water composition and other hydraulic parameters of porous media. Unfortunately, this relationship does not take into account the dependence of the EDL on the ionic concentration of the electrolyte [see the discussion in Jougnot et al., 2015].

The most rigorous approach to construct an accurate description of any phenomena at macroscopic scales is to start with a pore-scale description and then apply an upscaling method to the microscopic equations [Bernabé and Maineult, 2015]. In this study, the porous media is conceptualized as a bundle of capillary tubes with a fractal pore-size distribution. The porosity and permeability of the porous media are estimated using an equivalent medium theory. This geometrical description of the porous media and upscaling procedure have been successfully used to describe water flow in fractured media [Guarracino, 2006, Guarracino and Quintana, 2009], the evolution of multiphase flow properties during mineral dissolution [Guarracino et al., 2014] and relations between hydraulic parameters [Guarracino, 2007]. On the other hand, the effective excess charge in a single tube can be calculated from the radial distributions of excess charge and water velocity. Then the effective excess charge at the macroscopic scale is estimated using the flux-averaging technique proposed by Jougnot et al. [2012]. In that paper, they show that by combining these hydraulic and EK properties, a closed-form expression for the effective excess charge can be obtained in terms of permeability, porosity, tortuosity, ionic concentration, $\zeta$-potential, and Debye length. The dependence of the developed model on ionic concentration, permeability, porosity and also grain size is tested using different set of experimental data. It is also shown that the proposed model allows to derive from physical concepts the empirical relation between effective excess charge and permeability found by Jardani et al. [2007].

2 Theoretical development

The proposed model is based on the macroscopic description of effective excess charge in porous media from the upscaling of pore-size flow and electrokinetic phenomena. The porous medium is conceptualized as an equivalent bundle of water-saturated capillary tubes with a fractal law distribution of pore sizes. First, we derive expressions for the most important macroscopic hydraulic parameters (porosity and permeability). Then, we obtain an approximate expression for the excess charge in a single tube and, from this result, we estimate the effective excess in the porous media in terms of the hydraulic parameters.

To derive both hydraulic and electrokinetic properties we consider as representative elementary volume (REV) a cylinder of radius $R_{REV}$ (m) and lenght $L$ (m). The pores are assumed to be circular capillary tubes with radii varing from a minimum pore radius $R_{min}$ (m) to a maximum pore radius $R_{max}$ (m).

The cumulative size-distribution of pores whose radii are greater than or equal to $R$ (m) is assumed to obey the following fractal law [Tyler and Wheatcraft, 1990, Yu et al., 2003, Guarracino et al., 2014]:

[TABLE]

where $D$ is the fractal dimension of pore and $0<R_{min}\leq R\leq R_{max}<R_{REV}$. By using the Sierpinski carpet (a classical fractal object) Tyler and Wheatcraft [1990] show that the fractal dimension of (2) ranges from 1 to 2, with the highest values associated with the finest textured soils. Capillary tube models with this type of fractal pore-size distributions have been successfully used to describe macroscopic process and hydraulic properties in different porous media [e.g. Yu et al., 2003, Guarracino et al., 2014, Soldi et al., 2017].

Note that the cumulative number of pores given by (2) decreases with the increase of pore radius $R$, then differentiating (2) with respect to $-R$ we can obtain the number of pores whose radii are in the infinitesimal range $R$ to $R+dR$:

[TABLE]

2.1 Hydraulic properties

The porosity $\phi$ of the REV defined above can be straightforward computed from its definition:

[TABLE]

where $V_{p}(R)=\pi R^{2}l$ (m3) is the volume of a single tube of radius $R$ and length $l$ (m).

Substituting (3) into (4) and solving the definite integral we obtain

[TABLE]

where $\tau=l/L$ is the dimensionless hydraulic tortuosity of the capillary tubes [Scheidegger, 1958]. Note that the model assumes a single value of $\tau$ for all capillary tubes, this value must be considered as a mean tortuosity value of all tube sizes.

In order to obtain both the permeability and the excess charge of the REV we need to describe the water flow in a single tube. For a laminar flow rate, the velocity distribution inside the tube can be described by the Poiseuille model [Bear, 1988]:

[TABLE]

where $r$ (m) is the distance from the pore wall ($r=0$) to the center of the tube ($r=R$), $\rho_{w}$ the water density (kg/m3), $g$ the gravitational acceleration (m/s2), $\eta$ the dynamic viscosity (Pa s), and $\Delta h$ the pressure head drop across the REV (m). The average velocity $\overline{v}$ (m/s) in the capillary tube has the following expression:

[TABLE]

The total volumetric flow through the REV $V_{Q}$ (m3/s) is the sum of the volumetric flow rates of all individual tubes. According to (7) and (3) $V_{Q}$ can be computed as follows:

[TABLE]

On the basis of Darcy’s law (macroscopic scale), $V_{Q}$ can also be expressed as

[TABLE]

$k$ being the intrinsic permeability of the porous media (m2).

Then, combining (8) and (9) we obtain the following expression for permeability $k$ in terms of the geometrical parameters of the porous media:

[TABLE]

It is important to remark that a similar equations for $\phi$ and $k$ have been recently derived by Soldi et al. [2017] assuming constrictive capillary tubes. In the limit case of straight tubes both expressions for $\phi$ and $k$ are identical (see equations (11) and (15) of their paper).

For most porous media it can be assumed that $R_{min}<<R_{max}$ [Yu and Li, 2001]. Then, equations (5) and (10) can be reduced to

[TABLE]

[TABLE]

Finally, combining (11) and (12) we obtain a simple relationship to estimate permeabilty from porosity

[TABLE]

where $\gamma=\frac{R_{REV}^{2}}{8\tau(4-D)}\left(\frac{2-D}{\tau D}\right)^{\frac{4-D}{2-D}}$. Note that for $D=1$ the exponent of the porosity is 3 and (13) is equivalent to Kozeny’s equation [Kozeny, 1927].

2.2 Electrokinetic properties

Since electrokinetic phenomenon is caused by the coupling of fluid flow and charge distribution at pore scale, the magnitude of this phenomenom will be mainly determined by the macroscopic hydraulic and electrical properties of the porous medium. Based on the previous description of hydraulic properties we will compute the effective excess charge density $\hat{Q}_{v}^{REV}$ carried by the water flow in the REV (C/m3). The effective excess charge density, also called dynamic excess charge depending on the authors and symbolized by $\hat{Q}_{v}$ or $\bar{Q}_{v}^{eff}$, has to be distinguished from the other excess charge densities contained in the pore space: $Q_{v}$ the total excess charge density (C/m3), which includes all the charges from the Stern and diffuse layers of the EDL, and $\bar{Q}_{v}$ the excess charge located only in the diffuse layer (C/m3) (Fig. 1a) [see the discussion in Revil, 2017]. It is important to remark that in this study we do not consider the charges that are located in the Stern layer as they are fixed on the pore wall and do not contribute to the streaming current.

We start this derivation by defining the excess charge distribution $\bar{Q}_{v}$ in a capillary tube saturated by a binary symetric 1:1 electrolyte (e.g., NaCl) as follows

[TABLE]

where $N_{A}$ is the Avogadro’s number (mol*-1*), $e_{0}$ the elementary charge (C), $C^{0}$ the ionic concentration far from the mineral surface (mol/m3), $\psi$ the local electrical potential in the pore water (V), $k_{B}$ the Boltzman constant (J/K), and $T$ is the absolute temperature (K). For the thin double layer assumption (i.e., the thickness of the double layer is small compared to the pore size) the local electrical potential can be expressed [Hunter, 1981]:

[TABLE]

[TABLE]

where $\zeta$ (V) is the $\zeta$-potential on the shear plane (Fig. 1a), $l_{D}$ the Debye lenght (m), and $\epsilon$ the water dielectric permittivity (F/m).

As proposed by Jougnot et al. [2012], the effective excess charge density $\hat{Q}^{R}_{v}$ carried by the water flow in a single tube of radius $R$ is defined by

[TABLE]

being $A$ the cross sectional area of the tube (m2). Using a polar coordinate system with the pole located in the center of the tube, equation (17) becomes

[TABLE]

In order to obtain a closed-form analytical expression of $\hat{Q}^{R}_{v}$ we approximate the exponential terms of (14) by a four-term Taylor serie:

[TABLE]

Substituting (19) in (14) and solving (18) we obtain

[TABLE]

For the thin double layer assumption (i.e., the thickness of the double layer is small compared to the pore size) we consider $l_{D}<<R$ and (20) can be reduced to

[TABLE]

Figure 2 (a) shows the effective excess charge $\hat{Q}^{R}_{v}$ predicted by (20) and (21) for a ionic concentration $C^{0}=1$ mol/m3. Note that even though the number of terms of equation (21) is drastically reduced, both equations predict similar values of $\hat{Q}^{R}_{v}$. In order to test the general validity of (21)under the thin double layer assumption, we compare approximate values of $\hat{Q}^{R}_{v}$ with exact values obtained by the numerical solution of (18) assuming pore-sizes $R$ greater that 5 Debye lenghts. Figure 2 (b) presents the goodness of the fit for different values of ionic concentration $C^{0}$. From the analysis of Figure 2, we conclude that equation (21) predicts fairly well the effective excess charge in capillary tubes for a wide range of radius and ionic concentration values.

The effective excess charge $\hat{Q}^{REV}_{v}$ carried by the water flow in the REV is defined by

[TABLE]

where $v_{D}=\frac{\rho_{w}g}{\eta}k\frac{\Delta h}{L}$ is the Darcy’s velocity (m/s) (macroscopic scale). Substituting (21), (7) and (3) in (22) and assuming $R_{min}<<R_{max}$ yields

[TABLE]

Finally, combining (11), (12) and (23) we obtain the following expression for $\hat{Q}^{REV}_{v}$:

[TABLE]

The above equation constitutes the main result of this paper. Note that (24) predicts the effective excess charge density in terms of both macroscopic hydraulic parameters (porosity, tortuosity and permeability) and electrokinetic parameters (ionic concentration, $\zeta$-potential and Debye lenght). This equation gives insight into the role of macroscopic hydraulic parameters and it can be considered a starting point for designing non-invasive methods to monitoring groundwater flow using self-potential measurements.

In order to study the role of porosity $\phi$, permeability $k$, and ionic concentration $C^{0}$ on effective excess charge, we perform a parametric analysis of equation (24). The ionic concentration dependence of $\zeta$-potential is assumed to obey the relation proposed by Pride and Morgan [1991]:

[TABLE]

where $a$ and $b$ are fitting parameters. For this study we use the parameter values obtained by Jaafar et al. [2009] on silicate-based materials for NaCl brine $a$=-6.43 mV and $b$=20.85 mV.

Figure 3 summarizes the parametric analysis of $\log_{10}\left(\hat{Q}^{REV}_{v}\right)$ for the following ranges of variability of ionic concentration $C^{0}$, permeability $k$ and porosity $\phi$: $10^{-1}$ mol/m3 $\leq C^{0}\leq$ $10^{3}$ mol/m3, $10^{-20}$ m2 $\leq k\leq$ $10^{-8}$ m2 and $0.05$ $\leq\phi\leq$ $0.5$. Figure 3a shows the effect of porosity and permeability on $\hat{Q}^{REV}_{v}$ for a fixed value of ionic concentration ($C^{0}=100$ mol/m3). It can be observed that $\hat{Q}^{REV}_{v}$ is strongly determined by permeability while porosity only produces a slightly increase of $\hat{Q}^{REV}_{v}$ values. Figure 3b shows the effect of porosity and ionic concentration on $\hat{Q}^{REV}_{v}$. As shown in this panel, these parameters can change $\hat{Q}^{REV}_{v}$ values in two orders of magnitude for a fixed value of permeability ($k=10^{-14}$ m2). An increase of $\hat{Q}^{REV}_{v}$ can be observed when ionic concentration decreases and porosity increases. Finally, Figure 3c shows the role of permeability and ionic concentration on $\hat{Q}^{REV}_{v}$ for a fixed value of porosity ($\phi=0.2$). It can be observed that the effect of permeability on $\hat{Q}^{REV}_{v}$ is much more significant than ionic concentration.

From this parametric analysis we can conclude that effective excess charge $\hat{Q}^{REV}_{v}$ is highly sensitive to permeability values. However, porosity and ionic concentration can modify $\hat{Q}^{REV}_{v}$ values in two orders of magnitude for a given value of permeability [see discussions in Jougnot et al., 2012, 2015]. The increase of porosity or the decrease of ionic concentration produce an increase of $\hat{Q}^{REV}_{v}$ values.

3 Comparison with the empirical relationship proposed by Jardani et al. [2007]

The empirical relationship to estimate the effective excess charge $\hat{Q}^{REV}_{v}$ from permeability $k$ proposed by Jardani et al. [2007] reads as follows:

[TABLE]

where $A_{1}=-9.2349$ and $A_{2}=-0.8219$ are constant values obtained by fitting (26) to a large set of experimental data that includes various lithologies and ionic concentrations (Fig. 4). Note that this equation has been obtained regardless of pore water composition and other hydraulic parameters of porous media.

The empirical relationship (26) can be derived from the proposed equation (24). Using (13) we obtain the following expression for porosity in term of permeability: $\phi=(k/\gamma)^{\frac{2-D}{4-D}}$. Replacing this expression in (24) and taking the logarithm on both sides of the resulting equation, we obtain exactly (26) but with the following constants:

[TABLE]

According to our model, the $\text{log}_{10}(\hat{Q}^{REV}_{v})$-intercept $A_{1}$ mainly depends on chemical and interface parameters while the slope $A_{2}$ only depends on the fractal dimension of the pore size distribution ($1<D<2$). Note that the predicted range of the slope is $-1<A_{2}<-0.666$ and that the value obtained by Jardani et al. [2007] $A_{2}=-0.8219$ corresponds to the fractal dimension $D=1.567$.

Figure 4 shows the fit of (27) to an extensive set of $\hat{Q}^{REV}_{v}$ data determined by several authors [Ahmad, 1964, Cassagrande, 1983, Friborg, 1996, Pengra et al., 1999, Revil et al., 2005, 2007, 2012, Boleve et al., 2007, Jardani et al., 2007, Glover and Déry, 2010, Zhu and Toksöz, 2012, Jougnot and Linde, 2013]. The best fit of our model is obtained for fractal dimension $D=1.571$. Note that this single value of $D$ can fit a wide range of soil textures, then it can be considered a mean value of all pore size distributions.

It is worth mentioning that $\hat{Q}^{REV}_{v}$ data displayed on Fig. 4 are obtained from electrokinetic coupling coefficient $C^{EK}$ values (V/Pa) using the following equation (Revil and Leroy [2004], Revil et al. [2007]):

[TABLE]

where $\sigma$ (S/m) is the electrical conductivity, $\eta$ (Pa s) the dynamic viscosity, and $k$ (m2) the permeability. The values of $\sigma$ and $k$ of each sample have been measured independently.

If $D=1.527$ allows us to retrieve the exact same trend as Jardani et al. [2007], we are left with a cloud of values that spread around the proposed model. To test our model with more accuracy, we focus on a selection of experimental data where all the parameters have been measured in the following section.

4 Application to laboratory data

4.1 Effect of salinity

In order to analyze the effect of salinity on $\hat{Q}^{REV}_{v}$ we test the proposed model (24) with laboratory data obtained by Pengra et al. [1999]. These authors performed an exhaustive petrophysical characterization of a collection of rock and glass bead samples and they measured the electrokinetic coupling coefficient $C^{EK}$ for different NaCl brine concentrations. For this test $\hat{Q}^{REV}_{v}$ values are obtained from $C^{EK}$ using (28). The only parameter of(24) which is not measured by Pengra et al. [1999] is the hydraulic tortuosity $\tau$. Thus, we fit this parameter using a least square algorithm. The ionic concentration dependence of $\zeta$-potential is assumed to obey (25).

Figure 5 shows the fit of (24) to experimental values of effective excess charge $\hat{Q}^{REV}_{v}$ measured at different ionic concentrations $C^{0}$ for 3 sandstones samples of different permeabilities and 1 fused glass bead sample. The fitted values of tortuosity and measured values of porosity and permeability of each sample are listed in the figure caption. Experimental data show that $\hat{Q}^{REV}_{v}$ decreases with the increase of ionic concentration, and this behaviour can be adequately described by the proposed model. The decrease of $\hat{Q}^{REV}_{v}$ with the increase of $C^{0}$ was also predicted and discussed by Jougnot et al. [2015].

4.2 Effect of porosity

To test the effect of porosity on $\hat{Q}^{REV}_{v}$ estimates we fit the proposed model to experimental values of $\hat{Q}^{REV}_{v}$ in terms of both $k$ and $k/\phi$ obtained by Pengra et al. [1999] for $C^{0}=0.2$ mol/L in different porous media. Figure 6a shows the fit of (24) to experimental data in terms of $k$. The fitting parameter is the ratio $\phi/\tau^{2}$ and the RMS of the fit is 41.97 C/m3. Figure 6b shows the fit of the model to the same experimental values of $\hat{Q}^{REV}_{v}$ but in terms of $k/\phi$. In this case the fitting parameter is $\tau^{2}$ and the RMS is reduced to 19.64 C/m3. Even though more experimental evaluations of the model are needed, this test shows that the inclusion of porosity improves the estimate of $\hat{Q}^{REV}_{v}$.

4.3 Effect of grain size

Among the available data in the SP literature, Glover and Déry [2010] studied the effect of grain size on the electrokinetic coupling coefficient $C^{EK}$. These authors measured $C^{EK}$ on 12 packs of glass bead having different grain diameters $d$ ($d\in[1,990]\times 10^{-6}$ m) at two pore water salinities: $C^{0}=2\times 10^{-4}$ and $2\times 10^{-3}$ mol/L.

To test the performance of our model for different grain sizes, we used the measurements of Glover and Déry [2010] to predict the effective excess charge. For each grain size, we used the permeability and mercury porosity listed in their Table 1. For ionic concentrations $C^{0}=2\times 10^{-4}$ and $2\times 10^{-3}$ mol/L we use, respectively, the mean values $\zeta=$ -71.62 and -24.89 mV of $\zeta$-potentials calculated for all the samples (given in their Table 2). Note that the glass bead diameter should not influence the interface properties of the quartz mineral, this is an artifact due to surface conductivity effects [see also Leroy et al., 2012, Li et al., 2016].

In order to compare our model prediction with the measurements, we compute the coupling coefficient $C^{EK}$ using (28) and the electrical conductivity model proposed by Glover and Déry [2010]:

[TABLE]

where $\sigma_{w}$ is the electrical conductivity of pore water (calculated from $C^{0}$, see [Sen and Goode, 1992]); $\Sigma_{S}$, the glass beads surface conductance; $m$, the cementation index; and $F=\phi^{-m}$, the formation factor of the bead samples. The only parameter which is not provided in the study is the hydraulic tortuosity $\tau$.

Figure 7a shows the predicted effective excess charge for the two ionic concentrations as a function of the glass bead diameter and $\tau=1.2$ (this tortuosity has been chosen as it gives the best fit for all the data presented on Fig. 7b at once). We also plot in the figure the prediction of Jardani et al. [2007] empirical relationship which only depends on permeability $k$. Figure 7b illustrates the pretty good fit between the coupling coefficient measurements and predicted values of our model for both pore water salinities. Note that $\tau=1.20$ is the only fitting parameter of the model and this parameter allows to describe all measured data (different glass bead sizes and two ionic concentration values). It is worth mentioning that we could improve the fit between the proposed model and experimental data by adjusting $\tau$ for each glass bead diameter (especially for $d\ll 5\times 10^{-6}$ m), but we prefer keeping the number of fitting parameters as small as possible.

5 Discussion and conclusion

The present study is focused on the estimate of the effective excess charge $\hat{Q}^{REV}_{v}$ in fully saturated porous media, a key parameter to understand and model the streaming potential generation. Based on physical and geometrical concepts we derive a closed-form analytical expression to estimate $\hat{Q}^{REV}_{v}$ from electrokinetic and hydraulic parameters. The mathematical development involves up-scaling procedures at pore and REV scales, similar to the numerical up-scaling proposed by Jougnot et al. [2012] for partially saturated soils. The extension of the present model to unsaturated conditions is not straightforward as a relationship between permeability and saturation needs to be derived for the fractal pore size distribution.

The first step of this study consists in the up-scaling of EK properties for a binary symetric 1:1 electrolyte from the EDL scale to the capillary tube or pore scale (i.e., from few nm to possibly several mm). The approximate expression (21) relates effective excess charge $\hat{Q}^{R}_{v}$ to the radius of the tube $R$ and chemical parameters (ionic concentration $C^{0}$, $\zeta$-potential and Debye length $l_{D}$). The accuracy of this equation has been tested to be correct as long as the electrical double layer from the walls of the capillary do not overlap ($R>5l_{D}$), i.e. under the thin double layer assumption (Fig. 2). An extension of this model for EDL overlapping could be obtained by truncation of the diffuse layer [e.g. Gonçalvès et al., 2007].

The second step consists in the up-scaling of both hydraulic and EK properties from the pore scale (capillary tube) to the REV scale. The REV of porous media is conceptualized as a bundle of capillary tubes with a fractal pore-size distribution. Fractal models have been proven very useful to obtain different hydraulic properties of rocks and soils [Tyler and Wheatcraft, 1990, Yu et al., 2003, Guarracino, 2007, Guarracino et al., 2014, Soldi et al., 2017]. In this study, the fractal distribution approach allows us to obtain simple expressions for porosity $\phi$ (11) permeability $k$ (12) and effective excess charge $\hat{Q}_{v}^{REV}$ (23) in terms of the fractal dimension $D$, the maximum pore radius $R_{max}$ and the REV radius $R_{REV}$. Finally, combining these macroscopic properties we obtain the closed-form expression (24) to estimate $\hat{Q}_{v}^{REV}$ from hydraulic ($\phi$, $k$, $\tau$) and pore water chemistry ($C^{0}$, $l_{D}$, $\zeta$) parameters.

The proposed model is derived from physical and geometrical concepts and provides a mechanistic framework for understanding the role of hydraulic parameters on $\hat{Q}_{v}^{REV}$. In particular, the model corroborates the empirical relationship (26) of Jardani et al. [2007] which relies on an increasing number of data sets with different lithologies and it is used by many researchers. It is important to remark that the link between effective excess charge and permeability is the subject of debates in the scientific community [e.g. Jouniaux and Zyserman, 2016, Zyserman et al., 2017, Revil, 2017] and our model provides a theoretical justification to this link. However, in this study we show that it would be better to link $\hat{Q}_{v}^{REV}$ to the ratio of permeability and porosity ($k/\phi$) instead of $k$ (see Fig. 6).

Another limitation of the relationship proposed by Jardani et al. [2007] is that it does not explicitly depend on the pore-water chemistry, i.e. on the $\zeta$-potential, Debye length and ionic concentration. Indeed, the newly proposed model takes into account the ionic concentration and its $\zeta$-potential dependence and is therefore able to reproduce laboratory data from Pengra et al. [1999] where all the parameters but one (tortuosity) have been estimated independently.

The proposed model shows that $\hat{Q}^{REV}_{v}$ contains information on the lithology as predicted by Revil and Jardani [2013] (p. 64). Indeed, three petrophysical parameters are explicitly identified in the model: the permeability, $k$, the porosity, $\phi$, and the hydraulic tortuosity, $\tau$. It is worth to emphazise here that the relation between hydraulic and electrical tortuosity is not straightforward [e.g. Clennell, 1997]. In the present work, we test the simple model of Winsauer et al. [1952] to determine the hydraulic tortuosity used in Eq. (24) based on parameters that can be measured electrically [see Clennell, 1997, for a discussion]:

[TABLE]

where $\tau_{e}$ is the hydraulic tortuosity determined electrically from the formation factor and the porosity. Table LABEL:table1 show the comparison between the fitted tortuosity $\tau$ to the one predicted by (30). One can see that predicted and fitted tortuosities fall fairly close one to an other, which is a promising preliminary result. However, further work is needed to establish a better relation between these electrical and hydraulic properties, therefore to link the effective excess charge $\hat{Q}^{REV}_{v}$ to electrical conductivity.

The proposed model represents a major step forward in understanding the links between hydraulic and electrokinetic parameters in the framework of the excess charge approach. The model provides a theoretical basis for the empirical relationship with the medium permeability and describes the dependence of the excess charge on other petrophysical parameters. The simplicity of the model and its excellent fit to different set of experimental data opens-up new possibilities for a broader use of the SP method in hydrogeophysics studies based on the increasingly popular excess charge approach. The analytical development of a model for partially saturated porous media using the presented approach will be the natural next step to pursue this work.

Acknowledgments

The data used to test the proposed model is listed in the references. This research is partially supported by Universidad Nacional de La Plata and Consejo Nacional de Investigaciones Científicas y Técnicas (Argentina). The authors strongly thank the editor and the three reviewer for their nice and really constructive comments.