A new semi-analytic model for Stern-layer polarization in pore throats

SUMMARY

To explain induced polarization, membrane polarization is often referred to as a relevant process taking place in granular media – particularly, when narrow pore throats are present. This polarization effect is based on the membrane-like behaviour of pore throats caused by the presence of an usually negative charge on the pore surface, that influences charge transport in the pore fluid. Existing analytical, 1D models describe the pore system as a series of cylindrical pores with different radii and lengths. The polarization response is calculated by solving the Poisson–Nernst–Planck system for the current densities of one single anion and one single cation species representing the charge transport in the electrolyte and the diffuse layer at the pore surface. To include charge transport in the Stern layer, cations in the Stern layer have so far simply been considered by increasing the concentration of the diffuse layer cations. As we know from numerical modelling, this approach fails to predict the polarization response when the Stern layer is significantly charged. Here, we present a new semi-analytical model that treats the Stern-layer cations as a separate ion species and allows the Stern layer to polarize individually. To validate our new model, we compare it to the previously used analytical model and numerical simulations for different relative charges in Stern- and diffuse layer. We also use electrostatic surface-complexation models for two mineral surfaces (quartz and montmorillonite) to simulate the response of real geologic material under varying chemical conditions. This work is a step forward for considering realistic pore properties in induced-polarization modelling.

1 INTRODUCTION

Induced polarization (IP) is based on measuring the frequency-dependent behaviour of the complex conductivity. That frequency dependence is caused by polarization processes in geologic materials and is closely related to the structure of the pore space, particularly the pore surface (Börner & Schön 1991; Kruschwitz et al. 2010; Weller et al. 2010), and the electrochemical properties of the pore fluid (Revil 2012; Baierlein et al. 2016; Hördt et al. 2017). IP was first described by Schlumberger (1920) and is currently used in a broad field of applications, such as the prospection of minerals (e.g. Bleil 1953; Seigel et al. 2007; Revil et al. 2022), hydrogeology (e.g. Börner et al. 1996; Hördt et al. 2007; Slater 2007), permafrost research associated with climate warming (e.g. Maierhofer et al. 2022; Mudler et al. 2022) or for monitoring the transport and accumulation of contaminants in biogeophysics (e.g. Atekwana & Slater 2009; Wainwright et al. 2016).

Several empirical models have been proposed for the parametrization of measured complex conductivity spectra (e.g. Cole & Cole 1941; Dias 2000). However, only mechanistic microscale models can provide insight into the underlying physical processes.

In geological materials free of electron conductors, low-frequency polarization (typically between 1 mHz and 1 kHz) is commonly attributed to the electrical double layer (EDL), that forms due to electrochemical processes at the solid-electrolyte interface. These processes are usually described by surface-complexation models (SCM) that take the reactions at the mineral surface as well as the resulting electric potential in the electrolyte into account. Such models exist for a large number of mineral surfaces, including clays (Leroy & Revil 2004; Leroy et al. 2015) or amorphous silica and quartz (Leroy et al. 2008, 2013, 2022).

Early polarization models deal with the polarization of the Stern layer, the part of the EDL consisting of counter-ions (usually cations) directly adsorbed to the mineral surface, around a spherical particle embedded in an electrolyte solution (O’Konski 1960; Schwarz 1962; Schurr 1964). Lyklema et al. (1983) expand this particle-based model by including the coupling between charges of Stern layer and diffuse layer.

Polarization can also be understood in terms of pore constrictions acting as membranes. A membrane-like behaviour can be caused by the negative charge of the pore surface, which leads to an increased cation concentration and a decreased anion concentration and thus renders the transport of anions through the pore throat less effective than the transport of cations. This difference in transport numbers leads to a concentration polarization across the pore which is also known as membrane polarization.

Early membrane-polarization models consider the pore system as a sequence of “active” and “passive” zones (Marshall & Madden 1959), corresponding to a “wide” and “narrow” pore (compared to the thickness of the EDL), respectively. Membrane polarization is closely connected to the geometry of the pore space. This can be seen for example in a quadratic dependence of the relaxation time, the time the fully polarized system needs to relax, on the length of the wide pore (Blaschek & Hördt 2009; Volkmann & Klitzsch 2010) or the narrow pore, respectively (Titov et al. 2002). Bücker & Hördt (2013b) regard these two behaviours as limiting cases of the model of Marshall & Madden (1959).

Buchheim & Irmer (1979) discuss the influence of the EDL and of the pore’s cross-sectional area on membrane polarization. Gomaa & El-Diwany (2020) derive a simplified impedance equation for the case of thin diffuse layers (compared to pore size). Bücker & Hördt (2013a) expand the model of Marshall & Madden (1959) to include the influence of EDL and pore radius as they consider a radially varying ion concentration distribution. These authors also include the contribution of the Stern layer to the overall polarization in a largely simplified way by adding the positive charges in the Stern layer to the cations of the diffuse layer.

Hördt et al. (2017) study the influence of pore length and radius on the polarization model by Bücker & Hördt (2013a). On the basis of numerical simulations of the pore space, Bücker et al. (2019) investigate the coupled polarization of Stern layer and the diffuse layer in a cylindrical pore sequence. They show, that the model of Marshall & Madden (1959) together with the expansion by Bücker & Hördt (2013a) fails to reproduce the numerically calculated spectra in presence of a significant Stern-layer charge.

Here, we present a new semi-analytical model, that does not only include the ions in the diffuse layer, but also explicitly takes the cations in the Stern layer into account. The Stern-layer cations are treated as a separate ion species, that allow for an additional polarization of the Stern layer besides the classical membrane polarization of the diffuse layer. To include the Stern-layer cations, that are usually described by a surface charge density (e.g. in Schwarz 1962), into a 1D model, we define an effective Stern-layer ion concentration and specific boundary conditions to adapt for the strong binding of the cations to the surface. Additionally, we use surface-complexation models for the basal surface of montmorillonite and the (0001) surface of a quartz crystal to calculate the structure of the EDL for these two mineral surfaces under varying electrochemical conditions and use the results as an input for our new membrane polarization model.

2 STRUCTURE OF THE ELECTRICAL DOUBLE LAYER

The low-frequency polarization response strongly depends on the structure of the electrical double layer that in our model is described by the surface-charge density of the Stern-layer |$\Sigma _\mathrm{ S}$| and the potential |$\varphi _\mathrm{ d}$| at the d-plane as indicated in Fig. 1 (all symbols used in this manuscript are displayed in Tables 1 and 2). The values of these two parameters vary strongly depending on the mineral surface and the composition of the pore fluid, including pH and ionic strength, which equals the bulk ion concentration in the case of a monovalent, symmetric electrolyte. The respective structure of the EDL can be determined using an appropriate SCM.

