3D crustal stress state of Western Central Europe according to a data-calibrated geomechanical model-first results

The contemporary stress state in the upper crust is of great interest for geotechnical applications and basic research likewise. However, our knowledge of the crustal stress field from the data perspective is limited. For Western 10 Central Europe 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 (Heidbach et al., 2018) as well as a complementary compilation of stress magnitude data of Germany and adjacent regions (Morawietz et al., 2020). 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 Western 15 Central Europe. The model covers an area of about 1000 x 1250 km and extends to a depth of 100 km containing seven lithostratigraphic 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 modelled orientations of SHmax are almost entirely within the uncertainties of the WSM data 20 used and the Shmin magnitudes fit to various datasets well. Only the SHmax magnitudes show locally significant deviations, primarily indicating too low values 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. 25


Basics of the crustal stress state
The stress state at a given point can be described by a second rank tensor ( Fig. 2a) with Pascal (1 Pa = N m -2 ) as the basic unit. Due to its symmetry properties, only six out of nine components are independent from 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 85 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).

95
Assuming that the vertical stress (the overburden) in the Earth crust SV = g··z 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.

Data compilation of stress tensor components 100
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., 2018). 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;105 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;Reiter et al., 2015). The new data have been integrated in the latest WSM database release (Heidbach et al., 2018). For stress magnitude data, Morawietz et al. (2020) published a publicly accessible database with 568 data records including a quality assessment 110 of the data for Germany and adjacent regions. These two datasets, the SHmax orientation of the WSM with some additional data of Levi et al. (2019) and the stress magnitude database are used to calibrate the geomechanical model (Fig. 1a).

