Digital Rock Physics: A Geological Driven Workflow for the Segmentation of Anisotropic Ruhr Sandstone

Balcewicz, Martin; Siegert, Mirko; Gurris, Marcel; Ruf, Matthias; Krach, David; Steeb, Holger; Saenger, Erik H.

doi:10.3389/feart.2021.673753

ORIGINAL RESEARCH article

Front. Earth Sci., 29 June 2021

Sec. Solid Earth Geophysics

Volume 9 - 2021 | https://doi.org/10.3389/feart.2021.673753

This article is part of the Research Topic Rock Physics and Geofluid Detection View all 27 articles

Digital Rock Physics: A Geological Driven Workflow for the Segmentation of Anisotropic Ruhr Sandstone

Martin Balcewicz^1,2^*

Mirko Siegert^1,2

Marcel Gurris¹

Matthias Ruf³

David Krach³

Holger Steeb^3,4

Erik H. Saenger^1,2,5

¹Reservoir Engineering and Rock Physics, Bochum University of Applied Sciences, Bochum, Germany
²Institut für Geologie, Mineralogie und Geophysik, Ruhr-Universität Bochum, Bochum, Germany
³Institute of Applied Mechanics (CE), University of Stuttgart, Stuttgart, Germany
⁴Stuttgart Centre for Simulation Technology, University of Stuttgart, Stuttgart, Germany
⁵Fraunhofer-Einrichtung für Energieinfrastrukturen und Geothermie IEG, Bochum, Germany

Over the last 3 decades, Digital Rock Physics (DRP) has become a complementary part of the characterization of reservoir rocks due to the non-destructive testing character of this technique. The use of high-resolution X-ray Computed Tomography (XRCT) has become widely accepted to create a digital twin of the material under investigation. Compared to other imaging techniques, XRCT technology allows a location-dependent resolution of the individual material particles in volume. However, there are still challenges in assigning physical properties to a particular voxel within the digital twin, due to standard histogram analysis or sub-resolution features in the rock. For this reason, high-resolution image-based data from XRCT, transmitted-light microscope, Scanning Electron Microscope (SEM) as well as geological input properties like geological diagenesis, mineralogical composition, sample’s microfabrics, and estimated sample’s porosity are combined to obtain an optimal spatial segmented image of the studied Ruhr sandstone. Based on a homogeneity test, which corresponds to the evaluation of the gray-scale image histogram, the preferred scan sample sizes in terms of permeability, thermal, and effective elastic rock properties are determined. In addition, these numerically derived property predictions are compared with laboratory measurements to obtain possible upper limits for sample size, segmentation accuracy, and a geometrically calibrated digital twin of the Ruhr sandstone. The comparison corresponding gray-scale image histograms as a function of sample sizes with the corresponding advanced numerical simulations provides a unique workflow for reservoir characterization of the Ruhr sandstone.

1 Introduction

Digital Rock Physics (DRP) and its use of non-destructive methods to determine permeability or effective elastic rock properties have become a complementary part in reservoir characterization over the past 3 decades (e.g., Blunt and King, 1991; Iassonov et al., 2009; Andrä et al., 2013a; Andrä et al., 2013b). Whether in estimating the total reservoir volume, potential reservoir production rates, or interpreting seismic lines above reservoirs, DRP contributes in all major subfields (e.g., Fan et al., 2010; Aljamaan et al., 2017). To cover the increasing demand for energy supply (e.g., Scheer et al., 2013; Araújo, 2014) and to achieve the climate targets that have been set (e.g., Rogelj et al., 2016), currently non-fossil reservoirs are studied, such as in geothermal energy projects (e.g., Fridleifsson et al., 2008; Goldstein et al., 2011). For this purpose, the determination of the thermal conductivity, as well as the heat capacity, is becoming increasingly important in DRP. In the course of future regional explorations of the geothermal potentials in North Rhine-Westfalia, western Germany (e.g., Knutzen, 2017; Lippert et al., 2019; Kruszewski et al., 2020; Weiß, 2020; Balcewicz et al., 2021a; Jüstel, 2021; Lippert et al., 2021), results from DRP will become more important for this region. The Ruhr sandstone is to be the start of a comprehensive DRP investigation of the Upper Carboniferous Ruhr cyclothem, which occupies a predominant part in the Rhine-Ruhr metropolitan region, featuring the largest European district heating network (AGFW, 2009). One of the essential questions to be answered in the coming years is: Which geological horizons are suitable for geothermal projects? DRP will be used to further investigate the specific rock properties of the Ruhr sandstone, which so far can only be answered with great effort using classical laboratory methods.

Although high-resolution X-Ray Computed Tomography (XRCT) technology has opened up many possibilities through the visualization of the third spatial plane, it is not a stand-alone tool for material characterization (e.g., Holzer and Cantoni, 2012). This is especially true for fine-pored phases in carbonates or grain and phase boundaries in general (e.g., Saenger et al., 2016b). Nevertheless, in order to use the three-dimensional information provided by XRCT, it is necessary to utilize additional sources of information. An already established source of information is laboratory measurements that can provide rock properties such as permeability or porosity at ambient or elevated pressure and temperature conditions (e.g., Ahrens et al., 2018; Schepp and Renner, 2021). However, such laboratory measurements are resource and cost intensive, which is why the application of DRP is particularly interesting for small sample quantities. Nevertheless, basic laboratory measurements are helpful for appropriate XRCT interpretation, but they are not sufficient. Recent studies have already shown that by using additional imaging techniques (e.g., transmitted light microscope or Scanning Electron Microscope, SEM), features unrecognized by the XRCT can be visualized and implemented into the digital rock model (e.g., Saxena et al., 2017; Karimpouli et al., 2019). However, an additional source of information should always be the rock’s provenance. Especially in geological formations, diagenesis as well as the current geological setting is of crucial importance and can lead to internal changes in the rock at all scales. This study, therefore, focuses specifically on geologic features of the Ruhr sandstone formed during the diagenesis and the possibility of incorporating corresponding observations from thin sections into the segmentation process. In the context of the sandstone’s microstructure, this study takes a special role compared to similar studies of the Berea sandstone or the Fontainebleau sandstone (e.g., Bourbie and Zinszner, 1985; Churcher et al., 1991; Øren and Bakke, 2003), since eight individual phases were assigned to the digital twin of the Ruhr sandstone, showing indications of pressure dependencies, microporosities, and phase transformations. However, the workflow presented in this study implies more than just the assignment of gray-scale intensities within the XRCT volume. Rather, it is a process that begins before the actual segmentation and ends with assigning physical mineral properties to the corresponding segmented phases, which properties may have to be adjusted due to the prevailing microstructure. The incorporation of geologically based information for geophysical interpretations is already well established in areas of image processing, such as seismic reflection image interpretation, large-scale velocity modeling, or time reverse imaging interpretations (e.g., Ritzmann et al., 2007; Finger and Saenger, 2020). In DRP, a comparable geological interpretation of gray-scale images is not well established, yet. Despite this, the combination of all information sources is indispensable to perform an appropriate rock characterization-especially when it comes to the interpretation of reduced image information such as the gray-scale images of a XRCT volume. To obtain meaningful results, a Representative Elementary Volume (REV) must be used for the respective numerical studies (Bear and Braester, 1972; Ostoja-Starzewski, 1999). A studied rock sample may show different grain sizes (from $<$ 0.1 to $>$ 60 mm) due to geologically induced processes. Increasing grain sizes must be compensated with an increasing sample size in order to meet representative sample criteria. However, increased sample sizes can lead to difficulties when using a XRCT [see Lambert–Beer’s Law e.g., in Buzug (2011)]. Nevertheless, optimal spatial gray-scale images must be extensively processed before advanced numerical simulation results can be provided. This processing requires detailed knowledge of the material in order to assign physical mineral properties to the appropriate gray-scale intensities. Regarding the Ruhr sandstone, the basis for the correct assignment of the gray-scale intensities to the corresponding phases is geological property information such as phase and pore fractions, phase transitions, or the microstructures of the Ruhr sandstone. The technical assignment of gray-scale intensities to the corresponding phases, called segmentation, is a fast-developing aspect of DRP. Numerous mathematical approaches have already been used, such as gray-scale-filtering, gray-scale histogram processing, watershed algorithms, or the use of neural networks (e.g., Nehler, 2018; Karimpouli et al., 2019). The strengths of automated segmentation methods that are less prone to operator bias are increasingly being developed and applied (e.g., Schlüter et al., 2014; Khan et al., 2016). Nevertheless, advanced machine learning segmentation algorithms such as k-means clustering, fuzzy C-means clustering, or artificial neural networks need to be trained with a sufficient database (e.g., Chauhan et al., 2016; Alqahtani et al., 2021), which is not ensured in the present study of the Ruhr sandstone. At this stage of the study, no novel segmentation workflow can be presented that can reasonably prove itself against other segmentation methods. Rather, a possible segmentation approach is presented how a corresponding training dataset for possible unsupervized machine learning segmentation methods could look like. Due to the complex microstructure of the Ruhr sandstone, this training dataset requires extensive manual processing, which is presented in this manuscript based on gray-scale thresholding. The presented method is highly susceptible to operator bias, but also allows manual assignment of gray-scale intensities in the complex microfabrics of the present Ruhr sandstone. Sensitivity studies have shown that even small segmentation inaccuracies due to operator bias can lead to significant errors in numerical characterization approaches, such as permeability measurements (e.g, Leu et al., 2014). For this reason, geologic verifications were implemented along the segmentation process to minimize the extent of inaccuracies. Numerous DRP studies focus on the distinction between pore space and solid phase, which is sufficient for numerical permeability approaches. However, in this study, a raw XRCT volume is used for one specific segmentation and then provided for modeling permeability, thermal, and effective elastic rock properties. Thus, the limitations will be discussed in more detail, which XRCT resolution (e.g., Saxena et al., 2018; Guan et al., 2019; Feng et al., 2020), segmentation accuracy (e.g., Leu et al., 2014; Garfi et al., 2020) or numerical application (e.g., Andrä et al., 2013b) is sufficient for the Ruhr sandstone.