In this study, we focus on the basal surface of montmorillonite and on a quartz (0001) crystal surface in contact with an NaCl solution representing the pore fluid. Fig. 1 shows sketches of the resulting EDL configuration for these two mineral surfaces. We assume both surfaces to exhibit a negative surface-charge density |$\Sigma _0 < 0$| at the 0-plane (⁠|$x_0$|⁠). At a certain distance (⁠|$x_\mathrm{S}-x_0$|⁠) that can be calculated from the ratio of the water permittivity to the capacitance |$C_1$| of region (1) (Leroy et al. 2013), cations are specifically adsorbed to the surface, forming the Stern layer with a surface-charge density |$\Sigma _\mathrm{S}$|⁠. At distances beyond the Stern plane (⁠|$x > x_\mathrm{S}$|⁠), within region (2), anion concentrations increase and cation concentrations decrease with the distance to the surface. This layer is called diffuse layer. We consider a basic Stern model (Westall & Hohl 1980), where the inner boundary of the diffuse layer (d-plane) coincides with the Stern plane at |$x = x_\mathrm{S}$|⁠.

The difference in anion and cation concentration in the diffuse layer leads to a spatial charge density, that can be described in terms of a surface-charge density |$\Sigma _\mathrm{d}$|⁠. Electro-neutrality implies that the surface-charge density of the Stern plane |$\Sigma _\mathrm{S}$| and of the diffuse layer |$\Sigma _\mathrm{d}$| cancel out the surface-charge density at the 0-plane |$\Sigma _0$|⁠, i.e.

The Debye length

is a measure for the extension of the diffuse layer, with |$\varepsilon _\mathrm{ r}$| being the relative permittivity of the electrolyte solution, |$\varepsilon _0 = 8.85\times 10^{-12}\, \mathrm{F\, m^{-1}}$| the permittivity of a vacuum, the bulk ion concentration |$c_0$| (in mol L|$^{-1}$|⁠), Faraday’s constant |$F = 96\,485.33\, \mathrm{A\, s\, mol^{-1}}$|⁠, the elementary charge |$e = 1.602\times 10^{-19}\, \mathrm{C}$|⁠, Boltzmann’s constant |$k_\mathrm{ B} = 1.380\times 10^{-23}\, \mathrm{J\, K^{-1}}$|⁠, and the absolute temperature T (in K). At sufficiently large distances of several Debye lengths, i.e. in region (3), both the anion and the cation concentration approach the equilibrium concentration |$c_0$|⁠.

2.1 Surface-complexation models

For both minerals, we implement a 1D basic Stern model (Leroy et al. 2013, 2015, 2022) using Comsol Multiphysics^®. The model consists of three regions, as displayed in Fig. 1: In region (1) between the mineral surface and the Stern-plane, we solve the Laplace equation

for the electric potential |$\varphi$| (in V). At |$x_0 = 0$| (the so called 0-plane), we apply the boundary condition

(Bourg et al. 2007). Note that eq. (4) implies that the potential is constant within the mineral. The position |$x_\mathrm{S}$| of the Stern plane, i.e. the boundary between regions (1) and (2), is calculated from the capacitance |$C_1$| (in F m|$^{-2}$|⁠) as

The capacitance |$C_1$| can be determined experimentally using geochemical and electrokinetic methods (Hiemstra & Van Riemsdijk 2006) and is treated as an input parameter, which we take from literature (see Table 3 for reference). At the Stern plane, the gradient of the electric potential is given by

In the diffuse layer region (2), we solve the Poisson–Boltzmann equation, which reads

for the studied case of a monovalent, symmtetric electrolyte. Region (3) represents the bulk aqueous electrolyte, where the electric potential equals zero and the ion concentrations are equal to the bulk concentration |$c_0$|⁠.

2.1.1 Montmorillonite

For montmorillonite (Mt), we use the SCM presented in Leroy et al. (2015), considering only the adsorption of |$\mathrm{Na^{+}}$| at clay surface sites |$\mathrm{\gt X^{-}}$|⁠:

where |$K_\mathrm{Na,Mt}$| is the corresponding equilibrium constant, which is an input parameter that can also be estimated by geochemical and electrokinetic methods. The negative surface sites on the montmorillonite basal surface arise from isomorphic substitution of Al³⁺ and Fe³⁺ ions by Mg²⁺ and Fe²⁺ ions in the octahedral sheet (Leroy et al. 2015). The surface-charge density |$\Sigma _0$| at the 0-plane is also an input parameter of our model. For the basal surface of montmorillonite, it is possible to relate |$\Sigma _0$| to the measured cation exchange capacity (CEC, Leroy et al. 2015). After some algebraic manipulations given in Appendix A1, we obtain the surface-charge density at the Stern plane from the mineral surface-charge density:

The activity |$a_\mathrm{Na}$| of sodium ions in the bulk electrolyte can be approximated by

with |$\gamma _\mathrm{Na}$| being the sodium activity coefficient in the bulk electrolyte that can be calculated by using Davies’ equation (Davies 1962):

2.1.2 Quartz

The EDL of a quartz (0001) crystal surface can be modelled with the SCM presented in Leroy et al. (2022). In this model, a protonation reaction takes place at the 0-plane:

with |$K_\mathrm{H}$| being the respective equilibrium constant. Here, we neglect further protonation reactions due to their comparably low-equilibrium constants (Leroy et al. 2013). At the Stern plane, |$\mathrm{Na}$|-ions are adsorbed:

After some algebraic transformations given in Appendix A2, we obtain the surface-charge densities at the 0-plane

and at the Stern plane

with

|$\Gamma _\mathrm{S}$| is the total surface-site density of |$\mathrm{\gt SiO^-}$|- and |$\mathrm{\gt SiOH}$|-sites and is treated as an input parameter. While the activity of |$\mathrm{Na^+}$| in the bulk electrolyte can be calculated using Davies’ equation (eqs 10 and 11), we calculate the proton activity in the bulk electrolyte from pH as

2.2 Modelled double-layer structure

From our surface-complexation model, we extract the potential |$\varphi _\mathrm{d}$| at the d-plane and the surface-charge density |$\Sigma _\mathrm{S}$| of the Stern layer as well as the partition coefficient (⁠|$0 \le f_Q \le 1$|⁠), which is the fraction of the counter-charge compensated by the Stern layer

