The 3D stress state from geomechanical–numerical modelling and its uncertainties: a case study in the Bavarian Molasse Basin

Knowledge of the undisturbed stress state is a key parameter for borehole stability, productivity and induced seismicity hazard assessment. Due to the sparse and incomplete availability of data records, 3D geomechanical–numerical modelling is applied to estimate the stress state in a volume. To assess the quality of the model results, the model uncertainties have to be quantified. We present an approach that provides the uncertainties of the six independent stress tensor components at each location in the model volume by an average value and its standard deviation. We explore the uncertainties introduced to the model results with respect to the stress magnitude data used for model calibration. We test our approach on a model of the Bavarian Molasse Basin which includes the area around Munich with many geothermal projects. In the test area, we find large uncertainties in the modelled magnitude of the maximum horizontal stress (SHmax\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text{S}_\text {Hmax}$$\end{document}) in the order of 15–30%. The uncertainties in the magnitude of the minimum horizontal stress (Shmin\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text{S}_\text {hmin}$$\end{document}) are smaller between 5 and 20%. In connection with an adequate failure criterion, we compare our model results to the seismological observations at two neighbouring geothermal projects. Whilst Aschheim/Feldkirchen/Kirchheim remained seismically quiet, induced seismicity was recorded in Poing. Our modelled undisturbed stress state was confirmed by these observations by showing in average a stable stress state in Aschheim/Feldkirchen/Kirchheim and a critical stress state in Poing. If the uncertainties in terms of two standard deviations are considered, this does not change. This demonstrates the increase in model significance when uncertainties are available.

these events are above a site-specific value-after (Grünthal 2014) these events are labelled as "Seismic Events of Economic Concern" SEECo-geothermal exploitation is often stopped. Prominent examples of discontinued projects are, e.g. Basel, St. Gallen, and Pohang (Häring et al. 2008;Diehl et al. 2017;Grigoli et al. 2018). To identify a critical stress state, three components need to be determined: (a) The failure criterion of the rock, (b) the initial stress state, which is the undisturbed natural stress state before anthropogenic subsurface activities, and (c) the induced stress changes due to drilling, stimulation, production and re-injection. Only the superposition of these three components allows the determination of the spatio-temporal evolution of the absolute rock stability during the entire lifetime of a geothermal project.
Adequate description of rock failure exists (Jaeger and Cook 1971;Zang and Stephansson 2010;Hoek and Martin 2014) and the modelling of stress changes due to pore pressure changes is commonly performed (Van Wees et al. 2014;Altmann et al. 2010;Hornbach et al. 2015;Jacquey et al. 2016). Yet, without the knowledge of the initial stress state only relative stress changes can be estimated (Morris et al. 1996;Connolly and Cosgrove 1999;Henk and Nemčok 2008). To describe the initial stress state, 3D geomechanical-numerical models are used (Reiter and Heidbach 2014;Van Wees et al. 2003;Fischer and Henk 2013;Cacace et al. 2013). These models are calibrated against model-independent stress magnitude data and information on the stress regime, but since these data are in general sparse and incomplete in the model volume, the model results that describe the initial stress state have presumably large uncertainties. However, a rigorous quantification of this uncertainty has not been investigated. In using the term uncertainty we refer to the common definition (see, e.g. Begg et al. 2014) in contrast to the term variability. Uncertainty refers to the "likelihood of what the single, true value of the uncertain quality is" (Begg et al. 2014). Variability refers to the range of "multiple instances of the quantity, derived from observed data" (Begg et al. 2014).
In this paper, we focus on the rigorous estimation of the uncertainties of 3D geomechanical-numerical models that describe the initial stress field. We aim at provision of the average stress state out of all possible stress states supported by data records. Our approach is purely data driven, which means that we do not know a priori the individual error of each data record that is used for model calibration. This is in contrast to Lecampion and Lei (2010) who use reliable data records with known measurement error and thus, an inversion scheme can be applied to minimise the standard deviation. As a demonstration model, we use an updated and extended 3D geomechanical-numerical model of the Molasse Basin of the Alpine foreland , Figure 1). For the model calibration, we use 13 data records on the magnitudes of the minimum and maximum horizontal stress S hmin and S Hmax , respectively, from within the model volume. Instead of using the average of all data records to generate a single model scenario, we run 9 × 4 = 36 model calibration scenarios that cover all combinations of pairs of the 9 S hmin and 4 S Hmax data records. For each scenario, we compute the 3D stress tensor throughout the entire model to obtain a 3D stress field. From the 36 scenarios the average and its standard deviation are computed for the full stress tensor. In contrast to previous assessments of uncertainties , we focus herein on the variability of stress magnitude data records assuming that this is the only source of model uncertainties.
We obtain the average of each component of the full stress tensor and their standard deviations in 3D, i.e. the stress state and its uncertainties show lateral and vertical fluctuation due to lithological boundaries and material properties. Our model results have a prevailing strike-slip stress regime in general and in the relevant units for geothermal exploration in the Molasse Basin in particular. Furthermore, we show the modelled stress state at two geothermal sites of interest throughout the model area. Whilst one site remained seismically quiet (Aschheim/Feldkirchen/Kirchheim, AFK) during geothermal exploration, the other showed seismicity of up to Ml 2.1 (Poing) (Megies and Wassermann 2017). In combination with reasonable failure criteria, our modelled stress state and these observations are reciprocally confirmed.

