2.1. Geometry of the fracture

Measurements of fracture surface topographies demonstrate that they exhibit a roughness possessing long-range spatial correlations and a scale invariance characteristic of self-affinity (Schmittbuhl, Schmitt & Scholz Reference Schmittbuhl, Schmitt and Scholz1995). This means that the height distribution $h'(x',y')$ of an isotropic fracture surface, such as the fracture walls shown in figure 1, has the following scale invariance property:

(2.1)\begin{equation} \lambda^{H}f(\lambda^{H}\Delta h',\lambda r') = f(\Delta h',r'). \end{equation}

Here $f(\Delta h',r')$ is the probability density function (p.d.f.) of having a height difference $\Delta h'$ between two points belonging to the topography and separated by a horizontal distance $r'=\sqrt {(\Delta x')^2 + (\Delta y')^2}$, and $H$ denotes the Hurst exponent (Bouchaud, Lapasset & Planès Reference Bouchaud, Lapasset and Planès1990). Except for materials such as sandstone for which fracturing occurs in-between mineral grains (Boffa, Allain & Hulin Reference Boffa, Allain and Hulin1998), the value of the Hurst exponent has been measured to $0.8$ for a large range of different materials and length scales (Bouchaud et al. Reference Bouchaud, Lapasset and Planès1990); this value shall be used here. In addition, the p.d.f. in (2.1) has been measured to be Gaussian (Brown Reference Brown1995). Self-affinity imparts long-range correlations to the topography, which means that if the topography is isotropic, its autocorrelation function

(2.2)\begin{equation} \mathcal{C}(\Delta x',\Delta y') = \frac{1}{L'_x L'_y}\int_{0}^{L'x} \,{{\rm d} x}' \int_{0}^{L'_y} \,{{\rm d} y}' h'\left(x',y'\right)h'\left(x'+\Delta x',z'+\Delta y'\right) \end{equation}

is a function of $r'$ alone, and $C(0)-C(r')$ scales as a power law $r'^{2H}$ (Schmittbuhl et al. Reference Schmittbuhl, Schmitt and Scholz1995). Consequently, the power spectral density (PSD) $G(k')$, which is defined for an isotropic fracture as (Wang, Narasimhan & Scholz Reference Wang, Narasimhan and Scholz1988; Brown Reference Brown1995; Schmittbuhl et al. Reference Schmittbuhl, Schmitt and Scholz1995)

(2.3)\begin{equation} G(k') = \frac{1}{2{\rm \pi}}\int_{-\infty}^{\infty} \mathcal{C}(r') {\rm e}^{-{\rm i}k'r'}\, {\rm d}r' , \end{equation}

scales as $k'^{-2(H+1)}$.

Figure 1.

A synthetic single fracture with self-affine walls of Hurst exponent $0.8$ connects two reservoirs along the $x'$-axis. The aperture field $a'(x',y')$ has a mean $a'_{m} = 100\,\mathrm {\mu }$m and standard deviation $\sigma ' = 0.9\,a'_{m}$. The correlation length is $L_{c}= L'_x/2 = 10^3 a'_{m} = 10$ cm. The fracture is closed along its lateral boundaries (defined by $y'=0$ and $y'=L'_y$). The walls carry a surface charge distribution, represented by an equivalent zeta potential $\zeta '(x',y')$. The outlet and the inlet are subjected to different levels of pressure, concentration and electrical potential.

The two walls constituting a fracture possess the same self-affinity property, with identical Hurst exponent and standard deviation of $h'(x',z')$. The aperture field is defined as

(2.4)\begin{align} \left. \begin{array}{ll@{}} a'(x',y')= \max\left (h'_{t}(x',y') - h'_{b}(x',y')+a'_{m},0 \right) & \text{if } h'_{t}(x',y') + a'_{m} > h'_{b}(x',y'), \\ a'(x',y') = 0 & \text{otherwise}, \end{array} \right\} \end{align}

where the topographies for the upper and the lower walls, $h'_{t}$ and $h'_{b}$, are defined with a zero arithmetic mean and $a'_{m}$ is the fracture's mechanical aperture, which is the distance between the mean planes of the walls and is also the arithmetic mean of the aperture field $a'$ when the two fracture walls are not touching each other. At $(x,y)$ positions where $h'_{t}(x',y') - h'_{b}(x',y')+a'_{m}$ is negative, meaning that the two original wall topographies intersect each other, the aperture is set to zero, which is equivalent to considering perfect plastic closure of the fracture. Fracture walls pertaining to a freshly made fracture are identical down to very small scales, yielding $h'_{t}=h'_{b}$ and a constant aperture $a'(x,y)=a'_{m}\,\forall (x,y)$. But the walls of a geological fracture are only identical at scales larger than a characteristic correlation scale $L_{c}$. Below that scale they are both self-affine but uncorrelated to each other, so the aperture field is also self-affine (since self-affinity is a linear property). The PSD of the aperture field is then of the form

(2.5)\begin{equation} \left. \begin{array}{ll@{}} G(k') \propto ({k'_x}^2+{k'_y}^2)^{-(H+1)} & \text{if } \sqrt{{k'_x}^2+{k'_y}^2}\leq k'_{c}=\dfrac{\rm \pi}{L_{c}}, \\ G(k') = G(k'_{c}) & \text{otherwise}. \end{array} \right\} \end{equation}

