Freely available online through the author-supported open access option.

Abstract

In this study, spectral induced polarization (SIP) spectra were generated numerically to better understand how actual rock microstructure and electrolyte properties in rock pores affect the spectral pattern, i.e., the characteristic relaxation time of polarization as well as the polarization strength of a rock pore system. The dynamics of charge carriers in three-dimensional pore systems were simulated using a frequency-dependent formulation of the Nernst–Planck–Poisson (NPP) ion-transport equations. Basically, a pore-system model of alternating stacked cylinders of two different sizes was studied considering the electrical double layer (EDL). A reduced cationic mobility—resulting from increased adsorption of these cations at the rock–water interface—was assumed within the EDL. By solving the NPP equations using the finite element method, complex resistivity phase spectra were generated. Subsequently, the effect of pore structural properties and electrolyte conductivity on the magnitude and frequency position of the characteristic resistivity phase minimum of rocks was studied. The following results were found: First, regarding pore geometry, the characteristic frequency of the phase minimum fmin decreases with increasing pore length of the large pore. Second, both small pores having a radius of a few Debye lengths combined with larger pores are needed to ensure detectable phase amplitudes. Third, with regard to electrolyte concentration, the phase amplitude is inversely proportional to the concentration, whereas fmin remains constant. Because the studied model does not provide a direct and exclusive link between the simulated electrical properties and pore throat size, further research is needed here to specify a convincing SIP interpretation method for improved permeability estimation.

The influence of pore structure and water salinity on the frequency-dependent electrical properties (i.e., the spectral induced polarization response) of simple water-saturated pore systems was studied by simulating the underlying ion transport processes, diffusion and migration, at the microscale.

In geophysics, frequency-dependent complex resistivity measurements are a typical noninvasive method for field- and laboratory-scale applications. This method is called spectral induced polarization (SIP) or the complex resistivity method, first discovered by Schlumberger (1920). In other fields, it is known as impedance spectroscopy. Estimating permeability by such noninvasive electrical measurements has been extensively studied since the 1950s (e.g., Carman, 1956, Ch. 1) and is still a current field of research (e.g., Binley et al., 2005; Titov et al., 2010). To achieve this estimation, one can take advantage of the dependence of both hydraulic and complex or frequency-dependent electric conductivity on microstructural rock properties. A simple model with respect to permeability k is given by the Kozeny–Carman (Carman, 1956, Ch. 1) and related equations (e.g., Katz and Thompson, 1986). These state a dependence of k on porosity ϕ and the surface area to volume ratio Spor, which in turn can be linked to the electrical formation factor ℱ and a characteristic hydraulic length scale l of the pore space. There are several approaches to obtaining these quantities from electrical measurements (summarized by Slater, 2007). The value of ℱ (or ϕ) can be obtained using Archie's law (Archie, 1942), whereas Spor is related to the imaginary part of the electrical conductivity σ″ for a constant phase angle model (Börner et al., 1996). On the other hand, clear relationships between Spor and the SIP relaxation time τ as well as between the dominant pore throat size D0 and τ have been reported for samples—mostly from a single rock formation—that do not show a constant phase angle (e.g., Binley et al., 2005; Scott and Barker, 2005; Titov et al., 2010). Although such a relation does not necessarily apply to a compilation of rock samples from different formations—as shown for a data set comprising measurements on different lithologies (Kruschwitz et al., 2010)—some researchers have proposed a direct τ–k relationship. This approach has been confirmed by simple theoretical models predicting relationships between τ and a microstructural length property  
\[{\tilde{l}}\]
, which is possibly connected to the parameters l and Spor (Marshall and Madden, 1959; Schwarz, 1962; Titov et al., 2002).

Nevertheless, most of the relationships between the deduced hydraulic rock parameters and the electric quantities are empirical due to the lack of universal physical rock models. Consequently, measured spectra are often described by equivalent circuits that reproduce the spectral behavior by appropriately combining resistances and capacitances. One basic model of this type is the Cole–Cole model (Cole and Cole, 1941), for which there are many possible modifications. Dias (2000) summarized the typical equivalent circuits and their spectral behavior, which can be applied to the required experimental setup.

The empirical equivalent circuits are widely used for data interpretation, but as mentioned above, a few approaches exist that link the SIP quantities to grain sizes (Schwarz, 1962; Fixman, 1980) or to pore lengths and radii (Marshall and Madden, 1959; Titov et al., 2002) by means of microscopic models for simple systems of mineral grains or connected pores. This study focused on how microstructural and electrochemical parameters affect the SIP properties by means of numerical simulations taking into account earlier theories for simple model systems. Thus, the aim of this investigation was to advance understanding of which rock properties are the most relevant for the magnitude and shape of SIP phase spectra.

Theories

The frequency dependence of the electric impedance of rocks containing water (first described by Schlumberger, 1920) results from the buildup of an EDL at the solid–liquid interface caused by surface charges at the inner surface area. Thus, frequency-dependent effects of the diffusion and migration of charge carriers arise in the pore space. Recent theoretical studies have concentrated on approaches to characterize the electric behavior of rocks as a macroscopic multicomponent system for a wider frequency range, thus combining different influencing polarization mechanisms (e.g., de Lima and Sharma, 1992; Leroy et al., 2008; Leroy and Revil, 2009). In contrast, this study focused on the underlying microscopic effects that cause polarization at the scale of single grains or pores.

