# Quantitative Evaluation of the Effect of Pore Fluids Distribution on Complex Conductivity Saturation Exponents

Siyuan Qiang, Xiaoqing Shi, André Revil, Xueyuan Kang, Yalin Song, Kun Xing

► To cite this version:

Siyuan Qiang, Xiaoqing Shi, André Revil, Xueyuan Kang, Yalin Song, et al.. Quantitative Evaluation of the Effect of Pore Fluids Distribution on Complex Conductivity Saturation Exponents. Journal of Geophysical Research: Solid Earth, 2024, 129, 10.1029/2024jb028689. hal-04729125

HAL Id: hal-04729125

https://hal.science/hal-04729125v1

Submitted on 9 Oct 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
# JGR Solid Earth

## RESEARCH ARTICLE

10.1029/2024JB028689

**Key Points:**

* A micromodel setup allows for simultaneous pore fluid visualization and complex electrical conductivity measurement

* The two saturation exponents for the in-phase and quadrature conductivities are correlated with the ganglia size of the insulating phase

* A power-law-type relationship between the saturation exponents and the change rate of water connectivity with saturation is demonstrated

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

X. Shi,
shixq@nju.edu.cn

**Citation:**

Qiang, S., Shi, X., Revil, A., Kang, X., Song, Y., & Xing, K. (2024). Quantitative evaluation of the effect of pore fluids distribution on complex conductivity saturation exponents. *Journal of Geophysical Research: Solid Earth*, 129, e2024JB028689. https://doi.org/10.1029/2024JB028689

Received 8 JAN 2024
Accepted 3 AUG 2024

**Author Contributions:**

**Conceptualization:** Siyuan Qiang,
Xiaoqing Shi, André Revil

**Data curation:** Siyuan Qiang,
Xiaoqing Shi, André Revil, Xueyuan Kang

**Formal analysis:** Siyuan Qiang,
Xiaoqing Shi, André Revil

**Funding acquisition:** Xiaoqing Shi

**Investigation:** Siyuan Qiang,
Xiaoqing Shi, André Revil, Xueyuan Kang

**Methodology:** Siyuan Qiang,
Xiaoqing Shi, André Revil,
Xueyuan Kang, Yalin Song, Kun Xing

**Project administration:** Xiaoqing Shi

**Resources:** Siyuan Qiang, Xiaoqing Shi

**Software:** Siyuan Qiang, Xiaoqing Shi,
Xueyuan Kang

© 2024. The Author(s).
This is an open access article under the
terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

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

## Quantitative Evaluation of the Effect of Pore Fluids Distribution on Complex Conductivity Saturation Exponents

Siyuan Qiang¹ ID, Xiaoqing Shi¹ ID, André Revil² ID, Xueyuan Kang¹ ID, Yalin Song¹, and Kun Xing¹

¹Key Laboratory of Surficial Geochemistry of Ministry of Education, School of Earth Sciences and Engineering, Nanjing University, Nanjing, China, ²Université Grenoble Alpes, Université Savoie Mont Blanc, CNRS, UMR CNRS 5204, EDYTEM, Le Bourget du Lac, France

**Abstract** The induced polarization (IP) method holds a strong potential to better characterize the critical zone of our planet especially in areas characterized by multi-phase flow. Power-law relationships between the bulk, surface, and quadrature conductivities versus the pore water saturation are potentially useable to map the subsurface water content distribution. However, the saturation exponents *n* and *p* in these power-law relationships have been observed to vary with the texture of geomaterials and the wettabilities of pore fluids. Traditional experimental setups in the laboratory do not allow to independently visualize the pore fluid distribution. Therefore, the physical interpretations of the two saturation exponents have remained unclear. We developed a novel milli-fluidic micromodel using clay-coated glass beads that exhibit excellent visibility and high IP response. Through laboratory experiments, we simultaneously determined the micromodel complex conductivity and acquired the corresponding pore-scale fluid distributions generated by drainage and imbibition through such class of porous materials. Finite-element simulations of complex conductivity based on the upscaling of the complex surface conductance of grains were conducted to determine the saturation exponents under ideal pore fluid distributions. Results indicate that saturation exponents *n* and *p* vary depending on the ganglia size of the insulating fluids. The saturation exponents *n* and *p* exhibit power-law relationships with the change rate of pore water connectivity with saturation, which is calculated through the computation of the derivative of Euler characteristics. These findings provide a new physical explanation to the relationships between the saturation exponents and the microscopic fluid distributions within the geomaterials.

**Plain Language Summary** Water saturation of porous bodies can be related to the complex conductivity through power-law relationships. The existence of these power-law relationships has been clearly documented in the literature. They are critical in the realm of hydrogeophysics to better characterize the critical zone of the solid Earth. However, the values of the saturation exponents in these power-laws have been observed to vary with different material textures and no underlying mechanisms have been able to explain these variations to date. This lack of physical understanding could limit the applicability of the induced polarization method to characterize the critical zone especially when immiscible fluid phases are present. To tackle this, we developed a milli-fluidic pore model that allows investigating the saturation exponents while monitoring the pore-scale fluid distributions. Together with numerical simulations for ideal fluid distribution cases, we found a relationship between the saturation exponents and a microscopic pore parameter called “change rate of the pore water connectivity with saturation.” These findings suggest that when estimating subsurface water saturation from electrical parameters, taking the pore fluid distribution into account can significantly improve the estimation accuracy, therefore enhancing the efficiency of geo-electrical applications for a better characterization of hydrocarbon contaminated aquifers and oil reservoirs.

## 1. Introduction

In the last two decades, Induced Polarization (IP) has been extensively used for non-intrusively characterize and monitor the critical zone of the earth (e.g., Binley et al., 2015; Mendieta et al., 2021; Panwar et al., 2021). For field conditions, IP could be performed either in the time domain (square wave current pulses injected) or frequency domain (continuous alternating current injected, see Zonge et al., 1972), the two kinds of measured data can be converted to each other using Fourier transforms (Fiandaca et al., 2013).

We focus on frequency-domain IP which is also known in geophysics as spectral induced polarization (SIP). The application of SIP in the field allows to obtain tomogram of the complex conductivity $\sigma^*$ (in S/m) of the
**Supervision:** Xiaoqing Shi

**Validation:** Siyuan Qiang, Xiaoqing Shi,
André Revil, Xueyuan Kang

**Visualization:** Siyuan Qiang,
Xiaoqing Shi, André Revil

**Writing – original draft:** Siyuan Qiang,
Xiaoqing Shi

**Writing – review & editing:**

Siyuan Qiang, Xiaoqing Shi, André Revil,
Xueyuan Kang, Yalin Song, Kun Xing

subsurface at a given set of frequencies (typically inside the frequency range 0.1–100 Hz, e.g., Binley & Slater, 2020). Complex conductivity includes an in-phase component $\sigma'$ (in S/m) and a quadrature (out-of-phase) conductivity component $\sigma''$ (in S/m). They characterize the conduction and polarization of the porous material in a given thermodynamic state, respectively. The analysis of the spectra can be used to determine a main relaxation time or a distribution of relaxation times as well (e.g., Revil & Florsch, 2010).

The in-phase conductivity $\sigma'$ includes bulk and surface components, which represent the ability of a geomaterial to conduct electrical current along its pore network and internal interfaces, respectively (Revil et al., 2012). The quadrature conductivity $\sigma''$ describes the ability of geomaterial to store reversibly electric charges (Vinegar & Waxman, 1984). Both in-phase and quadrature conductivities are sensitive to the saturation of the liquid water phase and its salinity (Revil, 2013; Schmutz et al., 2010; Weller & Slater, 2012). They are also sensitive to the pore texture (Maineult et al., 2017), pore fluids morphology (Johansson et al., 2015), and electrochemical characteristics of the mineral surface (Revil et al., 2017).

For porous media saturated by two immiscible fluid phases, the petrophysical model developed by Waxman and Smits (1968) and Vinegar and Waxman (1984) establishes the existence of power-law relationships connecting their in-phase and quadrature conductivities to the water saturation $S_w$ (−). The so-called saturation exponents $n$ and $p$ in these power-laws (Archie, 1942; Vinegar & Waxman, 1984) for the in-phase and quadrature conductivity, respectively are usually determined by fitting the experimental complex conductivity data at different saturations (Han et al., 2021; Schmutz et al., 2012). A general trend of increasing conductivities with water saturation predicted by these power-law relationships agrees with most laboratory and field observations for water-wet geomaterials (water has a preferential wettability for the grain surface, see Archie, 1942; Schmutz et al., 2010; Vinegar & Waxman, 1984). However, researchers have observed that $n$ and $p$ (defined as the change rate of electrical conductivity with saturation) could vary with different pore textures (Revil, 2013; Sun & Torres-Verdín, 2022), wettabilities (Qiang et al., 2022; Revil et al., 2011), pore fractal dimensions (Chen et al., 2023), and pore fluids distributions at a given texture of porous material (Han et al., 2021; Maineult et al., 2018). As small errors in the values of these saturation exponents could lead to deviations in the estimated water saturation, these uncertainties in saturation exponents may lead to unrealistic interpretations of the geoelectrical data.

The variations of the saturation exponent $n$ for different geomaterials are attributed to the microscopic fluid distributions dictated in turn by the texture of the pore network and the history of drainage or imbibition through the pore networks when two immiscible fluid phases are present (Maineult et al., 2018; Sun & Torres-Verdín, 2022). For instance, Anderson (1986) finds that the saturation exponent $n$ becomes extremely large when the portions of brine are trapped, rendering them unable to contribute to electrical conductivity. Han et al. (2021) establish a relationship between the aspect ratio of non-conductive air bubbles in the pores and the saturation exponent $n$. Qiang et al. (2022) also observe a positive correlation between the size of the insulating non-aqueous phase ganglia in the pore space and the saturation exponent $n$. Tian et al. (2023) and Zhang et al. (2021) find that the saturation exponents $n$ is controlled by the pattern of the growth of the insulating hydrate phase in the pore space. Zhu et al. (2023) also attribute the high value of $n$ to the blockage of capillaries within the sediment by the development of hydrate. These findings indicate that saturation exponents are affected by the pore fluids distribution and how they contribute to the connectivity of the current path through the conductive fluid phase.

Furthermore, Glover (2017) proposes a model interpreting the saturation exponent $n$ as related to the rate of change of the factional connectedness with saturation and connectivity. Binley and Slater (2020) also recognize that the saturation exponent is attributed to the changes in pore water connectivity. Recently, Qiang et al. (2022) suggest a simpler positive correlation between the saturation exponent $n$ and the change rate of pore water connectivity with saturation through pore scale experiments and simulations. These theories further attribute the physical interpretation of the saturation exponent $n$ to the change of pore water connectivity, rather than the pore water connectivity itself. However, the concept of “connectivity” in previous studies is qualitatively defined, and none of these concepts are verified through independent laboratory experiments or numerical simulations. This is especially true for the saturation exponent $p$ of the quadrature conductivity.

One main obstacle in elucidating the physics of saturation exponents is the lack of technology to acquire the pore fluid distribution and the corresponding IP response simultaneously. Recently, micro- and milli-fluidic techniques have gained popularity in studying the geoelectrical method at the pore scale (Fernandez Visentini et al., 2023; Rembert et al., 2023). Compared with the opaqueness of traditional column or sandbox experiments, the excellent visibility and controllability of micromodels could help to elucidate the relationship between electrical signals and
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Sketch of the experimental setup. (a) Picture of the experimental micromodel without clamps. (b) Top view of the sketch of the micromodel. The dash lines correspond to the cross-section position for (c) the cross-section of the micromodel in line c' and (d) the cross-section of the micromodel in line d'.</figcaption></figure>

