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).

Typically, Ruhr sandstone consists mainly of quartz ca. 70 ± 5   %, feldspar ca. 10 ± 2   %, and clay and feldspar alterations ca. 15 ± 2   %, and minor constituents of iron oxides ( < 5 ± 2   %) (e.g., Kneuper, 1957; Rosenfeld, 1967; Duda and Renner, 2013; Stöckhert, 2015). The studied specimen consists mainly of quartz 70 ± 5   %, minor amount of feldspar phases like albite or orthoclase < 5 ± 2   %, and 20 ± 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 ± 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).

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.

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 × 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 × 2,940 × 2,141 voxel and shows a physical volume of 5.88 × 5.88 × 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.

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 × 2 × 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 × 4 × 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.

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.

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).

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 ⋅ 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 × c 11 ,  qtz and c 44 , watershed = 0.5 × c 44 ,  qtz . This factor was determined within the numerical studies as a result of a first calibration approach.

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 ± 1, 2,632 ± 1, and 2,647 ± 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 ± 1.87, 5.13 ± 1.90, and 5.37 ± 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 ± 0.78, 5.45 ± 0.77, and 5.10 ± 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.

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).

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)].

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 = [ 5.9710 0.00376 0.00519 0.00343 5.5045 0.00478 0.00294 0.00537 5.4968 ] ⋅ 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).

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.

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 × 2 × 2 and 4 × 4 × 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.

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.

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).

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.