In this regard, Schwarz (1962) studied the dielectric dispersion of spherical particles in an electrolyte solution and found a Debye-type dispersion with a relaxation time τp that is dependent on the particle radius R: 
\[\mathrm{{\tau}}_{\mathrm{p}}{=}\frac{R^{2}}{2uk_{\mathrm{b}}T}\]
[1]
where u is the mechanical surface mobility (in s kg−1) of the counterions, kb is Boltzmann's constant, and T is the temperature.
Another theory to explain the SIP effects in water-saturated rocks is the so-called “membrane polarization” theory. With regard to this theory, Marshall and Madden (1959) studied a simple one-dimensional pore system of alternating “active” zones and “passive” zones. The “active” zones, representing pore throat passages, are characterized by a reduced anion transport number. Therefore, because these active zones basically block anions, they behave as ion-selective membranes. For small excess electrolyte concentrations and electric fields, Marshall and Madden (1959) solved the dependence of the complex impedance Z on the angular frequency ω = 2πf, with f being the frequency, as follows: 
\begin{eqnarray*}&&Z{=}\\&&\frac{L_{1}}{\mathrm{{\mu}}_{\mathrm{p},1}CF}\left[t_{\mathrm{p},1}\right.{+}\frac{B}{A}t_{\mathrm{p},2}\\&&{+}\frac{(S_{2}{-}S_{1})^{2}}{\frac{X_{1}S_{1}}{(t_{\mathrm{p},2})^{2}t_{\mathrm{p},1}\mathrm{tanh}(X_{1})}{+}\frac{A}{B}\frac{X_{2}S_{2}}{t_{\mathrm{p},2}(t_{\mathrm{p},1})^{2}\mathrm{tanh}(X_{2})}}\right]\end{eqnarray*}
[2]
with 
\[X_{i}{=}\left(\frac{i\mathrm{{\omega}}}{2D_{\mathrm{p},i}t_{\mathrm{n},i}}\right)^{1/2}\frac{L_{i}}{2}\]
[3]
and 
\[A{=}\frac{L_{1}}{L_{2}},{\,}{\,}B{=}\frac{D_{\mathrm{p},1}}{D_{\mathrm{p},2}},{\,}S_{i}{=}\frac{t_{\mathrm{n},i}}{t_{\mathrm{p},i}}\]
[4]
where the subscript 1 denotes the quantities of the passive zone and subscript 2 those of the active zone. The subscript p stands for the cations and n for the anions, while Li indicates the lengths of the pores, μp,i and μni the ion mobilities, Dp,i and Dn,i their diffusion coefficients, tp,i and tn,i the ion-transport numbers [tp,i = μp,i/(μp,i + μn,i), tn,i = μn,i/(μp,i + μn,i)], C is the equilibrium concentration of both types of ions, and F is the Faraday constant. Because the Marshall and Madden model is a one-dimensional model, only the lengths of the pores contribute to the resistivity spectra. Similarly, Titov et al. (2002) developed the short narrow pore (SNP) model with the length of the active zones much smaller than that of the passive zones. They gave the following expression for the complex resistivity ρ(ω): 
\[\mathrm{{\rho}}(\mathrm{{\omega}}){=}\mathrm{{\rho}}_{0}\left(1{-}\mathrm{{\eta}}_{0}\left\{1{-}\frac{1{-}\mathrm{exp}\left[{-}2(i\mathrm{{\omega}{\tau}}_{\mathrm{snp}})^{1/2}\right]}{2(i\mathrm{{\omega}{\tau}}_{\mathrm{snp}})^{1/2}}\right\}\right)\]
[5]
with chargeability η0 being 
\[\mathrm{{\eta}}_{0}{=}\frac{4L_{\mathrm{small}}({\Delta}n)^{2}}{\left(\frac{L_{\mathrm{large}}}{A_{\mathrm{large}}}{+}\frac{L_{\mathrm{small}}}{\mathrm{{\alpha}}A_{\mathrm{small}}}\right)(A_{\mathrm{large}}{+}A_{\mathrm{small}})}\]
[6]
and the time constant τsnp being 
\[\mathrm{{\tau}}_{\mathrm{snp}}{=}\frac{L_{\mathrm{small}}^{2}}{4D}\]
[7]
where A is the pore section and α is a coefficient of efficiency of the active zones. As for the Schwarz model (Eq. [1]), a quadratic dependence on geometric properties is predicted by Eq. [7]. This is also approximately true for the one-dimensional Marshall–Madden model with respect to the active and passive zone lengths in a certain range, depending on the choice of model parameters (Eq. [2]). Thus, all the three aforementioned theories suggest a quadratic dependence of relaxation time on a geometric length scale of the porous medium, i.e., pore throat or grain sizes. This prediction can be qualitatively confirmed by some experimental data—cf., e.g., Kruschwitz (2007), where an exponent of 2.1 was given for the dependence of the Cole–Cole relaxation time τCC on the dominant pore throat diameter D0 (see Kruschwitz, 2007, Fig. 3.26). Nevertheless, from an experimental point of view, there is considerable uncertainty with respect to the influence of pore or grain geometry on the relaxation time. Binley et al. (2005) found a nearly linear behavior τCCD01.04. The recent work of Kruschwitz et al. (2010) has implied a relationship τGCCD03.0 with the relaxation time τGCC of a generalized Cole–Cole model. According to these researchers, this relation is valid for samples with large pore-throat diameters while there is no influence of the pore-throat diameter on the τGCC for small pore throats (<7 μm). Other experiments have suggested an exponential relation τpeak ∼ exp(D0/5.9) and τpeak ∼ exp(D0/5.17) (Scott and Barker, 2003, 2005), where τpeak is the peak relaxation time, which is essentially the inverse of the peak position frequency fpeak. This has been supported by the work of Titov et al. (2010) with a relation  
\[{\bar{{\tau}}}\]
∼ exp(D0/16.8), where  
\[{\bar{{\tau}}}\]
is a characteristic relaxation time, which is the mean value of modal relaxation times for different conductivities. The proposed τ–D0 relations are summarized in Table 1 . The definitions of the different relaxation times depend on the model that is used to fit SIP spectra but is always connected to the inverse of the characteristic phase-peak frequency. Hence, this phase-peak frequency was used in this study.