To validate the numerical implementation of our model in Comsol Multiphysics^®, we compare these results to the predictions of a corresponding PhreeqC simulation for both mineral surfaces. PhreeqC is a widely used free geochemical modelling software from USGS (Parkhurst & Apelo 2013). The parameters used for the calculations are given in Table 3. Fig. 2 shows the calculated EDL parameters of the montmorillonite basal surface for different surface-charge densities with varying |$\mathrm{NaCl}$| concentration |$c_0$| in the electrolyte. The results for the quartz (0001) crystal surface at different |$\mathrm{NaCl}$| concentrations with varying pH are displayed in Fig. 3.

Both numerical model predictions are in good agreement with the corresponding PhreeqC results. We therefore conclude that our new implementation of the model yields the correct structure of the EDL. In the range of parameter values investigated, we found that the EDL of the basal surface of montmorillonite does not behave like the EDL of the quartz (0001) surface. For both montmorillonite and quartz, the modelled d-plane potential magnitude decreases when salinity increases. However, the surface-charge density of the Stern layer and partition coefficient are high and relatively constant for montmorillonite whereas they vary over a broad range of values for quartz.

Therefore, in the case of quartz, a wide range of possible EDL parameters has to be taken into account when studying a modelled IP response. The partition coefficient of the montmorillonite basal surface is high (above 50 per cent) whatever the investigated salinity is. Therefore, the relative contribution of the Stern layer to the modelled IP response of the basal surface of montmorillonite is expected to be high. For quartz, varying the |$\mathrm{NaCl}$| concentration and/or the pH leads to a strong change in the partition coefficient (from about 0 per cent up to 80 per cent), indicating that the relative importance of Stern layer and diffuse layer for the modelled IP response may change drastically depending on the chemical conditions.

3 POLARIZATION MODEL

3.1 Basic equations

To model the polarization of the pore system, we use a Poisson–Nernst–Planck equation system for a symmetric and monovalent electrolyte consisting of one cation (subindex p) and one anion species (subindex n), e.g. NaCl. In the following equations, we will not explicitly consider the ion valence |$z_{p,n} = 1$|⁠. Neglecting the effect of advective flow on the electric current density of the ions, the electric current densities |$\mathbf {j}_{p,n}$| of cations and anions are given by

In this equation, |$c_{p,n}$| (in mol m|$^{-3}$| = 10|$^3$| mol L|$^{-1}$|⁠) is the ion concentration and U the electric potential (in V), |$D_{p,n}$| (in m|$^2$| s|$^{-1}$|⁠) and |$\mu _{p,n}$| (in m|$^2$| V|$^{-1}$| s|$^{-1}$|⁠) are the ion diffusion coefficient and the ion mobility, respectively. These two parameters are linked to each other by the Nernst–Einstein equation |$D_{p,n} = \mu _{p,n} k_\mathrm{B} T / (e)$|⁠. The vector |$\mathbf {r}$| represents the spatial coordinate and t is the time.

The divergence of the current density and the temporal variation of the ion concentration are connected by the continuity equation

We now consider a harmonic perturbation |$\propto \exp (i \omega t)$| of the electric potential and both ion concentrations, where |$\omega$| is the angular frequency (in rad s|$^{-1}$|⁠):

with the equilibrium electric potential |$U^{(0)}$| and the equilibrium ion concentrations |$c_{p,n}^{(0)}$|⁠. The perturbation terms |$\delta U$| for the electric potential and |$\delta c_{p,n}$| for the ion concentrations are assumed to be small compared to the respective equilibrium quantities.

The equilibrium ion concentrations can be calculated from

In this equation, |$c_{\infty }$| is the equilibrium ion concentration in the bulk electrolyte. For the case of an external electric field, combining eqs (21) and (22) with eqs (19) and (20) and neglecting products of two perturbation terms yields the linearized equation

The effect of non-vanishing charge densities resulting from unequal cation and anion concentrations on the electric potential can be expressed using Poisson’s equation:

For the equilibrium case, we can use eq. (23), which results in

Using the assumption of a harmonic perturbation for the case with an external electric field, we obtain

Since Stern-layer cations arrange within a thin layer at the mineral surface, we regard them as a surface-charge density |$\Sigma _\mathrm{S}$|⁠. Analogously to the ion concentrations of cations and anions, a harmonic perturbation is expected for the Stern-layer surface-charge density:

where |$\Sigma _{S}^{(0)}$| is the equilibrium value and |$\delta \Sigma _\mathrm{S}$| is the perturbation of the Stern-layer surface-charge density. Analogously to eq. (24), one can obtain a Nernst–Planck equation for the Stern-layer cations as well:

The operator |$\nabla _\mathrm{S}$| means that gradient and divergence are only formed along the pore surface.

3.2 New semi-analytical model

We calculate the low-frequency polarization response of a system consisting of two cylindrical pores with radii |$r_1$| and |$r_2$| and lengths |$2 L_1$| and |$2 L_2$| representing a pore constriction (Figs 4a and b). The wider pore is referred to as pore 1 and the narrow pore is referred to as pore 2. We apply periodic boundary conditions at the ends of the pore system, describing an infinite sequence of connected pores of the same geometry. The pore space is assumed to be filled with a symmetric, monovalent electrolyte with one cation and one anion species (e.g. NaCl).

The Poisson–Nernst–Planck equation system is now transformed into a 1D system of partial differential equations. We define a separate coordinate system for each of the two pores, ranging from |$[-L_1, L_1]$| and from |$[-L_2, L_2]$|⁠, respectively, and solve the equation systems individually. The individual solutions are then matched by suitable boundary conditions as described below.

To transform eq. (24) to the 1D, the gradient operators are exchanged by partial derivatives |$\partial /\partial x$|⁠. We assume the variation of the equilibrium electric potential |$U^{(0)}$| along the pores to be zero, because variations in the equilibrium electric potential and ion concentrations only occur in radial direction. The radial variation can be described by eqs (23) and (26) and suitable boundary conditions at the pore surface. We use the approach of Bücker & Hördt (2013a), in which the ion concentrations of anions and cations are averaged over the pore’s cross-section. The resulting mean concentration is then divided by the bulk ion concentration, which yield the dimensionless factor |$b_{p,n}$|⁠. Using this factor, an effective ion mobility of |$\mu _{p,n} = b_{p,n} \mu _{p,n}^{(0)}$| can be considered to account for the radial variation of concentrations in the equilibrium case. For a more detailed introduction of this approach, please refer to Bücker & Hördt (2013a).

In the non-equilibrium case, we assume an external potential |$\pm V$| to be applied at the two ends of the pore system.