Consequently, the aperture field of a geological fracture is well modelled as a stationary 2-D random process (Brown Reference Brown1995) defined by (i) the functional form of the p.d.f. for $a'$ – this includes specification of the mean $a'_{m}$ and standard deviation $\sigma$; (ii) the Hurst exponent $H$; (iii) the correlation length $L_{c}$; and (iv) the horizontal dimensions $L'_x$ (length) and $L'_y$ (width) of the fracture (Méheust & Schmittbuhl Reference Méheust and Schmittbuhl2003). The synthetic fractures of § 5.2 have been generated using a spectral method previously employed by Méheust & Schmittbuhl (Reference Méheust and Schmittbuhl2001, Reference Méheust and Schmittbuhl2003). The perfect plastic closure in contact zones, expressed by the notation $\max (\cdot, 0)$ in (2.4), is the strongest approximation made in this model of geological fracture.

2.2. Boundary conditions and assumptions of the model

We consider steady coupled Stokes flow and ion/charge transport through a synthetic fracture as shown in figure 1. The $x'$ and $y'$ axes run along the mean fracture plane, while the $z'$-axis is normal to that plane, hence running along the fracture's aperture. The top (respectively, bottom) surface has a surface charge density $q'_{t}(x',y')$ (respectively, $q'_{b}(x',y')$), which for the moment is considered to be an arbitrary function of $x'$ and $y'$. We assume that the fracture connects two reservoirs along the $x$-axis which contain solutions of the same $1:1$ symmetrical electrolyte, of respective concentrations $c'_{in}$ and $c'_{ex}$. A solution of the same electrolyte also saturates the fracture. The pressures at the inlet and outlet are maintained at $p'_{in}$ and $p'_{ex}$, respectively, while the electrical potential difference imposed between the two reservoirs is $\Delta V$. We further assume that there is no externally imposed electrical potential difference between the two fracture walls, as such a potential difference would not impact the flow and ion transport in the fracture. We further assume that the faces of the fracture lying perpendicular to the $y'$-axis, at $y'=0$ and $y'=L'_y$, are closed. The presence of surface charge on the fracture walls and electrolyte solution inside the fracture leads to the formation of EDLs in the vicinity of the walls. We assume that the EDLs are non-overlapping, so that for any $(x_0,y_0)$ position along the mean fracture plane, the region away from either of the walls lying on the line segment $(x=x_0,y=y_0)$ remains electroneutral. The fluid viscosity is denoted $\eta$, its density $\rho _f$ and the liquid's electrical permittivity $\gamma$. All these properties are assumed to remain constant throughout.

We present in the following a formulation for coupled fluid and charge transport in a fracture geometry based on the lubrication approximation. This model makes the following several assumptions.

1. (i) The topography varies slowly as a function of the horizontal coordinates, i.e. $\forall (x,y), \| \boldsymbol {\nabla } a \| \ll 1$. Note that the aperture field's self-affinity implies that the aperture gradient goes to infinity when computed over a vanishingly small horizontal length scale (see Méheust & Schmittbuhl Reference Méheust and Schmittbuhl2001); therefore, depending on the size of the computing mesh, the aperture topography may have to be smoothed below a length scale $\ell _c$ such that $\varepsilon = \Delta a_0 / \ell _c \ll 1$, where $\Delta a_0$ is the typical vertical variation in the aperture field over $\ell _c$, so as to ensure that the lubrication approximation be valid. However, Brown (Brown Reference Brown1987; Brown et al. Reference Brown, Stockman and Reeves1995) has shown that small wavelength features of the surface roughness only play a secondary role in controlling the flow heterogeneity at the fracture scale, so, if $\ell _c \ll L_{c}$, transport processes at the fracture scale will not be significantly impacted by this small-scale smoothing of the aperture field. The validity of the lubrication approximation will be verified by comparing the results of the lubrication-based model to the numerical solutions of the first principle equations in identical geometrical and electrical configurations in § 5.1.

2. (ii) For the sake of generality, we assume the surface charge to be variable as well, but its gradient is also much smaller than 1, which means that the characteristic length scales of variation $l_0$ (say) of the surface charges are of same order as those of the topography ($l_0\sim L_c$).

3. (iii) The Reynolds number is assumed to be small ($Re\ll 1$), which is a reasonable assumption for most configurations of fracture flow in the subsurface.

4. (iv) We assume that the magnitude of the surface charge and potential gradient imposed across the fracture are not asymptotically large (see details in § 4.1).

5. (v) Finally, since the theory is derived from the Poisson–Nernst–Planck equations, the usual assumptions are made for these equations (for example, point charges, no non-Coulombic interactions between ions, etc.). Although several modifications to the Nernst–Planck equations have been proposed (Bazant et al. Reference Bazant, Kilic, Storey and Ajdari2009), the Nernst–Planck equations are capable of capturing the essential physics with a reasonable accuracy, especially when investigating effective transport characteristics in subsurface permeable media (see § 5.2).

