Characterization of discontinuities in potential reservoir rocks for geothermal applications in the Rhine-Ruhr metropolitan area (Germany)

The importance of research into clean and renewable energy solutions has increased over the last decade. Geothermal energy provision is proven to meet both conditions. Therefore, conceptual models for deep geothermal applications were developed for different field sites regarding different local conditions. In Bavaria, Germany, geothermal applications were successfully carried out in carbonate horizons at depths of 4000 to 6000 m. Matrix permeability and thermal conductivity was mainly studied in karstified carbonates from the Late Jurassic reef facies. Similar to Bavaria, carbonates are located in the east of the Rhenohercynian Massif, in North Rhine-Westphalia (NRW), for which quantification of the geothermal potential is still lacking. Compared to Bavaria, a supraregional carbonate mountain belt is exposed at the Remscheid-Altena anticline (in NRW) from the Upper Devonian and Lower Carboniferous times. The aim of our study was to examine the potential geothermal reservoir by field and laboratory investigations. Therefore, three representative outcrops in Wuppertal, Hagen-Hohenlimburg, and Hönnetal were studied. During field surveys, 1068 discontinuities (139 open fractures without any filling, 213 joints, 413 veins filled with calcite, and 303 fractures filled with debris deposits) at various spatial scales were observed by scanline surveys. These discontinuities were characterized by trace length, true spacing, roughness, aperture, and filling materials. Discontinuity orientation analysis indicated three dominant strike orientations in NNW–SSE, NW–SE, and NE–SW directions within the target horizon of interest. This compacted limestone layer (Massenkalk) is approximately 150 m thick and located at 4000 to 6000 m depth, dipping northwards at a dip angle of about 30 to 40^∘. An extrapolation of the measured layer orientation and dip suggests that the carbonate reservoir could hypothetically extend below Essen, Bochum, and Dortmund. Our combined analysis of the field and laboratory results has shown that it could be a naturally fractured carbonate reservoir. We evaluated the potential discontinuity network in the reservoir and its orientation with respect to the prevailing maximum horizontal stress before concluding with implications for fluid flow: we proposed focusing on prominent discontinuities striking NNW–SSE for upcoming geothermal applications, as these (1) are the most common, (2) strike in the direction of the main horizontal stress, (3) have a discontinuity permeability that significantly exceeds that of the reservoir rock matrix, and (4) only about 38 % of these discontinuities were observed with a calcite filling. The remaining discontinuities either showed no filling material or showed debris deposits, which we interpret as open at reservoir depth. Our results indicate that even higher permeability can be expected for karstified formations related to the reef facies and hydrothermal processes. Our compiled data set, consisting of laboratory and field measurements, may provide a good basis for 3D subsurface modelling and numerical prediction of fluid flow in the naturally fractured carbonate reservoir.

1 Introduction

The Rhine-Ruhr metropolitan region in western Germany is one of the largest metropolitan areas in Europe. It has the largest European district heating network (AGFW, 2009) and is one of the largest energy systems in the world that is predominantly fed by fossil fuels (Klaus et al., 2010). In recent decades, the Ruhr metropolitan area has been subject to structural changes associated with the ending of traditional coal mining and steel industry, as well as accelerated urbanization. Considering effective climate change, these processes pose local to global challenges for a post-fossil energy future. The polycentric structure of the Rhine-Ruhr area is particularly promising for the requirements imposed by the energy transition (Wegener et al., 2019). To account for the increasing demand for energy supply (e.g., Scheer et al., 2013; Araújo, 2014), the successful implementation of sustainable non-fossil energy concepts, such as geothermal energy projects, is required (Fridleifsson et al., 2008; Goldstein et al., 2011).

In recent years, many economically successful geothermal energy projects in Bavaria, southern Germany, have been realized and could serve as a valuable model for the Rhine-Ruhr area. The Bavarian part of the South German Molasse Basin has a high hydrothermal energy potential (Schellschmidt et al., 2010; Böhm et al., 2010; Stober, 2014; Homuth, 2014), which can be extracted from pre-existing discontinuity and karst systems in deep carbonate rocks (e.g., limestone deposits Fritzer et al., 2010). Similar to deep Bavarian deposits, carbonate rocks were accumulated in western Germany. Geothermal systems aiming for comparable geological carbonate horizons in the Rhine-Ruhr area, located at depths relevant for deep geothermal applications (DEKORP Research Group, 1990), are promising solutions to initiate the energy transition (Knutzen, 2017). A feasibility study including extensive subsurface reservoir characterization and local geological site assessment is essential to describe the hydrothermal potential of presumably fractured and faulted carbonates in the Rhine-Ruhr area.

The estimation and description of individual discontinuities and entire discontinuity networks in deep reservoirs is of elementary importance when it comes to reservoir characterization. Discontinuities and discontinuity networks have a decisive influence on the fluid flow within the reservoir (e.g. Odling et al., 1999) and also on its stability (e.g. Cappa et al., 2005). They are typical observations at analogue outcrops and can be measured there (see review by Bonnet et al., 2001). While the origin of these typically tectonically driven discontinuities can be discussed extensively, the influence on reservoir permeability is obvious. Compared to the matrix permeability of the host rock, the permeability of discontinuities is generally several orders of magnitude higher (e.g. Nelson and Handin, 1977; Kranz et al., 1979; Evans et al., 1997; Ahrens et al., 2018). Detailed discontinuity analysis is essential to model reservoir quality, optimal well location, and well performance (Nelson, 2001). Field surveys are both the first and one of the most crucial steps in the exploitation of the reservoir (van Golf-Racht, 1982; Agosta et al., 2010), i.e. the investigation of the discrete fracture network (DFN). Furthermore, since direct information on in situ properties is limited and inaccessible (e.g. seismic tomography or well and borehole data), outcrop characterizations and laboratory measurements can provide the basis for reservoir characterization and the estimation of the potential reservoir behaviour under in situ conditions.

In this integrated study, we examine the geothermal potential of deep carbonates in the Rhine-Ruhr area by combining field and laboratory investigations of pre-existing discontinuities and discontinuity networks at outcrop and sample scales. The combined analysis of field and laboratory measurements allows us to describe the potential carbonate reservoir by characterizing its natural fracture systems. Three outcrop analogues on the eastern side of the northern Rhenohercynian Massif have been chosen for field survey and sample collection. All quarries are located within the Devonian Reef Complex. The dominant stratigraphic unit in those quarries is the Massenkalk facies from the Middle and Upper Devonian (Krebs, 1970; Paeckelmann, 1979; Schudack, 1993). Dominant discontinuity systems and different facies have been documented and sampled. Most of these measurements have been conducted by one-dimensional scanline surveys (Priest and Hudson, 1976, 1981). Representative fresh rock samples have been taken directly from the outcrop wall and petrophysically characterized by different laboratory methods.

In the following, we first introduce the geology of the Devonian Reef Complex. We proceed by describing the outcrops and the employed field and laboratory methods. Following the presentation of the measurements results, we conclude by discussing the geothermal potential of the geological subsurface model, i.e. the deep carbonate layers and their discontinuity systems.

2 Geology, outcrop investigations, and laboratory measurements

2.1 Geological setting: the Devonian Reef Complex

The Rhenohercynian Massif is the northernmost mountain belt of Europe, which shows the results of the Hercynian orogeny with its numerous northeast striking folds and thrusts at the northern edge. During the early Devonian a passive continental margin developed south of Laurussia, which was affected by crustal thinning. A shelf sea formed above the thin crust, which is called the Rhenish-Hercynian Basin. The shallow sea, which ranged from Brittany across Belgium to Germany, was filled with clastic shelf deposits and carbonates from Laurussia. During the Early Devonian, the sedimentation in the offshore belt decreased, and the horizons are characterized by a transition from siliceous, deltaic shallow marine environments (Rhenish facies) to open marine shale and purely limestone-dominated facies (Hercynian facies) (Krebs, 1970; Grabert, 1998; Pas et al., 2013; Meschede, 2018). In the late Middle Devonian, the continued northwestward movement favoured the growth of submarine volcanic mounds and thus also that of carbonates and coral and stromatoporoid reefs on the shelf (Jux, 1960; Grabert, 1998; Dallmeyer et al., 2013; Franke et al., 2017). Along the coast of Laurussia an elongated reef complex developed, the so-called Devonian Reef Complex (Krebs, 1970). In the Late Devonian, the first coral reefs disappeared in the course of the Late Devonian extinction (Kellwasser event), which was first favoured by a marine transgression and then by rapid regression (e.g. Meschede, 2018).

