Articles | Volume 12, issue 8
Solid Earth, 12, 1777–1799, 2021
Solid Earth, 12, 1777–1799, 2021

Research article 11 Aug 2021

Research article | 11 Aug 2021

3D crustal stress state of Germany according to a data-calibrated geomechanical model

3D crustal stress state of Germany according to a data-calibrated geomechanical model
Steffen Ahlers1, Andreas Henk1, Tobias Hergert1, Karsten Reiter1, Birgit Müller2, Luisa Röckel2, Oliver Heidbach3,4, Sophia Morawietz3,4, Magdalena Scheck-Wenderoth5,6, and Denis Anikiev5 Steffen Ahlers et al.
  • 1Engineering Geology, Institute of Applied Geosciences, TU Darmstadt, 64287 Darmstadt, Germany
  • 2Technical Petrophysics, Institute of Applied Geosciences, KIT, 76131 Karlsruhe, Germany
  • 3Seismic Hazard and Risk Dynamics, GFZ German Research Centre for Geosciences, 14473 Potsdam, Germany
  • 4Institute for Applied Geosciences, TU Berlin, 10587 Berlin, Germany
  • 5Basin Modelling, GFZ German Research Centre for Geosciences, 14473 Potsdam, Germany
  • 6Department of Geology, Geochemistry of Petroleum and Coal, Faculty of Georesources and Material Engineering, RWTH Aachen University, Aachen, Germany