Deriving permeability, thermal, or effective elastic rock properties from a digital twin require a five-step workflow: 1) Preparation of an XRCT volume, 2) tomographic reconstruction, 3) assessment and handling of the X-Ray artifacts 4) segmentation of pore and grain phases, respectively, and 5) solving equations due to the demanded properties. However, these strictly separated processes are increasingly becoming blurred by cross-process procedures, such as the process of reconstruction and segmentation (e.g., Yoon et al., 2010; Storath et al., 2015). This study will demonstrate that the derivation of a digital Ruhr sandstone twin cannot be reduced to just one of these five steps, but that the presented workflow is a cross-disciplinary method to produce reasonable numerical results for the complex microstructure of the Ruhr sandstone. Finally, the geological optimized digital twin can assist in minimizing the uncertainty of predicted rock properties at the mesoscopic and macroscopic scale in a non-destructive manner.

In the following, the geological formation of the Ruhr sandstone is presented first. Then, the Ruhr sandstone’s microfabrics and basic laboratory methods are described. Finally, after presenting the laboratory and numerical results, we discuss the workflow as a tool for deriving the digital twin of the Ruhr sandstone.

2 Materials and Methods

2.1 Geological Setting and Microfabrics Analysis of Ruhr Sandstone

The studied Ruhr sandstone from the Oberste quarry in Dortmund, Germany, is assigned to Namur B (Kaisberg) and Wesphalian A formation of the Upper Carboniferous. It is part of the Upper Carboniferous Ruhr cyclothem, which consists of clay-, silt-, mud-, sandstones and interbedded coal seams. These sedimentary deposits were chemically and mechanically compacted, folded, and faulted during the Hercynian orogeny. Therefore, paleo high burial depths of up to 6,000 m can be assigned to the Ruhr sandstone, which is located in a complex geologic setting in the Rhine-Ruhr metropolitan area (Drozdzewski, 1993; Karg et al., 2005).

The preparation of the studied thin sections was provided by the preparation facility of the Institute of Geology, Mineralogy, and Geophysics, at the Ruhr-Universität Bochum. Due to temporal tool limitations the thickness of the thin sections is 38 µm. This results in abnormal interference colors for the studied minerals (Figure 1A). The thin sections were digitized using the Keyence digital microscope, VHX-2000D, combined with a VH-Z100UR lens. The SEM images are open access and were acquired by the Ruhr-Universität Bochum using the high-resolution thermal-assisted field emission SEM Gemini2 Merlin HR-FESEM from ZEISS (Balcewicz et al., 2021c).

FIGURE 1

FIGURE 1. (A) Thin section of Ruhr sandstone shown in the upper half in crossed polars (XPL) and in the lower half in plane light (PPL). SEM images of Ruhr sandstone showing (B) broken quartz grains, (C) soiled pore space, and resulting (D) microporosities due to the sericitization process in feldspars. The abbreviations of the rock-forming minerals like Qtz (quartz), Ms (muskovite), Ab (albite), and Or (orthoclase) are used according to Kretz (1983).

Typically, Ruhr sandstone consists mainly of quartz ca. $70 \pm 5$ %, feldspar ca. $10 \pm 2$ %, and clay and feldspar alterations ca. $15 \pm 2$ %, and minor constituents of iron oxides ( $< 5 \pm 2$ %) (e.g., Kneuper, 1957; Rosenfeld, 1967; Duda and Renner, 2013; Stöckhert, 2015). The studied specimen consists mainly of quartz $70 \pm 5$ %, minor amount of feldspar phases like albite or orthoclase $< 5 \pm 2$ %, and $20 \pm 2$ % of feldspar altered phases like muskovite, sericite, kaolinite, as well as some illite, chlorite, and ankerite. Additionally, accessory minerals like zircon, tourmaline and pyrite are determined in the Ruhr sandstone (cf. Kneuper, 1957; Scherp, 1963; Rosenfeld, 1967). Furthermore, some grain conglomerates and a small amount of carbonates can be observed, which are most likely part of the feldspar altered matrix. The quartz grains show the largest shape with no systematic deviation from isometry or significant shape preferred orientation, and have a diameter of about 400 $\pm$ 200 µm. The clay-rich and sericite-rich matrix will be referred to as the sericitic phase in further discussions.

Quartz grains indicate both types of grain boundaries, high angle and low angle grain boundaries. Some grains show undulous extinction, other low angle boundaries have already transformed into quartz subgrains (Figure 1A). All intermediate steps of the transformation can be observed in the studied specimens. In addition, trans- and intragranular fluid inclusion trails can be observed in quartz (Figure 1A). Some high angle grain boundaries show sutured or concavo-convex grain boundaries. Further SEM analysis shows partially broken quartz grains in size between 10 and 100 µm (Figure 1B). The sericitic phase was pushed close to the quartz phase boundaries during diagenesis in the course of coalification, folding, overburdening, and alteration, indicating low porosity due to its cementing effect (e.g., Kneuper, 1957). In addition, the SEM images showed two types of pore spaces. On the one hand, free pore spaces corresponding to open cavities in the sandstone were found. On the other hand, pore spaces occupied by rock fragments such as debris or Quartzites could be identified (Figures 1B,C). The source of these fragments is mainly the host material. In the following, the pore spaces are referred to as clean pore (open cavities) and soiled pore (occupied by rock fragments). It is important to note that ankerite and bicarbonate solutions from mudstone deposits migrated into the Ruhr sandstone. Paleoenvironmental influences and precipitation of ankerite favored sericitization of the feldspars. This cementation by silica occurred at all depths and resulted in a quartzitic grain-to-grain structure observed in the microstructure analysis, which was found to be clearly temperature and depth dependent (e.g., Scherp, 1963; Hesemann, 2013).

2.2 Laboratory Measurements: Petrophysical Characterization of Samples