Since we are applying the lubrication theory, conditions for primitive variables (e.g. velocity) are not required at the flow domain boundaries. However, boundary conditions for the externally imposed potential ($\phi '$), concentrations and pressure are required to close the system of governing equations, as we show later. These conditions are given by

(2.6a)\begin{align} & \text{At}\quad x' = 0, \quad p'=p'_{in};\quad c'_i=c_{in};\quad \phi'=0, \end{align}

(2.6b)\begin{align} & \text{At}\quad x' = L'_x, \quad p'=p'_{ex};\quad c'_i=c_{ex};\quad \phi'={-}\Delta V, \end{align}

(2.6c)\begin{align} & \text{At}\quad y' = 0 \ \text{and}\ L'_y \quad \partial p'/\partial y'= \partial c'_{{\pm}}/\partial y'= \partial \phi'/\partial y' = 0. \end{align}

The last condition expresses the constraint that the diffusive fluxes of pressure, solute concentration and electrical potential, have no component along the $y'$-axis at the lateral boundaries of the fracture.

Note that the generalized lubrication equations derived in the following assume no specific relation between $h'_b$ and $h'_t$; thus, they remain valid for symmetric fractures (i.e. with wall topographies that mirror each other) as well as non-symmetric fractures. As we shall show below, the only required geometrical input to the generalized lubrication equations is the aperture field.

4.1. The rescaled equations

As already discussed, the essence of the theory states that (Brown Reference Brown1987; Leal Reference Leal2007): the typical length scale of variation of the relevant quantities (such as pressure, velocity, concentration, etc.) along the fracture plane is much larger than that across the fracture thickness. This may be summarized here as $\varepsilon = h_0/l_0 \ll 1$ (Leal Reference Leal2007), where here $l_0$ can be the correlation length ($L_c$) and $h_0$ can be the mean aperture ($a_{{m}}$). Such a physical paradigm necessitates the following rescaling of the relevant variables and fluxes (Leal Reference Leal2007): $p_* = \varepsilon ^{-2}p$; $w_* = \varepsilon w$; $z_* = \varepsilon z$ and $\kappa _* = \varepsilon ^{-1}\kappa$, while the rest of the variables remain the same and their rescaled versions are also expressed without the subscript ‘$_*$.’ The purpose of the rescaling is essentially to make all the variables $\mathcal {O}(1)$ quantities.

Now, following the principles of the lubrication theory (Batchelor Reference Batchelor2000; Leal Reference Leal2007), we expand all the rescaled variables in a regular asymptotic series in $\varepsilon$ as

(4.1)\begin{equation} \varTheta = \varTheta_0 + \varepsilon\varTheta_1+\varepsilon^2\varTheta_2 + \cdots, \end{equation}

where the $\varTheta$ can represent quantities such as $u$, $v$, etc. We substitute the above asymptotic series into the rescaled governing equations and retain the leading-order terms in $\varepsilon$ to deduce the following:

(4.2a)$$\begin{gather} \frac{\partial^2 \varphi}{\partial z^2} + \frac{1}{2}\kappa^2\rho = 0 \quad \text{and} \quad \frac{\partial^2 \phi}{\partial z^2} = 0, \end{gather}$$

(4.2b)$$\begin{gather}\frac{\partial^2 c}{\partial z^2}+\frac{\partial}{\partial z}\left(\rho\frac{\partial \varphi}{\partial z}\right)=0\quad \text{and}\quad \frac{\partial^2 \rho}{\partial z^2}+\frac{\partial}{\partial z}\left(c\frac{\partial \varphi}{\partial z}\right)=0, \end{gather}$$

(4.2c)$$\begin{gather}-\frac{\partial p}{\partial z} +\alpha\frac{\partial^2 \varphi}{\partial z^2}\frac{\partial \varphi}{\partial z}=0\quad \text{and}\quad {-}\boldsymbol{\nabla}_{H} p + \frac{\partial^2 \boldsymbol{v}_{H}}{\partial z^2}+\alpha\frac{\partial^2 \varphi}{\partial z^2}\boldsymbol{\nabla}_{H}\left(\varphi+\phi\right)=0, \end{gather}$$

(4.2d)$$\begin{gather}\boldsymbol{\nabla}_{H}\boldsymbol{\cdot}\boldsymbol{v}_{H} + \frac{\partial w}{\partial z} = 0. \end{gather}$$

Here $\boldsymbol {\nabla }_{H} \equiv \boldsymbol {\hat {e}}_x({\partial }/{\partial x})+\boldsymbol {\hat {e}}_y ({\partial }/{\partial y})$ and $\boldsymbol {v}_{H} = u\boldsymbol {\hat {e}}_x + v\boldsymbol {\hat {e}}_z$ are respectively the projection of the gradient operator onto the fracture plane and the velocity along that plane. In (4.2) we have made use of the fact that to leading order, $\phi$ is only a function of $x$ and $y$ ($\phi = \phi (x,y)$), as we shortly demonstrate. The equations are subject to the following boundary conditions, after application of the lubrication approximation. At $z= z_{t}(x,y) = -1/2+h_{b}(x,y)$ and $z = z_{b}(x,y) = 1/2+h_{t}(x,y)$,

(4.3a)$$\begin{gather} \left[\frac{\partial c}{\partial z} + \rho\frac{\partial \varphi}{\partial z}\right] _{z=z_{{t}\vert{b}}} = \left[\frac{\partial \rho}{\partial z} + c\frac{\partial \varphi}{\partial z}\right] _{z=z_{{t}{b}}}= 0, \end{gather}$$

(4.3b)$$\begin{gather}\varphi =\zeta_{t}(x,y)\quad\text{or}\quad\zeta_{b}(x,y);\quad \left(\frac{\partial \phi}{\partial z}\right)_{z=z_{{t}{b}}} = 0;\quad \text{and}\quad\boldsymbol{v} _{z=z_{{t}\vert{b}}}= 0. \end{gather}$$

Note that in (4.3b) the specified surface charge density condition from (3.5a–c) has been replaced by a Diritchlet boundary condition with a known potential (the zeta potential (Hunter Reference Hunter2013)) that defines the potential drop across the diffuse part of the EDL. It can be easily verified that these two forms are largely equivalent, since at leading order the EDL can effectively be considered at equilibrium (Hunter Reference Hunter2013; Ghosh et al. Reference Ghosh, Chaudhury and Chakraborty2016). We assume that $\max \{\zeta \} \sim \mathcal {O}(1)$, i.e. the surface potential is not asymptotically large as compared with the thermal potential. The boundary conditions at $x=0$ and $x=L_x$ remain the same as given in (3.6).

The governing equations (4.2) can be solved semi-analytically. Based on (3.4b) and (4.3b), it is straightforward to establish that $\phi$ does not depend on $z$ ($\phi = \phi (x,y)$). However, the functional form of $\phi (x,y)$ is not yet known and we elaborate on how to find this functional dependence in the next subsection. Solving the leading-order Nernst–Planck equations (4.2b) yields ionic concentrations that satisfy the Boltzmann distribution, which leads to (Kilic et al. Reference Kilic, Bazant and Ajdari2007)

(4.4a,b)\begin{equation} c = 2\tilde{c}(x,y)\cosh (\varphi)\quad\text{and}\quad\rho ={-}2\tilde{c}(x,y)\sinh (\varphi) , \end{equation}

where $\tilde {c}(x,y)$ is the bulk concentration of salt, assumed to remain neutral owing to non-overlapping EDLs; in general it will be a function of $x$ and $y$. This remains true even if there is no externally imposed concentration gradient along the fracture, as we demonstrate later on. Note that, just as $\phi (x,y)$, $\tilde {c}(x,y)$ is also an unknown and will be determined in the next subsection. The EDL potential ($\varphi$) then satisfies the Poisson–Boltzmann equation, given by

(4.5)\begin{equation} \frac{\partial^2 \varphi}{\partial z^2} = \kappa^2\tilde{c}(x,y)\sinh \varphi . \end{equation}

Note that this equation can be computed using a local vertical coordinate $\tilde {z}$ obtained from $z$ through a translation such that $\tilde {z}=0$ at the bottom wall, and thus $\tilde {z} \in [0, a(x,y)]$. The equations expressing the associated boundary conditions, (4.3a) and (4.3b), are not impacted by this variable change, but the walls are located at $\tilde {z} = 0$ and $\tilde {z} =a(x,y)$ instead of $z_{{t}\vert {b}}$. Therefore, the distribution of $\varphi$ across the fracture only depends on the local aperture $a(x,y)$ and not on the particular pair of topographies $z_{{t}\vert {b}}$ that result in the aperture $a(x,y)$.

The $z$-momentum equation in (4.2c) can then be solved for pressure to deduce

(4.6)\begin{equation} p = \alpha\kappa^2 \tilde{c}(x,y)\cosh (\varphi) + p_0(x,y) , \end{equation}

in which the first term, called the osmotic pressure, comes from the salt concentration in the fracture (Brunet & Ajdari Reference Brunet and Ajdari2004), and the second term ($p_0(x,y)$) denotes the contributions from the induced pressure as well as the externally imposed pressure gradient. Similar to $\phi (x,y)$ and $\tilde {c}(x,y)$, $p_0(x,y)$ is yet to be determined and forms part of a closure problem. Inserting the pressure from (4.6) into the right-hand side equation in (4.2c), we obtain the following simplified form:

(4.7)\begin{equation} {-}\boldsymbol{\nabla}_{H} p_0-\alpha\kappa^2 \cosh (\varphi)\boldsymbol{\nabla}_{H} \tilde{c}+\alpha\frac{\partial^2 \varphi}{\partial z^2}\nabla_{H}\phi+\frac{\partial^2 \boldsymbol{v}_{H}}{\partial z^2}=0. \end{equation}

From the linearity of this equation, we deduce that the velocity components along the mean fracture plane, $u(x,y,z)$ and $v(x,y,z)$, can be written in the following forms:

(4.8a)$$\begin{gather} \boldsymbol{v}_{H} = \bar{u}_p\boldsymbol{\nabla}_{H} p_0+\bar{u}_{c} \boldsymbol{\nabla}_{H} \tilde{c}+\bar{u}_e\boldsymbol{\nabla}_{H} \phi,\quad \text{where}, \end{gather}$$

(4.8b)$$\begin{gather}\bar{u}_p = \tfrac{1}{2}\left\{z^2 - \left(h_{b}(x,y) + h_{t}(x,y)\right) z- \left(\tfrac{1}{2}+ h_{t}(x,y)\right)\left(\tfrac{1}{2} -h_{b}(x,y)\right)\right\} , \end{gather}$$

(4.8c)$$\begin{gather}\bar{u}_e = \alpha\left\{g(x,y)+ f(x,y)z-\varphi\right\}\quad \text{and} \end{gather}$$

(4.8d)$$\begin{gather}\bar{u}_{c} = \alpha\kappa^2\int_{{-}1/2+h_{b}}^{1/2+h_{t}} \,{\rm d}Z \left[\frac{(Z_{>}-\frac{1}{2}-h_{t})(Z_{<}+\frac{1}{2}-h_{b})}{a(x,y)}\right] \cosh\left[\varphi(Z)\right] . \end{gather}$$

In (4.8), $Z_{>} = \max (Z,z)$, $Z_{<} = \min (Z,z)$, and the functions $g(x,y)$ and $f(x,y)$ have the following expressions: $f(x,y) = \alpha (\zeta _{t}-\zeta _{b})/a(x,y)$ and $g(x,y) = \alpha \{\zeta _{t}(1/2-h_{b})+\zeta _{b} (1/2+h_{t})\}/a(x,y)$. We can use the velocity profiles in (4.8) to obtain the expressions for the three main fluxes in the fracture, namely, the cross-sectional volumetric flux ($Q$), current flux ($I$) and salt flux ($J$). The central idea of the lubrication theory is that these three quantities are conservative; we elaborate more on this aspect in § 4.2. The volume flux is given by

(4.9)\begin{equation} \boldsymbol{Q}_{H} = Q_x\boldsymbol{\hat{e}}_x + Q_y\boldsymbol{\hat{e}}_y = Q_p(x,y)\boldsymbol{\nabla}_{H} p_0+Q_{c}(x,y)\boldsymbol{\nabla}_{H} \tilde{c}+Q_e(x,y)\boldsymbol{\nabla}_{H} \phi, \end{equation}

where the functions $Q_p$, $Q_{c}$ and $Q_e$ are given by

(4.10)\begin{align} \left. \begin{gathered} Q_p(x,y) = \int_{{-}1/2+h_{b}(x,y)}^{1/2+h_{t}(x,y)} \bar{u}_p \,{\rm d}z ={-}\frac{1}{12} a^3(x,y) ; \quad Q_{c}(x,y) = \int_{{-}1/2+h_{b}(x,y)}^{1/2+h_{t}(x,y)} \bar{u}_{c} \,{\rm d}z,\\ Q_e(x,y) = \int_{{-}1/2+h_{b}(x,y)}^{1/2+h_{t}(x,y)} \bar{u}_e \,{\rm d}z . \end{gathered} \right\} \end{align}