the microscopic properties of porous media. Related pore scale direct current (DC) studies focus on the sensitivity of in-phase conductivity to factors such as small-scale saline tracer heterogeneity (Jougnot et al., 2018), flow patterns and porosity efficiency (Sun et al., 2021), and diffusion-limited mixing of concentration (Fernandez Visentini et al., 2023). Nevertheless, obtaining accurate and repeatable IP measurements (including in-phase and quadrature conductivity) in micromodels faces significant challenges due to the minimal electrical signal-noise ratio compared to DC methods (see Kemna et al., 2012; Zarif et al., 2017). Izumoto et al. (2022) and Rembert et al. (2023) investigate the SIP signal of calcite precipitation and dissolution using well-designed milli- and micro-fluidic setups. However, the lack of a pore-scale SIP experimental setup that exhibits significant IP response of the porous medium itself has hindered quantitative SIP investigation for partially saturated conditions. Because decreasing saturation typically results in a decrease in the already inconspicuous quadrature conductivity (Vinegar & Waxman, 1984).

We explore below the relationships between pore fluids distribution and the saturation exponents *n* and *p*. We first develop a unique and novel milli-fluidic micromodel dedicated to complex conductivity measurements and real-time pore fluids imaging. A series of pore-network complex conductivity simulations are also conducted to check and extend the experimental results. Based on the change of Euler characteristics of pore water, we propose a novel quantitative relationship between the change rate of pore water connectivity with saturation $ΔC$, and the saturation exponents.

## 2. Materials and Methods

### 2.1. Experimental Setup

#### 2.1.1. Micromodels

In our experimental work, we aim to acquire different pore fluid distributions and determine their complex conductivity (see Text S1 in Supporting Information S1 for the introduction of spectral induced polarization method). Figure 1 shows the micromodel employed in the experiments. The bottom plate of the micromodel is a Plexiglas plate against a layer of soft acrylic acid (thickness 0.5 mm). The plate is uniquely shaped as a “spindle” to facilitate fixation. A 1 mm thick layer of silicone with a spindle-shaped cavity is then placed on the acrylic acid layer to define the in-plane boundaries of the porous medium (77.5 mm × 25 mm). On one of the long sides of the silicone, two 6 mm × 3 mm rectangle cavities (25 mm apart) are extended for the potential electrodes (see
Figure 1d). When the bottom setting is ready, glass beads coated with thin Ca-montmorillonite layer (thickness 30.83 ± 20.95 µm, which is obtained by estimating the thickness of highlighted dark Ca-montmorillonite layer surrounding the grains in 2D images) are randomly thrown into the cavity and attached upon contact with the acrylic acid on the exposed bottom plate. The total diameter of glass beads with clay coating is approximately 1.036 ± 0.033 mm (estimated directly from the 2D images). After that, the system is thoroughly shaken to ensure the glass beads form a “monolayer.” Finally, another Plexiglas plate with two holes for fluids inlet and outlet, four holes for electrodes, and two slots on the subsurface for the current electrodes is placed on the top of the setup. Clamps are used to clip all the layers and form a closed porous medium.

The pressure from clamps compresses the soft acrylic acid on the bottom plate and the surrounding silicon glue layer, making the monolayer of glass beads a tight “Hele-Shaw” structure with dual porosity. The macro porosity (i.e., the porosity inter glass beads) of the micromodel (0.5774 ± 0.042) is estimated by dividing the overall bulk volume of 1,403 clay coated glass beads (diameter 1.036 ± 0.033 mm) by the volume of the central rectangle area in the micromodel (77.5 × 25 × 1 mm). The verification of the porosity is made by dividing the maximum water volume that the central rectangle area of the micromodel could hold (1.16 ml) by the total volume, resulting in a value of 0.598. The contribution of the internal porosity of the clay coating on the glass beads $\phi_{ca}$ to the total porosity could be ignored considering the small volume proportion of the clay coatings, but note that $\phi_{ca}$ would affect the equivalent surface charge density of the clay coated glass beads, which will be discussed later in Section 2.3.3.

### 2.1.2. Porous Medium

Studying the relationship between the SIP signal and the pore fluid distributions requires both a good visibility and detectable IP response. Achieving simultaneously both requirements is challenging in itself. We finally addressed these issues with a porous medium consisting of clay-coated glass beads. As packs of pure glass beads exhibit very small SIP responses (phase angle smaller than 1.1 mrad when saturated with pore water conductivity $5 \times 10^{-3}$ S/m at 25°C, see Leroy et al., 2008). Such low phase value is expected to strongly reduce the signal-noise ratio of the SIP measurements in our micromodel. To amplify the SIP response of the porous medium, we utilize a modified coating technique (Jerez et al., 2006) by using highly polarizable clays in coating the surface of the glass beads. To this end, we first added 2 g of polyvinyl alcohol liquids (PVAL, Deli) to 40 g of glass beads in a beaker and stirred thoroughly to attach the glass beads surface with a layer of PVAL. Then we added 10 g of Ca-montmorillonite (purity 99%) to the mixture gram by gram while constantly stirring, as smectite is known to have a strong polarization because of its high cation exchange capacity (80–120 meq/100 g, see El Alam et al., 2023, and references therein). After which, the mixture of glass beads, PVAL, and clay was transferred to the oven, dried at 80°C for 12 hr, and washed with distilled water to remove the dissociative clay particles. The mixture was then dried at 80°C for another 12 hr, washed and dried again. Following these steps yields a polymer bridging that bonds the glass beads (silica surface) with the clay minerals through the polar functional groups in the non-ionic PVAL (Jerez et al., 2006). The thin clay coating layer is not uniformly distributed on the glass beads' surface (the opaque area on the glass beads exhibits the clay coating, see for instance Figure 1a), but stable to resist the water flow. Notably, the artificial configuration of this clay coating glass beads is consistent with some natural geomaterials such as chlorite-coated sandstones, where dispersed clay is distributed on the surface of the silica grains (Rabaute et al., 2003; Revil et al., 2015). Similar clay-coated glass beads have also been used as porous media for the studies of pollutant adsorption (Ake et al., 2003) and biodegradation (McBride et al., 1992).

During flow experiments, hydration and clay swelling could occur when the clay coatings contact with the fluid. However, we do not observe significant physical changes in either the clay coating or the pore structures during flooding. This is also proved by Aksu et al. (2015) that the porosity and permeability of clay-coated glass beads remain almost unchanged in the case of high porosity porous media. Thus, the pore structure change due to clay-related reactions is neglected for our experimental porous medium, which is characterized by a very high porosity (0.5774 ± 0.042).

Another important reaction between clay and pore fluid is the formation of an Electrical Double Layer (EDL) on their interface due to the surface electrical charge density at the surface of clay minerals (Miller & Low, 1990). The EDL includes the Stern layer of absorbed counterions and the diffuse layer (Revil et al., 2017). EDL polarization, as the dominating polarization mechanism in sand-clay mixtures (Revil et al., 2012), is expected to be the most important contribution of our measured SIP signals of the micromodel. Compared to the low SIP signal
of pure glass beads, the clay coatings exhibit significant SIP responses attributed to their high surface charge density and strong EDL polarization. This equivalent chargeable layer on the glass bead surface significantly improves the signal-noise ratio of the measurements (phase angle larger than 4 mrad when the model is saturated with pore water conductivity 0.02 S/m, 25°C) while preserving good media visibility for the camera.

The pore fluid distribution in this clay-coated porous medium may be different from that in the common fine-
grained geomaterials due to differences in the grain and pore size distribution. However, the fundamental
electrical mechanism of the clay-coated glass beads at play in our experiments is consistent with geomaterials
where clays-rich zones are coated on the surface of the rock matrix (Niu & Revil, 2016). The general electrical
behaviors of this porous medium with different pore fluid distributions could have implications for common
geomaterials. Similar porous media made by glass bead monolayers have also been used in monitoring the pore
scale motion of drainage fronts (Moebius & Or, 2014), dynamics of viscous entrapped saturated zones (Li
et al., 2018), and intermittent dynamics of slow drainage flows (Moura et al., 2020).

2.1.3. Electrode Setup

The impedance measurements are made with a four-electrode approach to avoid the effect of electrode polar-
zation. The potential (voltage) electrodes M and N (located on the top plate of the micromodel, see Figure 1) are
made of brass and have a diameter of 3 mm in contact with the electrolyte. To avoid electrode polarization in the
applied electrical field, the potential electrodes are retracted around 5 mm upwards and 3 mm away (Izumoto
et al., 2022; Zimmermann et al., 2008) from the boundary of the porous medium (Figure 1c). The space protruding
from the rectangle porous medium is filled with fine sands (diameter < 0.1 mm) to prevent the insulating non-
aqueous phase from contacting the potential electrodes and blocking the electrical currents. The point poten-
tial electrodes located at the single side of the micromodel may result in a localized voltage difference response of
the area between the potential electrodes, as a result, repeated experiments at similar conditions are needed (see
Text S2 in Supporting Information S1 for more details of the limitations associated with potential electrodes). The
current electrodes A and B are made by 24 mm long and 1 mm thick Ag electrodes connected with brass columns.
The subsurface of the rectangle Ag electrode is retracted to be parallel with the lower surface of the top plates to
prevent it from blocking the flow, as shown in Figure 1c. The geometric factor *K* of the micromodel in Figure 1 is
determined as 9.7 × 10⁻⁴ m by measuring the complex conductivity of micromodels filling with fluid of known
conductivity (Revil & Skold, 2011), which is close to the analytical approximation *K* = 1 × 10⁻³ m provided by
Binley and Slater (2020) for this case (*K* = *A*/L, where *A* = 2.5 × 10⁻⁵ m² is the cross-sectional area of the cell and
*L* = 0.025 m is the length between the potential electrodes). The electrode setup is tested by measuring the
magnitude and phase of the micromodel that was filled with only brine with a conductivity of 0.5 S/m, phase
values are all very close to zero at frequencies below 1 kHz (see Text S3 in Supporting Information S1 for more
details, together with a test of the performance of potential electrodes made by different materials).

## 2.2. Experimental Design

2.2.1. Experimental Procedure

We conduct two groups of experiments, including drainage and imbibition to acquire complex conductivity data for different pore fluid distributions. Here, we do not discuss the difference between the drainage and imbibition processes. Instead, we focus solely on the static fluid distributions when the flow system becomes steady. In the drainage experiments, we first slowly fill the micromodel with dyed tap water (dyed by red ink that contains brilliant crocein) with a syringe pump (LSP01-2A, Longerpump). Dying of the tap water cause a conductivity increase from 0.034 to around 0.04 S/m. Then we inject dimethyl silicone oil (CAS:63148-62-9, viscosity 10 mPa s, Macklin) as the non-aqueous phase liquid (NAPL) into the micromodel at different injection rates (0.25–20 ml/min) until the drainage became “stable” (i.e., the phase distribution in the micromodel did not change significantly). The SIP response of the static micromodel in a frequency range of 0.091 Hz–40 kHz is measured using SIP-LAB-IV (Radic-Research), and images of the micromodel are taken by CCD camera (Imaging, 4096 × 3000 pixels, RGB). The imbibition experiments begin with different initial phase distributions acquired by the drainages or another imbibition. A larger range of water injection rate (0.25–40 ml/min) is chosen to acquire a wider steady-state saturation value range.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** An example of the processing procedure of images taken from CCD camera. The original image taken from the CCD camera (a) is first separated into red and green (blue) channels (b and c). The difference between the binarized images of red and green (blue) channels (d and e) indicates the water phase distribution (f). Median filter is used to remove the isolated noise pixels in the images and the final image (g) provides the pore water distribution.</figcaption></figure>