For simplicity, in the following, we use the symbols for the different quantities themselves (⁠|$c_p$|⁠, |$c_n$|⁠, |$\Sigma _\mathrm{S}$| and U) rather than the symbols for the respective perturbation quantities (⁠|$\delta c_p$|⁠, |$\delta c_n$|⁠, |$\delta \Sigma _\mathrm{S}$| and |$\delta U$|⁠). However, all equations in this section refer to the respective perturbation of the respective quantities. Although the system of equations is solved for both pore types individually, we here only present the equations for the wider pore for brevity.

The resulting partial differential equations for the cation and anion concentration in the electrolyte have the following form (Marshall & Madden 1959):

Here, we use the assumption that in the equilibrium state, the monovalent, symmetric electrolyte is electroneutral, i.e. |$c_{p}^{(0)} = c_{n}^{(0)} = c_\infty$| and the diffusion coefficients |$D_{p,n}$| as well as the ion mobilities |$\mu _{p,n}$| refer to the effective parameters computed following the approach by Bücker & Hördt (2013a) as described above.

To include Stern-layer polarization into the model, we use an equivalent ion concentration. To this end, we multiply the surface-charge density by the pore’s circumference and divide by the pore’s cross-sectional area. Dividing the result by Faraday’s constant yields an equivalent ion concentration (in mol m|$^{-3}$|⁠):

Note that this approach is an assumption rather than a rigorous derivation. Since for the pore constriction |$r_2 < r_1$|⁠, this leads to a higher equivalent Stern-layer ion density in the narrow pore. We model the Stern-layer cations as an independent third ion species. A Nernst–Planck equation for the Stern-layer cations replaces equation (29):

where |$c_{S,0}$| is the equivalent equilibrium ion density of the Stern-layer ions, calculated from |$\Sigma _{S}^{(0)}$| via eq. (32). The factor M was introduced by Lyklema et al. (1983) as

with |$\kappa$| being the inverse Debye length given by |$\kappa = \lambda _\mathrm{d}^{-1}$|⁠. Lyklema et al. (1983) originally derived the coefficient M for a spherical particle embedded in an electrolyte solution. It accounts for the interaction between the Stern-layer cations and the ions in the diffuse layer and, for a positively charged Stern layer (⁠|$M \ge 1$|⁠), leads to a smaller relaxation time of the cations in the Stern layer. Although M was originally derived for a spherical geometry, we will see further below that the effect is the same in our cylindrical model geometry. To take the influence of Stern-layer charges on the variation of the electric potential along the pore into account, the equivalent ion concentration of the Stern-layer ions has to be included into Poisson’s equation as

Eqs (30), (31), (33) and (35) apply to the narrow pore as well. We use overlined quantities (⁠|$\overline{c}_p$|⁠, |$\overline{c}_n$|⁠, |$\overline{c}_\mathrm{S}$| and |$\overline{U}$|⁠) to distinguish the narrow-pore quantities from the ones in the wider pore (Fig. 4c). Note that for the Stern-layer cations, the equilibrium concentration |$\overline{c}_{S}^{(0)}$| is different in both pores due to the radius dependency expressed in equation (32).

Since all three ion flux densities have to be continuous between the two cylindrical pores, we have to account for the smaller cross-section of the narrow pore when defining the boundary conditions. Using the cross-section ratio |$a = r_2^2 / r_1^2$| these boundary conditions write:

The ion concentration of the ions in the electrolyte is assumed to be continuous at the boundary, as well:

As the Stern layer is a thin layer at the pore surface, the Stern-layer surface-charge density rather than the equivalent Stern-layer ion concentration should be assumed to be continuous between the two cylindrical pores. Therefore, we include the radius dependency of the equivalent ion density into the boundary condition by using the square root of the previously described cross-section ratio a:

Note that on the ion concentrations (eqs 37 and 38), the cross-section ratio a is not included into the boundary conditions. The equation system consisting of eqs (30), (31), (33) and (35) along with the boundary conditions in eqs (36), (37) and (38) is solved semi-analytically for the electric potential and the ion concentrations in both pores. We derive an eigenvalue equation which we solve numerically and continue calculating the solution analytically. A detailed description of the solution can be found in the Supplementary Material. Our semi-analytic approach differs from the model by Bücker & Hördt (2013a) not only with respect to the incorporation of the Stern layer as an additional ion species. It also does not require a limitation of the valid frequency range to low frequencies – one of the simplifications made in earlier analytic membrane polarization models (Marshall & Madden 1959; Bücker & Hördt 2013a).

To calculate the resulting effective conductivity of the 1D pore system, we insert the solution for the ion concentration into the corresponding Nernst–Planck equation to obtain the current density. The conductivity can then be determined by integrating the modelled current densities of the ion species over the length of the two pores, dividing it by the applied electric potential V, and multiplying the result with the total length of the pore system:

The resulting conductivty |$\sigma$| is a complex frequency-dependent quantity and can therefore be written as

with |$\sigma ^{\prime }$| being the real and |$\sigma ^{\prime \prime }$| the imaginary part of the conductivity.

3.3 Numerical model

To validate the new semi-analytical model, we compare the resulting conductivity spectra with a numerical model first presented by Bücker et al. (2019). The 3D model consists of two rotationally symmetric, cylindrical pores, representing the pore system (Fig. 4). The interior of the cylinders is assumed to be filled with a monovalent, symmetric electrolyte described by the respective ion concentrations of anions and cations. A surface charge density on the pore surface represents the Stern layer.

The narrow pore is enclosed by a non-conducting region representing the matrix material. In this region, the electric potential must fulfill the Laplace equation

The solution for the electric potential and the ion concentrations inside the pores are calculated in two steps: First, the static solution is calculated by solving eq. (26). As a boundary condition, the equilibrium electric potential on the surface is set to be equal to the d-plane potential, that is approximated by the Zeta potential (Bücker et al. 2019). Via eq. (23) the ion distribution in the diffuse layer and the free electrolyte is computed from the solution for the potential. For the Stern layer, the equilibrium surface-charge density is |$\Sigma _\mathrm{S}^{(0)}$|⁠.

With the static solution at hand, the frequency-dependent case of a harmonic external electric field applied at the boundaries can be calculated by solving eqs (24) and (27) for the ion concentrations and the electric potential, respectively, and eq. (29) for the Stern-layer surface-charge density.

The non-conducting matrix material, the Stern layer, and the electrolyte are coupled by using suitable boundary conditions as discussed in Bücker et al. (2019).

To calculate the effective conductivity of the pore system, both the contribution of the electrolyte and the Stern layer have to be taken into account. The total current density in the electrolyte is calculated by integrating the anion and cation current densities over the pore’s cross-section. For the Stern layer, we calculate the current density at one end of the pore system from the Nernst–Planck equation eq. (29). Both contributions are added and then normalized to the applied electric field to obtain the conductivity of the 3D model of the cylindrical pore sequence.