The three components in the expression of the volume flux are contributions resulting respectively from the pressure, concentration and the potential gradients. Note that the volume flux $\boldsymbol {Q}_{H}$ is a vector with components along $x$ and $y$ only. Likewise the cross-sectional current and salt fluxes along the fracture plane can be expressed as

(4.11)\begin{gather} \boldsymbol{I}_{H} = I_x\boldsymbol{\hat{e}}_x + I_y\boldsymbol{\hat{e}}_y = I_p(x,y)\boldsymbol{\nabla}_{H} p_0+I_{c}(x,y)\boldsymbol{\nabla}_{H} \tilde{c}+I_e(x,y)\boldsymbol{\nabla}_{H} \phi, \end{gather}

(4.12)\begin{gather} \boldsymbol{J}_{H} = J_x\boldsymbol{\hat{e}}_x + J_y\boldsymbol{\hat{e}}_y = J_p(x,y)\boldsymbol{\nabla}_{H} p_0+J_{c}(x,y)\boldsymbol{\nabla}_{H} \tilde{c}+J_e(x,y)\boldsymbol{\nabla}_{H} \phi. \end{gather}

The different subscripts in the flux components bear the same meaning as those in the expression of volume flow rate components in (4.10). The prefactors in the components of (4.11) and (4.12) have the following expressions (taking hint from (3.3)):