### 2.2.2. Estimating Water Saturation

Due to the conditions of the experiments (especially at high injection rates), we could not always get a good estimate of the water saturation directly from the injected volume of the fluids. The water saturation is then determined from the images of the steady phase distributions by first acquiring the 2D distribution of the water phase. Figure 2 shows the workflow of CCD image processing. During the process, the main difficulty is that the grain phase is not uniform in the pixel value due to the uniform clay coating, some part of the glass bead is not covered by dark clays and others were completely opaque. To deal with this, the RGB images taken from the steady state of drainage and imbibition are separated using MATLAB (Higham & Higham, 2016) into red and green (or blue) channels. Among them, the green (or blue) channel contains both the opaque grain and red pore water pixel information, while the red channel only identifies the opaque grains because the red signals representing the dyed water are filtered. We binarize the red and green channel images separately using the “im2bw” function based on Ostu's method (Ostu, 1979), and the difference between the two binarized images is regarded as the water distribution. We then use a median filtering function “medfilt2” based on the median value in a 5-by-5 neighborhood pixels around the corresponding pixel to remove the isolated noise pixels in the final image (Lim, 1990). For more details of the processing procedure, and potential errors from the median filter that merge some small water droplets or separated solid phase, please see Text S4 in Supporting Information S1.

Figure 2g provides the 2D water distribution used to estimate the water saturation. The saturation is then determined by calculating the ratio of the area of the water phase in each image to the area of pore space, the latter is obtained from the 2D image of a fully saturated micromodel. Similar saturation estimation technology has also been used by Moura et al. (2020). Note that using the 2D approximation of the micromodel may introduce errors in estimating saturation due to the intrinsic 3D spherical nature of the glass beads (refer to Text S5 in Supporting Information S1 for more details). Thus, we validate the saturation estimation technique by slowly injecting known volumes of water into the micromodel. Figure 3 shows the estimated saturation from 2D images versus the saturation determined as the ratio of the real injected water volume into the central rectangle area of the micromodel to the water volume that could just saturate the central rectangle area of the micromodel. The latter represents the pore space volume (approximately 1.16 ml). The saturation estimations data have a root mean square error (RMSE) of 0.031 in the saturation range from 0.21 to 0.84.

## 2.3. Numerical Methods

### 2.3.1. The Complex Conductivity of Porous Media

The measuring error and lack of experimental data in controlled conditions may prevent further discussion on the complex conductivity, numerical simulations are needed to reveal the fundamental physics of the saturation exponents. For complex conductivity simulation, we build a 2D porous medium by randomly adding circle grains
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Comparison of the estimated and calculated saturation. The y-axis denotes the estimated saturation from the ratio between the water phase area and pore space area in the 2D images. The x-axis denotes the saturation determined as the ratio between the known volumes of injected water into the central rectangle area of the micromodel to the maximum water volumes that the central rectangle area of the micromodel could hold. The data have a root mean square error (RMSE) of 0.031 in the saturation range from 0.21 to 0.84.</figcaption></figure>

(radius 0.5 mm) to a 77.5 mm × 25 mm domain. The grain positions are randomly selected but followed a minimum interval of 1 mm to each other to make the grain distributions consistent with the experiments. The addition of grains continued until the medium reached a porosity $\phi = 0.595$, which is close to the porosity of the experimental micromodel.

Here, we use a numerical method that calculates the complex conductivity of the grains by upscaling the complex conductivity of the EDL (Revil, 2012) on the grain surface to the whole grain (Niu & Zhang, 2017; Niu et al., 2020). The complex conductivity $\sigma^*$ of saturated porous medium could be written as:

$$ \sigma^* = \frac{\sigma_w}{F} + \sigma_s^* \quad (1) $$

where $\sigma_w$ (in S/m) denotes the conductivity of the pore water, $F = \phi^{-m}$ (-) denotes the (intrinsic) formation factor, $m$ (-) corresponds to the cementation or porosity exponent. The complex-valued quantity $\sigma_s^*$ (in S/m) denotes the complex surface conductivity associated with conduction and polarization of the EDL coating the grain surfaces. To acquire $\sigma_s^*$, we first calculate the surface conduction of the diffuse layer $C_S^d$ (in S) in EDL on the grain surface based on the simple assumption used by Revil (2013):

$$ C_S^d = \beta_{(+)}^d Q_d \quad (2) $$

where $\beta_{(+)}^d = 5.2 \times 10^{-8}$ (in m$^2$ s$^{-1}$ V$^{-1}$) (Revil, Florsch, & Camerlynck, 2014) denotes the ion mobility related to diffuse layer, $Q_d$ is the surface charge density and chosen to be 0.32 (C m$^{-2}$) (Woodruff & Revil, 2011), this value is much smaller than the one for the clay coated glass beads in the experiments (which will be discussed later in Section 2.3.3), but could be found in natural clay materials (Revil et al., 1998, 2013).

The complex conductance of the Stern layer $C_S^s$ (in S) is calculated based on the theory of Stern layer polarization (Leroy & Revil, 2009) following a Warburg distribution. Which is a special case of the Cole–Cole model (Tarasov & Titov, 2013) with the Cole–Cole exponent c fixed to 0.5 (Revil, Kessouri, & Torres-Verdín, 2014):

$$ C_S^s = \frac{i\omega\tau^{1/2}}{1 + i\omega\tau^{1/2}} \Sigma_S^s \quad (3) $$

where $\Sigma_S^s = \beta_{(+)}^s Q_s$ (in S) is the surface conductance related to Stern layer. As we are not discussing the detailed polarization mechanisms, we assume the ion mobility and surface charge density of Stern layer $Q_s$ and diffuse layer $Q_d$ are the same for simplification. $\omega = 2\pi f$ (in rad s$^{-1}$) is the angular frequency, where f is the measuring frequency. The relaxation time $\tau = r_d^2/2MD_{(+)}^s$ (Niu & Revil, 2016). Here we took the characteristic radius $r_d = 5 \times 10^{-4}$ m as identical to the grain radius in the micromodel, diffusion coefficient of the counterions in the Stern layer $D_{(+)}^s = 3.8 \times 10^{-12}$ m$^2$ s$^{-1}$ and coefficient $M = 3$ (-), which is the maximum value reached by M at high surface charge densities (>0.16 C m$^{-2}$, Niu & Revil, 2016). These parameters result in a relaxation time of 11 × 10$^4$ s for the porous medium we use, indicating that the peak frequency of the phase spectra of our micromodel should fall at around 10$^{-5}$ Hz using the equation $f_p = 1/2\pi\tau$ (Revil et al., 2013). Finally, the complex surface conductivity $\sigma_s^*$ (in S/m) of grain with EDL is electrically equivalent to a whole grain without the EDL by an upscaling procedure (Revil & Florsch, 2010):

$$ \sigma_s^* = \sigma_s' + i\sigma_s'' = \frac{2(C_S^d + C_S^s)}{\Pi} \quad (4) $$

where $\Pi$ is the characteristic particle size and is set to 0.5 mm.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Comparison of the complex conductivity calculated by the finite element method and differential effective medium (DEM) theory. The results are consistent. For the DEM model, the complex conductivity of pore water and grains are set to be consistent with the numerical simulation, the cementation exponent *m* of the 2D simulation domain (porosity *φ* = 0.595) is calculated to be 1.93, which results in a formation factor *F* = 2.68. The black cross and dots indicate the experimental measured complex conductivity under water saturated condition. The difference in the quadrature conductivity of experiments and simulations is related to the high equivalent surface charge density (*Q*<sub>s</sub> ≈ 630 C m<sup>-2</sup>) caused by the chargeable clay layer coating on the glass bead with respect to the surface charge density of a non-coated glass bead in the simulations (*Q*<sub>s</sub> = 0.32 C m<sup>-2</sup>).</figcaption></figure>

### 2.3.2. Simulation Implementations

The numerical simulations of the complex conductivity of the 2D porous medium are conducted thanks to the finite element code COMSOL Multi-physics® (version 5.4). The 2D porous medium domain is fine meshed, and the grids are assigned with grain, water, and NAPL phases according to the simulation schemes. For the numerical solution, we use apparent dielectric coefficients ε<sub>a</sub> to substitute the quadrature conductivity σ'' related to the IP process in the porous medium by ε<sub>a</sub> = σ''/ω. As the dielectric polarization mechanisms usually occur at higher frequencies (Revil et al., 2023), we neglect the dielectric polarization in this work (including the Maxwell-Wagner effect associated with the discontinuity of the current displacement at the interface between the different phases of the composite). So this apparent dielectric coefficient has nothing to do with the intrinsic dielectric permittivity of the porous medium at high frequencies, but the dielectric equivalent of the effect of the IP quadrature conductivity in Ampere's law (Vinegar & Waxman, 1984).

The in-phase conductivity of water phase is set to 0.04 S/m as in laboratory experiments and the apparent dielectric coefficient of water related to IP process is set to 80ε<sub>0</sub> where ε<sub>0</sub> = 8.854 × 10<sup>-12</sup> (in F m<sup>-1</sup>) denotes the dielectric coefficient in vacuum. The background in-phase conductivity and apparent dielectric coefficient of the grain phase were set to 10<sup>-5</sup> S/m and 4ε<sub>0</sub>, respectively (Niu & Zhang, 2017). While for the NAPL phase we use 10<sup>-6</sup> S/m and 2.5ε<sub>0</sub> (Gregory & Clarke, 2012; Jougnot et al., 2018). Those electrical properties remain unchanged over the evaluated frequency range. The add-on in-phase conductivity and apparent dielectric coefficient (quadrature conductivity) for the grain phase are acquired by Equation 4 using σ'<sub>add</sub> = σ'<sub>s</sub> and ε<sub>a-add</sub> = σ''<sub>s</sub>/ω (Xing et al., 2023), which yields frequency-dependent values.

The longer boundaries of the 2D porous medium domain are set as electric-insulation boundaries, and the shorter left boundary is specified as a complex voltage boundary while the shorter right one is the ground boundary. The complex conductivity of the domain is acquired by dividing the applied voltage ΔV by the electrical current I flow through the domain. The geometrical factor K of this domain is 0.3225 m when the effective thickness is fixed to 1 m, which is acquired by setting the domain with a uniform material with known conductivity. This value is aligned accurately with the analytical A/L approximation (0.32258) where the cross-sectional area A = 2.5 × 10<sup>-3</sup> m<sup>2</sup> and the length between the potential (also current) electrodes L = 0.0775 m. The numerical solution of a water-saturated case using this method is compared with the effective complex conductivity calculated by the differential effective medium (DEM) with a MATLAB routine (Bussian, 1983; Niu & Zhang, 2017) (see Figure 4):

$$ \sigma_{\text{eff}}^* = \frac{\sigma_1^*}{F} \left( \frac{1 - \sigma_2^*/\sigma_1^*}{1 - \sigma_2^*/\sigma_{\text{eff}}^*} \right)^m \quad (5) $$