In addition to the influence of pore geometry, an influence of electrolyte conductivity on the SIP spectra has been observed. Nevertheless, the classic model of Marshall and Madden (1959) and the SNP model of Titov et al. (2002) do not consider the effect of the concentration of charge carriers on the phase spectra. Some experiments, on the other hand, have suggested a general trend of decreasing phase-peak values (e.g., Lesmes and Frye, 2001) or chargeability (e.g., Titov et al., 2010) as a measure of polarization strength with increasing fluid conductivity. To eliminate this major influence of fluid conductivity on the phase-peak values, a normalized chargeability—namely, the product of chargeability and rock conductivity (e.g., Lesmes and Frye, 2001) or water conductivity (e.g., Titov et al., 2010)—has been introduced, which then showed minor variations with fluid conductivity. In addition, the data presented by Lesmes and Frye (2001) showed a maximum of the normalized chargeability at approximately 1 S m−1. The data of Scott and Barker (2005) showed an increasing normalized chargeability for conductivities <1 S m−1. In the same study, the data of Flath (1989) were presented, depicting a maximum of quadrature conductivity between 1 and 10 S m−1. In contrast, Titov et al. (2010) presented data with an increasing normalized chargeability even for conductivities >1 S m−1, thus assuming the maximum to be shifted to higher conductivities for their data set. Furthermore, variations of the relaxation time with variations in water conductivity have often been reported (Scott and Barker, 2005; Kruschwitz, 2007; Titov et al., 2010); however, these variations have typically been rather slight.

In conclusion:

  • a wide range of τ–D0 relations has been proposed in the literature (Table 1), e.g., relations τ ∼ D0a with exponents a ranging from 1 to 3 or relations τ ∼ exp(D0/b) with b ranging from 5.17 to 16.8;

  • increasing water conductivity → decreasing phase peak values;

  • relaxation time varies slightly with variations of water conductivity.

To improve the theoretical understanding of these effects, numerical approaches have been taken, e.g., by Blaschek and Hördt (2009) for one- and two-dimensional membrane polarization models. They applied the general concept of Marshall and Madden (1959), thus assuming a membrane effect in small pore throats, described, e.g., by a reduced anion mobility within these pores.

Unfortunately, the existing analytical and numerical polarization models cannot give any information about the size necessary to make a pore throat an active zone, thus it is difficult to obtain valuable results with respect to pore geometry. Because the dimensions of the EDL are not considered, the influence of fluid conductivity—defining the characteristic distance from the pore walls—is not described adequately, too. Consequently, in this study, a three-dimensional pore system model was developed that considers the double layer itself to examine the pore size with respect to the membrane polarization concept. Furthermore, this work investigated the influence of electrolyte properties (coupled to the thickness of the double layer) on the phase spectra. An effective numerical procedure was ensured by using a two-dimensional axial symmetric representation of the three-dimensional pore model and a frequency-dependent formulation of the governing equations.

The equations were solved using the finite elements software COMSOL Multiphysics (COMSOL AB, Stockholm, Sweden) to simulate the frequency-dependent resistivities with the general aim of achieving insight into the structural and electrochemical properties influencing SIP measurements, thus improving the prediction of permeability from SIP measurements.

Modeling Approach and Governing Equations

The microscale ion transport in rocks saturated with electrolyte solution was simulated and the frequency-dependent response of the model systems to a sinusoidal voltage was investigated. Therefore, three different models were used:

  1. The classic Marshall–Madden model enhanced to three dimensions (Model 1, Fig. 1a

    Fig. 1.

    Two-dimensional axial symmetric representation of three-dimensional cylindrical pore models with applied sinusoidal voltage U = U0sin(ωt) (with voltage amplitude U0, angular frequency ω, and time t): (a) Model 1, three-dimensional Marshall–Madden model with active (anion mobility μn < cation mobility μp) and passive (μn = μp) zones; (b) Model 2, mineral grain model with electrical double layer (μp < μn) and free electrolyte (μn = μp); (c) Model 3, simple pore system with electrical double layer (μp < μn) and free electrolyte (μn = μp).

    Fig. 1.

    Two-dimensional axial symmetric representation of three-dimensional cylindrical pore models with applied sinusoidal voltage U = U0sin(ωt) (with voltage amplitude U0, angular frequency ω, and time t): (a) Model 1, three-dimensional Marshall–Madden model with active (anion mobility μn < cation mobility μp) and passive (μn = μp) zones; (b) Model 2, mineral grain model with electrical double layer (μp < μn) and free electrolyte (μn = μp); (c) Model 3, simple pore system with electrical double layer (μp < μn) and free electrolyte (μn = μp).

    )

  2. A model of mineral grains suspended in electrolyte solution (Model 2, Fig. 1b)

  3. A pore-system model of alternating stacked cylinders of two different sizes taking into consideration the EDL (Model 3, Fig. 1c)

The surface charges at the pore walls affect the mobility of the charge carriers (gray areas in Fig. 1). We used an axial symmetric representation of the three-dimensional model system to reduce the computational effort.

The Model Systems

Enhanced Marshall–Madden Model (Model 1)

The original model of Marshall and Madden (1959) is a one-dimensional sequence of active and passive zones with reduced anion-transport numbers in the active zone representing narrow pore throats. The one-dimensional model of Marshall and Madden (1959) has been enhanced to three dimensions using an axial symmetric approach (Model 1, Fig. 1a). Thus, the pore throat is now characterized by reduced anion mobility and pore size. This enhanced model now has a pore length Lpassive = 1 μm and radius Rpassive = 1 μm for the passive zone and a length Lactive = 1 μm and radius Ractive = 0.1 μm for the active zone (Table 2 ). This Model 1 was mainly used to verify the modeling approach.

Mineral Grain Model (Model 2)

Similar to the concept of Schwarz (1962), a model of a mineral grain suspended in electrolyte solution was investigated (Model 2, Fig. 1b). This model is important because there is an ongoing controversy about whether grain size or pore size significantly affects polarization. To depict the adsorption of cations at the pore walls, a mobility reduced by one order of magnitude in the EDL was assumed. This reduced ion mobility represents the mean of the reduced mobility in the Stern layer (e.g., Zukoski and Saville, 1986; Revil and Glover, 1998) and the undiminished mobility in the diffuse layer. The thickness of the double layer is given by the Debye length (Debye and Hückel, 1923): 
\[\mathrm{{\lambda}}_{\mathrm{D}}{=}\sqrt{\frac{\mathrm{{\epsilon}}k_{\mathrm{b}}T}{2e^{2}X}}\]
[8]
with the ionic strength being 
\[X{=}\frac{1}{2}{{\sum}_{i}}z_{i}^{2}C_{i}\]
[9]
where e is the elementary charge, ε is the pore fluid permittivity, zi is the respective charge number and Ci is the ion concentration in the free electrolyte.