4 MODEL COMPARISON

4.1 Spectral response

We compare the spectra obtained from the new semi-analytic model and the numerical model by Bücker et al. (2019) to the response of the “old analytical” model by Bücker & Hördt (2013a). For the latter, we also use the correction for the Stern-layer cations, which Bücker & Hördt (2013a) in their largely simplifying approach simply add to the diffuse-layer cations. If not stated otherwise, the parameters used for the simulations are listed in Table 4.

Figs 5(a) and (c) show example spectra of the normalized complex conductivity |$\sigma (\omega )/\sigma _0$|⁠, where

is the bulk conductivity of the electrolyte. Fig. 5(b) shows the real part of the normalized complex conductivity relative to its low-frequency limit: |$\sigma (\omega )/\sigma ^{\prime }(\omega \rightarrow 0)$|⁠.

The real part of the normalized complex conductivity (Figs 5a and b) is different for all three models, but varies only over a rather small range. Normalizing the real part to the respective low-frequency limits illustrates that the overall spectral behaviour of the new semi-analytical model and the numerical model show a similar frequency dependence, while the old analytical model deviates.

For the imaginary part (Fig. 5c), the new semi-analytic model shows a fair agreement with the numerical simulation around the peak at lower frequencies, whereas there are significant differences at higher frequencies. Our semi-analytic model does not explicitly include Maxwell–Wagner polarization of the modelled material, which is caused by a difference in bulk electric properties. However, our semi-analytic approach involving two different types of 1D pores leads to an artificial difference in the electrical conductivity at the boundary between these pores. Therefore, our semi-analytic model exhibits a high-frequency peak, which is a mathematical artifact and should not be confused with the physical process of Maxwell–Wagner polarization in a real 3D pore sequence. In contrast, the numerical model includes a pore constriction in the form of a matrix material narrowing the pore channel. By doing that, the electrical properties of the matrix are included in the numerical model and therefore, the high-frequency peak can be interpreted as Maxwell–Wagner polarization (Bücker et al. 2019). Since the high-frequency increase of the new semi-analytic model is not related to a real physical process, we define an interpretable frequency range (white area in Fig. 5) as the regime where the low-frequency peak is not masked by the high-frequency increase and can clearly be distinguished from it.

In comparison to the other models, the old analytical model only produces a polarization response where the polarization peak is smaller and shifted to significantly higher frequencies for the chosen parameter set. Since the contribution at higher frequencies is omitted in the old analytical model, no additional high-frequency peak is visible in the spectra produced by that model.

The imaginary conductivity spectra for other parameter sets are displayed in Fig. 6. For spectra 1 and 3 (i.e. low d-plane potentials), the new semi-analytical model is able to reproduce the numerical model and exhibits the same non-symmetric spectral shape, while the old analytic model shows significant deviations. For spectrum 2 (i.e. high d-plane potential and Stern-layer surface-charge density), both the old and the new semi-analytic model are able to reproduce the primary peak at low frequencies, but fail to reproduce the behaviour at frequencies at around |$\sim 10^3$| rad s|$^{-1}$|⁠, where the numerical model shows a secondary peak. For spectrum 4 (i.e. a low d-plane potential and a high Stern-layer surface-charge density), the response of the new semi-analytic model merges into that of the old analytic model. Compared to the numerical model, the polarization magnitude is overestimated whereas the peak frequency is fairly well reproduced.

4.2 Maximum imaginary conductivity and relaxation time

To compare our new semi-analytic model with the numerical model and the old analytical model over a wider parameter range, we calculate the maximum imaginary conductivity |$\sigma ^{\prime \prime }_\mathrm{max}$| and the related relaxation time. We approximate the latter with the inverse of the peak angular frequency of the imaginary conductivity, that is |$\tau = \left(\omega (\sigma ^{\prime \prime }_\mathrm{max})\right)^{-1}$|⁠. If no low-frequency peak is detectable, no value is given for the chosen parameter set. Since our new semi-analytic model and the numerical model include the polarization of both the Stern layer and the diffuse layer (i.e. membrane polarization), two distinct peaks related to either process may eventually occur in the low-frequency range (e.g. Fig. 6b). In that case, we consider the peak with the higher maximum imaginary conductivity to be the primary peak.

Fig. 7 shows the maximum imaginary conductivity and the respective relaxation time for the numerical model (a and b), the new semi-analytic model (c and d) and the old analytic model (e and f) for varying d-plane potentials and Stern-layer surface-charge densities. All three models show similar patterns in terms of the maximum imaginary conductivity (Figs 7a, c and e) and cover the same range of values in the studied parameter space. For low d-plane potential magnitudes, the numerical model and the new semi-analytic model agree well, while the old analytic model shows slightly higher maximum imaginary conductivities for high Stern-layer surface-charge densities. At higher d-plane potential magnitudes, the new semi-analytic model rather fits the old analytic model with an approximately constant maximum imaginary conductivity for a changing surface-charge density of the Stern layer, while the numerical model shows some slight dependence on the Stern-layer surface-charge density. This behaviour can also be observed in the example spectra in Fig. 6: Spectra 1 and 3 show the same maximum imaginary conductivity for the numerical and the new semi-analytic model, while for Spectrum 4, the peak amplitude of the new semi-analytic model rather matches that of the old analytic model.

For the numerical model, the calculated relaxation time (Fig. 7b) exhibits two regimes: for low d-plane potential magnitudes, the relaxation time strongly varies with the surface-charge density of the Stern layer, while for sufficiently high d-plane potential magnitudes, no significant variation is apparent. This behaviour is related to the relative importance of Stern-layer and membrane polarization, that is which process dominates the overall polarization response. For low d-plane potential magnitudes, the dominating process is the polarization of the Stern layer, while for higher d-plane potential magnitudes, the relaxation time is dominated by a membrane-polarization process. The new semi-analytic model shows a similar behaviour, although the transition between the two regimes happens at slightly lower d-plane potential magnitudes. The old analytic model shows no strong variation in the relaxation time, since Stern-layer polarization is not included in the model and therefore membrane polarization is the dominating process in the whole parameter space. This behaviour is also apparent in the spectra displayed in Fig. 6: While we see a shift in the peak frequency for the numerical model and the new semi-analytic model, the peak occurs in the same frequency range for the old analytic model. Only in the parameter space where the polarization process is dominated by membrane polarization, all three models predict the same relaxation time.

Fig. 8 shows the relative difference of the maximum imaginary conductivity

and the logarithmic difference of the relaxation time

