2025 Volume 66 Issue 9 Pages 1258-1265
Understanding rock resistivity is essential for interpreting subsurface resistivity structures obtained from electromagnetic surveys. We can estimate porosity or water content from resistivity using rock physics models; however, predicted results are strongly dependent on the model with the assumption of an internal pore microstructure. In this study, we evaluated the relationship between the resistivity and pore microstructure of clay-free sandstones based on digital rock physics while focusing on tortuosity. Simulation results demonstrated an increase in resistivity and its anisotropy with decreasing porosity. Tortuosity values calculated from local electric current further explain the evolution of resistivity, which suggests that smaller pore volumes (i.e., porosity) prevent pore connectivity and enhance tortuosity, resulting in higher resistivity and anisotropy. The resistivity of high-porosity data can be fitted using Archie’s equation with empirical parameters, whereas the resistivity of low-porosity data cannot. This illustrates the difficulty in applying Archie’s equation to low-porosity rocks; however, the equivalent channel model reproduced the resistivity over a wide porosity range using the calculated tortuosity. Our results confirm that tortuosity is a key factor for determining electrical properties.
This Paper was Originally Published in Japanese in J. Soc. Mater. Sci., Japan 73 (2024) 232–239. The abstract and all figures are slightly modified.
Estimating the amount of subsurface water is crucial for understanding the structures of seismogenic zones and geothermovolcanic regions. Electromagnetic prospecting is a commonly used technique for estimating the amount of subsurface water [1]. The water distribution in high-resistivity bedrock can be estimated using resistivity because of the sensitivity of its resistivity structure to the presence of water. In geothermal and volcanic regions, hydrothermal systems closely related to geothermal resources and phreatic eruptions are often discussed based on resistivity structures [2]. In seismogenic zones, subsurface resistivity structures are widely used to determine the contribution of fluids to earthquake occurrence [3]. A rock physics model that expresses the relationship between the resistivity of water-bearing rocks and water content in an equation is necessary for quantitatively estimating the amount of subsurface water, not just its presence or absence, from resistivity. Rock physics models include theoretical models that describe analytical solutions by assuming the geometric shape of pores and empirical laws constructed to explain experimental values.
Empirical laws are preferred to describe rock resistivity because simple geometric shapes assumed in theoretical models cannot explain the experimental values and visualizing the actual internal pore structure is difficult. In recent years, the widespread use of microfocus X-ray computed tomography (CT) scanners with micrometer resolution made it possible to visualize pore structures of rocks at a high resolution. With the development of computer performance and computational techniques, digital rock physics, which can calculate macroscopic physical properties (such as elastic wave velocity and permeability) from digitized rocks (digital rocks) based on CT images, has been established [4–6]. Meanwhile, there are a few studies to estimate rock resistivity using digital rock physics. The few existing studies focus on evaluating its performance as an alternative to laboratory experiments [7], and no research has focused on the effect of the internal pore structure of porous rocks on rock resistivity.
Recent studies that utilize machine learning indicate that tortuosity is a dominant factor for determining resistivity [8]. In addition, the relationship between resistivity and porosity in rocks containing fractures changes when tortuosity exceeds a certain threshold [9]. Given this importance of tortuosity in controlling electrical properties, there is a need to systematically investigate how the tortuosity of electric current paths through pore networks influences rock resistivity. Therefore, this study calculated the electric current distribution and rock resistivity for digitized sandstones with various porosities using digital rock physics, and by focusing on the tortuosity of electric current paths, investigated the relationship between pore structure and resistivity in porous rocks. We did not consider the effects of clay minerals or other conductive minerals that could complicate the relationship between tortuosity and resistivity, and focused on sandstones that do not contain either.
Representative geometric models, which include spherical, tube, and film models [10], are respectively expressed as
| \begin{equation} \rho_{\text{eff}} \simeq \rho_{\text{R}}, \end{equation} | (1) |
| \begin{equation} \rho_{\text{eff}} = \frac{3\rho_{\text{W}}}{f\phi}, \end{equation} | (2) |
| \begin{equation} \rho_{\text{eff}} = \frac{3\rho_{\text{W}}}{2f\phi}, \end{equation} | (3) |
where ρeff, ρR, ρW, ϕ, and f represent the effective resistivity of a water-bearing rock, resistivity of a dry rock, resistivity of a liquid in the pores, porosity, and volume fraction of pores that belong to a connected network, respectively. These values are derived by approximating the pore structure of the rock as spheres, tubes, or films and assuming that ρW ≪ ρR and ϕ is small [11]. The Hashin–Shtrikman model, which assumes a two-phase spherical structure with inner and outer shells, corresponds to the geometric shapes of completely connected and isolated states [12]. The effective resistivity values of these structures are respectively expressed as
| \begin{equation} \rho_{\text{eff}} = \left(\frac{\rho_{\text{W}}^{-1}(3\rho_{\text{R}}^{-1} + 2\phi (\rho_{\text{W}}^{-1} - \rho_{\text{R}}^{-1}))}{3\rho_{\text{W}}^{-1} - \phi (\rho_{\text{W}}^{-1} - \rho_{\text{R}}^{-1})}\right)^{-1}, \end{equation} | (4) |
| \begin{equation} \rho_{\text{eff}} = \left(\frac{\rho_{\text{R}}^{-1}(\rho_{\text{W}}^{-1} + 2\rho_{\text{R}}^{-1} + 2\phi (\rho_{\text{W}}^{-1} - \rho_{\text{R}}^{-1}))}{\rho_{\text{W}}^{-1} + 2\rho_{\text{R}}^{-1} - \phi (\rho_{\text{W}}^{-1} - \rho_{\text{R}}^{-1})} \right)^{-1}. \end{equation} | (5) |
Other models include a percolation model [13], which expresses pores as a connection of randomly placed tubes, and an equivalent channel model [14], which considers the tortuosity τ2 of the flow path. These models are respectively expressed as
| \begin{equation} \rho_{\text{eff}} = \frac{4\rho_{\text{W}}}{f\phi}. \end{equation} | (6) |
| \begin{equation} \rho_{\text{eff}} = \frac{\tau^{2}}{\phi}\rho_{\text{W}}. \end{equation} | (7) |
Tortuosity represents the bending of a path; in porous media, it indicates the complexity of the pore path. When considering a certain path, tortuosity τ2 is defined as the ratio of the actual path l to the shortest path x as
| \begin{equation} \tau^{2} = \left(\frac{l}{x} \right)^{2}. \end{equation} | (8) |
Further, there are models that assume pores as truncated cone shapes [15] or consider the fractal dimension of pores [16]; however, there is a problem in that water content obtained from resistivity varies greatly depending on the assumed pore shape.
The most widely used empirical law constructed from the experimental values is Archie’s equation (eq. (9)) [17] obtained for sandstone saturated with brine.
| \begin{equation} \rho_{\text{eff}} = a\phi^{-m}\rho_{\text{W}}, \end{equation} | (9) |
where a and m represent the tortuosity coefficient and cementation exponent that change depending on the rock type, respectively. However, when the volume fraction of the liquid in the rock is very small, the effective resistivity of the water-bearing rock calculated from Archie’s equation becomes higher than that of the rock matrix. The modified Archie’s equation [18] that considers this has been proposed and is expressed as
| \begin{equation} \rho_{\text{eff}} = \left(\rho_{\text{R}}^{-1}(1 - \phi)^{\frac{\log (1 - \phi^{m})}{\log (1 - \phi)}} + \rho_{\text{W}}^{-1} \phi^{m} \right)^{-1}. \end{equation} | (10) |
Although empirical laws are useful, the physical image of fluid distribution and connectivity remains unclear, and therefore, it is necessary to investigate how these change with increasing fluids for various rocks [19].
We used digitized data from the Fontainebleau [20], Mt. Simon [21], and Bentheimer sandstones [22] for the calculations. Details on the binarization process of the CT images are presented in the cited references. CT images include three-dimensional data composed of voxels (cubes) binarized into pores and rock parts. The grid size of the Fontainebleau sandstone is 5.7 × 5.7 × 5.7 µm, and the porosity of the models ranges from 8–26%. Similarly, the resolution of Mt. Simon sandstone is 1.74 µm with porosities of 7–10%, and that of Bentheimer sandstone is 1.66 µm with porosities of 24–26%. The latter two datasets are downsampled by spline interpolation to unify the grid size to 5.7 µm for comparison with the Fontainebleau sandstone data. The model size must satisfy the representative elementary volume (REV), which is the minimum volume that can represent an entire medium. Sandstone has an REV ∼1 mm [23], and in this study, we performed calculations on models of an ∼1.4 mm cube. We prepared 56 sets of models that include 240 × 240 × 240 voxels from Fontainebleau sandstone and 8 sets each from Mt. Simon and Bentheimer sandstones. For the Fontainebleau sandstone, eight sets of models that include 240 × 240 × 240 voxels were created by changing the extraction location from large models with seven different porosities (8, 10, 13, 15, 18, 21, and 25%). For the Mt. Simon and Bentheimer sandstones, we created eight sets of models comprising 240 × 240 × 240 voxels by changing the extraction locations from large models with porosities of 13.4 and 24.9%. The x, y, and z axes defined in the large model were unified with those of the extracted model. Further, the coordinate axes were not defined along the bedding plane of the rock.
Figure 1 shows representative pore structures for each porosity of the Fontainebleau sandstone. The black areas indicate pores and confirm that the ratio of pores increase with an increase in the porosity of the model.
Three-dimensional digital models of Fontainebleau sandstones. The number in each axis denotes a grid distance, whereas φ represents the porosity of the model. (online color)
The resistivity of digital rocks was calculated using the finite element method [24, 25]. Assuming that charge conservation holds in the porous media, the current density vector J follows
| \begin{equation} \nabla \cdot \boldsymbol{J} = 0. \end{equation} | (11) |
| \begin{equation} \boldsymbol{J} = \sigma \boldsymbol{E}. \end{equation} | (12) |
| \begin{equation} \boldsymbol{E} = - {\boldsymbol{\nabla}} V. \end{equation} | (13) |
| \begin{equation} {\boldsymbol{\nabla}}\cdot (\sigma {\boldsymbol{\nabla}} \mathrm{V}) = 0. \end{equation} | (14) |
We solve eq. (14) by applying a potential difference to both ends of the model with periodic boundary conditions and minimizing the energy gradient through iterative calculations using the conjugate gradient method. The current density in eq. (12) was calculated for the entire digital rock from which the resistivity of the model was derived. The pores are assumed to be saturated with seawater (resistivity = 0.2 Ω·m), and a solid phase is assigned 105 Ω·m as a representative resistivity of the dry rock [26]. We calculated the resistivity in the three directions (x, y, and z axis) by changing the direction of the applied voltage for evaluating the anisotropy of the pore structure.
3.3 Tortuosity calculationWe analyzed the tortuosity τ2 of the current path from angle θ between the directions of the voltage application and current vector at each voxel using [27]
| \begin{equation} \tau^{2} = \left(\frac{l}{x} \right)^{2} {}= \frac{1}{\cos^{2}\theta}. \end{equation} | (15) |
There are two axes where the tortuosity is obtained for one axis on which the voltage is applied, and therefore, we calculated this for all voxels and considered the average tortuosity of the entire model. We used a value obtained by taking a weighted average according to the current magnitude for removing the effect of the tortuosity calculated from very weak currents, as shown by
| \begin{equation} \tau = \frac{\sum_{i=1}^{N}\tau_{i}w_{i}}{N}, \end{equation} | (16) |
where i and N represent the voxel number and total number of voxels, respectively. Weight wi is normalized by the maximum current. Further, θ = 90 ± 0.1° is excluded from the calculation to remove local errors.
Figure 2 illustrates an example of the calculation results for the local current distribution generated when a potential difference is applied in the y-axis direction. The current distribution in a color map indicates weak currents of less than 1% of the maximum current is shown as transparent. Further, the amount of current flowing increased with an increase in the porosity of the model, and the current flowed throughout the entire model.
Simulated electric currents for the digital models of Fontainebleau sandstones. φ represents the porosity of the model. (online color)
Figure 3 indicates the conversion of the current distribution shown in Fig. 2 into a vector distribution on the xy plane passing through the midpoint of the z axis. The current vectors, rock, and pores are presented in red, gray, and white, respectively. The vectors are scaled in length every eight points to prevent overlapping, and their lengths indicate the magnitude of the current. In models with high porosities, such as 21 and 25%, most current directions are aligned with the direction of voltage application (y axis). In models with low porosity, such as 8 and 10%, many currents perpendicular to the direction of voltage application (x axis) are confirmed.
Maps of simulated two-dimensional streamlines of electric currents for Fontainebleau sandstones. φ represents the porosity of the model. (online color)
Figure 4 presents a plot of the resistivity of the Fontainebleau sandstone calculated from the average current in the direction of voltage application against the porosity of the model. A difference in the markers indicates the direction of voltage application. The resistivity increases exponentially with decreasing porosity, and this trend is the same regardless of the direction of resistivity. In addition, the resistivity obtained in the models with low porosity showed large variations, while that obtained in the models with high porosity demonstrated small variations. A linear relationship was confirmed between porosity and resistivity on a double logarithmic graph; however, the slope differed on the low-and high-porosity sides with a boundary of ∼13% porosity.
Changes in the resistivity of Fontainebleau sandstones calculated in each direction as a function of porosity. (online color)
Figure 5 shows a plot of the anisotropy of the resistivity of Fontainebleau sandstone against porosity. Anisotropy is defined as the difference between the maximum value ρmax and minimum value ρmin of the resistivity in the three directions divided by the average ρmean of the resistivity in three directions.
| \begin{equation} \text{anisotropy} = \frac{\rho_{\text{max}} - \rho_{\text{min}}}{\rho_{\text{mean}}}. \end{equation} | (17) |
Changes in the resistivity anisotropy of Fontainebleau sandstones as a function of porosity.
As shown in Fig. 4, anisotropy increases with a decrease in porosity. The cases in which the value exceeded 1 were confirmed only in models with porosities below 10%. In models with porosities above 10%, most values were between 0.1 and 0.4.
4.3 Tortuosity of the current pathThe tortuosity τ2 of the current path was analyzed using eq. (15) based on the calculated current distribution (Fig. 3). Figure 6 presents a plot of the obtained results against the porosity. The tortuosity increased with a decrease in porosity in the direction of voltage application. The relationship between porosity and tortuosity is not a simple linear relationship, and the rate of increase in tortuosity increases with a decrease in porosity. In addition, the variations in tortuosity values were larger for models with lower porosity and smaller for models with higher porosity.
Changes in the electrical tortuosity of Fontainebleau sandstones calculated in each direction as a function of porosity. (online color)
Figure 3 confirms that current passing through rock flows along the pores. In models with high porosity, the current flowed along pore paths connected in the y-axis direction. In models with low porosity, many currents were tortuous, which can be interpreted as pores being insufficiently connected in the direction of the voltage application (y-axis). Moreover, the lower the porosity, the larger is the calculated resistivity and anisotropy (Figs. 4 and 5), and the larger is the tortuosity of the current path (Fig. 6). These results confirm that the connection of pores is hindered and the tortuosity of pores composed of connected pores increased with a decrease in porosity, thereby resulting in an increase in resistivity. The increase in the anisotropy of the resistivity with a decrease in porosity indicates that the connectivity at low porosity may be anisotropic. The anisotropy of tortuosity was also high at low porosities (Fig. 6).
We fitted the resistivity in each direction of the Fontainebleau sandstone (Fig. 4) using Archie’s equation to further examine the calculated resistivity values. For sandstone, a and m are reported to range from 0.5 to 2.5 and from 1.3 to 2.5, respectively [27, 28]. Applying the fitting to the porosity range of 17.2–25.5%, where the slope becomes linear on a double logarithmic graph. The resistivities in the x-, y-, and z-axis directions were fitted with a = 0.467, 0.432, and 0.456 and m = 2.216, 2.116, and 2.188, respectively. These values are consistent with those reported in previous studies [29].
Figure 7(a) shows a plot of the obtained resistivity against porosity, along with the fitted Archie’s equation, modified Archie’s equation with m obtained from fitting, and Hashin–Shtrikman model for the completely connected state (lower bound), tube model, film model, and percolation model. The Hashin–Shtrikman model for the completely isolated sphere model is not plotted here because their values deviate from the range of the vertical axis in Fig. 4, which is almost the same as that of the rock matrix. In addition, f = 1 is used for the tube (eq. (2)), film (eq. (3)), and percolation models (eq. (6)).
(a) Plots showing the resistivity of Mt. Simon, Bentheimer, and Fontainebleau sandstones in each direction and comparison with various rock physics models as a function of porosity. (b) Plots showing the tortuosity of Mt. Simon, Bentheimer, and Fontainebleau sandstones in each direction as a function of porosity. (online color)
When porosity falls below 10%, the resistivity values in all directions deviate from Archie’s law with fitted empirical parameters. This indicates the difficulty in applying Archie’s equation to low-porosity rocks and explaining the resistivity of porous media using Archie’s equation with a single set of a and m. The modified Archie equation with an empirical m does not reproduce the slope of the low-porosity resistivity on a double logarithmic axis. Comparing the calculated resistivity with models other than Archie’s equation, such as the film, tube, percolation, and Hashin–Shtrikman model (lower bound) indicated that they had slopes different from the resistivity at low porosity. This suggests that neither simple geometric shapes nor empirical laws derived from high-porosity sandstone can adequately explain the resistivity behavior across a wide range of porosity.
To verify that these findings are not specific to Fontainebleau sandstone but are reproducible across different sandstone types, calculated results for the Mt. Simon and Bentheimer sandstones are plotted in Fig. 7(a). These are all calculation results when the voltage is applied in the y-axis direction. Tortuosity is plotted in Fig. 7(b). Both resistivity and tortuosity agree well with the values of the Fontainebleau sandstone. The point at which the slope of the double logarithmic axis differs between low and high porosities is consistent; these results confirm that the discussion thus far is universal, beyond the Fontainebleau sandstone.
We estimated resistivity by substituting the calculated tortuosity of the current path into the equivalent channel model (eq. (7)) for verifying how well resistivity explained by tortuosity reproduces the resistivity shown in Fig. 4. Figure 8 presents a plot of the resistivity estimated from the equivalent channel model and resistivity shown in Fig. 4, along with Archie’s equation. Although rock physics models based on simple geometric shapes or empirical laws obtained under specific conditions cannot explain the resistivity values from high porosity to low porosity (Fig. 7(a)). The resistivity estimated using tortuosity and the equivalent channel model reproduces resistivity over a wide range from high porosity to low porosity. Further, it reproduces the point where the slope on the double logarithmic axis changes between low and high porosity sides. The non-linearity of the porosity-resistivity relationship (non-Archie behavior) is indicated in several studies, and attempts have been made for reproducing this non-linear behavior by increasing the parameters [30]. The results of this study can explain the resistivity of water-bearing porous rocks with various pore structures with a single parameter by utilizing the tortuosity of the pore path. In other words, the tortuosity of pores in the pore structure are an important factor to govern the resistivity of porous rocks.
Comparison of simulated and estimated resistivity values from the equivalent channel model. (online color)
Figure 9 shows the change in resistivity colored by the value of tortuosity to discuss the relationship between the deviation from Archie’s line in resistivity and tortuosity. The tortuosity at which the resistivity of Fontainebleau sandstone begins to deviate from Archie’s equation was ∼4.84–6.76 for the resistivity in the x- and y-axis directions, and ∼4.84–6.25 for the resistivity in the z-axis direction. The corresponding porosity is 13–15%. The deviation from Archie’s equation becomes more pronounced when the tortuosity reaches 9. Converting a tortuosity of 6.25 to an current angle, the angle of the current direction with respect to the direction of the potential difference is about 66°. Figure 3 shows that current vectors with an angle of 66° or more are dominant in models with a porosity of 8 or 10%.
Correlations between resistivity and porosity with respect to the tortuosity, shown as a color bar. (online color)
This implies that the tortuosity of the pores that become the current path in the water-bearing rock increases with a decrease in porosity, making it difficult to explain its resistivity with Archie’s equation using empirically obtained a and m.
We discuss the tortuosity at which the deviation from Archie’s equation begins by comparing it with geometric models. Substituting the tube and film models into the equivalent channel model (eq. (7)), revealed that the tortuosity in the connected state (f = 1) is 9 and 2.25 for the tube and film models, and the tortuosity of 4.84–6.76 is between these values. Figure 2 confirms that the current is connected in a tube-like manner when the porosity is lower than 13–15%, and the proportion of surface-like current increases when it is higher. These results suggest that the tortuosity of 4.84–6.76 may indicate the boundary at which the dominant pore shape in the resistivity of porous media transitions.
Previous calculations using models composed of combined elliptical particles [31] have examined the physical meaning of m in Archie’s equation, showing that larger m values indicate more rapid changes in pore complexity with varying porosity. The fact that m explains why the resistivity was different for low and high porosities in this study can be attributed to the difference in the tortuosity of the current path, resulting in different sensitivities of the complexity of the current path to changes in porosity. The results of this study suggest that tortuosity at which this sensitivity transitions is 4.84–6.76.
We calculated the electric current distribution flowing in digitized sandstones with a wide range of porosities and rock resistivities using digital rock physics. We found that resistivity and anisotropy increased with decreasing porosity. In addition, the slope of resistivity on the double logarithmic axis differed between low and high porosities, and the calculated current distribution confirmed that the connection of the pores was interrupted with a decrease in porosity. Accordingly, the current gradually deviated from the direction of the voltage application. The tortuosity calculated from this current distribution increased with a decrease in porosity. The results suggest that the lower the porosity, the more hindered is the connection of pores, which increases tortuosity, thereby resulting in higher resistivity and anisotropy. The resistivity of high-porosity samples could be fitted with Archie’s equation using empirical parameters; however, the resistivity of low-porosity samples could not be fitted, and the representative geometric models were similar. Thus, models that assume simple geometric shapes or empirical laws obtained from high-porosity sandstone cannot uniformly describe rocks with a wide range of resistivities.
In contrast, the equivalent channel model using tortuosity calculated from the current distribution could reproduce resistivity over a wide range of porosities and demonstrated that the tortuosity of the pore structure is an important factor to explain the electrical properties of porous rocks. The tortuosity at which the slope of resistivity on the double logarithmic axis changes was an intermediate value between the film and tube models, and the calculation results of this study suggest the process of the transition of the pore shape that dominates resistivity.
This research was supported by JSPS KAKENHI (Grant Numbers JP22K20383, JP22K14635, JP22H05303, and 24K01417).