The mineral grain model incorporates spheres with a radius Rsphere = 1.4 × 10−7 m. The whole system has a length of Ltot = 3 × 10−7 m and a radius of Rtot = 1.5 × 10−7 m and thus a porosity of ϕ = 0.5 (Table 2). This Model 2, compared with the model of Schwarz (1962), was used for further model verification.

Simple Pore System with Electric Double Layer (Model 3)

As illustrated in Fig. 1c, the membrane polarization concept is studied in a more fundamental way by modeling a pore system taking into consideration the EDL. This model consists of alternating stacked cylinders of two different sizes. In this pore space model, a membrane effect of small pore throats on ion transport (e.g., active zone in the Marshall Madden model) is not postulated for a certain range of pore sizes. This is advantageous, because the necessary pore size to make a pore throat an active zone is not known, whereas the dimension of the EDL is given by the Debye length (Eq. [8]). Again, a reduced cationic mobility in the EDL affects the dynamics of the charge carriers. An axial symmetric representation of the pore system, as shown in Fig. 1c, with larger pores (length Llarge = 5 × 10−6 m and radius Rlarge = 5 × 10−7 m) and smaller pores (length Lsmall = 5 × 10−8 m and radius Rsmall = 5 × 10−8 m) was simulated. Two different electrolyte concentrations of C = 1 mol m−3 (fluid conductivity σf ≈ 0.01 S m−1) and C = 0.1 mol m−3f ≈ 0.001 S m−1) were used (Table 2). Pore Model 3 was used to study the effects of changes in the pore geometry and electrolyte concentration on the characteristic phase shift of the electric resistivity.

In all three model concepts, a reduced ion mobility in certain areas is the decisive model parameter. The mobility of anions, and thus their transport number, is reduced in active zones for Model 1, whereas the mobility of cations is reduced in the EDL for the other two model concepts. These approaches can be explained in different ways. The reduced mobility of anions in the active zones of Model 1 represents their reduced concentration in small pore throats. According to Marshall and Madden (1959), these regions are strongly influenced by the EDL, which leads to a deficiency of anions in the pore throat. On the other hand, the reduced mobility of cations in Models 2 and 3 represents the influence of the EDL on the dynamic properties of the different ions. This reduced mobility represents the partial adsorption of ions to the pore walls, which is the mean of the mobility in the Stern layer and in the diffuse layer. The Debye length characterizes this region of the pore space, which is influenced by the charged solid surface.

Governing Equations

According to Marshall and Madden (1959), the following Nernst–Planck–Poisson system of equations can be used to describe the charge transport in an ideal one-component electrolyte: 
\begin{eqnarray*}&&\frac{{\partial}}{{\partial}{\,}t}C_{\mathrm{p}}(\mathbf{\mathrm{X}},t){=}D_{\mathrm{p}}{\Delta}C_{\mathrm{p}}(\mathbf{\mathrm{X}},t)\\&&{+}\mathrm{{\nabla}}{\cdot}\left[\mathrm{{\mu}}_{\mathrm{p}}C_{\mathrm{p}}(\mathbf{\mathrm{X}},t)\mathrm{{\nabla}}U(\mathbf{\mathrm{X}},t)\right]\end{eqnarray*}
[10]
 
\begin{eqnarray*}&&\frac{{\partial}}{{\partial}{\,}t}C_{\mathrm{n}}(\mathbf{\mathrm{X}},t){=}D_{\mathrm{n}}{\Delta}C_{\mathrm{n}}(\mathbf{\mathrm{X}},t)\\&&{-}\mathrm{{\nabla}}{\cdot}\left[\mathrm{{\mu}}_{\mathrm{n}}C_{\mathrm{n}}(\mathbf{\mathrm{X}},t)\mathrm{{\nabla}}U(\mathbf{\mathrm{X}},t)\right]\end{eqnarray*}
[11]
 
\[{\Delta}U(\mathbf{\mathrm{X}},t){=}\frac{F}{\mathrm{{\epsilon}}}\left[C_{\mathrm{n}}(\mathbf{\mathrm{X}},t){-}C_{\mathrm{p}}(\mathbf{\mathrm{X}},t)\right]\]
[12]
where Dp and Dn are the diffusion coefficients, μp and μn are the mobilities of the ions, and F is Faraday's constant. This system of equations can be solved for the concentration of cations Cp(x,t) and anions Cn(x,t) and the electric potential U(x,t) at position vector x and time t. For the calculations, a permittivity of ε = 80 was assumed. The undiminished mobility in the free electrolyte was μp/n,f = 5 × 10−8 m2 V−1 s−1 and reduced by a factor of 10 in the active zone and EDL, respectively, for one type of charge carriers. Equations [10] and [11] are continuity equations, which couple the time-dependent change in the concentrations to the divergence of their flux. This flux of ions consists of a diffusion term jD = −Dp/nCp/n(x,t), proportional to gradients of the concentration, and a migration term jM = −zp/nμp/nCp/n(x,t)∇U(x,t), proportional to the product of the concentration Cp/n(x,t) and the electrical field E(x,t) = −∇U (x,t). Generation of an electrical field due to space charges is considered in Poisson's Eq. [12]. In turn, the migration currents caused by this electric field tend to compensate for the charge buildup. Additionally, Einstein's relation couples the ion mobilities and diffusion coefficients, which reduces the number of independent constants: 
\[D_{p/n}{=}\frac{\mathrm{{\mu}}_{p/n}k_{\mathrm{b}}T}{e}\]
[13]
where kb is Boltzmann's constant and T is the temperature. Because frequency-dependent effects were investigated in this study, a time-harmonic approach was taken to convert the governing equations from the time domain to the frequency domain. Time-harmonic electric fields, e.g., an applied sinusoidal voltage, are typical for SIP. These fields can be regarded as a small perturbation of an equilibrium state with U(x,t) ≡ 0 and constant concentrations Cp(x,t) ≡ Cn(x,t) ≡ C ∈ R. It was assumed that the variations from this state are sine functions with small amplitudes compared with the equilibrium values. For the sake of simplicity, these functions are expressed as complex numbers.
Thus, in analogy to the corresponding one-dimensional approach (Marshall and Madden, 1959), the following time dependence was assumed in this study: 
\[C_{\mathrm{p}}(\mathbf{\mathrm{X}},t){=}c_{\mathrm{p0}}{+}c_{\mathrm{p}}(\mathbf{\mathrm{X}})\mathrm{exp}(i\mathrm{{\omega}}t)\]
[14]
 