Cored samples were drilled from a Ruhr sandstone cube with an edge length of 0.15 m. At the time of drilling, the orientation of the cube specimen in the field was not known. Nevertheless, macroscopic anisotropy was suspected, so the cube was cored in three direction x, y, and z. A total of 24 samples with a diameter of 20 mm and a mean length of 80 mm were extracted by diamond core drilling. Additionally, a cylindrical core sample with a diameter of 5 mm and a length of about 10 mm was extracted from x direction of the initial Ruhr sandstone cube for a high-resolution XRCT image. All samples were sawed perpendicular to the sample axis and their faces were ground square to the maximum possible length. All steps of the preparation were performed with water as the cooling medium. Basic petrophysical properties were determined on identically prepared samples at ambient conditions.

Bulk density $ρ_{geo}$ was derived by the geometrical volume of the cylindrical samples and their dry masses. The grain density $ρ_{grain}$ was determined by pycnometer measurements from rock powder according to the German standard (Normenausschuss Bauwesen (NABau), 2011). The total porosity results from the ratio between bulk density and grain density according to $ϕ_{tot} = 1 - (ρ_{geo} / ρ_{grain})$ . The connected porosity was determined from the mass difference of dry and distilled water saturated samples (see Duda and Renner, 2013).

The ultrasonic P-and S-wave velocities, $v_{P}$ and $v_{S}$ , were picked by hand on dry samples based on the first arrival times of the detected wavefront. Therefore, an ultrasonic device consisting of a waveform generator, two identical broadband ultrasonic sensors and a digital storage oscilloscope was used. The two identical broadband ultrasonic sensors were placed on the two ground head surfaces, resulting in measurements parallel to the drilling direction. P-and S-wave velocities were calculated by dividing the sample length by the determined arrival times less the internal device run time.

Thermal conductivity $λ_{dry}$ , thermal effusivity $e_{dry}$ and specific heat capacity $c_{dry}$ was determined with a thermal conductivity scanner from C-therm using a modified transient surface source (MTPS) at ambient conditions. The samples were installed in the direction of drilling, i.e., the bottom surfaces of the samples were treated with a thermal paste and then placed on the MTPS. Within the MTPS, a helical heating element is surrounded by a guard ring that supports one-dimensional heat transfer into the sample. The small amount of heat supplied results in a temperature rise at the interface between the sensor and the sample, which induces a change in the voltage drop of the sensing element. The rate of rise of the sensor voltage, which is factory calibrated to temperature, is used to determine the thermal properties of the sample. The precision and accuracy are specified by C-therm and correspond to better than 1% and better than 5% respectively.

2.3 Imaging and Processing of Ruhr Sandstone

The scan of the sample was performed in a self-built, modular µXRCT system employing an open micro-focus tube FineTec FORE 180.01C TT with a tungsten transmission target from Finetec Technologies GmbH, Germany in combination with a Shad-o-Box 6K HS detector with a CsI scintillator option from Teledyne DALSA Inc., Waterloo, Ontario, Canada. Latter provides a resolution of 2,940 $\times$ 2,304 pixel with a pixel size of 49.5 µm (see Ruf and Steeb, 2020). The geometric magnification for the scan was set to 24.78 which leads to the highest achievable spatial resolution of about 50 linepairs/mm of the system for 10% Modulation Transfer Function (MTF) (cf. Rossmann, 1969). The corresponding voxel size results in 2 µm. The X-ray tube voltage was set to 80 kVp and the tube flux to 120 µA. The resulting X-ray spectrum was additionally modified by the usage of a 0.5 mm thick Aluminium-filter. In total 1,800 projections angles were used with an exposure time of 3,000 ms each. The reconstruction was performed with the filtered back projection method using the software Octopus Reconstruction (Version 8.9.4–64 bit), Vlassenbroeck et al. (2007). The resulting data set consists of 2,940 $\times$ 2,940 $\times$ 2,141 voxel and shows a physical volume of 5.88 $\times$ 5.88 $\times$ 4.28 mm³. Consequently, the sample was scanned over the entire diameter of 5 mm and allows basically the definition and extraction of an appropriate REV. The XRCT dataset and associated metadata of the scan settings are open access (Ruf et al., 2021). The corresponding gray-scale intensities of the µXRCT volume and one characteristic 2D slice are shown in Figures 2A,B.

FIGURE 2

FIGURE 2. (A) Gray-scale intensities after a subsequent $2 \times 2 \times 2$ binning operation resulting from the XRCT volume and a corresponding slice in (B) with an edge length of 800 pixel. (C) Segmented gray-scale intensities color-coded with corresponding segmented phases. Note that the color-coded boundaries reflect the selected gray-scale intensity, but should not be interpreted as strict boundaries between two phases. Rather, the phase boundaries show overlapping characteristics that were adjusted by morphological post-processing operations (Sheppard et al., 2004; Jones et al., 2007). These critical areas between two phases are indicated by the critical boundary regions (CBR). (D) Segmented result with corresponding gray-scale intensities (C) of slice (B). Additionally, high angle quartz grain boundaries are implemented by the watershed-algorithm in a white color. The figures were created with Avizo (2019) and Matlab (2020).

The voxel size $v_{size}$ of the tomogram is the result of the employed geometric magnification M and the given pixel size $p_{size}$ of the detector, $v_{size} = p_{size} / M$ holds. The maximum spatial resolution of the system is about 50 linepairs/mm at 10% of the MTF (cf. Ruf and Steeb, 2020). This corresponds to a resolvable feature size of about 10 µm. Therefore, most features should be preserved and no significant information loss should occur if the data set is sub-sampled to a uniform voxel edge length of 4 µm. With this motivation, a subsequent $2 \times 2 \times 2$ binning of the raw tomogram data was performed to reduce the data set size. Before this, a 1,600³ voxel subvolume was extracted from the center of the original XRCT volume which corresponds to x = (650, 2,249), y = (750, 2,349), and z = (0, 1,599). Thus, the final voxel size after the performed binning of the 800³ voxel subvolume is 4 µm. This allows easier data handling and due to the limitation of the numeric simulation schemes, with respect to the maximum number of possible degrees of freedom, to cover a larger field of view. In addition, the 1,600³ voxel subvolume was reduced again by a factor of $4 \times 4 \times 4$ , taking into account the loss of information. Thus, the 400³ subvolume’s voxel size is 8 µm. Both binned subvolumes are published open access (Balcewicz et al., 2021b). Considering the number of largest grains per edge length of the initial 1,600³ voxel subvolume, a REV of the Ruhr sandstone can be assumed (Saxena et al., 2019). A detailed REV analysis should follow the approach of Ostoja-Starzewski (1999), which is beyond the scope of this study.

2.4 Initial Ruhr Sandstone Segmentation

Based on the microfabrics analysis, eight phases were segmented in the respective 800³ and 400³ voxel subvolumes: 1) clean pore, 2) soiled pore, 3) quartz, 4) sericitic phase, 5) albite, 6) orthoclase, 7) pyrite, and 8) quartz high angle phase boundaries. All segmentation was done in Avizo (2019) with basic installed modules. Both segmented subvolumes are open access (Balcewicz et al., 2021b). All applied segmentation processes are based on simplified gray-scale thresholding methods which are probably the simplest way to binarize an image. Nevertheless, some adjustments were made to avoid known problems such as overlapping phase boundaries (e.g., Sheppard et al., 2004). The applied segmentation can thus basically be assigned to the supervised multistage approaches (e.g., Sezgin and Taşaltın, 2000; Wu and Amin, 2003). Starting with a global threshold value that roughly divides the Ruhr sandstone into the individual phases, followed by a local threshold to incorporate additional information such as geological input parameters into the segmentation. In the course of adaptive methods, morphological post-processing operators such as “remove small holes” or logical difference operators were used (Table 3).

After the initial gray-scale thresholding, a geologic review was performed before further steps could be applied. During the geological review, the segmented image was compared and evaluated with thin section images, SEM images, and porosity measurements. The geological verification may result in some modifications to the segmented image, such as fitting the pore space or elimination of incorrectly assigned phases due to X-ray artifacts. The geological verification was based on qualitative and quantitative evaluation criteria such as adjacent mineral phases or calculation of pore space (Table 1; c.f. Table 3). This process of adjusting the initial segmentation results was repeated until the geologic input parameters, such as phase or porosity fractions, approximately matched the binarized phases.