Figure 1 (a) Simplified geological map of the Rhine-Ruhr metropolitan area, Germany (modified from Jansen et al., 1986; Drozdzewski et al., 2007). Here, three stone pits have been chosen for field surveys: Osterholz (Wuppertal), Oege (Hagen), and Asbeck (Hönnetal), which are part of the Devonian Reef Complex (Massenkalk). The red lines indicate the three geological cross sections as shown in Fig. 2. In panel (b), the geological map from panel (a) was further simplified, and in addition the rose diagrams of the discontinuity sets as measured by scanline investigations at the quarries in Wuppertal (blue), Hagen-Hohenlimburg (green), and Hönnetal (orange) are added (for more detail, see Figs. 4 to 6). The arrows indicate the direction of maximum horizontal stress σ_H according to the World Stress Map (Heidbach et al., 2016). The specific directions and angles of incidence of the investigated, exposed limestone layers are also shown in panel (b).

Figure 2 Simplified cross sections of the Remscheid-Altena anticline (modified after Jansen et al., 1986; Drozdzewski et al., 2007). The exact position and orientation of the three cross sections is illustrated in Fig. 1a.

The Devonian Reef Complex extends over the area from Neandertal (Düsseldorf) across Dornap, Wuppertal, Schwelm, Hagen-Hohenlimburg, Letmathe, Sundwig, Iserlohn, Hemer, and Balve (Fig. 1a). Furthermore, these horizons can also be observed from drilling in northern Germany (e.g. Fig. 2; Hesemann, 1965, Münsterland 1). During the Early Carboniferous, the sedimentation changed during the closure of the Rhenohercynian basin. Primarily caused by tipping blocks and faults, first greywacke turbidites and second a quartz-rich sandstone was deposited in the foreland. Marine regressions occurred due to the mentioned sediment deposits. This was followed by a gradual subsidence of the bedrock. These periodic marine transgressions (i.e. sea level rises) and regressions (i.e. sea level falls) can be followed in the Upper Carboniferous Ruhr sequences (cyclothems) in today's Rhine-Ruhr area. From the Lower Carboniferous onwards, the deposited sediments were compacted and folded by the Hercynian orogeny. In the meantime, the thick sediment accumulations (3 to 12 km) were thrusted by reactivated listric normal faults during the NW movement (Oncken, 1997; Franke et al., 1990; Engel et al., 1983; Franke and Engel, 1982; Holder and Leversidge, 1986; Meschede, 2018). Between the Early Triassic and the Late Cretaceous at least two further extension sequences developed, which caused a NE–SW-directed shortening (Brix et al., 1988; Drozdzewski and Wrede, 1994). Since the Eocene, the complex tectonic environment is affected by another extension regime (Kley, 2013).

2.2 Study area

Three stone pits have been chosen for field studies within the Devonian Reef Complex: (1) Osterholz stone pit, located in Wuppertal; (2) Oege quarry, located in Hagen-Hohenlimburg; and (3) Asbeck quarry, located in Hönnetal (Fig. 1a). The outcrops in Osterholz, Oege, and Asbeck are active open-pit mines and are driven by the following independent companies: Kalkwerke H. Oetelshofen GmbH & Co. KG, Hohenlimburger Kalkwerke GmbH, and Rheinkalk GmbH, respectively. The entrances of the three quarries are located at (1) 51.248219^∘ N, 7.042342^∘ E, (2) 51.347811^∘ N, 7.581534^∘ E, and (3) 51.400452^∘ N, 7.842271^∘ E (Fig. 3). In the following, the outcrops will be referred according to their geographical location: Wuppertal (WOH), Hagen-Hohenlimburg (HKW), and Hönnetal (HLO).

