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

D Kreith, Philippe Leroy, M Bücker

► To cite this version:

D Kreith, Philippe Leroy, M Bücker. A new semi-analytic model for Stern-layer polarization in pore throats. Geophysical Journal International, 2024, 239 (3), pp.1910 - 1927. 10.1093/gji/ggae370. hal-04788967

HAL Id: hal-04788967

https://brgm.hal.science/hal-04788967v1

Submitted on 18 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A new semi-analytic model for Stern-layer polarization in pore throats

D. Kreith<sup>1</sup>,<sup>1</sup> P. Leroy<sup>2</sup> and M. Bücker<sup>1,3</sup>

Accepted 2024 October 9. Received 2024 May 2; in original form 2024 September 10

## 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.

**Key words:** Electrical properties; Microstructure; Induced polarization; Numerical modelling; Hydrogeophysics.

## 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

<footnote n="1"><sup>1</sup>Institute of Geophysics and Extraterrestrial Physics, TU Braunschweig, Mendelssohnstr. 3, 38106 Braunschweig, Germany.
E-mail: d.kreith@tu-braunschweig.de</footnote>

<footnote n="2"><sup>2</sup>BRGM, French Geological Survey, 3 Avenue Claude Guillemin, 45060 Orléans Cedex 2, France</footnote>

<footnote n="3"><sup>3</sup>Institute of Geosciences, Christian-Albrechts-Universität zu Kiel, Otto-Hahn-Platz 1, 24118 Kiel, Germany</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Schematic representation of the basic Stern model of (a) the basal surface of montmorillonite (after Leroy *et al.* 2015) and (b) the main quartz surface, quartz (0001) crystal surface (after Leroy *et al.* 2013) in contact with a binary monovalent electrolyte. The surface-charge density Σ₀ of the mineral surface at *x* = *x*₀ is partly counterbalanced by the surface-charge density Σₛ of the Stern layer (region 1) at *x* = *x*ₛ (determined by the capacitance *C*₁). The electric potential drops from φ₀ at the surface to φ₉ at the Stern plane. The remaining charge is counterbalanced by the surface-charge density Σ₉ of the diffuse layer (region 2). The characteristic length scale of the diffuse layer is the Debye length λ₉, which is determined by the exponential drop of the electric potential in the diffuse layer. Outside the double layer, the electrolyte is undisturbed and the electric potential approaches zero (region 3).</figcaption></figure>

