Predicting Electrokinetic Coupling and Electrical Conductivity in Fractured Media Using a Fractal Distribution of Tortuous Capillary Fractures

: Electrokinetics methods have attracted increasing interest 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 coefﬁcient 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 ﬂow 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 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). That parameter 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, properties of pore water 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 in case of negligible surface electrical conductivity, 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 theliterature 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. 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 given by [62] where r (no units) is the relative permittivity of the fluid, 0 (F/m) is the dielectric permittivity in vacuum, ζ is the zeta potential, η (Pa s) is the dynamic viscosity, and σ w (S/m) is the electrical conductivity of the fluid. Equation (2) is valid when the surface electrical conductivity of porous media is neglected, e.g., [38]. When the surface electrical conductivity is taken into account, one can apply the modified HS equation given by, e.g., [36,40,41] where Σ s (S) is the specific surface conductance and Λ (m) is a characteristic length scale [63]. Addition to the HS or modified HS equation, one can express the SPCC via the excess charge density Q v (C/m 3 ) that is dragged by pore water in porous media as, e.g., [42][43][44] where k (m 2 ) and σ (S/m) are permeability and electrical conductivity of porous materials at fully saturated conditions, respectively.

Geometry of Fractured Media
It has been shown that the fractal geometry can be applied in the analysis of flow and transport properties in porous media in general and specifically in fractured media, e.g., [64][65][66][67][68]. Fractured media are assumed to be composed of the fractures and the surrounding matrix, e.g., [69]. 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 Figure 1. The fracture is approximately considered as plane with rectangular cross-section, whose widths follow the fractal scaling law, e.g., [70][71][72]. The aperture and width of parallel fractures are 2a (m) and w (m), respectively. In order to derive electrokinetic properties at macroscale, a representative elementary volume (REV) as a cube with the length of L o and the cross-sectional area A REV that is perpendicular to the flow direction is considered as shown in Figure 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., [71][72][73][74]: 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 minus sign in Equation (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., [71,72] 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., [71,72,75,76]: where β is the fracture aspect ratio.

Hydraulic Properties
The porosity of the REV can be calculated 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., [77]. 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 [68,78,79], 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 Figure 2) inside the fracture are, respectively, given by, e.g., [57,80] 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 Figure 2. The flow rate through a fracture is therefore given by It is noted that Equation (11) is the well-known cubic law, e.g., [81,82]. Consequently, the total volumetric flow through the REV is the sum of volumetric flow rates through all fractures and is given by: Combining Equations (5), (7), (11) and (12) yields: According to Darcy s law (macroscopic scale) for Newtonian fluid flow in porous media, V REV is expressed as Combining Equations (8), (13) and (14), the following is obtained In case of w max >> w min (α → 0), Equation (15) reduces to

Electrical Conductivity
Adapting the reported approach for the electrical conductivity in porous media, e.g., [46,83], 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 Figure 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 Equations (5), (19) and (20), the following is obtained Substituting A REV from Equation (8) into Equation (21), the electrical conductivity of porous media σ under saturated conditions is obtained as Equation (22) indicates that σ of fractured media is related to porosity, electrical conductivity of pore water, specific surface conductance and microstructural parameters of porous media (D f , φ, α, β, w max ). It has a similar form to published models for porous media, e.g., [37,[84][85][86][87][88][89]. From Equation (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, Equation (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 . Equation (22) can be rewritten as:

Streaming Current in the REV
Adapting the strategy given by [90], 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 Figure 2, v(y) (m/s) is the velocity profile given by Equation (9) andQ v (y) (C/m 3 ) is the charge distribution in a fracture. Under Debye-Hückel approximation and the fluid of a binary symmetric 1:1 electrolyte, e.g., [51,91],Q v (y) is given by, e.g., [57,92] Q where λ is the Debye length depending on properties of the fluid and not on properties of the solid surfaces, e.g., [36,93]. Combining Equations (9), (26) and (27), one obtains 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 [94]. Reservoirs can be saturated with brine with much higher ionic concentrations. Therefore, the Debye length is typically less than 10 nm at 25 • C [93]. In addition, typical characteristic size of pore and fracture aperture in geological media is tens of micrometer, e.g., [95]. Therefore, the Debye length is normally much smaller than the pore sizes (thin EDL assumption). As seen in Figure 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 Figure 3). Under that condition, Equation (28) is simplified to Therefore, the total streaming current through the REV is given by, e.g., [53,56] Combining Equations (5), (29) and (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 Equations (21), (32) and (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 Equations (31), (34) and (35) yields Consequently, the SPCC is obtained as Equation (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, Equation (37) reduces to the HS equation as shown by Equation (2). It is indicated that for negligible surface conductivity, the SPCC is independent of fracture size distribution and geometrically shaped pore structures.  (3) and (37), the characteristic length scale Λ (m) is obtained for fractured media as

Comparing Equations
It should be noted that Thanh et al. (2020) [56] and Rembert et al. (2020) [96] 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 Equations (9), (10) and (27), respectively. Combining Equations (9), (10), (27) and (39), one obtains Under a thin EDL condition, Equation (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 Equations (42) and (43), the following is obtained: Substituting A REV from Equation (8) into Equation (44), Q v in the REV is written as Equation (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,94]. Note that Equation (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 Equations (15) and (45) to get the following: Comparing Equations (21), (37) and (44), one can find the SPCC from Q v as Equation (47) is the same as Equation (4) that has been developed for porous media using the volume averaging procedure, e.g., [43,97]. However, Equation (47) is here explicitly developed from the assumptions made for fractured media. This confirms the findings of [47,48,94], 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 . Figure 4 shows the variation of the permeability of fractured media k with the fracture aspect ratio β and porosity φ predicted from Equation (15) for three values of the ratio of minimum to maximum apertures α (0.0001, 0.001 and 0.01) with representative values of w max = 200 × 10 −6 m, τ = 1.2: (a) φ is fixed at 0.15 and (b) β is fixed at 0.001. Note that we determine D f from φ and α using Equation (6). It is seen that the permeability increases with an increase of φ as reported in literature for fractured rocks, e.g., [74,[98][99][100]. It can be explained by the fact that larger porosity means larger space occupied by fractures leading to larger permeability. Additionally, Figure 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 [74].

Results and Discussion
(a)     (6) with the knowledge of φ and α. It is seen that the SPCC in magnitude increases with an increase of β for a given value of α. Additionally, one can see that when α increases, the SPCC also increases. The reason is that σ decreases with an increase of β and with increase of α. As deduced from Equations (25) and (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 Equation (22).  Figure 7 shows the variation of the Q v with the fracture aspect ratio β predicted from Equation (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 Figure 4 the permeability k increases with increasing β or α. That is in good match with results reported in [19] or [101] in which the Q v decreases with an increase of permeability k for porous media. Due to lack of experimental data for fractured media, the model result is only compared with the experimental data published by [102] 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. Figure 8 shows the comparison between the model result from Equation (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 [76]. Besides the 1:1 line (the solid line), the 1:10 and the 10:1 lines are also shown in Figure 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. Figure 8. Comparison between measured permeability and calculated one for three fractured limestone samples obtained from [102] using Equation (16).

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.