Figure 3 Map view of all three study areas in which scanline surveys were carried out: (a) the Osterholz stone pit, which is operated by Kalkwerke H. Oetelshofen GmbH & Co. KG in Wuppertal; (b) the Oege stone pit, operated by Hohenlimburger Kalkwerke GmbH in Hagen-Hohenlimburg; (c) the Asbeck rock mine, managed by Lhoist Deutschland, Rheinkalk GmbH in Hönnetal. Red lines indicate the predominant rock walls, i.e. the direction of the mining. In each of the three outcrops there were two dominant mining directions, which were approximately perpendicular to each other. This figure was created using ArcGIS^® software by Esri (Basemap, World Imagery: http://www.arcgis.com/home/item.html?id=10df2279f9684e4a9f6a7f08febac2a9, last access: 24 March 2020). © ArcGIS and ArcMap™ are the intellectual property of Esri and are used herein under license. © Esri. All rights reserved. For more information about Esri™ software, please visit http://www.esri.com (last access: 24 March 2020).

The Osterholz quarry in Wuppertal is located in the Osterholz graben, which is limited by the Herzkamper syncline to the north and the Wupper river to the south (Fig. 1a). The regional geology is dominated by WSW–ENE-striking folds, which are associated with the Hercynian orogeny (Paeckelmann, 1979). This local tectonic regime was mainly affected by the Hercynian orogeny and became further compressed by recent folding processes in the south. The stratigraphic unit of the quarry Osterholz in Wuppertal is Upper Givetian (approx. 382 Ma), primarily the Schwelm facies (according to Krebs, 1970). The studied limestone layers, with a mean thickness of 1 to 5 m, showed greyish, compacted layers of well-bedded carbonates with corals, stromatoporoids, and bioclastic materials. This approximately 150 m thick horizon sits on top of dolomitic carbonates and is related to the brachiopoda- and coral-rich series or stromatoporoid series. The limestone bed boundaries show mechanical and fracture stratigraphy.

The regional tectonic structure is characterized by the Remscheider anticline, the Ennepe thrust, and the Großholthausener fault in Hagen-Hohenlimburg (Fig. 1a). In the Remscheider anticline's northern limb, where the studied quarry was set, the Hasper anticline and Voerder syncline are located. Both folds are monoclines, and their steep limbs rapidly end when intersecting the regional strike–slip fault, the Großholthausener fault. The Oege quarry in Hagen-Hohenlimburg is based in the eastern part of the Großholthausener fault. Krebs (1970), Flügel and Hötzl (1976), and Koch-Früchtel and Früchtel (1993) characterized the local compacted limestone in the quarry Oege as Schwelm facies. The quarry is enriched in dolomite due to hydrothermal veins (von Kamp and Ribbert, 2005). These dolomitic carbonates manifest in a high degree of fractured blocks and karst formations, which were filled by Tertiary and Cretaceous sediments, respectively. The carbonates, which were untouched by the hydrothermal veins, contain brachiopoda, bivalves, crinoids, corals, gastropoda, stromatoporoids and micrite. These carbonates are identified as Middle Devonian (Givetian according to Basse et al., 2016).

The Asbeck quarry in Hönnetal is located to the north of the Remscheid-Altena anticline (Fig. 1a). The exposed compacted limestone deposits dip gently to the north and are surrounded by Lower Devonian sediments in the south and Upper Devonian to Lower Carboniferous sediments to the north. The regional setting is dominated by NNW–SSE-striking strike–slip faults, which are associated with the folding mechanics during Hercynian orogeny like the Großholthausener fault in Hagen-Hohenlimburg (Paeckelmann, 1938). Those strike–slip faults can be tracked along with the major Remscheid-Altena anticline structure. At the most eastern part of the Remscheid-Altena anticline, the fold's tip plunges in the NNE direction. The local stratigraphy of the Asbeck quarry in Hönnetal consists of an initial bank stadium (Schwelm facies) that developed into a true bioherm complex, i.e. Dorp facies (Becker et al., 2016). Schudack (1993) has investigated the different facies and cyclic depositions of the reef and back-reef lagoon facies and subdivided these facies into six subfacies, which were investigated in more detail by Kloke (2007). The greyish, highly compacted limestones (approx. 600 m thickness) containing corals, stromatoporoids, and brachiopods are associated with late Middle Devonian (Givetian) or early Late Devonian (Frasnian). Except for some minor parts that were mentioned by Kloke (2007), no dolomitic rocks were documented.

The mean bedding orientation in the quarries Wuppertal, Hagen-Hohenlimburg, and Hönnetal is approximately 327, 345, and 18^∘ in dip direction, with an average dip angle of about 37, 42, and 28^∘. Each stone pit's working levels were approximately oriented perpendicular to each other and strike NNW–SSE and NE–SW (Fig. 3).

2.3 Field method and theoretical estimation: 1D outcrop scanline surveys and 2D discontinuity connectivity

The term discontinuity and its origin can be looked up in many textbooks on structural and engineering geology (e.g. Priest, 1995). In the following, we define a discontinuity as a weak zone within a rock mass that can be referred to as a fracture, joint, vein, or fault. Unlike the International Society for Rock Mechanics (ISRM, 1978), we neglect weak bedding planes or schistosity planes in the field surveys. Nevertheless, we distinguish between natural discontinuities, which are associated with geological origins, and artificial discontinuities, which result from anthropogenic processes (e.g. drilling, blasting, or excavation).

Table 1 Overview of the documented discontinuity properties and their classification as defined for scanline surveys (modified from Markovaara-Koivisto and Laine, 2012). In our observations we distinguish between open, closed, and filled discontinuities. Open discontinuities are those that are neither filled with calcite nor with debris. Closed discontinuities are in fact joints. Filled discontinuities are not categorized as closed ones, since both discontinuity surfaces are not in direct contact with each other. Filled discontinuities can either be filled with calcite or with debris (i.e. veins). The roughness of the discontinuities was classified on field scale (> 10⁻¹ m) and on a mesoscopic scale (⩽ 10⁻¹ m).

The discontinuities were recorded by a scanline sampling technique according to Priest and Hudson (1981). The advantages and disadvantages of this method have already been discussed in detail in several studies over the last few decades (e.g. Attewell and Farmer, 1976; Priest and Hudson, 1981). Based on our field observations, we simply assume that spatially homogeneous discontinuity networks dominate in the exposures, depending on the spatial observation scale. In regions with spatially heterogeneous discontinuity networks, however, more complex methods than the scanline sampling technique are required for discontinuity characterization (e.g. Watkins et al., 2015). Following the study of Priest and Hudson (1981), our paper focuses only on the most important parts of this technique. The setting of a scanline survey is a tape line, which is fixed firmly to the exposed rock wall and shows a start and an endpoint. The length of a scanline can range from a few decimetres to several metres. Estimates of all possible error sources have shown that the number of measured and recorded discontinuities is more significant than the length of the scanline (Priest and Hudson, 1981). Therefore, we selected representative, safe rock walls for each quarry where at least 30 measurable discontinuities (number of discontinuity observations n ≥ 30) could be recorded. After the tape line was fixed, the trend and plunge angles of the tape line were measured. The trend angle reflects the orientation of the tape line (${\mathit{\alpha }}_{\mathrm{S}}=\mathrm{0}{}^{\circ }\equiv$ north; ${\mathit{\alpha }}_{\mathrm{S}}=\mathrm{180}{}^{\circ }\equiv$ south), which corresponds to the strike of the rock wall. The tape's line levelling is defined by the angle of plunge (approx. $\mathrm{0}{}^{\circ }\le {\mathit{\beta }}_{\mathrm{S}}\le \mathrm{15}{}^{\circ }$). After the tape line was installed, each discontinuity that intersected the tape line was recorded and categorized according to the following six properties (see also Table 1):

1. Intersecting distance. The intersecting distance corresponds to the distance between each intersecting discontinuity and the fixed starting point of the tape line.

2. Discontinuity length. The installed tape line and the intersecting discontinuity meet at an intersecting point. The trace length above and below the tape towards the discontinuity tips indicates the length of the discontinuity. If both tips exceed the rock wall height, the discontinuity is recorded as “through-going”. Accordingly, discontinuities of which one or both ends are visible are categorized as “one end visible” or “both ends visible”. On rare occasions, neither end can be recorded, and thus the discontinuity is documented as “neither end visible”.

3. Orientation. The dip and dip direction of the intersecting discontinuities were recorded.

4. Roughness. Roughness is the inherent alignment and waviness of the surface at the mean level of a discontinuity. A large-area waviness may also alter the dip locally. All roughness measurements were recorded manually. Therefore, we distinguish between “smooth”, “slightly rough”, and “rough”, both on mesoscopic and field scales. On the mesoscopic scale, the discontinuity surfaces were analysed for mineral steps, stylolites, or plumose structures. On the field scale, the discontinuity roughness was determined by subjective (qualitative) and visual evaluation of wavelength measurements with a tape line. The specified roughness classes correspond to the average roughness on the respective scale.

5. Aperture. Aperture is the perpendicular distance between adjacent rock walls of a discontinuity where the space between them is filled with a material, i.e. filling. The determined values correspond to the mean apertures determined by geometric (qualitative) analysis of the respective total discontinuity lengths.

6. Filling. Filling is the material that separates the adjacent rock wall of a discontinuity and is usually weaker than the parent rock. In our study, typical filling materials are calcite, clay, or debris. These include thin mineral coatings and healed discontinuities such as quartz and calcite veins.

In spite of the above-mentioned records, it must be emphasized again how important the calculation of the intersecting angle is. Due to an irregular outcrop wall, the distance between the starting point and the recorded discontinuity (i.e. apparent spacing) is measured incorrectly. This apparent distance distorts the true orientation of the discontinuity. The true angle δ between the scanline orientation (trend and plunge) and the normal to recorded discontinuity was calculated according to

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathit{\delta }& =\arccos \left(\cos \left({\mathit{\alpha }}_{\mathrm{n}}-{\mathit{\alpha }}_{\mathrm{s}}\right)\cos \left({\mathit{\beta }}_{\mathrm{n}}\right)\cos \left({\mathit{\beta }}_{\mathrm{s}}\right)\right.\\ & \left.+\sin \left({\mathit{\beta }}_{\mathrm{n}}\right)\sin \left({\mathit{\beta }}_{\mathrm{s}}\right)\right),\end{array}\end{array}$$

where α_s, β_s, α_n, and β_n denote the trend direction of the tape line (scanline), the plunge direction of the tape line with regard to trend orientation, and the dip direction and dip of the normal to the discontinuity, respectively. Further, the true angle δ and the apparent spacing χ_s between two discontinuities were used to calculate the true spacing assuming a homogeneous discontinuity pattern χ_r=cos (δ)χ_s. By evaluating the ratio of true discontinuity spacing χ_r and the number of discontinuity observations, the mean discontinuity spacing $\stackrel{\mathrm{‾}}{\mathit{\chi }}={\mathit{\lambda }}^{-\mathrm{1}}=(\sum {\mathit{\chi }}_{\mathrm{r}})n^{-\mathrm{1}}$ was determined, where λ denotes the discontinuity frequency. The evaluation of the mean discontinuity track length with respect to the required accuracy has already been discussed in detail by Priest and Hudson (1981). The authors have made various suggestions when the end of a discontinuity trace on an investigated rock wall may not be visible, for example, due to geological events or anthropogenic processes. In our study, the proposed solution for censored semi-trace length sampling was adopted to evaluate the mean discontinuity track length μ_i. If the distribution of trace length over the recorded rock wall follows a power-law distribution, then the frequency f(l) is given by $f\left(l\right)=\mathit{\lambda }{\mathrm{e}}^{-\mathit{\lambda }\stackrel{\mathrm{‾}}{\mathit{\chi }}}$. Following the approach of Priest and Hudson (1981), for censored semi-trace length analysis, the mean discontinuity trace length was determined according to

$$\begin{array}{}\text{(2)}& {{\mathit{\mu }}_i}^{-\mathrm{1}}={\mathit{\mu }}^{-\mathrm{1}}-\left(c{\mathrm{e}}^{-\mathit{\mu }c}\right){\left(\mathrm{1}-{\mathrm{e}}^{-\mathit{\mu }c}\right)}^{-\mathrm{1}},\end{array}$$

where μ and c denote mean population frequency and concealed trace length of the recorded discontinuity, respectively. Typically, the concealed trace length is an order of magnitude smaller than the recorded trace length. Note that discontinuities exceeding the height of the outcrop wall could be detected at the nearest level above and/or below the quarry.

In addition, the simplified 2D discontinuity connectivity of discontinuities aligned perpendicular to the examined outcrop face, i.e. the average number of discontinuity intersections per discontinuity P, was estimated for each opencast mine according to Ozkaya (2011)

$$\begin{array}{}\text{(3)}& P=\sum _k^n{d}_k^{-\mathrm{1}}\sum _k^n{d}_k\sum _{j\ne k}^n{d}_j{L}_j\sin \left(\mathrm{\Delta }{\mathit{\theta }}_{jk}\right),\end{array}$$

where Δθ_jk, d_j, and L_j denote the angle between the average striking directions of two discontinuity sets j and k, the number of discontinuities per total scanline length (i.e. discontinuity density), and the mean discontinuity length of each discontinuity set, respectively. In each outcrop we encountered three sets of discontinuities (n=3) in the respective investigated carbonate layers. Therefore, the average discontinuity connectivity P of each exposed carbonate layer was estimated on the basis of three average discontinuity densities and mean discontinuity lengths.

The recorded discontinuities were evaluated with MATLAB (2018) following Hudson (2005) and Markovaara-Koivisto and Laine (2012). Markovaara-Koivisto and Laine (2012) developed an open-source MATLAB code for the visualization of scanline survey results that was adapted to the purpose of this study. All measurements were recorded according to the metric system of measurement.

2.4 Laboratory measurements: petrophysical characterization of samples

During the fieldwork, loose rock blocks were taken in each quarry for further petrophysical characterization. Compacted limestones, dolomites, and red dolomitic rocks were named MKB (black Massenkalk), MKY (yellow Massenkalk), and MKR (red Massenkalk), respectively. After sampling, the rocks were named according to the sampled quarry, sampled working level, rock type, and sample number. Representative cylindrical cores with a diameter of 40 mm were extracted from each of these blocks by diamond core drilling perpendicular to the bedding direction. All samples were saw-cut perpendicular to the core axis and their end faces were ground square to the maximal possible length l. All steps of the preparation were conducted with water as coolant and for rinsing removed material. Following preparation, samples were oven-dried at 60 ^∘C for about 48 h. Basic petrophysical properties were determined on identically prepared samples at ambient conditions, except permeability, which was derived under elevated pressures.

Bulk density ρ_geo was calculated from the geometrical volume of the cylindrical samples and their dry masses. Grain density ρ_grain was gained from pycnometer measurements on rock powder, produced by crushing and grinding of leftover rock fragments, in compliance with the German (NABau, 2011) standard. By evaluating the ratio between bulk and grain density, the total porosity was determined according to ${\mathit{\varphi }}_{\mathrm{tot}}=\mathrm{1}-({\mathit{\rho }}_{\mathrm{geo}}/{\mathit{\rho }}_{\mathrm{grain}})$. The connected porosity, i.e. the externally accessible and connected pore volume, was determined using the difference of the masses of dry and saturated samples with distilled water (see Duda and Renner, 2013).

Ultrasonic P- and S-wave velocities, v_P and v_S, were determined on dry and saturated samples from first-arrival measurements using an ultrasound benchtop unit composed of a waveform generator, two identical broadband ultrasound sensors (1 MHz centre frequency and 0.5 in., i.e. 12.7 mm, diameter), and a digital storage oscilloscope (200 MHz sampling frequency). Measurements were performed parallel to the cylinder axis, i.e. the drilling direction. Velocities were calculated by dividing the sample length by the determined arrival times, less the assembly parts' travel time. Drained dynamic Young's modulus ${\stackrel{\mathrm{̃}}{E}}_{\mathrm{d}}$ and Poisson's ratio ${\stackrel{\mathrm{̃}}{\mathit{\nu }}}_{\mathrm{d}}$ were calculated by dry P- and S-wave velocities, v_P,dry and v_S,dry, and bulk density assuming isotropy (see Mavko et al., 2020). Moreover, saturated P- and S-wave velocity, v_P,sat and v_S,sat, and the density of the fluid-saturated sample was used to calculate undrained Poisson's ratio ${\stackrel{\mathrm{̃}}{\mathit{\nu }}}_{\mathrm{ud}}$ and Young's modulus ${\stackrel{\mathrm{̃}}{E}}_{\mathrm{ud}}$, assuming isotropy and employing Gassmann's hypothesis (Gassmann, 1951) that S-wave velocity remains unaffected by the presence of a fluid ($v_{\mathrm{S},\mathrm{dry}}=v_{\mathrm{S},\mathrm{sat}}$).

Thermal conductivity κ_dry of the dry sample cores was determined by a thermal conductivity scanner using an optical scanning method at ambient conditions (Popov, 1997). The scanner consists of an emitter and a measuring unit that is moved along the sample at a fixed distance. The emitted light and heat radiation are focused on the sample's surface, which heats up the sample pointwise. To ensure absolute absorption of the energy, a part of the sample, usually a strip, was painted black. Furthermore, it was ensured that the samples for thermal conductivity measurements meet the geometrical requirements for sample size in order to reduce boundary effects. Infrared temperature sensors are located at a fixed distance from the emitter (lead sulfide infrared receiver) and measure the sample's temperature difference before and after heating. The thermal conductivity was determined by comparison with known standards. In this study, we used the standards of quartz (κ_qtz=1.35 W (m K)⁻¹) and titanium alloy (κ_Ti=6.05 W (m K)⁻¹). According to the manufacturer (Lippmann and Rauen GbR; TCS no. 2010-013), the determination of thermal conductivity is subject to an absolute error of approximately 3 %.

Permeability k of the samples was determined using Darcy's law (Darcy, 1856) and a conventional Hoek cell (Hoek and Brown, 1997) to apply the axial load, pore-fluid pressure, and confining pressure. Axial pistons and the pre-saturated samples were jacketed by a rubber tube to prevent oil from penetrating the sample, i.e. a connection between confining and pore pressure. In addition, confining and pore-fluid pressure were kept below the axial pressure. A hand pump was operated, supplying the axial load (12 MPa). A computer-controlled, high-pressure metering pump was used to apply confining pressure on the sample by compressing distilled water (10 MPa). The axial pore-fluid flow was ensured through central bores in the axial loading pistons. Distilled water was pumped from a water reservoir through the lower end face of the samples either by applying constant flow rates (0.001 to 0.15 L h⁻¹). The lower axial loading piston was equipped with an outlet pipe discharging fluid pressure to the atmosphere (i.e. 100 kPa). The temperature and pressure of the pore fluid were measured to calculate the temperature- and pressure-dependent fluid viscosity, according to Wagner (2009). Each permeability of the samples corresponds to the arithmetic mean value of all experiments conducted with various flow rates for steady-state fluid flow conditions.

3 Results

3.1 Field discontinuity observations and discontinuity connectivities

A total of 1068 discontinuity observations were recorded and classified (e.g. discontinuity type, filling, and roughness) by using the field surveys in the three outcrops in Wuppertal, Hagen-Hohenlimburg, and Hönnetal (Table 2). To map all the features of the outcrop, we decided to carry out surveys at different levels of the quarry. During the field observations in Wuppertal, scanline surveys were taken in two orientations: NNW–SSE and NE–SW (Fig. 3a). In contrast, the scanline investigations in Hagen-Hohenlimburg and Hönnetal were conducted on two exposed rock faces, each with the orientations NNW–SSE and ENE–WSW (Fig. 3b, c). Comparable scanline surveys over total lengths of 66.60, 76.15, and 62.76 m were measured in Wuppertal (461 discontinuities), Hagen-Hohenlimburg (361 discontinuities), and Hönnetal (246 discontinuities), respectively. It was found that in all quarries the majority of the discontinuities were filled (Table 2). In Wuppertal most of the discontinuities were filled with debris, whereas in the other two outcrops calcite was the dominant filling material. In general, most of the ends of all discontinuities were visible on the outcrop wall. On the field scale the discontinuities had a smooth roughness and on a smaller (mesoscopic) scale they were slightly rough. The mean discontinuity trace lengths clearly differed between the outcrops, so that the largest one was determined in Wuppertal with 6.65 m and the smallest in Hagen-Hohenlimburg with 1.23 m. The mean discontinuity spacing of the quarries were all in the same order of magnitude, amounting to 0.15, 0.24, and 0.27 m in Wuppertal, Hagen-Hohenlimburg, and Hönnetal, respectively.

Table 2 Overview of the recorded properties and their classification of the measured discontinuities in the three studied outcrops in Wuppertal, Hagen-Hohenlimburg, and Hönnetal. See the text and Table 1 for a detailed description of the classes.

Table 3 Overview of the settings and the results of the scanline surveys performed on rock walls in the outcrops Wuppertal (WOH), Hagen-Hohenlimburg (HKW), and Hönnetal (HLO). The abbreviations MKY, MKB, and MKR correspond to the coloured rock types Massenkalk yellow (dolomite), Massenkalk black (limestone), and Massenkalk red (dolomite), respectively. Three main families or sets of discontinuity orientations were found to be predominant in all outcrops.

l_scan: scanline length; n: number of discontinuities; ${\stackrel{\mathrm{‾}}{l}}_{\mathrm{dis}}$: mean trace length; ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{dis}}$: mean discontinuity aperture; ${\stackrel{\mathrm{‾}}{d}}_{\mathrm{dis}}$: mean discontinuity spacing. The resolution of the length measurements amounts to 0.01 m.

