Petrophysical characterization, BIB-SEM imaging, and permeability models of tight carbonates from the Upper Jurassic (Malm ß), SE Germany

Tight carbonate rocks are important hydrocarbon and potential geothermal reservoirs, for example, in CO 2 -Enhanced Geothermal Systems. We report a study of outcrop samples of tectonically undeformed tight carbonates from the upper Jurassic “Malm ß” formation in Southern Germany near the town of Simmelsdorf (38 km NE of Nurem-berg) to understand bulk petrophysical properties in relation to microstructure and to compare models for permeability prediction in these samples. We applied Archimedes isopropanol immersion, Helium pycnometry, mercury injection, gamma density core logging, and gas permeability measurements, combined with microstructural investigations


Introduction
The energy transition requires exploration for unconventional renewable energy sources, such as Enhanced Geothermal Systems (EGS), where hydraulically active fracture corridors are created by hydraulic and chemical stimulation techniques (Stober and Bucher 2013). Here, tight but thermally highly conductive carbonate rocks bear the potential of playing a vital role in the energy transition (e.g., Gosnold et al. 2010;Hofmann et al. 2014). Such reservoirs are characterized by very low matrix permeabilities, typically in the range of 0.001-10 mD (Akanji et al. 2013), implying that fluid transport in low-permeable or tight carbonate rocks is focused to fractures and faults (Al-Obaid et al. 2005;Dimmen et al. 2017;Litsey et al. 1986;O'Neill 1988;Zeybeck and Kuchuk 2002) with only minor fluid flow through the rock matrix (Bohnsack et al. 2020(Bohnsack et al. , 2021. However, microporosity and structures may impact fracture pattern, fracture density, fracture propagation, well logging (e.g., PHI, acoustic log, neutron log), and fluid losses, and thus the effectivity of carbonate matrix stimulation treatments (Barri et al. 2021;Ziauddin and Bize 2007).
A major aquifer in southern Germany is the Upper Jurassic (Malm) limestone reservoir from which geothermal energy for electricity and/or heat supplies is currently (August 2022) produced by > 22 geothermal power plants in the South German Molasse Basin (SGMB) Weber et al. 2019). As geophysical borehole measurements or cores from this low-permeable thermal aquifer are scarce, incomplete, or even non-existent (Bohnsack et al. 2020(Bohnsack et al. , 2021, petrophysical investigations of rock samples from stratigraphically equivalent reservoir analogs are common practice. The tight carbonates of the Upper Jurassic Malm β carbonates exposed in the Frankenalb of Northern Bavaria (Fig. 1) are part of the same stratigraphic unit as the producing geothermal aquifer in the SGMB. In terms of lithofacies, the Malm β carbonates in South Germany predominantly classify as mud-to wackestones (Koch et al. 2005). These units of which maximum burial depth of c. 1100 m has been recently determined represent an excellent opportunity to acquire and compare petrophysical macro-and microstructural parameters of reservoir analogs to geothermal target lithologies in the Munich area of the SGMB 100-200 km to the south, where a similar burial depth was reached (Freitag et al. 2022). From a number of available natural and artificial outcrops we chose a study area NE of Nuremberg (Fig. 1).
The primary goal of this study is to characterize the tight carbonate rocks in terms of their structural anisotropy, their petrophysical and microstructural properties, and the lateral and vertical variability of these parameters. Secondly, by understanding heterogeneities we aim to improve estimates of the contribution of fluid transport between matrix and fracture system at the scale of our field lab (c. 32 × 34 m) and for modeling purposes. As pore systems in tight carbonates are primarily controlled by lithofacies and later modified by diagenetic and tectonic processes (Haines et al. 2016), the petrophysical properties and pore characteristics of the matrix of faulted carbonate rocks are assessed. We also investigate the hydraulic influence of stylolites, which tend to form either barriers (Rashid et al. 2017) or conduits for fluid flow (Bruna et al. 2019). Consequently, stylolites are thought being responsible for most of the variability in petrophysical properties in these rock types (Bruna et al. 2019).
Previous studies focused on a lithofacies-dependent thermophysical and simple petrophysical characterization (porosity and permeability) of core and cutting samples of the Upper Jurassic reservoir (Malm β) in the SGMB and their outcrop analogs in the Northern and Southern Franconian Alb (Homuth et al. , 2015. A similar study, though including geomechanical measurements as well as a petrographic study, was conducted by Mraz et al. (2018) on outcrop reservoir analog samples from the Southern Franconian Alb. Bohnsack et al. (2020Bohnsack et al. ( , 2021 conducted extensive studies investigating potential lithofacies-dependent porosity-permeability relationships of Upper Jurassic limestones (i.e., Malm α and β) and their stress sensitivity (Malm ζ) within the SGMB. Potten et al. (2019) and Potten (2020), on the other hand, characterized Upper Jurassic Fig. 1 Geological map of northern Bavaria, showing major stratigraphic units, structural features, and the test field location ~ 38 km NE of Nuremberg. Background data source: Bayerisches Landesamt für Umwelt, www. lfu. bayern. de. The city data (transparent crosses) were provided by https:// mapcr uzin. com/ (downloaded: 19 July 2021 at 9:45). The administrative region data belong to © GeoBasis- DE/BKG (2021) limestones derived from both the subsurface (243-5225 m TVD) of the SGMB and from outcrops ~ 40 km N of our test area in terms of their geomechanical properties. Taking this a step further our paper additionally evaluates the strength of Franconian Alb reservoir analog samples and we discuss to which degree their petrophysical properties are applicable to buried reservoir rocks of the same stratigraphic units in the SGMB. By applying a combination of state-of-the-art direct, indirect, and graphical methods, we try to achieve more comprehensive understanding of factors that are primarily controlling fluid flow within the rock matrix so that more reliable permeability predictions can be made.
A second aim is more technical. As for many geothermal reservoirs only a limited number of samples is available, petrophysical data of the target rock formation determined from cuttings are essential for subsequent reservoir quality evaluation and modeling. Plug ends, that result from preparing the plug cores in the quarry compared to cuttings in terms of their relatively small size, are investigated by combined mercury intrusion capillary pressure (MICP) measurements and broad ion beam scanning electron microscopy (BIB-SEM). Based on MICP measurements, bulk information on quantitative pore characteristics are obtained, including the pore throat size distribution, down to a pore throat size of only 3 nm (Clarkson et al. 2013;Giesche 2006;Okolo et al. 2015;Webb 2001;Xu et al. 2018a, b;Zhao et al. 2018). A variety of models enable permeability estimation based on MICP data. Therefore, this method is considered as a standard for the investigation of tight reservoir rocks (Gao et al. 2016;Okolo et al. 2015). We evaluate various permeability models, mainly based on MICP measurements, by comparing them to steady-state Darcy flow/permeability measurements. Although for tight rocks typically transient methods (e.g., pulse decay) are used for measuring the permeability, their benefit in terms of receiving more reliable results in shorter time is still a matter of debate (Sander et al. 2017), particularly for rocks with permeabilities > 10E−20 m 2 . The BIB-SEM application enables microstructural investigation of relatively small rock pieces and quantification of pore sizes down to c. 5 nm (Klaver et al. 2012). Additional capillary tube modeling then provides a solid basis for comparison with permeability measurements in tight rocks (Philipp et al. 2017;Sinn et al. 2017).

Geological background
During the Jurassic, the study area was occupied by an epicontinental sea as part of the northern to northwestern Tethys shelf (e.g., Koch et al. 2005;Meyer 1996;Pieńkowski et al. 2008). Stress-induced lithospheric deflections related to far-field compression (Scheck-Wenderoth et al. 2008) and a wrench-dominated tectonic regime at the southern end of the North Sea rift system (Pharaoh et al. 2010) led to rapid shallowing of the South German shelf areas during the latest Jurassic to earliest Cretaceous. This resulted in the deposition of peritidal carbonates and anhydrites, and ultimately led to their exposure and pronounced erosion (Bachmann et al. 1987;Schröder 1968;Vejbaek et al. 2010;Voigt et al. 2007Voigt et al. , 2008Ziegler 1990). Widespread post-Jurassic, Mid-to Late-Cretaceous erosion and karstification due to compressional stresses associated with the reactivation of NW-SE and NNE-SSW striking normal faults (e.g., Peterek et al. 1996;Schröder 1968Schröder , 1987Voigt et al. 2007;Wagner et al. 1997) were followed by a sea-level rise, recorded by a northward marine transgression (Bachmann et al. 1987;Scheck-Wenderoth et al. 2008). After a phase of tectonic quiescence, N-S to NW-SE oriented far-field compressive stresses within the Central European lithosphere related to the Alpine Orogeny (e.g., Scheck-Wenderoth et al. 2008;Ziegler 1987;Ziegler et al. 1995) and/or the onset of Africa-Iberia-Europe convergence (Kley and Voigt 2008) caused the Late-Cretaceous inversion and the formation of N-S to NNW-SSE trending strike-slip faults accompanied by stylolites with orientations vertical to bedding, bedding parallel, and sub-parallel to the strike of normal faults (Koehler et al. 2022) as well as the reversed reactivation of NW-SE striking normal faults (Zulauf 1993). Ultimately, Late-Cretaceous inversion not only led to a cessation of Cretaceous sedimentation but also to the erosion of several hundreds of meters of Mesozoic sediments and widespread reverse reactivation of normal faults (Bachmann et al. 1987;von Eynatten et al. 2021;Fazlikhani et al. 2022;Scheck-Wenderoth et al. 2008;Voigt et al. 2008Voigt et al. , 2021. A second major uplift phase between latest Late Cretaceous and Palaeocene was caused by the combined effects of the Alpine continental collision Reicherter et al. 2008;Schröder 1987;Wagner et al. 1997;Ziegler 1987) and mantle-induced domal uplift in the area of the Upper Rhine Graben Rift (URG) to the west of the study area (von Eynatten et al. 2021) (Fig. 1 inset). Consequently only erosional remnants of Cretaceous and Cenozoic sediments of < 100 m thickness are preserved in the Frankenalb of N Bavaria (Meyer 1996) but significantly higher thicknesses of > 5000 m are preserved within the SGMB to the south of the study area (Bachmann and Müller 1992;Meyer 1996).
A detailed description of limestone facies and petrography in the study area is given by Koch et al. (2005). Lower Kimmeridgian mud-to wackestones are overlying well-bedded Oxfordian limestones, separated by thick marl beds of the Platynota zone (Koch et al. 2005;Meyer 1974;Zeiss 1977). Thick-bedded Middle Kimmeridgian limestones form the top of the section.

Sampling
To minimize surface weathering effects, our test field is located in an active quarry near the town of Simmelsdorf, about 38 km to the northeast of Nuremberg (Fig. 1), exposing Upper Jurassic (Malm) limestones. Exposed sub-horizontal beds vary in thickness from 0.1 to 0.6 m and were consecutively labeled B-2 to B10 ( Fig. 2B and C). A total of 40 carbonate Malm β samples were collected for this study's workflow ( Fig. 2A) from a NNE-SSW outcrop section ( Fig. 2B and C), which is oriented sub-perpendicular to two NW-SE trending normal faults forming a graben structure. The offset of the normal faults is relatively small, approximately 0.5 m. Blocks (~ 50 cm × 50 cm × 40 cm) of each bed ( Fig. 2B and C) were systematically sampled in a vertical section at approximately six meters distance to the fault and numbered 0 to 10. For a horizontal section, blocks of similar size were extracted from bed B1 every few meters and labeled 1E-H (Fig. 2B) and 1A-D (Fig. 2C). The blocks were carefully removed by an excavator, layer by layer. They were cut on-site into smaller pieces and a selection of samples from the vertical section was drilled with diameters of c. 6 cm and 3 cm to get cores (6 × 20 cm) and plugs (3 × 4 cm) suitable for further lab-based Multi-Sensor Core Logger (MSCL), gas permeability (plugs), helium Fig. 2 A Workflow from sampling to analysis and data usage; B and C are photos of the exposed quarry walls (see Fig. 1 for quarry location) illustrating consecutive labeling of beds B-2 to B10 (white background) and positions of the sampled rock blocks (yellow background), labeled 0 to 10 (vertical sampling) and 1A to 1H. (horizontal sampling of B1). Traces of NW-SE striking normal faults and inferred polarity of offset in red

Microstructural investigations via BIB-SEM
The plug ends of each sample from the vertical section were sub-sampled for microstructural analyses using BIB-SEM. The 5 sub-samples have maximum dimensions of 4 × 10 × 6 mm 3 (height × width × depth) and were taken from a macroscale relatively homogenous, representative part of the plug ends. The BIB cross-section was polished using the JEOL SM-09010 cross-section polisher producing a 1 mm 2 planar Gaussian-shaped cross-section on the sample by removing approximately 100 μm of material in 8 h (Klaver et al. 2012). Prior to SEM image mapping by a Zeiss Supra 55 Field Emission SEM, samples were prepared with approximately 7.5 nm of tungsten coating. In the SEM, the pore space was imaged with the secondary electron (SE2) detector systematically in a raster pattern throughout the whole BIB section at about 50 locations without overlap. At each location 5 images were acquired with different magnifications of 2500, 5000, 10,000, 20,000, and 40,000 times. This approach enables unbiased sampling of the whole section and scans a wide range of pore sizes simultaneously. The visible BIB-SEM porosity was quantified by calculating the fraction of visible pore space using image segmentation (Klaver et al. 2012) for 5 selected crosssections on several single images. The other BIB cross-sections were used for qualitative investigations as the pore microstructures were similar to each other according to visual inspection.
Two additional sub-samples were injected with Wood's Metal (LMI-BIB-SEM) following the workflow of Klaver et al. (2015) at 100 and 200 MPa, respectively. As Wood's Metal has non-wetting properties comparable to mercury, it is a similar principle as MICP, though the Wood's Metal is solid at room temperature, which enables visualizing the metal-filled pore space. Subsequently, after BIB low angle polishing in a Leica TIC3X BIB, the metal-filled pore space was imaged in the SEM to qualitatively evaluate the pore connectivity in the carbonate matrix.

Porosity
Multi-Sensor Core Logger (MSCL) The Geotek Multi-Sensor Core Logger (MSCL) system includes an assembly of tools that log the geophysical and geochemical properties of cores. For this study, only data acquired by the gamma density tool were used, and a total of 0.94 m core lengths were measured at 1 cm intervals. The samples were exposed to a focused beam of gamma rays (energies principally at 0.662 MeV) that become attenuated by Compton scattering while they pass through the core with the degree of attenuation directly relating to the diameter and the electron density of the core. The bulk density (ρ bulk ) of the measured core section was then calculated by measuring the number of transmitted gamma photons that passed through the core unattenuated (I), considering the core thickness (d): where µ is the Compton attenuation coefficient and I 0 the gamma source intensity (Weber 1997). The porosity Φ GD (%), in the following termed gamma ray (GD) porosity, is calculated by applying gamma ray-derived bulk density (ρ bulk , g/cm 3 ), matrix density (ρ matrix = 2.71 g/cm 3 for CaCO 3 ), and the density of air (ρ Air = 0.001225 g/cm 3 ): Archimedes (buoyancy) isopropanol immersion method The Archimedes method was applied to a total of 40 drilled cylindrical core plugs. Advantages, disadvantages, operating principle, and potential measurement errors of the method are thoroughly discussed by Hall and Hamilton (2016). The Archimedes porosity Φ (%) was calculated via Eq. 3: where V P (m 3 ) is the open pore space volume filled with isopropanol, derived from subtracting the dry mass m dry (g) from the saturated mass m sat (g). Saturating the open pore space with isopropanol was achieved by fully submerging the sample in an isopropanol bath and placing it in a vacuum chamber until no visible air bubbles were exiting the sample. V tot (m 3 ) is the total volume of the sample, including the open pore space and the solid rock matrix volume, which is equal to the difference between the saturated mass m sat and the weight m im (g) of the sample submerged in isopropanol. The weights were determined by a Sartorius ED2245 working at an accuracy of 0.1 mg. Only open pores that are connected to the open pore network of the sample are determined by this method, as the isopropanol cannot access closed pores (Zinszner and Pellerin 2007).
He pycnometry Helium porosity measurements were carried out on 40 selected plugs from the horizontal and vertical sections. The porosities were calculated based on the difference between the total dry plug volume (calculated from diameter and length measures by a high precision gauge) and the matrix volumes determined by a Micromeritics AccuPyc 1330 pycnometer. The instrument measures the skeletal volume of a sample at an accuracy of 0.03% (i.e., matrix volume) by the gas displacement technique, based on the ideal gas law. The use of helium gas enables the filling of connected pores as small as 0.1 nm in diameter.

Mercury injection capillary pressure (MICP)
Mercury intrusion porosimetry is one of the most widely used methods to determine the bulk porosity and the pore (throat) size distribution by utilizing the property of non-wetting liquids that only intrude capillaries under pressure. Washburn (1921) described this relationship between pressure and capillary diameter: where P is the pressure (Pa), γ the surface tension of the liquid (in this case mercury) (mN/m), θ the contact angle of the liquid (for mercury θ = 140°), and d cap the diameter of the capillary (m). The intruded volume of mercury entering the pores at each pressure increment is recorded and from that, the pore (throat) size distribution is derived, whereas the total porosity can be calculated from the total intruded mercury volume (Abell et al. 1999). The r 35 values correspond to the pore diameter (µm) at 35% mercury saturation, while r Main (µm) gives the pore diameter at which the largest amount of mercury intruded, with both values listed in Appendix 1: Table 4. We used the PoreMaster 60 by Quantachrome, operating at an accuracy of ± 1% fso (full Scale Output) of sample cell stem volume on 24 samples.

Permeability measurements
Permeability k Ar measurements (in m 2 ) were carried out on 12 samples, using argon gas. The 1-inch diameter plugs were placed inside a stainless-steel cylinder (autoclave) and permeability determinations were carried out at increasing gas pressure steps of 1 bar with a confining pressure P conf of 15 bar. For each plug, the flow rate through the sample was measured at six different injection pressure steps P Inj , from 3 to 8 bar to correct for the so-called "Klinkenberg Effect" (Klinkenberg 1941). After each pressure increment, the flow was measured once measurements were stable for at least ten minutes before proceeding to the next step. A similar method was applied to determine the permeability (k Air ) of all samples (34), using compressed air instead of argon. Once a stable flow after a minimum of 10 min was established at a particular pressure increment, the gas (air) flow rate exiting the sample was measured every second over a period of 30 s. An "Aarberg" mass flow meter operating at an accuracy of 1% at flow rates between 0 to 50 ml/min was used as logging device. Due to mechanical limitations, only a confining pressure of 8 bar could be applied. Therefore, the maximum flow rate through the sample was reached at a pressure of 5.8 bar, before the through-flowing air might have bypassed the sample due to a too low difference between confining and injection pressure. A Klinkenberg correction for the permeability measurements with air was not possible, as we received negative slippage factors. Possible reasons for that are discussed later. We, therefore, also used the uncorrected, mean measured permeabilities for the comparison to other applied methods and models.

Permeability models
Many different permeability models that were calibrated or validated on different sample sets over the last decades were applied to our data set but only the permeability estimation models that performed best will be treated and discussed in this study.

Models based on percolation theory
Permeability estimation based on MICP data was introduced by Thompson (1986, 1987). Their model is based on the percolation theory (K-T Model) and relates the pore diameter l (m) to the intrinsic permeability k KT (m 2 ). When l is optimally selected, k KT (m. 2 ) can be derived via Eq. (5): where Φ (%) represents the porosity of the rock, l c (m) the critical length, l h max (m) the maximum hydraulic length, and f(l h max ) the fraction of the whole rock that is filled by mercury at l h max . While critical length l c is defined as the critical pore diameter at which mercury can finally percolate through the sample (equal to the steepest slope of the capillary pressure vs. cumulative porosity curve after cut-off ), l h max corresponds to the capillary pressure, where the product of the mercury saturation and the cubic pore throat diameter, f(l h max ) * l 3 , are at its maximum (Nishiyama and Yokoyama 2014;Rashid et al. 2015). Originally, the theoretical consideration that the pore's diameter is equal to its length led to the constant 1/89 (Nishiyama and Yokoyama 2014;Rashid et al. 2015). However, this constant was empirically determined for porous rocks. As this study focuses entirely on tight carbonates, a value of 1013/89 for C was used instead, based on the work of Rashid et al. (2015) recommending this value for tight carbonates. RASHID et al. (2015) apply eight permeability models, all based on the Poiseuille Model. We tested a selection of these models, too. The Winland Model was originally introduced in various unpublished reports between 1972 and 1976, which we could not obtain; we therefore reference published studies by Comisky et al. (2007), Gunter et al. (2014), and Rashid et al. (2015). The Winland Model uses the radius r 35 (μm), which is calculated using the Washburn equation (Eq. 4) at a mercury saturation of 35% (Rashid et al. 2017) and relates it to permeability k W (in m 2 ) according to where C W , a 1, and a 2 are empirically determined variables (−), Φ is the porosity (%), and cf (= 9.86923*E−16) is the factor for converting milliDarcy (mD) to m 2 . These variables were derived from the calibration of Winland's equation on a data set consisting of 82 samples, 56 of which were sandstones and 26 carbonates (Klinkenberg-corrected permeabilities), as well as 240 samples, where only uncorrected air permeabilities were known. The calibration resulted in the following values: C = 49.4, a 1 = 1.70, and a 2 = 1.47. Dastidar et al. (2007) introduced another Poiseuille-based permeability model (Dastidar Model), calibrated on tight gas sandstones. The authors suggest taking the entire pore throat spectrum into account when estimating the permeability from MICP data. They introduced a length scale based on the geometric mean of the pore throat radius (r wgm ) which is calculated:

Poiseuille-based models
with the pore throat radius R i at the ith capillary pressure (m), the total number of incremental pressure steps n, and w i the ratio of the incremental mercury volume intruded into the sample at the specific capillary pressure p i (Pa) and the total mercury volume intruded (m 3 ). With this, we can calculate permeability k D (in m 2 ) after Dastidar et al. (2007), where cf is the factor for converting mD to m 2 : An alternative Poiseuille-based model is the CT Model (or CTM) that can be directly applied to the pore geometries determined from segmented BIB-SEM images. It assumes that the flow through the rock is analog to laminar flow through a bundle of pipes. As pore networks in rocks are never perfectly straight round tubes, but follow a tortuous path, a modified version of the Hagen-Poiseuille equation, taking into account the tortuosity factor τ (−) (Philipp et al. 2017;Sinn et al. 2017), is employed to determine permeability k H-P (in m 2 ) according to the CT Model: where r i (m) equals the hydraulic radius (pore area divided by its perimeter), Φ i is the porosity (%) of each segmented pore in the BIB-SEM image, and cf is the factor for converting mD to m 2 . The tortuosity value can be a fitting parameter or taken from literature data. For this study, a tortuosity of 2.0 was assumed initially, as tight carbonates are slightly less tortuous than tight sandstones (Cai et al. 2019) with typical values of 2.1 (Du 2019).
The Kozeny-Carman Model (K-C Model) is an extension by Carman (1937) that bases on the permeability model developed by Kozeny (1927). He used the specific surface area related to the rock volume S 0 (1/m) and the effective porosity Φ eff (%), hence the pore space contributing to fluid flow (Fens 2000): with the Kozeny constant c and τ as the tortuosity factor, the permeability k K-C (m 2 ), based on the K-C Model), and cf the factor for converting mD to m 2 . S 0 was obtained from MICP data. The tortuosity factor τ is derived from the optimized CT Model and the Kozeny constant c for cylindrical capillaries is 1.57 (Carman 1937).

Empirical models
Although various empirical equations exist, we only applied the models developed by Bohnsack et al. (2020), Saki et al. (2020), Lucia (2001), and Jennings and Lucia (2003) in this study.
Based on the porosity-permeability relationship measured on a subset of ~ 50 mudsupported limestones out of a set of 363 Upper Jurassic limestone core samples, Bohnsack et al. (2020) inferred the following power law (termed Bohnsack Model) where k B is the permeability (m 2 ), Φ the effective water porosity (%), and cf the factor for converting mD to m 2 . Saki et al. (2020), on the other hand, established the following relationship (termed Saki Model) between gas permeability, porosity, and pore/throat diameters from 187 Freitag et al. Geothermal Energy (2022) 10:30 sandstone, limestone, and dolostone samples derived from 8 different Iranian oil and gas fields: where k S is the gas permeability (m 2 ), Φ the porosity (%), r 35 is the smallest pore throat radius (μm) filled by mercury at 35% mercury saturation, and cf is the factor for converting mD to m 2 . An extensive study on a variety of limestones (n = 416) was conducted by Lucia (2001) and Jennings and Lucia (2003). They related rock-fabric petrophysical classes and interparticle porosity to permeability via a multilinear regression, termed global porosity-permeability transform (GPPT). Each rock-fabric petrophysical class represents a different type of pore distribution and interconnection (Lucia 1995). Three classes are thereby distinguished and assigned a specific rock-fabric number (rfn) (0.5-4.0): class 1 represents grainstones and coarse-crystalline dolostones with an rfn of 0.5-1.5, class 2 includes grain-supported packstones and medium-crystalline dolostones with an rfn between 1.5 and 2.5, and class 3 comprises mud-supported limestones and fine-crystalline dolostones with an rfn of 3.5-4.0 (Lucia 2001). The GPPT Model is given by where k GPPT is the rock permeability (in m 2 ) based on Lucia (2001) and Jennings and Lucia (2003), A = 9.7982, B = 12.0838, C = 8.6711, D = 8.2965, Φ ip the fractional interparticle porosity (effective porosity), and cf the factor for converting mD to m 2 .

Geomechanical properties
Geomechanical parameters, such as the dynamic Poisson's number ν dyn (−), dynamic Young's Modulus E dyn (GPa), the dynamic Bulk Modulus K dyn (GPa), and the dynamic Shear Modulus G dyn (GPa), are calculated from P-(V P ) and S-wave (V S ) velocities (m/s) which were measured on 33 limestone samples using an automated core logger equipped with a Geotron.UKS12 ultrasonic device. This comprises conical-shaped piézoelectric stainless-steel p-wave transmitter and receiver probes, operating with pulse transmission at a frequency of 80 kHz and measuring the time taken by the first p-wave to cross the sample. Only plugs that were drilled vertically to bedding were measured. For each sample, five measurements were conducted for reproducibility reasons. For a detailed description of the operating principle of the ultrasonic core logging device we refer to the study of Filomena and Stollhofen (2011). The dynamic geomechanical parameters were derived from the following equations (Schön 2015): A detailed description of the derivation of the relative porosity change ΔΦ (%) with increasing depth, hence effective pressure P e (Pa), is given by Bohnsack et al. (2021) and can be calculated after Cheng (2016) as follows: with Φ i representing the initial porosity (%), measured with the Archimedes method. However, from Cheng (2016) it does not emerge whether he refers to the dynamic (G dyn ) or the static (G stat ) shear modulus. Therefore, we calculated the relative porosity change for both shear moduli. As the static shear modulus cannot easily be calculated from the dynamic shear modulus, we used the empirical correlation by Bastos et al. (1998) for limestone samples: Further models were introduced by Bohnsack et al. (2021) and Hu et al. (2020), applying solely porosity stress sensitivity coefficients that were adjusted to the specific data set and the differential pressure P e -P i (atmospheric pressure). As the stress sensitivity values were not available for our data set, we did not include these relationships in our study.
The relative change in permeability Δk with increasing effective stress can be calculated from various equations, integrating different empirically determined parameters. Applying β as dimensionless constant for the stress sensitivity exponent, provides the following exponential law (Bohnsack et al. 2021;David et al. 1994;Xu et al. 2018a, b): For β, β Min (= 28.3) and β Max (= 46.3) values were determined by Bohnsack et al. (2021) for limestones of the Upper Jurassic (Malm Zeta) from two drill cores in the SGMB that were subsequently used to calculate Δk Min and Δk Max , respectively. Again, two further models by Shi and Wang (1986) and Katsube et al. (1991) describe the relative permeability change with increasing effective stress, solely applying stress sensitivity-describing parameters adjusted to the specific data set and the differential pressure P e -P i . As the stress sensitivity parameters were also not available in this study, we excluded these models from our study.

Results
The investigated Oxfordian to Lower Kimmeridgian limestones classify as mud-to wackestones (Cohen et al. 2013;Koch et al. 2005) and can be described as mud-supported limestones, organized in horizontally layered, continuous beds (Fig. 2) which only rarely contain fossil remnants or vugs filled with pyrite. These beds are separated by thin marlstone layers. The alternating lime-marlstone layers are crosscut by three normal faults, each of which showing an offset of ~ 2 m. In general, these limestones are finegrained and have a very homogeneous appearance. Only rarely they are bearing fossil 4G dyn/stat . remnants and occasionally horizontal, oblique, and vertical pressure solution seams or stylolites are developed with increasing frequency toward the faults.

Pore throat sizes and pore connectivity
A summary of all petrophysical measurement results is listed in Appendix 1: Table 3,  Table 4, and Appendix 1: Table 5. In total, 24 samples were analyzed using MICP. The results for the vertical section (Fig. 3A) show a unimodal pore diameter distribution for most samples. The samples 0, 1, 7, 8, and 9 have pore throat diameters of around 11-15 nm, and samples 2, 3, 4, 5, 6, and 10 had larger pore throats, with diameters of 31-39 μm. The pore throat size distribution of most samples from the horizontal section are clearly conform to each other (Fig. 3b). All samples had pore throat diameters of 9-15 nm, except sample 1F which had slightly larger pore throat diameters averaging 23 nm. However, some samples (1G, 4, and 10) indicate continued mercury intrusion even at the device's maximum pressure, corresponding to the smallest pore diameter of 3 nm (Fig. 3). All samples were conformance corrected (> 4 µm) for surface irregularities (Newsham et al. 2004;Sigal 2009). In conclusion, the samples are relatively homogeneous and reflect a narrow pore throat size distribution with small variations between 8 and 40 nm and no trends related to fault proximity or bed thicknesses (see Fig. 1b for corresponding sample positions relative to the fault).  The MICP breakthrough fits well with BIB-SEM observations on the Wood's Metal injected (LMI-BIB-SEM) samples at 100 and 200 MPa, equivalent to pore throat diameters of approximately 15 and 9 nm (Fig. 4). The sample injected at lower pressure (100 MPa) showed both filled and unfilled interparticle pores (Fig. 4A), indicating a percolating network that was not fully reached by the injection at this low pressure. Instead, the samples injected at twice the pressure (200 MPa) showed that virtually all visible interparticle pores were filled with metal (Fig. 4B).

Microstructure and pore geometry
All plugs prepared for the gas (Ar) permeability tests were macroscopically investigated regarding the presence and orientation of stylolites (Appendix 1: Table 3*). A quarter of the plugs contained stylolites identifiable at macroscale, both sub-parallel to the bedding (Fig. 5A) as well as sub-perpendicular to bedding (Fig. 5B). Such pervasive stylolites were clearly visible in the SEM and were examined at very high detail in the BIB cross-section to determine and characterize different types of microporosity (pores < 1 micron) and nanoporosity (pores < 100 nm) (Fig. 5C-F). Figure 5D illustrates the presence of nanoporosity along the stylolite plane in contrast to the almost tight nature of the limestone's matrix containing only few clusters of interparticle nanopores. Insoluble minerals, like clay, quartz, dolomite, and denser minerals, accumulated as solution residues at the stylolites' peak tops contributing to interparticle/intercrystalline nano-/microporosity.  (Fig. 6B), but at higher resolution a matrix with clusters of moldic pores (Fig. 6C), and nanoporosity associated with larger fossil remnants (Fig. 6D). Besides the moldic pores, the matrix also contains interparticle pores, both being usually smaller than 1 micron ( Fig. 6C and D). However, these pores that are associated with partly dissolved fossil remnants and which are recognizable by their specific pore arrangement resembling the shape of fossils are not common in the investigated BIB cross-sections (Fig. 6D). The BIB cross-sections of C5V and C6-2 V (Appendix 1: Table 3) also contain some much larger interparticle (macro)pores (e.g., Figure 6E) which significantly contribute to the visible porosity. In general, very few but spatially distant macropores were observed. Figure 6F shows the typical microstructure at high resolution showing typical (tri-) angular interparticle and intercrystalline pores relatively close to each other. The smallest visible pores that could be segmented are a few pixels in size.
Overall, the average matrix is characterized by pores mostly in the 100 to 1000 nm pore size range in BIB-SEM (Fig. 7A) indicating a general increase of visible nano-and microporosity with improved magnification. However, this increase stagnates at higher magnification, suggesting that not much additional nanoporosity becomes visible beyond 40.000 × magnification (pixel size 7.3 nm), or below the practical pore resolution (PPR). This is also validated by uniform pore diameter distributions (Fig. 7B), independent from the applied magnification. For the permeability estimation via the CT Model, the measurements at 14.7 magnification (equal to 20,000×) were applied to improve pore detection.

Porosity
The MSCL-logged cores show gamma density porosities varying between 3.7 and 5.9% (mean 5.0 ± 0.7%) that are, on average, slightly higher than porosities derived from He pycnometry, ranging from 2.5 to 6.0% (mean 4.3 ± 1.0%) ( Fig. 8A and B). The average Archimedes porosities of the different limestone beds in the vertical section range from  Table 3), with an overall average of 3.0 ± 0.7%. Similar ranges and mean porosities were determined by both, the BIB-SEM and MICP methods.
The reported porosities of all applied methods are derived from the vertical section. Complementary Φ Ar , Φ He , and Φ Hg porosities from the same beds and within the horizontal section showed comparable porosity variability and no trend toward increased or decreased porosities with distance from the fault (Fig. 8C). Hence, we conclude that the bulk matrix porosity is almost homogeneous with an average MICP porosity of 2.9% (1.8 to 4.4%). The visible BIB-SEM porosities Φ BIB-SEM are also on average 2.9%, varying from 2.1% to 3.5%. Macropores are present; however, they are rather distant from each other in the tight matrix so do not affect the fluid transport properties significantly. Additionally, (micro-) fractures were purposely excluded from the BIB-SEM analyses as our focus was on the microporosity of the matrix.

Permeability and permeability models
The measured air permeabilities K Air vary between 1.9E−18 m 2 and 10.8E−18 m 2 with a mean of 6.8E−18 ± 2.1E−18 m 2 . The lowest determined permeability equals the detection limit of the device, at least for air injection pressures up to ~ 6.8 bar. A mean of 1.4 3) E−17 m 2 with a more confined range between 1.2E−17 m 2 and 2.1E−17 m 2 was determined for the same samples applying argon gas (Fig. 9A). On average the determined air permeabilities are a factor approximately two times lower than the argon permeabilities with no obvious correlation (Fig. 9A). No correlation between the measured air permeabilities K Air and the maximum hydraulic length l h max nor the r 35 (Appendix 1:   (Fig. 9B), which also applies to the correlation between the measured argon permeabilities K Ar and the maximum hydraulic length l h max . In contrast, a positive correlation between measured argon permeabilities K Ar and the r 35 is indicated in Fig. 9B. Correlating air permeabilities with porosities derived from various methods shows no distinct relationship between these properties (Fig. 9C). This is also true when comparing the porosity results of all four applied methods with the permeabilities measured with argon gas (Fig. 9C).
The highest K Air values were found in layers B0, B4, and B7, though, similar to the porosity measurements, no vertical (Fig. 9D) or horizontal (Fig. 9E) trends toward the normal fault can be identified from the matrix permeabilities of the sampled beds. Plug permeabilities from the horizontal section (Bed B1, Fig. 9E) consistently show very low values in the range between 7.6E−18 m 2 and 14.3E-18 m 2 (Appendix 1: Table 4). Furthermore, no permeability trends in N-S, E-W directions, and vertical sampling successions were determined (Appendix 1: Table 4), suggesting that the rocks are behaving relatively isotropically regarding fluid flow. All permeability models applied in this study (Appendix 1: Fig. 11) show at least moderate correlations with measured permeability values. With exception of the CT Model (Appendix 1: Fig. 11E), which is solely based on the BIB-SEM analysis, we only display permeability prediction models based on Φ Hg , Φ He , and Φ Ar for better clarification, as Φ Gd and Φ He , as well as the Φ Ar , Φ Hg , and Φ BIB-SEM values, show very low variation among each other (Fig. 8A).
Comparing the model-based permeability estimations to air permeability measurements (K Air ), decent fits, hence low mean residual square error (MRSE) ( Table 1), were achieved for the Katz-Thompson (K-T) Model (Appendix 1: Fig. 11A) and the Winland Model using Φ Ar and Φ Hg (Appendix 1: Fig. 11B). The same applies to the Kozeny-Carman Model (K-C) using the best-fit tortuosity factor of τ = 1.57 that was achieved by optimizing the CT Model applying He porosity (Appendix 1: Fig. 11D), and the Bohnsack Model using Φ Ar and Φ Hg (Appendix 1: Fig. 11F). Great fits are achieved for the CT Model (Appendix 1: Fig. 11E) and the Saki Model (Appendix 1: Fig. 11H), both being well within the 2.5 variance factor from the perfect fit. The Winland Model using Φ He (Appendix 1: Fig. 11B), the Dastidar Model (Appendix 1: Fig. 11C), the Bohnsack Model also using Φ He (Appendix 1: Fig. 11F), and the GPPT Model (Appendix 1: Fig. 11G) strongly overestimate the measured permeabilities. In contrast, the Kozeny-Carman Model (K-C) using Φ Ar and Φ Hg (Appendix 1: Fig. 11D) underestimates the air permeability K Air .
In general, similar observations are made when comparing the permeability models to the measured argon permeabilities (K Ar ), although with a better correlation. The latter applies except the K-T (Appendix 1: Fig. 11A) and the K-C Models (Appendix 1: Fig. 11D). Particularly well fitting permeability estimations are achieved by applying the Capillary tube (Appendix 1: Fig. 11E) and the Saki Model (Appendix 1: Fig. 11H).

Mechanical properties
We calculated the dynamic elastic parameters (Poisson's number, the dynamic Young's Modulus, and the Bulk Modulus) from sonic velocity measurements using Eqs. (5)-(9), respectively ( Table 2). For samples where the shear velocities were not available and hence, the Poisson's number could not be calculated, and the average value of ν = 0.27 Table 1 For determining the best permeability estimation model, the mean residual square error (MRSE) in relation to the corresponding air K     derived from all other samples was used. As we achieved very similar velocity values despite the samples were analyzed in two different laboratories, the measurements are considered as reliable. All limestone beds are rated as very homogeneous in terms of their geomechanical properties, with the dynamic E-, K, G dyn , and G stat moduli ranging between 65 and 85 GPa with a mean of 73 ± 5 GPa, 41-66 GPa with a mean of 53 ± 5 GPa, 24-33 GPa with a mean of 27 ± 2 GPa, and 14-20 GPa with a mean of 16 ± 1 GPa, respectively.
The relative porosity and permeability changes at a depth of ~ 2000 m TVD (≈ 25 MPa), corresponding roughly to the burial depth of the Malm reservoir in the SMGB according to Drews et al. (2020), are listed in Appendix 1: Table 5. Based on Eq. (9), a relative porosity change of 0.06% to 0.08% with a mean of 0.07 ± 0.01% based on the G dyn is to be expected, while calculations using the static shear module G stat suggest a relative porosity change of 0.10% to 0.14%. with a mean of 0.12 ± 0.01%. A larger relative change at same burial depth is predicted for the permeability by applying Eq. (10), using two different values β min (= 28.3) and β max (= 46.3) (Bohnsack et al. 2021). The results suggest a relative change Δk Min of 1,65% to 2.27% with a mean of 2.01 ± 0.15% and Δk Max of 2.69% to 3.69% with a mean of 3.27 ± 0.24% using G dyn . In contrast, a permeability change of 2.77% to 3.88% (mean 3.41 ± 0.27%) for Δk Min and of 4.50% to 6.26% (mean 5.52 ± 0.43%) for Δk Max based on the static shear module is to be expected.

Pore throat size and pore connectivity
The BIB-SEM porosity, as well as the LMI-BIB-SEM observations regarding the filling of almost the entire visible open pores at high injection pressure (Fig. 6), indicate that a significant part, if not the entire BIB-SEM porosity, is connected. This is in good agreement with the similarity of the average BIB-SEM porosity and the average MICP porosity values ( Fig. 8C and D). Moreover, comparing the mean pore size from BIB-SEM (~ 500 nm, Fig. 7) with mean pore throat sizes from MICP (~ 11 nm and ~ 33 nm, Fig. 3) results in pore body to pore throat (aspect) ratios of about 25:1. According to Zhao et al. (2018), who also apply MICP, such ratios for tight carbonates are typically in the range of c. 40:1 to 480:1, with small ratios, indicating increased conductivity to fluid flow compared to large ratios. From our ratio, we propose higher permeabilities for tight carbonates when compared to matrix permeability values of 2.3E−19 to 2.0E−20 m 2 , determined for similar rock types by Hu et al. (2020). However, the ratio only allows for a qualitative statement. For example, Cai et al. (2019) have shown that tortuosity and connectivity, among others, play a significant role in controlling hydraulic flow in porous media. Using the ratio as the sole factor for permeability estimation is therefore not advisable. Overall, the negligible variations in pore volumes, pore throat sizes, and pore size distribution with varying distances to the fault (Figs. 3B and 8C) imply that faulting and fracturing are essentially post-diagenetic.

Microstructure and pore geometry
Submicron interparticle and intercrystalline pores of sub-angular shape are the most abundant types of pores within the matrix of Malm ß limestones. Also, elongated interparticle pores are present, as well as moldic pores, which resulted from the dissolution of the still recognizable fossil remnants but are relatively rare though can be relatively large in size (Fig. 6). The elongated open void spaces are associated with stylolites that are interpreted to have re-opened due to stress release and/or desiccation. The pores at the stylolite interfaces alternate with the accumulation of insoluble, dense minerals mainly at their peaks ( Fig. 5E and F), coinciding with the observations of, e.g., Heap et al. (2014) and Toussaint et al. (2018). These authors assign a complex internal structure of varying thickness to these stylolites, which are interpreted to be a product of the horizontal linkage and vertical coalescence of multiple pressure solution seams (Nenna and Aydin 2011;Toussaint et al. 2018). Hence, our findings confirm the microporous nature of stylolites, which can significantly contribute to the total porosity of low-porous carbonates, at least under unconfined stress conditions. However, due to their elongated shape, it is reasonable to assume that these structures will successively be closed with increasing effective stress (Si et al. 2018). This applies to cases where the principal stress is directed perpendicular to the longitudinal axis of the stylolite and hence to the elongated pores, which are typically arranged parallel to the stylolites. Such an assumption is emphasized by the continued mercury intrusion at the device's maximum pressure (Fig. 3), reflecting either the closure of pores or microfractures or suggesting that some of the open pore space still acts as flow paths for the intruding fluid.

Porosity
A comparison between the various porosity measurement methods shows general conformance among the results with values between 2 and 5% (Fig. 8A), and also with published data for mud-supported limestones (mud-to wackestones) of the same Malm ß stratigraphy from outcrops on the Southern and Northern Franconian Alb (Homuth et al. , 2015Mraz et al. 2018) and the subsurface of the SGMB (Beichel et al. 2014;Böhm et al. 2010;Bohnsack et al. 2020). Neither across the vertical sampled limestone section nor along an individual layer (B1) in a horizontal section, significant changes in the mean matrix porosity were visible (Fig. 8C). This in turn suggests that the petrophysical properties can be considered as rather constant in the investigated section. We think that a missing matrix porosity gradient toward faults is mostly related to both the prevailing stress conditions and the rock's high strength, which promoted dilatant fractures/faults. This, in turn, resulted in localized fluid flow through these structures and prevented the alteration of the tight matrix. Particularly well matching are Φ Gd and Φ He values with only minor method-related deviations ( Fig. 8A and B). In contrast, Φ Ar , Φ Hg , and Φ BIB-SEM values for the same samples are (Fig. 8A) uniformly lower.
Deviations in porosity values are likely related to the capability of the measurement device to register total porosity (e.g., including isolated pores) or only effective porosity (selectively recording connected pores). Also, injection/saturation methods apply different fluids with pore-scale fluid occupancy and connectivity being controlled by surface roughness, intrinsic (wetting) contact angle of the fluid, the wetting state, and the spatial distribution of wettability (Armstrong et al. 2021). Physical limitations in the injection or a negative pressure may prevent the detection of pores below a diameter of 3 nm (MICP) and slightly smaller pores (Archimedes method) (Clarkson et al. 2013;Giesche 2006;Kiula et al. 2014;Okolo et al. 2015;Webb 2001). Full saturation even at the devices' technical limits cannot be ensured with absolute certainty, as air might still be trapped at isolated and dead-end pores. In contrast, the bulk gamma density and the BIB-SEM method are not affected by this problem, with the latter method being only limited by its maximum resolution (≥ 10 nm). However, this should not play a significant role, as it is clear from the pore size distributions of our samples that the majority of pore space is in the 0.1-1 µm range (Figs. 3 and 7B).
Keeping the method drawbacks in mind, our observations indicate a low ratio of isolated pores within the matrix, inferred from only slight differences between Φ Gd and Φ He values (Fig. 8A) and the rare presence of unconnected pores in the LMI-BIB-SEM visualization. The He pycnometry provides the most reliable values for the effective porosity, while the bulk gamma density method detects additional isolated pore space. Regarding the determination of effective porosity for geothermal applications, the Archimedes method provides the most practical method.

Permeability and permeability models
Although fluid flow in tight carbonates is primarily controlled by dissolution and structural features such as faults, fractures, and fracture corridors (Al-Obaid et al. 2005;Dimmen et al. 2017;Litsey et al. 1986;O'Neill 1988;Zeybek and Kuchuk 2002), matrix porosity can be a significant reservoir for fluid recharge. Therefore, different permeability measurement methods and permeability prediction models will be discussed regarding their reliability and limitations.

Permeability
A large spread in air permeabilities (1.9E−18 m 2 to 10.8E−18 m 2 ) most likely results from device-related issues, such as relatively low confining and injection pressures and the missing Klinkenberg correction of the air permeabilities. This also leads to a deviation from the argon permeabilities (1.2E−17 m 2 to 2.1E−17 m 2 , Fig. 9A). Surprisingly, a missing Klinkenberg correction should rather result in increased than in lower permeabilities, as it is the case here. No laminar flow may have been established during the air permeability measurements, despite measurements were only initiated after fluid flow did no longer show large variations. A similar problem is reported by Bohnsack et al. (2020) for mud-supported limestone (mud-to wackestone) samples which are comparable to ours. As this problem seems to be common for low-to ultralow-permeable reservoir rocks, unsteady-state permeability measurement methods (e.g., pulse decay, oscillating pressure, GRI method) on these rock types are frequently preferred over the steady-state methods (Sander et al. 2017). Methodological studies by Boulin et al. (2012), Chenevert and Sharma (1993), and Bertoncello and Honarpour (2013), however, on low-to ultralow-permeable (unconventional gas) reservoir rocks have concluded that the steady-state method is still more reliable than for instance the pulse decay method. Thus, the best practical solution for permeability measurements on low-permeable rocks is still matter of debate (Sander et al. 2017). We consider our data set as reliable as permeability values of a similar range were determined for the same sample set in different laboratories and by different methods.
The missing correlation between porosity and air/argon permeability values (Fig. 9C) suggests that fluid flow through the matrix primarily depends on pore connectivity and pore throat size distribution rather than on the effective porosity. This is in accordance with results of Smodej et al. (2020), showing that the permeability is mostly controlled by the smallest throats. The slightly positive correlation between the argon permeabilities and the r 35 also points to the pore throat size distribution, particularly the r 35 , as the matrix permeability-controlling factor, thereby agreeing with the findings of Saki et al. (2020). However, a statement regarding the importance of the pore throat size distribution as a whole or particular pore throat diameters (e.g., l h max or the r 50 ) for the matrix permeability cannot be made due to a strongly confined range of permeabilities and pore throat diameters (Fig. 9B).
No influence on the matrix permeability was observed in case of the presence of stylolites, regardless of their orientation with respect to the measurement direction. This is in contrast to the findings of Heap et al. (2014) or Hu et al. (2020), who related increased permeabilities to the presence of stylolites. Other authors (e.g., Mehrabi et al. 2016;Rashid et al. 2017;Vandeginste and John 2013) consider stylolites primarily as barriers to fluid flow, depending on their roughness profile (Koehn et al. 2016). This highlights the ambivalent nature of these structural features, as already described by various authors (e.g., Bruna et al. 2019;Burgess and Peter 1985;Korneva et al. 2014;Toussaint et al. 2018). Nevertheless, it is most likely that the stylolite fillings have a lower permeability as the carbonate matrix, as the microstructure of the stylolite filling (Fig. 5D) corresponds to that of porous clays or claystones with permeabilities in the same order of magnitude.
Overall, faulting and fracturing had no noticeably effect on the permeability of the surrounding protolith, at least at the scale of this investigation. We relate this to the post-diagenetic timing of faulting and the limestones overall low permeability, which inhibits intensive fluid-rock interactions, such as karstification or leaching due to the negligible exchange of percolating fluids. Ziauddin and Bize (2007) experimental results support this, as they found limestones of the Khuff formation with similar petrophysical properties compared to the investigated samples to be insensitive to matrix acidizing due to their too low permeabilities.

Permeability models
As illustrated in Appendix 1: Fig 11, most applied permeability prediction models correlate relatively well with measured permeabilities (K Ar and K Air ). However, except for the CT and the Saki Models, the variation of the predicted permeabilities is at least one order of magnitude. As permeabilities measured with argon gas are about twice as high as permeabilities measured with air (Fig. 8A), the difference between these methods is low with a slight trend of K Air toward lower permeability values. Still, some models, such as the CT Model and the Saki Model, show a better correlation to measured permeabilities than others, particularly when applying argon for permeability measurements.
The very good fit of the CT Model (particularly with τ = 1.57) with the measured permeabilities (AAppendix 1 : Fig 11) is primarily related to the involved parameters, namely, the tortuosity τ, the hydraulic radius r i , and the porosity of each segmented pore Φ i , all derived from BIB-SEM analysis. In contrast to other permeability prediction models (except the K-C Model), the CT Model does neither integrate additional empirically determined constants nor porosity values from any of the applied porosity measurement methods. It is most likely for this reason that permeabilities calculated by the CT Model correlate well with the measured permeabilities (Appendix 1: Table 5). The fact that the K-C Model performs worse in terms of the match of predicted with measured permeabilities, despite avoiding the integration of empirical parameters, is most likely related to two facts: firstly, the CT Model has been continuously refined along with improving technical possibilities over the last decades. And secondly, the specific surface area S 0 was measured with the MICP, which might lead to an underestimation at high injection and confining pressures due to the closure of microfractures and stylolites, all resulting in inaccurate permeability predictions. We, therefore, think that the three parameters (τ, r i , and Φ i ) are quite important regarding the assessment of fluid flow within the matrix.
A very good correlation between predicted and measured permeabilities is also achieved by the Saki Model (Appendix 1: Fig 11H). The fact that this model (Saki et al. 2020) applies the smallest pore/throat radius r 35 among other parameters points to the high importance of the r 35 , as a decisive factor in controlling the matrix permeability. Even though the Winland Model (Appendix 1: Fig 11B) also applies the r 35 parameter this results in a worse prediction/measurement match. This suggests that either the calibration samples used by Saki et al. (2020) are more similar to the samples investigated in our study than the ones used by Winland (Comisky et al. 2007;Gunter et al. 2014;Rashid et al. 2015), or rather that this parameter's impact on the Saki equation (Eq. 15) is larger than in the Winland equation (Eq. 9). Nevertheless, the first reason for a worse fit of the model is, in our opinion, the main controlling factor for the deviation of most model-based predictions (K-T, Dastidar, Bohnsack, and the GPPT Models) from measured permeabilities as all of them apply parameters that were empirically determined on the specific data set. Using different petrophysical rock characterization methods for the model calibration might lead to further deviations (e.g., Bohnsack Model). Whether the Saki Model represents the best permeability estimation model for tight carbonates, in general, is therefore questionable and has to be investigated in further studies on other tight carbonates, as this model also incorporates various constants calibrated on their specific data set.
Other models, such as the K-T and the Dastidar Models, which also use pore throat radii at a specific saturation (K-T Model) or capillary pressure (Dastidar Model) correlate significantly worse with measured permeabilities compared to the CT and the Saki Models, suggesting once more that the r 35 parameter is an important matrix permeability-controlling factor, at least in the case of the samples investigated in this study. Although other studies propose that the permeability in tight limestones is primarily controlled by the pore throat radii (Cai et al. 2019;Zhao et al. 2018) or the smallest pore body sizes (Smodej et al. 2020) rather than by the entire effective porosity range, this does not contradict our results. Instead, our results are rather complementary to these findings and specify that one particular pore throat radius, namely, the r 35 , has a particularly strong control on fluid flow in the matrix.
Implications for matrix permeability-controlling factors Relationships between the various petrophysical parameters lead to the following conclusions. As indicated by the overestimation of predicted permeabilities that use the more accurate He porosity values, microporosity does not have a significant impact on the rock's matrix permeability. As the He porosity, which is higher compared to the Ar, MICP, and BIB-SEM porosities due to the registration of microporosity, is integrated in permeability calculations as a multiplication factor (K-T model), as an exponential factor (Winland, Dastidar, K-C, and Bohnsack Models), or as logarithmic factor (Saki and GPPT Models), the predicted permeabilities are also elevated. Resulting permeability predictions deviating from measured permeabilities, therefore, implies a negligible effect of microporosity, but a higher importance of the pore throat size distribution and the pore network connectivity, conforming to the conclusions by Lala and El-Sayed (2017) and Philipp et al. (2017). In particular, the pore throat radius at a mercury saturation of 35% (r 35 ) seems to be the parameter controlling most to the fluid flow in the pore network, as both best-performing models (Winland and Saki) integrate this measure. Also, a correlation between argon permeabilities and the r 35 could be found, though without a clear trend due to the confined value range (Fig. 9B).
We can further conclude that most permeability prediction models deviate from the measured permeabilities due to the involvement of empirically determined variables that were calibrated on rock sample sets of slightly different lithologies, derived from locations that experienced a different diagenetic history compared to our samples. Also misleading appears the application of effective porosity as the sole physical input parameter, as this is not the main permeability-controlling factor according to Cai et al. (2019) and Zhao et al. (2018). Fluid flow in pore networks also strongly depends on the wettability of the permeating fluid/gas, on the spatial distribution of the wettability, which is related to the type of minerals present in the sample, as well as on the surface roughness (Armstrong et al. 2021). Both the spatial wettability distribution and the surface pore roughness are difficult to quantify and were not available in this study. This fact and the complex interrelation of these factors impede reliable permeability predictions, as most permeability prediction models applied in this study-except for the Capillary Tube and the Bohnsack Model-try to approximate the permeability based on MICP-derived pore throat diameter distributions, which in turn depend on the above mentioned permeability-controlling factors proposed by Armstrong et al. (2021). As these factors also concern the permeability measurement itself, the permeability prediction derived from measurements that use different permeating fluids will always be prone to inaccuracies unless detailed information on the wettability distribution, the fluid's wettability, and the surface roughness of the sample is considered.
Model rating To quantify the performance of the various models for permeability estimation, the Mean Residual Square Error (MRSE) was calculated for each model and each porosity-determining method (Appendix 1: Table 5) as illustrated in Fig. 10. We will only compare our rankings of the models' performances in the context of argon permeabilities to the rankings of other authors from here on, as the argon permeabilities proved a greater reliability compared to air permeabilities.
Based on the MRSE, a ranking for each method and an overall comparative ranking is determined (Table 1). The CT Model could not be taken into account as it does not use a bulk porosity method. Relating the permeability predictions to the measured air permeabilities K Air (Fig. 10A) provides satisfying results for the K-T (overall rank 1), the Saki (overall rank 2), and the Bohnsack Models (overall rank 3). A slightly different rating results from models using Φ Ar and Φ Hg , where both the K-T (rank 1) and the Bohnsack Models (rank 1) perform equally well, and the Saki Model (rank 3) a little worse (Fig. 10b).
While the comparison to argon permeabilities K Ar (Fig. 10B) shows best correlations with the Saki Model (overall rank 1), both the Winland (overall rank 2) and the K-T Models (overall rank 3) achieve slightly lower rankings. Again, a slightly different ranking order results when only models applying Φ Ar and Φ Hg are considered. The Saki Model (rank 1) performs best and the Winland (rank 2) as well as the Bohnsack and the K-T Models (both rank 3) show good correlations with the measured argon permeability. The CT Model is always among the top four prediction models. Hence, we recommend using either the Saki or the CT Model for estimating permeabilities of similar sample lithologies.
The comparison of our permeability estimation model ranking to rankings of Comisky et al. (2007) and Rashid et al. (2015) shows some similarities, although the spectrum of their models differs from ours. Comisky et al. (2007) for instance found a good correlation between measured and modeled permeabilities based on the K-T method (their rank 2) and the Winland method (their rank 5), whereas these models achieved ranks 3 an 2 in our comparison.
Those models that were applied by both our study and that of Rashid et al. (2015) are the Winland (their rank 4), the K-T (their rank 9), and the Dastidar Models (their rank 12). Although a comparison of their ranking results to our ranking is hardly representative, their relative ranking order differs from ours, as the K-T Model performed very well in the case of our data set. Nevertheless, even the best-performing permeability estimation model (generic model) in the ranking of Rashid et al. (2015) still has a mean residual square error of 0.402 and therefore shows a poorer correlation to measured permeabilities than the CT and the Saki Models.  Table 5 Implications for a "best practice" in the petrophysical analyses of tight carbonates We conclude from our data set that porosity data for permeability models are best acquired by mercury intrusion porosimetry, or even better, LMI-BIB-SEM, as this method provides not only reliable information on the porosity (the connected pore volume relative to the bulk rock volume) but, compared to other methods (e.g., NMR, SEM, X-ray CT), also on a larger range of pore throat sizes (Wu et al. 2019). Pore throat sizes (particularly the r 35 ) exercise primary control over matrix permeability and are therefore an indispensable prerequisite for a reliable permeability estimation. Although a larger pore size range is accessed and measured through Helium pycnometry, Mercury intrusion porosimetry offers a wider range of applications as also the permeability can be estimated quite accurately, even from small sample volumes, such as drill cuttings. Hence, even when the availability of sample material is limited, the rock's matrix permeability can be estimated, which is particularly useful for subsequent modeling purposes. However, special care should be taken regarding the presence of microfractures induced by sample handling. Ideally one should investigate the sample after injection, for example by employing LMI-BIB-SEM, as this method is well suited for the detection of microfractures and injection-related artifacts.

Mechanical properties and the relation of petrophysical properties to buried samples
All samples that were investigated in terms of their rigidity revealed similar values in the range of 65 to 86 GPa (Young's Modulus, Table 1). Geomechanical measurements on rocks from the same stratigraphic units (Malm α and β) and of comparable lithology (mud-to wackestone), but buried 240-5200 m TVD in the SGMB were carried out by Potten (2020) showing values of a similar range (70 to 74 GPa), confirming the comparability of outcrop analog and buried tight carbonates, at least in this particular case. The rock's stiffness depends mostly on the pore type, pore geometry, and the total porosity (Eberli et al. 2003;Li et al. 2018;Weger et al. 2004). These parameters, together with pre-existing zones of weakness, such as stylolites (e.g., Agosta et al. 2015;Antonellini et al. 2008Antonellini et al. , 2014Graham et al. 2003;Micarelli et al. 2005;Tondi 2007) and the in situ effective stress conditions control the deformation mode and the resulting fault zone structure and consequently its hydraulic capacity to conduct fluids (e.g., Antonellini et al. 2008;Michie 2015;Sagi et al. 2016). The majority of the samples' pores were found to be of the submicron interparticle, intercrystalline, and occasionally also of the moldic small mesopore type, and of (tri-)angular pore geometry. All determined properties in combination with the low total porosity explain the very high mechanical strength of the investigated rocks (Eberli et al. 2003;Potten 2020).
Relating the mechanically very competent limestones from Frankenalb outcrop analogs to reservoir depths, which in case of the SGMB would be > 2000 m TVD toward the south of Munich (Mraz 2019), a relative volume change of only 0.06% to 0.08% or 0.10% to 0.14% is calculated (see Appendix 1: Table 5), depending on whether G dyn or G stat (Eqs. 17 and 19) are used for calculations. In contrast, Bohnsack et al. (2021) determined a relative porosity change of ~ 2% when buried ~ 2 km TVD, hence one order of magnitude larger. This large discrepancy is most likely related to the fact that porosities and geomechanical parameters of their investigated reservoir rocks (Grainstones from the Malm ζ 4-5) differ significantly from our samples (higher porosity, but lower rock stiffness/strength), which strongly increases their sensibility to effective stress changes. This highlights the importance of pore geometry and pore connectivity in relation to geomechanical properties and the stress sensitivity, a dependency has been thoroughly investigated by various authors (e.g., Zoback and Byerlee 1975;Pei et al. 2014;Xu et al. 2018a, b).
Comparing our relative permeability changes Δk Min (1.65% to 2.27%) and Δk Max (2.69% to 3.69%) calculated with G dyn or Δk Min (2.77% to 3.88%) and Δk Max (4.50% to 6.26%) calculated with G stat to those of Bohnsack et al. (2021) (33.0% to 56.7%) again shows that for our rocks, the permeability change is about one order of magnitude smaller than the values determined by Bohnsack et al. (2021) or other authors investigating similar rock types (e.g., Bakhtiari et al. 2011;David et al. 1994;Moosavi et al. 2014;Hu et al. 2020). As the stress sensitivity of permeability is a function of porosity (Cheng 2016) and mainly dependent on the aspect ratio, type and size of pores, and pore throats (Rashid et al. 2017) as well as the presence of fracture-like pores or cracks (Bohnsack et al. 2021), permeability change with increasing effective stress is typically higher than porosity change. The permeability reduction is thereby caused by the size reduction of voids, pore closure, pore throat collapse, and/or clogging of pores (Selvadurai and Głowacki 2008). Low-porous carbonate rocks with nano-intercrystalline pore types and high aspect ratio pores should therefore be much more sensible to effective stress-induced permeability reduction (c.f., Rashid et al. 2017).
As the rocks investigated in this study are already strongly lithified and consequently have very high elastic moduli, low porosities, and low pore throat diameters in relation to the samples investigated by Bohnsack et al. (2021), they are much more resistant to effective stressinduced deformation. Together with the comparatively low aspect ratios for these types of rocks found by Zhao et al. (2018), the significantly lower sensitivity to effective stress associated with low porosity/permeability changes with increasing effective stress can be explained. It has to be considered that the static Shear Modulus (G stat ) was calculated based on an empirical relationship identified by Bastos et al. (1998) and these values are therefore prone to inaccuracies due to lithological and petrophysical differences. Studies by  as well as Pei et al. (2014) on rocks lithologically and stratigraphically similar to those of our study showed that other than the effective stress, rather the prevailing temperature at reservoir level (≥ 60 °C) appears to control the permeability reduction due to thermal expansion. As the rocks investigated in this study are all low porous and primarily composed of rock matrix, thermal expansion with increasing depth and temperature might have a much larger impact on the petrophysical properties of these rocks at depth than the effective stress.
Consequently, stress conditions at reservoir level should hardly influence the petrophysical properties of the limestone matrix. Rock matrix petrophysical properties measured at surface conditions should therefore be well relatable to properties at reservoir level. Lower porosities measured by Bohnsack et al. (2021) at 25 MPa confining stress (~ 2.0%) compared to the porosities measured at atmospheric conditions (~ 3.3%) by  and by us (2.9-3.2%) of the same stratigraphic unit likely relate to changes in lithology (mud-supported limestones), different degrees of dolomitization related to regional heterogeneities of the depositional system, and variation in diagenetic alteration processes (Bohnsack et al. 2020). As the differences in porosities are very low, these factors are thought to play a minor role in their influence on the petrophysical properties, at least in the case of Malm β carbonates. Also, as according to our findings, the pore throat distribution is controlling the matrix permeability rather than the effective porosity, these small differences in porosity should be negligible and are considered as intrinsic sample heterogeneities. Thus, we are confident that our results reflect a good estimation of in situ petrophysical properties due to the rocks' high stiffness.

Conclusion
In this study, petrophysical properties and microstructures of tight carbonates were investigated with various methods, aiming for characterizing the matrix properties of the Upper Jurassic (Malm ß) limestones in Northern Bavaria for modeling purposes and deriving correlations between the different measurements using various models. Furthermore, their interdependence of the protolith's porosity with stylolitization and normal fault distance were examined. The overall average porosity measured by MICP, He, GD, BIB-SEM, and Archimedes methods is 3.1 ± 1.0% and the overall average gas permeability is 1.4E−17 m 2 , defining them as tight carbonates. Nevertheless, it has to be considered that the effective porosity, e.g., for geothermal applications is expected to be much lower than the measured values as the percolating fluids are thermal waters or brines, for which the microporosity can be neglected in terms of its importance for subsurface fluid transport within the matrix. No significant impact of stylolites on the permeability and porosity under confined conditions was observed. Furthermore, no improved petrophysical properties in a preferred direction were recorded, indicating a homogeneous, isotropic behavior of the rock. Normal faulting had no observable impact on the matrix's poro-perm properties at the meter scale. We, therefore, assume postdiagenetic fracturing and faulting, as otherwise faulting-related processes would have most likely altered the petrophysical properties of the rock matrix. Few leached pores were observed, but with negligible effect on porosity and permeability. The majority of the porosity is in the submicron range. Two models, the CT Model and the Saki Model for estimating permeability from porosity and pore size distribution data, obtained from either the MICP or the BIB-SEM method, gave best results compared to measured permeabilities (air and argon) for the examined rock types. The CT Model applies a tortuosity factor of 1.57 fitted from BIB-SEM analyses, while the Saki Model integrates the pore throat size at 35% mercury saturation. We found that the pore throat diameter is a main factor controlling fluid flow in the rock matrix rather than the effective porosity or microporosity. The application of other models developed by Katz-Thompson (K-T), Winland, and Bohnsack et al. using Ar and MICP porosities also results in reasonable permeability estimations. Due to the relative high rock stiffness of ~ 73 GPa for the dynamic Young's Modulus, ~ 53 GPa for the dynamic Bulk Modulus, 27 GPa for the dynamic and 16 GPa for the calculated static Shear Modulus our findings can be reliably transferred to reservoir conditions (even at up to 2000 m depth), as the volumetric change of porosity (maximum 0.12%) and of permeability (maximum 5.52%) due to effective stress at this depth are negligibly small.*Poisson's number calculated from samples, where the p-wave (V P ) and the s-wave (V S ) velocity were available     Table 5 Relative porosity (ΔΦ/Φ) and permeability (Δk/k) changes for the investigated rock at a depth of ~ 2 kmTVD with an assumed vertical effective stress of dp ≈ 25 MPa (according to Drews et al. 2020), calculated after Eqs. (18) and (20)