Correspondence: Steffen Ahlers (


The contemporary stress state in the upper crust is of great interest for geotechnical applications and basic research alike. However, our knowledge of the crustal stress field from the data perspective is limited. For Germany basically two datasets are available: orientations of the maximum horizontal stress (SHmax) and the stress regime as part of the World Stress Map (WSM) database as well as a complementary compilation of stress magnitude data of Germany and adjacent regions. However, these datasets only provide pointwise, incomplete and heterogeneous information of the 3D stress tensor. Here, we present a geomechanical–numerical model that provides a continuous description of the contemporary 3D crustal stress state on a regional scale for Germany. The model covers an area of about 1000×1250 km2 and extends to a depth of 100 km containing seven units, with specific material properties (density and elastic rock properties) and laterally varying thicknesses: a sedimentary unit, four different units of the upper crust, the lower crust and the lithospheric mantle. The model is calibrated by the two datasets to achieve a best-fit regarding the SHmax orientations and the minimum horizontal stress magnitudes (Shmin). The modeled orientations of SHmax are almost entirely within the uncertainties of the WSM data used and the Shmin magnitudes fit to various datasets well. Only the SHmax magnitudes show locally significant deviations, primarily indicating values that are too low in the lower part of the model. The model is open for further refinements regarding model geometry, e.g., additional layers with laterally varying material properties, and incorporation of future stress measurements. In addition, it can provide the initial stress state for local geomechanical models with a higher resolution.

1 Introduction

Knowledge about the stress state in the upper crust is of great importance for many economic and scientific questions. Some examples are wellbore stability (Bell, 2003; Kristiansen, 2004), operation and stimulation of hydrocarbon and geothermal reservoirs (Altmann et al., 2014; Azzola et al., 2019; Henk, 2009; Smart et al., 2014), slip and dilation tendency of existing faults and fractures (Hettema, 2020; Konstantinovskaya et al., 2012), underground mining (Brady and Brown, 2004) and deep tunneling (Diederichs et al., 2004). Furthermore, it plays a decisive role in the search for a disposal site for high-level radioactive waste, since it is crucial for the short and long-term safety of a possible repository (StandAG, 2017; Nagra, 2008; BGR, 2015). For all these applications the contemporary stress state is a key parameter, and thus the quantification of the complete 3D stress tensor is essential.

However, from the data perspective, our knowledge of the stress state in western central Europe is limited in particular regarding stress magnitude information. Public stress information is provided by the World Stress Map (WSM) project, which supplies a global database of the orientation of the maximum horizontal stress (SHmax) and the stress regime (Heidbach et al., 2016) and by a compilation of stress magnitude data for Germany and adjacent regions of Morawietz and Reiter (2020). However, these two datasets contain only pointwise information, which is incomplete as only a subset of the stress tensor components is provided and their spatial distribution is sparse and irregular (Fig. 1a).

To provide a continuous description of the 3D stress tensor in the upper crust on a regional scale, we developed the first 3D geomechanical model covering Germany (Fig. 1). Our model comprises seven units, with specific material properties and laterally varying thicknesses: a sedimentary unit, four different units of the upper crust, the lower crust and the lithospheric mantle. The finite element method (FEM) is used to solve the partial differential equation which describes the equilibrium of body and surface forces within an inhomogeneous medium. Our input parameters are density and elastic material properties (Young's modulus and Poisson's ratio). The model is calibrated using appropriate initial conditions and displacement boundary conditions to find a best fit with respect to the stress orientation and magnitude datasets described above. This modeling approach has been used for a wide range of scales and different tectonic settings (Buchmann and Connolly, 2007; Heidbach et al., 2014; Hergert and Heidbach, 2011; Hergert et al., 2015; Reiter and Heidbach, 2014).

2 Fundamentals and state of the art

2.1 Geology and tectonic setting of the study area

The crustal and lithospheric structure in the model domain reflects the complex geodynamic evolution of central Europe since Precambrian times (McCann, 2008; Meschede and Warr, 2019) (Fig. 1c and d). The north-eastern part of the study area belongs to the cratonic unit of Baltica and, more specifically, the East European Craton (EEC). This unit consists mainly of high-grade magmatic and metamorphic rocks of Precambrian and early Paleozoic age. Crustal thickness in the area is about 50 km, and a thick mantle lithosphere down to depths of 200 km has been observed (Mazur et al., 2015). The EEC is separated from the Avalonia microplate to the south-west by the Tornquist Suture and the Thor Suture, respectively (e.g Linnemann et al., 2008, and various references therein). At this boundary, a sharp transition to the significantly thinner crustal and lithospheric thicknesses typical for Paleozoic and Mesozoic Europe can be observed (Ziegler and Dèzes, 2006). Western Baltica and eastern Avalonia got into contact during closure of the Tornquist Ocean during Ordovician to Silurian times and the Caledonian orogeny, respectively. At this stage, Laurussia (composed of Laurentia, Baltica and Avalonia) was formed, whose continental crust makes up the northern and eastern part of the study area.

The central part of the model domain comprises the southern part of Avalonia as well as Armorica – microplates and terrane assemblages – which collided during the Variscan orogeny in Late Paleozoic times (Franke, 1989, 2006; Meschede and Warr, 2019). The low-grade metamorphic rocks of the Rhenohercynian Zone (sensu Kossmat, 1927) represent passive margin sediments which were deposited on thinned crust of south-eastern Avalonia. South-eastward-directed subduction and closure of the Rheic Ocean led to the formation of an active margin at the northern rim of Armorica, which nowadays comprises the medium-grade metamorphic and magmatic rocks of the so-called Mid German Crystalline High (Oncken, 1997). Further to the south, the Saxothuringian and Moldanubian Zone represent the remnants of the internal zone of the Variscan orogen with medium- to high-grade metamorphic rocks and abundant granitoids presently exposed at surface. Crustal thickness in this part of the model domain and outside the areas affected by Cenozoic rifting and mountain building is typically on the order of 30 km (Ziegler and Dèzes, 2006). The late Cretaceous to Paleogene evolution was influenced by NE-directed Africa–Iberia–Europe convergence which led to intraplate contraction and inversion of NW–SE-striking structural elements (Kley and Voigt, 2008). The final stage of this phase coincides with W–E to NW–SE directed extension and the onset of rifting in the Upper Rhine Graben (URG) and Eger Graben, among others (Kley et al., 2008).

The southernmost parts of the study area are located in the so-called ALCAPA (Alps–Carpathians–Pannonian) unit or terrane (e.g., Brückl et al., 2010; Schmid et al., 2004). Its geodynamic evolution is closely related to the collision between Europe and the Adriatic–Apulian microplate leading to the Alpine orogen. Since Eocene times, its northern foreland has been characterized by N–S to NW–SE-directed compression and thrusting, respectively (Reicherter et al., 2008).

Figure 1Maps of western central Europe. The red polygon indicates the model area. (a) Overview of calibration data used. Color-coded lines indicate the orientation of SHmax and the stress regime of the WSM (Heidbach et al., 2016) and additional data of Levi et al. (2019). Grey dots show the positions of stress magnitude data of Morawietz and Reiter (2020). (b) Topography and mean SHmax orientations on a regular 0.5 grid derived from the WSM. Each grid point requires at least 10 data points within a fixed search radius of 200 km (details in Sect. 4.1). The topography is based on Smith and Sandwell (1997). (c) Tectonic framework of the model area based on Asch (2005) and Kley and Voigt (2008). EG – Eger Graben, FL – Franconian Line, LRB – Lower Rhine Basin, M – Massif, URG – Upper Rhine Graben. (d) Overview of the crustal units in central Europe (modified after Kroner et al., 2010; Brückl et al., 2010). Black titles show tectonic units and white titles sutures and Variscan units. Coastlines and borders used in the figures are based on the Global Self-consistent Hierarchical High-resolution Geography (GSHHG) of Wessel and Smith (1996).

2.2 Basics of the crustal stress state

The stress state at a given point can be described by a second rank tensor (Fig. 2a) with the pascal (1 Pa = N m−2) as the basic unit. Due to its symmetry properties, only six out of nine components are independent of each other (e.g., Jaeger et al., 2011). In the principal axis system, the off-diagonal components vanish and the remaining three components are the principal stresses σ1, σ2 and σ3. Their orientations and magnitudes describe the absolute stress state (Fig. 2b). Assuming that the vertical stress (SV) is one of the three principal stresses (Fig. 2c), the orientation of this so-called reduced stress tensor is determined by the orientation of SHmax. Given that SV can be approximated by depth and density of the overburden, the remaining unknowns are the magnitudes of SHmax and the minimum horizontal stress (Shmin).

Figure 2(a) The nine components of the stress tensor define the stress state at an arbitrary point and enable to compute the stress vector on any surface through that point. To describe the stress tensor components an infinitely small cube with uniform surfaces is used. (b) Due to the conservation of momentum, the stress tensor is symmetric and thus a coordinate system exists where shear stresses vanish along the faces of the cube. In this principal axis system, the remaining three stresses are the principal stresses. (c) Assuming that the vertical stress (the overburden) in the Earth crust SV=gρcotz is a principal stress (g is gravitational acceleration, ρ is the rock density, z is depth), the two horizontal stresses Shmin and SHmax, the minimum and maximum horizontal stress, respectively, are principal stresses as well. This so-called reduced stress tensor is fully determined with four components: the SHmax orientation and the magnitudes of SV, Shmin and SHmax (Heidbach et al., 2018).


The stress state of the continental crust is influenced by stress sources on different scales from several meters up to several thousand kilometers: first-order stress sources (>100 km) related to plate boundary forces, e.g., ridge push or slab pull, second-order stress sources ( 100 km) related to large volume forces, e.g., lithospheric flexure due to mountain ranges or deglaciation and third-order stress sources (<100 km) related to local density or stiffness contrasts, e.g., faults or diapirs. Second- and third-order stress sources are able to disturb the overall stress orientation trend from regional through local to reservoir scale (Heidbach et al., 2007).

2.3 Data compilation of stress tensor components

The orientation of the stress tensor in the Earth's crust is provided by the WSM database, which is a global compilation providing data on the SHmax orientation and the stress regime (Heidbach et al., 2016). This stress information is derived from a variety of methods, primarily earthquake focal mechanism solutions, borehole breakouts and drilling-induced tensile fractures (from borehole image or multi-arm caliper log data), in situ stress measurements (overcoring or hydraulic fracturing) and geologic indicators, such as fault slip and volcanic vent alignment (Amadei and Stephansson, 1997; Ljunggren et al., 2003; Schmitt et al., 2012). The stress information in the WSM database is compiled in a standardized format and quality-ranked for reliability and comparability on a global scale (Heidbach et al., 2010; Zoback, 1992). For Germany and adjacent regions, the SHmax orientations have been re-evaluated recently (Reiter et al., 2016, 2015). The new data have been integrated in the latest WSM database release (Heidbach et al., 2016). For stress magnitude data, Morawietz et al. (2020) published a publicly accessible database with 568 data records including a quality assessment of the data for Germany and adjacent regions. These two datasets (Fig. 1a) – the SHmax orientation of the WSM with some additional data of Levi et al. (2019) and the stress magnitude database (Morawietz and Reiter, 2020) – are used to calibrate the geomechanical model.

2.4 Previous models

Modeling the contemporary crustal stress in western central Europe has been addressed by various authors since the mid-1980s. However, except the model of Buchmann and Connolly (2007), who provide a 3D model of the broader URG, all models are 2D and with a strong emphasis on the SHmax orientation and little regarding the stress magnitudes. Table 1 gives a short overview of their key technical characteristics. If several model versions are published by one author, the most current one is listed. In general, different plastic and elastic material properties have been tested so far and also various boundary conditions have been applied. For a detailed overview we refer to, e.g., Cacace (2008), Heidbach et al. (2007) or Jarosiński et al. (2006).

The results indicate different main factors influencing the contemporary stress state. The majority of the studies have found lateral stiffness contrasts in the lithosphere, such as the Bohemian Massif, the Elbe Fault Zone or the Avalonia–EEC boundary (Cacace, 2008; Grünthal and Stromeyer, 1994; Jarosiński et al., 2006; Marotta et al., 2002) and isostatic effects (Bada et al., 2001; Kaiser et al., 2005; Jarosiński et al., 2006) to be the main cause of stress perturbations. In addition, faults or fault zones are held responsible for third- or second-order stresses, respectively (Jarosiński et al., 2006; Kaiser et al., 2005). An indirect influence of the depth of the asthenosphere lithosphere boundary (LAB) and the resulting temperature contrasts and changed mechanical properties are described by Cacace (2008).

Table 1Overview of regional-scale stress models within the model area. If several model versions are published, the most current one is listed. The (X) is used for Buchmann and Connolly (2007), since the boundary conditions are not derived from the plate boundary forces but still represent them.

Download Print Version | Download XLSX

3 Model setup

3.1 Conceptual modeling approach

To model the contemporary 3D stress field of the upper crust, we assume linear elasticity and neglect thermal stresses and pore pressure effects. With these assumptions, the partial differential equation of the equilibrium of forces has to be solved (Jaeger et al., 2007). The two contributing forces are volume forces from gravitational acceleration and surfaces forces that are mainly attributed to plate tectonics. The latter are key drivers for the tectonic stress that we observe, and they are parameterized with displacement boundary conditions that are chosen in a way that the resulting stresses deliver a best fit with respect to the model-independent stress data. Although these displacement boundary conditions are mainly representing the tectonic stresses, they are not derived from these. Accordingly, our results do not allow any conclusions regarding the sources of crustal stress in the model area. This process is called model calibration, which can also be used to estimate model uncertainties by means of standard deviation (Ziegler et al., 2016; Ziegler and Heidbach, 2020). The technical procedure is presented in Fig. 3 with a schematic general workflow. The individual text boxes are color-coded indicating the four major steps.

The model geometry reflects the contemporary distribution of rock properties such as density and stiffness and Poisson's ratio. An appropriate initial stress equilibrates the gravitational stresses and resembles a reference stress state (Fischer and Henk, 2013; Hergert et al., 2015; Reiter and Heidbach, 2014). The orientations of the lateral model boundaries where the displacement boundary conditions are applied are chosen in such a way that the mean SHmax orientation (Fig. 1b) is perpendicular or parallel to them.

For the solution of the partial differential equation of the equilibrium of forces, we use the FEM to estimate an approximated numerical solution. The FEM is appropriate as it allows discretizing complex geometries with unstructured meshes. The commercial FEM software package Abaqus v2019 is used. For post-processing we are using Tecplot 360 enhanced with the GeoStress add-on (Stromeyer and Heidbach, 2017). For the construction and discretization of the 3D model geometry, GOCAD and HyperMesh are used.

Figure 3General workflow of 3D geomechanical–numerical modeling. White boxes: assembly of model geometry and rock properties. Left figure: 3D view of the discretized model volume. Grey boxes: initial stress field and kinematic boundary conditions, gravity load and numerical solution. An appropriate initial stress state and kinematic boundary conditions are determined and applied as well as gravity load. Right figure: discretized model volume including boundary conditions used. The partial differential equation of the equilibrium of forces in 3D is solved using the FEM (σij stress tensor, xj Cartesian coordinates, ρ density, and Xi body forces). Orange boxes: model results are calibrated against model-independent observations. Yellow box: once the fit to the model-independent observations is acceptable, i.e., within their uncertainties, an interpretation and analysis of the model results can be performed.


3.2 Model geometry

The model geometry extends over 1250 km in east–west direction from eastern Poland to western France and by 1000 km in north–south direction from southern Scandinavia to northern Italy, covering an area of about 1.25 million km2. This area was chosen with regard to the orientation of SHmax to simplify the definition of boundary conditions later on and with regard to important crustal structures which may affect the recent stress field in Germany, e.g., the Bohemian Massif, the Avalonia–EEC suture and the European Cenozoic Rift System. Additionally, model boundaries are selected distal to the German border to avoid possible boundary effects in the area of main interest (Fig. 1).

The model geometry contains seven units: a sedimentary cover, the upper crust subdivided into four units, the lower crust and parts of the lithospheric mantle. The units are bounded by five surfaces: the topography, the top of the crystalline basement, the top of the lower crust, the Mohorovičić discontinuity (Moho) and the model base at 100 km depth. The bottom of the model is thus not defined as the LAB, and therefore the thickness of the lithospheric mantle can deviate from its real thickness. The Moho was chosen as the deepest surface since almost all calibration data are from above, and also the depth interval of greatest interest is the upper 10 km of the crust. Although deeper structures may exert a long-wavelength effect on the stress state in the upper crust, we expect that the primary contributions to the stress field are captured by the considered interfaces. The upper crust is laterally sub-divided into four parts: the EEC, Avalonia, the Armorican Terrane Assemblage and the ALCAPA unit referring to the tectonic units displayed in Fig. 1d.

Figure 4Database and depth maps of the top of the crystalline basement, the top of the lower crust and the Mohorovičić discontinuity. We used the color map “roma” of Crameri (2021). Sources of used data are as follows. Crystalline basement: Anikiev et al. (2019) (dark green), Diebold et al. (1991) (brown), GeORG-Projektteam (2013) (light blue), Geothermieatlas Bayern (2004) (pink), Hurtig et al. (1992) (blue green), Kirsch et al. (2017) (red), Korsch and Schäfer (1995) (light green), Lindner et al. (2004) (grey), Maystrenko and Scheck-Wenderoth (2013) (dark blue), Reinhold (2005) (orange), Rupf and Nitsch (2008) (cyan), Sommaruga (1999) (yellow), Tašárová et al. (2016) (light red), and black profiles (Behr et al., 1994; Bokelmann and Bianchi, 2018; Cazes et al., 1985; Freeman and Mueller, 1992; Grad et al., 2009a; Heinrichs et al., 1994; Hirschmann, 1996; Janik et al., 2011; Meschede and Warr, 2019; Oncken et al., 2000; Reinhold, 2005; Schintgen, 2015; Wenzel and Brun, 1991). A labeled, larger size map of this database is available in the Supplement. Top of the lower crust: Anikiev et al. (2019) (dark green), Maystrenko and Scheck-Wenderoth (2013) (dark blue), Tašárová et al. (2016) (light red), and Valasek and Mueller (1997) (cyan). Mohorovičić discontinuity: Anikiev et al. (2019) (dark green), Grad et al. (2009b) (red), Maystrenko and Scheck-Wenderoth (2013) (dark blue), and Wagner et al. (2012) (pink). Coastlines and borders used in the figures are based on the Global Self-consistent Hierarchical High-resolution Geography (GSHHG) of Wessel and Smith (1996).

Figure 4 shows the depth maps of the top of the crystalline basement, the top of the lower crust and the Moho with the corresponding database used. The model is mainly based on three existing models: the 3D Deutschland model (Anikiev et al., 2019), the Central European Basin model (Maystrenko and Scheck-Wenderoth, 2013) and the central Europe model of Tašárová et al. (2016).

The key challenge was the construction of the top of the crystalline basement. In all the models used and also in most other datasets, the base of the sedimentary layer is defined as the top of the basement regardless of whether the basement consists of crystalline or low-grade metamorphic rocks. This is an assumption which is not sufficient to represent the stiffness contrast correctly. The main reason for this assumption is the lack of data due to the usually great depths and the lack of economic interest in these units. Especially in the Rhenohercynian and Saxothuringian Zone (Fig. 1c), only a few seismic profiles exist from research projects like DEKORP (Meissner and Bortfeld, 1990), EGT (Freeman and Mueller, 1992) or ZENTROSEIS (Bormann et al., 1986). Despite the uncertainties due to this poor amount of data, the use of the sediment–crystalline boundary is necessary for a geomechanical–numerical model, because of the strongly different mechanical properties. An extreme example within our model area is the western part of the Rhenohercynian Zone. Here the basement is outcropping, e.g., in the Rhenish Massif, but the top of the crystalline crust is suspected to be at about 20 km depth (Schintgen, 2015; Oncken et al., 2000). Therefore, in those areas where the definition of the basement does not correlate to the top of the crystalline basement, we constructed this surface to obtain a mechanically uniform surface; data used are shown in Fig. 4. The boundaries between the Variscan basement units are simplified as vertical due to the poor knowledge.

3.3 Model discretization

Our final mesh shown in Fig. 5 comprises 1.32 million hexahedral elements with a lateral homogenous resolution of approximately 6×6 km2. The vertical resolution decreases with depths from 800 m near the surface up to 7500 m at the base of the model. An exception is the uppermost element layer, which is only 50 m thin to reduce the impact of free surface effects in the uppermost units. Due to the complex geometry of our model we decided not to use the common approach in the upper units in which each unit is meshed individually. Only the mantle and the crust are meshed as whole. Then we use the tool ApplePy (Ziegler et al., 2019) to assign each finite element to the respective subunits and the appropriate rock properties.

Figure 5Six different views of the discretized model showing the internal model structure. The sedimentary unit is colored in yellow, the upper crust in different red shades regarding to different tectonic units (Fig. 1c), the lower crust in light grey and the lithospheric mantle in dark grey. The dimension of the model is 1000×1250×100 km3 comprising 1.32 million hexahedral elements. ALCAPA – Alps–Carpathians–Pannonian, ATA – Armorican Terrane Assemblage, EEC – East European Craton.


3.4 Rock properties

The material properties used in the model and corresponding references are shown in Table 2. The assignment of mean rock properties to the sediment unit is a difficult task, due to the large number of different rock types represented by unconsolidated rocks, claystones, sandstones, salt or limestones. Therefore, the values are approximate mean values. For the upper crust we applied a different density for each tectonic unit in the range of 2750 to 2820 kg m−3. Young's modulus and Poisson's ratio are values for granodiorite as the characteristic rock of the upper crust. For the lower crust density we use the results of Maystrenko and Scheck-Wenderoth (2013) and Tašárová et al. (2016), but since, unlike them, we have only one uniform unit, we use an average of 3000 kg m−3. Young's modulus and Poisson's ratio are again values for the characteristic rock of the unit, in this case gabbro.

Table 2Overview of the parameters used for the parametrization.

a Turcotte and Schubert (2014). b Maystrenko and Scheck-Wenderoth (2013). c Tašárová et al. (2016). d Przybycin et al. (2015).

Download Print Version | Download XLSX

3.5 Initial stress state

Before applying displacement boundary conditions to the model, an initial stress state is generated representing a reference stress state. We use a simple semi-empirical function by Sheorey (1994) for the stress ratio k depending on depth (z) and Young's modulus (E) which can be considered as being representative of tectonically inactive regions with low lateral density contrasts:

(1) k = 0.25 + 7 E 0.001 + 1 z .

To achieve our initial stress state we compare the k values defined as

(2) k = S Hmean S V = S Hmax + S hmin 2 S V

from 29 synthetic profiles with the stress ratio calculated for a Young's modulus of 30 and 70 GPa, representing the sedimentary and upper crust units.

In order to establish the initial stress state, an underburden and a sideburden are added and this extended model is implemented in a conic shell (Fig. 6a and b). Then, the model has to settle down frictionless within this conical shell. During that procedure, Young's modulus in the underburden as well as Poison's ratio in all units is varied until the virtual wells fit the Sheorey equation (Eq. 1, Fig. 6c). This procedure has been used and described several times (Buchmann and Connolly, 2007; Hergert, 2009; Hergert and Heidbach, 2011; Reiter and Heidbach, 2014). The resulting stress state represents the initial stress state, which is subsequently perturbed by applying displacement boundary conditions that impose the tectonic stress.

Figure 6(a) Top view of the model implemented in the shell (green) and the sideburden (dark blue). Blue dots indicate the synthetic calibration wells. (b) Side view of the model implemented in the shell, the sideburden and the underburden (bright blue). (c) k values of the calibration wells (blue curves) in comparison with k values calculated with a semi-empirical function by Sheorey (1994) for a Young's modulus of 30 and 70 GPa representing the sedimentary and the upper crust units (red curves).


3.6 Displacement boundary conditions

The base of the model is fixed vertically, lateral movements are allowed and the model surface is free. At the five lateral boundaries of the model, displacement boundary conditions are applied to parametrize past and ongoing tectonic kinematics. The orientations of the model boundaries are chosen parallel or perpendicular to the mean SHmax orientation (Fig. 1b). The eastern and western lateral model boundaries are aligned parallel and the northern and southern boundaries perpendicular to the mean SHmax orientation. Accordingly, extension is applied to the eastern and western boundaries and shortening to the northern and southern ones (Fig. 3).

We use a two-stage approach to find a good agreement with the stress orientation and stress magnitude datasets. First a best fit with respect to a mean SHmax orientation (see details in Sect. 4.1) is estimated by an appropriate ratio between the extension and shortening applied. In a second step we vary the magnitude of these displacements on the model boundaries while keeping the ratio constant so that a best fit with the stress magnitude data is achieved as well. The calibration is mainly based on the Shmin magnitude due to the larger amount of data from the compilation of Morawietz and Reiter (2020) and the fact that SHmax magnitudes are often calculated and not measured and therefore less reliable. For the best-fit model a total extension of 465 m in east–west direction and a total shortening of 325 m in north–south direction are applied.

4 Results

4.1 Orientation of the maximum horizontal stress (SHmax)

We compare our model results with the stress orientation from the WSM database (Heidbach et al., 2016) and some additional data by Levi et al. (2019) from western Austria (Fig. 1a). From the WSM database we use only SHmax orientations that have a WSM quality A to C. However, we do not use individual data records but a mean SHmax orientation on a regular 0.5 grid (Figs. 1b and 7c). Using a mean SHmax orientation avoids effects of data clustering which is often the case in the WSM database and it filters the data for a wavelength of the stress pattern that is representative of the resolution of the model. For the estimation of the mean SHmax values, we use the tool stress2grid from Ziegler and Heidbach (2017). The SHmax data records are weighted according to their quality and their distance to the grid points. Each grid point requires at least 10 data points within a fixed search radius of 200 km. The resulting mean orientation of SHmax has a median standard deviation of  28 using the statistics of bi-polar data (Mardia, 1972). The model results are interpolated linearly on a plane at 5 km depth, and then the nearest value to each grid point is chosen for the comparison with the mean WSM data.

Figure 7a displays the SHmax orientation of the model at 5 km depth, whereas Fig. 7c shows the calculated mean SHmax orientation of the WSM data within the model area. The modeled SHmax orientations at the model boundaries are controlled by the assigned boundary conditions; thus the orientations are perpendicular to the northern and southern boundaries and parallel to the eastern and western edges. Within the model area the orientation of SHmax shows a homogenous pattern with a dominant NNW–SSE orientation which rotates slightly to a north–south orientation at the eastern boundary. Figure 7b visualizes the deviation of the model results from the mean WSM data. Blue indicates regions where the model results are rotated anti-clockwise with respect to the mean WSM data and orange regions with a clockwise rotation. There are three areas with larger deviations: one with primarily clockwise rotation in the area of Belgium and the two other areas, located in the northern and south-western part of the model, including the North German Basin (NGB), the eastern part of the Alps and western part of the Carpathians showing an anti-clockwise rotation. Apart from these two areas, the dominant color is orange, conterminous with a slight clockwise rotation. This trend is also visible in Fig. 7d, where the histogram of the deviation between the mean SHmax orientation derived from the WSM data and the modeled orientation is shown with a median deviation of 5.6 and a mean deviation, calculated from the absolute differences, of 15.6.

Figure 7Comparison of SHmax orientation of the model results with the mean SHmax orientation derived from WSM data (Heidbach et al., 2016). (a) SHmax orientation of the model at 5 km depth. (b) Deviation of the model result relative to the mean SHmax orientation derived from WSM data. (c) Orientation of the mean SHmax of WSM data (details are described in the text). (d) Histogram of the deviation of the modeled SHmax orientation to the mean SHmax orientation derived from WSM data. Coastlines and borders used in the figures are based on the Global Self-consistent Hierarchical High-resolution Geography (GSHHG) of Wessel and Smith (1996).

4.2 Stress magnitudes

4.2.1 Minimum horizontal stress (Shmin) magnitudes

The modeled magnitudes of Shmin in comparison to stress magnitude data of Morawietz and Reiter (2020) are shown in Fig. 8. The figure is divided into three subfigures displaying the differences depending on depth and quality (Fig. 8a), the spatial distribution of the calibration data (Fig. 8b) and a histogram showing the distribution of the differences between the modeled and observed Shmin magnitudes (Fig. 8c). The differences are calculated as interpolated model results minus data; thus positive differences correspond to model values that are too large and negative ones to model values that are too low. We use only data from Morawietz and Reiter (2020) with a quality of A, B and C and from depths >200 m, to avoid topographic effects. Thus, we use 74 Shmin magnitude data records from a depth of 200 to 4600 m, most of them from the upper 1000 m. As shown in Fig. 8b the data are mainly located within the south-western part with the exception of one measurement from the NGB. With 42 data records, more than half of all data records originate from three localities: from Falkenberg near the German–Czech border (Baumgärtner et al., 1987), from Benken in Switzerland (Nagra, 2001) and from Wittelsheim in eastern France (Cornet and Burlet, 1992). Due to the calibration process described in Sect. 3.6 a median difference of 0 MPa is achieved. The differences are, with two exceptions, in a range of 10 to 10 MPa and seem to be independent of depth and quality. This together with a mean difference of 3.3 MPa indicates a very good fit with the data of Morawietz and Reiter (2020).

Figure 8Shmin magnitudes of the model in comparison to the data of Morawietz and Reiter (2020). The differences are calculated as model results minus calibration data. (a) Depth depending differences. Color of dots indicates the quality of the calibration data. (b) Spatial distribution of the calibration data used, numbers indicating localities with multiple data used. (c) Histogram of the differences displayed in (a). Coastlines and borders used in the figures are based on the Global Self-consistent Hierarchical High-resolution Geography (GSHHG) of Wessel and Smith (1996).

4.2.2 Maximum horizontal stress (SHmax) magnitudes

For the model calibration regarding the SHmax magnitudes, 57 data records are used from the database of Morawietz and Reiter (2020). Again, only data with a quality of A to C and from a depth of >200 m are used. Similar to the Shmin data they are mainly located in the south-western part of the model area (Fig. 9b). The data are from seven different localities, whereby the data from Falkenberg near the German–Czech border (Baumgärtner et al., 1987) and from Benken in Switzerland (Nagra, 2001) with 25 data records make up almost half of the comparison data used. The mean of the absolute difference is 20.6 MPa, and the median difference is 19.3 MPa. This difference can be explained by the asymmetric depth distribution of the values (Fig. 9a). There are significantly more data records from shallower depths (200 to 1000 m), which indicate model results that are too large, than from the greater depths (>1000 m), which indicate results that are too low. Regardless of this, a trend is visible from positive to negative stress differences with increasing depth; i.e., the model seems to predict values of SHmax that are too large in the upper part of the model and values of SHmax that are too low in the deeper part. Furthermore, it is striking that the differences of quality A data are almost all negative and almost all of quality B and C are positive.

Figure 9SHmax magnitudes of the model in comparison to the data of Morawietz and Reiter (2020). The differences are calculated as model results minus calibration data. (a) Depth depending differences. Color of dots indicates the quality of the calibration data. (b) Spatial distribution of the calibration data used, numbers indicating localities with multiple data used. (c) Histogram of the differences displayed in (a). Coastlines and borders used in the figures are based on the Global Self-consistent Hierarchical High-resolution Geography (GSHHG) of Wessel and Smith (1996).

4.2.3 Stress gradients and stress regime

In addition to the calibration and comparison of the model with stress magnitude data of Morawietz and Reiter (2020), the absolute stress magnitudes of SV, Shmin and SHmax for three hypothetical wells up to 10 km depth are shown in Fig. 10. We have chosen these three locations (Fig. 11) due to the availability of stress data for a comparison later on, the quite uniform distribution over Germany (north, south-west and south-east) and the different depths of the crystalline surface. The hypothetic well 1 is entirely within the crystalline basement, well 2 entirely within the sedimentary unit and well 3 partly within the sediment unit and partly within the crystalline basement. As with the previous results we do not show the results of the upper 200 m. The depths are relative to the model surface and do not correspond to the z values of the model.

Figure 10Gradients of regime stress ratio (RSR), SV, Shmin and SHmax for three hypothetic wells in comparison with data. The orientation of SHmax (not shown here) is constant over the entire depth (well 1: 161, well 2: 163, well 3: 162). The colored lines show the RSR (red), the stresses of SV (yellow), Shmin (blue) and SHmax (green). Blue stars and green rectangles show measured respectively calculated magnitude of Shmin and SHmax. The uncertainties of the magnitudes if specified are displayed as error bars. Yellow dotted lines are calculated SV from density data. Well 1: comparison data from the “Kontinentale Tiefbohrung” (KTB) of Brudy et al. (1997). Well 2: comparison data from the NGB of Röckel and Lempp (2003). Well 3: comparison data from Soultz-sous-Forêts and Rittershoffen. Measured Shmin magnitudes of Valley and Evans (2007), calculated SHmax values of Klee and Rummel (1993), and calculated SV magnitudes based on density values of Azzola et al. (2019).


The SV gradients of well 1 and 2 are constant with the exception of a slight increase in the upper 1000 m, which is related to free surface effects. Well 3 also shows this effect but with an additional gradient change at 1500 to 3500 m depth. Above 1500 m depth the gradient corresponds to well 2 and below 3500 m to well 1. Overall well 2 shows the lowest SV gradient of about 22.5 MPa km−1, well 1 the highest SV gradient of about 27 MPa km−1 and well 3 is with  25.5 MPa km−1 in between. The horizontal stresses of Shmin and SHmax have almost constant gradients in well 1 and 2; only the absolute stresses differ. In well 1 we have a gradient of about 17 MPa km−1 resulting in 170 MPa at 10 km depth for Shmin and about 205 MPa at 10 km for SHmax. Due to the identical gradients the differential stress between SHmax and Shmin is constant 35 MPa all over the well path. Well 2 shows a similar pattern with a Shmin and SHmax gradient of about 15 MPa km−1 and a constant differential stress between Shmin and SHmax of about 15 MPa. The gradients of well 3 are, as with the SV magnitudes, a combination of well 1 and 2. This can be seen particularly clearly by the Shmin values of well 3. From the surface to a depth of  1500 m the gradient is quite similar to well 2 and below  3500 m to well 1. In between, the gradient increases to  25 MPa km−1, which is the highest gradient of all horizontal stresses displayed. As a result, the differential stress in well 3 between SHmax and Shmin also changes with depth. It amounts to 20 MPa at 1500 m depth, increasing with depth to about 40 MPa at 4000 m and then remains constant leading to 165 MPa for Shmin and 205 MPa for SHmax at 10 km depth. Well 3 thus shows the only significant change of horizontal differential stress with depth of all three wells shown and also the highest differential stress with a maximum of about 40 MPa below 4000 m depth.

All three wells show a change of the stress regime from strike-slip to normal faulting, with SV becoming greater than SHmax. In well 1 the transition is at about 3500 m, in well 2 at about 2500 m and in well 3 at about 4000 m depth. But despite these minor differences in depth, there are almost no differences between the stress regimes.

As an additional result the regime stress ratios (RSRs) (Simpson, 1997) for four model sections and for the three wells are shown in Figs. 10 and 11. The RSR (Eq. 3) is a unitless value combining the stress regime index n (Eq. 4) of Anderson (1905) and the ratio of stress differences ϕ (Eq. 5):

(3)RSR=n+0.5+-1n(ϕ-0.5)(Simpson, 1997)(4)n=0Shmin<SHmax<SV1Shmin<SV<SHmax2SV<Shmin<SHmax(Anderson, 1905)(5)ϕ=σ2-σ3σ1-σ3(Angelier, 1979).

The resulting value between 0 and 3 is the RSR indicating the stress regime divided into six classes from radial extension over extension, transtension, transpression and compression to constriction. The results in Fig. 11 show with the exception of peripheral areas, displaying some boundary effects, a rather uniform change of the stress regime from strike-slip to normal faulting with increasing depth. Starting with a RSR of 1 to 2 at 1000 m resulting in a RSR of 1 to 0.25 at 10 km depth. The two sections in between at 2000 and 4000 m show the transition with a dominant RSR of 1 to 1.5 and 0.75 to 1.25, respectively. The RSRs of the three wells displayed in Fig. 10 confirm this observation. In the shallower parts the RSR lies between 1.75 and 2 then decreasing with depth to values smaller than 0.5. A special aspect is visible in well 3, where the RSR is almost constant over 2 km along the transition between the sedimentary and the upper crust unit. The lowest RSR occurs in well 2, which is located entirely within the sedimentary unit. This correlation is also visible in Fig. 11, where the lower RSR is related to areas with high sediment thicknesses, e.g., the NGB. On the other hand, a higher RSR seems to correlate with basement areas like the Bohemian Massif or the Mid German Crystalline High, well visible at 2000 and 4000 m depth.

Figure 11RSR indicating the stress regime for four model sections for different depths. The three black dots show the locations of the three hypothetical wells in Fig. 10. The high RSR values at the model edges in the upper 4000 m, representing a constriction, are edge effects due to the applied boundary conditions. We used a color map based on “lajolla” of Crameri (2021).Coastlines and borders used in the figures are based on the Global Self-consistent Hierarchical High-resolution Geography (GSHHG) of Wessel and Smith (1996).

5 Discussion

5.1 Orientation of the maximum horizontal stress (SHmax)

The results of the SHmax orientation are in comparison to the mean WSM data quite good with a median deviation of 5.6 and an absolute mean deviation of 15.6 (Fig. 7). Therefore, the results are within the error range of the used WSM A to C quality data records which have uncertainties of 15 to 25 (Heidbach et al., 2018). Apart from some reorientation at the model edges in the upper 100 s of meters, the orientations are almost constant over the entire model depth. For example, the SHmax orientations of our three hypothetical wells are 161, 163 and 162 (Fig. 10).

However, there are three areas with noticeable lateral deviations: in the north-eastern part of Germany, in Belgium and along the Carpathians. The region in north-eastern Germany belongs to the NGB in which there are thick salt deposits. Salt can act as a mechanical decoupling horizon between the layers above and below (Ahlers et al., 2019; Bell, 1996; Cornet and Röckel, 2012; Heidbach et al., 2007; Hillis and Nelson, 2005; Tingay et al., 2011). In such cases, the stress state below represents the regional trend transferred through the crust, while the stress state above is only affected by local sources often controlled by local density and strength contrasts. More than 20 % of the data from this region are above the salt and mostly E–W oriented, in contrast to the data below, which are more N–S oriented (e.g., Cornet and Röckel, 2012; Grote, 1998; Röckel and Lempp, 2003; Roth and Fleckenstein, 2001). However, since we do not distinguish between the data from these different layers the derived mean SHmax values are influenced by these data above the salt layer. Possibly the misfit in this area can also be explained by the Pritzwalk anomaly, a positive gravity anomaly due to high-density lower crust (Krawczyk et al., 2008).

The deviations in Belgium and adjacent areas can have several reasons. This region is the border between two massifs, the Rhenish and the Brabant Massif (Pharaoh, 2018), and the strength contrast between these two massifs may play a role. Such contrasts are often considered to be responsible for reorientations in the stress field (e.g., Adams and Bell, 1991; Heidbach et al., 2007; Rajabi et al., 2017). Another reason could be the tectonically active Lower Rhine Basin nearby or the uplift of the Rhenish Massif (Reicherter et al., 2008). Possible boundary effects can be excluded, since the orientation of SHmax is uniform along the entire western boundary of the model. The deviations in the south-eastern part of the model are located along the Carpathians and the adjacent Pannonian Basin. This is possibly an area with low far-field or first-order horizontal stress sources resulting in a near-isotropic stress state (Heidbach et al., 2007), and thus the topography contrast between the mountain range and the Pannonian Basin probably has a dominant influence (Bada et al., 2001, 1998). Furthermore, the NE–SW-oriented SHmax indicated by the WSM data implies an NW–SE extension in this area, which is in agreement with the orientation of back-arc extension arising from the retreating slab beneath the Carpathians in Romania (Sperner et al., 2001), which the model does not account for.

In general, our modeled SHmax orientations show a rather simple stress pattern without local perturbations. This is to a certain extent an expectable result since our model is in equilibrium with gravitation. Therefore neither isostatic effects as described by, e.g., Kaiser et al. (2005), Bada et al. (2001) or Jarosiński et al. (2006) nor local perturbations due to faults or fault zones (Kaiser et al., 2005; Jarosiński et al., 2006) can be considered since such structural elements are not implemented. Nevertheless, our model results also show no impact of mechanical contrasts on the orientation of SHmax, e.g predicted by Grünthal and Stromeyer (1994), Marotta et al. (2002) or Cacace (2008) despite mechanical contrast, e.g., a Young's modulus difference of 40 GPa between the sedimentary (30 GPa) and the upper crust units (70 GPa). Probably a lateral stiffness difference and a weak unit seems to be necessary to get some perturbation due to a stiffness contrast (Reiter, 2021). To test this thesis, we defined an unrealistic low Young's modulus of 30 GPa to the upper crust unit of Avalonia (Fig. 5). In this case we could see perturbations at the border between the upper crust units of Avalonia and the EEC.

Although our kinematic boundary conditions applied are not derived from plate tectonic forces, they fit in general the tectonic setting of the model area. A shortening in N–S direction can be related to the alpine orogeny in the south, and an extension in E–W direction correlates with the evolution of several extensional structures like the Cenozoic Rift System or the Eger Graben since the Paleogene (Kley et al., 2008).

5.2 Magnitudes of Shmin, SHmax and SV

The Shmin magnitudes (Fig. 8a) in general show a very good correlation with the data of Morawietz and Reiter (2020) with a mean difference of 3.3 MPa, a median difference of 0 MPa and an almost even distribution independent of data quality and depth. However, the model was calibrated with these values, so the almost perfect match should not be overrated.

The comparison of the SHmax magnitudes does not show such a good match with a mean difference of 20.6 MPa and a median difference of 19.3 MPa (Fig. 9). Due to the calibration process described a much better fit should be achieved. But we have decided not to force a median of 0 MPa for various reasons. Compared to the Shmin magnitude data, the scattering is significantly larger and the distribution is less even. This is an expectable result since SHmax values are usually calculated and not measured and therefore SHmax magnitudes have a lower reliability compared to Shmin values (Morawietz et al., 2020), but additionally there seems to be a dependency on depth and quality. A major part of the data indicating SHmax magnitudes that are too large are from shallow depths (200 to 1000 m) (Fig. 9a). It can be assumed that the median and mean difference would be significantly better for a uniform depth distribution of the SHmax values since the results below 1000 m show a good match. A reduction of the SHmax magnitudes in our model and thus a statistically better fit would therefore only lead to a better fit of the model result in the upper part of the model. Whether there is a dependency of the results on quality is difficult to assess, although data of quality B and C tend to show larger deviations than data of quality A. But most quality B and C data are also from shallow depths. Therefore, the depth dependency may overlay the quality dependency.

The stress magnitudes of the three hypothetical wells displayed in Fig. 10 show by and large the expected results. Since SV is only dependent on the density the gradients of well 1 and 2 located in a single unit are constant all over the total well depth. The gradient of well 3 changes in between 1500 and 3500 m depth due to the change of units. The transition zone is quite large because of the vertical element resolution of about 800 m. Based on the sum of the overburden, the maximum SV at 10 km depth is the highest in well 1, followed by well 3 and 2. The stress differences between SHmax and Shmin are dependent on the elastic rock properties. Therefore, these results again are mainly based on the unit the well is located in. Since the Poisson's ratio is constant for the units involved (0.25), Young's modulus is probably the decisive parameter. This explains the constant horizontal stress differences within well 1 and 2 and the variations over the length of well 3. In addition, the maximum stress differences seemed to be mainly dependent on Young's modulus. Well 2 shows the smallest stress differences of 15 MPa, and well 1 shows differences of 35 MPa according to Young's modulus of 30 and 70 GPa, respectively. An exception is well 2 with the highest total stress differences of up to 40 MPa.

Also the RSRs displayed in Figs. 10 and 11 indicate a strong dependency on Young's modulus since the highest values occur usually in the units of the upper crust characterized by a Young's modulus of 70 GPa and the lowest values are visible in areas with a high sediment thicknesses with a Young's modulus of 30 GPa. This correlation is almost perfectly visible in 4000 m depth in Fig. 11 in comparison to the depth of the top of the crystalline basement shown in Fig. 4. This stiffness difference between these two units is also responsible for the constant RSR in well 3 in the transition zone between these units in 1500 to 3500 m depth (Fig. 10). The explanation for this correlation between the RSR and Young's modulus is larger horizontal stresses due to a higher Young's modulus, which lead to a more compressive regime and vice versa. In addition, our results indicate a change of the stress regime with depth for the whole model area from a dominant strike-slip regime (1 < RSR < 2) to a normal faulting regime (RSR < 1) (Figs. 10 and 11). This change occurs, with few exceptions, between 2000 and 4000 m depth. The stress regime and thus in particular a change with depth is a decisive factor, e.g., for the wellbore stability, especially in case of directional or deviated drilling (Rajabi et al., 2016) or the stimulation of enhanced geothermal reservoirs (Azzola et al., 2019). Such depth-dependent stress regimes are for example described by Brooke-Barnett et al. (2015), Cornet et al. (2007), Rajabi et al. (2016) and Rajabi et al. (2017).

To get a more detailed insight we compare our hypothetic wells 1, 2 and 3 with local magnitude data (Fig. 10). The model results of the hypothetic well 1 are displayed in comparison to values of the “Kontinentale Tiefbohrung” (KTB), a major scientific drilling project in Germany (Brudy et al., 1997). Our results of SV are in a very good agreement with the SV calculated from a mean density value. Only at greater depths does the difference increase to about 5 MPa, and in the uppermost 750 m our results are too large, possibly due to free surface effects caused by our model resolution. The results of Shmin and SHmax show significantly larger differences to the data of Brudy et al. (1997). Except for the values at 3000 m depth the Shmin magnitudes of Brudy et al. (1997) are at least 15 MPa larger than the model results. The maximum difference of about 35 MPa is at 6400 m depth. The results of SHmax show even greater deviations. All values of Brudy et al. (1997) are at least 15 MPa larger than the model results. The maximum difference is about 180 MPa at 7800 m depth and thus larger than our model results with 160 MPa. Remarkable is the change in the horizontal stress magnitudes of Brudy et al. (1997) at 3000 m depth. The Shmin magnitudes increase from 50 to 70 MPa within 200 m, and SHmax increases even by 30 MPa from 100 to 130 MPa. At the same time the inaccuracies also increase significantly. This can be attributed to the fact that the values between 3000 and 7000 m depth are only calculated and not directly measured (Brudy et al., 1997), which is why the values only got a quality of worse than C (Morawietz et al., 2020) and are not used by us for calibration. A remarkable difference between our model and the geomechanical properties at the KTB site is the values of Young's modulus. The calculated values of Brudy et al. (1997) are about 90 GPa on average between 3000 and 8000 m depth, which is about 20 MPa larger than our model assumption of 70 GPa in this area. Furthermore, in our model a normal faulting regime is established from about 3500 m depth downwards (Fig. 11), which is contrary to the stress regime of Brudy et al. (1997) showing a strike-slip regime from 1 km depth downwards. This indicates that our SHmax values are possibly too low within this model area. In general, our model results show a constant differential stress between Shmin and SHmax of 35 MPa, whereas the data of the KTB indicate an increasing differential stress with depth. The fact that our model does not include faults can also have an effect. The KTB is located above the Franconian Line and even intersecting it (Wagner et al., 1997). The Franconian Line is a major fault zone at the south-western margin of the Bohemian Massif with a polyphase development from late Paleozoic to Neogene times (Zulauf, 1993; Peterek et al., 1997). In the end it is probably a combination of the lower Young's modulus in the model, the very large uncertainties of the calculated values and SHmax values that are too low in our model, which may explain the differences of up to 180 MPa for the SHmax magnitudes.

The stress magnitudes of the hypothetic well 2 are shown in comparison with SV and Shmin data from the NGB of Röckel and Lempp (2003). The SV values are in good agreement with our results down to  2000 m depth, whereas the difference increases below this level. This shows that the density chosen for the sedimentary unit is at least appropriate for the upper part of this unit. A larger density would lead to better results at greater depths, but since our calibration data mainly comes from depths shallower than 3500 m (Figs. 8 and 9), we consider assuming a density of 2300 kg m−3 for the sedimentary unit is reasonable. The Shmin values indicate a good fit with our results across the entire depth range to 7000 m, which agrees with our general comparison shown in Fig. 8. Due to missing data, a comparison is not possible for the SHmax values. However, Röckel and Lempp (2003) mention that the actual stress regime in the NGB can be characterized as normal faulting for the sub-salt level. At an average depth of the salt layer in the NGB of 4 km (Scheck-Wenderoth and Lamarche, 2005), our results show a normal faulting regime beneath 4 km depth too (Figs. 10 and 11).

For the comparison of our hypothetic well 3, we use data from the geothermal project in Soultz-sous-Forêts and some values from Rittershoffen, another geothermal project nearby, both located at the western edge of the URG. The URG is part of the major European Cenozoic Rift System located in south-western Germany and eastern France (Ziegler and Dèzes, 2006). The dashed SV gradients in Fig. 10 are calculated on the basis of density values of Rittershoffen (Azzola et al., 2019) correlated to the stratigraphic column of Soultz-sous-Forêts based on Aichholzer et al. (2016) up to 5080 m, the total depth of the deepest well in Soultz-sous-Forêts. Despite the density change between the sedimentary unit and the crystalline basement, the SV magnitudes of the model are in good agreement. An exception is the upper 750 m, where modeled SV magnitudes are slightly too low. The data of Valley and Evans (2007) show measured Shmin magnitudes between 1500 and 4500 m depth. They are in very good agreement with our model results within the upper 3000 m. In between 3500 and 4500 m depth the agreement is slightly worse but with a maximum difference of about 10 MPa at 4500 m depth still good. For the validation of the SHmax magnitudes, we use four calculated values of Klee and Rummel (1993) between 2200 and 3500 m. Within this depth interval our model results show a good correlation with a maximum deviation of about 5 MPa. Due to the small amount of data available for SHmax, a comparison with some calculated stress gradients can be helpful. Assumptions for stress gradients of Heinemann (1994), Klee and Rummel (1993), and Valley and Evans (2007) result in SHmax magnitudes of 35 up to 55 MPa at a depth of 2000 m, of 95 up to 143 MPa at 5000 m and of 195 up to 310 MPa at 10 km. Even though these values show a quite wide range, the comparison of the model results allows the conclusion that our SHmax values show a quite good agreement in the upper 5000 m but tend to be too low with increasing depth. This is also supported by the observation of seismic events, which show a slight dominance of strike-slip versus normal faulting focal mechanisms at depths of 8 to 10 km in the URG (Cornet et al., 2007). In contrast, our results show a normal faulting regime, which implies SHmax values that are too low.

In general, the model results show a good agreement with real magnitude data. The SV magnitudes show a good correlation for all three described cases of the KTB, the NGB and the Soultz-sous-Forêts site (Fig. 10). The modeled SV magnitudes appear to be only slightly too low with increasing depth, but in general the densities seem to be quite well chosen. A better agreement would probably only be possible with a higher stratigraphic resolution in the sedimentary unit and a density gradient within the upper crust. Such a simple gradient, which could also reduce the differences in the sedimentary unit, would be less useful in this case, because densities between two sedimentary units can differ considerably, independent of their depth. Overall, the Shmin magnitudes show a very good correlation. This can be seen in the general location-independent comparison in Fig. 8 as well as with regard to the Soultz-sous-Forêts and NGB location in Fig. 10. Only the results for the KTB site show some considerable differences, with the greatest deviations from calculated Shmin values (3000 to 7000 m depth) that have not been directly measured and are of rather questionable quality. The SHmax magnitudes show the largest deviations, both in the general comparison (Fig. 9) and in the local comparisons (Fig. 10). The general comparison shows that our values in the upper part are rather too large and at greater depths rather too small. The results of SHmax at the KTB and Soultz-sous-Forêts sites only confirm the trend of values that are too low in the deeper parts but not the trend of values that are too low in the shallower parts. In contrast, the SHmax magnitudes at Soultz-sous-Forêts show a good agreement down to 2000 m depth and the values at the KTB site even indicate magnitudes down to 800 m depth that are too low. An indication of SHmax magnitudes that are generally too small with increasing depth is also the RSR values at 10 km depth (Fig. 10), which show larger areas of values lower than 0.5 indicating a radial extension, an uncommon stress regime in the upper crust. In the upper part of the model up to 4000 m, the rather uniform tectonic regime, between normal faulting and transtensional, corresponds mainly to the tectonic conditions expected (Röckel and Lempp, 2003; Cornet et al., 2007). However, in detail, e.g., for the KTB site, the model cannot reflect differing local conditions. This could simply be a consequence of the simplifications made, which cannot resolve all local conditions, e.g., differing rock properties or nearby faults.

6 Conclusions

The model presented is the first 3D geomechanical model for Germany predicting the first-order 3D stress tensor. The model is calibrated with SHmax orientations from the WSM database (Heidbach et al., 2016) and compilation of Shmin and SHmax stress magnitude data from Morawietz and Reiter (2020). Overall, our model shows good results regarding the orientation of SHmax and Shmin magnitudes despite the necessary simplifications due to the model resolution and rock property distributions as well as the highly irregular spread of the calibration data and their varying quality. The SHmax orientations of the model are to a large extent within the uncertainty of the mean SHmax orientations that are derived from the A to C quality data of the WSM database. Furthermore, the Shmin magnitudes show a quite good fit to various datasets (Röckel and Lempp, 2003; Valley and Evans, 2007), but the SHmax magnitude results show locally significant differences. Modeled SHmax magnitudes are too small in the lower part of the model, whereas some results indicate values that are too high in the upper part. But in general, our model describes the regional 3D contemporary stress state quite well. Some larger deviations due to local structures are expectable. Therefore, the model results cannot be used for stress prediction on a local or reservoir scale as the resolution is not sufficient, but it can deliver initial stress conditions for smaller-scale models that contain little or no stress magnitude data at all.

To improve our large-scale model a better stratigraphic resolution of the sedimentary unit and thus a better representation of the lithologies has to be implemented. This would increase the reliability of the comparison between measured stress magnitude data and the modeled ones. In addition to a vertical refinement, resolving lateral variations of the rock properties would be useful and these potentially account for lateral variability of the stress tensor.

Data availability

The results of our model have been published and are publicly available at (see Ahlers et al., 2021).


The supplement related to this article is available online at:

Author contributions

Conceptualization of the project was done by AH, TH, KR, OH and BM. Construction and discretization of the model were done by SA. Data for the model and its calibration were collected and provided by SA, TH, LR, SM, MSW and DA. Evaluation of the model results and their interpretation were performed by SA with the support of AH, TH, KR, BM, LR, OH and SM. SA wrote the paper with help from all coauthors. All authors read and approved the final paper.

Competing interests

The authors declare that they have no conflict of interest.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


This study is part of the SpannEnD Project (, last access: 2 August 2021), which is supported by Federal Ministry for Economic Affairs and Energy (BMWI) and managed by Projektträger Karlsruhe (PTKA) (project code: 02E11637A).

Financial support

This research has been supported by the Bundesministerium für Wirtschaft und Energie (grant no. 02E11637A).

Review statement

This paper was edited by Patrice Rey and reviewed by two anonymous referees.


Adams, J. and Bell, J. S.: Crustal stresses in Canada, in: Neotectonics of North America, edited by: Slemmons, D. B., Engdahl, E. R., Zoback, M. D., and Blackwell, D. D., Geol. Soc. Am. USA, 1, 367–386,, 1991. 

Ahlers, S., Hergert, T., and Henk, A.: Numerical Modelling of Salt-Related Stress Decoupling in Sedimentary Basins–Motivated by Observational Data from the North German Basin, Geosciences, 9, 19,, 2019. 

Ahlers, S., Henk, A., Hergert, T., Reiter, K., Müller, B., Röckel, L., Heidbach, O., Morawietz, S., Scheck-Wenderoth, M., and Anikiev, D.: Crustal stress state of Germany – Results of a 3D geomechnical model, TUdatalib [data set],, 2021. 

Aichholzer, C., Duringer, P., Orciani, S., and Genter, A.: New stratigraphic interpretation of the Soultz-sous-Forêts 30-year-old geothermal wells calibrated on the recent one from Rittershoffen (Upper Rhine Graben, France), Geothermal Energy, 4, 1–26,, 2016. 

Altmann, J. B., Müller, B., Müller, T. M., Heidbach, O., Tingay, M., and Weißhardt, A.: Pore pressure stress coupling in 3D and consequences for reservoir stress states and fault reactivation, Geothermics, 52, 195–205,, 2014. 

Amadei, B. and Stephansson, O.: Rock Stress and Its Measurement, Springer Netherlands, Dordrecht, 490 pp.,, 1997. 

Anderson, E. M.: The dynamics of faulting, Trans. Edin. Geol. Soc., 8, 387–402,, 1905. 

Andeweg, B.: Cenozoic tectonic evolution of the Iberian Peninsula: Effects and causes of changing stress fields, Ph. D. thesis, Faculty of Earth Sciences, Vrije Universiteit, Amsterdam, 178 pp., 2002. 

Angelier, J.: Determination of the mean principal directions of stresses for a given fault population, Tectonophysics, 56, T17–T26,, 1979. 

Anikiev, D., Lechel, A., Gomez Dacal, M. L., Bott, J., Cacace, M., and Scheck-Wenderoth, M.: A three-dimensional lithospheric-scale thermal model of Germany, Adv. Geosci., 49, 225–234,, 2019. 

Asch, K.: The 1:5 Million International Geological Map of Europe and Adjacent Areas (IGME5000), Bundesanstalt für Geowissenschaften und Rohstoffe, Hannover, 2005. 

Azzola, J., Valley, B., Schmittbuhl, J., and Genter, A.: Stress characterization and temporal evolution of borehole failure at the Rittershoffen geothermal project, Solid Earth, 10, 1155–1180,, 2019. 

Bada, G., Horváth, F., Cloetingh, S., Coblentz, D. D., and Tóth, T.: Role of topography-induced gravitational stresses in basin inversion: The case study of the Pannonian basin, Tectonics, 20, 343–363,, 2001. 

Bada, G., Cloetingh, S., Gerner, P., and Horváth, F.: Sources of recent tectonic stress in the Pannonian region: Inferences from finite element modelling, Geophys. J. Int., 134, 87–101,, 1998. 

Baumgärtner, J., Rummel, F., and Zhaotan, C.: Wireline hydraulic fracturing stress measurements in the Falkenberg granite massif, Geol. Jb., 39, 83–99, 1987. 

Behr, H. J., Duerbaum, H. J., Bankwitz, P., Bankwitz, E., Benek, R., Berger, H. J., Brause, H., Conrad, W., Foerste, K., Frischbutter, A., Gebrande, H., Giese, P., Goethe, W., Guertler, J., Haenig, D., Haupt, M., Heinrichs, T., Horst, W., Hurtig, E., and Kaempf, H.: Crustal structure of the Saxothuringian Zone; results of the deep seismic profile MVE-90(East), Z. Geol. Wissenschaft., 22, 647–770, 1994. 

Bell, J. S.: Petro geoscience 2, In situ stresses in sedimentary rocks (part 2): Applications of stress measurements, Geosci. Can., 23, 135–153, 1996. 

Bell, J. S.: Practical methods for estimating in situ stresses for borehole stability applications in sedimentary basins, J. Petrol. Sci. Eng., 38, 111–119,, 2003. 

BGR: Abriss der Standortauswahl und Darstellung der angewandten geowissenschaftlichen Kriterien bei den Endlagerprojekten in den Ländern Schweiz, Frankreich, Schweden, Belgien und USA, Hannover, 126 pp., 2015. 

Bokelmann, G. and Bianchi, I.: Imaging the Variscan suture at the KTB deep drilling site, Germany, Geophys. J. Int., 213, 2138–2146,, 2018. 

Bormann, P., Bankwitz, P., Apitz, E., Bankwitz, E., and Franzke, H. J.: Komplexinterpretation des Profilnetzes ZENTROSEIS – G4 Bericht der SAG Tiefenerkundung – Abschlussbericht, Zentralinstitut für Physik der Erde, Potsdam, 162 pp., 1986. 

Brady, B. H. and Brown, E. T.: Rock Mechanics for underground mining, 3th, Springer Netherlands, Dordrecht, 628 pp.,, 2004. 

Brooke-Barnett, S., Flottmann, T., Paul, P. K., Busetti, S., Hennings, P., Reid, R., and Rosenbaum, G.: Influence of basement structures on in situ stresses over the Surat Basin, southeast Queensland, J. Geophys. Res., 120, 4946–4965,, 2015. 

Brückl, E., Behm, M., Decker, K., Grad, M., Guterch, A., Keller, G. R., and Thybo, H.: Crustal structure and active tectonics in the Eastern Alps, Tectonics, 29, TC2011,, 2010. 

Brudy, M., Zoback, M. D., Fuchs, K., Rummel, F., and Baumgärtner, J.: Estimation of the complete stress tensor to 8 km depth in the KTB scientific drill holes: Implications for crustal strength, J. Geophys. Res., 102, 18453–18475,, 1997. 

Buchmann, T. J. and Connolly, P. T.: Contemporary kinematics of the Upper Rhine Graben: A 3D finite element approach, Glob. Planet. Change, 58, 287–309,, 2007. 

Cacace, M.: Stress and Strain modelling of the Central European Basin System, Ph. D. thesis, Freie Universität Berlin, Berlin, 167 pp.,, 2008. 

Cazes, M., Torreilles, G., Bois, C., Damotte, B., Galdeano, A., Hirn, A., Mascle, A., Matte, P., van Ngoc, P., and Raoult, J. F.: Structure de la croute hercynienne du Nord de la France; premiers resultats du profil ECORS, B. Soc. Geol. Fr., 8, 925–941,, 1985. 

Cornet, F. H. and Burlet, D.: Stress field determinations in France by hydraulic tests in boreholes, J. Geophys. Res.-Sol. Ea., 97, 11829–11849,, 1992. 

Cornet, F. H. and Röckel, T.: Vertical stress profiles and the significance of “stress decoupling”, Tectonophysics, 581, 193–205,, 2012. 

Cornet, F. H., Bérard, T., and Bourouis, S.: How close to failure is a granite rock mass at a 5 km depth?, Int. J. Rock Mech. Min, 44, 47–66,, 2007. 

Crameri, F.: Scientific colour maps, Zenodo,, 2021. 

Diebold, P., Naef, H., and Ammann, M.: NTB 90-04: Zur Tektonik der zentralen Nordschweiz - Interpretation aufgrund regionaler Seismik, Oberflächengeologie und Tiefbohrungen, Nagra, Wettingen, 277 pp., 1991. 

Diederichs, M., Kaiser, P., and Eberhardt, E.: Damage initiation and propagation in hard rock during tunnelling and the influence of near-face stress rotation, Int. J. Rock Mech. Min, 41, 785–812,, 2004. 

Fischer, K. and Henk, A.: A workflow for building and calibrating 3-D geomechanical models – a case study for a gas reservoir in the North German Basin, Solid Earth, 4, 347–355,, 2013. 

Franke, W.: Tectonostratigraphic units in the Variscan belt of central Europe, GSA Special Papers, 230, 67–90,, 1989. 

Franke, W.: The Variscan orogen in Central Europe: Construction and collapse, in: European Lithosphere Dynamics, edited by: Gee, D. R. and Stephenson, R., Geological Society of London, London, 333–343,, 2006. 

Freeman, R. and Mueller, S. (Eds.): A continent revealed: The European geotraverse, Cambridge Univ. Pr, Cambridge, 275 pp., 1992. 

GeORG-Projektteam: Geopotentiale des tieferen Untergrundes im Oberrheingraben: Fachlich-Technischer Abschlussbericht des INTERREG-Projekts GeORG, Teil 4, Freiburg i. Br., 104 pp., 2013. 

Geothermieatlas Bayern: Grundgebirge (Prä-Perm) (Verbreitung und Tiefenlage), Bayerisches Staatsministerium für Wirtschaft, Landesentwicklung und Energie, available at: (last access: 2 August 2021), 2004. 

Goelke, M. and Coblentz, D.: Origins of the European regional stress field, Tectonophysics, 266, 11–24,, 1996. 

Grad, M., Brückl, E., Majdański, M., Behm, M., and Guterch, A.: Crustal structure of the Eastern Alps and their foreland: seismic model beneath the CEL10/Alp04 profile and tectonic implications, Geophys. J. Int., 177, 279–295,, 2009a. 

Grad, M., Tiira, T., and ESC Working Group: The Moho depth map of the European Plate, Geophys. J. Int., 176, 279–292,, 2009b. 

Grote, R.: Die rezente horizontale Hauptspannungsrichtung im Rotliegenden und Oberkarbon in Norddeutschland, Erdöl-Erdgas-Kohle, 114, 478–483, 1998. 

Grünthal, G. and Stromeyer, D.: The recent crustal stress field in Central Europe sensu lato and its quantitative modelling, Geol. Mijnbouw, 73, 173–180, 1994. 

Heidbach, O., Reinecker, J., Tingay, M., Müller, B., Sperner, B., Fuchs, K., and Wenzel, F.: Plate boundary forces are not enough: Second- and third-order stress patterns highlighted in the World Stress Map database, Tectonics, 26, TC6014,, 2007. 

Heidbach, O., Tingay, M., Barth, A., Reinecker, J., Kurfeß, D., and Müller, B.: Global crustal stress pattern based on the World Stress Map database release 2008, Tectonophysics, 482, 3–15,, 2010. 

Heidbach, O., Hergert, T., Reiter, K., and Giger, S.: NAB 13-88: Local Stress field sensitivity analysis – Case study Nördlich Langen, Wettingen, 50 pp., 2014. 

Heidbach, O., Rajabi, M., Reiter, K., Ziegler, M., and WSM Team: World Stress Map Database Release 2016 v1.1, GFZ Data Services [data set],, 2016. 

Heidbach, O., Rajabi, M., Cui, X., Fuchs, K., Müller, B., Reinecker, J., Reiter, K., Tingay, M., Wenzel, F., Xie, F., Ziegler, M. O., Zoback, M.-L., and Zoback, M.: The World Stress Map database release 2016: Crustal stress pattern across scales, Tectonophysics, 744, 484–498,, 2018. 

Heinemann, B.: Results of scientific investigations at the HDR test site Soultz-sous-Forêts: Alsace (1987–1992), SOCOMINE report, 126 pp., 1994. 

Heinrichs, T., Giese, P., Bankwitz, P., and Bankwitz, E.: Dekorp 3/MVE-90(West) – preliminary geological interpretation of a deep near-vertical reflection profile between the Rhenish and the Bohemian Massifs, Germany, Z. Geol. Wissenschaft., 22, 771–801, 1994. 

Henk, A.: Perspectives of Geomechanical Reservoir Models – Why Stress is Important, Oil Gas: European Magazine, 125, OG20–OG24, 2009. 

Hergert, T.: Numerical modelling of the absolute stress state in the Marmara region – a contribution to seismic hazard assessment, Dissertation, Universität Karlsruhe,, 2009. 

Hergert, T. and Heidbach, O.: Geomechanical model of the Marmara Sea region-II, 3-D contemporary background stress field, Geophys. J. Int., 185, 1090–1102,, 2011. 

Hergert, T., Heidbach, O., Reiter, K., Giger, S. B., and Marschall, P.: Stress field sensitivity analysis in a sedimentary sequence of the Alpine foreland, Northern Switzerland, Solid Earth, 6, 533–552,, 2015. 

Hettema, M.: Analysis of mechanics of fault reactivation in depleting reservoirs, Int. J. Rock Mech. Min, 129, 104290,, 2020. 

Hillis, R. R. and Nelson, E. J.: In situ stresses in the North Sea and their applications: Petroleum geomechanics from exploration to development, in: Petroleum Geology: North-West Europe and Global Perspectives – Proceedings of the 6th Petroleum Geology Conference, 551–564,, 2005. 

Hirschmann, G.: KTB – The structure of a Variscan terrane boundary: seismic investigation – drilling – models, Tectonophysics, 264, 327–339,, 1996. 

Hurtig, E., Cermak, V., Haenel, R., and Zui, V.: Geothermal atlas of Europe, Haack, Gotha, Germany, 156 pp., 1992. 

Jaeger, J. C., Cook, N. G. W., and Zimmerman, R. W.: Fundamentals of rock mechanics, 4th Edn., Blackwell Publ, Malden, MA, 475 pp., 2011. 

Janik, T., Grad, M., Guterch, A., Vozár, J., Bielik, M., Vozárova, A., Hegedűs, E., Kovács, C. A., Kovács, I., and Keller, G. R.: Crustal structure of the Western Carpathians and Pannonian Basin: Seismic models from CELEBRATION 2000 data and geological implications, J. Geodyn., 52, 97–113,, 2011. 

Jarosiński, M., Beekman, F., Bada, G., and Cloetingh, S.: Redistribution of recent collision push and ridge push in Central Europe: insights from FEM modelling, Geophys. J. Int., 167, 860–880,, 2006. 

Kaiser, A., Reicherter, K., Huebscher, C., Gajewski, D., Marotta, A. M., and Bayer, U.: Variation of the present-day stress field within the North German Basin; insights from thin shell FE modeling based on residual GPS velocities, Tectonophysics, 397, 55–72,, 2005. 

Kirsch, M., Kroner, U., Hallas, P., and Stephan, T.: 3D Model of the Erzgebirge – Crustal-Scale 3D Modelling of the Allochthonous Domain of the Saxo-Thuringian Zone, available at: (last access: 30 April 2019), 2017. 

Klee, G. and Rummel, F.: Hydrofrac stress data for the European HDR research project test site Soultz-Sous-Forets, Int. J. Rock Mech. Min, 30, 973–976,, 1993. 

Kley, J. and Voigt, T.: Late Cretaceous intraplate thrusting in central Europe: Effect of Africa-Iberia-Europe convergence, not Alpine collision, Geology, 36, 839–842,, 2008. 

Kley, J., Franzke, H.-J., Jähne, F., Krawczyk, C., Lohr, T., Reicherter, K., Scheck-Wenderoth, M., Sippel, J., Tanner, D., and van Gent, H.: Strain and Stress, in: Dynamics of complex intracontinental basins: The Central European Basin System, edited by: Littke, R., Bayer, U., Gajewski, D., and Nelskamp, S., Springer, Berlin Heidelberg, 97–124,, 2008. 

Konstantinovskaya, E., Malo, M., and Castillo, D. A.: Present-day stress analysis of the St. Lawrence Lowlands sedimentary basin (Canada) and implications for caprock integrity during CO2 injection operations, Tectonophysics, 518-521, 119–137,, 2012. 

Korsch, R. J. and Schäfer, A.: The Permo-Carboniferous Saar-Nahe Basin, south-west Germany and north-east France: basin formation and deformation in a strike-slip regime, Geol. Rundsch., 84, 293–318,, 1995. 

Kossmat, F.: Gliederung des variszischen Gebirgsbaus, Abh. Sächs. Geol. Landesamtes, 1, 1–39, 1927. 

Krawczyk, C. M., Rabbel, W., Willert, S., Hese, F., Götze, H.-J., Gajewski, D., and SPP-Geophysics Group: Crustal structures and properties in the Central European Basin system from geophysical evidence, in: Dynamics of complex intracontinental basins: The Central European basin system, edited by: Littke, R., Bayer, U., Gajewski, D., and Nelskamp, S., Springer, Berlin, Heidelberg, 67–95,, 2008. 

Kristiansen, T. G.: Drilling Wellbore Stability in the Compacting and Subsiding Valhall Field, IADC/SPE Drilling Conference, 2-4 March, Dallas, Texas,, 2004. 

Kroner, U., Romer, R. L., and Linnemann, U.: The Saxo-Thuringian Zone of the Variscan Orogen as part of Pangea, in: Pre-Mesozoic geology of Saxo-Thuringia: From the Cadomian active margin to the Variscan orogen, edited by: Linnemann, U. and Romer, R. L., Schweizerbart, Stuttgart, 3–16, 2010. 

Levi, N., Habermueller, M., Exner, U., Piani, E., Wiesmayr, G., and Decker, K.: The stress field in the frontal part of the Eastern Alps (Austria) from borehole image log data, Tectonophysics, 769, 228175,, 2019. 

Lindner, H., Scheibe, K., Seidel, K., and Hoffmann, N.: Berechnung von Relief, Tiefenlage und Magnetisierung des magnetisch wirksamen Kristallins für das Norddeutsche Becken, Z. Angew. Geol., 50, 65–74, 2004. 

Linnemann, U., D'Lemos, R., Drost, K., Jeffries, T., Gerdes, A., Romer, R. L., Samson, S. D., and Strachan, R. A.: Cadomian tectonics, in: The Geology of Central Europe Volume 1: Precambrian and Palaeozoic; Volume 2: Mesozoic and Cenozoic, edited by: McCann, T., Geol. Soc. Lond., 103–154,, 2008. 

Ljunggren, C., Chang, Y., Janson, T., and Christiansson, R.: An overview of rock stress measurement methods, Int. J. Rock Mech. Min, 40, 975–989,, 2003. 

Mardia, K. V.: Statistics of Directional Data: Probability and Mathematical Statistics, Academic Press, London, 380 pp., 1972. 

Marotta, A. M., Bayer, U., Thybo, H., and Scheck, M.: Origin of the regional stress in the North German Basin – results from numerical modelling, Tectonophysics, 360, 245–264,, 2002. 

Maystrenko, Y. P. and Scheck-Wenderoth, M.: 3D lithosphere-scale density model of the Central European Basin System and adjacent areas, Tectonophysics, 601, 53–77,, 2013. 

Mazur, S., Mikolajczak, M., Krzywiec, P., Malinowski, M., Buffenmyer, V., and Lewandowski, M.: Is the Teisseyre-Tornquist Zone an ancient plate boundary of Baltica?, Tectonics, 34, 2465–2477,, 2015. 

McCann, T. (Ed.): The Geology of Central Europe Volume 1: Precambrian and Palaeozoic; Volume 2: Mesozoic and Cenozoic, Geol. Soc. Lond., 1449 pp.,, 2008. 

Meissner, R. and Bortfeld, R. K.: DEKORP-Atlas: Results of Deutsches Kontinentales Reflexionsseismisches Programm, Springer Berlin Heidelberg, Berlin, Heidelberg, 21 pp.,, 1990. 

Meschede, M. and Warr, L. N.: The Geology of Germany, Springer, 304 pp.,, 2019. 

Morawietz, S. and Reiter, K.: Stress Magnitude Database Germany v1.0, GFZ Data Services [data set],, 2020. 

Morawietz, S., Heidbach, O., Reiter, K., Ziegler, M., Rajabi, M., Zimmermann, G., Müller, B., and Tingay, M.: An open-access stress magnitude database for Germany and adjacent regions, Geothermal Energy, 8,, 2020. 

Nagra: Sondierbohrung Benken: Technical Report NTB 00-01, Nagra, 288 pp., 2001. 

Nagra: Vorschlag geologischer Standortgebiete für das SMA- und das HAA-Lager, Begründung der Abfallzuteilung, der Barrieresysteme und der Anforderungen an die Geologie, Bericht zur Sicherheit und technischen Machbarkeit: NTB 08-05, Nagra, Wettingen, 2008. 

Oncken, O.: Transformation of a magmatic arc and an orogenic root during oblique collision and it's consequences for the evolution of the European Variscides (Mid-German Crystalline Rise), Geol. Rundschau, 86, 2–20,, 1997. 

Oncken, O., Plesch, A., Weber, J., Ricken, W., and Schrader, S.: Passive margin detachment during arc-continent collision (Central European Variscides), in: Orogenic Processes: Quantification and Modelling in the Variscan Belt, edited by: Franke, W., Haak, V., Oncken, O., and Tanner, D., London, 199–216,, 2000. 

Peterek, A., Rauche, H., Schröder, B., Franzke, H.-J., Bankwitz, P., and Bankwitz, E.: The late-and post-Variscan tectonic evolution of the Western Border fault zone of the Bohemian massif (WBZ), Geol. Rundsch., 86, 191–202,, 1997. 

Pharaoh, T.: The Anglo-Brabant Massif: Persistent but enigmatic palaeo-relief at the heart of western Europe, P. Geol. Assoc., 129, 278–328,, 2018. 

Przybycin, A. M., Scheck-Wenderoth, M., and Schneider, M.: Assessment of the isostatic state and the load distribution of the European Molasse Basin by means of lithospheric scale 3D structural and 3D gravity modelling, Int. J. Earth Sci., 104, 1405–1424,, 2015. 

Rajabi, M., Tingay, M., and Heidbach, O.: The present-day state of tectonic stress in the Darling Basin, Australia: Implications for exploration and production, Mar. Petrol. Geol., 77, 776–790,, 2016. 

Rajabi, M., Tingay, M., Heidbach, O., Hillis, R., and Reynolds, S.: The present-day stress field of Australia, Earth-Sci. Rev., 168, 165–189,, 2017. 

Reicherter, K., Froitzheim, N., Jarosinski, M., Badura, J., Franzke, H.-J., Hansen, M., Hubscher, C., Müller, R., Poprawa, P., Reinecker, J., Stackebrandt, W., Voigt, T., von Eynatten, H., and Zuchiewicz, W.: Alpine tectonics north of the Alps, in: The Geology of Central Europe Volume 1: Precambrian and Palaeozoic, Vol. 2, Mesozoic and Cenozoic, edited by: McCann, T., Geol. Soc. Lond., 1233–1285,, 2008. 

Reinhold, K.: Tiefenlage der ”Kristallin-Oberfläche” in Deutschland – Abschlussbericht, Bundesanstalt für Geowissenschaften und Rohstoffe, Hannover, 89 pp., 2005. 

Reiter, K.: Stress rotation – impact and interaction of rock stiffness and faults, Solid Earth, 12, 1287–1307,, 2021. 

Reiter, K. and Heidbach, O.: 3-D geomechanical–numerical model of the contemporary crustal stress state in the Alberta Basin (Canada), Solid Earth, 5, 1123–1149,, 2014. 

Reiter, K., Heidbach, O., Reinecker, J., Müller, B., and Röckel, T.: Spannungskarte Deutschland 2015, Erdöl-Erdgas-Kohle, 131, 437–442, 2015. 

Reiter, K., Heidbach, O., Müller, B., Reinecker, J., and Röckel, T.: Stress Map Germany 2016,, 2016. 

Röckel, T. and Lempp, C.: Der Spannungszustand im Norddeutschen Becken, Erdöl-Erdgas-Kohle, 119, 73–80, 2003. 

Roth, F. and Fleckenstein, P.: Stress orientations found in NE Germany differ from the West European trend, Terra Nova, 13, 289–296,, 2001. 

Rupf, I. and Nitsch, E.: Das Geologische Landesmodell von Baden-Württtemberg: Datengrundlagen, technische Umsetzung und erste geologische Ergebnisse, Freiburg i. Br., LGRB-Informationen, 21, 81 pp., 2008. 

Scheck-Wenderoth, M. and Lamarche, J.: Crustal memory and basin evolution in the Central European Basin System – new insights from a 3D structural model, Tectonophysics, 397, 143–165,, 2005. 

Schintgen, T.: The Geothermal Potential of Luxembourg – Geological and thermal exploration for deep geothermal reservoirs in Luxembourg and the surroundings, Ph. D. thesis, Universität Potsdam, Potsdam, 313 pp., 2015. 

Schmid, S. M., Fügenschuh, B., Kissling, E., and Schuster, R.: Tectonic map and overall architecture of the Alpine orogen, Ecl. Geolog. Helv., 97, 93–117,, 2004. 

Schmitt, D. R., Currie, C. A., and Zhang, L.: Crustal stress determination from boreholes and rock cores: Fundamental principles, Tectonophysics, 580, 1–26,, 2012. 

Sheorey, P. R.: A theory for In Situ stresses in isotropic and transverseley isotropic rock, Int. J. Rock Mech. Min, 31, 23–34,, 1994. 

Simpson, R. W.: Quantifying Anderson's fault types, J. Geophys. Res., 102, 17909–17919,, 1997. 

Smart, K. J., Ofoegbu, G. I., Morris, A. P., McGinnis, R. N., and Ferrill, D. A.: Geomechanical modeling of hydraulic fracturing: Why mechanical stratigraphy, stress state, and pre-existing structure matter, AAPG Bull., 98, 2237–2261,, 2014. 

Smith, W. H. F. and Sandwell, D. T.: Global Sea Floor Topography from Satellite Altimetry and Ship Depth Soundings, Science, 277, 1956–1962,, 1997. 

Sommaruga, A.: Décollement tectonics in the Jura forelandfold-and-thrust belt, Mar. Petrol. Geol., 16, 111–134,, 1999. 

Sperner, B., Lorenz, F., Bonjer, K., Hettel, S., Müller, B., and Wenzel, F.: Slab break-off – abrupt cut or gradual detachment? New insights from the Vrancea Region (SE Carpathians, Romania), Terra Nova, 13, 172–179,, 2001. 

StandAG: Gesetz zur Suche und Auswahl eines Standortes für ein Endlager für hochradioaktive Abfälle, 2017. 

Stromeyer, D. and Heidbach, O.: Tecplot 360 Add-on GeoStress, GFZ Data Services,, 2017. 

Tašárová, Z. A., Fullea, J., Bielik, M., and Środa, P.: Lithospheric structure of Central Europe: Puzzle pieces from Pannonian Basin to Trans-European Suture Zone resolved by geophysical-petrological modeling, Tectonics, 35, 722–753,, 2016. 

Tingay, M., Bentham, P., Feyter, A. de, and Kellner, A.: Present-day stress-field rotations associated with evaporites in the offshore Nile Delta, GSA Bull., 123, 1171–1180,, available at:, 2011. 

Turcotte, D. L. and Schubert, G.: Geodynamics, 3rd Edn., Cambridge Univ. Press, Cambridge, 623 pp., 2014. 

Valasek, P. and Mueller, S.: A 3D tectonic model of the Central Alps based on an integrated interpretation of seismic refraction and NRP 20 reflection data, in: Deep structure of the Swiss alps: results of NRP 20, edited by: Pfiffner, O. A., Lehner, P., Heitzmann, P., Mueller, S., and Steck, A., Birkhauser Verlag, Basel, 302–325, 1997. 

Valley, B. and Evans, K. F.: Stress State at Soultz-Sous-Forêts to 5 km Depth from wellbore failure and hydraulic observations, in: Thirty-Second Workshop on Geothermal Reservoir Engineering, 22–24 January 2007. 

Wagner, G. A., Coyle, D. A., Duyster, J., Henjes-Kunst, F., Peterek, A., Schröder, B., Stöckhert, B., Wemmer, K., Zulauf, G., Ahrendt, H., Bischoff, R., Hejl, E., Jacobs, J., Menzel, D., Lal, N., van den haute, P., Vercoutere, C., and Welzel, B.: Post-Variscan thermal and tectonic evolution of the KTB site and its surroundings, J. Geophys. Res., 102, 18221–18232,, 1997. 

Wagner, M., Kissling, E., and Husen, S.: Combining controlled-source seismology and local earthquake tomography to derive a 3-D crustal model of the western Alpine region, Geophys. J. Int., 191, 789–802,, 2012. 

Warners-Ruckstuhl, K. N., Govers, R., and Wortel, R.: Tethyan collision forces and the stress field of the Eurasian Plate, Geophys. J. Int., 195, 1–15,, 2013. 

Wenzel, F. and Brun, J. P.: A deep reflection seismic line across the Northern Rhine Graben, Earth Planet. Sc. Lett., 104, 140–150,, 1991. 

Wessel, P. and Smith, W. H. F.: A global, self-consistent, hierarchical, high-resolution shoreline database, J. Geophys. Res., 101, 8741–8743,, 1996. 

Ziegler, P. A. and Dèzes, P.: Crustal evolution of Western and Central Europe, Geological Society, London, Memoirs, 32, 43–56,, 2006. 

Ziegler, M. O. and Heidbach, O.: Matlab Script Stress2Grid v1.1, GitHub [code],, 2017.  

Ziegler, M. O. and Heidbach, O.: The 3D stress state from geomechanical–numerical modelling and its uncertainties: a case study in the Bavarian Molasse Basin, Geothermal Energy, 8,, 2020. 

Ziegler, M. O., Heidbach, O., Reinecker, J., Przybycin, A. M., and Scheck-Wenderoth, M.: A multi-stage 3-D stress field modelling approach exemplified in the Bavarian Molasse Basin, Solid Earth, 7, 1365–1382,, 2016. 

Ziegler, M. O., Ziebarth, M., and Reiter, K.: Python Script Apple PY v1.0, GFZ Data Services [code],, 2019. 

Zoback, M. L.: First- and second-order patterns of stress in the lithosphere: The World Stress Map Project, J. Geophys. Res.-Sol. Ea., 97, 11703–11728,, 1992. 

Zulauf, G.: Brittle deformation events at the western border of the Bohemian Massif (Germany), Geol. Rundsch., 82, 489–504,, 1993. 

Short summary
Knowledge about the stress state in the upper crust is of great importance for many economic and scientific questions. However, our knowledge in Germany is limited since available datasets only provide pointwise, incomplete and heterogeneous information. We present the first 3D geomechanical model that provides a continuous description of the contemporary crustal stress state for Germany. The model is calibrated by the orientation of the maximum horizontal stress and stress magnitudes.