In situ stress state
The stress state can be mathematically described by the symmetric second-order stress tensor The stress tensor has six independent components (three normal components σ ij , i = j and three shear components τ ij , i = j ). In geomechanics, the stress tensor is usually transformed into its principal axis system in which the shear stresses vanish and only the three principal stresses σ 1 , σ 2 , and σ 3 , that are perpendicular to each other, remain It is assumed that in the upper crust one of the principal stresses is the vertical stress S v . This assumption reduces the unknowns to four: the magnitudes of S v and two horizontal principal stresses ( S Hmax and S hmin ), and the orientation of one of the horizontal principal stresses, usually S Hmax . The relation between the three principal stress components' magnitudes defines the stress regime. Three main stress regimes are distinguished: normal faulting ( S v > S Hmax > S hmin ), strike-slip ( S Hmax > S v > S hmin ), and thrust faulting ( S Hmax > S hmin > S v ).
The orientation of S Hmax is derived by methods such as the analysis of borehole breakouts or drilling-induced fractures, the principal strain axes of earthquake focal mechanisms, or geological indicators, amongst other methods (Amadei and Stephansson 1997;Schmitt et al. 2012). Data records on the orientation of S Hmax are compiled and quality ranked in the World Stress Map (WSM) Zoback 1992;Zoback et al. 1989). In the model area, mainly data records from hydrocarbon exploratory boreholes are available . Their data provide a good coverage of the S Hmax orientation that shows a prevailing pattern of almost N-S orientation Reiter et al. 2015;Reinecker et al. 2010). The information of the S Hmax orientation is important for the choice of the orientation of the model boundaries as they have to be parallel or perpendicular, respectively, to the horizontal stress orientations Heidbach et al. 2018).
Information on the vertical stress magnitude S v is most easily obtained by with the depth z, the depth-dependent rock density ρ , and the gravitational acceleration g. Estimates on the rock density ρ are available from various sources. Lithological borehole logs can be correlated with laboratory density measurements of analogue rock materials. If borehole cores are available, the density can be directly measured. Finally, density borehole logs may be available. Information on the horizontal stress magnitudes are usually sparse since these observations are required to be made in situ which accounts for high cost. Hydraulic fracturing (HF) tests provide comprehensive information on the magnitude of S hmin (Haimson and Fairhurst 1969;Haimson and Cornet 2003). Some further data on the S hmin magnitude can be derived from borehole tests that are often conducted by default to test borehole stability. During a Leak-Off Test (LOT), the fluid pressure is increased in the borehole until the fluid leaks into the formation (Bell 1996;Addis et al. 1998). The fluid pressure at which the leak-off occurred is used as a proxy for the S hmin magnitude, even though it is not as reliable as an HF test. During a Formation Integrity Test (FIT) the fluid pressure is increased, but only up to a pre-defined pressure. Reaching this level without any leak-off proves the boreholes integrity (White et al. 2002). The maximum pressure achieved during an FIT is used as a lower bound for S hmin (Zoback et al. 2003). These three methods (HF, LOT, FIT) only give estimates of the S hmin magnitude when S v is not the minimum principal stress. Therefore, in thrust faulting stress regime these tests are not useful.
The S Hmax magnitude can be derived from HF tests, when details on the fluid pressure curve are given (Haimson and Cornet 2003;Zoback 2010). However, different approaches exist and lots of assumptions are required which reduces the significance of the resulting S Hmax magnitude (Ito et al. 1999;Zoback 2010;Zoback et al. 2003). Another approach is to use the frictional limit assumption (Jaeger and Cook 1969). Hence, it is assumed that the Earth's crust contains pre-existing faults that are optimally oriented in the prevailing stress field and that these faults are at their frictional limit (Jaeger and Cook 1971;Sibson 1974). Then, the ratio of the maximum and minimum principal stress is determined with σ 1 /σ 3 = ((µ 2 + 1) + µ) 2 . With a friction coefficient of µ = 0.6 this ratio would be 3.1 (Jaeger and Cook 1971;Jaeger et al. 2007;Sibson 1974). When the stress regime is known and with the assumption that S v is a principal stress, upper bounds for S Hmax in a thrust faulting and lower bounds for S hmin in normal faulting stress regime, respectively, can be made. In a strike-slip stress regime where S v is the intermediate principal stress, the ratio is determined by S Hmax /S hmin and further assumptions or information is needed. With these the stress state can be narrowed with a so-called stress polygon (Schoenball and Davatzes 2017;Zoback 2010;Zoback et al. 2003). Here, additional information, e.g. from FITs or LOTs can be introduced to further constrain the upper and lower boundaries of the principal horizontal stresses as shown by Seithel et al. (2015) for the borehole Sauerlach south of Munich in the Bavarian Molasse.

