A digital rock physics approach to effective and total porosity for complex carbonates: pore-typing and applications to electrical conductivity

Recent advances in micro-CT techniques allow imaging heterogeneous carbonates at multiple scales and including voxel-wise registration of images at different resolution or in different saturation states. This enables characterising such carbonates at the pore-scale targeting the optimizing of hydrocarbon recovery in the face of structural heterogeneity, resulting in complex spatial fluid distributions. Here we determine effective and total porosity for different pore-types in a complex carbonate and apply this knowledge to improve our understanding of electrical properties by integrating experiment and simulation in a consistent manner via integrated core analysis. We consider Indiana Limestone as a surrogate for complex carbonate rock and type porosity in terms of macroand micro-porosity using micro-CT images recorded at different resolution. Effective and total porosity fields are derived and partitioned into regions of macro-porosity, micro-porosity belonging to oolithes, and micro-porosity excluding oolithes’ rims. In a second step we use the partitioning of the micro-porosity to model the electrical conductivity of the limestone, matching experimental measurements by finding appropriate cementation exponents for the two different micro-porosity regions. We compare these calculations with calculations using a single cementation exponent for the full micro-porosity range. The comparison is extended to resistivity index at partial saturation, further testing the assignment of Archie parameters, providing insights into the regional connectivity of the different pore types.


Introduction
Electrical resistivity measurements of partially saturated rock are commonly used to estimate the productivity and size of oil and gas reservoirs. While the empirical Archie's laws [1] describe the resistivity behavior of clean sandstones -- where FF is Formation Factor, RI is Resistivity Index, Rw is fluid resistivity (water), R0 and Rt are the resistivity of the fully and partially water-saturated rock respectively, is porosity, Sw is water saturation, a is a constant, m is cementation exponent and n is saturation exponent -carbonate rocks frequently display non-Archie RI behavior attributed to porosity at different length scales, e.g. "micro-porosity" and/or heterogeneity [2][3][4][5][6][7][8][9][10]. A good summary of different behavior types has been given by Fleury [5] including 'type II' (Fig. 1) with an initially lower slope bending up towards lower water saturation, and then flattening again. The author introduced a triple porosity model to explain the behavior where a series of conductors consisting of micro-porosity and macroporosity are in parallel to a micro-porous conductor. Sen [4] modeled micro-porous grain assemblies using effective medium approaches, where different RI shapes with n being a function of saturation resulted, but no 'type II' curve was evident. Fig. 1. Schematic of different RI curve types after Fleury [5] with type II bending at intermediate saturation, then flattening again with lower water saturation.
In recent years digital core analysis has made significant advances in particular for the analysis of heterogeneous rocks with complex micro-structure like carbonates. Numerical capabilities to carry out electrical conductivity calculations directly on tomographic images were established [11][12][13][14][15] and morphological transforms applied to carry out resistivity index calculations [16][17][18]. Integrated approaches combining experiments with tomogram-based numerical modeling were presented for non-Archie behavior of clean sandstones [19][20][21] including for pressure-dependent return to Archie behavior of Fontainebleau sandstone at high confining pressures [21] and more recently also carbonates [22], where the effect of films on the RI behavior of dualporosity carbonate rock is assessed with network models derived from tomograms acquired via micro-CT.
The analysis of transport properties like electrical conductivity for heterogeneous carbonates via digital core analysis requires balancing field of view (FOV) and resolution. Given the heterogeneity of carbonates, we consider core plugs of 1in diameter. [a] [b] [a] SEM image illustrating the main porosity types and solid phases; Ach: allochems, detrital skeletal remains of marine organisms, CaCm: athigenic calcite cement. Squares 1-4 have side length of 40 µm and depict various states of the environments of outer shells of ooliths (rims, 1-3) and isolated vugs (4). [b] Slice (1000x1000 voxel) through a 10003 subdomain of a tomogram of a 1in diameter core with cubic voxel size of 10.7 µm.
Calculation of resistivity index behavior from tomographic images then requires the application of effective medium type approximations for unresolved porosity at image resolution and the ordering of microporosity in terms of fluid invasion; e.g. where different micro-porosity types exist one of the micro-porosity types may be invaded preferentially (at lower pressure) by a non-wetting fluid. In this work, we distinguish between two micro-porosity types that are prominent in terms of their X-ray density and geometrical features in order to assess their relative impact on resistivity index behavior (see Fig. 2). In particular, we test the hypothesis that type II RI behavior can be reproduced directly from tomographic images by carefully addressing associated numerical problems originating from the choice of resolution; namely strong partial volume effects at partial saturations.