\[C_{\mathrm{n}}(\mathbf{\mathrm{X}},t){=}c_{\mathrm{n0}}{+}c_{\mathrm{n}}(\mathbf{\mathrm{X}})\mathrm{exp}(i\mathrm{{\omega}}t)\]
[15]
 
\[U(\mathbf{\mathrm{X}},t){=}u(\mathbf{\mathrm{X}})\mathrm{exp}(i\mathrm{{\omega}}t)\]
[16]
with equilibrium concentrations cp0,cn0 ∈ R and complex amplitudes cp(x),cn(x),u(x) ∈ C of the excess concentrations and electric potential.

Insertion of the approach of Eq. [14], [15], and [16] into Eq. [10], [11], and [12] leads to the following significant simplifications:

  1. The space and time derivatives of the equilibrium concentrations vanish: 
    \[\frac{{\partial}}{{\partial}t}c_{\mathrm{p0}}{=}\frac{{\partial}}{{\partial}t}c_{\mathrm{n0}}{=}\mathrm{{\nabla}}c_{\mathrm{p0}}{=}\mathrm{{\nabla}}c_{\mathrm{n0}}{=}0\]
  2. Time derivatives become products: 
    \[\frac{{\partial}}{{\partial}{\,}t}\mathrm{exp}(i\mathrm{{\omega}}t){=}i\mathrm{{\omega}exp}(i\mathrm{{\omega}}t)\]
  3. For small excess concentrations cp/n(x) and electrical fields |E(x)| = |−∇u(x)|, any product of these quantities can be neglected in Eq. [10] and [11]: 
    \[c_{p/n}(\mathbf{\mathrm{X}})\mathrm{{\nabla}}u(\mathbf{\mathrm{X}}){\rightarrow}0\]
    For cp0 = cn0 = C (e.g., NaCl solution), the differences in equilibrium concentrations vanish in Eq. [12]: 
    \[c_{\mathrm{n}0}{-}c_{\mathrm{p}0}{=}0\]

With these simplifications, Eq. [10], [11], and [12] are modified to 
\[i\mathrm{{\omega}}c_{\mathrm{p}}(\mathbf{\mathrm{X}}){=}D_{\mathrm{p}}{\Delta}c_{\mathrm{p}}(\mathbf{\mathrm{X}}){+}\mathrm{{\nabla}}{\cdot}\left[\mathrm{{\mu}}_{\mathrm{p}}C\mathrm{{\nabla}}u(\mathbf{\mathrm{X}})\right]\]
[17]
 
\[i\mathrm{{\omega}}c_{\mathrm{n}}(\mathbf{\mathrm{X}}){=}D_{\mathrm{n}}{\Delta}c_{\mathrm{n}}(\mathbf{\mathrm{X}}){-}\mathrm{{\nabla}}{\cdot}\left[\mathrm{{\mu}}_{\mathrm{n}}C\mathrm{{\nabla}}u(\mathbf{\mathrm{X}})\right]\]
[18]
 
\[{\Delta}u(\mathbf{\mathrm{X}}){=}\frac{F}{\mathrm{{\epsilon}}}\left[c_{\mathrm{n}}(\mathbf{\mathrm{X}}){-}c_{\mathrm{n}}(\mathbf{\mathrm{X}})\right]\]
[19]

The boundary condition must be transformed accordingly. The alternating voltage was applied to the upper and lower ends of the system, represented in the time-harmonic approach by constant values for the complex voltage amplitude u with a phase shift of π between the upper and lower boundaries [e.g., u(z = 0) = +i and u(z = Ltot) = −i]. Periodic boundary conditions were implemented for the complex amplitudes of the excess concentrations at these boundaries [cp/n(z = 0) = cp/n(z = Ltot)]; hence, the structure repeats itself in the z axis. A symmetry boundary condition was used for the symmetry axis. For the boundaries at the rock–electrolyte interface, no normal flux of ions or electric field normal to the boundary are allowed, because of which there are no chemical reactions at this interface. Furthermore, surface charges were not considered as a boundary condition.

The system of Eq. [17], [18], and [19] was solved parametrically with respect to the angular frequency ω. From this solution, the total current Itot through the pore system was determined as a volume integral of the z component of the current density jz over the pore volume Vpore divided by the length of the system Ltot: 
\[I_{\mathrm{tot}}{=}\frac{I}{L_{\mathrm{tot}}}{{\int}_{V_{\mathrm{pore}}}}j_{z}(\mathbf{\mathrm{X}})\mathrm{d}V\]
[20]
where the current density is the sum of the current of cations and anions: 
\begin{eqnarray*}&&j_{z}(\mathbf{\mathrm{X}}){=}F\left[{-}D_{\mathrm{p}}{\Delta}c_{\mathrm{p}}(\mathbf{\mathrm{X}}){-}\mathrm{{\mu}}_{\mathrm{p}}C\mathrm{{\nabla}}u(\mathbf{\mathrm{X}})\right.\ \\&&\left.\ {+}D_{\mathrm{n}}{\Delta}c_{\mathrm{n}}(\mathbf{\mathrm{X}}){-}\mathrm{{\mu}}_{\mathrm{n}}C\mathrm{{\nabla}}u(\mathbf{\mathrm{X}})\right]_{z}\end{eqnarray*}
[21]