(4.13a)$$\begin{gather} I_p = \int_{{-}1/2+h_{b}}^{1/2+h_{t}} \rho \bar{u}_p \,{\rm d}z; \quad I_{c} = \int_{{-}1/2+h_{b}}^{1/2+h_{t}} (\rho \bar{u}_{c}+2Pe^{{-}1}\sinh \varphi) \,{\rm d}z; \end{gather}$$

(4.13b)$$\begin{gather}I_e = \int_{{-}1/2+h_{b}}^{1/2+h_{t}} (\rho \bar{u}_e-2Pe^{{-}1}\tilde{c}\cosh \varphi) \,{\rm d}z; \quad J_p = \int_{{-}1/2+h_{b}}^{1/2+h_{t}} c \bar{u}_p \,{\rm d}z, \end{gather}$$

(4.13c)$$\begin{gather}J_{c} = \int_{{-}1/2+h_{b}}^{1/2+h_{t}} (c \bar{u}_{c}- 2Pe^{{-}1}\cosh \varphi) \,{\rm d}z;\quad J_e = \int_{{-}1/2+h_{b}}^{1/2+h_{t}} (c \bar{u}_e+ 2Pe^{{-}1}\tilde{c}\sinh \varphi) \,{\rm d}z. \end{gather}$$

Note that the concentration $c$ and density $\rho$ are obtained from $\tilde {c}$ and $\varphi$ through (4.4a,b). The final step is to write the equations enforcing that the fluxes $\boldsymbol {Q}_{H}$, $\boldsymbol {I}_{H}$ and $\boldsymbol {J}_{H}$ are conservative.

4.2. Closure: conservation of fluxes

In the classical lubrication theory of hydrodynamics (Batchelor Reference Batchelor2000; Leal Reference Leal2007) (i.e. without any electrical body forces), the volume flow rate conservation is a consequence of the continuity equation (Leal Reference Leal2007) for velocity, which implies that $\boldsymbol {Q}_{H}$ – in that case, defined by the sole first term on the right-hand side of (4.9) – is conservative ($\boldsymbol {\nabla }_{H}\boldsymbol {\cdot }\boldsymbol {Q}_{H} = 0$), which yields the Reynolds equation for pressure. In the present case the fluxes for fluid mass ($\boldsymbol {Q}_{H})$, electrical charge $(\boldsymbol {I}_{H})$ and solute mass $(\boldsymbol {J}_{H})$, defined respectively by (4.9), (4.11) and (4.12), are all conservative (i.e. divergence free), as we show in detail in Appendix A. The conservation of these three fluxes, combined with the expressions (4.9), (4.11) and (4.12), yields a set of three independent conservation equations for the fluid mass, electrical charge and salt mass, which constitute a generalized lubrication theory,

(4.14a)$$\begin{gather} \boldsymbol{\nabla}_{H}\boldsymbol{\cdot}\left(Q_p\boldsymbol{\nabla}_{H} p_0\right) + \boldsymbol{\nabla}_{H}\boldsymbol{\cdot}\left(Q_{c}\boldsymbol{\nabla}_{H} \tilde{c}\right) + \boldsymbol{\nabla}_{H}\boldsymbol{\cdot}\left(Q_e\boldsymbol{\nabla}_{H} \phi\right) = 0, \end{gather}$$