where σ<sub>1</sub><sup>*</sup> and σ<sub>2</sub><sup>*</sup> (in S/m) are the complex conductivity of pore water and grains, respectively, which are consistent with the numerical simulation. *m* is estimated to be 1.93 for this 2D simulated domain, which results in *F* = 2.68 with a porosity *φ* = 0.595. To compare with the experiments, this simulated formation factor is close to the apparent formation factor σ'/σ<sub>w</sub> = 2.31 for the experiment under saturated condition with a pore water conductivity σ<sub>w</sub> = 0.04 S/m. However, the acquisition of the intrinsic experimental formation factor should consider excluding the surface conductivity σ'<sub>s</sub> from the in-phase conductivity σ', which will be discussed later in the result section. As Figure 4 shows, the numerical and analytical solutions in the range of 10<sup>-8</sup> Hz to 10<sup>4</sup> Hz are consistent.

### 2.3.3. Comparisons Between the Simulations and Experiments

In Figure 4, we compare the complex conductivity acquired by the experiments and simulations under water saturated conditions. Due to the prolonged SIP measuring time at low frequencies, the experiment only measures
the complex conductivity between 0.091 Hz and 5 kHz (the measurements at 46.8 Hz are close to the 50 Hz and prone to uncertainty as corresponding to the powerline frequency). The experimental in-phase conductivity is close to both the numerical simulated and DEM calculated results in the log-log plots. Despite the similarity in the frequency dependence observed between the experimental and numerical quadrature conductivities (depression in the experimentally measured frequency range), there exists a notable discrepancy in their absolute value. Specifically, the experimental quadrature conductivity is three orders of magnitude higher than the corresponding numerical results. This extreme magnitude of experimental quadrature conductivity is caused by the dense clay coatings on the glass beads, which form an equivalent chargeable layer with high equivalent surface charge density. The equivalent surface charge density of the clay coated glass beads could be determined as the ratio of total charge of the clay coating layer ($4\rho_g(1 - \phi_{ca})CECh\pi r^2$) to the glass bead surface area ($4\pi r^2$), which is written as:

$$Q_s = \rho_g(1 - \phi_{ca})CECh \quad (6)$$

where CEC (in C kg⁻¹) denotes the cation exchange capacity of the Ca-montmorillonite, $\phi_{ca}$ and h denote the “equivalent” porosity and thickness of Ca-montmorillonite coating, respectively. Given CEC = 77,056 C kg⁻¹ as for the smectite (El Alam et al., 2023), $\phi_{ca}$ = 0.9 (Revil et al., 2018, note that this “equivalent” porosity is only a reference value used for $Q_s$ estimation and do not correspond to the real porosity of the clay coatings) and h = 30.83 µm, we have a surface charge density for the clay coated glass beads around 630 C m⁻², which is far greater than the value used in the simulations (0.32 C m⁻²). This significant surface charge density and resulting IP response are expected for the experiments to make the micromodel detectable by the SIP monitoring setups, as it is difficult to ensure a good signal-noise ratio for SIP monitoring with small IP response. Meanwhile, we chose to keep the lower surface charge density for the simulations to make the numerical situation schemes more consistent with natural geomaterials (Revil et al., 1998, 2013). Despite the difference in surface charge density, there could also be differences in a mass of parameters such as pore geometry, relaxation time, electrode polarization, and the Cole-Cole exponent between the experiments and simulations. These potential differences interact complicatedly with each other and are difficult to determine experimentally, which prevents the complete reproduction of the parameters and absolute IP response of the micromodel experiments by the simulations. However, the numerical simulations could be used to reveal the reasonable relative changes in the complex conductivity with saturation, and further helped us to study the relationships between the saturation exponents and pore fluids distributions of natural geomaterials.

# 3. Results

## 3.1. Experimental Results

With different injected phases, rates, and pre-distributions of the water/NAPL phase, we acquire a total of 23 (6 for drainage and 17 for imbibition) IP spectra and micromodel images with different fluid distributions. Figure 5 shows a part of the in-phase (a) and quadrature (b) conductivity frequency spectra at 0.091 Hz to 5 kHz and all the in-phase (c) and quadrature (d) conductivity at the lowest measure frequency 0.091 Hz versus water saturation.

Figures 5c and 5d show that the in-phase and quadrature conductivities versus saturation at high saturations fall on a straight line in the log-log plots, but the conductivity decrease seems to slow down at lower saturations and deviate from the power-law straight line. Previous works on in-phase conductivity have also observed changes in the slope of the conductivity with saturation in log-log plots (e.g., Deng et al., 2018; Han et al., 2009; Li et al., 2015; Montaron, 2009). Some of them are explained by the percolation theory (Ewing & Hunt, 2006; Kozlov et al., 2012) or residual water films around the grains at low saturations (Han et al., 2009; Sun & Torres-Verdín, 2022). However, for the determination of saturation exponent n of these data, care should be taken in checking that surface conductivity $\sigma'_s$ (Revil, Kessouri, & Torres-Verdín, 2014) has been properly accounted as the in-phase conductivity comprised surface conductivity components (Waxman & Smits, 1968):

$$\sigma'(S_w) = \frac{1}{F} S_w^n (\sigma_w + BQ_v/S_w) \quad (7)$$
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Experimental results. (a) In-phase and (b) quadrature conductivity frequency spectra at 0.091 Hz to 5 kHz, lighter colors indicate high water saturation (the lines are just guides for the eyes). (c) The complete in-phase and (d) quadrature conductivity at 0.091 Hz versus saturation in the drainage and imbibition experiments. The dashed lines in (c) and (d) show the power-law fitting of the imbibition data. Gray data points in (c) indicate the bulk component of in-phase conductivity (the surface component is excluded for saturation exponents estimation). The saturation exponents shown in the lower right corner of (c) and (d) were estimated from solely the drainage, imbibition, and all data series, respectively.</figcaption></figure>

where B (in m² s⁻¹ V⁻¹) denotes the equivalent conductance of the counterions (Revil et al., 2017). Q_v (in C m⁻³) denotes the volumetric charge density (Revil et al., 2013; Waxman & Smits, 1968). The term BQ_v/S_w multiplicated by S_w^n/F describes the saturation dependence of surface conductivity. Here with empirical relationships between the quadrature conductivity σ'' and the in-phase component of surface conductivity σ'_s (Revil et al., 2017; Weller et al., 2013), we could distinguish the surface conductivity from the in-phase conductivity. The bulk component of in-phase conductivity data series in Figure 5c shows the result of excluding the surface components from the in-phase conductivity using the empirical relationship σ''/σ'_s = l = 0.042 (fit value by multiple types of geomaterials, Weller et al., 2013). This value of l yields an experimental (intrinsic) formation factor F = 3.50, which decreases to 2.71 when l = 0.097 (value for saprolite, Revil et al., 2013).

The saturation exponent p for quadrature conductivity is directly related to the slope of the quadrature conductivity with saturation in the log-log plots as the quadrature conductivity is only influenced by the quadrature conductivity under saturated conditions σ''(S_w = 1) and saturations (Vinegar & Waxman, 1984):

$$ \sigma''(S_w) = \sigma''(S_w = 1)S_w^p. \qquad (8) $$
<figure><img src="image_8.png" /><figcaption>**Figure 6.** The in-phase (a) and quadrature (b) conductivity at 0.091 Hz versus the Euler characteristics of pore water acquired by the drainage and imbibition experiments. The dashed lines indicate the power-law fits for all data. R² describes the goodness of fit.</figcaption></figure>

Finally, the estimated saturation exponents in Figures 5c and 5d show a difference between the drainage (low saturation) and imbibition (high saturation) data series. In Figure 6, we plot the in-phase and quadrature conductivity at 0.091 Hz versus the corresponding Euler characteristics of water phase. The Euler characteristic of pore water is defined as the number of connected water phases minus the number of disconnected NAPL and solid phases in an image, which represents the connectivity of the water phase (Armstrong et al., 2019; Mecke & Arns, 2005; Scholz et al., 2012). The Euler characteristic could be calculated expediently through the MATLAB function "bweuler" from the binarized images where the water phase is assigned with "1" and the NAPL and solid phase are assigned to "0." As Figure 6 shows, the in-phase and quadrature conductivity are both fitted with the Euler characteristics by a power-law with an R² > 0.6, which agrees with and extends the relationship between the in-phase conductivity and the Euler characteristics of water $\chi_w(-): \sigma' \sim \chi_w^\beta$ to complex conductivity case, where $\beta(-)$ is the exponent related to phase morphology (Qiang et al., 2022; Scholz et al., 2012).

## 3.2. Simulation Scheme and Results

The experimental results show that the saturation exponents change with the saturation of the micromodel. However, the lack of controllability and accuracy in the experimental data impedes further discussions on the inherent factor that controls these changes. For instance, the Representative Elementary Volume (REV) of the low saturation data points in Figure 5 may be larger than the micromodel itself because the water distribution at low saturation is highly heterogeneous (Jougnot et al., 2018). Which means that the complex conductivity measured at low saturations could vary for different water distributions even with the same saturation. We do not test these potential errors in experiments systematically because we could not control the final saturation and fluids distribution acquired by the drainage or imbibition process. Alternatively, we construct ideal numerical simulations to implement uncontrolled pore fluids distributions and calculate their corresponding complex conductivity at $10^{-5}$ Hz (close to the peak frequency of the simulated SIP phase spectra). The simulations could not completely restore the absolute value of the experiments, but could provide us with the potential evolution pattern of complex conductivity under well-controlled saturation and pore fluids distributions, which help us to further understand and explain the data measured in the experiments.

The simulations change the pore fluids distribution by adding different sizes of NAPL ganglia into the pore space (see Figures 7a and 7b) using the patchy saturation model (Toms et al., 2006), which has become in the last 15 years as a classical method to simulate the pore fluids distribution in porous media (e.g., Lebedev et al., 2009; Mustofa et al., 2022; Xing et al., 2021). The detailed steps are as follows: (1) First, the 2D pore space is initially saturated with water. (2) Subsequently, 625 circular NAPL ganglia (radius 0.2 mm) are randomly added to the pore space. During this process, some NAPL ganglia may be covered by the grains, or overlap with each other and merge into larger ganglion. The complex conductivity of the 2D domain with this water-NAPL distribution is then simulated. (3) Another 625 NAPL ganglia (radius 0.2 mm) are randomly added to the pore space (which has already contained 625 NAPL ganglia after step 2), and the corresponding complex conductivity is simulated again. (4) Step 3 is repeated 9 times, together with the complex conductivity under saturated condition, we acquire
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Sketch of pore space in the simulations. (a) Pore space with small NAPL ganglia. (b) Pore space with large NAPL ganglia. The two pore spaces have a similar saturation.</figcaption></figure>

10 data points for the NAPL ganglia radius of 0.2 mm case at different sat-
urations, as shown in Figure 8 by the line on the top. (5) Steps 1–4 are
repeated for NAPL ganglia radius ranging from 0.2 to 1.4 mm in 0.1 mm
increments. The number of added NAPL ganglia each time in steps 2 and 3 is
adjusted to be negatively proportional to the square of NAPL ganglia radius.
For example, when the added NAPL ganglia radius is 0.5 mm, the number of
NAPL ganglia added to the pore space for each time is 625/(2.52) = 100. This
adjustment ensures the same water saturation range (0.75 to 1) for different
NAPL ganglia radius cases. (6) Finally, a water film (thickness 0.1 mm) is
applied to all the grain surfaces to mimic a water-wet situation as observed in
the experimental porous media (where the contact angles of water to clay-
coated glass bead surface are smaller than 20° as observed in the CCD im-
ages). The phenomenon that all the grain surface is covered by water may not
be observed in the natural water-wet materials because a small fraction of the
solid phase could still be covered by the NAPL phase (e.g., Molnar
et al., 2020). While this small inconsistency would not significantly affect the