(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 Σₛ and the potential φ₉ 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
<table><caption>**Table 1.** List of all symbols (latin characters) used in alphabetical order.</caption><thead><tr><th>Symbol</th><th>Meaning</th><th>Unit</th></tr></thead><tbody><tr><td>a</td><td>Pore cross-section ratio</td><td>dimensionless</td></tr><tr><td>a<sub>i</sub></td><td>Ion activity</td><td>mol L<sup>-1</sup></td></tr><tr><td>b<sub>p,n</sub></td><td>Effective ion-mobility ratio</td><td>dimensionless</td></tr><tr><td>c<sub>0</sub></td><td>Bulk ion concentration (surface-complexation model)</td><td>mol L<sup>-1</sup></td></tr><tr><td>c<sub>∞</sub></td><td>Bulk ion concentration (polarization model)</td><td>mol m<sup>-3</sup></td></tr><tr><td>c</td><td>Ion concentration (polarization model)</td><td>mol m<sup>-3</sup></td></tr><tr><td>C<sub>1</sub></td><td>Capacitance between mineral surface and Stern plane</td><td>F m<sup>-1</sup></td></tr><tr><td>D<sub>i</sub></td><td>Ion diffusion coefficient</td><td>m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>e</td><td>Elementary charge</td><td>C</td></tr><tr><td>F</td><td>Faraday's constant</td><td>A s mol<sup>-1</sup></td></tr><tr><td>f<sub>Q</sub></td><td>Partition coefficient</td><td>dimensionless</td></tr><tr><td>j<sub>i</sub></td><td>Current density</td><td>A m<sup>-2</sup></td></tr><tr><td>k<sub>B</sub></td><td>Boltzmann constant</td><td>J K<sup>-1</sup></td></tr><tr><td>K<sub>i</sub></td><td>Equilibrium constant</td><td>dimensionless</td></tr><tr><td>L<sub>i</sub></td><td>Pore length</td><td>m</td></tr><tr><td>M</td><td>Correction coefficient (Lyklema et al. 1983)</td><td>dimensionless</td></tr><tr><td>R</td><td>Particle radius</td><td>m</td></tr><tr><td>r</td><td>Spatial coordinate (3D)</td><td>m</td></tr><tr><td>r<sub>i</sub></td><td>Pore radius</td><td>m</td></tr><tr><td>r<sub>S</sub></td><td>Spatial coordinate (along surface)</td><td>m</td></tr><tr><td>T</td><td>Absolute temperature</td><td>K</td></tr><tr><td>t</td><td>Time</td><td>s</td></tr><tr><td>U</td><td>Electric potential (polarization model)</td><td>V</td></tr><tr><td>V</td><td>Initial electric potential (polarization model)</td><td>V</td></tr><tr><td>x</td><td>Spatial coordinate (1D)</td><td>m</td></tr><tr><td>x<sub>0</sub></td><td>Coordinate of the mineral surface</td><td>m</td></tr><tr><td>x<sub>S</sub></td><td>Coordinate of the Stern plane</td><td>m</td></tr><tr><td>z<sub>i</sub></td><td>Ion valence</td><td>dimensionless</td></tr></tbody></table>

<table><caption>**Table 2.** List of all symbols (greek characters) used in alphabetical order.</caption><thead><tr><th>Symbol</th><th>Meaning</th><th>Unit</th></tr></thead><tbody><tr><td>γ<sub>Na</sub></td><td>Sodium ion activity coefficient</td><td>dimensionless</td></tr><tr><td>Γ<sub>i</sub></td><td>Surface site density</td><td>m<sup>-2</sup></td></tr><tr><td>ε<sub>0</sub></td><td>Permittivity of vacuum</td><td>F m<sup>-1</sup></td></tr><tr><td>ε<sub>r</sub></td><td>Relative permittivity</td><td>dimensionless</td></tr><tr><td>κ</td><td>Inverse Debye length</td><td>m<sup>-1</sup></td></tr><tr><td>λ<sub>d</sub></td><td>Debye length</td><td>m</td></tr><tr><td>μ<sub>i</sub></td><td>Ion mobility</td><td>m<sup>2</sup> V<sup>-1</sup> s<sup>-1</sup></td></tr><tr><td>σ</td><td>Total complex conductivity</td><td>S m<sup>-1</sup></td></tr><tr><td>σ'</td><td>Real conductivity</td><td>S m<sup>-1</sup></td></tr><tr><td>σ''</td><td>Imaginary conductivity</td><td>S m<sup>-1</sup></td></tr><tr><td>σ''<sub>max</sub></td><td>Maximum imaginary conductivity</td><td>S m<sup>-1</sup></td></tr><tr><td>Σ<sub>0</sub></td><td>Surface-charge density of the mineral surface</td><td>C m<sup>-2</sup></td></tr><tr><td>σ<sub>0</sub></td><td>Bulk electrolyte conductivity</td><td>S m<sup>-1</sup></td></tr><tr><td>Σ<sub>d</sub></td><td>Surface-charge density of the diffuse layer</td><td>C m<sup>-2</sup></td></tr><tr><td>Σ<sub>S</sub></td><td>Surface-charge density of the Stern layer</td><td>C m<sup>-2</sup></td></tr><tr><td>τ</td><td>Relaxation time</td><td>s</td></tr><tr><td>φ</td><td>Electric potential (surface-complexation model)</td><td>V</td></tr><tr><td>φ<sub>d</sub></td><td>Electric potential at the d-plane</td><td>V</td></tr><tr><td>ω</td><td>Angular frequency</td><td>rad s<sup>-1</sup></td></tr></tbody></table>

surfaces to exhibit a negative surface-charge density $\Sigma_0 < 0$ at the 0-plane ($x_0$). At a certain distance ($x_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_S$. At distances beyond the Stern plane ($x > x_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_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_d$. Electro-neutrality implies that the surface-charge density of the Stern plane $\Sigma_S$ and of the diffuse layer $\Sigma_d$ cancel out the surface-charge density at the 0-plane $\Sigma_0$, i.e.

$$ \Sigma_0 + \Sigma_S + \Sigma_d = 0. \qquad (1) $$

The Debye length

$$ \lambda_d = \sqrt{\frac{\epsilon_0 \epsilon_r k_B T}{2 \times 10^3 e c_0 F}} \qquad (2) $$
<table><caption>**Table 3.** Model parameters used for surface-complexation modelling, if not stated otherwise.</caption><thead><tr><th>Parameter</th><th>Symbol [unit]</th><th>Value</th></tr></thead><tbody><tr><td>Absolute temperature</td><td><i>T</i> [K]</td><td>293</td></tr><tr><td>Relative electrolyte permittivity</td><td><i>ε</i><sub>r</sub> [-]</td><td>80</td></tr><tr><td>Stern-layer capacitance (Mt)<sup>a</sup></td><td><i>C</i><sub>1</sub> [F m<sup>-2</sup>]</td><td>0.71</td></tr><tr><td>Stern-layer capacitance (quartz)<sup>b</sup></td><td><i>C</i><sub>1</sub> [F m<sup>-2</sup>]</td><td>3.43</td></tr><tr><td>Equilibrium constant (Na-adsorption at Stern plane, Mt)<sup>c</sup></td><td><i>K</i><sub>Na,Mt</sub> [-]</td><td>10<sup>0.88</sup></td></tr><tr><td>Equilibrium constant (Na-adsorption at Stern plane, quartz)<sup>b</sup></td><td><i>K</i><sub>Na,Qtz</sub> [-]</td><td>10<sup>0.58</sup></td></tr><tr><td>Equilibrium constant (protonation)<sup>b</sup></td><td><i>K</i><sub>H</sub> [-]</td><td>10<sup>7.28</sup></td></tr><tr><td>Surface-site density (quartz)<sup>d</sup></td><td><i>Γ</i><sub>S</sub> [m<sup>-2</sup>]</td><td>4.6 × 10<sup>18</sup></td></tr></tbody></table>

is a measure for the extension of the diffuse layer, with $ε_r$ being the relative permittivity of the electrolyte solution, $ε_0 = 8.85 × 10^{-12}$ 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$ A s mol$^{-1}$, the elementary charge $e = 1.602 × 10^{-19}$ C, Boltzmann's constant $k_B = 1.380 × 10^{-23}$ 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

$$ \Delta\varphi = 0 \qquad (3) $$

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

$$ \nabla\varphi|_{x_0} = -\frac{\Sigma_0}{\varepsilon_0\varepsilon_r}. \qquad (4) $$

(Bourg *et al.* 2007). Note that eq. (4) implies that the potential is constant within the mineral. The position $x_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

$$ x_S = \frac{\varepsilon_0\varepsilon_r}{C_1}. \qquad (5) $$

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

$$ \nabla\varphi|_{x_S} = -\frac{\Sigma_0 + \Sigma_S}{\varepsilon_0\varepsilon_r}. \qquad (6) $$

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

$$ \Delta\varphi = \frac{Fc_0}{\varepsilon_0\varepsilon_r} \left( \exp\left(\frac{e\varphi}{k_B T}\right) - \exp\left(-\frac{e\varphi}{k_B T}\right) \right) \qquad (7) $$

for the studied case of a monovalent, symmetrical 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 Na$^+$ at clay surface sites $> X^-$:

$$ > X^- + Na^+ \rightleftharpoons > X^- \cdots Na^+, \quad K_{Na,Mt}, \qquad (8) $$

where $K_{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$^{3+}$ and Fe$^{3+}$ ions by Mg$^{2+}$ and Fe$^{2+}$ 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:

$$ \Sigma_S = - \frac{K_{Na,Mt} a_{Na} \exp\left(-\frac{e\varphi_d}{k_B T}\right)}{1 + K_{Na,Mt} a_{Na} \exp\left(-\frac{e\varphi_d}{k_B T}\right)} \Sigma_0. \qquad (9) $$

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

$$ a_{Na} = \gamma_{Na} c_0 \qquad (10) $$

with $\gamma_{Na}$ being the sodium activity coefficient in the bulk electrolyte that can be calculated by using Davies' equation (Davies 1962):

$$ \log \gamma_{Na} = -\frac{1}{2} \left( \frac{\sqrt{c_0}}{1 + \sqrt{c_0}} - 0.3c_0 \right). \qquad (11) $$

### 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:

$$ > \text{SiO}^- + \text{H}^+ \rightleftharpoons > \text{SiOH}, \quad K_{\text{H}}, \qquad (12) $$

<footnote><sup>a</sup>The SCM for the basal surface of montmorillonite is not sensitive to $C_1$. We arbitrarily choose $C_1$ to relate it to a Stern-layer thickness of 1 nm using eq. (5)</footnote>

<footnote><sup>b</sup>Leroy *et al.* (2022)</footnote>

<footnote><sup>c</sup>Tournassat *et al.* (2009)</footnote>

<footnote><sup>d</sup>Hiemstra *et al.* (1989)</footnote>
with $K_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, Na-ions are adsorbed:

$$> \text{SiO}^{-} + \text{Na}^{+} \rightleftharpoons > \text{SiO}^{-} \cdots \text{Na}^{+}, \quad K_{\text{Na,Qtz}}. \quad (13)$$

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

$$\Sigma_0 = -e \left( 1 + K_{\text{Na,Qtz}} a_{\text{Na}} \exp \left( -\frac{e\varphi_d}{k_B T} \right) \right) \frac{\Gamma_S}{A} \quad (14)$$

and at the Stern plane

$$\Sigma_S = -e K_{\text{Na,Qtz}} a_{\text{Na}} \exp \left( -\frac{e\varphi_d}{k_B T} \right) \frac{\Gamma_S}{A} \quad (15)$$

with

$$A = 1 + K_H a_H \exp \left( -\frac{e\varphi_0}{k_B T} \right) + K_{\text{Na,Qtz}} a_{\text{Na}} \exp \left( -\frac{e\varphi_d}{k_B T} \right). \quad (16)$$

$\Gamma_S$ is the total surface-site density of > SiO⁻ and > SiOH-sites and is treated as an input parameter. While the activity of 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

$$a_H = 10^{-pH}. \quad (17)$$

## 2.2 Modelled double-layer structure

From our surface-complexation model, we extract the potential $\varphi_d$ at the d-plane and the surface-charge density $\Sigma_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

$$f_Q = \frac{\Sigma_S}{|\Sigma_0|}. \quad (18)$$

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 NaCl concentration $c_0$ in the electrolyte. The results for the quartz (0001) crystal surface at different 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

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Calculated d-plane potential (upper panel), Stern-layer surface-charge density (middle panel), and partition coefficient (lower panel) with varying NaCl concentrations at different surface-charge densities for the montmorillonite basal surface. Solid lines show the results of our numerical simulation, dots show the PHREEQC results.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Calculated d-plane potential (upper panel), Stern-layer surface-charge density (middle panel), and partition coefficient (lower panel) with varying pH at different salinities for a quartz (0001) crystal surface. Solid lines show the results of our numerical simulation, dots show the PHREEQC results. Model parameters are given in Table 3.</figcaption></figure>
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 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 $j_{p,n}$ of cations and anions are given by

$$ \mathbf{j}_{p,n}(\mathbf{r}, t) = -F D_{p,n} \nabla c_{p,n}(\mathbf{r}, t) \\ \mp F \mu_{p,n} c_{p,n}(\mathbf{r}, t) \nabla U(\mathbf{r}, t). \quad (19) $$

In this equation, $c_{p,n}$ (in mol m⁻³ = 10³ mol L⁻¹) is the ion concentration and *U* the electric potential (in V), $D_{p,n}$ (in m² s⁻¹) and $\mu_{p,n}$ (in m² V⁻¹ s⁻¹) 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_BT/e$. The vector **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

$$ \nabla \mathbf{j}_{p,n} + F \frac{\partial}{\partial t} c_{p,n}(\mathbf{r}, t) = 0. \quad (20) $$

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⁻¹):

$$ U(\mathbf{r}, t) = U^{(0)}(\mathbf{r}) + \delta U(\mathbf{r}, \omega) \exp(i\omega t), \quad (21) $$

$$ c_{p,n}(\mathbf{r}, t) = c_{p,n}^{(0)}(\mathbf{r}) + \delta c_{p,n}(\mathbf{r}, \omega) \exp(i\omega t), \quad (22) $$

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

$$ c_{p,n}^{(0)}(\mathbf{r}) = c_{\infty} \exp \left( \mp \frac{e}{k_B T} U^{(0)}(\mathbf{r}) \right). \quad (23) $$

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

$$ i\omega \delta c_{p,n}(\mathbf{r}, \omega) = \nabla\{D_{p,n} \nabla\delta c_{p,n}(\mathbf{r}, \omega) \\ + \mu_{p,n} [c_{p,n}^{(0)}(\mathbf{r}) \nabla\delta U(\mathbf{r}, \omega) \\ + \delta c_{p,n}(\mathbf{r}, \omega) \nabla U^{(0)}(\mathbf{r})]\}. \quad (24) $$

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:

$$ \nabla^2 U(\mathbf{r}, t) = -\frac{F}{\varepsilon_r \varepsilon_0} (c_p(\mathbf{r}, t) - c_n(\mathbf{r}, t)). \quad (25) $$

<figure><img src="image_5.png" /><figcaption>**Figure 4.** (a) Sketch of an electrolyte-filled pore channel with a negatively charged mineral surface. The darker grey area indicates the location of the pore constriction. (b) Corresponding cylindrical model with the red line indicating the Stern layer. (c) 1D-representation of the cylindrical model with effective values for ion concentration *c* and ion mobility *μ*.</figcaption></figure>

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

$$ \nabla^2 U^{(0)}(\mathbf{r}) = -\kappa^2 \frac{k_B T}{e} \sinh \left( \frac{e}{k_B T} U^{(0)}(\mathbf{r}) \right). \quad (26) $$

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

$$ \nabla^2 \delta U(\mathbf{r}, \omega) = -\frac{F}{\varepsilon_r \varepsilon_0} (\delta c_p(\mathbf{r}, \omega) - \delta c_n(\mathbf{r}, \omega)). \quad (27) $$

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

$$ \Sigma_S(\mathbf{r}_S, t) = \Sigma_S^{(0)} + \delta \Sigma_S(\mathbf{r}_S, \omega) \exp(i\omega t), \quad (28) $$

where $\Sigma_S^{(0)}$ is the equilibrium value and $\delta \Sigma_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:

$$ i\omega \delta \Sigma_S(\mathbf{r}_S, \omega) = \nabla_S(D_S \nabla_S \delta \Sigma_S(\mathbf{r}_S, \omega) \\ + \mu_S \Sigma_S^{(0)} \nabla_S \delta U_S(\mathbf{r}_S, \omega)). \quad (29) $$

The operator $\nabla_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 $2L_1$ and $2L_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 yields 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_S$ and $U$) rather than the symbols for the respective perturbation quantities ($\delta c_p$, $\delta c_n$, $\delta \Sigma_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):

$$i\omega c_p = D_p \frac{\partial^2 c_p}{\partial x^2} - \mu_p c_0 \frac{\partial^2 U}{\partial x^2}, \quad (30)$$

$$i\omega c_n = D_n \frac{\partial^2 c_n}{\partial x^2} + \mu_n c_0 \frac{\partial^2 U}{\partial x^2}. \quad (31)$$

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}$):