TABLE 1

TABLE 1. Matrix properties determined by laboratory measurements on Ruhr sandstone samples from the Oberste quarry in Dortmund. All measurements were conducted at ambient conditions. The abbreviations RSST and e.g., X indicate the rock sample Ruhr sandstone and the drilled direction from the parent cube, respectively.

After all individually segmented phases passed the geologic verification, a location-dependent model of the digital Ruhr sandstone twin was generated. This includes the merging of seven individually segmented phases into a single binary object. Afterwards, a watershed algorithm was applied to the segmented volume according to Beucher and Meyer (1992) to automatically assign unassigned voxels to the appropriate phases. The watershed algorithm is based on the idea of a gray-scale image interpreted as a topographic relief, provided that the pixel values are considered as altitudes. Based on this idea, the watershed algorithm calculates the catchment basin located between two maxima. These catchment basins define the image partitions. This results into a simplified digital twin of the Ruhrsandstone.

In a final step, the quartz grain boundaries were added to the Ruhr sandstone model. For this purpose, the quartz phase was isolated from the rest of the location-dependent digital twin and processed individually. The “separate objects” module was used to implement grain boundaries. This module is a high-level combination of watershed, distance transform, and numerical reconstruction algorithms which computes the watershed lines of a binary image. The selected properties are 3D interpretation, voxels with a least one common vertex are considered connected, marker extent is set to seven, output type is set to line, and a repeatable algorithm method for a conservative chosen method. Followed by an erosion module to set the quartz boundaries to a maximum size of one voxel. Finally, all segmented phases were merged resulting into the digital twin of the Ruhr sandstone (Figure 2D).

2.5 Modelling Permeability, Thermal, and Effective Elastic Rock Properties

Smoothed Particle Hydrodynamics (SPH) is a fully-Lagrangian mesh-free collocation method, which has proven feasible for the calculation of flows through porous materials. In SPH the continuous weakly compressible Navier–Stokes equations are transformed in discrete summation of short-range interaction forces.

For the calculation of the permeabilities the $400^{3}$ voxel subvolume was used which corresponds to a number of about $16 \cdot 10^{6}$ degrees of freedom (DOF) after subtracting non-essential solid fractions. For the simulation, the HOOSPH code was used in which the specific SPH-particle-particle interactions are implemented and that is based on HOOMD-blue (Anderson et al., 2008; Glaser et al., 2015). This is a rather computationally expensive method but has shown to be highly suitable for effective permeability computations of porous rocks in the Darcy and weak inertia regime as well as with onset of more pronounced inertia effects (Sivanesapillai et al., 2014; Osorno et al., 2020).

In order to determine the effective thermal conductivity, the steady-state heat conduction equation is solved using the finite volume method (Siegert et al., 2021). To limit the numerical error, the interfacial conductivity is modelled with two different averaging approaches, namely the harmonic mean (HM) and the arithmetic mean (AM). Each phase of the numerical model is assigned a fixed thermal conductivity to determine the effective thermal conductivity of the entire sample, these are: $λ_{air} = 0.026 W {(m K)}^{- 1}$ , $λ_{sericitic} = λ_{muscovite} = 2.28 W {(m K)}^{- 1}$ , $λ_{albite} = 2.14 W {(m K)}^{- 1}$ , $λ_{orthoclase} = 2.31 W {(m K)}^{- 1}$ , and $λ_{pyrite} = 19.21 W {(m K)}^{- 1}$ (e.g., Clauser and Huenges, 2013). Considering the observed reduced physical properties of quartz in Saenger et al. (2016a), the original literature value of the thermal conductivity of quartz was reduced by 27.75%, resulting in $λ_{quartz} = 5.56 W {(m K)}^{- 1}$ .

Effective elastic P- and S-wave velocity were derived by using the rotated staggered finite difference grid technique (Saenger and Bohlen, 2004). This numerical scheme allows modeling of seismic waves in complex viscoelastic and anisotropic media with strong contrast interfaces of arbitrary shape. All calculations were performed with second-order spatial finite-difference operators and with a second-order time update. The method was successfully compared with other approaches (see in Andrä et al., 2013b; Saxena et al., 2019). The segmented phases (A) clean pore, (B) soiled pore, (C) quartz, (D) sercitic phase, (E) albite, (F) orthoclase, and (G) pyrite were assigned with the corresponding elastic moduli c₁₁, c₄₄, and $ρ_{solid - phase}$ : (A,B) 0.0 and 0.0 GPa, (C) 68.0 GPa, 33.5 GPa, and 2,648 kgm^-3, (D) 184.3 GPa, 16.0 GPa, and 2,844 kgm^-3, (E) 74.0 GPa, 17.3 GPa, and 2619 kgm^-3, (F) 67.0 GPa, 14.3 GPa, and 2,575 kgm^-3, (G) 361.0 GPa, 105.2 GPa, and 5,016 kgm^-3, respectively (Bass, 1995; Saenger et al., 2016a). Additionally, for the 800³_watershed and 400³_watershed subvolumes the quartz grain boundaries were assigned with half the elastic moduli of the host mineral, i.e., $c_{11, watershed} = 0.5 \times c_{11, qtz}$ and $c_{44, watershed} = 0.5 \times c_{44, qtz}$ . This factor was determined within the numerical studies as a result of a first calibration approach.

3 Results

3.1 Laboratory Characterizations

Twenty-four samples were obtained from one Ruhr sandstone cube with an edge length of 0.15 m. Samples were drilled parallel to the three directions x, y, and z. All ten samples drilled in the x direction showed macroscopic indications of weathering, which decreased gradually from top to bottom. Except for this, no large fractures were observed in the Ruhr sandstone. The mean grain densities of the samples retrieved from x, y, and z direction are 2,676 $\pm$ 1, 2,632 $\pm$ 1, and 2,647 $\pm$ 1 kgm^-3, respectively (Table 1). Considering the measurement uncertainty and standard deviation, all studied samples showed comparable densities.

Mean total and connected porosity results show slightly higher total and connected porosities for samples drilled in x direction. The derived mean total porosities are 6.58 $\pm$ 1.87, 5.13 $\pm$ 1.90, and 5.37 $\pm$ 1.88% in x, y, and z direction, respectively (Table 1). Maximum and minimum derived total porosities show results of 7.89 and 6.07%, respectively. Connected porosity results show a similar trend, revealing 5.53 $\pm$ 0.78, 5.45 $\pm$ 0.77, and 5.10 $\pm$ 0.72% for samples taken from x, y, and z direction, respectively. However, a maximum connected porosity was determined in sample RSST-Z-05 with 6.42%. Corresponding minimum result was determined for sample RSST-X-01 sampled from x direction.

Ultrasound P- and S-wave velocities of dry Ruhr sandstone range between 4,667 and 2,303 ms^-1, respectively (Table 2). P-wave velocities measured in x direction appear to be up to 300 ms^-1 lower than corresponding measurements in y and z direction. P-wave velocities measured in y direction show the highest values of about 4,667 ms^-1. Equivalent measurements in z direction show a small variation from P-wave velocities determined in y direction. The mean P-wave velocities for x, y, and z directions are 4,363, 4,644, and 4,542 ms^-1, respectively. Similar to the P-wave velocities, determined S-wave velocities in x direction appeared to be the lowest among all measurements. Mean S-wave velocities for x, y, and z direction are 2,368, 2,507, and 2,558 ms^-1, respectively.

TABLE 2

TABLE 2. Ultrasound velocities, thermal conductivity, thermal effusivity, and specific heat capacity derived by laboratory measurements for Ruhr sandstone from Oberste quarry in Dortmund. All measurements were conducted at ambient conditions. The abbreviations RSST and e.g., X indicate the rock sample Ruhr sandstone and the drilled direction from the parent cube, respectively.