(4.14b)$$\begin{gather}\boldsymbol{\nabla}_{H}\boldsymbol{\cdot}\left(I_p\boldsymbol{\nabla}_{H} p_0\right) + \boldsymbol{\nabla}_{H}\boldsymbol{\cdot}\left(I_{c}\boldsymbol{\nabla}_{H} \tilde{c}\right) + \boldsymbol{\nabla}_{H}\boldsymbol{\cdot}\left(I_e\boldsymbol{\nabla}_{H} \phi\right) = 0, \end{gather}$$

(4.14c)$$\begin{gather}\boldsymbol{\nabla}_{H}\boldsymbol{\cdot}\left(J_p\boldsymbol{\nabla}_{H} p_0\right) + \boldsymbol{\nabla}_{H}\boldsymbol{\cdot}\left(J_{c}\boldsymbol{\nabla}_{H} \tilde{c}\right) + \boldsymbol{\nabla}_{H}\boldsymbol{\cdot}\left(J_e\boldsymbol{\nabla}_{H} \phi\right) = 0 . \end{gather}$$

The three unknowns to these three coupled equations are the functional dependences of $p_0$, $\tilde {c}$ and $\phi$ on $x$ and $y$. The expressions for the coefficients $Q_p$, $Q_{c}$, $I_{c}$, $J_e \ldots\, $, etc., which effectively act as diffusion coefficients in the equations, may be evaluated using (4.10) and (4.13), except for $\bar {u}_e$, whose definition includes the EDL potential $\varphi$. Note that the quantities $\varphi$, $c$ and $\rho$ explicitly depend on $\tilde {c}$, which indicates that the coefficients $Q_c,\;I_p,\;I_c,\ldots\, $, etc. all depend on $\tilde {c}$ as well. This makes the coupling between the (4.14) strongly nonlinear. Based on the conditions stated at the end of § 2.2, i.e. equations (2.6), we infer that the (4.14) are subject to the following boundary conditions for the said unknown functions:

(4.15a)$$\begin{gather} \text{At } x = 0, \forall y, \quad p_0 = p_{in}-\alpha\kappa^2;\quad \tilde{c}=1;\quad \phi=0, \end{gather}$$

(4.15b)$$\begin{gather}\text{At }x = L_x, \forall y, \quad p_0=p_{ex}-\alpha\kappa^2\tilde{c}_{ex};\quad \tilde{c}=c_{ex};\quad \phi={-}\beta L , \end{gather}$$

(4.15c)$$\begin{gather}\text{At } y = 0 \quad\text{and}\quad y=L_y, \forall x, \quad \frac{\partial p_0}{\partial y}= \frac{\partial \tilde{c}}{\partial y}= \frac{\partial \phi}{\partial y} = 0. \end{gather}$$

In the absence of any electrical charge and concentration gradients in the dissolved salt, (4.14) trivially simplify to the Reynolds equation (Brown Reference Brown1987)

(4.16)\begin{equation} \boldsymbol{\nabla}_{H}\boldsymbol{\cdot}\left(Q_p\boldsymbol{\nabla}_{H} p_{0}\right) = \boldsymbol{\nabla}_{H}\boldsymbol{\cdot}(a^3(x,y)\boldsymbol{\nabla}_{H} p_{0}) = 0, \end{equation}

where $a(x,z) = 1-h_{b}+h_{t}$ is the fracture's aperture field. It may be verified that the above flux coefficients (such as $Q_c,\;Q_e,\;I_p,\ldots\, $, etc.), when expressed in respective dimensional forms, satisfy the Onsager reciprocal relations (Brunet & Ajdari Reference Brunet and Ajdari2004) (see the seminal work of Onsager (Reference Onsager1931a,Reference Onsagerb) for the meaning of these relations in the general framework of non-equilibrium thermodynamic systems). Note that to satisfy the Onsager relation, one has to consider the chemical potential gradient, $\boldsymbol {\nabla }_H\tilde {\mu } = \tilde {c}^{-1}\boldsymbol {\nabla }_H \tilde {c}$, rather than the bulk concentration gradients $\boldsymbol {\nabla }_H \tilde {c}$. The validity of the Onsager relation herein stems from the fact that the EDL essentially remains in equilibrium, to leading order (Brunet & Ajdari Reference Brunet and Ajdari2004).

In summary, (4.14) are a generalization of the Reynolds equation to the conditions of electrohydrodynamic coupling. Instead of only solving one linear equation for pressure, as is the case for the Reynolds equation, we have to solve three coupled nonlinear equations for three unknowns: $p_0$, $\phi$ and $\tilde {c}$. This formalism, which extend the Reynolds equations to electrokinetic and diffuso-osmotic flows, has not been reported in the literature so far (see related discussion in paragraph 4.3 below). It can be applied to a wide range of flow pathways and surface potential distributions, provided that the assumptions of the lubrication theory (see § 2.2) are valid.

Note that detailed analytical solutions to the generalized lubrication equations in case of weakly charged 2-D geometries are included in § S2 of the supplementary material document.

4.3. The generalized lubrication equations as compared with the previous lubrication-based models of electrohydrodynamics

4.3.2. Analytic comparison to previous lubrication based models

We now proceed to show that the aforementioned approaches (Brown Reference Brown1989; Ghosal Reference Ghosal2002; Park et al. Reference Park, Russo, Branton and Stone2006) may be derived as special cases of the formulation carried out herein, when one ignores either one or both of the two additional conservation laws.