change rate of complex conductivity with saturation (especially the values of the saturation exponents), which is
very a true and accurate assumption for the high saturations (S<sub>w</sub> > 0.75 for all simulations) where water dominates
the pore space. Furthermore, the thin water layer here also ensures that the EDL surrounding the grains is un-
affected by the presence of NAPL, thus the upscaled complex conductivity assigned to the grain phase calculated
by Equation 4 remains unchanged with different saturations. This is crucial because it allows us to attribute the
evolution of complex conductivity of the 2D domain entirely to the changes in the current path controlled by pore
fluids distribution, thereby simplifying the physical mechanisms. Meanwhile, the thin water layer around the
grain could prevent the excessive electrical current flow through the grains when NAPL blocks the pore throats,
and ensures a stable numerical solution. This simulation scheme restores a possible way for the pore water
connectivity to change in real conditions and could help us evaluate how saturation exponents change with pore
fluids distribution. We test the REV effect of the simulation scheme by repeated simulations using the case where
the added NAPL ganglia radius equals 0.8 mm. The complex conductivities acquired in three close saturations
(0.823 ± 0.003) exhibit a very small variation (2.4%), indicating that the scale of the simulation domain is close to
or larger than the REV for similar partially saturated case (saturation range from 0.75 to 1).

Figure 8 shows the complex conductivity at 10⁻⁵ Hz versus water saturation acquired by the simulations. To distinguish the bulk component of in-phase conductivity from the surface component in Figure 8b, the surface conductance in Equation 4 is set as zero when doing the specific simulations. At the same saturation, the in-phase and quadrature conductivity of the larger NAPL ganglia series are lower than those with smaller ganglia sizes. It is expected because larger insulating NAPL ganglia size usually corresponds to poor water connectivity as they can more easily block the pore throats and current path (see Figure 7). Similar situations could also occur when air

<figure><img src="image_10.png" /><figcaption>**Figure 8.** The normalized in-phase (exclude surface conductivity) (a) and quadrature (b) conductivity versus water saturation acquired by simulations and experiments. The lighter color for the simulation series indicates larger NAPL ganglia sizes. At the same saturation, the in-phase and quadrature conductivity of the larger NAPL ganglia series are lower than those with smaller ganglia sizes. The black dash lines denote the normalized complex conductivity predicted by the saturation exponents (n = 2.00 and p = 1.17) fitted from the experimental imbibition data in Figure 5.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** The estimated saturation exponents *n* and *p* acquired by simulations versus the setting NAPL ganglia size, the dashed lines are linear fittings. *R*² describes the goodness of fit.</figcaption></figure>

bubbles obstruct the pore space, which also has an effect on the electrical conductivity (Han et al., 2021; Rembert et al., 2023). We also plot the normalized complex conductivity predicted by the saturation exponents obtained from the experimental imbibition data in Figure 5. The experimental fitting line shows similar slopes (equal to the value of saturation exponents) to the simulated data series, particularly when the added ganglia size equal to 0.4–0.5 mm. This similarity indicates the accuracy of the numerical simulations to provide a physical picture consistent with the results of the experiments.

Figure 9 shows the estimated saturation exponents using all 10 data points for each fixed NAPL ganglia size acquired by the simulations (an extended data set with NAPL ganglia radius range from 0.2 to 2 mm is shown). Due to the randomness of the position of NAPL ganglia and uncertainties in the resulting saturations, simulations with the same NAPL size may give different values of saturation exponents. However, both the saturation exponents for the in-phase and quadrature conductivity increase significantly with the NAPL ganglia size. Qiang et al. (2022) argue that the size of NAPL ganglia in pore space could be a representation of the change rate of the pore water connectivity with saturation, which indicates the relationship between the saturation exponents and the dynamic of pore fluids distributions.

# 4. Discussion

## 4.1. Influence of the Change Rate of Pore Water Connectivity With Saturation

The results section as well as previous studies have pointed out that saturation exponents are related to the distribution of the oil/water phase in pore space (Johansson et al., 2019; Schmutz et al., 2012), or more specifically, the change in the connectivity of pore water as conductive fluid phase (Binley & Slater, 2020; Glover, 2017). However, there are still no universal theories to quantitatively explain the variation of the values of saturation exponents with different geomaterials, especially for the saturation exponent *p* for quadrature conductivity.

Here we attempt to derive a universal interpretation of the saturation exponents from the microscopic fluids distributions. We make a key deduction that the factors that affect the saturation exponents such as the size of the NAPL ganglia added in our simulations can be summarized by the “change rate of pore water connectivity with saturation.” For example, when the water saturation is relatively high in our simulations (*S*<sub>w</sub> = 0.75–1), the size of the NAPL ganglia added to the pore space is proportional to the probability that it will cut off the connected water phase (i.e., the change rate of pore water connectivity with saturation) (Qiang et al., 2022). Thus, we can infer further that the “change rate of pore water connectivity with saturation” could be used to infer the saturation exponents.

Then, inspired by Euler characteristics *χ* (−) as a simple topology parameter to reflect the connectivity of a specific phase in pore space (Armstrong et al., 2019), we propose a parameter Δ*C* (−) to quantify the change rate of pore water connectivity with saturation to some extent:

$$ \Delta C = -(\partial\chi_w/\partial S_w) \qquad (9) $$

this equation calculates the partial derivative of Euler characteristic of pore water *χ*<sub>w</sub> with water saturation *S*<sub>w</sub>. Meanwhile, a derivative of *χ*<sub>w</sub> with *S*<sub>w</sub> could relieve the overdependence of *χ*<sub>w</sub> by taking the NAPL ganglia size into concern, and reduce the effects caused by the sensitivity of Euler characteristics itself to pore textures (Li et al., 2015). As an increase in Euler characteristic usually indicates a decrease in pore water connectivity (McClure et al., 2020), a minus is added to redirect the equation. It is worth noting that some other parameters such as coordination number (Bernabé et al., 2010) or fractal dimensions (Chen et al., 2023) could also be used here to describe the pore water connectivity, which may exhibit similar numerical characteristics.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** The estimated saturation exponents versus the change rate of pore water connectivity with saturation $\Delta C$ acquired by the simulations (a) and experiments (b). The dashed lines are guides for the eyes that calculated by the models $n = a \times \Delta C_{eq}^b + n_{emp}$ and $p = a \times \Delta C_{eq}^b + p_{emp}$, where $\Delta C_{eq} = \Delta C - \min(\Delta C) + 1$ denotes the positive equivalent of $\Delta C$, a and b are power parameters to be fit. $n_{emp}$ (equals 2 for simulations and 1.5 for experiments) and $p_{emp}$ (equals 1 for both simulation and experiments) are the initial empirical values of the saturation exponents given to the model. The experimental data exhibit a scattered nature but the saturation exponents still show an upward trend with the increase of $\Delta C$.</figcaption></figure>

Figure 10 plots the estimated $n$ and $p$ versus $\Delta C$ for the simulations and experiments. For the experimental data, $\Delta C$ is calculated by their central finite difference approximations: $-(\chi_{i+4} - \chi_{i-4})/(S_{wi+4} - S_{wi-4})$, where $i$ indicates the evaluated data point. For simulations, $\Delta C$ of each data series is acquired by $-(\chi_1 - \chi_{10})/(S_{w1} - S_{w10})$, where $S_{w1}$ denotes the data points with the highest saturation. Here the Euler characteristic and saturation of the simulated data points are determined directly from the assigned phase distributions in the simulation schemes.

The saturation exponents in Figure 10 follow power-law-type relationships with $\Delta C$ ($n = a \times \Delta C_{eq}^b + n_{emp}$ and $p = a \times \Delta C_{eq}^b + p_{emp}$, where $\Delta C_{eq} = \Delta C - \min(\Delta C) + 1$ is the positive equivalent of $\Delta C$, a and b are power parameters to be fit, $n_{emp} = 2$ for simulations and 1.5 for experiments and $p_{emp} = 1$ for both simulations and experiments denote the initial empirical values of the saturation exponents). The saturation exponents increased dramatically with $\Delta C$ in the power-law-type relationships. This feature is less significant for the experimental results due to the scattered nature of the experimental data, but it is still clear that the saturation exponents increased with $\Delta C$.

Figure 8 shows that the saturation exponents of the high saturation imbibition data are consistent with the simulation schemes with small NAPL ganglia size (0.4–0.5 mm). This is because the imbibition data is obtained in a water-wet micromodel at relatively high saturations, where the corresponding water phase is distributed widely and evenly in the pore space. While the simulation schemes with smaller NAPL ganglia size also correspond to more water-wet scenarios (Qiang et al., 2022; Yang et al., 2007) and exhibit high saturation, the pore fluids evolution pattern of the imbibition experiments and simulation with small NAPL ganglia would be somehow consistent. These simultaneous consistencies in both the pore fluids distribution and saturation exponents indicate the accuracy of the numerical simulations to provide a similar electrical response mechanism to the experiments, even if their absolute values are different. With this consistency, we further attribute the variation of $\Delta C$ and saturation exponents within the experimental and simulated data themselves respectively in Figure 10 to the different evolution patterns of pore fluids distributions. For example, in simulations, $\Delta C$ and saturation exponents increase with the size of the NAPL ganglia added to the pore space, because larger NAPL ganglia added to the pore space are more easily to cut off the connected water phase and reduce the water connectivity. While for experiments, when saturation is relatively high (imbibition data), the isolated pore water phase in experiments is easy to be connected by the add-on saturation as water distributed widely in the pore space, leading to a larger $\Delta C$. Then at low saturations (drainage data), pore water is distributed very dispersively in the form of a limited number of narrow bands along the solid phases. The add-on saturation could not improve the pore water connectivity significantly because they could hardly connect the isolated water phase and form additional current paths. Thus the $\Delta C$ and saturation exponents at this stage are relatively small. However, although the pore fluids evolution pattern in the experiments and simulations, respectively is different, the power-law-type relationships between the saturation exponents and $\Delta C$ for the experiments and simulations are still fixed.
$\Delta C$ could directly describe the final pore fluid distributions determined by different fluid and porous matrix properties, such as wettability (Montaron, 2009; Qin et al., 2023), pore textures (Sun & Torres-Verdín, 2022), viscosity of the fluids (Freytes et al., 2018), and capillary number (Dong et al., 1998; Guo et al., 2022). $\Delta C$ indeed serves as a summary of the factors that potentially control how fluids are distributed in the pore space and how fluids distribution changes with saturation. For example, the wettability of the geomaterials controls the size of the residual NAPL ganglia in pore space (larger in oil-wet media, see Montaron, 2009; Qiang et al., 2022), leading to different $\Delta C$ and saturation exponents for geomaterials with different wettabilities. While we do not aim to figure out the detailed relationships between these factors and pore fluid distribution, as it is intuitively complicated due to numerous controlling factors and possible correlations which is beyond this work's scope. However, the complexity in the controlling factors of pore fluid distribution highlights the value of the parameter $\Delta C$, which could help us to directly describe the state of pore fluids and uncover the physics of saturation exponent.

## 4.2. Comparison With Existing Theories