Sample description
The samples used in this study originate from an outcrop block of Indiana limestone (ILS) sourced from the Salem formation near Bedford, Indiana. We carried out XRF elemental analysis showing that ILS is basically a monomineralic rock with 98.8% calcite and spurious amounts of quartz (0.5%) and clay minerals (0.7%). SEM images of ILS were acquired of resin impregnated polished samples after bubbles were carefully removed by an applied vacuum and surfaces polished to roughness below 40 nanometres, using an FEI Nova NanoSEM230 system. A series of SEMs at different resolution were acquired and one SEM at intermediate resolution (200x magnification) is depicted in Fig. 2a. In addition to macroscale intergranular porosity different micro-porosity types are clearly visible. Due to their spatial extend, we select the most striking low density boundaries of ooliths, sometimes partially dissolved or collapsed, as separate micro-porosity class termed 'rims' in the following. Other micro-porosity types are more difficult to distinguish in the micro-CT images (see Fig. 2b) and therefore combined as 'non-rims'. Fig. 2a shows regions of cementation extending significantly into the macroporosity, and which may make some of the macroporosity inaccessible through larger porosity pathways connected to the outer boundary of the sample. Mercury intrusion capillary pressure (MICP) measurements (Fig.  3a) indicate a bi-modal pore accessibility distribution function with significant micro-porosity accessible at invasion radii below 2 micrometre [23].
To distinguish between total and effective porosity we measure porosity on 12 Indiana limestone cores via the so-called resaturation method (a setup utilising a desiccator and a vacuum pump) resulting in an effective porosity of 16.26% ± 0.25%. The SEM image suggests that there may also be very poorly connected or disconnected porosity present. Consequently we measured porosity both on intact and crushed cores using essentially a thermogravimetric method to evaporate the connate water (capillary bound and surface adsorbed water); comparing methods we find connate water of 0.28% ± 0.02%. Crushing the core frees another 0.32% of water that was not dried even at 160C; we therefore assume that fluid inclusions may exist and total porosity may indeed be larger than effective porosity for this rock (Fig. 3b). [a] [b] Poorly connected or isolated porosity quantified by thermogravimetric means, derived by drying with temperature increments over 24h for two sets of intact cores and crushed core.

Electrical measurements
Cores of 25.4mm diameter and length were drilled from the Indiana limestone block and cleaned by Soxhlet extraction [24] for 9-12h with a 50/50 solvent mixture of methanol and toluene, then placed in a humiditycontrolled oven at 60C. All cleaned and dried cores were vacuumed in a desiccator at 0.005 kPa over 12 hours and then submerged and soaked in brine (2% NaCl by weight), followed by pressure saturation at 2,000 psi to assure 100% saturation. Electrical resistivity measurements were carried out with a modified Coretest AEP-710 Electrical Properties System and resistance measured by a 1920 Precision LCR Meter at frequencies of 100Hz, 1kHz and 10kHz. We report measurements at 1kHz.
Measurements on fully saturated cores were repeated every 6-12h until resistance values varied less than ±3%. [a] [b] The resultant Formation Factor measurements are shown in Figure4a and compared to the data of [24]. While the actual FF values are different, they are on trend with similar cementation exponent. Resistivity index measurements were performed using the same electrical property measurement system where desaturation was achieved using a gravimetric capillary pressure system. To ensure stable water saturation, electrical measurements at each desaturation point were performed once observed pipette volume of the produced fluid did not change over an 8h period. A switch to the next desaturation point occurred once subsequent desaturation values were within ±2% saturation after another 12h. Water saturation at each desaturation point was determined by weight to avoid dead volume effects. The RI measurements including a repeat run are reported in Fig. 4b and compared with literature [8]. We find repeatable results which agree well with literature. Type II behavior [5] is clearly evident with an inflection point of the RI curves at a water saturation of SW = 45 %. This approximately coincides with the partition of porosity into macro-and micro-porosity (e.g. their relative fractions from MICP measurements).