Three sets of discontinuity orientations were predominant in all outcrops (Table 3). In Wuppertal the sets were grouped into the arithmetic mean discontinuity azimuth directions of strike 172^∘ (181 discontinuities), 55^∘ (142 discontinuities), and 105^∘ (138 discontinuities) (Fig. 4c, d; Table 6). In contrast, the discontinuity sets in Hagen-Hohenlimburg and Hönnetal were oriented towards 135^∘ (134 discontinuities), 12^∘ (95 discontinuities), and 75^∘ (132 discontinuities) and 176^∘ (140 discontinuities), 65^∘ (24 discontinuities), and 87^∘ (82 discontinuities), respectively (Figs. 5c, d; 6c, d; Table 6). The dip angles of the discontinuities were documented as very steep, i.e. between 80 and 90^∘ (dip) in all three quarries (Table 3). By further differentiating the strike directions of the classified discontinuity sets according to their filling material, no strike rotation between calcite- and debris-filled veins and fractures was identified (Fig. 7). Besides, several generations of repeated cracking and sealing were observed that were located within the recorded veins in Hagen-Hohenlimburg. A precise analysis of the individual vein generations is currently the subject of recent studies and could not be provided within the scope of this study.

Figure 4 Results of all scanline investigations at the quarry in Wuppertal: (a) stereogram of measured discontinuity sets and a (b) rose diagram of rotated discontinuity orientations according to the bedding layer in panel (a) for a horizontal position, (c) discontinuity spacing diagram, and (d) discontinuity trace length diagram.