Using Ohm's law, the complex resistance was obtained from both the current Itot and the applied voltage U; thus, the complex resistivity was obtained when considering the dimensions of the model system.

The permittivity ε of the pore fluid can be taken into account by adding a term jε = −F[iωε∇u]z when integrating the current density. Because this effect would dominate other effects of interest, especially at higher frequencies (see Fig. 2 ), this term was neglected in our simulation.

Preliminary Discussion about Model Verification and Concept

The numerical approach was verified by comparison of a one-dimensional sequence of “pores” (a line with varying ion mobilities) with the analytical solution (Eq. [2]) of Marshall and Madden (1959). Because the results of Marshall and Madden (1959) can be reproduced (Fig. 3 ), the one-dimensional Marshall–Madden model was further enhanced to three dimensions using an axial symmetric approach (Model 1, Fig. 1a). The axial symmetric results were verified by three-dimensional modeling with absolute agreement of the phase spectra (Fig. 4 ).

The frequency-dependent Eq. [17], [18], [19], and [13] allow a more direct and effective calculation of the resistivity's amplitude and phase spectra than the time-dependent formulation of Eq. [10], [11], [12], and [13]. For the time-dependent equations, the resulting sine current response of the pore system to an applied sinusoidal voltage must be analyzed time dependently for each frequency step, whereas in the frequency-dependent formulation the frequency-dependent complex amplitudes cp, cn, and u already contain the information about amplitude and phase. As expected, the approximation on neglecting small product terms in the derivation of Eq. [17], [18], and [19] did not make a noticeable difference in the spectra compared with the time-dependent calculation (Fig. 5 ). The computing time was reduced from hours to minutes for one SIP spectrum (on a standard personal computer) depending on the required accuracy, number of data points, and geometry. A second advantage of this approach is the decreased scattering of the resulting phase spectra (see Fig. 5).

Starting from the basic geometry of Model 1 (Fig. 1a) with Lpassive = 1 μm, Rpassive = 1 μm, Lactive = 1 μm, and Ractive = 0.1 μm (Table 2), the whole geometry (all radii and lengths) was increased by a factor s between 1 and 100. For the different geometries, the phase spectra kept a constant amplitude of the phase minimum but changed their frequency position according to fmins−2 (Fig. 6 ). Because the relaxation time is coupled to this characteristic frequency by τ ∼ fmin−1, a quadratic dependence on geometric scaling in agreement with the corresponding fundamental theoretical results (e.g., Schwarz, 1962; Titov et al., 2002) was obtained.

Nevertheless, Model 1 is conceptually problematic with respect to the influence of absolute pore dimensions. For example, scaling the geometry by a factor of 10 would lead to new active zones with the dimensions of the former passive zones. This demonstrates the inability of this model concept to provide useful information on absolute geometric pore properties because the size of pores to be in active or passive zones must be postulated. Therefore, the EDL itself was modeled—with the Debye length as the corresponding length scale—in the other two model concepts (Models 2 and 3). For the spherical particle model (Model 2, Fig. 1b), the influence of the radius of the sphere Rsphere was investigated and compared with the model of Schwarz (1962). The volume fraction of electrolyte solution and rock matrix was kept constant. Thus, in this case all geometric properties of the model system (basic values according to Table 2) were varied in equal measure with the exception of the double layer thickness given by Eq. [8], which remained constant. Of course, the volume fraction of the EDL changed due to changes in the radius of the sphere. Logarithmically equidistant values between Rsphere = 1.4 × 10−7 and Rsphere = 1.4 × 10−5 were used. Because the porosity was kept constant at ϕ = 0.5, the model system started with a length of Ltot = 3 × 10−7 m and a radius of Rtot = 1.5 × 10−7 m. Two characteristic phase minima were obtained (see Results and Discussion), which were sometimes difficult to differentiate. With a combined Maxwell–Wagner and Davidson–Cole fit, it was possible to determine two corresponding relaxation times. The influence of this grain size variation on the Davidson–Cole relaxation time is shown in Fig. 7 . The Davidson–Cole fit was used because of the unbalanced shape of the low-frequency peak. The Maxwell–Wagner fit was done according to an explicit expression given by Hanai (1960), neglecting the permittivity effect at very high frequencies as was done in the simulations (see Governing Equations). An approximately quadratic behavior τ ∼ Rsphere1.95 was found (Fig. 7) for this low-frequency (high relaxation time) peak. Although a different model approach was used to describe the double layer, the simulated dependence of the relaxation time on the particle radius qualitatively agreed with the result of Schwarz (1962; Eq. [1]).

Hence, the modeling concept can reproduce known theoretical results and is in agreement with different mathematical formulations of the problem. Although it significantly reduces the computing time with improved numerical stability, the time-harmonic formulation of the SIP problem has some major constraints. First, because the studied quantities are assumed to be time harmonic, it is impossible to implement constant boundary conditions (e.g., surface charges) at the pore walls, which is necessary for studying more realistic pore systems. Additionally, the equilibrium electrolyte concentrations (in contrast to small excess concentrations) must be equal for both kinds of ions. Thus, different concentration characteristics between the different kinds of ions in the EDL as those present in real pore systems cannot be considered within the direct frequency-dependent calculation. Consequently, for an enhanced coupling among electrolyte concentration, surface charges, and ion mobilities, the more laborious, time-dependent formulation has to be used in future studies.

Results and Discussion

The pore space model (Fig. 1c) was used to assess the influence of pore structure and electrolyte composition on SIP phase spectra. Starting with initial values of Llarge = 5 × 10−6 m, Rlarge = 5 × 10−7 m, Lsmall = 5 × 10−8 m, and Rsmall = 5 × 10−8 m (Table 2), the geometric properties were varied separately in the range of 0.1 to 10 times the corresponding basic value. For each study about the influence of one geometric property, e.g., Rlarge, all other properties were kept constant at their basic values according to Table 2. These dimensions were varied as far as possible without principally changing the geometric configuration. The studies on geometric influence were conducted using two different electrolyte concentrations of C = 1 mol m−3 (fluid conductivity σf ≈ 0.01 S m−1) and C = 0.1 mol m−3 (fluid conductivity σf ≈ 0.001 S m−1). Furthermore, the effect of changing the electrolyte concentration systematically in the range from 0.1 mol m−3 to 100 mol m−3f ≈ 0.001–1 S m−1) was studied for the basic model geometry.