Thermal conductivity results show similar values for samples oriented in y and z direction. The derived maximum and minimum thermal conductivity for samples oriented in y and z direction are 3.64 and 4.03 W (mK)^-1, respectively. Mean thermal conductivity results for samples oriented in y or z direction are 3.96 and 3.87 W (mK)^-1 (Table 2). Compared to this, samples oriented in x direction show a lower thermal conductivity on average. Here, the mean thermal conductivity is 3.44 W (mK)^-1. The mean measured thermal effusivity for samples in x direction is 2,659 Ws^-1/2(m²K)^-1 with a minimum and maximum thermal effusivity of 2,601 and 2,678 Ws^-1/2(m²K)^-1, respectively. Samples oriented in y and z direction show a mean thermal effusivity of 2,902 and 2,865 Ws^-1/2(m²K)^-1, respectively. The maximum thermal effusivity was measured for sample RSST-Z-01 with a result of 2,941 Ws^-1/2(m²K)^-1. Average derived thermal heat capacities were measured for the same amount of samples as thermal conductivity and thermal effusivity. The mean thermal heat capacities for samples taken from x, y, and z direction are 822, 852, and 844 J (kg K)^-1, respectively (Table 2).

3.2 Digital Rock Physics Results

The presented segmentation resulted in the corresponding fractions for both 800³ voxel subvolume phases (A $-$ G) and 400³ voxel subvolume phases (see Table 3; A’ $-$ G’): (A) clean pore phase 0.38%, (B) soiled pore phase 6.11%, (C) quartz 72.71%, (D) sericitic phase 16.71%, (E) albite 1.38%, (F) orthoclase 1.48%, and (G) pyrite 0.05%. The implemented high angle quartz grain boundary fractions are 1.16 and 1.58% for both subvolumes 800³ and 400³ voxel, respectively. The total porosity showed 6.50 and 8.51% for 800³ and 400³ voxel subvolumes, respectively. For connected porosity further distinctions can be made like the considered connectivity between measured voxels. For example, considering that at least one common vertex should be analyzed for the segmented pore voxels, the connected porosity for the 800³ voxel subvolume is 4.65 percent. In case of the 400³ voxel subvolume, the connected porosity is 6.59% (Table 4). Considering the mentioned porosities, this results in a higher total and connected porosity for the 400³ voxel subvolume by 2.01 and 1.94%, respectively. All segmented label analysis was derived by Avizo (2019) modules [e.g., volume fraction, axis connectivity, and label analysis; see in Thermo Fisher Scientific (2018)].

TABLE 3

TABLE 3. Phase segmentation matrix with applied settings and derived phase quantities for the subvolumes of Ruhr sandstone 800³ and 400³ voxel.

TABLE 4

TABLE 4. Derived porosity measurements considering three cases of segmented pore voxels obtained for the 800³ and 400³ voxel subvolume, respectively.

Numerical permeability was calculated using the SPH method for the $400^{3}$ voxel subvolume with a total and connected porosity of 8.51 and 6.59%, respectively. Using a body force b as the driving force, all components of the intrinsic permeability tensor can be determined as follows

k^{I} = [\begin{matrix} 5.9710 & 0.00376 & 0.00519 \\ 0.00343 & 5.5045 & 0.00478 \\ 0.00294 & 0.00537 & 5.4968 \end{matrix}] \cdot 10^{- 12} m^{2} . (1)

The entries that do not lie on the main diagonal can be computed using the entries in the velocity vector perpendicular to the direction of the driving force. In this case, however, these are in the small single-digit per mille range relative to the main diagonal. Thus, a consideration of $k_{11}^{I}$ , $k_{22}^{I}$ , and $k_{33}^{I}$ is sufficient to make a statement about anisotropy of the investigated subvolume.

The segmented 800³ and 400³ voxel subvolumes were further subdivided into: 1) subvolume_watershed, corresponding to the segmented digital twin with high angle quartz grain boundaries and 2) subvolume_solid, without watershed segmented quartz grain boundaries. This results in a total of 24 simulations for the thermal conductivity. Results based on the AM-approach for the 800³_watershed voxel subvolume in x, y, and z direction are 4.049, 4.121, and 4.123 W (mK)^-1, respectively. Corresponding results for 800³_solid show 4.139, 4.246, and 4.247 W (mK)^-1 for x, y, and z directions, respectively. Lower results are derived for the simulations based on the HM-approach. Here, the results of the 800³_watershed voxel subvolume are 3.721, 3.696, and 3.706 W (mK)^-1 for the x, y, and z directions. The corresponding results for the 800³_solid voxel subvolume show higher values than for 800³_watershed and read as 3.922, 4.068, and 4.068 W (mK)^-1 for x, y, and z directions, respectively.

Derived thermal conductivity results for the corresponding 400³_watershed and 400³_solid show comparable tendencies for the individual measurements. Simulations based on the AM-approach vary in the range of 3.902 and 4.115 W (mK)^-1 for 400³_watershed and 400³_solid. Results of the HM-approach vary in the range of 3.512 and 3.880 for the corresponding 400³_watershed and 400³_solid subsamples. This results in total lower values for 400³ voxel subvolumes than for 800³ voxel subvolumes and generally lower values for segmented digital twins with high angle quartz grain boundaries than for those without (cf. Table 5 and Figure 3).

TABLE 5

TABLE 5. Correlation between laboratory and numerical derived results.

FIGURE 3

FIGURE 3. Correlation between (A) connected and total porosity, (B) P-wave velocity and total porosity, (C) thermal conductivity and total porosity, and between (D) permeability and confining pressure [see permeability measurements in Duda (2012), Brenne (2016), Lippert (2017)]. Filling and coloring of markers indicate sample and sample direction, respectively. Measurement uncertainties are indicated by error bars. Where error bars do not exceed marker size, uncertainties appear to be small. The blue line in panel (A) indicates identity, i.e., connected porosity corresponds exactly to total porosity. Filled and unfilled red triangles indicate different numerical approaches for 800³ and 400³ voxel subvolumes, respectively. Measurements in panel (A–C) were conducted at ambient conditions. Laboratory measurements in panel (D) have not been determined in the scope of this study and can be looked up for detailed information in Duda (2012), Brenne (2016), Lippert (2017). The abbreviations RSST and e.g., X indicate the rock sample Ruhr sandstone and the drilled direction from the parent cube, respectively.

Effective P-and S-wave velocities were modelled for 800³ and 400³ voxel subvolumes of the Ruhr sandstone. Similar to thermal conductivity models, the segmented subvolumes were subdivided into subvolume_watershed and subvolume_solid. Derived P- and S-wave velocities for the 800³_solid voxel subvolume were modelled without the segmented high angle quartz grain boundaries and result in 4,735, 4,845, and 4,865 ms^-1 for x, y, z directions, respectively. Determined S-wave velocities read as 3,210, 3,259, and 3,240 ms^-1 for the corresponding subvolume_solid in x, y, and z directions, respectively. Modeled P- and S-wave velocity in x direction with segmented quartz grain boundaries 800³_watershed result in 4,713 and 3,194 ms^-1, respectively.

Computed P-wave velocity results for 400³_watershed show 4,591, 4,740, and 4,729 ms^-1 for the three directions. The 400³_solid P-wave velocity in x direction is 4,620 ms^-1. Additionally, the P-wave velocity was determined for the 400³ voxel subvolume with assigned quartz grain boundaries of either 10% of quartz moduli or equal to pore moduli. Corresponding results are 4,540 and 4,500 ms^-1 for 10% of quartz moduli and equal to pore moduli, respectively.

4 Discussion