Brown (Reference Brown1989) considers electroneutral fractures. In § 5.2.4 we have established that our generalized lubrication equations boil down to those used by Brown in the special case of electroneutral fractures and in the absence of any imposed concentration gradients.

Park et al. (Reference Park, Russo, Branton and Stone2006) consider a very specific case of flow through a cylindrical bottleneck. Hence, their configuration is not directly applicable to a fracture geometry as ours is. In addition, Park et al. (Reference Park, Russo, Branton and Stone2006) disregard concentration variations, while their model is valid in the low surface charge ($\zeta \ll 1$) and weak field limits ($\beta \ll 1$). The low surface charge limit for 2-D conduits is also analysed in detail in the supplementary material (see § S2 therein). We shall therefore omit the details and simply note that, for $\zeta (x)\ll 1$, disregarding concentration variations ($\tilde {c}(x) = 1$), integrating (4.14) leads to the following asymptotic forms for the induced pressure ($p_0$) and the potential ($\phi$) in a 2-D confinement:

(4.17a)$$\begin{gather} \frac{{\rm d}p_0}{{\rm d} x} = \frac{\mathcal{K}_1}{a^3(x)} + \zeta_0\left[\frac{\mathcal{K}_2}{a(x)}+\frac{\mathcal{K}_3\bar{\zeta}(x)}{a^3(x)}\left(1-\frac{\tanh\omega}{\omega}\right)\right] + O(\zeta_0^2), \end{gather}$$

(4.17b)$$\begin{gather}\frac{{\rm d}\phi}{{\rm d} x} = \frac{\mathcal{K}_4}{a(x)} + O(\zeta_0). \end{gather}$$

Here $\zeta _{t} = \zeta _{b} = \zeta _0\bar {\zeta }(x)$ with $\zeta _0\ll 1$ and $\zeta \sim O(1)$. Additionally, $\omega = \kappa a(x)/2$ and $\mathcal {K}_1 - \mathcal {K}_4$ are constants. For further details, see equations (S10)–(S13) in the supplementary material. Although the above equations cannot be directly compared with those reported in Park et al.'s (Reference Park, Russo, Branton and Stone2006) work because of the critical difference in geometry (planar vs cylindrical), it is nevertheless apparent that the pressure and the potential gradients have similar forms in terms of their dependence on the geometry. However, the fact that the model by Park et al. (Reference Park, Russo, Branton and Stone2006) cannot consider a spatially varying concentration is a strong limitation in comparison to our model's capability.

Ghosal (Reference Ghosal2002) assumes the EDL to be very thin (Debye length $\kappa \gg 1$, see §§ 3 and 4.1 for its definition) and derives a modified Reynolds’ equation for the flow along a channel whose cross-section is of arbitrary shape and where flow occurs predominantly along the channel's axis. No variations in the solute (salt) concentration are considered by Ghosal. Introducing these additional limitations on the flow geometry, solute concentration field and EDL size, our generalized lubrication equations do reduce to those proposed by Ghosal (Reference Ghosal2002). Indeed, let us first recall that for a 2-D conduit (we disregard any flow along the $y$-direction, hence $\boldsymbol {\nabla }_{{H}}\equiv \hat {\boldsymbol {e}}_x\,\textrm {d}/{\textrm {d} x}$), the local cross-sectional area per unit width may be written as $\mathcal {A} = a(x)$. Then the coefficients appearing in (4.14) take the following form: $Q_p = -a^3(x)/12 = \int _{\mathcal {A}} \bar {u}_p \,\textrm {d}\mathcal {A}$. For thin EDLs (Ajdari Reference Ajdari1996), $\varphi = 0$ in the channel (except in a very thin region of size $O(\kappa ^{-1})$ close to the walls) and hence $\rho = 0$ and $c(x) = 2\tilde {c}(x)$. This yields $Q_e = \alpha \zeta (x) a(x) = \int _{\mathcal {A}} \bar {u}_e \,\textrm {d}\mathcal {A}$, wherein $u_e = \alpha \zeta (x)$ and $\zeta _b = \zeta _t = \zeta (x)$, whereas $Q_c = \alpha \kappa ^2Q_p$. Therefore, the electro-osmotic flow is essentially driven by the slip velocity at the edge of the EDL (Ajdari Reference Ajdari1996). The current and the salt fluxes (see (4.13)) take the following forms (Brunet & Ajdari Reference Brunet and Ajdari2004): $I_p,\; I_c \sim O(\kappa ^{-1}) \ll 1$, $I_e = -2Pe^{-1}\tilde {c}(x)a(x)+O(\kappa ^{-1})$, $J_p = 2\tilde {c}Q_p + O(\kappa ^{-1})$, $J_c = -2Pe^{-1}a(x)+\alpha \kappa ^2\tilde {c}(x)Q_p(x)+O(\kappa ^{-1})$ and $J_e = 2\alpha \tilde {c}(x)a(x)\zeta (x)+O(\kappa ^{-1})$. Substituting the above expressions for the various coefficients in (4.14), we deduce the following simplified version of the general lubrication equations:

(4.18a)$$\begin{gather} \frac{{\rm d}}{{\rm d} x}\left[Q_p(x)\frac{{\rm d}p}{{\rm d} x} + Q_e(x)\frac{{\rm d}\phi}{{\rm d} x}\right] = 0, \end{gather}$$

(4.18b)$$\begin{gather}\frac{{\rm d}}{{\rm d} x}\left[\tilde{c}(x)a(x)\frac{{\rm d}\phi}{{\rm d} x}\right] = 0 , \end{gather}$$

