Digital Carbonate Rock Physics

Modern estimation of rock properties combines imaging with advanced numerical simulations, an approach known as Digital Rock Physics (DRP). In this paper we suggest a specific segmentation procedure of X-Ray micro-Computed Tomography data with two different resolutions for 15 two sets of carbonate rock samples. These carbonates were already characterized in detail in a previous laboratory study which we complement with nano-indentation experiments. In a first step a non-local mean filter is applied to the raw image data. We then apply different thresholds to identify pores and solid phases. Because of a non-neglectable amount of unresolved micro-porosity (“micritic phase”) we also define intermediate phases. Based on this segmentation we determine porosity-dependent values 20 for Pand S-wave velocities as well as for the intrinsic permeability. The porosity measured in the laboratory is then used to predict the effective rock properties for a comparison with experimental data. Anisotropy is observed for some sub-samples, but seems to be insignificant in our case. Because of the complexity of carbonates we suggest to use DRP as a complementary tool for rock characterization in addition to classical experimental methods. 25 Solid Earth Discuss., doi:10.5194/se-2016-45, 2016 Manuscript under review for journal Solid Earth Published: 9 March 2016 c © Author(s) 2016. CC-BY 3.0 License.

measured in the laboratory. This conclusion is mostly based on Berea sandstone although carbonates are considered in this study as well. However, also Jouini et al. (2015) reports about an overestimation of effective elastic properties of carbonates by DRP. Therefore we conclude here that the digital rock images themselves and/or the computational workflow has to be improved to provide better estimate of effective properties of rocks. In this paper we consider in detail a carbonate dataset and suggest 5 techniques to achieve a better agreement between numerical predictions and laboratory measurements.

Rock Samples and Laboratory Characterization
2.1 Carbonate samples 10 measured at room pressure and temperature on core plugs 1'' (2.5 cm) in length and diameter (Vialle et al., 2013). The associated errors did not exceed 0.5%, 1%, and 2%, respectively. P-and S-wave velocities, at 1 and 0.7 MHz frequency, respectively, were acquired on dry samples under increasing (up to 30 MPa) and decreasing hydrostatic stress. Velocities were measured by using a pulse-transmission technique. The errors in V P and V S are about 1%. The results are listed in Tables 1 and 2. The obtained 5 grain densities (2.69 ± 0.01g/cm 3 and 2.70 ± 0.01g/cm 3 , for samples Carb-A and Carb-B, respectively) are in agreement with a mineralogy of pure calcite (Mavko et al., 2009).