of both the new semi-analytic model and the old analytic model compared to the results of the numerical model. Although both models exhibit deviations compared to the numerical model, the differences are significantly larger for the old analytical model. In particular, the prediction of the relaxation time at low d-plane potential magnitudes is improved by the new semi-analytic model.

Overall, the new semi-analytic model is able to reproduce the main features in terms of polarization time and maximum imaginary conductivity, while the old analytic model only agrees well when the polarization process is dominated by membrane polarization. Therefore, the new semi-analytic model is a significant improvement, particularly in the regime dominated by Stern-layer polarization.

5 RESULTS AND DISCUSSION

5.1 Polarization response for different mineral surfaces

We investigate the polarization response of a montmorillonite and a quartz surface predicted by the new semi-analytic model. We vary the NaCl concentration along with the permanent negative surface-charge density due to isomorphic substitution (montmorillonite) or the pH (quartz), respectively.

In a first step, we use our surface-complexation models to calculate the Stern-layer surface-charge density and the d-plane potential for each parameter combination. The results are displayed as isolines in Fig. 9. We then use the resulting EDL parameters to calculate the IP response of the new semi-analytic membrane-polarization model. To analyse the results, we consider the maximum imaginary conductivity and the relaxation time as previously defined.

We will use the same pore geometry for both types of surface (see Table 3) in order to investigate the effect of the different surface-complexation models separately from the effect of geometry. Additionally, we will use the same pore geometry as in Bücker et al. (2019) to make our results comparable to previous results. However, the choice of the same pore geometry for quartz and clay is not realistic in comparison to realistic geologic materials. While the chosen length scales are reasonable for quartz grains (e.g. Koch et al. 2011; Weller et al. 2011), the microstructure of clay typically exhibits smaller length scales (e.g. Leroy et al. 2015; Tournassat et al. 2015). We will discuss our results considering more realistic pore geometries for montmorillonite in Section 5.4.

5.1.1 Montmorillonite surface

For the basal surface of montmorillonite (Figs 9a and c), the Stern-layer surface-charge density mostly depends on the surface-charge density at the clay surface. The d-plane potential magnitude is largely controlled by the NaCl concentration but also shows a slight dependence on the surface-charge density. In the parameter space considered here, the maximum imaginary conductivity strongly varies with the d-plane potential, exhibiting the largest values for high d-plane potential magnitudes. This is due to higher d-plane potential magnitudes corresponding to a higher charge density in the diffuse layer, leading to a larger polarization effect (e.g. Hördt et al. 2016).

The relaxation time decreases with increasing surface-charge density of the Stern-layer, especially for high NaCl concentrations and thus low d-plane potential magnitudes. A similar effect is described by the Stern-layer polarization model by Lyklema et al. (1983), where the relaxation time of the Stern layer depends on the parameter M:

where R is the radius of the polarizing grain. Because M increases with increasing Stern-layer surface-charge density (eq. 34), the relaxation time |$\tau _\mathrm{ L}$| decreases accordingly. This indicates, that for montmorillonite, Stern-layer polarization is the dominant polarization process due to the comparably high partition coefficients of around 70 per cent in the studied parameter space (see Fig. 2).

For higher d-plane potential magnitudes, that is at lower NaCl concentrations, the dependence on the Stern-layer surface-charge density becomes weaker. This can be explained by the contribution of membrane polarization becoming more important for the overall polarization response, which is, however, still dominated by Stern-layer polarization. Consequently, in our membrane polarization model, Stern-layer polarization is the more important polarization process for clay minerals for the chosen pore lengths and radii. Note that the pore geometry might have an influence on this result, as it also has a strong effect on the polarization behaviour.

5.1.2 Quartz surface

For a quartz surface (Figs 9b and d), the surface-charge density of the Stern layer strongly increases with increasing pH and with increasing NaCl concentration, while the d-plane potential magnitude increases with an increasing pH and decreases with an increasing NaCl concentration. As it is the case for montmorillonite (Fig. 9a), the maximum imaginary conductivity increases for an increasing d-plane potential magnitude due to the larger number of ions in the diffuse layer. We can identify two regimes in the relaxation time (Fig. 9d): For an ion concentration smaller than 10 mol m|$^{-3}$|⁠, the relaxation time is mostly controlled by the d-plane potential, indicating that membrane polarization is the dominant process. For ion concentrations larger than 10 mol m|$^{-3}$|⁠, the relaxation time is dominated by the surface-charge density of the Stern layer, indicating that Stern-layer polarization is the dominant process. In both regimes, a decrease of the relaxation time for increasing pH can be observed. For the membrane-polarization regime, this is in accordance with the results from Hördt et al. (2016) (their fig. 14). For the regime dominated by Stern-layer polarization, the pH dependency of the relaxation time is even more pronounced.

We now want to investigate the transition between these two regimes more closely. For that purpose, we consider selected spectra of the imaginary conductivity for a quartz surface at pH 5.5 for different NaCl concentrations around 10 mol m|$^{-3}$|⁠. The resulting spectra are displayed in Fig. 10. The picked peaks of the imaginary conductivity for each spectrum are indicated by black circles in Fig. 10. For concentrations above the threshold of 10 mol m|$^{-3}$|⁠, the spectrum is dominated by a low-frequency peak at |$\sim$|0.3 rad s|$^{-1}$| caused by Stern-layer polarization. For decreasing concentration, a second peak caused by membrane polarization at a higher frequency of |$\sim$|2 Hz builds up and becomes dominant below 10 mol m|$^{-3}$|⁠. As soon as the membrane polarization becomes dominant, the calculated relaxation time (i.e. the inverse of the peak frequency) shows a sudden increase, which leads to the transition between regimes visible in Fig. 9(d).

This observation indicates that for quartz surfaces, depending on the ion concentration, both membrane polarization and polarization of the Stern layer make a significant contribution to the overall polarization response. In addition, pH also influences the respective contribution of membrane polarization and Stern-layer polarization to the overall polarization. Note that at the limit between the two regimes, both processes have a similar magnitude. Therefore, the maximum imaginary conductivity in Fig. 9(b) exhibits a smooth transition around the threshold concentration.

5.2 Influence of pore geometry

Besides the structure of the EDL, also the pore geometry, that is the pore lengths and radii, have a large influence on the IP-response. We therefore run a parameter study varying the geometric parameters of the pore system and compare our results to those Hördt et al. (2017) obtained for the old analytic model.