The proposed power-law-type relationship between $\Delta C$ and saturation exponents qualitatively validates the theory summarized by Binley and Slater (2020) and Qiang et al. (2022). This finding enhances our understanding of the dependence of the saturation exponents from that on the water connectivity (Anderson, 1986; Montaron, 2009) to the change rate of water connectivity with saturation. Nevertheless, these two concepts may be aligned under certain conditions. Porous media with poor water connectivity often exhibit a small $\Delta C$, given that the water connectivity usually increases slowly with saturation in low saturation scenarios.

Interestingly, Yue (2019) derives a relationship between the cementation exponent $m$ and “the change rate of pore complexity with porosity” for saturated porous media, where the concept of “pore complexity” is also associated with the pore tortuous (connectivity). Similarly, Glover (2009) proposes a definition of the cementation exponent $m$ by the rate of change of the connectedness with porosity. These findings support our study well as the saturation exponent is found to be essentially equivalent to the cementation exponent for the unsaturated media (Han et al., 2021). Meanwhile, the cementation exponent is proven to be related to the fractal building process (Chen et al., 2023; Wang & Revil, 2020) of the pore, which is also related to the dynamic process (change) of the pores structure, rather than the pore structure itself.

The relationship between $\Delta C$ and saturation exponents may also offer insights into the percolation behavior of electrical conductivity (Jougnot et al., 2010; Kozlov et al., 2012; Montaron, 2009), which suggests that the porous media do not conduct when saturation is below a certain percolation threshold. We argue that when the saturation is below the percolation threshold, the current electrodes are indeed not connected, and hence the change rate of water “connectivity” between electrodes with saturation is zero. It is noted that the “connectivity” between electrodes here is different from the pore water connectivity we defined based on the Euler characteristic, instead, it emphasized the connectivity between the electrodes to pass the electrical currents. The evolution of electrical conductivity below the percolation threshold could be regarded as a case where both the change rate of water “connectivity” between electrodes with saturation and the saturation exponents equal zero. Meanwhile, the different evolution slope of electrical conductivity with saturation (saturation exponents) near the percolation threshold in log-log scales (Ewing & Hunt, 2006) could also be explained by the difference in the change rates of water “connectivity” between electrodes with saturation.

Han et al. (2009) and Sun and Torres-Verdín (2022) also observed the bend-up behavior of the conductivity curve and reduction of saturation exponents at low saturations as shown in Figure 5. They attribute this deviant behavior of the conductivity curve to the water films on the grain surface, which reserves the percolating path at low saturations. Their theories could be explained by the proposed relationships between $\Delta C$ and saturation exponents. When water is distributed in the form of films along the grains, it usually exhibits a low $\Delta C$ because the add-on saturation could hardly form new percolating paths that are independent of the existing paths. Thus, the pore water connectivity does not increase significantly with saturation (or vice versa), which results in a smaller saturation exponent as in the low saturation cases for our experiments.

## 4.3. Practical Implications for Geophysical Application and Future Directions

The proposed power-law-type relationship between $\Delta C$ and saturation exponents establishes a bridge between the geoelectrical observations and the microscopic properties of the subsurface porous media. This unique relationship could help to refine the physical interpretation of the saturation exponents and enhance the interpretations
of the geoelectrical data under partially saturated conditions. Notably, this work first links the saturation exponent *p* for the quadrature conductivity to the pore fluid distribution. Which has rarely been studied in previous geophysical works that focus solely on the bulk components of in-phase conductivity (Han et al., 2021). The new findings on quadrature conductivity could provide an alternative path for the application of geoelectrical methods and data interpretations in partially saturated clay-rich areas. Because clay-rich areas usually exhibit significant IP response and surface conduction that would disturb the interpretation of in-phase conductivity data alone.

Moreover, given some information on the porous matrix and fluid properties, the proposed relationship between ΔC and saturation exponents can assist in refining the saturation exponent estimations in the absence of laboratory measurements. For example, when estimating subsurface hydrocarbon distribution using Archie's law, a higher saturation exponent for oil-wet geomaterials can be directly assumed for oil-wet media. This is because oil-wet media are known to have higher ΔC due to the susceptible internal pore water connectivity (Montaron, 2009; Qiang et al., 2022). This larger expectation of the saturation exponents may significantly improve the accuracy of the hydrocarbon volume estimation than simply assigning them with empirical values. Conversely, knowledge of the saturation exponents for specific geomaterial, can aid in predicting the properties of porous matrix and fluid. For instance, the value of saturation exponents can help to estimate the amount or size of the air bubbles generated in the pore space by reactive transport processes such as calcite dissolution (Rembert et al., 2023). Thereby improving our understanding of the status of subsurface reactions.

One limitation of this work is that the simulation work does not consider low saturation cases. This is because the simulation schemes require the non-conductive NAPL phase to have a very small effective conductivity (10⁻⁶ S/m) for the convergence of the numerical solutions. Consequently, this leads to unrealistically complex conductivity at low saturation and fails the fitting of saturation exponents. Meanwhile, the “upscaling” numerical method used in this study, which transforms the EDL polarization on the grain surface to the entire grain, cannot quantitatively account for the detailed change in EDL structures. Therefore, a more generalized and practical mathematical framework that incorporates the Poisson-Nernst-Planck equations (Bücker et al., 2019) is needed for future studies. With which we could fully model the EDL and the low-saturation behaviors of the complex conductivity with different pore fluid distributions.

Future studies should also aim to establish the linkage between saturation exponent and more generalized physical parameters. Although the acquisition of parameter ΔC is intuitive, it requires extremely high-resolution images of the porous media, which poses great challenges in practical applications. Fractal dimensions, recognized as a powerful tool in characterizing microstructures (Chen et al., 2023; Katz & Thompson, 1985), are more widely used than Euler characteristics in describing porous morphology. Its relationship to electrical conductivity has been extensively studied (Thanh et al., 2020; Wang & Revil, 2020; Wei et al., 2015). Most importantly, the fractal building process inherently corresponds to changes in the pore structure (Wang & Liu, 2017), which is aligned with the current finding that saturation exponents are related to the change rate of pore water connectivity. Thus, future studies could focus on the dependence of saturation exponents on the fractal properties of the porous media for a more generalized relationship. Furthermore, the internal physic of the empirical relationship p = n − 1 between the saturation exponents (Revil, 2013; Vinegar & Waxman, 1984) should also be explored in future studies, which may help to elucidate the difference in the physical interpretation of the saturation exponents.

**5. Conclusions**

We have developed a visual milli-fluidic micromodel to investigate the saturation exponents of porous media. The micromodel allows us to acquire the microscopic distributions of pore fluids while measuring complex conductivity. Pore network complex conductivity simulations are also conducted for the uncontrolled ideal fluids distribution.

The experimental and simulated results show a variation of saturation exponents value with different distributions of pore fluids in the same porous medium. We then further introduced a parameter ΔC defined as the change rate of pore water connectivity with saturation to capture and quantify the fluid and porous matrix properties associated with saturation exponents. ΔC is calculated by the partial derivative of pore water Euler characteristics with saturation. Notably, our results reveal that the saturation exponents for both in-phase and quadrature conductivity vary with the change rate of pore water connectivity with saturation ΔC according to a power-law-type relationship. These findings establish a more physical and quantitative relationship between the saturation exponents for induced polarization and the textural properties of the porous material and pore fluid properties. For future
studies, it is recommended to focus on the physics of the electrical double layer under different wetting conditions, as well as solving the Poisson-Nernst-Planck equations. The dependence of saturation exponents on the pore fluid fractal properties may also be worth studying for more generalized relationships.

## Data Availability Statement

The experimental raw data and numerical models used for this manuscript are available on Zenodo via https://doi.org/10.5281/zenodo.12819973 (Qiang, 2023).

## Acknowledgments

The work was supported by the National Key Research and Development Program of China (2022YFC3703101) and the National Natural Science Foundation of China (42272276, 41977157). Xiaoqing Shi was partly supported by the Fundamental Research Funds for the Central Universities (020614380217 and 020614380181). Siyuan Qiang thanks the China Scholarship Council (202306190219) and the Postgraduate Research and Practice Innovation Program of Jiangsu Province (KYCX24_0209) for their support. We appreciate Professor Lanchang Xing for the discussion on the numerical methods. We are grateful to Editor Yves Bernabé, the anonymous Associate Editor, and two anonymous reviewers for their insightful comments and suggestions, which significantly improve the quality of this work.

## References

Ake, C., Wiles, M., Huebner, H., McDonald, T., Cosgriff, D., Richardson, M., et al. (2003). Porous organoclay composite for the sorption of polycyclic aromatic hydrocarbons and pentachlorophenol from groundwater. *Chemosphere*, *51*(9), 835-844. https://doi.org/10.1016/S0045-6535(03)00040-7

Aksu, I., Bazilevskaya, E., & Karpyn, Z. (2015). Swelling of clay minerals in unconsolidated porous media and its impact on permeability. *GeoResJ*, *7*, 1-13. https://doi.org/10.1016/j.grj.2015.02.003

Anderson, W. G. (1986). Wettability literature survey-part 3: The effects of wettability on the electrical properties of porous media. *Journal of Petroleum Technology*, *38*(12), 1371-1378. https://doi.org/10.2118/13934-PA

Archie, G. E. (1942). The electrical resistivity log as an aid in determining some reservoir characteristics. *Transactions of the AIME*, *146*(01), 54-62. https://doi.org/10.2118/942054-G

Armstrong, R. T., McClure, J. E., Robins, V., Liu, Z., Arns, C. H., Schlüter, S., & Berg, S. (2019). Porous media characterization using Minkowski functionals: Theories, applications and future directions. *Transport in Porous Media*, *130*(1), 305-335. https://doi.org/10.1007/s11242-018-1201-4

Bernabé, Y., Li, M., & Maineult, A. (2010). Permeability and pore connectivity: A new model based on network simulations. *Journal of Geophysical Research*, *115*(B10), B10203. https://doi.org/10.1029/2010JB007444

Binley, A., Hubbard, S. S., Huisman, J. A., Revil, A., Robinson, D. A., Singha, K., & Slater, L. D. (2015). The emergence of hydrogeophysics for improved understanding of subsurface processes over multiple scales. *Water Resources Research*, *51*(6), 3837-3866. https://doi.org/10.1002/2015WR017016

Binley, A., & Slater, L. (2020). *Resistivity and induced polarization: Theory and applications to the near-surface Earth*. Cambridge University Press. https://doi.org/10.1017/9781108685955

Bücker, M., Flores Orozco, A., Undorf, S., & Kemna, A. (2019). On the role of Stern-and diffuse-layer polarization mechanisms in porous media. *Journal of Geophysical Research: Solid Earth*, *124*(6), 5656-5677. https://doi.org/10.1029/2019JB017679

Bussian, A. (1983). Electrical conductance in a porous medium. *Geophysics*, *48*(9), 1258-1268. https://doi.org/10.1190/1.1441549

Chen, X., Thanh, L. D., Luo, C., Tahmasebi, P., & Cai, J. (2023). Dependence of electrical conduction on pore structure in reservoir rocks from the Beibuwan and Pearl River Mouth Basins: A theoretical and experimental study. *Geophysics*, *88*(2), MR35-MR53. https://doi.org/10.1190/GEO2021-0682.1

Deng, Y., Shi, X., Revil, A., Wu, J., & Ghorbani, A. (2018). Complex conductivity of oil-contaminated clayey soils. *Journal of Hydrology*, *561*, 930-942. https://doi.org/10.1016/j.jhydrol.2018.04.055