Figure 5 Results of all scanline investigations at the quarry in Hagen-Hohenlimburg: (a) stereogram of measured discontinuity sets and a (b) rose diagram of rotated discontinuity orientations according to the bedding layer in panel (a) for a horizontal position, (c) discontinuity spacing diagram, and (d) discontinuity trace length diagram.

Figure 6 Results of all scanline investigations at the quarry in Hönnetal: (a) stereogram of measured discontinuity sets and a (b) rose diagram of rotated discontinuity orientations according to the bedding layer in panel (a) to a horizontal position, (c) discontinuity spacing diagram, and (d) discontinuity trace length diagram.

Figure 7 Rose diagrams of measured discontinuity sets at the quarries (a, b) Wuppertal, (c, d) Hagen-Hohenlimburg, and (e, f) Hönnetal as a function of their discontinuity filling, i.e. whether the discontinuities are filled with (a, c, e) calcite or (b, d, f) debris. The discontinuity strike azimuths given next to the figures correspond to the arithmetic mean values and the corresponding standard deviations.

The statistical analysis of the three subdivided sets from Wuppertal and Hagen-Hohenlimburg showed that although the mean discontinuity lengths differed significantly, the discontinuity densities were comparable (Table 6). In Wuppertal and Hagen-Hohenlimburg the discontinuity densities ranged from about 2.00 to 2.75 and 1.25 to 1.75 m⁻¹, respectively. Only in Hönnetal were clear differences in both mean discontinuity lengths and discontinuity densities found between the three sets of discontinuities. Here, the mean discontinuity length and the discontinuity densities amounted to 4.51, 5.57, and 2.85 m and 2.23, 0.38, and 1.31 m⁻¹, respectively. The identified 2D discontinuity connectivity was about 6 to 8 times higher in Wuppertal (P=21.42) than in the other outcrops at Hagen-Hohenlimburg (P=3.71) and Hönnetal (P=2.49).

3.2 Laboratory characterizations

Samples could be obtained from all outcrops of dolomitic carbonates and limestone rocks, although it was challenging to obtain samples from the Hönnetal due to the lack of revealed outcrops. Most of the exposed Hönnetal limestone rocks were purely calcitic and showed a low grade of dolomitization. Nevertheless, some dolomitic carbonates have been sampled. Dolomitic carbonate samples from Hagen-Hohenlimburg occurred as highly fractured rocks with karst and highly porous formations, ranging from a few millimetres to several centimetres.

Table 4 Matrix properties, permeabilities, and thermal conductivities as determined by laboratory measurements on limestone and dolomite samples from Wuppertal (WOH), Hagen-Hohenlimburg (HKW), and Hönnetal (HLO). The abbreviations MKY, MKB, and MKR indicate the coloured rock types Massenkalk yellow (dolomite), Massenkalk black (limestone), and Massenkalk red (dolomite), respectively.

l: sample length; ρ_geo: bulk density; ρ_grain: grain density; ϕ_tot: total porosity; ϕ_con: connected porosity; k: permeability; κ_dry: thermal conductivity. Quoted uncertainties reflect accuracy of the measurements.

The mean bulk densities of the samples retrieved from the outcrops in Wuppertal, Hagen-Hohenlimburg, and Hönnetal were 2704 ± 39, 2537 ± 128, and 2620 ± 46 kg m⁻³, respectively (Table 4). These mean values were determined independently of whether they were limestone rocks or dolomitic carbonates. In Wuppertal and Hagen-Hohenlimburg, one dolomitic carbonate sample was investigated in each case that differed significantly from the other samples and was not included in the calculation of the mean bulk densities. Taking into account the measurement uncertainty and standard deviation, all samples examined showed comparable densities. A similar trend was observed for the mean grain densities, as these values are also very close together. The mean grain densities of the samples from Wuppertal, Hagen-Hohenlimburg, and Hönnetal amounted to 2801 ± 51, 2777 ± 41, and 2848 ± 32 kg m⁻³, respectively.

Figure 8 Correlation between (a) connected and total porosity and between (b) permeability and connected porosity. Filling and colouring of markers indicate sample and outcrop, 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 (b) indicates identity, i.e. connected porosity corresponds exactly to total porosity.

Laboratory measurements to determine the total and connected porosities indicated that the dolomitic carbonates investigated were significantly more porous than the limestones. This trend was evident at all quarries. On average, the volumetric proportion of pores in the dolomitic carbonates was twice as high as in the limestones among almost all samples. Our uncertainty analysis showed no significant volume of isolated pores, and pre-existing cracks could be determined in the investigated samples. Considering the measurement uncertainties, the total and connected porosities overlapped (Fig. 8a).

Ultrasound P- and S-wave velocities of the dry samples ranged between 3000 and 6500 and 1300 and 3500 m s⁻¹, respectively (Table 4). Regardless of the outcrop, it was shown that the dry P-wave velocities of the limestones tended to be 1000 to 1500 m s⁻¹ higher than those of the dolomitic carbonates. The Wuppertal samples have shown that the dolomitic carbonates and limestones cannot be distinguished according to the S-wave velocities (Table 4). Both rock types had comparable S-wave velocities. The samples of the remaining outcrops showed a clear difference in S-wave velocities between dolomitic carbonates and limestone. In Hagen-Hohenlimburg and Hönnetal the S-wave velocities of the limestones were about 500 to 1000 m s⁻¹ higher than those of the dolomitic carbonates. These differences were also apparent in the drained dynamic Young's moduli (Table 4). The mean drained Young's moduli of the limestone and the dolomitic carbonates amounted to about 60 and 70 GPa, respectively. The mean drained Poisson's ratio, that is, the relative relationship between P- and S-wave velocities, is about 0.35 for both rock types. Upon saturation, P-wave velocities of the samples collected in Wuppertal, Hagen-Hohenlimburg, and Hönnetal increased by up to 4 %, 11 %, and 14 % on average. This increase was also apparent in the undrained Poisson's ratios. The undrained Young's moduli hardly differ from the drained ones (Table 4).