When changing the length of the pores, Hördt et al. (2017) observed two regimes for the relaxation time (their figs 3 and 4): For sufficiently short wide pores, the relaxation time quadratically depends on the length of the wide pore, while for longer wide pores, only a quadratic dependence on the length of the narrow pore can be observed. This behaviour can be explained as the relaxation mainly takes place in one of the two pores, leading to a relaxation time that only depends on one pore length. These two regimes are called long and short narrow regime, depending on the pore dominating the relaxation process (Bücker & Hördt 2013b). Our new semi-analytic model shows a similar behaviour and is therefore not further discussed (see Supplementary Material for the corresponding figures).

Fig. 11 shows the maximum imaginary conductivity (a and c) and the relaxation time (b and d) for the new (a and b) and the old analytic model (c and d). Note that results for |$r_2 > r_1$| are not shown, because for those parameter combinations, the narrow and the wide pore simply switch their roles. For the new semi-analytic model, no value is given for high radii of the wider pore, because in that regime, the polarization peak is masked by the high-frequency polarization and can therefore not be identified. We observe a similar behaviour of both models for most parameters, indicating that membrane polarization is the dominant process in the parameter space studied. Only for small radii, the relaxation times of the new semi-analytic model deviate from those of the old analytic model, because Stern-layer polarization becomes the dominating process. Hence, not only the structure of the EDL, but also the pore geometry, or more precisely, the radii of the two pores, influence the relative importance of the two competing polarization processes.

5.3 Comparison to grain-based models

Due to their complex geometrical structure, natural geologic materials are only represented in a highly simplified manner by both, pore-based models like those discussed in this study and grain-based models considering spherical particles (e.g. Schwarz 1962; Lyklema et al. 1983; Leroy et al. 2008; Bücker et al. 2019). A comparison to grain-based models can help identifying differences and similarities of the predictions of both model families.

The IP relaxation time for non-conducting geologic materials typically exhibits a quadratic dependence on a characteristic length scale. Just as existing models (e.g. Volkmann & Klitzsch 2010; Bücker & Hördt 2013a, b), for our new semi-analytic model this length scale is related to the pore length of either the wide or the narrow pore. In contrast, the relaxation time of grain-based models is controlled by the grain radius (e.g. Schwarz 1962; Dukhin & Shilov 1974; Fixman 1980; Lyklema et al. 1983).

Lyklema et al. (1986) and de Lima & Sharma (1992) observed that the diffuse layer polarization model involving polarization of neutral ion clouds outside the EDL better reproduces low-frequency dielectric dispersion measurements on colloidal particles than a Stern-layer polarization model. Delgado et al. (1998) showed that the surface diffusion mechanism associated with Stern-layer polarization dominates the polarization response compared to the volume-diffusion mechanism associated with diffuse-layer polarization when particle concentration increases in such a way that the space allowing ion cloud formation is not available anymore.

Studies of the coupled polarization of Stern layer and diffuse layer (Lesmes & Morgan 2001; Bücker et al. 2019) show that for grain-based models, typically the polarization of the Stern layer is the dominating polarization process as it exhibits a larger magnitude (de Lima & Sharma 1992), while the polarization of the diffuse layer only plays a minor role. In our new semi-analytic model, the polarization of the Stern layer and the membrane polarization can lead to comparably large polarization amplitudes (see Fig. 10). Whether one or the other process dominates depends on the structure of the EDL, which in turn depends on electrolyte concentration and pH, and pore geometry. Therefore, Stern-layer polarization should no longer be neglected in membrane-polarization modelling.

5.4 Comparison to experimental data

Now we want to compare the predictions of our new semi-analytic model to some experimental results from previous studies. Mendieta et al. (2021) performed spectral induced-polarization measurements on unconsolidated clays with varying fluid conductivity. They observed an increasing imaginary conductivity with fluid conductivity at low values of fluid conductivity and a decreasing imaginary conductivity with fluid conductivity at higher fluid conductivity. A similar behaviour was already observed in the modelling results of Hördt et al. (2016) for membrane polarization with sufficiently small radii.

For the comparison, we use a pore geometry similar to the microstructure of montmorillonite. For the pore constriction, we use a pore cross-section of approximately four Debye lengths, which is a typical distance between two montmorillonite particles (Leroy et al. 2015). Note that the Debye length changes with fluid conductivity (eq. 2). For the wide pore, we assume a pore cross-section which is 2.5 times larger than that of the narrow pore (i.e. ten Debye lengths). Typically, montmorillonite clay platelets have a length between 50 and 1000 nm (Tournassat et al. 2015). We therefore choose a length of 500 nm for the narrow pore, such that the pore constriction represents the space between two montmorillonite particles. For the larger pore, we chose a length of 5 |$\mu$|m.

For this clay-like pore geometry, we calculate the maximum imaginary conductivity for different surface-charge densities and varying ion concentration, which are displayed in Fig. 12. The results of the full parameter space (i.e. for varying |$\Sigma _0$| and |$c_\infty$| as shown in Fig. 9 for the larger pore geometry) are given in the Supplementary Material. The maximum imaginary conductivity shows increasing values at low ion concentrations and decreasing values at higher ion concentrations. The concentration of the peak tends to higher concentrations for larger surface-charge densities. This peak-like behaviour of the maximum imaginary conductivity is similar to the one observed by Mendieta et al. (2021).

For glass beads, Leroy et al. (2008) observe a decreasing maximum phase shift with increasing fluid salinity and Weller et al. (2011) see a similar behaviour for the chargeability of sandstones. In agreement with these experimental observations, our model predicts a decreasing maximum imaginary conductivity (Fig. 9b) for quartz surfaces. For sandstone Weller et al. (2011) report a constant or slightly decreasing relaxation with increasing fluid salinity. Fig. 9(d) shows such a behaviour for some parameter combinations. However, the variation of the relaxation time in Fig. 9(d) is highly dynamic and therefore no clear prediction of the direction of change of the relaxation time with fluid conductivity can be made.

5.5 Limitations of the model

The new semi-analytic model exhibits various limitations. Although the predicted main relaxation time and the amplitude of the imaginary-conductivity peak of the new semi-analytic model are in agreement with those of the numerical model, the exact spectral response including minor peaks can not be reproduced well for some parameter sets. This limitation is probably due to the use of the effective concentration of Stern-layer cations instead of describing the Stern-layer cations in terms of a surface-charge density.