To develop a digital Ruhr sandstone twin a three-step interdisciplinary geological driven workflow is proposed, starting with the 1) pre-segmentation process, followed by the 2) assignment of gray-scale intensities to the resolved phases, and finally a 3) numerical characterization approach of the digital Ruhr sandstone twin in terms of permeability, thermal conductivity, and effective elastic wave propagation (Figure 4). In the pre-segmentation process, geological properties serve as input, such as Ruhr sandstone’s diagenesis, or mineralogical composition. Additionally, the raw XRCT data set must be technically evaluated based on the gray-scale histogram, existing X-ray artifacts, and the agreement with a REV. The gray-scale assignment to the corresponding mineralogical phases presented in this study is based on the simple gray-scale thresholding method, which allows a maximum of individual manipulation to ensure a first approximation of a meaningful segmented Ruhr sandstone model. The manipulation of the assigned phases is contrasted with four geological verification steps, which are intended to ensure a meaningful result during the segmentation process (Figure 4). The first geological verification is the number of expected phases based on the mineralogical composition and the accepted, high-quality XRCT image. The second geological review is during the individual phase segmentation based on visual evaluation of wrongly assigned voxels and the inspection of segmented phase fractions. The latter criterion results from the geological input properties. After individual phase segmentation, all phases are merged into a location-dependent model of the digital twin, which is supported technically by watershed algorithms and assisted by the third geological verification based on the geological property input. The last geological driven verification is the additional implementation of features, which cannot be resolved by the XRCT but might show a significant impact on further numerical characterizations. For Ruhr sandstone, this applies to high angle quartz grain boundaries.

FIGURE 4

FIGURE 4. (A) Segmented 800³ voxel subvolume with additional implemented high angle quartz grain boundaries in a white color. (B) Corresponding connected pore network. The ball sizes are indicating the pore size volumes. Visualized pore throats show the connectivity between two pore volumes. (C) Visualized connected porosity of panel (B) indicating critical pore throats in a subvolume. All figures were created with Avizo (2019).

4.1 The Pre-Segmentation Process: Evaluating the Reconstructed Gray-Scale Intensities

The reconstructed gray-scale intensities result from the materials’ mass attenuation coefficients which, in case of this study, can be defined as the phase attenuation coefficient $μ_{mineral}$ divided by the phase density $ρ_{mineral}$ . This leads to the fact that minerals with similar attenuation coefficients are located very close to each other in the gray-scale histogram. Comparable observations are made when fine structures cannot be resolved by the XRCT (e.g., micritic phases in carbonates, see in Saenger et al., 2016b). At perfect conditions, each gray-scale intensity peak represents a single phase. However, in XRCT-based characterizations, image noise and artifacts are always present and should somehow not be included in the segmentation interpretation (e.g., ring artefacts or beam hardening; Münch et al., 2009; Hsieh, 2009). For this reason, numerous image denoizing algorithms are evaluated in the literature (e.g., linear and nonlinear filtering methods; Mery, 2015). However, first applications of a 3D-based nonlocal-mean filter have shown that individual voxels associated with the sericitic phases became filtered. One common answer could be the achieved technical resolution of the used acquisition system, which could not sufficiently reproduce the grain contacts. Another explanation could be the tight grain-to-grain boundaries observed in the Ruhr sandstone’s microstructure due to the diagenesis. As a result, no distinct edges appear in the XRCT volume, making the edge-preserving smoothing by the nonlocal-mean filter give false results. This is especially true for those voxels whose size roughly corresponds to the resolution of the subvolume to be segmented. Due to the already applied binning procedure, we decided not to extract any further gray-scale intensities from the respective subvolumes and thus not to apply any of the established filtering methods (e.g., Schlüter et al., 2014; Saxena et al., 2018). However, the X-ray artifacts within the subvolumes could be handled and reduced to a minimum by the subsequent morphological post-processing operations during the segmentation process. Nevertheless, we would like to point out that the application of linear or nonlinear filtering methods during reconstruction or in the course of segmentation handles these artifacts much better than manual processing of the voxels.

The applied binning procedure increased the contrast of the 400³ voxel subvolume compared to the 800³ voxel subvolume at the sacrifice of imaged small particles in the XRCT volume. This results in a higher total and connected porosity for the 400³ than for the 800³ voxel subvolume (cf. Table 4 and Figure 3A). Results of the modeled thermal conductivity as well as the effective elastic wave propagation show that the chosen coarsening of the data set is not in disadvantage for the mentioned physical investigations. This suggests that both the $2 \times 2 \times 2$ and $4 \times 4 \times 4$ binned data sets are basically sufficient in terms of the resolution provided for thermal conductivity and effective elastic wave propagation simulations. This statement does not apply if the permeability is determined numerically. The resolved pore throats show that the binning procedure resulted in incorrect permeability simulations (cf. Figures 3D, 5C, Supplementary Material S1). At this state of the study, it is not possible to say whether the original data set (Ruf et al., 2021) has the necessary spatial resolution to provide sufficient resolution in the segmented rock model. We therefore advise a critical evaluation of the physical quantities to be studied and the resolution offered by the segmented data set, especially for small resolved features like micro porosities.

FIGURE 5

FIGURE 5. Schematic representation of the geologically driven workflow. It should be noted that geological input parameters such as geological diagenesis, mineral composition, existing microstructure, or estimated sample porosity are repeatedly consulted in the segmentation process in order to obtain meaningful results. The so-called geological verification thus takes place a total of four times during the threshold segmentation.

4.2 The Segmentation Process: Evaluating the Mineralogical Segmentation Based on Reconstructed Gray-Scale Intensities

Primarily, the applied segmentation workflow is based on the observed mineral structure in the examined thin sections. Minor differences may occur within the studied microstructures of the Ruhr sandstone due to local geologically influenced processes (cf. Kneuper, 1957). Subsequently, the reconstructed XRCT volume must be compared to the given microfabrics and the observations of the available reconstructed gray-scale histogram. The gray-scale histogram already gives a good indication of which phases could be resolved by the XRCT (see Figure 2A).

The presented workflow is based on a repeated application of gray-scale thresholding. In this approach, the minima and maxima of the gray-scale intensities of a phase are segmented, resulting in strict boundaries between the two phases. However, studies have shown that, e.g., in terms of noise in the gray-scale histogram due to multiple phase attenuation coefficients, this is not necessarily correct (e.g., Ketcham and Carlson, 2001; Sheppard et al., 2004; Jones et al., 2007; Wildenschild and Sheppard, 2013). Therefore, critical boundary regions arise during the segmentation of the individual phases (Figure 2C). Critical boundary regions close to two peaks can be approximated by a watershed algorithm. However, the main idea of the presented segmentation workflow is to reduce critical boundary regions between two peaks to minimize possible watershed errors, i.e., exact numerical assignment of gray-scale intensities, which, however, is not geologically meaningful. Furthermore, the presented segmentation method allows control over the segmentation process, which in our opinion is not necessarily the case with e.g., watershed-based segmentations. Jones et al. (2007) presented a method for minimizing critical boundary regions by the use of intensity–gradient histograms. The extent to which the described segmentation also provides meaningful results for the example of the Ruhr sandstone must be further analyzed. First segmentation attempts using the intensity-gradient-histogram did not yield satisfactory results. The reason for this could be that the described phase contrasts of the Ruhr sandstone’s minerals, i.e., the mean of total observed phase attenuation coefficients, is much lower than the phase contrasts between the scaffold of a bone and air in the gray-scale histograms of the examined biomaterials. Based on the presented segmentation idea, only those phases can be assigned which phase attenuation coefficients is sufficient to be detected by the given XRCT. Furthermore, these distinguishable phases must be large enough to exceed the occurring signal-to-noise ratio (SNR). The presented segmentation workflow has shown that phases that have very similar gray-scale intensities (e.g., both feldspar phases) can be separated clearly from each other, by manipulating false-assigned gray-scale intensities. The presented manipulation does not only imply the separation of phase contrasts at possible grain boundaries, but also assigns gray-scale intensities within grain boundaries to the correct phases - without a necessarily application of watershed algorithms. Due to the subsequent correction of mismatched gray-scale intensities, we expect a clear quality improvement of the segmented images compared to the simple threshold based phase separation e.g., multi-thresholding segmentation. Nevertheless, the segmentation workflow does not work entirely without watershed algorithms. Primarily, it is used to assign individual unassigned voxels to the appropriate ambient phase in the final process of segmentation. This step is mandatory to ensure that the resulting segmentation model does not contain any unassigned voxels that may have been created during processing. Although the critical boundary regions were previously reduced to a minimum, we cannot exclude the possibility of voxels being incorrectly assigned at this point due to the automation of the watershed algorithm. However, with regard to the corresponding model sizes of 400³ and 800³ voxel subvolumes, these incorrectly assigned voxels are assumed to be neglected.