The thermal conductivity measurements showed similar results for limestones as for dolomitic carbonates. However, a slight scattering can be observed. The dolomitic carbonate samples from Wuppertal were an exception, as these samples showed relatively high thermal conductivities that exceeded those of the limestones. On average, the highest and lowest conductivities could be determined on the Wuppertal and the Hagen-Hohenlimburg samples, respectively. Measured thermal conductivities of all samples varied between 1.10 and 3.40 W (m K)⁻¹.

The permeabilities measured correlated as expected with the corresponding porosities, i.e. a permeability increase with an increasing proportion of connected pore volume (Fig. 8b; Table 4). The measured permeabilities of the dolomitic carbonates from Hagen-Hohenlimburg are the highest and averaged to about 4.27 × 10⁻¹⁶ m². The Hagen-Hohenlimburg limestones were, on average, about 3 orders of magnitude less permeable than the limestones. The differences in permeability between dolomitic carbonates and limestone samples amounted to approximately 2 orders of magnitude in Wuppertal. In contrast, dolomitic carbonates and limestones in Hönnetal showed similar permeabilities of the same order of magnitude, that is, permeabilities in the range of 10⁻¹⁸ m². Only one sample from Hönnetal was low permeable (permeability of about 4 × 10⁻²⁰ m²), but here the associated connected porosity was also significantly smaller than that of the comparative samples.

4 Discussion

The discontinuities of deep carbonates in the Rhine-Ruhr area were characterized and described using outcrop scanline surveys. In addition, the results of laboratory investigations of density, porosity, dynamic elastic moduli, thermal conductivity, and permeability of the collected representative rock samples of the Devonian Reef Complex were presented. In the following, we correlate the field and laboratory investigations and address the geothermal potential of these deep carbonates in the Rhine-Ruhr region. We focus our results of the DFN on discontinuity length, density, orientation, and aperture. Characterizing the DFN requires access to the reservoir rock at reservoir-related depth, which is not present. However, we have reasonable cause to expect comparable discontinuity orientations in our geological subsurface model, as in our studied outcrops (Lorenz et al., 1991; Narr, 1991). For example, Narr (1991, 1996) presented how to derive DFN from drill core samples of reservoirs.

4.1 Estimating the geological subsurface model of the carbonate reservoir

The geological setting in the Rhine-Ruhr metropolitan area is highly complex, with details still being debated (e.g. Brix et al., 1988). However, we can name large tectonic events that have influenced local formations to a varying extent:

1. crustal thinning that enhanced the formation of a shelf sea during the Early Devonian (e.g., Dallmeyer et al., 2013).

2. an overall NW movement that reduced clastic sedimentation within the sea and enabled the formation of reef carbonates on the clastic shelf; besides those clastic shelf carbonates, some reef deposits formed on volcanic mounds within the hemipelagic realm (the isolated carbonate complexes on the outer shelf are the geothermal horizon of interest in this study, e.g. Franke et al., 2017; Salamon and Königshof, 2010);

3. the ongoing NW movement that changed the material deposits and thus finally created the Ruhr coal district during the Carboniferous (e.g., Meschede, 2018);

4. at least two different extension sequences during the Early Triassic and Early Cretaceous (Drozdzewski and Wrede, 1994);

5. a NE–SW-directed shortening during the Late Cretaceous (e.g., Brix et al., 1988);

6. a further extensional regime since the Eocene (Kley, 2013).

Due to the geological formation's complexity, our underground model has been simplified concerning the most important tectonic features, such as folds and faults. The foundation of our model is an approx. 150 m thick carbonate layer, dipping northwards at a dip angle of about 30 to 40^∘ (see Fig. 1b; Jansen et al., 1986; Drozdzewski et al., 2007). We expect this layer to be at depths between 4000 to 6000 m (DEKORP Research Group, 1990). This layer forms the Devonian basement, overlaid by interbedded sequences of sand, clay, silt, and coal layers from the Carboniferous period. These sediments can be mapped to the surface (e.g. Brix et al., 1988). Many folds and thrusts are found within these extremely complex interbedded sequences (for more details, see Brix et al., 1988; Drozdzewski and Wrede, 1994).

All studied outcrops are located in the large-scale Remscheid-Altena anticline formation. However, there is a difference between the three outcrop results, which could be influenced by the regional Remscheid-Altena anticline, different stress conditions, and/or a different time of origin (Table 3). The discontinuity orientations were rotated according to the corresponding bedding layers (Figs. 4a, 5a, 6a) to a horizontal position (Figs. 4b, 5b, 6b). The discontinuity unfolding procedure indicates a similar tectonic origin for fracture sets from Wuppertal and Hagen-Hohenlimburg. Discontinuities striking NE–SW are consistent with folding mechanism and are oriented parallel to fold axes that have been studied within the Rhine-Ruhr area (Drozdzewski, 1985; Brix et al., 1988; DEKORP Research Group, 1990; Drozdzewski and Wrede, 1994). The dominant discontinuity strike directions in Wuppertal and Hagen-Hohenlimburg (NNW–SSE) agree with the structure of the regional Remscheid-Altena anticline (Fig. 1b) and the overall assumed mean horizontal stress direction according to the World Stress Map (Heidbach et al., 2016) and additional available stress data (Rummel and Weber, 1993). In western Germany, or to be more precise in North Rhine-Westphalia, the World Stress Map contains a wide variability of mean horizontal stress directions (Heidbach et al., 2016), which can be explained by shallow stress measurements, local anomalies due to weak coal seams, or regional NE–SW-striking thrusts. The observed strong scattering of the vein strike directions in the dolomitic carbonates exposed in Hagen-Hohenlimburg is most probably due to their hydrothermal veins' formation during the Hercynian orogeny (Gillhaus et al., 2003) but can be also explained by the recent study of Michie (2015). Michie (2015) found a crucial correlation between fracture propagation and petrological properties of carbonate rocks, specifically grain size and sorting, porosity and uniaxial compressive strength. Furthermore, Gillhaus et al. (2003) explain that the existing NNW–SSE-striking discontinuities are of post-Hercynian orogeny origin. This leads to the conclusion that the investigated outcrops in Wuppertal and Hagen-Hohenlimburg have been further deformed by later tectonic events. The cause of the slightly different discontinuity strike directions in Hönnetal cannot be specified according to the current state of scientific knowledge and is therefore unknown. However, the discontinuity formation could be explained by various local and temporal stress anomalies and different formation times.

Brudy et al. (1997) showed no significant stress orientation changes have to be expected in the brittle crust with depth. Consequently, it is reasonable to assume that we can expect quite similar discontinuity directions in our target horizon of the deep Devonian limestone, i.e. the potential geothermal reservoir. However, if the carbonate layers exposed in the investigated opencast mines are extrapolated in dip direction of bedding, the carbonate reservoir of interest is approximately located below the cities of Essen, Bochum, and Dortmund at a depth of 4000 to 6000 m (Fig. 2). This extrapolation of the carbonates in dip direction corresponds in this case approximately to the direction of the main horizontal stress (Heidbach et al., 2016; Rummel and Weber, 1993) with respect to the simplified geological setting (Fig. 1a). This assumption allows us to predict local DFN in the deep Devonian limestone (i.e. naturally fractured carbonate reservoir), whose exact depth and characteristics should be verified by additional geophysical prospecting techniques to further describe the geothermal potential of this reservoir (e.g. Hirschberg et al., 2015).

We carried out laboratory experiments under ambient and elevated pressure conditions to gain insights into the potential reservoir's petrophysical properties. The derived porosities of the limestone samples are in agreement with literature values between 1 % and 6 % (Fig. 8a; Gebrande, 1982). The porosity of carbonate rocks may increase via a wide variety of processes. Next to dissolution, cementation, and recrystallization, the process of dolomitization is the most common. Dolomitization describes the geochemical process of replacing Ca ions by Mg ions, forming dolomite from calcite: $\mathrm{2}{\mathrm{CaCO}}_{\mathrm{3}}+{\mathrm{Mg}}^{\mathrm{2}+}\to \mathrm{CaMg}({\mathrm{CO}}_{\mathrm{3}}{)}_{\mathrm{2}}+{\mathrm{Ca}}^{\mathrm{2}+}$. Lucia et al. (2007) pointed out that dolomitization may increase the carbonate's porosity by 13 %. This correlates fairly well with our results (Fig. 8a). Consequently, porosity measurements are in accordance with the field study's outcome: porosity decreases towards the core of the reef formation, in agreement with Homuth et al. (2015b).