Nano-indentation
Nanoindentation tests were performed to obtain stiffness (Young's modulus) of the carbonates at the micrometer-scale. These tests were performed on room-dried samples consisting of two small irregular 10 pieces, about 5 mm thick and with a surface of a few cm 2 , taken from the cuttings of the 1'' core plugs of Carb-A and Carb-B. Prior to testing, the surface of each sample was polished with carbide paper (grit 120). Roughness (Sq) of the surface measured by DS 95 AFM system (Semilab) on 10µm x 10µm areas was 1.4µm for Carb-A and not measured for Carb-B (RMS values). The IBIS nano-indentation system (Model B, Fisher-Crips Laboratories Pty.Ltd.) is equipped with a Berkovich-type diamond indenter 15 (Lebedev et al, 2014) and was used in a static mode: the tests consist in continuously recording the load, P, and the displacement, h, of the indenter as it pushes into and withdraws from the surface of the sample. A constant maximum loading force of 10mN and an initial contact force of 0.15mN were used.
In total, 961 (31 x 31) measurements were performed on a 300 x 300µm surface with a spacing of 10µm between measurement points. 20 Typically, the extraction of the mechanical properties is achieved by using the P-h curves and by applying a continuum scale mechanical model to obtain the indentation modulus M: where S the is unloading indentation stiffness S = dP dh ( ) h=h max and A c the contact area, extrapolated from the maximum penetration depth h max and using the relation A c =24.5.h max 2 according to the geometry of Berkovitch-type indenters.
Data was further corrected considering deviation of the indenter tip from ideal geometry, initial penetration into the rock below a load threshold and compliance of the loading column, leading to a 5 nominal uncertainty of indentation moduli of < 2GPa.
Young's moduli , , can be calculated from the indentation moduli according to Indenter properties are i= 1220GPa and i = 0.06, according to Klein & Cardinale (1992) and Fischer-Cripps, (2004) for diamond material. Each performed measurement covers a projected surface of about 10 40µm² on average (which corresponds to an equilateral triangle with a side of 6µm), and contains both pores and solid grains. A Poisson's ratio, , has also to be assumed for each nano-indentation measurement: even though we can compute it from the laboratory ultrasonic P-and S-wave measurements at the core scale, there is no reason to assume that this value is constant for each individual measurement. Figure 1 displays the distribution of the indentation moduli for both  and Carb-B (right).

Procedure to get CT-Images
Two samples were prepared for imaging with micro-XRCT from the cuttings, one from Carb-A and one from Carb-B. A cylindrical-shape sample of 1.5cm in height and 2mm in diameter was achieved by 20 gently grinding the cuttings, first on the side on a rock saw blade, and then by hand using sand paper (grit 120). This procedure allows obtaining very thin cylinders while minimizing mechanical damage that classical drilling would produce. These cylindrical samples were then glued with Crystalbond509 (SPI suppliers) on a 2mm diameter flat-head metal pin, which was itself inserted in the core holder of the micro-tomograph. The 3D X-ray Microscope Versa XRM 500 (Zeiss -XRadia) was used with a X-25 Solid Earth Discuss., doi:10.5194/se-2016-45, 2016 Manuscript under review for journal Solid Earth Published: 9 March 2016 c Author(s) 2016. CC-BY 3.0 License. ray energy of 60keV. Two different settings of source-to-sample distance were used to achieve two nominal voxel sizes of (3.4µm) 3 and (0.6µm) 3 for Carb-A, and of (3.4µm) 3 ("low" resolution) and (1.1µm) 3 ("high" resolution) for Carb-B. X-ray microscope and image acquisition settings are summarized in Table 3 for each of the two samples.
The number of radiographic projections acquired during sample imaging with "low" and "high" 5 resolutions were 3001 and 5001, respectively. The total scanning time for one sample was about 8h.
Initial cone-beam 3D image reconstruction was performed using internal "Reconstructor" software (Zeiss-XRadia). To illuminate geometrical artefacts during reconstruction a secondary reference was acquired for samples image with maximum resolution.

Segmentation procedure 10
In addition to solid grains and pore space different micritic phases are visible in the raw images of the scans entailing an advanced segmentation procedure. For our segmentation (Figure 2) we select a region of interest (ROI) from the raw data of the two types of carbonates with two different resolutions (Tab 3). The ROI is subdivided in eight partly overlapping sub-volumes, each of a size of 400 3 voxels ( Figures 3 and 4). For the low-resolution images, it gives thus eight cubes with a side of 1.37mm for 15 both samples, and for the low-resolution images, 8 cubes with a side of 0.25mm for Carb-A and of 0.46mm for Carb-B.
Our segmentation workflow is applied to the full ROI including all subvolumes. Image enhancement and segmentation steps were carried out using the software package Avizo Fire 9 (FEI Visualization Sciences Group). Before actual segmentation the image noise and scan artefacts are reduced while 20 preserving interfaces using a non-local mean image filter calibrated for the appropriate dimensions and kernel window sizes. Note that every step of image enhancement changes the original data set affecting subsequent steps required for data analysis.
The image-enhanced datasets were segmented into classes using global thresholds for the covered range • high-confidence mineral (illustrated with dark red color), • and five classes in between the mineral class and the pore class.
A non-negligible part of the pore space is below the resolution limit of the µ-CT scans. cf. results of mercury intrusion, Figure 6 of Vialle et al., 2013, showing pores and pore throats as low as 0.06µm.

Numerical Results 5
In order to numerically calculate the effective intrinsic permeability k of the digitized rock sample we calculate the fluxes under creeping flow condition based on a parallelized Stokes-solver. The parallelized Finite Differences-based Stokes-solver is suitable for the calculation of effective hydraulic parameters for low and high porous materials (cf. Osorno et al., 2015). Using volume averaging technique, we coarse-grain the local velocity field u(x) obtaining the global velocity component u m in 10 flow direction.
The intrinsic permeability is calculated with Darcy's law: ( The pressure gradient Δ is imposed with pressure boundary conditions in the numerical simulations.
The dynamic viscosity of the pore fluid is . 15 To obtain effective P-and S-wave velocities of the digitized rock samples we use a technique described in detail in Saenger et al. (2011) and references therein. The basic idea of this approach is to study speeds of elastic waves through heterogeneous materials in the long wavelength limit (pore size ≪ wavelength) using the RSG FD algorithm (Saenger et al., 2000) for the simulation of elastic wave  From the results of the high-resolution samples, Carb-A and Carb-B, it could be observed that the variation in pores channels arrangement is significant and the permeability in the different subsamples of the same material does not necessarily increase with the porosity increment.

Elasticity
Several micritic phases have been identified in the raw images of the carbonate rock. The porosity of 20 these regions cannot be determined exactly, as some pores are below the resolution of the scans: typically, micrites exhibit pore sizes with a maximum diameter of 1µm (Moshier, 1989;Cantrell and Hagerty, 1999), and pore sizes as low as 0.06µm have been measured by MICP for the samples under investigation (Vialle et al., 2013). To account for the unresolved pore space we perform a number of two-phase simulations. For these simulations we assign vacuum properties to the pore phase, while the 25 rest of the digital image including the micritic phase will be assumed to be solid with the mineral properties of calcite (e.g. Andrae et al., 2013b). For the second simulation, we assign vacuum properties to the pores and the first micritic phase, while the rest will be assumed to be solid. We continue this way for all micritic phases, so that the last simulation assigns the mineral properties of calcite only to the high-confidence mineral phase. By this technique we obtain a porosity-velocity trend (Figure 8) for any random selection of high-resolution sub-samples for Carb-A and Carb-B: This porosity-velocity trend is exactly the same as has been observed for a carbonate dataset from a different location used in Saenger et al. (2014). In their paper this trend has been observed for three different resolutions (65nm, 1µm, 4µm).

Permeability
Similar to the procedure of the numerical simulation for elasticity we vary the sample porosity. This way we get six possible domains for each subsample with different porosities. To reduce computational times for the Stokes flow simulation we eliminate the disconnected pores. Some sub-samples solid-pore configurations with the lower porosities do not present connected pores, and we assume for the effective 15 permeability k=0.  Figure 11 (left hand side) shows that the Carb-A sample permeability is anisotropic with a variation between directions of up to 80%. In some sub-samples the 20 permeability value varies by up to two orders of magnitude.
The largest sample available for the low resolution scan of Carb-A is 2.4mm x 2.4mm x 2.4mm (689 3 voxels). The permeability calculated for this sample is 13D for a porosity of 0.173.

Elasticity
For the low-resolution scans we repeat the two-phase simulations for Carb-A and Carb-B as described 25 in section 4.1.2 ( Figure 10). Interestingly, the two-phase trend given by Equations (4) and (5) is Especially in the low-resolution case we expect to have images with a large amount of unresolved 5 porosity, mainly due to micritic phases. Therefore we perform multi-phase simulations and vary the porosity by assigning effective elastic properties to an interval of micritic phases (always starting with the class closest to the high-confidence pore phase). As effective medium approach we use the trend given by the simulations using two single phases only (Equations (4) and (5)). Despite the interval of micritic phases we assign vacuum values for the high-confidence pores and use for the remaining 10 phases the known elastic moduli of calcite (e.g. Andrae et al., 2013b). The results are displayed with green dots for Carb-B in Figure 10. We repeat the procedure with different intervals of the micritic phases. There are three interesting observations: (1) The resulting effective velocities are always significantly below the observed two-phase trend, (2), the curves for different intervals will intersect each other, and (3), the experimental determined velocities for high confining pressures are between the 15 multi-phase results and the two-phase trend as described above.

Discussion
In this paper we compare results from laboratory investigations with numerical estimates based on digital images. Note that in laboratory experiments we use samples on the cm-scale for the determination of permeability and ultrasonic velocities and compare it with DRP-predictions based on 20 images on the mm-scale. Especially because of the known heterogeneity of carbonates there is always a risk that the selected scanned area is not representative compared to the full sample size used for laboratory characterizations.

Discussion Experimental Characterization 25
Even with the highest resolution currently available in micro-XRCT imaging there will be a significant amount of unresolved pore features which need to be treated in the DRP workflow (Saenger et al., Solid Earth Discuss., doi:10.5194/se-2016-45, 2016 Manuscript under review for journal Solid Earth

2015)
. On the grey-scale intensity level histograms of the low-and high-resolution images of the micro-CT scanning ( Figure 5) this is reflected in a continuum in the intensity levels between the phase identified as pores and the phase identified as calcite grains. In this paper we have dealt with these micritic phases by replacing, step-by-step, and in a cumulative way, each of the micritic phases by void, and establishing a porosity-velocity trend. The nano-indentation technique provides a measure for the 5 distribution of effective elastic properties at the micrometer-scale, and can thus potentially constrain the input parameters for the different phases identified during the segmentation. To be able to do so, nanoindentation needs to provide bulk and shear moduli from each of the measurements (load-displacement curves) and we need to obtain effective bulk and shear moduli values for each of the identified phases in the micro-tomography (pores, calcite grains and the five micritic phases). However, if nano-indentation 10 technique is a well-established technique in material sciences, which deals with homogeneous, purely elastic materials, this is, as of today, not yet the case for rocks, which are heterogeneous materials with both elastic and non-elastic behaviour (creep). Though nano-indentation tests provide significant insights into elastic properties of heterogeneous rocks such as carbonates (Lebedev et al., 2014;Vialle & Lebedev, 2015) or shale (Ulm & Abousleiman, 2006;Abousleiman et al., 2007) consistent with the existence of micritic phases identified with X-ray tomography. However, we did not observe two peaks in the histogram for the nano-indentation results (Figure 1) in contrast to the histograms of the scanned micro-XRCT-images ( Figure 5). Therefore the direct translation of moduli derived from nano-indentation remains also to be difficult. The "resolution" of nano-indentation used in 25 this study allows for determining effective elastic properties at slightly bigger scales than that used here for the micro-XRCT. Solid Earth Discuss., doi:10.5194/se-2016-45, 2016 Manuscript under review for journal Solid Earth Published: 9 March 2016 c Author(s) 2016. CC-BY 3.0 License. Regnet et al. (2014) showed that there is a relationship between micrite microstructure and laboratory ultrasonic velocities on core samples, with samples with higher content of 'tight micrite' exhibiting higher velocities and samples with higher content of 'microporous micrite' exhibiting lower velocities.
Studied core samples were though a mixture of different types of micrite and the measured velocities represent effective properties at the core scale. This observation is reflected in the established porosity-5 velocity trends (Equations 4 and 5): micritic phases with density closer to that of calcite ('tight micrites') have higher velocities than micritic phases with lower density closer to that of pores ("microporous micrites').

Discussion Porosity: Experiments vs. Digital Rock Physics
After the segmentation it is also possible to estimate the porosities of the samples. Based on the 10 suggested workflow described in section 3.2 there will be a lower and upper bound. For the lower bound we will treat only the "high-confidence pores" as pores; for the upper bound we count only "high-confidence minerals" as minerals.
We obtain a porosity range between 25% and 35% for the high-resolution data of Carb-A and a range between 7.5% and 31% for the low-resolution data. We observe that the mean value is in rough 15 agreement with the experimentally determined porosity of 16.7% (Vialle et a. 2013) only for the lowresolution case. Although also the experimental value comes with an error we conclude that the highresolution dataset for Carb-A is maybe not representative for the full sample. In case of Carb-B the intervals range from 13% to 45% and 7% to 48% for the high-resolution and low-resolution case, respectively. Here the mean value is in both cases closer to the experimentally determined porosity of 20 29.4%.
We conclude that the porosity values of carbonates using micro-XRCT-data will only provide estimates with a relatively high uncertainty due to the significant amount of unresolved pore features in the images. An indication is the result of the mercury-intrusion experiments presented by Vialle et al.

Discussion Permeability: Experiments vs. Digital Rock Physics
Permeability numerically estimated for Carb-B presents an error of 97% on average with respect to the experimental value. In some cases the error is as low as 55%. The numerical error in comparison with the experimental values is within the expected range for the numerical method at these porosities, cf. Table 1 in Osorno et al. (2015). 5 Experimental results for Carb-A sample are below the measurement error tolerance. This could imply a sample with no connected pores between the inlet and outlet defined for the experiment. The numerical estimation of the permeability for the Carb-A low resolution sample (on average 7 D) is four orders of magnitude higher than the experimental estimation. In the high resolution case the numerical estimation is closer to experimental results (on average 100 mD), but the porosity presents a large numerical error, 10 therefore we do not take this domain as representative of the sample. However the numerically calculated permeability does not differ much from values found in the literature for porous materials with alike porosity (cf. Andrä et al., 2013b). On the other hand it is observed that for a porosity below 25% permeability values of carbonates can span several orders of magnitudes (e.g. Figure 3 of Vialle et al. 2013). Therefore we suggest considering a statistically significant number of samples to characterize 15 a formation.

Discussion Elasticity: Experiments vs. Digital Rock Physics
There are two important observations. The two-phase trend (displayed with solid and dashed-dotted lines) seems to be an upper bound for velocities. This data-driven upper bound is much stricter than the bound given by Hashin-Shtrikman (see Figure 8) and is now confirmed for several carbonates using 20 several resolutions (this study and Saenger et al. 2014). Only for Carb-B we observe a slightly different trend for the low-resolution case for P-waves (equation (6) and Figure 10 right hand side).
The trend given by the envelope of the multi-phase simulations (displayed by green dots in Figure 10) is not a strict lower bound, because the shape will strongly depend on the applied method to determine effective elastic properties for areas which are below the resolution limit of the used  The best choice to our knowledge is the two-phase trend discussed above. This trend can be regarded as "carbonate-data-driven" effective medium approach. We suggest implementing here in the future also Solid Earth Discuss., doi:10.5194/se-2016-45, 2016 Manuscript under review for journal Solid Earth Published: 9 March 2016 c Author(s) 2016. CC-BY 3.0 License. the findings of the nano-indentation experiments. However, we observe that the velocities obtained for the multi-phase simulations are in a reasonable agreement with laboratory measurements. This is the case for a known porosity determined in complementary laboratory studies (see also section 5.1). For carbonates the distribution of the micritic phases and their effective elastic behavior is crucial to predict the effective wave speeds. 5

Discussion Anisotropy: Elasticity vs. Permeability
In general we don't observe any significant anisotropy for permeability and for velocities of the considered samples. However, a few samples are out of this general trend. One example is a sub-sample of Carb-A (low-resolution case), for which we show the results for P-wave velocities and permeabilities in Figure 11. Interestingly, the significant anisotropy for the permeability is not present for the velocity. 10

Summary
With the current imaging techniques it remains difficult to resolve microstructures (on sub-micrometer scale) and image a representative volume at the same time, which is essential to understand the effective material properties of rocks. The porosity of the rock samples is the most relevant parameter. To overcome this problem, we have conducted a careful calibration of DRP estimates with laboratory data. 15 For carbonate samples it is difficult to estimate the porosity from raw-CT data. Therefore, we use our presented numerical results in an inverse way. We suggest using the porosity determination from the laboratory and go back to our low-resolution results given in Figure 9 and 10. With a given porosity we can estimate the permeability and the effective wave velocities.
In case of the studied samples Carb-A and Carb-B, we can predict P-and S-wave velocities with a good 20 agreement to laboratory results. The predicted permeability values are only in good agreement for Carb-B. Most probably the low-resolution image of Carb-A is not representative for the sample used in the laboratory.
However, for carbonate rock the resolution of the XRCT is not sufficient for a more exact estimate because the micritic phases cannot be resolved. The effective elastic properties have to be 25 approximated. Our suggestion is to use the trend of the two-phase simulations. The implemented Solid Earth Discuss., doi:10.5194/se-2016-45, 2016 Manuscript under review for journal Solid Earth Published: 9 March 2016 c Author(s) 2016. CC-BY 3.0 License. workflow in this paper can be applied in general for numerical estimates of mechanical and transport properties of carbonates. Because of the known strong heterogeneity of carbonates we suggest to use a statistically significant amount of digital images to characterize a formation. Fig. 1: Nano-Indentation results for Carb-A (left hand side) and Carb-B (right hand side). In blue we 5

Figures
illustrate the corresponding moduli range from ultrasonic experiments on dry samples from 0 to 30MPa confining pressure, and in green we illustrate the moduli range given by the solid anisotropic calcite crystal. Overall we observe that the medium effective indentation module M is slightly stiffer for Carb-A (26GPa vs. 25GPa).  Tables   Table 1: Helium bulk and grain density (in g.cm -3 ), helium porosity (in PU) and air-permeability (in mD) measured at benchtop conditions, for Carb-A and Carb-B. Permeability of plug Carb-A is below the sensitivity level of the apparatus used (0.1 mD). Values are from Vialle et al. (2013