$$c_S = \frac{2 \Sigma_S}{F r_1}, \quad \bar{c}_S = \frac{2 \Sigma_S}{F r_2}. \quad (32)$$

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):

$$i\omega c_S = M D_S \frac{\partial^2 c_S}{\partial x^2} - \mu_S c_{S,0} \frac{\partial^2 U}{\partial x^2}, \quad (33)$$

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

$$M = 1 + \frac{\kappa \Sigma_S}{2 F c_\infty \cosh\left(\frac{e \varphi_d}{2 \kappa_B T}\right)}. \quad (34)$$

with $\kappa$ being the inverse Debye length given by $\kappa = \lambda_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

$$\frac{\partial^2 U}{\partial x^2} = \frac{F}{\varepsilon_0 \varepsilon_r} (c_p - c_n + c_S). \quad (35)$$

Eqs (30), (31), (33) and (35) apply to the narrow pore as well. We use overlined quantities ($\bar{c}_p$, $\bar{c}_n$, $\bar{c}_S$ and $\bar{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 $\bar{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:

$$\begin{align}
j_p(L_1) &= a \bar{j}_p(-L_2), \nonumber \\
j_n(L_1) &= a \bar{j}_n(-L_2), \nonumber \\
j_S(L_1) &= a \bar{j}_S(-L_2). \tag{36}
\end{align}$$

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

$$c_p(L_1) = \bar{c}_p(-L_2), \quad c_n(L_1) = \bar{c}_n(-L_2). \quad (37)$$

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$:

$$c_S(L_1) = \sqrt{a} \bar{c}_S(-L_2). \quad (38)$$

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:

$$ \sigma = \frac{2(L_1 + L_2)}{V} \left[ \int_{-L_1}^{L_1} (j_p - j_n + j_s) \, dx + a \int_{-L_2}^{L_2} (\bar{j}_p - \bar{j}_n + \bar{j}_s) \, dx \right]. \quad (39) $$

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

$$ \sigma(\omega) = \sigma'(\omega) + i\sigma''(\omega) \quad (40) $$

with $\sigma'$ being the real and $\sigma''$ 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

$$ \nabla^2 \delta U(\mathbf{r}, \omega) = 0. \quad (41) $$

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_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.

<table><caption>**Table 4.** Model parameters used for the calculation of the IP response, if not stated otherwise.</caption><thead><tr><th>Parameter</th><th>Symbol [unit]</th><th>Value</th></tr></thead><tbody><tr><td>Equilibrium NaCl concentration</td><td>c<sub>∞</sub> [mol m<sup>-3</sup>]</td><td>1</td></tr><tr><td>Absolute temperature</td><td>T [K]</td><td>293</td></tr><tr><td>Relative electrolyte permittivity</td><td>ε<sub>r</sub> [-]</td><td>80</td></tr><tr><td>Equilibrium ion mobility</td><td>μ<sub>0</sub> [m<sup>2</sup> V<sup>-1</sup> s<sup>-1</sup>]</td><td>5 × 10<sup>-8</sup></td></tr><tr><td>Equilibrium Stern-layer cation mobility</td><td>μ<sub>S</sub> [m<sup>2</sup> V<sup>-1</sup> s<sup>-1</sup>]</td><td>5 × 10<sup>-9</sup></td></tr><tr><td>Pore length of wide pore</td><td>2L<sub>1</sub> [m]</td><td>9 × 10<sup>-5</sup></td></tr><tr><td>Pore length of narrow pore</td><td>2L<sub>2</sub> [m]</td><td>1 × 10<sup>-5</sup></td></tr><tr><td>Pore radius of wide pore</td><td>r<sub>1</sub> [m]</td><td>2 × 10<sup>-6</sup></td></tr><tr><td>Pore radius of narrow pore</td><td>r<sub>2</sub> [m]</td><td>2 × 10<sup>-7</sup></td></tr></tbody></table>

## 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

$$ \sigma_0 = 2F\mu_0c_\infty \quad (42) $$

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'(\omega \to 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
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Complex-conductivity spectra. (a) Real part normalized to the bulk conductivity $\sigma_0$, (b) real part normalized to the respective low-frequency limit $\sigma_{DC}$ and (c) imaginary part of the complex conductivity normalized to the bulk conductivity, produced by the new semi-analytical model, the numerical model by Bücker *et al.* (2019) and the old analytical model by Bücker & Hördt (2013a) for $\Sigma_S = 1$ mC m$^{-2}$ and $\varphi_d = -25$ mV. The yellow lines indicate the maximum imaginary conductivity $\sigma_{max}$ and the relaxation time $\tau$ of the new semi-analytical model. The grey area indicates the frequency range, in which our new semi-analytical model is not interpretable.</figcaption></figure>

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 ~ $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''_{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 =$

$(\omega(\sigma''_{max}))^{-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
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Imaginary part of the normalized complex conductivity of the three models for (a) $\Sigma_S = 57.5$ mC m$^{-2}$ and $\varphi_d = -46$ mV, (b) $\Sigma_S = 57.5$ mC m$^{-2}$ and $\varphi_d = -165$ mV, (c) $\Sigma_S = 2.5$ mC m$^{-2}$ and $\varphi_d = -46$ mV and (d) $\Sigma_S = 2.5$ mC m$^{-2}$ and $\varphi_d = -165$ mV.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Maximum imaginary conductivity (a, c and e) and relaxation time (b, d and f) for the numerical model (a and b), the new semi-analytical model (c and d) and the old analytical model (e and f). The red star indicates the position of the example spectrum in Fig. 5 in the parameter space, the white dots correspond to the spectra in Fig. 6.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** (a and c) Relative Difference of the maximum imaginary conductivity and (b and d) logarithmic difference of the relaxation time for (a and b) the new semi-analytical model and (c and d) the old analytical model.</figcaption></figure>

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

$$ \Delta_{\text{rel}}\sigma''_{\text{max}} = \frac{\sigma''_{\text{max}} - \sigma''_{\text{max, num}}}{\sigma''_{\text{max, num}}} \qquad (43) $$

and the logarithmic difference of the relaxation time

$$ \Delta_{\log} \tau = \log(\tau) - \log(\tau_{\text{num}}) \qquad (44) $$

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
<figure><img src="image_10.png" /><figcaption>**Figure 9.** (a and b) Maximum imaginary conductivity and (c and d) relaxation time of the new semi-analytic model for (a and c) the basal surface of montmorillonite at varying NaCl concentration and absolute surface-charge density and (b and d) the quartz (0001) surface at varying NaCl concentration and pH. The black dashed lines display isolines of the d-plane potential (in V), while the white dashed lines display isolines of the logarithm of the Stern-layer surface-charge density (in C m⁻²).</figcaption></figure>

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*:

$$ \tau_L = \frac{R^2}{2D_S M} \propto \frac{1}{M}, \qquad (45) $$

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_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⁻³, 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⁻³, 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.
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Imaginary conductivity spectra for different NaCl concentrations for a pore system with a quartz surface at pH 5.5. Black circles indicate the picked peak position for each spectrum.</figcaption></figure>

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 \text{ 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 \text{ mol m}^{-3}$, the spectrum is dominated by a low-frequency peak at $\sim 0.3 \text{ 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 \text{ Hz}$ builds up and becomes dominant below $10 \text{ 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
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Maximum imaginary conductivity (a and c) and relaxation time (b and d) for the the new (a and b) and the old analytic model (c and d). The length of the wide pore is fixed to $2L_1 = 2 \times 10^{-5}$ m and that of the narrow pore to $2L_2 = 2 \times 10^{-4}$ m. The EDL parameters are $\Sigma_S = 7.8$ mC m$^{-2}$ and $\phi_d = -75$ mV. The grey area indicates the parameter space with $r_2 > r_1$ that has not been considered. The white area (a and b) indicates the parameter combination where no distinct low-frequency peak was observable.</figcaption></figure>

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 µm.

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Variation of the maximum imaginary conductivity with NaCl concentration for a pore system with a montmorillonite surface with different surface-charge densities. The wide pore has a length of 5 µm and a radius of 5 Debye lengths, the narrow pore has a length of 500 nm and a radius of 2 Debye lengths. The dots indicate the position of the data points.</figcaption></figure>

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. 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. 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 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 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 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⁻² and NaCl concentrations between 1 and 100 mol m⁻³, 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 ~10 mol m⁻³) and a dominant polarization of the diffuse layer (for NaCl concentrations above ~10 mol m⁻³), 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 × 10⁻⁸ nm for the wide and below 5 × 10⁻⁹ 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.

## SUPPORTING INFORMATION

Supplementary data are available at *GJIRAS* online.
*stern layer pore throats suppl mat revision1.pdf*

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

## DATA AVAILABILITY

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

## REFERENCES

* Atekwana, E.A. & Slater, L., 2009. Biogeophysics: A new frontier in earth science research, *Rev. Geophys.*, **47**, RG4004. doi: 10.1029/2009RG000285.

* Bairlein, K., Bücker, M., Hördt, A. & Hinze, B., 2016. Temperature dependence of spectral induced polarization data: experimental results and membrane polarization theory, *Geophys. J. Int.*, **205**, 440–453.

* Blaschek, R. & Hördt, A., 2009. Numerical modelling of the IP effect at the pore scale, *Near Surf. Geophys.*, **7**, 579–588.

* Bleil, D.F., 1953. Induced polarization: A method of geophysical prospecting, *Geophysics*, **18**, 636–661.

* Buchheim, W. & Irmer, G., 1979. Zur Theorie der induzierten galvanischen Polarisation in Festkörpern mit elektrolytischer Porenfüllung, *Gerl. Beitr. Geophys.*, **88**, 53–72.

* Bücker, M. & Hördt, A., 2013a. Analytic modelling of membrane polarization with explicit parametrization of pore radii and the electrical double layer, *Geophys. J. Int.*, **194**, 804–813.

* Bücker, M. & Hördt, A., 2013b. Long and short narrow pore models for membrane polarization, *Geophysics*, **78**, E299–E314. doi: 10.1190/geo2012-0548.1.

* Bücker, M., Flores-Orozco, A., Undorf, S. & Kemna, A., 2019. On the role of Stern- and diffuse-layer polarization mechanisms in porous media, *J. Geophys. Res. Solid Earth*, **124**, 5656–5677.

* Börner, F. & Schön, J.H., 1991. A relation between the quadrature component of electrical conductivity and the specific surface area of sedimentary rocks, *Log anal.*, **32**, 612–613.

* Börner, F., Schöpper, J. & Weller, A., 1996. Evaluation of transport and storage properties in the soil and groundwater zone from induced polarization measurements, *Geophys. Prospect.*, **44**, 583–601.

* Bourg, I.C., Sposito, G. & Bourg, A.C. M., 2007. Modeling the acid-base surface chemistry of montmorillonite, *J. Colloid Interface Sci.*, **312**, 297–310. doi: 10.1016/j.jcis.2007.03.062.

* Chuprinko, D. & Titov, K., 2017. Influence of mineral composition on spectral induced polarization in sediments, *Geophys. J. Int.*, **209**, 186–191.

* Cole, K.S. & Cole, R.H., 1941. Dispersion and absorption in dielectrics. I. Alternating current characteristics, *J. Chem. Phys.*, **9**, 341–351.

* Dias, C.A., 2000. Developments in a model to describe low-frequency electrical polarization of rocks, *Geophysics*, **65**, 437–451.

* Davies, C.W., 1962. *Ion Association*, pp. 37–53, Butterworths, London

* Delgado, A.V., Arroyo, F.J., González-Caballero, F., Shilov, V.N. & Borkovskaya, Y.B., 1998. The effect of the concentration on the mechanisms of low-frequency dielectric dispersion (LFDD) in colloidal systems, *Colloids Surf. A: Physicochem. Eng. Asp.*, **140**, 139–149.

* Dukhin, S.S. & Shilov, V.N., 1974. *Dielectric Phenomena and the Double Layer in Disperse Systems and Polyelectrolytes*, Wiley, New York.

* Fixman, M., 1980. Charged macromolecules in external fields. I. The sphere, *J. Chem. Phys.*, **72**, 5177–5186.

* Gomaa, M.M. & El-Diwany, E.A., 2020. A new generalized membrane polarization frequency-domain impedance formula, *J. Appl. Geophys.*, **177**, 104023. doi: 10.1016/j.jappgeo.2020.104023.

* Hiemstra, T. & Van Riemsdijk, V.H., 2006. On the relationship between charge distribution, surface hydration, and the structure of the interface of metal hydroxides, *J. Colloid Interface Sci.*, **301**, 1–18.

* Hiemstra, T., De Wit, J.C. M. & Van Riemsdijk, W.H., 1989. Multiside proton adsorption modeling at the solid/solution interface of (hydr)oxides: A new approach: II. Application to various important (hydr)oxides, *J. Colloid Interface Sci.*, **133**, 105–117.

* Hördt, A., Druiventak, A., Blaschek, R., Binot, F., Kemna, A., Kreye, P. & Zisser, N., 2007. Case histories of hydraulic conductivity estimation with induced polarization at the field scale, *Near Surf. Geophys.*, **7**, 529–545. doi: 10.3997/1873-0604.2009035.

* Hördt, A., Bairlein, K., Bielefeld, A., Bücker, M., Kuhn, E., Nordsiek, S. & Stebner, H., 2016. The dependence of induced polarization on fluid salinity and pH, studied with an extended model of membrane polarization, *J. Appl. Geophys.*, **135**, 408–417.

* Hördt, A., Bairlein, K., Bücker, M. & Stebner, H., 2017. Geometrical constraints for membrane polarization, *Near Surf. Geophys.*, **15**, 579–592.

* Kirchhoff, G., 1845. Ueber den Durchgang eines elektrischen Stromes durch eine Ebene, insbesondere durch eine kreisförmige, *Ann. Phys.*, **140**, 497–514.

* Koch, K., Kemna, A., Irving, J. & Holliger, K., 2011. Impact of changes in grain size and pore space on the hydraulic conductivity and spectral induced polarization response of sand, *Hydrol. Earth Syst. Sci.*, **15**, 1785–1794.

* Kruschwitz, S., Binley, A., Lesmes, D. & Elshenawy, A., 2010. Textural controls on low-frequency electrical spectra of porous media, *Geophysics*, **75**, WA113–WA123. doi: 10.1190/1.3479835.

* Leroy, P. & Revil, A., 2004. A triple-layer model of the surface electrochemical properties of clay minerals, *J. Colloid Interface Sci.*, **270**, 371–380.

* Leroy, P., Revil, A., Kemna, Cosenza, P. & Ghorbani, A., 2008. Complex conductivity of water-saturated packs of glass beads, *J. Colloid Interface Sci.*, **321**, 103–117.

* Leroy, P., Devau, N., Revil, A. & Bizi, M., 2013. Influence of surface conductivity on the apparent zeta potential of amorphous silica nanoparticles, *J. Colloid Interface Sci.*, **410**, 81–93.

* Leroy, P., Tournassat, C., Bernard, O., Devau, N. & Azaroual, M., 2015. The electrophoretic mobility of montmorillonite. Zeta potential and surface conductivity effects, *J. Colloid Interface Sci.*, **451**, 21–39.

* Leroy, P., Weigand, M., Mériguet, G., Zimmermann, E., Tournassat, C., Fagerlund, F., Kemna, A. & Huisman, J.A., 2017. Spectral induced polarization of Na-montmorillonite dispersions, *J. Colloid Interface Sci.*, **505**, 1093–1110.

* Leroy, P., Maineult, A., Li, S. & Vinogradov, J., 2022. The zeta potential of quartz. Surface complexation modelling to elucidate high salinity measurements, *Colloids Surf. A: Physicochem. Eng. Asp.*, **650**, 129507.

* Lesmes, D.P. & Morgan, F.D., 2001. Dielectric spectroscopy of sedimentary rocks, *J. Geophys. Res.*, **106**, 13261–13816. doi: 10.1029/2000JB900298.
de Lima, O.A. L. & Sharma, M.M., 1992. A generalized Maxwell-Wagner theory for membrane polarization in shaly sands, *Geophysics*, **57**, 431–440. doi: 10.1190/1.1443257.

Lyklema, J., Dukhin, S. & Shilov, V., 1983. The relaxation of the double layer around colloidal particles and the low-frequency dielectric dispersion. Part I. Theoretical considerations, *J. Electroanal. Chem. Interfacial Electrochem.*, **143**, 1–21.

Lyklema, J., Springer, M.M., Shilov, V.N. & Dukhin, S.S., 1986. The relaxation of the double layer around colloidal particles and the low-frequency dielectric dispersion. Part III. Application of theory to experiments, *J. Electroanal. Chem.*, **198**, 19–26.

Maierhofer, T., Hauck, C., Hilbich, C., Kemna, A. & Flores-Orozco, A., 2022. Spectral induced polarization imaging to investigate an ice-rich mountain permafrost site in Switzerland, *Cryosphere*, **16**, 1903–1925.

Maineult, A., Revil, A., Camerlynck, C., Florsch, N. & Titov, K., 2017. Upscaling of spectral induced polarization response using random tube networks, *Geophys. J. Int.*, **209**, 948–960.

Marshall, D.J. & Madden, T.R., 1959. Induced polarization, a study of its causes, *Geophysics*, **24**, 790–816.

Mendieta, A., Jougnot, D., Leroy, P. & Maineult, A., 2021. Spectral induced polarization characterization of non-consolidated clays for varying salinities–an experimental study, *J. Geophys. Res. Solid Earth*, **126**, e2020JB021125. doi: 10.1029/2020JB021125.

Mudler, J., Hördt, A., Kreith, D., Sugand, M., Bazhin, K., Lebedeva, L. & Radic, T., 2022. Broadband spectral induced polarization for the detection of permafrost and an approach to ice content estimation–a case study from Yakutia, Russia, *Cryosphere*, **16**, 4727–4744.

Niu, Q., Zhang, C. & Prasad, M., 2020. A framework for pore-scale simulation of effective electrical conductivity and permittivity of porous media in the frequency range from 1 mHz to 1 GHz, *J. Geophys. Res. Solid Earth*, **125**, e2020JB020515. doi: 10.1029/2020JB020515.

O'Konski, C.T., 1960. Electrical properties of macromolecules. V. Theory of ionic polarization in polyelectrolytes, *J. Phys. Chem.*, **64**, 605–619.

Parkurst, D.L. & Apelo, C.A., 2013. Description of input and examples for PhreeqC version 3–A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations, *US Geol. Surv. Tech. Methods*, **2013**(6), 497. URL: http://pubs.usgs.gov/tm/06/a43

Revil, A., 2012. Spectral induced polarization of shaly sands: Influence of the electrical double layer, *Water Resour. Res.*, **48**, W02517. doi: 10.1029/2011WR011260.

Revil, A., Vaudelet, P., Su, Z. & Chen, R., 2022. Induced polarization as a tool to assess mineral deposits: A review, *Minerals*, **12**, 571. doi: 10.3390/min12050571.

Schlumberger, C., 1920. *Étude sur la prospection électrique du sous-sol*, Gouthier-Villars, Paris.

Schurr, J.M., 1964. On the theory of the dielectric dispersion of spherical colloidal particles in electrolyte solution, *J. Phys. Chem.*, **68**, 2407–2413.

Schwarz, G., 1962. A theory of the low-frequency dielectric dispersion of colloid particles in electrolyte solution, *J. Phys. Chem.*, **66**, 2636–2642.

Seigel, H., Nabighian, M., Parasnis, D.S. & Vozoff, K., 2007. The early history of the induced polarization method, *Leading Edge*, **26**, 312–321.

Skold, M., Revil, A. & Vaudelet, P., 2011. The pH dependence of spectral induced polarization of silica sands: Experiment and modeling, *Geophys. Res. Lett.*, **38**, L12304. doi: 10.1029/2011GL047748.

Slater, L., 2007. Near surface electrical characterization of hydraulic conductivity: From petrophysical properties to aquifer geometries–A review, *Surv. Geophys.*, **28**, 169–197.

Stebner, H., Halisch, M. & Hördt, A., 2017. Simulation of membrane polarization of porous media with impedance networks, *Near Surf. Geophys.*, **15**, 563–578.

Titov, K., Komarov, V., Tarasov, V. & Levitski, A., 2002. Theoretical and experimental study of time domain-induced polarization in water-saturated sands, *J. Appl. Geophys.*, **50**, 417–433.

Tournassat, C., Chapron, Y., Leroy, P., Bizi, M. & Boulahya, F., 2009. Comparison of molecular dynamics simulations with triple layer and modified Gouy-Chapman models in a 0.1 M NaCl-montmorillonite system, *J. Colloid Interface Sci.*, **339**, 533–541.

Tournassat, C., Bourg, I.C., Steefel, C.I. & Bergaya, F., 2015. Surface properties of clay minerals, *Developments in Clay Science, Vol. 6*, Elsevier B.V., pp. 5–31

Volkmann, J. & Klitzsch, N., 2010. Frequency-dependent electric properties of microscale rock models for frequencies from one millihertz to ten kilohertz, *Vadose Zone J.*, **9**, 858–870.

Wainwright, H.M., Flores-Orozco, A., Bücker, M., Dafflon, B., Chen, J., Hubbard, S.S. & Williams, K.H., 2016. Hierarchical Bayesian method for mapping biogeochemical hot spots using induced polarization imaging, *Water Resour. Res.*, **52**, 533–551.

Weller, A., Slater, L., Nordsiek, S. & Ntarlagiannis, D., 2010. On the estimation of specific surface per unit pore volume from induced polarization: A robust empirical relation fits multiple data sets, *Geophysics*, **75**, WA105–WA112.

Weller, A., Breede, K., Slater, L. & Nordsiek, S., 2011. Effect of changing water salinity on complex conductivity spectra of sandstones, *Geophysics*, **76**, F315–F327.

Weller, A., Slater, L. & Nordsiek, S., 2012. On the relationship between induced polarization and surface conductivity: Implications for petrophysical interpretation of electrical measurements, *Geophysics*, **78**, D315–D325. doi: 10.1190/geo2013-0076.1.

Westall, J. & Hohl, H., 1980. A comparison of electrostatic models for the oxide/solution interface, *Adv. Colloid Interface Sci.*, **12**, 265–294.

# APPENDIX A: CALCULATION OF SURFACE CHARGE DENSITIES

## A1 Montmorillonite

The equilibrium constant of Na-adsorption can be calculated as

$$K_{\text{Na,Mt}} = \frac{\Gamma_{\text{X}^{-\dots\text{Na}^{+}}}}{\Gamma_{\text{X}^{-} - a_{\text{Na}}}} \exp\left(\frac{e\varphi_{\text{d}}}{k_{\text{B}}T}\right) \quad (\text{A1})$$

with $\Gamma_i$ being the surface-site density of species i. In particular, $\Gamma_X^{-}$ is the surface-site density of the clay sites > X$^{-}$ without a sodium cation attached. The total surface-charge density is connected to the surface-site densities by

$$\Sigma_0 = -e (\Gamma_{X^{-}} + \Gamma_{X^{-} \dots Na^{+}}). \quad (\text{A2})$$

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

$$\Gamma_{X^{-} \dots Na^{+}} = - \frac{K_{\text{Na,Mt}} a_{\text{Na}} \exp\left(-\frac{e\varphi_{\text{d}}}{k_{\text{B}}T}\right)}{1 + K_{\text{Na,Mt}} a_{\text{Na}} \exp\left(-\frac{e\varphi_{\text{d}}}{k_{\text{B}}T}\right)} \Sigma_0 \quad (\text{A3})$$

and with

$$\Sigma_S = e \Gamma_{X^{-} \dots Na^{+}} \quad (\text{A4})$$

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

$$K_H = \frac{\Gamma_{\text{SiOH}}}{\Gamma_{\text{SiO}^{-}} - a_H} \exp\left(\frac{e\varphi_0}{k_B T}\right) \quad (\text{A5})$$
and

$$K_{\text{Na,Qtz}} = \frac{\Gamma_{\text{SiO}^{-\dots\text{Na}^{+}}}}{\Gamma_{\text{SiO}^{-}} a_{\text{Na}}} \exp\left(\frac{e\varphi_{\text{d}}}{k_{\text{B}}T}\right). \qquad (\text{A6})$$

These equations can be rearranged to

$$\Gamma_{\text{SiOH}} = K_{\text{H}} a_{\text{H}} \exp\left(-\frac{e\varphi_{\text{0}}}{k_{\text{B}}T}\right) \Gamma_{\text{SiO}^{-}} \qquad (\text{A7})$$

and

$$\Gamma_{\text{SiO}^{-\dots\text{Na}^{+}}} = K_{\text{Na,Qtz}} a_{\text{Na}} \exp\left(-\frac{e\varphi_{\text{d}}}{k_{\text{B}}T}\right) \Gamma_{\text{SiO}^{-}}. \qquad (\text{A8})$$

The total surface-site density is

$$\Gamma_{\text{S}} = \Gamma_{\text{SiO}^{-}} + \Gamma_{\text{SiOH}} + \Gamma_{\text{SiO}^{-\dots\text{Na}^{+}}} = A\Gamma_{\text{SiO}^{-}} \qquad (\text{A9})$$

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

$$\Gamma_{\text{SiO}^{-}} = \frac{\Gamma_{\text{S}}}{A} \qquad (\text{A10})$$

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

$$\Sigma_{0} = -e(\Gamma_{\text{SiO}^{-}} + \Gamma_{\text{SiO}^{-\dots\text{Na}^{+}}}) \qquad (\text{A11})$$

and at the Stern plane as

$$\Sigma_{\text{S}} = -e\Gamma_{\text{SiO}^{-\dots\text{Na}^{+}}} \qquad (\text{A12})$$

yielding eqs (14) and (15).