(4.18c)$$\begin{gather}\frac{{\rm d}}{{\rm d} x}\left[\tilde{c}(x)Q_p(x)\frac{{\rm d}p}{{\rm d} x} + \tilde{c}(x)Q_e(x)\frac{{\rm d}\phi}{{\rm d} x}+Pe^{{-}1}a(x)\frac{{\rm d}\tilde{c}}{{\rm d} x}\right] = 0 . \end{gather}$$

Note that in (4.18a) and (4.18c), we have used the total pressure $p = p_0+\alpha \kappa ^2\tilde {c}$. In the absence of any concentration gradient, $\tilde {c}(x) = 1$ and $\textrm {d}\tilde {c}/{\textrm {d} x} = 0$, so it follows that (4.18c) becomes identical to (4.18a) and, hence, the above equations further simplify to

(4.19a,b)\begin{equation} \frac{{\rm d}}{{\rm d} x}\left[Q_p(x)\frac{{\rm d}p}{{\rm d} x} + Q_e(x)\frac{{\rm d}\phi}{{\rm d} x}\right] = 0\quad \text{and}\quad \frac{{\rm d}}{{\rm d} x}\left[a(x)\frac{{\rm d}\phi}{{\rm d} x}\right] = 0 , \end{equation}

which are identical to (3.54) and (3.43) in Ghosal (Reference Ghosal2002). Hence, the model of Ghosal becomes a special case of the configurations that our model can address, namely those for which (i) there is a single dominant direction of flow, (ii) the solute concentration is homogeneous, and (iii) the EDL is very thin. In the configurations investigated in the results section below, assumptions (i) and (ii) are not valid while assumption (iii) is only valid to a limited extent.

However, (4.19a,b) can be extended to configurations for which the fracture geometry, concentrations, pressure and electrical potential depend both on the longitudinal coordinate $x$ and on the transverse in-plane coordinate $y$. It is thus possible to apply this model to one of the configurations otherwise addressed with our model in § 5 below, and compare their outputs. Such a comparison is presented in § S3 of the supplementary material, where we show that ignoring the equations for solute mass and charge conservations may lead to local relative errors larger than 90 % on the pressure field and larger than 70 % on the magnitude of the solute flux.

4.4. Numerical solutions to the generalized lubrication equations

We solve (4.14) numerically, using an iterative finite volume scheme (Patankar Reference Patankar1980). The algorithm is briefly discussed in what follows. We start with an initial guess for $p$, $\tilde {c}$ and $\phi$ (linear variation with $x$) and subsequently evaluate $\varphi$, $\bar {u}_p$, $\bar {u}_{c}$ and $\bar {u}_e$, from which $c$ and $\rho$ are evaluated at every vertical section in the fracture. Thereafter, we compute the axial flux coefficients, $Q_p$, $Q_e$, $I_p$, $I_c$, $J_e,\, \ldots\, $, etc. at every node of the horizontal mesh. These coefficients are first inserted in the (4.14a) to compute an updated $p_0$, while the older guesses for $\tilde {c}$ and $\phi$ are used to compute the last two terms therein. Subsequently, the fluxes and the updated $p_0$ are inserted in (4.14b) to update $\phi$ and finally, following the same procedure $\tilde {c}$ is updated from (4.14c). We then compute the relative errors in all the quantities as follows: for any variable $\xi$ (this may represent $p_0$, $\tilde {c}$, etc.), the relative error is defined as $\mathcal {E} = \max (|\xi -\xi _0|)/\max (|\xi _0|)$, where $\xi _0$ is the value of the relevant quantity in the previous iteration. If the maximum error falls below $10^{-5}$, the iterations are assumed to have converged. Otherwise, the above steps are repeated to compute the next corrections, until convergence occurs.

At every position $(x,y)$ in the mean fracture plane, the aperture has been discretized into 500 grid points to solve (4.5) for the EDL potential ($\varphi$), using a central finite difference scheme, and then compute $\bar {u}_{c}$. The cross-sectional fluxes are subsequently computed by numerically integrating the corresponding flux densities across these 500 grid points. The fracture plane was divided into a $100\times 100$ equally spaced grid, to solve the generalized equations; it has been verified that further refinement of the mesh does not alter the transport characteristics. For ease of computation, in places where the fracture is closed, i.e. $a = 1 + h_{t} - h_{b} = 0$, the non-dimensional aperture has been assigned a very small value (here $a = 10^{-6}$) to avoid having singular matrices after discretization. We have verified that lowering the limiting value below its current magnitude does not alter the final results. It is also important to note that below a certain value of the aperture, $a_{lim} = 10^{-2}$, the effect of surface charge has been neglected, for two main reasons: (i) below a certain value of $a$, the EDLs start to overlap and the assumption of zero centreline potential becomes invalid, which would call for special treatment of these parts in the fracture space; and (ii) since the aperture field in these areas have very small values (1 % of the average aperture), they would contribute very little to the overall transport anyway, so neglecting their contribution is a reasonably accurate approximation. We have also verified that lowering the limiting value ($a_{lim}$) does not change the transport characteristics in the fracture.

Note that $Q_{c}$, $Q_p$, $I_e,\ldots\, $, etc. can be negative (contrary to actual molecular diffusitivies) and, hence, in order to avoid a singular coefficient matrix, we evaluate these coefficients at all the faces of the control volumes, using a staggered grid (Patankar Reference Patankar1980). Once the solutions for $p_0$, $\tilde {c}$ and $\phi$ are known, the fluxes $\boldsymbol {Q}_{H}$, $\boldsymbol {J}_{H}$ and $\boldsymbol {I}_{H}$ can be easily deduced from (4.9), (4.11) and (4.12).