Previous models
Modelling the contemporary crustal stress in Western Central Europe has been addressed by various authors since the mid-80s. However, except the model of Buchmann and Connolly (2007) which provide a 3D model of the broader Upper Rhine 115 Graben (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 with a wide range have https://doi.org/10.5194/se-2020-199 Preprint. Discussion started: 14 December 2020 c Author(s) 2020. CC BY 4.0 License. been tested so far and also various boundary conditions have been applied. For a detailed overview we refer to Cacace (2008) and Jarosiński et al. (2006). 120 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 Tornquist-Teisseyre Zone (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 local variation or second-order stresses, respectively (Jarosiński et al., 2006;Kaiser et al., 125 2005). An indirect influence of the depth position of the asthenosphere lithosphere boundary (LAB) and the resulting temperature contrasts and changed mechanical properties are described by Cacace (2008).

Conceptual modelling 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 135 parametrized 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. 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. 140 https://doi.org/10.5194/se-2020-199 Preprint. Discussion started: 14 December 2020 c Author(s) 2020. CC BY 4.0 License.
The model geometry reflects the contemporary distribution of rock properties such as density and stiffness and Poisson's ratio. An appropriate initial stresses 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 observed SHmax orientation is perpendicular or parallel to them (Fig. 1a). 145 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. 150    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 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 195 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, 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 200 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 205 boundaries between the Variscan basement units are simplified as vertical.

Model discretization
Our final mesh shown in Fig

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

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), which can be considered as being representative for tectonically inactive regions with low lateral density contrasts. The resulting stress state has to be in equilibrium with gravity and represents the initial stress state, which is subsequently perturbed by applying displacement boundary conditions that 240 impose the tectonic stress. (Buchmann and Connolly, 2007;Hergert and Heidbach, 2011;Reiter and Heidbach, 2014)

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 perpendicular to these (Fig. 3) to parametrize past and ongoing tectonic kinematics. The orientations of the model boundaries are chosen parallel or perpendicular to the 245 observed SHmax orientation. The eastern and western lateral model boundaries are aligned parallel and the northern and southern boundaries perpendicular to the SHmax orientation. Accordingly, extension is applied to the eastern and western boundaries and shortening to the northern and southern ones.
We use a two-stage approach to find a good agreement with the stress orientation and stress magnitude datasets. First a bestfit with respect to the SHmax orientation is estimated by an appropriate ratio between the extension and shortening applied. In 250 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 et al. (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 is applied. 255

Orientation of the maximum horizontal stress (SHmax)
We compare our model results with the stress orientation (WSM quality A to C) from the WSM database (Heidbach et al., 2018) and some additional data by Levi et al. (2019) from western Austria (Fig. 1a). However, we do not use individual data records for this comparison, but a mean SHmax orientation is used on a regular 0.5° grid. Using a mean SHmax orientation 260 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 for 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 https://doi.org/10.5194/se-2020-199 Preprint. Discussion started: 14 December 2020 c Author(s) 2020. CC BY 4.0 License. distance to the grid points. Each grid point requires at least ten data points within a fixed search radius of 200 km.
Furthermore, we only use mean SHmax orientations that have a standard deviation <25° using the statistics of bi-polar data 265 (Mardia, 1972). The grid contains a total number of 619 data points. 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. Carpathians show 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. 6d where the histogram of the deviation between the mean SHmax orientation derived from the WSM data and the modelled orientation is shown with a median deviation of 6.9° and a mean deviation, calculated from the absolute differences, of 15.4°. 280

Minimum horizontal stress (Shmin) magnitudes
The modelled magnitudes of Shmin in comparison to stress magnitude data of Morawietz et al. (2020) are shown in Fig. 7.
The figure is divided into three subfigures displaying the differences depending on depth and quality (Fig. 7a), the spatial distribution of the calibration data (Fig. 7b) and a histogram showing the distribution of the differences between the 290 modelled and observed Shmin magnitudes (Fig. 7c). The differences are calculated as interpolated model results minus data, thus positive differences correspond to too large model values and negative ones to too low model values. We use only data from Morawietz et al. (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.   7b the data are mainly located within the south-western part with the exception of one measurement from the NGB. With 42 295 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 et al. (2020). 300 https://doi.org/10.5194/se-2020-199 Preprint. Discussion started: 14 December 2020 c Author(s) 2020. CC BY 4.0 License.

Maximum horizontal stress (SHmax) magnitudes
For the model calibration regarding the SHmax magnitudes 57 data records are used from the database of Morawietz et al. (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. 8b). The data are from seven different localities, whereby 310 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. 8a). There are significantly more data records from shallower depths (200 to 1000 m) which indicate too large model results than from the greater depths (>1000 m) which indicate too low results. Regardless of this, a trend is visible 315 from positive to negative stress differences with increasing depth, i.e. the model seems to predict too large values of SHmax in the upper part of the model and too low values of SHmax 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.

Stress gradients and stress regime
Additional to the calibration of the model with stress magnitude data, the absolute stress magnitudes of SV, Shmin and SHmax for three hypothetical wells up to 10 km depth are shown in Fig. 9. We have chosen these three locations (Fig. 10) 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 330 basement, well 2 entirely within the sedimentary unit and well 3 partly within the sediment unit and partly within the   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 shows thus the only significant 355 https://doi.org/10.5194/se-2020-199 Preprint. Discussion started: 14 December 2020 c Author(s) 2020. CC BY 4.0 License. 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.

of SV (yellow), Shmin (blue) and SHmax (red). Blue stars and red 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 Rittershofen. Measured Shmin magnitudes of Valley and Evans (2007), calculated SHmax values of Klee and Rummel (1993) and calculated SV magnitudes
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. 360 As an additional result the regime stress ratio (RSR) (Simpson, 1997) for four model sections are shown in Fig. 10. The RSR (Eq. 1) is a unitless value combining the stress regime index n (Eq. 2) of Anderson (1905) and the ratio of stress differences ϕ (Eq. 3): = ( + 0.5) + (−1) ( − 0.5) (Simpson, 1997) (1) n = { 0 S hmin < S Hmax < S V 1 S hmin < S V < S Hmax 2 S V < S hmin < S Hmax (Anderson, 1905) (2) 365 (Angelier, 1979) (3) The resulting value between 0 and 3 is the RSR indicating the stress regime divided into 6 classes from radial extension over extension, transtension, transpression and compression to constriction. The results in Fig. 10 show with the exception of peripheral areas 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 at 10 km depth. The two sections in between at 2000 and 370 4000 m show the transition with a dominant RSR of 1 to 1.5 and 1 to 0.5, respectively. https://doi.org/10.5194/se-2020-199 Preprint. Discussion started: 14 December 2020 c Author(s) 2020. CC BY 4.0 License.

Orientation of the maximum horizontal stress (SHmax)
The results of the SHmax orientation are in comparison to the mean WSM data quite good (Fig. 6). With a median deviation of 6.9° and an absolute mean deviation of 15.4° 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). However, there are three areas with noticeable deviations: In the north-eastern part of Germany, in Belgium and along the Carpathians. The region in north-eastern 380 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. 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. Heidbach et al., 2007;Roth and Fleckenstein, 2001;Röckel and Lempp, 2003;Grote, 1998;Ahlers 385 et al., 2019). However, since we do not distinguish between the data from these different layers the derived mean SHmax https://doi.org/10.5194/se-2020-199 Preprint. Discussion started: 14 December 2020 c Author(s) 2020. CC BY 4.0 License. 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 390 contrasts are often considered to be responsible for reorientations in the stress field (e.g. Grünthal and Stromeyer, 1994;Heidbach et al., 2007;Marotta et al., 2004;Reiter, 2020). 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 395 horizontal stresses resulting in a near isotropic stress state (Heidbach et al., 2007) and thus the topography contrast between the mountain range and the Pannonian Basin has probably 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. 400 In general, our modelled 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 perturbations are not implemented. Nevertheless, our model results also show no impact of mechanical contrasts on the orientation of SHmax, despite lateral mechanical 405 contrast, e.g. a Young's modulus difference of 40 GPa between the sedimentary (30 GPa) and the upper crustal units (70 GPa).