Dong, M., Dullien, F. A. L., & Zhou, J. (1998). Characterization of waterflood saturation profile histories by the 'complete' capillary number. *Transport in Porous Media*, *31*(2), 213-237. https://doi.org/10.1023/a:1006565621860

El Alam, J., Revil, A., & Dick, P. (2023). Influence of the water content on the complex conductivity of bentonite. *Engineering Geology*, *322*, 107183. https://doi.org/10.1016/j.enggeo.2023.107183

Ewing, R., & Hunt, A. (2006). Dependence of the electrical conductivity on saturation in real porous media. *Vadose Zone Journal*, *5*(2), 731-741. https://doi.org/10.2136/vzj2005.0107

Fernandez Visentini, A., de Anna, P., Jougnot, D., Le Borgne, T., Méheust, Y., & Linde, N. (2023). Electrical signatures of diffusion-limited mixing: Insights from a milli-fluidic tracer experiment. *Transport in Porous Media*, *146*(1), 435-461. https://doi.org/10.1007/s11242-021-01607-0

Fiandaca, G., Ramm, J., Binley, A., Gazoty, A., Christiansen, A. V., & Auken, E. (2013). Resolving spectral information from time domain induced polarization data through 2-D inversion. *Geophysical Journal International*, *192*(2), 631-646. https://doi.org/10.1093/gji/ggs060

Freytes, V. M., Rosen, M., & D'Onofrio, A. (2018). Capillary film and breakup mechanism in the squeezing to dripping transition regime at the mesoscale between micro and milli-fluidics. *Chaos: An Interdisciplinary Journal of Nonlinear Science*, *28*(10). https://doi.org/10.1063/1.5033451

Glover, P. (2009). What is the cementation exponent? A new interpretation. *The Leading Edge*, *28*(1), 82-85. https://doi.org/10.1190/1.3064150

Glover, P. (2017). A new theoretical interpretation of Archie's saturation exponent. *Solid Earth*, *8*(4), 805-816. https://doi.org/10.5194/se-8-805-2017

Gregory, A. P., & Clarke, R. (2012). Tables of the complex permittivity of dielectric reference liquids at frequencies up to 5 GHz.

Guo, H., Song, K., & Hilfer, R. (2022). A brief review of capillary number and its use in capillary desaturation curves. *Transport in Porous Media*, *144*(1), 3-31. https://doi.org/10.1007/s11242-021-01743-7

Han, M., Youssef, S., Rosenberg, E., Fleury, M., & Levitz, P. (2009). Deviation from Archie's law in partially saturated porous media: Wetting film versus disconnectedness of the conducting phase. *Physical Review E*, *79*(3), 031127. https://doi.org/10.1103/PhysRevE.79.031127

Han, T., Yan, H., & Fu, L. Y. (2021). A quantitative interpretation of the saturation exponent in Archie's equations. *Petroleum Science*, *18*(2), 444-449. https://doi.org/10.1007/s12182-021-00547-0

Higham, D. J., & Higham, N. J. (2016). *MATLAB guide* (3rd ed.). Society for Industrial and Applied Mathematics. https://doi.org/10.1137/1.9781611974669

Izumoto, S., Huisman, J. A., Zimmermann, E., Heyman, J., Gomez, F., Tabuteau, H., et al. (2022). Pore-scale mechanisms for spectral induced polarization of calcite precipitation inferred from geo-electrical millifluidics. *Environmental Science & Technology*, *56*(8), 4998-5008. https://doi.org/10.1021/acs.est.1c07742

Jerez, J., Flury, M., Shang, J., & Deng, Y. (2006). Coating of silica sand with aluminosilicate clay. *Journal of Colloid and Interface Science*, *294*(1), 155-164. https://doi.org/10.1016/j.jcis.2005.07.017
Johansson, S., Fiandaca, G., & Dahlin, T. (2015). Influence of non-aqueous phase liquid configuration on induced polarization parameters: Conceptual models applied to a time-domain field case study. *Journal of Applied Geophysics*, **123**, 295–309. https://doi.org/10.1016/j.jappgeo.2015.08.010

Johansson, S., Rossi, M., Hall, S. A., Sparrenbom, C., Hagerberg, D., Tudisco, E., et al. (2019). Combining spectral induced polarization with x-ray tomography to investigate the importance of DNAPL geometry in sand samples. *Geophysics*, **84**(3), E173–E188. https://doi.org/10.1190/geo2018-0108.1

Jougnot, D., Ghorbani, A., Revil, A., Leroy, P., & Cosenza, P. (2010). Spectral induced polarization of partially saturated clay-rocks: A mechanistic approach. *Geophysical Journal International*, **180**(1), 210–224. https://doi.org/10.1111/j.1365-246X.2009.04426.x

Jougnot, D., Jiménez-Martínez, J., Legendre, R., Le Borgne, T., Méheust, Y., & Linde, N. (2018). Impact of small-scale saline tracer heterogeneity on electrical resistivity monitoring in fully and partially saturated porous media: Insights from geoelectrical milli-fluidic experiments. *Advances in Water Resources*, **113**, 295–309. https://doi.org/10.1016/j.adwwatres.2018.01.014

Katz, A. J., & Thompson, A. (1985). Fractal sandstone pores: Implications for conductivity and pore formation. *Physical Review Letters*, **54**(12), 1325–1328. https://doi.org/10.1103/PhysRevLett.54.1325

Kemna, A., Binley, A., Cassiani, G., Niederleithinger, E., Revil, A., Slater, L., et al. (2012). An overview of the spectral induced polarization method for near-surface applications. *Near Surface Geophysics*, **10**(6), 453–468. https://doi.org/10.3997/1873-0604.2012027

Kozlov, B., Schneider, M., Montaron, B., Lagues, M., & Tabeling, P. (2012). Archie's law in microsystems. *Transport in Porous Media*, **95**, 1–20. https://doi.org/10.1007/s11242-012-0029-6

Lebedev, M., Toms-Stewart, J., Clennell, B., Pervukhina, M., Shulakova, V., Paterson, L., et al. (2009). Direct laboratory observation of patchy saturation and its effects on ultrasonic velocities. *The Leading Edge*, **28**(1), 24–27. https://doi.org/10.1190/1.3064142

Leroy, P., & Revil, A. (2009). A mechanistic model for the spectral induced polarization of clay materials. *Journal of Geophysical Research*, **114**(B10), B10202. https://doi.org/10.1029/2008JB006114

Leroy, P., Revil, A., Kemna, A., Cosenza, P., & Ghorbani, A. (2008). Complex conductivity of water-saturated packs of glass beads. *Journal of Colloid and Interface Science*, **321**(1), 103–117. https://doi.org/10.1016/j.jcis.2007.12.031

Li, M., Tang, Y., Bernabé, Y., Zhao, J., Li, X., Bai, X., & Zhang, L. (2015). Pore connectivity, electrical conductivity, and partial water saturation: Network simulations. *Journal of Geophysical Research: Solid Earth*, **120**(6), 4055–4068. https://doi.org/10.1002/2014JB011799

Li, S., Liu, M., Hanaor, D., & Gan, Y. (2018). Dynamics of viscous entrapped saturated zones in partially wetted porous media. *Transport in Porous Media*, **125**(2), 193–210. https://doi.org/10.1007/s11242-018-1113-3

Lim, J. S. (1990). *Two-dimensional signal and image processing*. Prentice-Hall.

Maineult, A., Jougnot, D., & Revil, A. (2018). Variations of petrophysical properties and spectral induced polarization in response to drainage and imbibition: A study on a correlated random tube network. *Geophysical Journal International*, **212**(2), 1398–1411. https://doi.org/10.1093/gji/ggx474

Maineult, A., Revil, A., Camerlynck, C., Florsch, N., & Titov, K. (2017). Upscaling of spectral induced polarization response using random tube networks. *Geophysical Journal International*, **209**(2), 948–960. https://doi.org/10.1093/gji/ggx066

McBride, J., Brockman, F., Szecsody, J., & Streile, G. (1992). Kinetics of quinoline biodegradation, sorption and desorption in a clay-coated model soil containing a quinoline-degrading bacterium. *Journal of Contaminant Hydrology*, **9**(1–2), 133–154. https://doi.org/10.1016/0169-7722(92)90055-J

McClure, J. E., Ramstad, T., Li, Z., Armstrong, R. T., & Berg, S. (2020). Modeling geometric state for fluids in porous media: Evolution of the euler characteristic. *Transport in Porous Media*, **133**(2), 229–250. https://doi.org/10.1007/s11242-020-01420-1

Mecke, K., & Arns, C. (2005). Fluids in porous media: A morphometric approach. *Journal of Physics: Condensed Matter*, **17**(9), S503–S534. https://doi.org/10.1088/0953-8984/17/9/014

Mendieta, A., Jougnot, D., Leroy, P., & Maineult, A. (2021). Spectral induced polarization characterization of non-consolidated clays for varying salinities—An experimental study. *Journal of Geophysical Research: Solid Earth*, **126**(4), e2020JB021125. https://doi.org/10.1029/2020JB021125

Miller, S. E., & Low, P. (1990). Characterization of the electrical double layer of montmorillonite. *Langmuir*, **6**(3), 572–578. https://doi.org/10.1021/la00093a010

Moebius, F., & Or, D. (2014). Pore scale dynamics underlying the motion of drainage fronts in porous media. *Water Resources Research*, **50**(11), 8441–8457. https://doi.org/10.1002/2014WR015916

Molnar, I. L., Gerhard, J. I., Willson, C. S., & O'Carroll, D. M. (2020). Wettability effects on primary drainage mechanisms and NAPL distribution: A pore-scale study. *Water Resources Research*, **56**(1), e2019WR025381. https://doi.org/10.1029/2019WR025381

Montaron, B. (2009). Connectivity theory-a new approach to modeling non-Archie rocks. *Petrophysics-The SPWLA Journal of Formation Evaluation and Reservoir Description*, **50**(02).

Moura, M., Måløy, K. J., Flekkøy, E. G., & Toussaint, R. (2020). Intermittent dynamics of slow drainage experiments in porous media: Characterization under different boundary conditions. *Frontiers in Physics*, **7**, 217. https://doi.org/10.3389/fphy.2019.00217

Mustofa, M. B., Fauzi, U., Latief, F. D. E., & Warsa, W. (2022). Experimental and modeling of electrical resistivity changes due to micro-spatial distribution of fluid for unconsolidated sand. *Journal of Petroleum Science and Engineering*, **208**, 109472. https://doi.org/10.1016/j petrol.2021.109472

Niu, Q., & Revil, A. (2016). Connecting complex conductivity spectra to mercury porosimetry of sedimentary rocks. *Geophysics*, **81**(1), E17–E32. https://doi.org/10.1190/GEO2015-0072.1

Niu, Q., & Zhang, C. (2017). Pore-scale modelling of complex conductivity of saturated granular materials. *Near Surface Geophysics*, **15**(6), 593–602. https://doi.org/10.3997/1873-0604.2017055

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. *Journal of Geophysical Research: Solid Earth*, **125**(10), e2020JB020515. https://doi.org/10.1029/2020JB020515

Ostu, N. (1979). A threshold selection method from gray-level histograms. *IEEE Trans SMC*, **9**(1), 62–66. https://doi.org/10.1109/TSMC.1979.4310076

Panwar, N., Revil, A., Sharma, R., Schmutz, M., Duvillard, P., Garcia, B., et al. (2021). Induced polarization of carbonates. *Journal of Geophysical Research: Solid Earth*, **126**(6), e2021JB022029. https://doi.org/10.1029/2021JB022029