Geological setting and seismotectonics of the model area
The Northern Alpine Molasse Basin extends along the mountain front approx. 1100 km from Lake Geneva in the West to Vienna in the East. Its largest N-S extend with approx. 140 km is in Bavaria/Germany. The lithological units dip south towards the Alps. At the southern border of the Molasse Basin in front of the mountain range, a maximum thickness of the sediments of approx. 5000 m is recorded. There, the basement consists of crystalline crust which is overlain by a lithological succession of Mesozoic units. These units are themselves overlain by the sedimentary Folded Molasse in front of the mountain range and the Foreland Molasse farther north (Lemcke 1978). Regarding the contemporary tectonic regime in the Bavarian Molasse Basin, older studies proposed that it is between strike-slip and thrust faulting (Illies and Greiner 1978;Lemcke 1988). More recent studies, however, favour normal faulting (Bachmann et al. 1987;Drews et al. 2018Drews et al. , 2019Reinecker et al. 2010), normal faulting to strike-slip (Budach et al. 2017;Backers et al. 2017) or strike-slip (Seithel et al. 2015;Megies and Wassermann 2014).
In the Upper Jurassic karst lithology-formerly referred to as Malm-hot water is present and provides a large potential for geothermal applications (Agemar et al. 2014). With a geothermal gradient of 30-40 °C/km, a temperature of > 100 °C is reached at a depth of 3000 m (Agemar et al. 2012). Even though this temperature is not sufficiently high for extensive electric power generation, the water is used for applications such as district heating and balneology. This is reflected in the current usage of geothermal energy in the Molasse Basin, which has an installed thermal capacity of more than 323 MW th , but an electric capacity of only little more than 35 MW e as of 2019 (Bundesverband Geothermie 2019).
Before the onset of geothermal exploitation, no natural seismicity has been recorded in the area (Grünthal 2011(Grünthal , 2014Grünthal and Wahlström 2012). This is partly due to the limited sensitivity of the seismological stations in Bavaria before an improvement in station density in 2000 and the following years. With this network and additional installations at geothermal sites, low-magnitude seismicity can be recorded. So far, at eight geothermal projects in Bavaria an onset of seismicity was recorded after production had started (Megies and Wassermann 2017). The largest events were recorded in Unterhaching (Ml 2.4), Poing (Ml 2.1), Sauerlach (Ml 1.2) (Megies and Wassermann 2017), and Dürrnhaar (Ml 2.0) (Wassermann and Megies 2019). At the other geothermal project sites, no seismicity has been recorded or felt by the population so far. Furthermore, 48 seismic events could not be associated to a particular project site since the location could not be determined accurately enough since their magnitudes were all below Ml ≤ 0 (Megies and Wassermann 2017).

Model area and basic assumptions
Our model focus is on the part of the Molasse Basin that currently experiences most activity in geothermal energy applications and exploration. It is located between the cities of Augsburg and Salzburg (E-W) and between Freising and the Alps (N-S). The area extends 145 × 70 km 2 (E-W, N-S, Fig. 1) and reaches from the ground surface to a depth of 11 km below sea level. The lithology is based on the geological model of Przybycin et al. (2017Przybycin et al. ( , 2014 that uses freely available lithological information from boreholes and seismic lines to describe the subsurface geometry in the Molasse Basin. Furthermore, the geometry has been constrained by 3D gravity modelling (Przybycin et al. 2017). A total of 13 geologic units are distinguished (Table 1) with a focus on the six thin southward dipping Jurassic units of Malm alpha to Malm zeta, which are target formations for geothermal applications (Lemcke 1988;Bachmann et al. 1987;Fritzer et al. 2012). The model area is approximately E-W/N-S oriented with a small clockwise rotation by 1.7° to align the model boundaries parallel to the prevailing regional S Hmax orientation .
To model the initial stress state, we follow the two-step approach of applying gravity followed by lateral displacement as explained in detail in Hergert et al. (2015). The model rheology is linear elastic and rock properties are the same as given in Ziegler et al. (2016) and listed in Table 1. To solve the resulting 3D partial differential equation that describes the equilibrium of the volume forces (gravity) and surface forces (displacement boundary conditions), we use the finite element method. Thus, the geologic model is discretised with 2.95 × 10 6 8-node hexahedron elements. Each element has a lateral extend of 1500 × 1500 m 2 which approximately corresponds to the resolution of the geologic model (Przybycin et al. 2017(Przybycin et al. , 2014. Despite this comparably The resulting high vertical resolution allows resolving the mechanical differences between the six thin Malm units. The assignment of elements to lithological units is done automatically by the algorithm Apple PY (Ziegler et al. 2019). Afterwards, the elements are populated with the rock properties, i.e. the Young's modulus, Poisson ratio, and density according to the lithological units (Table 1).

Model calibration
Whilst the vertical stress component is adequately initiated by the density of the rock, the horizontal stresses depend on the uniform displacements assigned perpendicular to the lateral model boundaries. Since the model boundaries are oriented parallel to the mean S Hmax orientation, the displacement in each of the two orientations (here, approx. E/W and approx. N/S) is used individually to adjust S Hmax and S hmin , respectively (Hergert et al. 2015;Reiter and Heidbach 2014). A variety of displacement boundary conditions exist which calibrate the S hmin magnitude (Fig. 2a) or the S Hmax magnitude. However, there is only one pair of displacement boundary conditions that fit both, the S Hmax and S hmin magnitude data records (Fig. 2b). If more than one data record for the calibration of each S hmin or S Hmax is available, boundary conditions that result in a model that fits the observations best are sought (Fig. 2c). The boundary conditions that minimise the mean difference between modelled stress state and the calibration data records are called best-fit boundary conditions of the best-fit model scenario (Fig. 3). This entire process is referred to as model calibration and is assisted by the FAST Calibration tool which is documented in detail in Ziegler et al. (2016) and Ziegler (2018). At the end of the calibration procedure, the best-fit stress state is available at each node in the model.
A prerequisite for this modelling approach is the assumption that the available magnitude data are consistent in first order, i.e. two data records that are located at different lateral locations but in the same true vertical depth and lithology are expected to have slightly but not systematically different stress magnitudes. These slight deviations are due to lithological inhomogeneities and/or measurement errors related to the stress indication methods. Errors in the estimation or derivation of stress magnitudes (especially S Hmax ) may further blur their spatial and depth distribution. The calibration data availability in the model area is indicated in Table 2 and Fig. 1. In Freiham and Geretsried (Backers et al. 2017), Unterhaching (Budach et al. 2017;Backers et al. 2017), and Sauerlach (Seithel et al. 2015), estimates on the magnitude of S Hmax and S hmin are available. Budach et al. (2017) and Backers et al. (2017) provide two completely different sets of stress magnitude data for either a strike-slip ( S Hmax > S v > S hmin ) or transtensional ( S Hmax = S v > S hmin , there called normal faulting) stress regime in Freiham, Geretsried, and Unterhaching. Due to the observed induced seismicity, Megies and Wassermann (2014), Budach et al. (2017) and Backers et al. (2017) prefer the strike-slip stress regime. In addition, the S hmin magnitude is constrained by two LOTs in Unterhaching (Drews et al. 2019), a suspected LOT during an FIT in Sauerlach which is treated as an LOT (pers. comm. T. Fritzer, Ziegler et al. 2016), an LOT in Unterföhring, and an S hmin magnitude data record from an undisclosed method in Kirchweidach (Backers et al. 2017), assumed to be an LOT. Lower row shows the associated implementation in the model with stress magnitude data that are used for calibration. Red ( S Hmax ) and blue ( S hmin ) indicate the horizontal stress magnitudes and the associated possible boundary conditions that fit their respective value. a For a single S hmin data record the combinations of boundary conditions that satisfy the calibration data record is described by a linear function. One possible set of boundary conditions is indicated here. b An S hmin and S Hmax calibration data record can be modelled individually again by boundary conditions on a straight line each. However, if both S Hmax and S hmin should be satisfied the intersection of these two lines, indicated by a solid white dot in the top row, is the unique pair of boundary conditions that fits both. c If several S Hmax and S hmin data records are available, for each S Hmax and S hmin the linear function with the least mean difference between observation and model is derived. Then, the intersection of these two functions is used as the valid set of boundary conditions that minimises the differences between modelled and observed stress state Fig. 3 Best-fit modelling approach (left) vs. average stress state with uncertainty quantification approach (right) of geomechanical-numerical modelling. In a best-fit approach the difference between the modelled and observed stress magnitudes is minimised. The uncertainty quantification approach combines each S hmin with each S Hmax magnitude and thus generates several model scenarios that are each adjusted to fit an individual pair of calibration stress data records. From all of these model scenarios, the average and its standard deviation of the full stress tensor is computed at each node throughout the model volume Table 2 Disclosed stress magnitude data from the model area The asterisk signifies that the magnitudes are assumed for a strike-slip regime. For the corresponding magnitudes for a transtensional stress regime, see Table 3 Stress component

Estimation of model uncertainty
To estimate the uncertainty of the modelled stress tensor, we have to investigate the influence of individual stress magnitude data records on the resulting stress tensor of the model. Herein, we call this the uncertainties of the modelled stress state since it describes the likelihood of the stress state at a certain location and depth in the model. Neglecting the uncertainties that may origin from our assumption that the rock volume can be represented by homogeneous, isotropic, linear-elastic continuum, the uncertainties in the modelled stress state are only due to variabilities in the a priori unknown stress data records that are used for the model calibration. Thus, to estimate the uncertainties in the modelled stress state, each S hmin data record is combined with each S Hmax data record and used for the calibration of an individual calibration scenario (Fig. 2b). This permutation results in 36 model scenarios. The boundary conditions for each model scenario are adjusted so that the data records agree with the observations (Fig. 2b). Thereby we explore the parameter space defined by the variations in the calibration data. Then, an average value for each component of the stress tensor is computed for each element in the model by with n = 36 as the number of model scenarios. Furthermore, the corresponding standard deviation of each component of the stress tensor is computed at each element throughout the model (Fig. 3).

Results and discussion
In the following four sub-sections, we present and discuss the model results from four different perspectives. As a reference, we first show and explain the key features of the model result and their uncertainties. In the second sub-section, we analyse the sources of the estimated model uncertainty. In the third sub-section, we show the results that we receive for the stress regime and discuss them in reflection of the published literature and data. In sub-section four, we link the model results to the key motivation, i.e. the stress state criticality and to what extent this model and the estimated uncertainties can explain the occurrence or absence of induced seismicity at different geothermal sites.

Average modelled stress state with uncertainty
The average modelled stress state σ ij is presented in Fig. 4 by the distribution of the S hmin magnitudes and in Fig. 5 with the S Hmax magnitudes. In both figures, the influence of the southward dipping Malm units on the stress state is indicated in the profile of A-A′. The increasing stress magnitude with depth is depicted in the map-view of the dipping horizon top Malm delta. The standard deviation of S Hmax is approximately 10 times the standard deviation of S hmin . The influence of the Landshut-Neuöttinger-Hoch on the stress state is indicated by a decrease in magnitudes (in the north east of the map view). Both S Hmax and S hmin show smaller averages and smaller standard deviations in the area of the Landshut-Neuöttinger-Hoch. Generally, the standard deviation is significantly higher in the basement compared to the sediments (Fig. 6). In both, basement and sediments, a general trend of an increase in standard deviation southward towards the mountains is observed (Fig. 6). This trend, differentiating the Foreland Molasse from the basement, stops within the Folded Molasse towards the mountain front which is also close to the boundary of the model. The average model result with its standard deviation is provided in the data publication (Ziegler and Heidbach 2020). The following variables are available: the coordinates and the vertical stress S v as well as the average and standard deviation of S Hmax and S hmin magnitudes, the regime stress ratio (RSR) after Simpson (1997), and the differential stress S 1 -S 3 where S 1 and S 3 are the maximum and the minimum principal stress, respectively.

Sources of uncertainties
The uncertainties in the modelled S hmin magnitude are comparably small (Fig. 4) in contrast to the uncertainties in the modelled S Hmax which are very high (Figs. 5,  6). The magnitude of S hmin is with nine data records clearly better constrained than the magnitude of S Hmax with only four data records. The reason is that S hmin data records used for calibration are more easily obtained. S Hmax calibration data records,  however, cannot be estimated from measurements and their derivation requires several assumptions on material properties and the stress state. Hence, the available data records differ significantly even for the same units and depths, which leads to a high data variability. The modelled magnitude of S v is not subject to uncertainties in this approach. S v is only a result of the overburden whose density is defined by the material properties. However, the variability of material properties of the individual units, including their density, is not regarded here.
An increase in model uncertainties may be observed if variability in material properties as well as lateral or vertical inhomogeneities within the same unit, that may affect all modelled stress magnitudes, are regarded. Available studies from other areas indicate large standard deviations in material properties (Reyer and Philipp 2014). In particular, this is significant for the locations of the calibration data points. Due to the modelling approach that uses displacements as boundary conditions to induce the stress state, it is imperative that the calibration data points are mechanically in the same lithology and depth in the model and in reality. Therefore, information on the stress state (HF, LOT, FIT) combined with triaxial laboratory tests to estimate the material properties from the same depth and lithology would be preferable but are rarely available. The numerical error due to the choice of linear finite elements and the discretisation has been checked by increasing the resolution to second-order finite elements. The deviation to the linear finite elements is in average < 1% with local maxima < 3%. In addition, the stress magnitude data records used for model calibration do not have a (reported) measurement error since usually only one estimation or measurement was conducted per location. Neither are measurement errors reported, which are likely present, and could be used to assigned an uncertainty to the data records. Knowledge of these would further improve the significance of the model even though the modelled uncertainties are likely to increase at the same time.
In addition to the estimation of the average stress state and its uncertainties, a bestfit model scenario is estimated following the approach of minimisation of differences between modelled and observed stress state. The best-fit approach uses all available data records without weighting them and thus follows the calibration procedure for a minimisation of the differences between model and observation displayed in Fig. 2c. In general, compared to the average stress state, the magnitude of S hmin is overestimated and the magnitude of S Hmax is underestimated by the best-fit model scenario. In the Molasse sediments, S hmin in the best-fit model scenario is approximately 0.3 MPa larger than the average S hmin . In the basement, the best-fit scenario overestimates S hmin by between 1 and 2 MPa. An underestimation of S Hmax by the best-fit model of up to 6 MPa is indicated. This combination of over-and underestimation of S Hmax and S hmin in particular in a strike-slip setting is meaningful since it decreases the differential stress S 1 -S 3 (Fig. 7). Thus, a Mohr circle derived from the best-fit model indicates a less critical stress state than that provided by the average of the 36 model scenarios. Thus, the reactivation of faults seems less likely if the stress state is derived from a best-fit model.
The applied method of uncertainty quantification founds on the assumption that the stress data records available for calibration cannot be graded according to their reliability in terms of a standard deviation. Even though the degree of believe in each stress data record is variable but can be quantified (Morawietz et al. in review), no measurement error or variability derived from several consecutive measurements is available. In particular, information on the data variability is important if the data records originate from a wide range of sources. Thus, the presented approach cannot aim at minimising the standard deviation (see, e.g. Lecampion and Lei 2010) but instead provides a data-based average stress state with uncertainties that is supported by all the available data records.

Tectonic stress regime
The tectonic stress regime characterises the stress state using the relation of the magnitudes of the three principal stresses. In the Bavarian Molasse Basin, the stress regime is highly debated (Drews et al. 2019;Reinecker et al. 2010;Seithel et al. 2015). Seithel et al. (2015) assume a strike-slip stress regime but Budach et al. (2017) and Backers et al. (2017) provide magnitude estimations of S hmin and S Hmax for either a strike-slip or transtensional stress regime for three locations. The authors prefer a strike-slip tectonic regime due to the observed induced seismicity in the Molasse Basin of which the ones with focal mechanism solutions showed a strike-slip mechanism Wassermann 2014, 2017). This indicator has to be regarded with caution since the recorded seismicity is most likely induced events. Thus, an unknown alteration of the undisturbed stress state took place before the events occurred, which limits the informative value in terms of the initial, i.e. the undisturbed stress state.
We tested the stress magnitude data for both a strike-slip and a transtensional ( S Hmax = S v > S hmin ) stress regime. Therefore, we applied a method comparable to that used by Drews et al. (2019) in comparing the stress state implied by the provided data records with Formation Integrity Tests (FITs) throughout the model area. Since an FIT does not fracture the rock, the maximum pressure achieved is a lower boundary for the least principal stress magnitude. Hence, if the modelled least principal stress at a certain location is lower than the recorded FIT fluid pressure, the model result is probably not correct in that it shows a stress state in which the FIT should have developed into a hydraulic fracture and loosing fluid into the formation. We used three stress states both for strikeslip (Table 2) and transtensional (Table 3) individually for calibration of the model. For each of the three transtensional stress states, the majority of recorded FIT fluid pressures are larger than the modelled σ 3 magnitude. In contrast, barely any FIT fluid pressure is larger than the modelled σ 3 magnitude for a modelled strike-slip stress regime. This indicates that the data reflecting a strike-slip regime have to be preferred over those reflecting a transtensional stress regime.
These findings of a preferred strike-slip regime are in contrast to Drews et al. (2019), who expect a normal faulting regime in the Molasse Basin. This discrepancy can have different reasons. (1) It can be due to the fact that we use more diverse density data for the modelling of S v . Hence, Drews et al. (2019) might have overestimated S v which changes a strike-slip regime to a normal faulting regime.
(2) They only address a shalerich lithology whilst our model includes a large part of the Molasse Basin with a variety of lithological units. Depth and lithology-dependent changes of the tectonic stress regime are likely to be present (Hergert et al. 2015). Hence, the findings of Drews et al. (2019) are not necessarily in contradiction to ours. At least a transtensional stress regime in the depth of the Malm units may closely be present within the model uncertainties at the Malm units which is reflected by the standard deviation ( 1σ ) in Fig. 8.

Geothermal sites and induced seismicity
We apply our model to investigate the design and operation of geothermal power plants. Therefore, we compare two neighbouring geothermal sites Aschheim/Feldkirchen/Kirchheim (AFK) and Poing that are approximately at 4 km distance. They are in operation since 2009 and 2011, respectively, and provide hot water for district heating (geothermie.de). Both geothermal projects have two boreholes, one for production and one for reinjection. AFK has a flow rate of 75 l/s and a target depth of 2630 m, and Poing has a slightly higher flow rate of 100 l/s and a deeper target depth of 3014 m (geothermie.de). The boreholes in Poing reach the top basement whilst the boreholes in AFK only reach the Malm units. Whilst no seismicity has been associated with AFK, yet seismic events of up to Ml 2.1 have been recorded in Poing close to the reinjection target in December 2016 (Megies and Wassermann 2017).
From our model results, we extract the average stress state and the standard deviation along two depth profiles at the two reinjection sites, and discuss the context of  (Fig. 9), but is relatively stable at AFK (Fig. 10). Yet, in light of the large uncertainties, this is not a significant finding. However, significantly different cohesions and frictions are widely reported for the Malm units (Hedtmann and Alber 2018;Hergert et al. 2015;Tondera et al. 2013). Thus, at the AFK site we use the mean cohesion and friction for a saturated Malm dolomite (C = 13 MPa, µ = 1.4) or limestone (C = 11.9 MPa, µ = 1.3) reported by Hedtmann and Alber (2018). These failure envelopes indicate a stable rock mass below AFK for the average stress state, even considering the uncertainties (Fig. 10). Slip tendency is a value between 0 and 1 which indicates the stability of a fault of a given orientation in a stress field (Morris et al. 1996). Slightly adapted from its original form, it is computed by Fig. 8 The average regime stress ratio (RSR) (top) and its standard deviation (bottom) on the profile A-A′ with a vertical exaggeration of 1:2. The southward dipping Malm units are clearly delineated by an increase in standard deviation. The RSR is a continuous scale from 0 to 3 to discuss the stress regime (Simpson 1997). A value of 0.5 is pure normal faulting (NF), 1.5 is pure strike-slip faulting (SS), and 2.5 is pure thrust faulting (TF) with the maximum shear stress on the fault τ max , the cohesion C, the fault normal stress σ n , and the friction coefficient µ . A value of ST = 0 indicates a lithostatic stress state with no possibility to reactivate a fault. A value of ST = 1 indicates that a fault is critical and about to fail. The stability of the two compared geothermal sites is indicated by a value of ST = 0.99 ± 0.48 in Poing and ST = 0.36 ± 0.21 in AFK.
Thus, the modelled stress state in combination with reasonable failure criteria is in agreement with the observed occurrence of induced seismicity. This is an encouraging result which shows the value of stress models for an induced seismic hazard assessment. Even with high uncertainties, beneficial information is provided. This indicates the importance of further research directed towards the quantification of additional uncertainties introduced by the variability of input parameters such as material properties and subsurface geometry.

ST =
τ max − C σ n µ −1 Fig. 9 Modelled stress state at the geothermal project Poing where seismicity up to Ml 2.1 has been recorded. The left panel shows a depth vs. stress plot with S Hmax (red), S v (green), and S hmin (blue). The average values are the bold lines, a standard deviation of one sigma (strongly shaded) and two sigma (lightly shaded, approximately 95% probability of the stress state to be within these bounds) is shown for S Hmax and S hmin . The target depth at − 2503 m below sealevel is indicated by the bold black horizontal line. The right panel shows a Mohr-Coulomb diagram at target depth. The differential Mohr circle ( S 1 -S 3 ) is indicated in blue and S 1 -S 2 and S 2 -S 3 are indicated in grey. The averages are indicated by bold lines, the standard deviation of one sigma by a strong shaded area and of two sigma by a lightly shaded area. A standard failure envelope (black) is included for reference with a friction coefficient of µ = 0.6 and a cohesion of C = 10 MPa. Failure envelopes for Malm dolomite (orange, C = 13 MPa, µ = 1.4) and limestone (blue, C = 11.9 MPa, µ = 1.3) reported by Hedtmann and Alber (2018) are also shown

Fig. 10
Modelled stress state at the geothermal project Aschheim/Feldkirchen/Kirchheim where no seismicity has been recorded yet. The target depth at − 2116 m below sealevel is indicated by the bold black horizontal line. For a detailed explanation refer to Fig. 9 Simultaneously, the possibilities of reduction of uncertainties introduced by stress data records should be investigated.

Conclusion
The presented work aims at the quantification of uncertainties in 3D geomechanicalnumerical models for the initial stress state demonstrated on a high-resolution 3D model of the Bavarian Molasse Basin. We estimate the uncertainties in the modelled stress state due to variabilities in the stress magnitude data records used for calibration. Large uncertainties in the modelled S Hmax magnitude and significantly smaller uncertainties in the S hmin magnitude are the result. The model results are supported by the observed seismicity at the geothermal project Poing, which has a critical stress state in the model. In contrast to this, the seismically quiet project Aschheim/Feldkirchen/Kirchheim at 4 km distance to Poing is uncritical in the model both for the average stress state and considering its standard deviation of 2σ . This shows the value of the model for induced seismic hazard assessment. Since the uncertainties are still large, the model encourages further research towards strategies for quantifying and reducing uncertainties of 3D geomechanical-numerical models. This includes further uncertainties due to the variability of (inhomogeneous) material properties or the geological model. Still, an average stress approach is preferred to a best-fit approach since it provides and average stress state and uncertainties instead of only a single stress state. The availability of uncertainties allows to make more transparent, significant, and reliable statements on the initial stress state.