A much more critical process of segmentation occurs during the implementation of grain boundaries (i.e., high angle quartz grain boundaries). Madonna et al. (2012) already presented a possible implementation of unresolved features within a XRCT image, almost ten years ago. We have adopted the described process by Madonna et al. (2012), which is based on watershed algorithms and assumed pressure-dependent quartz contact zones (Saenger et al., 2016a). However, numerical P- and S-wave velocity modelling results indicated further pressure-dependent features. These numerical results are in a good agreement with the observed quartz microfabrics (e.g., low angle quartz grain boundaries). Based on the geological setting, we assume at least a third pressure-dependent phase - namely the sericitic phase.

The used threshold segmentation can be applied to a single XRCT dataset with fairly inexpensive effort. This differs from the approach used for advanced machine learning algorithms, which require a critical amount of training data to produce satisfactory initial results. Consequently, the presented segmentation workflow should present a possible approach how future neural network based segmentation methods can be trained. The obvious advantages of unsupervised segmentation methods are the expected consistency, repeatability, and comparability of the segmented results (e.g., Chauhan et al., 2016; Alqahtani et al., 2021). These arguments gain importance especially when a segmented model of the digital twin is already accepted. However, these arguments can only apply to a limited extent when it comes to the first approximations to create a new verified 3D binary. This is where geological verification comes in for DRP. The steps presented in this study are intended to illustrate that operator bias can be well estimated by microfabric analyses (phase fractions or phase transformations) or by simple laboratory experiments (e.g., determination of porosity) and does not have to be characterized as arbitrariness of the operator. Finally, with respect to the results from the microfabric analysis or laboratory measurements, it can be considered that the raw XRCT data set has been reasonably segmented. In the future, these initial observations need to be applied to a larger data set and verified with appropriate methods. Here, the application of neural network based segmentation methods can make a major contribution to future Ruhr sandstone segmentations.

4.3 The Capabilities of the Geological Driven Workflow

Diagenesis reconstructions have shown that the Ruhr sandstone was covered at depths of up to 6,000 m (Karg et al., 2005). Assuming a normal gradient, this results in a pressure of about 150 MPa, which was applied to the Ruhr sandstone. Karg et al. (2005) also reconstructed possible temperatures of above 120°C that must have acted on the Ruhr sandstone. These P-T conditions fit both the presented microfabrics and the generally accepted highly complex geologic setting (e.g., Balcewicz et al., 2021a, and references therein). Saenger et al. (2016a) have already shown that pressure-dependent indicators of a Berea sandstone sample at 20 MPa confining pressure cannot be resolved in the XRCT. This resulted in pressure-dependent grain contact zones. The presented microfrabric analyses also show highly stressed quartz grains. Implemented high angle quartz grain boundaries decrease the results for the effective elastic wave propagation as well as for the numerically derived thermal conductivity approach. These results fit to the insulating character of thin segmented phases assigned with elastic moduli properties of air. To be specific regarding thermal conductivity results, comparing the two averaging methods, it also becomes clear that the HM-approch reflects the insulating character of the contact phase much more strongly than the AM-approach. The reason for this is the low resolution of the contact phase. In such a case, the AM-approach cannot fully capture the isolating effect. However, additional low angle quartz grain boundaries, subgrains, fluid inclusion trails, and sutured grain boundaries could be identified for quartz grains compared to the observations in Saenger et al. (2016a). For this reason, we conclude that it is mandatory to assign reduced elastic properties to the digital rock twin in order to achieve meaningful agreements with laboratory results. This results in reduced elastic properties for the quartz phase by ca. 30% (cf. Bass, 1995; Saenger et al., 2016a). It must be pointed out that the degree of reduction derived in Saenger et al. (2016a) is not necessarily correct for the Ruhr sandstone (Figure 3B). An explicit numerical calibration of the elastic quantities must follow this study of geometrical calibration to answer open questions. The purpose of this study is to demonstrate that the geologically driven workflow potentially opens up the possibilities to answer subsequent questions. A corresponding numerical calibration should therefore also be investigated for all further segmented phases. Due to the geological diagenesis of the Ruhr sandstone and the resulting segmentation possibilities we want to inspire other research groups to focus on the soiled pore phase and the sericitic phase for future studies.

SEM images (Figures 1B–D) have shown that the pore space in the Ruhr sandstone occurs in at least two different types. The process of segmentation can be approximated as a function of the applied time vs. the location-dependent accuracy of the segmented phases. For this reason, only two pore space observations were segmented: 1) clean pore and 2) soiled pore phase. While the explanation of the clean pore space and further handling for numerical purpose is obvious, more extensive aspects have to be considered for the soiled pore space. Particles located in this soiled pore space can have different influence on the rock properties. Mineral fragments from the Ruhr sandstone are likely to be present within the soiled pore space. However, these particles may seal off the connected porosity during permeability measurements, which would effectively result in lower permeabilities being measured in the laboratory. This observation should be taken into account, especially when evaluating numerical permeability results with comparable microfabrics. Besides sealing the connected porosity, the mineral particles can also influence hydrochemical reactions. Tran et al. (2020) studied mine waters in the Rhine-Ruhr metropolitan area and the influence on the local rheological units (e.g., Ruhr sandstone). First results showed that mine waters were weakly acidic to neutral, resulting in significant water-rock interaction correlations. Transport or hydrochemical induced clogging indicate that numerical permeability results tend to be too high compared to laboratory measurements (Figure 3D). In this respect, we recommend explicit investigations of particle transport or mineral scaling in the Ruhr sandstone based on XRCT and other imaging techniques (e.g, Weinhardt et al., 2021). Until then, the segmented soiled pore phase provides an opportunity for other research groups to adapt their numerical approaches to the observed Ruhr sandstone’s microfabric and incorporate possible reduced permeabilities.

Another phase that was segmented due to the workflow, but needs to be studied in more detail as part of a numerical calibration, is the sericitic phase. The segmented sericite phase consists of a mixture of different mineral phases. Basically, this is due to the feasibility of XRCT resolution of microporosities and extremely fine-grained mineral phases (e.g., pore diameter $<$ 1 µm and mean grain size $<$ 5 µm). In addition, however, the sericitic phase also exhibits distinct paragenetic sequences that can be attributed to the metasedimentary background of the Ruhr sandstone. Blatt (1967) has already pointed out the importance of provenance studies on recycled sediments and that identical mineral phases may occur at different shapes and with different properties. Thus, considering the described Ruhr sandstone’s diagenesis of Kneuper (1957) and Scherp (1963) and the large number of phases combined in the segmented sericitic phase, the digital rock twin can only be described as an approximation of the real Ruhr sandstone. SEM images show enclosed fine pore spaces by the sericitic phase and quartzitic-sericitic grain-to-grain structures which show further possible pressure-dependent zones (Figures 1B,C). Due to simplifications, we have assumed for the elastic properties of the sericitic phase the most represented mineral, i.e., muscovite. This means that observed phases such as kaolinite, illite, or chlorite are not considered for the numerical approaches. Therefore, we have reasonable concerns to suspect that the numerical results of the thermal conductivities as well as P- and S-wave velocities are too high for the digital rock twin. Reduced elastic properties and thermal conductivities could be a solution to the discrepancy between laboratory and numerical results - which, however, we consider to be the wrong approach for this phase. Instead, the individual phases must be segmented in an additional process (based on thin section, SEM or XRCT) and then implemented as high-resolution subvolumes into the digital rock twin. Karimpouli et al. (2019) has already shown that a similar method can be used to implement pore spaces in a XRCT volume which resolution is insufficient for small features (e.g., XRCT volumes from medical tomography scanners).

4.4 The Numerical Characterization Approach: Verification of Observed Anisotropy