X-ray micro-CT image acquisition and processing
A micro-CT image of Indiana limestone was acquired at a resolution of 10.67 m [22]. The tomogram of 2394 x 2385 x 3800 voxel covers the full core plug; here we consider a 1000 3 central cubic subdomain (see Fig. 2b). Image segmentation is carried out using an active contour method [25] seeded by initial thresholds followed by an additional segmentation of the non-resolved porosity into n micro-porous phases (Fig. 5b). Given this n-phase segmentation, we select a high porosity range of the micro-porous voxels as seed voxel set to detect voxels corresponding to the 'rim'-type micro-porosity. The rims are made compact by applying a morphological closing operation on the selected seeds (Fig. 5a). Such defined micro-porosity classes and n-phase micro-porosity field are overlaid in Fig. 5c and relabeled to continuous labels (Fig. 5d). The resulting labeling scheme is given in Table  1. To introduce partial saturations, we model drainage into water-wet ILS using the morphological capillary drainage transform (CDT) [15,16]; boundary effects are avoided by calculating the CDT on a larger 1200 3 subdomain which symmetrically enclosed the 1000 3 subdomain. The light blue in Fig. 5e highlights macro-porosity not connected to the 1200 3 domain boundary through macroporosity. Fig. 5f depicts a cross-section of the domain at intermediate water-saturation. Motivated by the high porosity of the rims and larger-scale porosity showing on the SEM images we model still lower water saturation by desaturating class 1 micro-porosity for the final water-saturation point. This allows considering saturations below the measured RI inflection point (Figure 4b). [a] [b] [c] [d] [e] [f]

Formation factor calculations
We solve for electrical conductivity using a Laplace solver previously described in [11,12,15,20]. In particular, we use a MPI-parallel conjugate gradient finite element method directly on the segmented image (labels given in Table 1) to solve the Laplace equation with fixed voltage gradient aligned to the z-axis of the domain. For the assignment of voxel-scale electrical conductivity we apply Archie's law (Eqn. 1) with independent cementation factors m1 and m2 for each micro-porosity class. The prefactor a in Eqn. (1) is set to 1. Due to a mismatch between imaged porosity (for the full imaged core) and experimental porosity, we scale the image porosity of each voxel proportionally to its solid fraction, considering an additional porosity of the macro-solid of Φm (voxel porosity Φv (v=micro) is rescaled according to Φv' = Φv + Φm (1-Φv). A value of Φm =0.02% is in agreement with the difference measured between effective and total porosity. Scenarios for the resulting voxel-scale electrical conductivity with fully effective or partially disconnected porosity are given in Fig. 6a, while Fig. 6b depicts the respective calculations on the 1000 3 domain for different cementation exponent pairs. Fig. 6b shows that without connected background porosity the experimental Formation factor measurements cannot be matched. For RI calculations we select reasonable combinations of m1-m2 pairs, e.g. the cases (m1 = 1.75, m2 = 2.00), (m1 = m2 = 2), and (m1 = 2.00, m2 = 1.75) with background porosity Φm = 0.02%.

Resistivity index calculations
Consider now the results of RI calculations given in Fig.  7a. Qualitatively there is a good match to experiment for the higher saturation range. It is also notable that the cementation exponent m1 of the rims has a relatively small effect on RI behavior with m2 playing a much more important role. This is not surprising given that the relative fraction of class 2 micro-porosity is significantly larger than for class 1. However, the inflection point of RI is not present for any of the of m1-m2 pairs. The reason for this behavior can be found by close inspection of Figure  5d: the bluish color is associated with boundaries between the macro-porosity and micro-porosity of class 2 and indicates the presence of strong partial volume effects. The latter do not matter significantly in FF calculations since neighboring voxels in the macro-porosity offer higher conductivity and there is some continuity. However, for increasing depletion of water with increasing non-wetting fluid saturation these highporosity boundary layers offer an effective shortcut for electrical current and thus result in artificially low RI values (Fig. 7a). [a] [b] Fig. 6. [a] Voxel-scale effective conductivity for different cementation factors against original porosity-bin values (see Table 1), m=mi (i=1,2): cementation exponent, and [b] resulting Formation factor of Indiana limestone, considering both different cementation exponent scenarios for each microporosity type and connected background porosity; σmicro/σfl=1/F: conductivity via Archie's law normalised to fluid conductivity mi: cementation exponent of class i.
We correct this effect by selecting voxels with high porosity of micro-porosity class 2 and set those voxels to zero conductivity for the lowest water saturation. This is justified since the micro-porosity of class 1 is already nonconducting for that saturation point as it is invaded by non-wetting fluid; the voxels to be invaded next are logically the ones on the surface of micro-porosity regions of class 2, which have an artificially high water saturation. Removing those high conductivity values while respecting the non-wetting fluid invasion order results in the RI values of Fig. 7b, which exhibits a transition of type II in RI behavior according to [5]. Finally, we compare in Fig. 8 the computed RI values to the experimental data, resulting in an excellent match to experiment. From this we can conclude that an excellent match is possible with the inflection point essentially located at the right water saturation value and slopes matched both at high and lower water saturation. [a] [b]