To achieve a suitable description of Stern-layer polarization, the mobility of the cations in the Stern layer is an important parameter. However, this value is not well constrained (Weller et al. 2013). Modelling of the surface conductivity of clay materials indicates that the Stern-layer cation (e.g. |$\mathrm{Na^+}$|⁠) mobility might be several orders of magnitude lower than the cation mobility in the bulk electrolyte (Revil 2012; Weller et al. 2013), whereas cation (e.g. |$\mathrm{Na^+}$|⁠) mobilities in the Stern layer and the bulk electrolyte are more likely to have similar values for silica minerals (Leroy et al. 2008; Revil 2012). In the case of dispersed clay suspensions, Leroy et al. (2017) successfully reproduced measured low-frequency complex conductivity spectra considering that |$\mathrm{Na^+}$| mobility in the Stern layer of montmorillonite is half its value in bulk electrolyte. Compared to Leroy et al. (2017), increasing clay-particle concentration in the experiments used in Revil (2012) may explain the very low-effective |$\mathrm{Na^+}$| mobility in the Stern layer found by Revil (2012) due to interacting clay particles.

Furthermore, the spectral response at high frequencies can not be modelled suitably, limiting the model’s interpretable frequency range to the frequency range, at which the low-frequency polarization is not masked by the steep high-frequency increase in imaginary conductivity.

At this stage, our model is limited to a monovalent, symmetric electrolyte and the surface-complexation model is only valid for a NaCl solution. Since naturally occuring pore fluids contain many different ion species with some of them having higher valences, this limitation restricts the practical applicability of our model to the modelling of realistic geological materials. We do not take the contribution of |$\mathrm{H^+}$|-ions to the overall conductivity into account, which might cause deviations at low pH (Skold et al. 2011).

Additionally, our model describes the pore channel as a system consisting of two types of cylindrical pores. This represents a strong simplification compared to the pore space of real geologic materials. In addition, our model can only account for two pore radii, which is in contrast to the wide range of pore radii occuring in natural rocks. Therefore, a suitable upscaling technique is necessary to model realistic geologic materials. Such models have been proposed, for example by Stebner et al. (2017) or Maineult et al. (2017) in the form of random networks of different pores, where the overall polarization response is calculated from Kirchhoff’s law (Kirchhoff 1845). An alternative hybrid approach is presented by Niu et al. (2020), where small-scale polarization models are convoluted with the pore-size distribution of a rock sample. The result is used to calculate an effective fluid conductivity, which is then used to predict the overall response of the rock geometry taken from digital images. Our new model can be used with either upscaling method, of course always considering the limitations discussed above.

To this point, our model relates the difference in transference numbers causing the polarization effect to a difference in pore radii. Chuprinko & Titov (2017) studied a model consisting of two pores of the same radius, but exhibiting different EDL structures, which they relate to different mineral surfaces. Their work indicates that such a situation might lead to another membrane polarization mechanism in addition the polarization of the Stern layer and membrane polarization caused by pore constrictions (Chuprinko & Titov 2017). An extension of our model to different mineral surfaces in the two pores might enable modelling this kind of polarization, too.

6 CONCLUSION

We present a new model for membrane polarization that does not only include the polarization of the diffuse layer but also the polarization of the cations in the Stern layer. This extension is implemented by considering a third Nernst–Planck differential equation describing the Stern-layer cations besides those for cations and anions in the electrolyte and the diffuse layer. Comparison with numerical simulations using the model by Bücker et al. (2019) shows a good agreement and a great improvement relative to the analytical model by Bücker & Hördt (2013a).

We use surface-complexation models for both a clay and a silica surface to model the structure of the EDL for different chemical conditions in the electrolyte solution (NaCl concentrations, pH). We then use these results in our polarization model: For montmorillonite, in the studied parameter space of surface charge densities between |$-0.1$| and |$-0.15$| C m|$^{-2}$| and NaCl concentrations between 1 and 100 mol m|$^{-3}$|⁠, the overall polarization response is dominated by Stern-layer polarization. In the case of a quartz surface, two regimes can be identified that correspond to a dominant Stern-layer polarization (for NaCl concentrations below |$\sim 10$| mol m|$^{-3}$|⁠) and a dominant polarization of the diffuse layer (for NaCl concentrations above |$\sim 10$| mol m|$^{-3}$|⁠), respectively. It is noteworthy that both polarization processes contribute with similar magnitudes in the pore constriction geometry, which differs significantly from the relative contributions of Stern- and diffuse-layer polarization predicted by grain-based models.

Beside the structure of the EDL, as in earlier membrane-polarization models, the pore geometry has a strong influence on the predicted polarization behaviour. Overall, our new semi-analytic model shows a similar dependence on geometrical parameters (pore lengths and radii) as the model of Bücker & Hördt (2013a). The relaxation time systematically depends on the length of either the wider or the narrow pore, depending on which pore dominates the relaxation process. Similar to the structure of the EDL, the pore radius has an influence on the respective contribution of Stern-layer polarization and the polarization of the diffuse layer: For small radii (below |$5 \times 10^{-8}$| nm for the wide and below |$5 \times 10^{-9}$| nm for the narrow pore), Stern-layer polarization is the dominant process, while for larger radii, membrane polarization dominates the overall polarization response.

Since the inclusion of Stern-layer polarization in the membrane-polarization model is a significant improvement compared to previous analytic models, our model allows further investigation of the polarization taking place at constrictions of electrolyte-filled pores. Future investigations of the modeled polarization response will give deeper insight into the processes taking place at the pore scale.

ACKNOWLEDGMENTS

PL warmly thanks BRGM for funding a one-month stay of D. Kreith in Orléans (2023 June). We thank L. Slater and D. Jougnot for their insightful comments, which contributed significantly to the improvement of the manuscript. We acknowledge support by the Open Access Publication Funds of Technische Universität Braunschweig. The work of Philippe Leroy was also supported by the IMAGE project, funded by the French National Research Agency (grant ANR-21-CE04-0013).

DATA AVAILABILITY

Data and model files associated with this research are available at Zenodo (https://doi.org/10.5281/zenodo.13929230).

REFERENCES

APPENDIX A: CALCULATION OF SURFACE CHARGE DENSITIES

A1 Montmorillonite

The equilibrium constant of |$\mathrm{Na}$|-adsorption can be calculated as

with |$\Gamma _i$| being the surface-site density of species i. In particular, |$\Gamma _{\mathrm{X}^-}$| is the surface-site density of the clay sites |$\gt X^-$| without a sodium cation attached. The total surface-charge density is connected to the surface-site densities by

Inserting eq. (A2) into eq. (A1) and rearranging leads to

and with

we obtain eq. (9).

A2 Quartz

The two equilibrium constants for protons and sodium cations can be expressed in terms of the surface-site densities |$\Gamma _i$| as

and

These equations can be rearranged to

and

The total surface-site density is

with A defined in eq. (16). We can therefore make the substitution

in eqs (A7) and (A8). Now the surface-charge densities at the 0-plane can be expressed as

and at the Stern plane as

yielding eqs (14) and (15).