For each parameter variation, 30 phase spectra were taken, consisting of 100 logarithmically equidistant frequency steps in the range from 10−3 to 108 Hz. Figure 2 shows an exemplary spectrum for one set of model parameters. The frequency of the absolute minimum fmin and the corresponding phase value ϕmin were studied as a function of the altered model parameter. Hence, a measure of both the characteristic relaxation time and of the polarization strength were obtained. The results were compared with the theories of Marshall and Madden (1959; Eq. [2]) and Titov et al. (2002; Eq. [5]).

Characteristic Phase Shifts in the Frequency Range from 10−3 to 108 Hertz

The simulated electric resistivity spectra in the frequency range from 10−3 to 108 Hz show two different characteristic phase minima (Fig. 2). It is possible to associate them with different physical mechanisms by comparison with a modified model. By replacing Poisson's Eq. [19] with the electroneutrality condition cp = cn in the system of equations, the buildup of charges is completely stopped. Until then, charge buildup has only been limited by the development of a counteractive internal electric field. If this charge buildup is stopped, the remaining difference with an equilibrium of cp = cn = 0 is now the equal buildup of concentration gradients for cations and anions (Fig. 8 ). The effect that replacing Poisson's equation with the electroneutrality condition has on the resulting phase spectra is the disappearance of the main “high-frequency” minimum. Thus, the charge buildup using the unmodified model produces one large minimum in the higher frequency range. The so-called “low-frequency” minimum is usually smaller, however, may have a more complicated shape, and depending on pore geometry, may have a high-frequency part superimposed by the larger space charge effect. Figure 9 shows a phase spectrum with and without space charges allowed for the same model properties as in Fig. 2. In this special case, a small “high-frequency” effect remains after switching off the space charge buildup. Whereas the “high-frequency” minimum can mainly be associated with the charge buildup at the interface between materials with different electric properties (here, ion mobilities and thus conductivities), the “low-frequency” part is not affected by the space-charge effects at all, thus representing a diffusion current effect. The charge buildup using the unmodified model, observed at the interfaces between areas with different ion mobilities, is known as the Maxwell–Wagner polarization, described by the MWBH theory, named after Maxwell (1873), Wagner (1914), Bruggeman (1935), and Hanai (1960). The studies of this “high-frequency” minimum are not discussed here. In this study, the focus was on the diffusion-controlled “low-frequency” minimum (10−3−104 Hz). Diffusion determines the ion transport if a space charge does not exist. This behavior is illustrated in Fig. 8 for the one-dimensional Marshall–Madden model and can be explained as follows: For low frequencies, additional diffusion currents limit the resulting total current. Above a characteristic frequency, there is not enough time to build up the concentration gradients. These concentration gradients were studied by Blaschek and Hördt (2009) for a similar one-dimensional system in the time domain. Their results concur with those of this study, showing similar curves to those of Fig. 8 for different times after the current was “switched on.” Thus, a higher total current density was obtained in the pore system at higher frequencies because of decreasing diffusion currents. Minimal phase values were found, however, in the frequency region where the current density increased most.

It should be mentioned here that for the pore model, a large difference (that means at least a factor of 10–100) in the dimensions of the larger and the smaller pores with a size of a few Debye lengths was necessary to visually separate the two minima. Especially for geometric values closer together (e.g., values of Llarge = 10−6 m, Lsmall = 10−6 m, Rlarge = 10−6 m, and Rsmall = 10−7 m), the “low-frequency” minimum became very small and was shifted to higher frequencies, so that it was superimposed on the “high-frequency” minimum. The large differences in the pore dimensions combined with very small smaller pores might possibly be associated with recent studies concerning a contribution from the roughness of pores to SIP spectra found at intermediate frequencies (Leroy et al., 2008). For the spherical particle model (Model 2), the effect was observed as well (see Preliminary Discussion about Model Verification and Concept).

Regarding the pore space model (Model 3, Fig. 1c), the results of the influence of pore and electrolyte properties on the phase spectra are summarized in Table 3 .

Dependence of the Phase Shift on Pore Geometry

For the pore space model (Model 3, Fig. 1c), Llarge and Lsmall as well Rlarge and Rsmall were varied separately to study their effect on the characteristic phase spectra. Studying the dependence of the frequency fmin (the position of the minimal phase shift; see Fig. 2) on the geometric pore properties, the relationship fminLlarge−1.7…−1.8 was obtained. This behavior is very similar to the quadratic dependence of the relaxation time on geometric properties proposed by theory (Eq. [1] and [7]). Nevertheless, it was associated with the large pores instead of the small ones (Fig. 10 and 11 ). The frequency of the minimal phase shift fmin did not significantly depend on the pore radii (Fig. 12 ), which is in good agreement with the SNP model's prediction. Consequently, an important result of the simulations is that the only significant parameter influencing fmin was Llarge (Fig. 10). As pointed out above, different experiments have suggested different relations between relaxation time and pore throat size, such as exponents between 1 and 3, logarithmic relationships, or a vanishing influence for very small pore throats (Binley et al., 2005; Scott and Barker, 2005; Kruschwitz, 2007; Kruschwitz et al., 2010; Titov et al., 2010), or even constant phase angle behavior (Börner et al., 1996). In this study, the simulations indicated an importance of the size of larger structures connected to very small pore throats of a few Debye lengths, which by themselves do not influence the phase peak position. On the other hand, these very small pore throat passages linked to much larger pores are necessary to cause noticeable polarization effects in the model system. The corresponding minimal phase shifts ϕmin, as a measure of polarization strength (for definition, see Fig. 2), did not agree with the theoretical models of Marshall and Madden (1959) or of Titov et al. (2002), who predicted large phase values for relatively larger “small pores” (active zones, Fig. 13 ) and relatively smaller “large pores” (passive zones, Fig. 14 ). For the simulations in our study, there was an absolute minimum of the characteristic phase shift at a very large pore length Llarge relative to the other geometric properties (Llarge > 10−5 m, Fig. 14). The amplitude of ϕmin depended on the Rlarge as well. Polarization was maximal for a ratio Rlarge/Rsmall between 5 and 10 (Fig. 15 ) if the larger radius was varied. While the larger pores must be comparatively large for maximal polarization effects, a decrease in the absolute phase value was obtained for an increasing Lsmall (Fig. 13). Considering the radius of the smaller pore, an absolute phase minimum (maximal polarization) was found for a very small radius Rsmall of approximately 1.5 Debye lengths (where there still is a small “channel” of free electrolyte in the center of the smaller pore). For a ratio of RsmallD > 10, the amplitudes fell below 1 mrad and were hardly measurable. This is in contrast to the SNP model, predicting that larger pore throats would cause noticeable phase values (Fig. 16 ). Hence, maximum polarization was obtained for very small structures linked to comparatively large pore bodies.