During preparation, macroscopic indications of weathering were observed assuming macroscopic anistropy. Laboratory results show a weak anistropy for samples drilled in x direction. This is especially true for derived P- and S-wave velocities (cf. Figure 3B and Table 2). Similar observations were made for laboratory derived thermal conductivity results. Samples taken from x direction tend to show lower thermal conductivity results than samples from y and z direction, respectively (Figure 3C). However, this observation cannot be made for the comparison between total and connected porosities (Figure 3A). Advanced numerical approaches to determine permeability, thermal conductivity, and P- and S-wave velocity were applied for a cylindrical core sample parallel to laboratory determined sample results in x direction. Due to the cylindrical shape of the scanned sample, true y and z directions could not be reconstructed for the digital rock twin. The numerical results have shown a fairly good agreement with the results obtained in the laboratory, considering the XRCT resolution as well as the assignment of the elastic moduli to the segmented phases. In addition, weak anisotropies were determined independently from the respective numerical approaches.

The thermal conductivities in y and z direction show almost identical results, while the thermal conductivity in x direction seems to be lower. In the case of implemented high angle quartz grain boundaries, this effect seems to be suppressed for the thermal conductivity results. The numerically derived P- and S-wave velocity results do also show weak but clear anisotropy for x direction. Modeled P- and S-wave velocity show up to 100 ms^-1 lower results for x direction of the 800³ voxel subvolume and up to 200 ms^-1 lower results for x direction of the 400³ voxel subvolume compared to y and z directions, respectively. These observations are in good agreement with the mean velocity differences determined in the laboratory of approx. 300 ms^-1, which are smaller for the x direction than for the y and z directions. Löer et al. (2018) have demonstrated the degree of anisotropy using installed seismic arrays in the Parkfield area. This verified anisotropy, which could be determined from surface wave measurements, is in the range of about 100 ms^-1. Thus, because of the explicit alignment of our samples with respect to the assumed anisotropy, we can demonstrate a comparable degree of anisotropy. The influence of an anisotropy of about 300 ms^-1 is shown by the interpretation of a velocity model in a geothermal field in Los Humeros (Mexico). The combination of a geological subsurface model and a geophysical velocity model shows a good agreement between both interpretations, which provide either a resolution of tens of meters or hundreds of meters per second, respectively (Löer et al., 2020).

Due to the observed macroscopic anisotropy of the samples, further investigations were considered. Thus, it was initially assumed that the anisotropy occurring in the color of the rock samples was caused by the intrusion of the pyrites. A comparison of the modeled P-wave velocities, in which all phases except the pyrites were assigned to quartz, did not provide clear evidence for this hypothesis. In a second set of tests, all phases were set equal to quartz except for the sericitic phase, which was shown to be the result of alteration of the feldspars. This comparison of the modeled P-wave velocities also did not confirm any anisotropy due to the solid phases. In a final numerical run, all phases were assigned equal to quartz moduli except for the two pore phases, i.e., clean and soiled pore. The result for this analysis reads as 4,484, 4,615, and 4,626 ms^-1 for x, y, and z directions, respectively.

The applied SPH based permeability solver confirms the observed anisotropy in the x direction. Nevertheless, the comparison between the numerical and laboratory derived results reveals that not all physical quantities can be derived with the same degree of sophistication using the identical segmented model. In particular, the resolution of micro porosities and the associated pore throats make, for example, the segmented 400³ voxel subvolume an inappropriate choice for such permeability simulations. However, the same segmented data set gives reasonable results for both numerical thermal conductivity and effective elastic wave propagation results. From this, we conclude that valuable numerical thermal conductivity or effective elastic wave propagation results can be expected from segmented volumes after binning processes as described in this study. We consider this observation to be equivalent to the lower resolution in XRCT images. Whereas numerical permeability approaches at this reduced resolution can no longer provide reasonable results. This implies high resolutions in microtomography processing are necessary to visualize critical volume zones in the pore throats for very low porosity rocks like Ruhr sandstone (see supplementary materials).

Thus, we summarize that the main anisotropy is due to the Ruhr sandstone’s pore network (Figure 5). This result indicates that the apparent seismic anisotropy at the field-scale (Löer et al., 2018) can not only stem from fault systems but from the pore network, too. Future studies may analyze the implications regarding fluid flow directions in reservoirs.

5 Conclusion

Based on the example of the Ruhr sandstone, we presented a simple segmentation approach that takes into account geological properties such as mineralogy, microfabrics, and basic laboratory measurements like total and connected porosity. To create the digital twin of the Ruhr sandstone sample we have prepared a cylindrical core sample with a diameter of 5 mm and a length of about 10 mm to provide an appropriate REV. The created XRCT volume resulted in a 800³ voxel subvolume with a voxel size of 4 µm. Due to the limitation of the numeric simulation schemes a second 400³ voxel subvolume with a voxel size of 8 µm was prepared. The geological driven workflow was applied on both subvolumes, primarily based on gray-scale threshold segmentation. By the additional use of e.g., logical difference operators or watershed segmentation, a digital Ruhr sandstone twin with eight phases could be acquired. The presented segmentation workflow is intended to serve as a first approach to train future neural network-based segmentation approaches with sufficiently accurate segmentation phases.

The laboratory results indicate a transversely isotropic medium whose properties vary only with the x direction studied, but not with the y or z directions studied. After comparing the laboratory measurements with the numerical results, a weak anisotropy was confirmed by the numerical approaches, i.e., modeling permeability, thermal, and effective elastic rock properties. Further investigations suggest that the observed anisotropy is due to the existing pore network in the Ruhr sandstone. Moreover, we present the complex microstructure of the Ruhr sandstone with its altered and clay-rich matrix on the basis of thin section and SEM images. On the one hand, the assigned properties are intended to reflect the investigated rock properties, such as low porosity, and on the other hand, they are intended to provide suggestions for matching the numerical results more closely to the laboratory measurements. We encourage other research groups to take advantage of our segmented datasets and use the elaborated phases to address open questions, such as how to assign elastic moduli for pressure- and temperature-dependent phases or how to simulate the effect of clogging or hydrochemical reactions during flow modeling.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author Contributions

Seven authors have contributed to this paper. MB was the lead author and carried out laboratory measurements, sample characterization, segmentation, and analytical calculations. MS helped with numerically derived thermal conductivity. MG was involved in interpreting numerical derived results. MR provided high-resolution XRCT images. DK was involved in numerical permeability calculations and their interpretations. HS provided helpful background knowledge and laboratory as well as numerical resources. ES was involved in effective elastic wave propagation simulations, provided helpful background knowledge, and was involved in the data interpretation. All authors contributed to writing the manuscript.

Funding

MB and ES acknowledge generous funding by the Federal Ministry of Education and Research and geomecon GmbH for the project 3D-RuhrMarie (“FHprofUnt2016”; Project No. 13FH279PX6). MR and HS acknowledge funding from the German Science Foundation (DFG) through Project No. STE 969/13-1. HS thanks the DFG for supporting this work under Grant No. SFB 1313 (Project No. 327154368). The authors declare that this study received funding from geomecon GmbH. The funder was not involved in the study design, collection, analysis, interpretation of data, the writing of this article or the decision to submit it for publication.

Conflict of Interest

The authors declare that this study received funding from geomecon GmbH. The funder was not involved in the study design, collection, analysis, interpretation of data, the writing of this article or the decision to submit it for publication.

Acknowledgments

We thank all students and researchers involved in the laboratory measurements, in particular Mandy Duda. We especially thank Claudia Finger for fruitful discussions and helpful comments on an earlier draft of the manuscript. Furthermore, we would like to thank the editor Pascal Audet and the two reviewers Steffen Berg and Arne Jacob for their comments and discussions on the manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2021.673753/full#supplementary-material

References

AGFW (2009). Arbeitsgemeinschaft für Wärme und Heizkraftwirtschaft: Technisches Handbuch Fernwärme Technik und Normung (Frankfurt am Main: AGFW-Projektges. für Rationalisierung, Information und Standardisierung).

Ahrens, B., Duda, M., and Renner, J. (2018). Relations between Hydraulic Properties and Ultrasonic Velocities during Brittle Failure of a Low-Porosity sandstone in Laboratory Experiments. Geophys. J. Int. 212, 627–645. doi:10.1093/gji/ggx419