Predicting electrokinetic coupling and electrical conductivity in fractured media using a bundle of tortuous capillary fractures

The electrokinetics methods have a great potential to characterize hydrogeological processes in geological media, especially in complex hydrosystems such as fractured formations. In this work, we conceptualize fractured media as a bunch of parallel capillary fractures following the fractal size distribution. This conceptualization permits to obtain analytical models for both the electrical conductivity and the electrokinetic coupling in water saturated fractured media. We explore two different approaches to express the electrokinetic coupling. First, we express the streaming potential coupling coefficient as a function of the zeta potential and then we obtain the effective charge density in terms of macroscopic hydraulic and electrokinetic parameters of porous media. We show that when the surface electrical conductivity is negligible, the proposed models reduces to the previously proposed one based on a bundle of cylindrical capillaries. This model opens up a wide range of applications to monitor the water flow in fractured media.

Even if fractured media constitute a large fraction of hydrosystems, it remains largely under studied. The interest of geo-electrical methods to study fractured media is clearly established in the literature using electrical conductivity [e.g., [21][22][23] or streaming potential [e.g., 5,24,25], but the quantitative use of these methods remains scarce in the literature.
The nature of fractured media makes it a challenge for integrative methods such as electrical conductivity imaging and streaming potential [26]. Indeed, fractures usually present a huge contrast in terms of properties with the surrounding matrix making it difficult to be clearly identified through diffusion-based methods such as electrical resistivity tomography [e.g., 27,28]. Different numerical methods have been proven useful for both electrical conductivity and self-potential simulations, among which the very fine discretization in finite-element modelling [e.g., 29,30] or discrete-dual-porosity [e.g., 23,31]. Nevertheless, analytical models and petrophysical relationships can be beneficial to the development and use of geo-electrical methods through for faster and simpler approaches.
Petrophysical relationships relating electrical conductivity and structural properties in fractured media are not abundant in the literature. Among other works, Bernabe [32] uses crack networks and percolation to study the transport properties, while Roubinet et al. [33] revisited the classic model of Archie [34] to relate topological and electrical properties. To the best of the authors knowledge, no analytical petrophysical relationships for fractured media exist in the literature.
Regarding the modeling of electrokinetic coupling, the most crucial parameter is the streaming potential coupling coefficient (SPCC). The SPCC controls the coupling between the fluid flow and the electrical field in porous media. It can be obtained through either the Helmholtz-Smoluchoski (HS) equation or the effective excess charge density approach. The HS equation is commonly used to determine the SPCC via the zeta potential at the solid liquid interface, fluid properties for a single cylindrical capillary [e.g., 35,36]. Applying a volume averaging procedure, Pride [37] showed that the HS equation is still valid for porous media as long as the surface conductivity can be neglected [see also , 38]. Modified HS equations have been proposed in the literature when the surface conductivity cannot be neglected [e.g., [39][40][41]. Then, an alternative approach allows the determination of the SPCC via the excess charge that is effectively dragged by the water flow in porous media [e.g., [42][43][44]. Capillary-based models have been successfully applied to describe SP in porous media, i.e., the porous medium is described as a bundle of parallel cylindrical tubes [e.g., [45][46][47][48][49][50][51][52][53][54][55][56]. However, microstructure of porous media is normally very complex with tortuous and non-circular pores. Recently, Shi et al. (2018) have theoretically studied the dependence of dynamic electrokinetic coupling coefficient on the electrical double layer thickness for a single cylindrical capillary and a single capillary fracture [57]. Then, they directly extended the obtained electrokinetic coupling coefficient from a single capillary to porous media without applying the upscaling technique. It should be noted that the electroosmosis that is a opposite effect to the SP in porous media has also been studied using capillary bundle models with different capillary geometry such as rectangular, cylindrical and annular geometries [e.g., [58][59][60][61]. Nevertheless, to the best of our knowledge, no model has been published to describe both the electrical conductivity an the electrokinetic coupling using a distribution of rectangular (i.e., fracture-like) capillaries.
Therefore, in this work, we conceptualize porous media as a bunch of tortuous parallel fractures or rectangular cross-sections following the fractal pore size distribution to study electrical conductivity and electrokinetic coupling in fractured media. Namely, we develop the model for the SPCC that is expressed in terms of both the zeta potential and the effective charge density. We also propose an analytical model for the effective charge density in terms of macroscopic hydraulic and electrokinetic parameters of porous media. Incidentally, we also obtain an expression for the electrical conductivity of saturated fractured media. All the proposed models are then compared with experimental data available in literature and published models.

Theory of the electrokinetic coupling
The streaming current is controlled by the relative movement between the charged solid surfaces and pore fluid and is directly related to an electric double layer (EDL) existing at the interface between the fluid and solid surfaces [e.g., 35]. This streaming current is then balanced out by a conduction current, leading to a so-called streaming potential. In other words, the water flowing through porous media creates a streaming potential [e.g., 62]. The important parameter for SP is the SPCC. At the steady state condition, the SPCC is defined as [e.g., 39]: where ∆V (V) and ∆P (Pa) are the generated streaming potential and applied fluid pressure difference, respectively. There normally exist two expressions to determine the SPCC at fully saturated conditions in literature. The first one is the HS equation that links the SPCC with the zeta potential ζ (V) at the solid-liquid interface and fluid properties and is given by [63] where Σ s (S) the specific surface conductance and Λ (m) is a characteristic length scale [64]. Addition to the HS equation or modified HS equation, the SPCC can be expressed via the effective excess charge density dragged by pore water in porous media as [e.g., [42][43][44] where Q v (C/m 3 ), k (m 2 ) and σ (S/m) are the effective excess charge density, permeability and electrical conductivity of porous media at fully saturated conditions, respectively. Note that the electrical conductivity plays a very important role in the electrokinetic coupling.

Geometry of fractured media
It has been shown that the fractal geometry can be applied in analysis of flow and transport properties in porous media in general and specifically in fractured media [e.g., [65][66][67][68][69]. Fractured media are assumed to be composed of the fractures and the surrounding matrix [e.g., 70]. The matrix permeability is normally much smaller than that of the fractures and thus the matrix can be considered as impermeable and no fluid exchange through the fracture walls. In this work, fractures are conceptualized as a bunch of parallel capillary slits as shown in Fig. 1. The fracture is approximately considered as plane with rectangular cross-section, whose widths follow the fractal scaling law [e.g., [71][72][73]. The aperture and width of parallel fractures are 2a (m) and w (m), respectively. In order to derive electrokinetic properties at macroscale, we consider a representative elementary volume (REV) as a cube with the length of L o and the cross-section area of the REV perpendicular to the flow direction of A REV (see Fig. 1). As mentioned, the REV is conceptualized as a bundle of fractures with width varying from w min to w max . It is shown that the number of fractures whose widths are in the range from w to w + dw in the REV is given by [e.g., [72][73][74][75]: where D f (no units) is the fractal dimension for pore space, 0 < D f < 2 in two-dimensional space and 0 < D f < 3 in three dimensional space. The negative sign in Eq. (5) implies that the number of fractures decreases with an increase of fracture width. The fractal dimension for fracture space (D f ) can be estimated from properties of porous media using the following relation [e.g., 72,73] where φ is the porosity of porous media and α is the ratio of the minimum width to the maximum width of the fractures (α = w min /w max ). It has been shown that the aperture a (m) is related to the width w (m) of the fracture by the following linear scaling law [e.g., 72,73,76,77]: where β is the fracture aspect ratio.

Hydraulic properties
The porosity of the REV defined above can be computed as the ratio of the total pore volume V p and the total volume V REV of the REV: where L o is the length of the REV, L τ is the real length of the fracture and τ=L τ /L o is the dimensionless hydraulic tortuosity of the microfracture [e.g., 78]. For the sake of simplicity, the length of the fractures is assumed to be independent of the fracture width, so an average tortuosity τ is considered in the model. It is trivial to take into account the correlation between the tortuosity and the capillary size as performed by [69,79,80], for example. However, this extra complexity does not really make the model any more representative of real porous media or affect the key results of the model [e.g., 46]. For a laminar flow, the velocity profile and average velocity (see Fig. 2) inside the fracture are, respectively, given by [e.g., 57,81] and v av = a 2 3ητ where ∆P is the fluid pressure difference, η is the dynamic viscosity of the fluid, a is the half aperture and y is the coordinate along the aperture as shown in Fig. 2.
The flow rate through a fracture is therefore given by It is noted that Eq. (11) is the well-known cubic law [e.g., 82,83] Consequently, the total volumetric flow through the REV is the sum of volumetric flow rates through all fractures and is given by: Combining Eq. (5), Eq. (7), Eq. (11) and Eq. (12) According to Darcy s law (macroscopic scale) for Newtonian fluid flow in porous media, V REV is expressed as Combining Eq. (8), Eq. (13) and Eq. (14), the following is obtained In case of w max >> w min (α → 0), Eq. (15) reduces to

Electrical conductivity
Adapting the reported approach for the electrical conductivity in porous media [e.g., 46,84], the bulk conductivity within a single fracture is given by The surface conductivity within a single fracture is given by recall that 2aw and 2w + 4a are the cross sectional area and the perimeter of the fracture (see Fig. 1). Consequently, the total electrical conductivity in a fracture is given by The electrical conductivity of the REV σ under saturated conditions is obtained by integrating over all saturated capillaries as Combining Eq. (5) Substituting A REV from Eq. (8) into Eq. (21), the electrical conductivity of porous media σ under saturated conditions is obtained as Eq. (22) indicates that the electrical conductivity of fractured media under saturated conditions is explicitly related to the porosity, electrical conductivity of the pore liquid, specific surface conductance and microstructural parameters of porous media (D f , φ, α, β, w max ). It has the similar form to published models for porous media [e.g., 37,[85][86][87][88][89][90]. From Eq. (22), the formation factor F is deduced as [34]: Roubinet et al. (2018) conducted a systematic numerical analysis for modeling electrical current flow in complex and heterogeneous fractured rocks [33]. For simple case of a simple fracture-matrix system with N f horizontal fractures embedded in a matrix having domain size L (τ = 1), fracture aperture b f and matrix porosity φ b , the formation factor of such a domain is given by [33] where m b is the cementation factor for the matrix. If a matrix is impervious to electrical current, that is, φ b = 0, for which the porosity only depends on the fractures as φ = N f b f /L, Eq. (24) now becomes for fractured media F = φ −1 . This is expected from simple geometry principles, as for parallel straight and horizontal capillaries with no surface conductivity: σ = φσ w . Eq. (22) can be rewritten as:

Streaming potential coupling coefficient 3.4.1. Streaming current in the REV
Adapting the strategy given by [91], the streaming current in a fracture due to transport of charge by the fluid under a fluid pressure difference is given by where a is the half aperture and y is the coordinate along the aperture as shown in Fig.  2, v(y) (m/s) is the velocity profile given by Eq.
Under most environmental conditions, ionic strengths (i.e., a proxy for ionic concentration) in potable water are normally between 10 −3 and 10 −2 mol/L [95]. Reservoirs can be saturated with brine with much higher ionic concentrations. Therefore, the Debye length is typically less than 10 nm at 25 • C [94]. In addition, typical characteristic size of pore and fracture aperture in geological media is tens of micrometer [e.g., 96]. Therefore, the Debye length is normally much smaller than the pore sizes (thin EDL assumption). As seen in Fig.  3, when a is 20 times larger than λ, 1 − λ a tanh( a λ ) can be approximated as 1 with around 5% difference (see the red dot in Fig. 3). Under that condition, Eq. (28) is simplified to Therefore, the total streaming current through the REV is given by [e.g., 53,56] Combining Eq. (5), Eq. (29) and Eq. (30) yields recall that α =w min /w max .

Conduction current in the REV
The streaming current is responsible for the streaming potential ∆V that is set up across porous media under a fluid flow. This streaming potential will generate an electric conduction current that is opposite in direction with the streaming current. According to Ohm s law, the electrical conduction current density is given by recall that σ is the electrical conductivity of the REV that has been previously presented. Therefore, the total electrical conduction current through the REV is given by Combining Eq. (21), Eq. (32) and Eq. (33), the following is obtained

Streaming potential coupling coefficient
At steady state, the total electrical current through the REV is zero and one obtains Combining Eq. (31), Eq. (34) and Eq. (35) yields Consequently, the SPCC is obtained as Eq. (37) indicates that the SPCC for fractured media under fully saturated conditions is related to the zeta potential, fluid properties, specific surface conductance at the solid water interface and microstructural properties of the media (D f , α, w max and β). When Σ s = 0, Eq. (37) reduces to the HS equation as shown by Eq. (2). It is indicated that for the negligible surface conductivity, the SPCC is independent of fracture size distribution and geometrically shaped pore structures. Comparing Eq. (3) and Eq. (37), the characteristic length scale Λ (m) is obtained for fractured media as It should be noted that Thanh et al. (2020) [56] and Rembert et al. (2020) [97] also obtained expressions for the characteristic length scale in porous media.

Effective excess charge density
Following the formalism presented by [52] or [54], we determine the effective excess charge density Q v (C/m 3 ) carried by the water flow in the REV. The effective excess charge density Q w v (C/m 3 ) carried by the water flow in a single fracture with the width of w and aperture of 2a is defined by [e.g., 50,52,54]: recall that v(y), v av andQ v (y) is the velocity profile, the average velocity and the charge distribution in the fracture, respectively. Their expressions are presented in Eq.
Under a thin EDL condition, Eq. (40) may be simplified as The effective excess charge Q v (C/m 3 ) carried by the water flow in the REV is defined by [e.g., 50,52]: where A REV is the cross-section area of the REV perpendicular to the flow direction and v D is the Darcy s velocity (m/s) that is given by recall that k is the permeability of porous media. Combining Eq. (42) and Eq. (43), the following is obtained: Substituting A REV from Eq. (8) into Eq. (44), Q v in the REV is written as Eq. (45) shows the relationship between Q v , macroscopic hydraulic parameters of porous media (the permeability, porosity, tortuosity) and electrokinetic parameters (the zeta potential). It is seen that the fracture size distribution and geometrically shaped fracture structures (w max , α and β) do not directly appear in the closed-form equation for Q v as reported in literature [e.g., 52,95]. Note that Eq. (45) shows a behavior that is completely consistent with the empirical relationship proposed by [19]. To obtain the dependence of Q v on the parameters of w max , α and β, we combine Eq. (15) and Eq. (45) to get the following: Comparing Eq. (21), Eq. (37) and Eq. (44), the relationship between the streaming potential coupling coefficient and effective charge density is obtained as Eq. (47) is the same as Eq. (4) that has been developed for porous media using the volume averaging procedure [e.g., 43,98]. However, Eq. (47) is here explicitly developed from the assumptions made for fractured media. This confirms the findings of [47], [48] and [95], the coupling coefficient dependence to the pore space geometry (pore size distribution and pore shapes) can be taken into account in the effective excess charge density Q v .   . It is seen that the permeability increases with an increase of φ as reported in literature for fractured rocks [e.g., 75,[99][100][101]. It can be explained by the fact that larger porosity means larger space occupied by fractures leading to larger permeability. Additionally, Fig. 4 also shows that the permeability increases with an increase of β and α. The reason for the increase of k with an increase of β may be that larger the fracture aspect ratio leads to higher flow rate and therefore higher permeability. The increase of k with an increase of α can be explained by the increase of the average fracture aperture with α as shown by [75]. Figure 5 shows the variation of the electrical conductivity of fractured media σ with the fracture aspect ratio β and the fractal dimension of the fracture sizes D f for three values of α (0.0001, 0.001 and 0.01) with representative values of σ w = 10 −2 S/m, Σ s = 1.5 × 10 −9 S, F = 20 and w max = 200.10 −6 m: (a) D f is fixed at 1.7 and (b) β is fixed at 0.001. It is seen that the σ is very sensitive with β, D f and α. For given values of D f and α, σ decreases with an increase of β. It can be explained by the decrease of the surface conductivity when β increases. Additionally, σ increases with an increase of D f and with a decrease of α. That can be explained by an increase of the total number of fractures in the REV and therefore, the increase of surface conductivity with the increase of D f and with the decrease of α (e.g., see [75] for more details). Figure 6 predicts the variation of the SPCC with the fracture aspect ratio β and the maximum width w max predicted from Eq.   an increase of β and with increase of α. As deduced from Eq. (25) and Eq. (37), when σ decreases then the SPCC increases. The prediction also shows that the SPCC increases with an increase of w max . This may be attributed to the fact that a larger fracture size leads to a lower surface electrical conductivity and therefore an increase of the SPCC. When w max exceeds a certain value, the SPCC becomes independent of fracture sizes. The reason is that the surface electrical conductivity can be negligible at large fracture sizes as indicated by Eq. (22). Figure 7 shows the variation of the Q v with the fracture aspect ratio β predicted from Eq. (46) for three values of α (0.0001, 0.001 and 0.01) with representative values of σ w = 2.10 −2 S/m, Σ s = 1.10 −9 S, ζ = -30 mV, φ = 0.15 and w max = 200.10 −6 m. It is seen that the Q v decreases with an increase of β and increase of α. As shown in Fig. 4 the permeability k increases with increasing β or α. That is in good match with results reported in [19] or [102] in which the Q v decreases with an increase of permeability k for porous media.

Results and discussion
Due to lack of experimental data for fractured media, the model result is only compared with the experimental data published by [103] for three fractured limestones from different depths (samples A, B and C). The permeability, porosity and maximum width of fractured samples A, B and C are reported to be 1×10 −17 m 2 , 0.007, 80 µm; 4.69×10 −16 m 2 , 0.0107, 200 µm and 4.8×10 −17 m 2 , 0.006, 150 µm, respectively. Fig. 8 shows the comparison between the model result from Eq. (16) with the measured data for samples A, B and C. For the comparison, the model parameters β, τ and D f are determined by seeking a minimum root-mean-square error (RMSE) through the "fminsearch" function in the MATLAB. We found β = 0.009, τ = 1.5 and D f = 1.8 to give the minimum RMSE for fractured samples. Note that Ghanbarian et al. (2019) found the values of β between 0.1 and 0.001 by fitting their model to experimental data of tensile fractures in the Krafla fissure swarm of northeast Iceland [77]. Besides the 1:1 line (the solid line), the 1:10 and the 10:1 lines are also shown in Fig. 8 (the dashed lines). It is seen that the model predicts well the experimental data within  less than 1 order of magnitude. Nevertheless, additional experimental laboratory studies or well controlled field measurements should be conducted to test our model further.

Conclusions
In this work, we conceptualize porous media as a bunch of parallel capillary fractures following a fractal pore size distribution to obtain the model for both electrical conductivity and SPCC. The later is expressed in terms of both the zeta potential and the effective charge density. We also propose an analytical model originating from the effective charge density approach that depends on macroscopic hydraulic and electrokinetic parameters of porous media. This work shows that when the surface electrical conductivity can be neglected, the proposed models reduces to ones obtained by a volume averaging procedure or a cylindrical capillary bundle model. It means that, in saturated conditions, the proposed models do not depend on the pore size distribution and geometrically shaped pore structures (w max , α and β). We believe that the models proposed in this work can open-up new possibilities to study and monitor water flow in fractured media.