Dependence of the Phase Shift on the Electrolyte Concentration

Again regarding the pore space model (Model 3, Fig. 1c), the effect of varied electrolyte concentrations (fluid conductivities) was studied. Variations in the electrolyte concentration did not change fmin in this model system, which is in agreement with theoretical predictions (Fig. 17 ). The frequencies themselves were highly varied due to a different geometric influence for the different models (see Dependence of the Phase Shift on Pore Geometry), so unfortunately, the (often slight) variations in relaxation time with changing fluid conductivity, which were observed in experimental studies (Scott and Barker, 2005; Kruschwitz, 2007; Titov et al., 2010) cannot be explained with the present model concept.

The corresponding phase value ϕmin, on the other hand, shows a clear dependence of ϕmin ∼ −C−1 in the pore space model (Fig. 18 ). This dependency agrees with neither the Marshall–Madden nor the SNP model nor the enhanced Marshall–Madden model (Model 1, Fig. 1a), which did not predict any dependence on the concentration at all. A decreasing SIP effect with increasing electrolyte concentration does concur, however, with most experimental results. Thus, many researchers (Lesmes and Frye, 2001; Scott and Barker, 2005; Titov et al., 2010) have studied a normalized chargeability—the product of chargeability and fluid or rock conductivity—instead of the chargeability itself to eliminate this major influence. Although further minimum phase changes cannot be simulated by altering the fluid conductivity, this major dependency can be reproduced. The dependence of fmin and ϕmin on all studied pore model properties is summarized in Table 3.

Conclusions and Outlook

The frequency-dependent complex resistivity of three different microscale rock models was studied by numerically describing ion transport in the pore fluid and at the fluid–rock interface in an electric field. In this study, the EDL was considered by means of a reduced cation mobility within a Debye length distance from the solid–liquid interface. A time-harmonic approach allowed a direct frequency-dependent calculation of the resistivity phase spectra, and an axial symmetric representation of the three-dimensional model geometry led to further reductions of the computational effort. An enhanced Marshall–Madden model and a simple spherical grain model were used for model verification by comparison with theoretical and earlier numerical results. The quadratic dependence of the relaxation time on geometric dimensions proposed by these results was confirmed. Because in this study the pore sizes were assumed to be the decisive parameter, a pore space model was used to study the influence of pore size and electrolyte properties on the SIP phase spectra.

Two characteristic phase minima in a larger frequency range from 1 mHz to 100 MHz were found. They can be associated with different physical mechanisms, which are diffusion effects for the low-frequency behavior and charge buildup in the high frequency range. The low-frequency minimum below 104 Hz was studied here.

This study confirms the dependence of the frequency of the minimal phase shift on a geometric length scale. This length scale was identified as the pore size, more precisely the length in the direction of the electric current of pores linked to very small pore throat passages in the current path. The characteristic frequency fmin decreased with an increasing larger pore length with exponents of −1.7 and −1.8. Furthermore, this length Llarge was the only significant geometric parameter influencing fmin. This relation qualitatively concurred with earlier results, but was associated with the larger pores instead of the small pore throat. These small pore throats, on the other hand, were decisive to cause considerable polarization strength. A ratio of RsmallD ≈ 1.5 between the radius of very short pore constrictions and the Debye length combined with very long passages with a radius Rlarge with values between five and 10 times Rsmall caused maximal SIP phase values. For Rsmall > 10λD (≈0.01 μm for an electrolyte concentration of C = 1 mol m−3), the phase amplitudes fell below 1 mrad, thus becoming hardly noticeable. Consequently, the frequency of the phase minimum can be regarded as a measure for the length of pores that are in contact with sufficiently small narrow passages in the direction of an electric field.

The effect of varied electrolyte concentrations (fluid conductivities) was studied with the important result that the phase magnitude ϕmin showed a clear dependence of ϕmin ∼ −C−1 on electrolyte concentration C, thus showing the major influence of fluid conductivity on polarization strength known from experimental observations. Second, widely concurring with experimental results and in agreement with earlier theories, the electrolyte concentration did not affect the characteristic frequency fmin at all.

Consequently, combining the information obtained by the characteristic frequency and polarization magnitude, the inner pore structure might be inferred from SIP measurements. An improved permeability estimation will be aimed at by considering these relations in further studies. Moreover, more realistic pore geometries will be simulated. The parameterization of the double layer by reduced ion mobilities will be replaced with a more realistic description of the processes close to the surface by an implementation of surface charges. A rock model composed of spherical particles as a transition between pore space and mineral grain model concepts will be simulated to separate the influence of grain and pore sizes. The model concept will also be applied to study the influence of additional liquid phases, such as air or oil in the pore space, on the SIP properties.

This study was conducted within the Transregional Collaborative Research Centre 32 (SFB TR 32). We thank the German Research Foundation for funding. Present and future studies are supported by the Deutsche Gesellschaft für Erdöl, Erdgas und Kohle e.V. (DGMK).