The P- and S-wave velocities and the matrix permeabilities and dynamic mechanical properties derived in the laboratory provide statistical information for further numerical simulations of the reservoir. The thermal conductivity results are appropriate for the analysed rock formations (cf. Čermák and Rybach, 1982; Clauser and Huenges, 2013; Jorand et al., 2015). The low connected porosities match with the determined permeabilities. Compacted limestone (Massenkalk), which is characterized by a low connected porosity, showed the lowest permeability, whereas the dolomitic carbonates showed an inverse relation (Fig. 8b). However, the results of the determined matrix permeability seem to be comparable with similar geothermal reservoirs (e.g. Homuth et al., 2015a); we do not expect sufficient reservoir permeability in 4000 and 6000 m depth. Therefore, we have focused on potential DFN within the geothermal reservoir, whose influence on fluid flow in fractured reservoirs has already been discussed in detail (e.g. Guerriero et al., 2013).

4.2 Evaluating a discrete fracture network in the carbonate reservoir: from outcrop and laboratory measurements to implications for fluid flow on a reservoir scale

The uncertainties associated with the description of naturally fractured systems and DFN are strongly dependent on sampling effects such as truncation or censoring (Bonnet et al., 2001; Baecher and Christian, 2005). Truncation effects explain the influence of underestimating small characteristics due to the resolution of the sampling method, i.e. outcrop and survey size. Censoring effects occur when broad characteristics are incompletely observed due to the outcrop size mentioned or when characteristics are selected due to subjective choice (Priest and Hudson, 1981; Lacazette, 1991). While the truncation effects can be reduced by a sophisticated set of data and resolution, the censoring effect can be reduced depending on the method applied, adapted technology, and the amount of data obtained (Santos et al., 2015). In the literature, these uncertainties are typically described as aleatory and epistemic (Der Kiureghian and Ditlevsen, 2009). Various survey methods and their corrections concerning the effects mentioned are described in more detail in Bonnet et al. (2001).

The approach we have chosen to use, the scanline survey method described by Priest and Hudson (1981), is a widely accepted method for the characterization of discontinuity networks (e.g. La Pointe and Hudson, 1985; Lacazette, 1991; Narr, 1996). The advantages of the method are especially the acquisition time and the associated verification as well as the acquisition costs compared to other methods (e.g. lidar outcrop survey; Wilson et al., 2011). The uncertainties due to scanline surveys are generally known and have been critically questioned for decades (e.g. Terzaghi, 1965; Cruden, 1977; La Pointe and Hudson, 1985; Ganoulis, 2008). However, due to the long period of application, the method explained by Priest and Hudson (1981) has been adapted for a variety of geological and outcrop settings, which leads to various error corrections regarding the application case. This allows us to keep the limitations and sources of error as low as possible. We pursued the concerns about epistemic uncertainties and identified the main problems in defining representative areas' locations to evaluate a discrete fracture network in the carbonate reservoir in the Rhine-Ruhr area.

4.2.1 Orientation, roughness, and filling of the discontinuity network

The architecture of each open-pit mine studied is directly related to the directions of the dominant discontinuities, i.e. NNW–SSE and NE–SW (Fig. 3). In addition, the direction NW–SE was documented as another dominant discontinuity direction. These findings fit with the general observations documented in the literature (e.g. Brix et al., 1988; Drozdzewski and Wrede, 1994). Nickelsen and Hough (1967) and Roberts (1974) investigated discontinuity orientations in sedimentary layers with respect to the maximum horizontal stress and their corresponding filling material. The orientations of the respective structures in their area are almost identical to the values we measured. They concluded that the discontinuity directions can be related to residual stresses and stress release during past unloading regimes, which might be transferable to our studied areas (Fig. 7). Besides, our roughness analyses revealed predominating smooth discontinuity surfaces on the field scale (> 10⁻¹ m) and mainly slightly rough discontinuity surfaces on the mesoscopic scale (⩽ 10⁻¹ m) (Table 2). Considering the numerous tectonic events and the repetitive reactivation of pre-existing discontinuities, smooth discontinuities tend to be reasonable in the deep carbonate horizon of interest. Since other studies have already observed that the orientations of regional discontinuities are continuous and consistent over large areas of several thousand square metres (Kelley and Clinton, 1960; Stearns and Friedman, 1972; Hancock and Bevan, 1987), we assume comparable consistency and distribution for the discontinuity sets recorded.

The separate consideration of the three respective discontinuity sets per quarry (Figs. 4b, 5b, 6b; Table 6), as well as the distinction between discontinuities filled with calcite or debris, shows us that the mean orientation of one discontinuity set (set 1; Fig. 7) matches remarkably well with the prevailing orientation of the principal horizontal stress (i.e. NNW–SSE orientation) (Fig. 1; Heidbach et al., 2016; Rummel and Weber, 1993). We could not detect any rotation between the strike direction of calcite- and debris-filled discontinuities (Fig. 7). Hence, the maximum horizontal stress may have changed only slightly, if at all, in the time between the formation of the calcite- and the debris-filled discontinuities. Therefore, we assume that the discontinuity orientation in the reservoir is comparable to the ones observed in the representative outcrops. In addition to the discontinuities' relative orientation to the direction of the main principal stress, the filling is also decisive for whether the discontinuities are potentially open or closed in the subsurface (Laubach et al., 2004). More than half of all our observations showed paleo filling materials such as calcite and, in addition, approximately 17 % showed calcite-enriched slickensides on the discontinuity surface (Fig. 7). Thus, apart from the orientation, no differences in occurrence between calcite- and debris-filled discontinuities can be detected, and we must further consider the orientation of the discontinuities. Even if we cannot predict it with certainty, pre-existing discontinuities oriented parallel to the principal horizontal stress could be observed openly, while discontinuities oriented perpendicular to this stress tend to be closed (Lorenz et al., 1991). Accordingly, the highest discontinuity permeability can be expected in NNW–SSE direction. But there might be open discontinuities that are not necessarily aligned in the principal stress direction and are still open. This is particularly true for discontinuities that are partially filled, for example, with cement (Laubach et al., 2004). In general terms, a large number of veins might be paleo fluid-flow paths striking in NNW–SSE direction, which could be exploited by advanced drilling methods (e.g. hydraulic fracturing; Dahi Taleghani and Olson, 2013) for future geothermal applications.

Before utilization of geothermal energy from the potential deep fractured carbonates in the Rhine-Ruhr Area, the material changes between host rock and vein material might predict discontinuity orientations in depth by further geophysically driven studies (e.g., density changes or reflection coefficients). Further studies have to prove if the material change between the encountered filling materials and the host rock is sufficient for seismic interpretation. Considering the typical seismic resolution concerning the depth of investigation and discontinuity aperture, individual discontinuities will probably not be resolved (Martin et al., 2006). However, the discontinuity density may be sufficient to demonstrate macroscopic anisotropy using seismic methods (Löer et al., 2018).

4.2.2 Connectivity of the discontinuity network

Comparing the discontinuity trace lengths with the true discontinuity spacing of the different discontinuity sets gives a first idea of the connectivity of the discontinuity network. Predominantly discontinuities with short trace lengths occur on the field scale (Figs. 4d, 5d, 6d). This observation was made for all exposed outcrop faces in all three quarries. Almost 70 % of the recorded discontinuities classified as “both ends visible” tended to be between 1 and 2 m long (Table 3). It is conceivable that the observed trace lengths increase with increasing outcrop length and height. The recorded end types suggest that 40 % of all discontinuities investigated exceeded the observed outcrop height (Table 3). In other words, 40 % of all measured trace lengths are probably greatly underestimated by the method. This source of error was attempted to be minimized by the censored semi-trace length analysis (Eq. 2). In contrast, 60 % of all joints counted showed a start and end tip (Table 2). However, a full 3D discontinuity analysis of the entire quarry will likely reflect better trace length results depending on the accuracy. At this point, economic considerations must be made regarding the required survey accuracy and associated costs.