Discussion
The above derivation is based on an analysis of a subdomain of the full 1in diameter plug and as such still subject to heterogeneity effects at a scale greater than (10.7mm) 3 , which was the sample size considered here. This analysis is ongoing as well as a further tightening of model parameters utilizing a registered 8mm diameter sample at higher resolution. Here we base our analysis using effective parameters to work with the different scales of the sample, but constrain these effective parameters to reasonable values both in terms of background porosity (measured) and effective conductivities (extracted from sensitivity analysis). High resolution modeling will give additional clarity/confirmation on both the argument about the order of fluid invasion of the micro-porosity by drainage (including for the water-wet case considered here) and the effect of partial volume effects, which can be reduced by increasing resolution. Furthermore, voxels with strong partial volume effects were located here on the basis of a micro-porosity cut-off (>50%). While this is a reasonable approach, the set of voxels affected can be constrained further by selecting only voxels a distance smaller than √2 away from the macro-porosity boundary saturated by the non-wetting fluid which are also part of one of the micro-porous regions. An indication of such a filter based on Euclidean distance is given in Fig. 9. Visible in Fig. 9 are thin boundaries between the macro-porosity saturated by the invading fluid and the micro-porous layers, which are treated as extra phases. While in Fig. 9a these are combined with the non-wetting phase (since they do not contribute to conductivity and thus conductivity calculations are possible on this basis) in Fig. 9b these are represented as separate micro-porosity types with an associated expansion of the label range by an additional 202 phases. [a] [b] Fig. 9. Illustration of micro-porosity typing including the selection of boundary voxels of the invading fluid saturating macro-porosity.
[a] labeling using 205 labels (all invaded microporosity boundary voxels relabeled to same label) and [b] 405 labels (allowing to keep porosity scales consistent for convenience). This boundary assignment is a function of invading fluid saturation. Without going into detail about the color scheme, the main point to observe here is the boundary around the yellow-([a]) and red-toned ([b]) macroporosity, which is the distinguishing factor to Fig. 5. In particular, that boundary is located within the micro-porosity, and does not exist in the water-saturated macro-porosity; macro-porosity invaded by simulated capillary drainage at invasion radius of rCDT=1.5; slices depict a 1000x1000 voxel cross-section (10.7 µm resolution).

Conclusion
A careful analysis of resistivity index behavior of Indiana limestone was achieved by separately considering the properties of two different micro-porosity types, which show a distinct behavior in saturation order. This may be thought as being akin to the parallel/serial interpretation of 'type II' RI behavior [5], but taking into account the spatial distribution of these elements. The match to experiment was excellent and achieved through sensitivity studies varying both cementation factors of the two micro-porosity regions and global background porosity (proportional to voxel-scale solid fraction). We can draw the following conclusions: 1. The occurrence of partial volume effects has a significant effect on RI calculations and if not addressed is giving erroneous results. 2. The cementation exponent of the majority microporosity type (class 2, excluding rims) has a significant impact on RI behavior. 3. Adaption to mixed-wet systems will require a capability to model fluid distributions at high resolution over a large field of view. 4. The previous point may require the introduction of super-resolution techniques for resistivity modeling to capture both necessary detail and field of view.