Qiang, S. (2023). Experimental raw data and numerical model [Dataset]. Zenodo. https://doi.org/10.5281/zenodo.12819973

Qiang, S., Shi, X., Revil, A., Kang, X., Liu, Y., & Wu, J. (2022). Residual NAPL morphology effects on electrical resistivity: Insights from micromodel Fiandaca ement experiments and pore network simulations. *Water Resources Research*, **58**(12), e2022WR033233. https://doi.org/10.1029/2022WR033233
Qin, X., Xia, Y., Qiao, J., Chen, J., Zeng, J., & Cai, J. (2023). Modeling of multiphase flow in low permeability porous media: Effect of wettability and pore structure properties. *Journal of Rock Mechanics and Geotechnical Engineering*, **16**(4), 1127–1139. https://doi.org/10.1016/j.jrmge.2023.06.007

Rabaute, A., Revil, A., & Brosse, E. (2003). In situ mineralogy and permeability logs from downhole measurements: Application to a case study in chlorite-coated sandstones. *Journal of Geophysical Research*, **108**(B9), 2414. https://doi.org/10.1029/2002JB002178

Rembert, F., Stolz, A., Soulaine, C., & Roman, S. (2023). A microfluidic chip for geoelectrical monitoring of critical zone processes. *Lab on a Chip*, **23**(15), 3433–3442. https://doi.org/10.1039/d3lc00377a

Revil, A. (2012). Spectral induced polarization of shaly sands: Influence of the electrical double layer. *Water Resources Research*, **48**(2), W02517. https://doi.org/10.1029/2011WR011260

Revil, A. (2013). Effective conductivity and permittivity of unsaturated porous materials in the frequency range 1 MHz–1 GHz. *Water Resources Research*, **49**(1), 306–327. https://doi.org/10.1029/2012WR012700

Revil, A., Binley, A., Mejus, L., & Kessouri, P. (2015). Predicting permeability from the characteristic relaxation time and intrinsic formation factor of complex conductivity spectra. *Water Resources Research*, **51**(8), 6672–6700. https://doi.org/10.1002/2015WR017074

Revil, A., Cathles, L. M., Losh, S., & Nunn, J. (1998). Electrical conductivity in shaly sands with geophysical applications. *Journal of Geophysical Research*, **103**(B10), 23925–23936. https://doi.org/10.1029/98JB02125

Revil, A., Coperey, A., Deng, Y., Cerepi, A., & Seleznev, N. (2018). Complex conductivity of tight sandstones. *Geophysics*, **83**(2), E55–E74. https://doi.org/10.1190/geo2017-0096.1

Revil, A., Coperey, A., Shao, Z., Florsch, N., Fabricius, I. L., Deng, Y., et al. (2017). Complex conductivity of soils. *Water Resources Research*, **53**(8), 7121–7147. https://doi.org/10.1002/2017WR020655

Revil, A., & Florsch, N. (2010). Determination of permeability from spectral induced polarization in granular media. *Geophysical Journal International*, **181**(3), 1480–1498. https://doi.org/10.1111/j.1365-246X.2010.04573.x

Revil, A., Florsch, N., & Camerlynck, C. (2014). Spectral induced polarization porosimetry. *Geophysical Journal International*, **198**(2), 1016–1033. https://doi.org/10.1093/gji/ggu180

Revil, A., Ghorbani, A., Jougnot, D., & Yven, B. (2023). Induced polarization of clay-rich materials—Part 1: The effect of desiccation. *Geophysics*, **88**(4), MR195-MR210. https://doi.org/10.1190/GEO2022-0510.1

Revil, A., Kessouri, P., & Torres-Verdín, C. (2014). Electrical conductivity, induced polarization, and permeability of the Fontainebleau sandstone. *Geophysics*, **79**(5), D301–D318. https://doi.org/10.1190/GEO2014-0036.1

Revil, A., Koch, K., & Holliger, K. (2012). Is it the grain size or the characteristic pore size that controls the induced polarization relaxation time of clean sands and sandstones? *Water Resources Research*, **48**(5), W05602. https://doi.org/10.1029/2011WR011561

Revil, A., Schmutz, M., & Batzle, M. L. (2011). Influence of oil wettability upon spectral induced polarization of oil-bearing sands. *Geophysics*, **76**(5), A31–A36. https://doi.org/10.1190/GEO2011-0006.1

Revil, A., & Skold, M. (2011). Salinity dependence of spectral induced polarization in sands and sandstones. *Geophysical Journal International*, **187**(2), 813–824. https://doi.org/10.1111/j.1365-246X.2011.05181.x

Revil, A., Skold, M., Hubbard, S. S., Wu, Y., Watson, D. B., & Karaoulis, M. (2013). Petrophysical properties of saprolites from the Oak Ridge integrated field research challenge site, Tennessee. *Geophysics*, **78**(1), D21–D40. https://doi.org/10.1190/geo2012-0176.1

Schmutz, M., Blondel, A., & Revil, A. (2012). Saturation dependence of the quadrature conductivity of oil-bearing sands. *Geophysical Research Letters*, **39**(3), L03402. https://doi.org/10.1029/2011GL050474

Schmutz, M., Revil, A., Vaudelet, P., Batzle, M., Viñao, P. F., & Werkema, D. (2010). Influence of oil saturation upon spectral induced polarization of oil-bearing sands. *Geophysical Journal International*, **183**(1), 211–224. https://doi.org/10.1111/j.1365-246X.2010.04751.x

Scholz, C., Wirner, F., Götz, J., Rüde, U., Schröder-Turk, G. E., Mecke, K., & Bechinger, C. (2012). Permeability of porous materials determined from the Euler characteristic. *Physical Review Letters*, **109**(26), 264504. https://doi.org/10.1103/PhysRevLett.109.264504

Sun, Z., Mehmani, A., & Torres-Verdín, C. (2021). Pore-scale investigation of the electrical resistivity of saturated porous media: Flow patterns and porosity efficiency. *Journal of Geophysical Research: Solid Earth*, **126**(12), e2021JB022608. https://doi.org/10.1029/2021JB022608

Sun, Z., & Torres-Verdín, C. (2022). The role of pore-shape and pore-space heterogeneity in non-Archie behavior of resistivity index curves. *Journal of Geophysical Research: Solid Earth*, **127**(9), e2022JB024792. https://doi.org/10.1029/2022JB024792

Tarasov, A., & Titov, K. (2013). On the use of the Cole–Cole equations in spectral induced polarization. *Geophysical Journal International*, **195**(1), 352–356. https://doi.org/10.1093/gji/ggt251

Thanh, L. D., Jougnot, D., Van Do, P., Tuyen, V. P., Ca, N. X., & Hien, N. T. (2020). A physically based model for the electrical conductivity of partially saturated porous media. *Geophysical Journal International*, **223**(2), 993–1006. https://doi.org/10.1093/gji/ggaa307

Tian, H., Liu, L., Zhu, L., Ge, X., Ding, P., & Cai, J. (2023). Relationship between normalized permeability and resistivity index in hydrate-bearing sediments: Fractal model and numerical simulation. *Geophysical Journal International*, **234**(1), 684–698. https://doi.org/10.1093/gji/ggad090

Toms, J., Müller, T., Ciz, R., & Gurevich, B. (2006). Comparative review of theoretical models for elastic wave attenuation and dispersion in partially saturated rocks. *Soil Dynamics and Earthquake Engineering*, **26**(6–7), 548–565. https://doi.org/10.1016/j.soildyn.2006.01.008

Vinegar, H., & Waxman, M. (1984). Induced polarization of shaly sands. *Geophysics*, **49**(8), 1267–1287. https://doi.org/10.1190/1.1441755

Wang, H., & Liu, T. (2017). Derivation of Archie's law based on a fractal pore volume. *Geophysical Journal International*, **209**(3), 1403–1407. https://doi.org/10.1093/gji/ggx095

Wang, H., & Revil, A. (2020). Surface conduction model for fractal porous media. *Geophysical Research Letters*, **47**(10), e2020GL087553. https://doi.org/10.1029/2020GL087553

Waxman, M. H., & Smits, L. (1968). Electrical conductivities in oil-bearing shaly sands. *Society of Petroleum Engineers Journal*, **8**(02), 107–122. https://doi.org/10.2118/1863-A

Wei, W., Cai, J., Hu, X., & Han, Q. (2015). An electrical conductivity model for fractal porous media. *Geophysical Research Letters*, **42**(12), 4833–4840. https://doi.org/10.1002/2015GL064460

Weller, A., & Slater, L. (2012). Salinity dependence of complex conductivity of unconsolidated and consolidated materials: Comparisons with electrical double layer models. *Geophysics*, **77**(5), D185–D198. https://doi.org/10.1190/GEO2012-0030.1

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

Woodruff, W., & Revil, A. (2011). CEC-normalized clay-water sorption isotherm. *Water Resources Research*, **47**(11), W11502. https://doi.org/10.1029/2011WR010919

Xing, L., Qi, S., Xu, Y., Wang, B., Lao, L., Wei, W., et al. (2021). Numerical study on complex conductivity characteristics of hydrate-bearing porous media. *Journal of Natural Gas Science and Engineering*, **95**, 104145. https://doi.org/10.1016/j.jngse.2021.104145
Xing, L., Zhang, H., Wang, S., Wang, B., Lao, L., Wei, W., et al. (2023). Pore-scale modelling on complex-conductivity responses of hydrate-bearing clayey sediments: Implications for evaluating hydrate saturation and clay content. *Geoenergy Science and Engineering*, 221, 211356.
https://doi.org/10.1016/j.geoen.2022.211356

Yang, Z., Shao, C., Zhou, G., & Qiu, C. (2007). Modeling non-aqueous phase liquid displacement process. *Petroleum Science*, 4(3), 40–46. https://doi.org/10.1007/s12182-007-0007-8

Yue, W. (2019). Pore-scale explanation of the Archie's cementation exponent: Microstructure, electrical anisotropy, and numerical experiments. *Geophysical Research Letters*, 46(11), 5799–5807. https://doi.org/10.1029/2019GL082585

Zarif, F., Kessouri, P., & Slater, L. (2017). Recommendations for field-scale induced polarization (IP) data acquisition and interpretation. *Journal of Environmental & Engineering Geophysics*, 22(4), 395–410. https://doi.org/10.2113/JEEG22.4.395

Zhang, Z., Liu, L., Li, C., Cai, J., Ning, F., Meng, Q., & Liu, C. (2021). Fractal analyses on saturation exponent in Archie's law for electrical properties of hydrate-bearing porous media. *Journal of Petroleum Science and Engineering*, 196, 107642. https://doi.org/10.1016/j.petrol.2020.107642

Zhu, L., Wu, S., Zhou, X., & Cai, J. (2023). Saturation evaluation for fine-grained sediments. *Geoscience Frontiers*, 14(4), 101540. https://doi.org/10.1016/j.gsf.2023.101540

Zimmermann, E., Kemna, A., Berwix, J., Glaas, W., Münch, H., & Huisman, J. (2008). A high-accuracy impedance spectrometer for measuring sediments with low polarizability. *Measurement Science and Technology*, 19(10), 105603. https://doi.org/10.1088/0957-0233/19/10/105603

Zonge, K., Sauck, W., & Sumner, J. (1972). Comparison of time, frequency, and phase measurements in induced polarization. *Geophysical Prospecting*, 20(3), 626–648. https://doi.org/10.1111/j.1365-2478.1972.tb00658.x