The results of the 2D discontinuity connectivity analysis clearly show that the reservoir in Wuppertal has the highest average 2D discontinuity connectivity, i.e. an average number of discontinuity intersections per discontinuity of 21.42 (Table 6). Both the discontinuity density and the mean discontinuity length of each discontinuity set are significantly higher in Wuppertal than in the other outcrops (Table 6). Accordingly, we can assume that this reservoir probably has considerable geothermal potential due to its discontinuity network. For the Hönnetal reservoir, a discontinuity intersection per discontinuity of 3.71 was determined. One reason for this relatively low value is the reservoir's location on the anticline axis and the associated tectonics (Fig. 1). Although the recorded discontinuities in Hagen-Hohenlimburg exhibit the greatest variety of strike directions (Fig. 5), the estimated 2D discontinuity connectivity is the lowest at this location due to the short discontinuity lengths and the associated low discontinuity densities (Table 6). Based on the observations of many karst formations and altered host rocks by hydrothermal veins in Hagen-Hohenlimburg, we conclude that there might be a higher permeability in the Hagen-Hohenlimburg reservoir and the surrounding areas than in the other studied areas. However, the distribution of this reservoir in the subsurface is still under debate (e.g. Salamon and Königshof, 2010), and its conditions due to karstification still have to be proven.

The densities of the dolomitic carbonates and limestones determined in the laboratory did not show significant differences, which was to be expected (Table 4). However, the comparatively high porosity of the dolomitic carbonates indicated a moderate pore and crack volume. Consequently, a higher density of potentially hydraulically stimulated pores and cracks (or discontinuity network) can be expected in the dolomitic carbonate layers on the field scale. On the other hand, the comparable densities pose a warning regarding the unreflected interpretation of, for example, density measurements in naturally fractured carbonate reservoirs. Such measurements do not appear to be very useful for differentiating between low porous limestone and relatively high porous dolomitic carbonates.

Table 5 Ultrasound velocities and dynamic elastic moduli derived by laboratory measurements for limestone and dolomite samples from Wuppertal (WOH), Hagen-Hohenlimburg (HKW), and Hönnetal (HLO). The abbreviations MKY, MKB, and MKR indicate the coloured rock types Massenkalk yellow (dolomite), Massenkalk black (limestone), and Massenkalk red (dolomite), respectively.

v_P,dry, v_P,wet: P-wave velocity of dry and wet samples; v_S,dry: S-wave velocity of dry samples; ${\stackrel{\mathrm{̃}}{\mathit{\nu }}}_{\mathrm{d}}$: dynamic drained Poisson’s ratio; ${\stackrel{\mathrm{̃}}{\mathit{\nu }}}_{\mathrm{ud}}$: dynamic undrained Poisson’s ratio; ${\stackrel{\mathrm{̃}}{E}}_{\mathrm{d}}$: dynamic drained Young’s modulus; ${\stackrel{\mathrm{̃}}{E}}_{\mathrm{ud}}$: dynamic undrained Young’s modulus. Quoted uncertainties reflect accuracy of the measurements.

Table 6 Overview of the 2D discontinuity connectivity results.

θ: discontinuity strike azimuth; d: discontinuity density; L: mean discontinuity length; P: 2D discontinuity connectivity.

The influence of pore and crack volume in the dolomitic carbonate samples was also reflected in the measured ultrasonic velocities (Table 5). The dolomitic carbonate samples showed the lowest P-wave velocities. The S-wave velocities of dolomitic carbonates and limestones hardly differ, which align with the given densities. Therefore, a warning must be given at this point because it is crucial to observe the temporal and spatial variation of elastic wave velocities during a potential reservoir stimulation of a dolomite and limestone reservoir in order to infer its properties and composition. In the next stage, these findings can be of great importance for interpreting flow rates in the deep carbonate reservoir. Ahrens et al. (2018) demonstrated a uniform correlation between P-wave velocity and the key parameter for transient flow processes experimentally, i.e. hydraulic diffusivity, in the direction of fluid flow during inelastic deformation. Their observation allows the conclusion that monitoring variations in P-wave velocities in field studies can indicate changing fluid flow rates, for example, in developing reservoirs. Consequently, these measurements have the potential to provide insights into the deformation state, connectivity, and hydraulic properties of pore and discontinuity networks in the reservoir.

The low impact of the fluid saturation on the samples' dynamic moduli indicates that heterogeneous pore and crack networks prevail in the investigated samples (Table 5). In the samples' measuring or drilling direction, only a small increase of the P-velocities and thus also of the undrained dynamic moduli occurs. This indicates that the connected pore and crack space tends to be locally distributed (upper Hashin–Shtrikman bound for a simplified two-phase carbonate material; Hashin and Shtrikman, 1963), since ultrasonic velocities only characterize the momentary elastic behaviour along a wave path, which was chosen by the waves due to its above-average rigidity. On the scale of wavelength (1 to 7 mm), the samples are thus probably heterogeneous. Moreover, a velocity anisotropy should also prevail, which is particularly relevant for seismic surveys. However, under the assumption of constant porosity, a possibly existing heterogeneous and local distribution of the pore and crack volume can favour the interconnection of these volumes to large pores or opened cracks more than a homogeneous distribution of many small pores could. In this case heterogeneity possibly favours the fluid flow.

Through hydraulic stimulation, however, the fractures prevailing there and their networks could be reactivated, thus increasing the permeability even though this prediction depends on the discontinuity network's orientation in the stress field and the fluid pressure. However, it is indispensable for geothermal projects to carry out a seismic risk analysis before planning to manipulate the local stress conditions, i.e. the reactivation of fractures.

5 Conclusion

We investigated the geothermal potential of a carbonate reservoir in the Rhine-Ruhr area, Germany, by combining outcrop scanline surveys of pre-existing discontinuities and laboratory measurements of petrophysical properties on a sample scale. The target horizon of interest is an approximately 150 m thick and widely distributed compacted limestone layer (Massenkalk) from the Late Devonian that dips northwards at 30 to 40^∘ and likely occurs at a depth of 4000 to 6000 m. If we extrapolate the course of the exposed Devonian limestone layers, which were investigated at three outcrop analogues, in the direction of the maximum horizontal stress, the carbonate reservoir could hypothetically extend below Essen, Bochum, and Dortmund.

Our petrophysical laboratory measurements on representative outcrop samples showed insufficient matrix porosity, permeability, and thermal conductivity for geothermal energy production. However, outcrop investigation suggests that these units could act as a naturally fractured reservoir at depth. Generally, three discontinuity sets striking in NNW–SSW, NW–SE, and NE–SW directions with dipping angles between 80 and 90^∘ were measured. Considering the numerous tectonic events in the Rhine-Ruhr region and the resulting repeated reactivation of pre-existing discontinuities and thus surface wear and/or degradation, the occurrence of slightly smooth discontinuity surfaces in the reservoir is more likely. We suggest that geothermal exploration efforts focus on discontinuities that are approximately NNW–SSE oriented, bearing in mind that this direction is approximately parallel to the main horizontal stress according to the World Stress Map. Although the discontinuities may not be open, more than half of them showed paleo-filling materials such as calcite beside calcite-enriched slicksides, giving them the potential to be reactivated and improve fluid flow in the carbonate reservoir.

Based on the mean discontinuity length, spacing, and orientation observed at the surface and our discontinuity connectivity analysis, the reservoir at depth most likely has a well-connected network of discontinuities. The hypothetical fractured carbonate formation below Essen could most likely have the comparatively best connected discontinuity network, since a high discontinuity connectivity could be determined for the corresponding reference reservoir in Wuppertal. Despite this, we conclude that the permeability is more significant in the Devonian reef regions facies than in the other areas investigated because of the presence of many karst formations and altered host rocks by hydrothermal veins. This assumption was also confirmed by our laboratory measurements, whose results suggest the presence of heterogeneities that may favour fluid flow.

Our study may provide the basis for improved subsurface characterization concerning the structure and characteristics of the naturally fractured Devonian limestone reservoir in the Rhine-Ruhr area. More sophisticated geophysical prospecting techniques and a combination of pressure- and temperature-dependent laboratory measurements have to be elaborated to verify the reservoir depth and orientation, as well as its discontinuity density and connectivity. It must be demonstrated whether the prevailing discontinuity network can be reactivated, thereby improving the fluid flow in the carbonate reservoir. In addition, the local stress field should be further verified in follow-up studies. The understanding of the discontinuity network's structural characteristics and physical properties and their impact on reservoir properties that, in turn, define the geothermal potential will benefit from, for instance, seismic surveys and in situ measurements in a pilot borehole.