Magnitudes of Shmin, SHmax and SV
The Shmin magnitudes (Fig. 7a) in general show a very good correlation with the data of Morawietz et al. (2020) independent of data quality and depth. The statistical evaluation supports this with a mean difference of 3.3 MPa, a median difference of 410 0 MPa and an almost even distribution. 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. 8). 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 415 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 too large SHmax magnitudes are from shallow depths (200 to 1000 m). 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 420 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. 9 show by and large the expected results. Since SV is 425 only dependent on the density the gradients of well 1 and well 3 located in a single unit are constant all over the total well depth. The gradient of well 2 changes in between 1500 and 3500 m depth due to the change of the lithological unit. The transition zone is quite large because of the vertical element resolution of about 800 m. Based on the sum of the overburden, https://doi.org/10.5194/se-2020-199 Preprint. Discussion started: 14 December 2020 c Author(s) 2020. CC BY 4.0 License. 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 430 located in. Since the Poisson's ratio is constant for the units involved (0.25) the 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 the Young's modulus. Well 2 shows the smallest stress differences of 15 MPa and well 1 shows differences of 35 MPa according to a 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. 435 To get a more detailed insight we compare our hypothetic wells 1, 2 and 3 with local magnitude data (Fig. 9). 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 the difference increases 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 440 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)  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 are the values of the Young's modulus. The calculated values of 450 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. 10), which is contrary to the stress regime of Brudy et al. (1997) showing a strike-slip regime from 1 km depth downwards. This indicate that our SHmax 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 455 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 too low horizontal stresses in this model region, which may explain the 460 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 465 from depths shallower than 3500 m ( Fig. 7 and 8), 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. 7. Due to missing data, a comparison is not possible for the SHmax values. However, Röckel and Lempp (2003)  https://doi.org/10.5194/se-2020-199 Preprint. Discussion started: 14 December 2020 c Author(s) 2020. CC BY 4.0 License.
For the comparison of our hypothetic well 3 we use data from the geothermal project in Soultz-sous-Forêts and some values from Rittershofen, 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. 9 are calculated on the base of density values of Rittershofen (Azzola et al., 2019) 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 480 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)  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 to low SHmax 490 values.
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. 9). The modelled 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 be probably only possible with a higher stratigraphic resolution in the sedimentary unit and a density 495 gradient within the upper crust. Such a simple gradient, which could also reduce the differences in the sedimentary unit, would be in this case less useful, 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. 7, but also with regard to the Soultz-sous-Forêts and NGB location in Fig. 9. Only the results for the KTB site show some considerable differences, with the greatest deviations from calculated Shmin values (3000 to 7000 m depth) 500 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. 8) and in the local comparisons (Fig. 9). 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 too low values in the deeper parts, but not the trend of too large values in the shallower parts. In contrast, the SHmax magnitudes at Soultz-sous-Forêts show a good agreement down to 2000 m depth 505 and the values at the KTB site even indicate too low magnitudes down to 800 m depth. An indication for generally too small values with increasing depth are 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 510 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. https://doi.org/10.5194/se-2020-199 Preprint. Discussion started: 14 December 2020 c Author(s) 2020. CC BY 4.0 License.

Conclusion
The model presented is the first 3D geomechanical model for Western Central Europe predicting the first order 3D stress tensor. The model is calibrated with SHmax orientations from the WSM database and compilation of Shmin and SHmax stress 515 magnitude data from Morawietz et al. (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 extend 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 520 Evans, 2007), but the SHmax magnitudes result show locally significant differences. Modelled SHmax magnitudes are too small in the lower part of the model, whereas some results indicate too high values 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 525 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 modelled ones. In addition to a vertical refinement, resolving lateral variations of the rock properties would be useful as well as these potentially account for lateral variability of the stress tensor. 530