A multi-stage 3-D stress field modelling approach exemplified in the Bavarian Molasse Basin

The knowledge of the contemporary in situ stress state is a key issue for safe and sustainable subsurface engineering. However, information on the orientation and magnitudes of the stress state is limited and often not available for the areas of interest. Therefore 3-D geomechanical– numerical modelling is used to estimate the in situ stress state and the distance of faults from failure for application in subsurface engineering. The main challenge in this approach is to bridge the gap in scale between the widely scattered data used for calibration of the model and the high resolution in the target area required for the application. We present a multi-stage 3-D geomechanical–numerical approach which provides a state-of-the-art model of the stress field for a reservoir-scale area from widely scattered data records. Therefore, we first use a large-scale regional model which is calibrated by available stress data and provides the full 3-D stress tensor at discrete points in the entire model volume. The modelled stress state is used subsequently for the calibration of a smaller-scale model located within the large-scale model in an area without any observed stress data records. We exemplify this approach with two-stages for the area around Munich in the German Molasse Basin. As an example of application, we estimate the scalar values for slip tendency and fracture potential from the model results as measures for the criticality of fault reactivation in the reservoir-scale model. The modelling results show that variations due to uncertainties in the input data are mainly introduced by the uncertain material properties and missing SHmax magnitude estimates needed for a more reliable model calibration. This leads to the conclusion that at this stage the model’s reliability depends only on the amount and quality of available stress information rather than on the modelling technique itself or on local details of the model geometry. Any improvements in modelling and increases in model reliability can only be achieved using more high-quality data for calibration.


Introduction
The contemporary in situ upper crustal stress field is of key importance for our understanding of geodynamic processes such as natural and induced seismicity (Häring et al., 2008;Gaucher et al., 2015;Scholz, 2002;Heidbach and Ben-Avraham, 2007;Townend and Zoback, 2004;Zang et al., 2014).The stress field also provides critical a priori information for safe and sustainable underground engineering such as wellbore planning and stability, reservoir management, tunnelling, mining, and underground waste storage (Altmann et al., 2014;Cornet et al., 1997;Fuchs and Müller, 2001;Moeck and Backers, 2011;Tingay et al., 2008;Zang et al., 2013;Ziegler et al., 2015;Zoback, 2010).The quantification of the criticality of the in situ stress state in terms of fault reactivation in advance of any underground treatment is essential for identifying areas of low criticality for safe and efficient utilization of the subsurface (Hornbach et al., 2015;Zoback et al., 1985;Häring et al., 2008;Kohl and Mégel, 2007).In particular, the enhancement of permeability through hydraulic fracturing should be achieved without reactivation of sealing faults or inducing seismic events of economic concern (Deichmann and Ernst, 2009;Yoon et al., 2015;Zoback et al., 1985;Townend and Zoback, 2000).
The main focus of current research is to quantify stress changes due to anthropogenic underground usage (McClure and Horne, 2014;Jeanne et al., 2014;Orlecka-Sikora, 2010;Gaucher et al., 2015;Magri et al., 2013).Induced changes of the 3-D stress state in geo-reservoirs are simulated with thermo-hydro-mechanical (THM) models since the treatment of the underground, e.g. the rate of injected fluid or the amount of mass removal, is well known (Kohl and Mégel, 2007;Gaucher et al., 2015;Van Wees et al., 2014;Jeanne et al., 2014;Cacace et al., 2013;Rutqvist et al., 2013;Magri et al., 2013).However, to assess whether the subsurface engineering pushes the system into a critical stress state in terms of absolute values, knowledge of the contemporary in situ stress, i.e. the undisturbed stress state, is essential (Hergert et al., 2015;Häring et al., 2008).
The 3-D in situ stress state can be described with a symmetric tensor of second degree with six independent components (Jaeger et al., 2007;Zang and Stephansson, 2010).Assuming that the vertical stress S v is one of the principal stresses in the upper crust, the number of independent unknowns reduces to four (Zoback, 2010).In the principal axis system these are the orientation of one of the two principal horizontal stresses, i.e. the maximum and minimum horizontal stress, S H max and S h min , as well as the magnitudes S v , S H max and S h min (Zoback, 2010;Schmitt et al., 2012).Thus, the orientation of this so-called reduced-stress tensor is described by the S H max orientation, which is systematically compiled by the World Stress Map (WSM) project (Heidbach et al., 2010(Heidbach et al., , 2008;;Sperner et al., 2003;Zoback, 1992).
Figure 1 shows a stress map with a typical density of S H max orientation data records with 172 reliable data records for the 82 000 km 2 large part of the Alpine Foreland Molasse (Reiter et al., 2015;Reinecker et al., 2010;Heidbach and Reinecker, 2013).This results in an average data density of 0.21 data records per 100 km 2 , which is the typical claim size for exploration.In general, the orientation of the stress field does not change with depth in the upper crust (Rajabi et al., 2016;Pierdominici and Heidbach, 2012;Heidbach et al., 2007).Laterally the stress field in the Alpine Molasse rotates only gently anticlockwise from east to west (Reinecker et al., 2010).Thus, the available stress orientation data allows the determination of the orientation of the reduced-stress tensor to a relatively high degree of reliability (Heidbach et al., 2007;Ziegler et al., 2016;Reiter et al., 2015).
More important for assessing criticality is the estimation of the differential stress between the magnitudes of the largest and smallest principal stresses and their changes during stimulation and production.The S v magnitude can be derived from the vertical-density distribution.In contrast to this, the horizontal stress magnitudes originate from geological history and ongoing tectonic evolution and cannot be determined directly from rock properties (Brown and Hoek, 1978;Zang et al., 2012;Zang and Stephansson, 2010).Furthermore, the increase of horizontal stress magnitude with depth is often described with a linear gradient, which is only justified when rock strength and density do not change significantly with depth (Brudy et al., 1997;Lund and Zoback, 1999).In sedimentary basins this linear increase cannot always be assumed.Competent layers, e.g. from the Malm and Muschelkalk, alternate with weaker layers with high clay content such as Dogger and Keuper and result in a sudden deviation of the stresses from a linear trend across these layers (Warpinski, 1989;Hergert et al., 2015;Cornet and Röckel, 2012;Gunzburger and Cornet, 2007;Zang et al., 2012).Moreover, the density of stress magnitude data records is, in general, up to 2 magnitudes lower than that of the orientation data (Fig. 1).
To summarize, our knowledge of the 3-D in situ stress state is based on sparsely distributed and incomplete information.Only the orientation of the reduced-stress tensor and, to a lesser extent, information about the stress regime are relatively well estimated from stress indicators.The crustal in situ stress magnitudes are underdetermined, since they vary laterally and vertically.To determine the full stress tensor for every point in a volume, a 3-D geomechanical-numerical model workflow that uses the available stress information for model-independent constraints for calibration is essential.Moreover, at reservoir scale, often no stress information is available for model calibration (Fig. 1).Thus, it is necessary to enlarge the model area until sufficient stress data are within the model volume.In the Bavarian Molasse Basin, which we use as an example, this enlargement of the model area leads to an increase in model size from a 10 × 10 km 2 reservoir-sized model to 70 × 70 km 2 regional model (Fig. 1).Considering a constant resolution, this enlargement would lead to a higher number of model degrees of freedom by a factor of 50.In most cases of THM reservoir modelling, this is beyond feasibility due to the time required for iterations and limitations in computation power.One option for avoiding a high degree of freedom is to refine the structure only in the area of interest (Jeanne et al., 2014;Westerhaus et al., 2008) (Fig. 2a).However, this becomes challenging when local structures have to be integrated.An alternative option is to use nested modelling, which is applied in various scientific disciplines such as meteorology, climate simulations, and the simulation of seismic cycles (Warner and Hsu, 2000;Cacas et al., 2001;Giorgi et al., 1998;Hergert and Heidbach, 2011).Essentially, a nested modelling approach can be (1) a local high-resolution grid inside a coarse grid where the variables are matched at the boundaries (Fig. 2b) (Oey and Chen, 1992) or (2) a multi-stage approach of two or more individual models which increase their resolution within the same area or subarea (Warner and Hsu, 2000) (Fig. 2c).In contrast to the previously named nested  (Heidbach et al., 2008(Heidbach et al., , 2010) ) and additional data from Reiter et al. (2015) and Heidbach and Reinecker (2013).Lines show the S H max orientation with line length proportional to WSM data quality (Heidbach et al., 2010).Colour coding of the data shows the stress regime with red for normal faulting (NF), green for strike-slip (SS), blue for thrust faulting (TF), and black if the regime is unknown (U).The star marks the location of the Sauerlach project where information on the S h min magnitude is available (Seithel et al., 2015).The orange box shows the lateral boundaries of the 3-D geomechanical-numerical model area (70 × 70 km 2 ) and the small black box indicates the typical dimensions of a reservoir model (10 × 10 km 2 ).The cross section A-A' (based on Przybycin, 2015) shows a 1 : 2.5 exaggeration of the area with red lines being the borehole sections and stress indicators within the model area.A complete remesh and re-evaluation is required in case of geometry or input data changes.(c) A multi-stage approach has the easiest mesh generation since the differently sized models are generated independently.Furthermore, several reservoir models can be based on the same regional model.approaches, the multi-stage procedure is most favourable in terms of required workload (fast and simple mesh generation) and quality of results (high spatial resolution in the area of interest).Furthermore, it may serve several individual reservoir model locations within the regional-area model volume (Fig. 2c).However, so far this procedure has not been applied in 3-D geomechanical-numerical modelling of the crustal stress field.
In this paper we demonstrate the applicability of the multistage nesting workflow for the 3-D geomechanical modelling of the stress tensor.We exemplify our approach with a 3-D model of the Greater Munich area in the northern Alpine Molasse Basin and a generic reservoir model (Fig. 1).We demonstrate the conceptual advantages of the multi-stage approach as a detailed, yet fast workflow for exploration from planning to exploitation.Furthermore, we quantify the impact of the uncertainties of the model parameters and the limited amount of calibration data on the model results and discuss the reliability of the 3-D geomechanical-numerical modelling.

Geological setting
The northern Alpine Molasse Basin is a typical asymmetric foreland basin which extends over 1000 km from Lake Geneva in the west to Lower Austria in the east (Bachmann et al., 1987).Its widest N-S extent is 130 km in southern Germany (Lemcke, 1988).The basin mainly consists of Tertiary sediments on top of Mesozoic successions and a Variscan basement with Permo-Carboniferous troughs (Lemcke, 1988;Bachmann et al., 1987).In the Foreland Molasse these sediments dip towards the south where a maximum thickness of approximately 6000 m is reached in front of and beneath the Alpine mountain chain and the Folded Molasse (Fig. 1) (Bachmann et al., 1987).Due to the Molasse Basin's close link to the Alpine orogeny (Schmid et al., 2008) most of the main faults in the Bavarian Foreland Molasse are steeply dipping (> 60 • ) and strike at least subparallel to the Alpine front approximately E-W (Reinecker et al., 2010;Bachmann et al., 1987;Lemcke, 1988).They are considered mostly inactive at the moment (Reinecker et al., 2010;Bachmann et al., 1987;Lemcke, 1988).
For our model geometry we use the upper 9 km of the 3-D structural model of the northern Alpine Foreland Basin by Przybycin (2015), which covers the entire German part of the Molasse Basin.It provides 12 stratigraphical units in total with a focus on the Malm and Purbeck, two target horizons for geothermal exploration (Lemcke, 1988;Bachmann et al., 1982;Fritzer et al., 2012).The lateral resolution of the structural model (1 × 1.7 km 2 ) is sufficient to provide the geometry for the generation of our regional-scale model of the Greater Munich area.The structure is based on freely available data on the depth and thicknesses of stratigraphic units from wells and seismic lines as well as 3-D gravity mod-elling as a further constraint (Przybycin, 2015).The part of the structural model used for the geomechanical model has a size of 70 × 70 km 2 and is referred to as the root model.It includes the sediments in the Molasse Basin in their entire vertical extent.The bottom of the model is situated at a depth of 9 km entirely within the upper crust.The generic reservoir model located within the root model volume is called a branch model.It has a size of 10 × 10 km 2 with more detailed structural information, e.g.provided by a 3-D seismic survey.
3 In situ stress data

Orientation of S H max
Within the root model area (Fig. 1, orange box) 18 reliable S H max orientation data records are located, while there are none in the branch model area (Fig. 1, black box).These data are exclusively from borehole measurements using drillinginduced tensile fractures (Aadnoy, 1990) and borehole breakouts (Bell and Gough, 1979;Bell, 1996) as indicators for the S H max orientation (Reinecker et al., 2010).In 15 wells in the model area, borehole breakouts are found with a combined length of 7.7 km.In 3 wells drilling-induced fractures are found with a combined length of 0.3 km.The stress indicators are found mainly between the surface and a depth of 2-3 km even though some are located at greater depth (Fig. 1).No significant stress rotation or perturbation with depth is observed in the available data (Reinecker et al., 2010;Heidbach et al., 2016).The quality of the data is exceptionally good according to the WSM quality ranking (Heidbach et al., 2010;Sperner et al., 2003;Zoback, 1992) with eight A-quality data records (i.e. an uncertainty of ±15 • ), six B-quality data records (±15-20 • ), and four C-quality data records (±20-25 • ).Under the assumption that S v is a principal stress component, the reduced 3-D stress tensor within the model area has a mean S H max orientation of 1.7 • ± 19.2 • which is approximately perpendicular to the Alpine front (Fig. 1).

Stress magnitudes
The magnitude of S v can be estimated with a relatively high reliability from the thickness of the different overlying units (z) in the structural model provided by Przybycin (2015), the density of the corresponding rock material (ρ rock , Table 1) and the gravitational acceleration (g) given by S v = σ zz = ρ rock gz.
(1) However, information on the horizontal stress magnitudes is sparse even within the root model area.The magnitude of S h min is usually derived from hydraulic fracturing (Haimson and Fairhurst, 1969;Hubbert and Willis, 1972), but such data are not available publicly for the Bavarian Molasse Basin.Alternatively, leak-off tests (LOTs), which rely on a cheaper  Lama and Vutukuri (1978), e Koch and Clauser (2006) and faster method, are more frequently used for the estimation of S h min .They provide information on the break-down pressure of the tested formation, which is then used as an approximation for the S h min magnitude (Haimson and Fairhurst, 1969;Bell, 1990;Zang et al., 2012).Further information on the S h min magnitude can be derived from a formation integrity test (FIT).It does not fracture the rock but provides a minimum pressure value at which the rock is stable, which in turn provides a lower bound for the S h min magnitude (Zoback et al., 2003).Even though no hydraulic fracturing was done in the model area a LOT has been conducted in the Unterhaching Gt 1/1a borehole which is used for calibration (T.Fritzer, personal communication, 2016).Furthermore, several FITs have been performed in the borehole Sauerlach (Fig. 1) that is in the root model area (Seithel et al., 2015).In contrast to the LOTs FITs are not used for calibration since the difference between the FIT pressure and the actual magnitude of S h min is not known.However, during one of the FITs in the Sauerlach borehole bore fluid was lost into the formation (T.Fritzer, personal communication, 2016).Hence a leak-off occurred and this FIT is treated as a LOT and used for the model calibration.
The direct estimation of the S H max magnitude would only be possible with overcoring measurements (Hast, 1969;Sjöberg et al., 2003).In addition, reasonable values for the S H max magnitude can be derived on the basis of the frictional equilibrium theory (Zoback et al., 2003) or physics-based relations for which reliability is largely dependent on the quality of the S h min magnitude estimation (Zoback, 2010;Cornet, 2015).Seithel et al. (2015) use the friction equilibrium approach and derive a single S H max magnitude between 112 and 116 MPa at a depth of 4 km.We use an S H max magnitude of 112 MPa in Sauerlach for the calibration even though the uncertainties introduced by the derivation are high.The impact of these high uncertainties on the model results is discussed later on.

Stress regime
In areas with a low number of magnitude stress data records, the stress regime provides information on the relative magnitudes of S v , S H max , and S h min .The stress regime is mainly derived from focal mechanisms of seismic events and, to a small extent, from geological indicators or hydraulic fracturing experiments (Zoback, 1992;Sperner et al., 2003).In the Swiss part of the northern Alpine Molasse Basin a strike-slip (S H max > S v > S h min ) and, to a smaller extent, extensional (S v > S H max > S h min ) stress regime is mainly observed (Heidbach and Reinecker, 2013).However, in the Bavarian Molasse Basin north of the Alpine front, no natural seismicity has been recorded (Grünthal, 2011;Grünthal and Wahlström, 2012) to derive the stress regime from focal mechanisms.
Information from structural geology observing steeply dipping faults in the Bavarian Molasse Basin (Bachmann et al., 1987;Lemcke, 1988) indicates an extensional tectonic faulting regime (Anderson, 1905(Anderson, , 1951)).In contrast to this Illies and Greiner (1978); Lemcke (1988), andReinecker et al. (2010) propose a compressional (S H max > S h min > S v ) or strike-slip stress regime in the Molasse Basin.Seithel et al. (2015) also propose a strike-slip stress regime at a depth of 4 km for the Sauerlach project according to their analysis based on the frictional equilibrium theory.However, without further estimations of the stress magnitudes in other depth sections and locations, the regional tectonic stress regime setting is subject to large uncertainties. of the equilibrium of forces.Furthermore, we assume a linear elastic rheology and solve for absolute stresses (no pore pressure).The general model procedure follows the technical workflow explained in detail by Hergert et al. (2015) and Reiter and Heidbach (2014).
The root model extends 70 × 70 × 10 km 3 in the entire Greater Munich area (Fig. 1).It consists of six different stratigraphic layers (Table 1) based on the 3-D structural model by Przybycin (2015).The generic branch model of a potential geothermal site has a size of 10 × 10 × 10 km 3 and includes six different stratigraphic units (Table 1).For both models the boundaries are oriented perpendicular and parallel to the orientation of S H max and S h min respectively (Fig. 1).Both models are populated with the Young's modulus, the Poisson ratio and the density according to the stratigraphic units (Table 1).
An exact fit of the overburden S v is achieved by applying gravity, provided that the density of the stratigraphic units is correctly chosen.We implement an equilibrated initial stress state close to lithostatic (S H max ≈ S h min ≈ S v ).Dirichlet boundary conditions (i.e.displacements) are applied to the sidewalls of the model to create horizontal differential stresses.The boundary conditions are adjusted in a way that the modelled magnitude of S H max and S h min at the calibration points fit the observed magnitudes.
Due to the complex topology of the stratigraphy and inhomogeneous rock properties of the different units the finite element method (FEM) that allows unstructured meshes is used to solve the partial differential equation of the equilibrium of forces at discrete points.Thus, both models are discretized into finite element meshes.The root model is constructed with approximately 10 6 hexahedral elements resulting in approximately 400 m of horizontal and between 15 and 700 m of vertical resolution (Fig. 3).A vertically refined resolution is created in the units of interest for geothermal exploration, namely the Malm and Purbeck formation.The Cretaceous and the Triassic (pre-Malm) are only preserved in parts of the root model.Compared to the root model a significantly finer resolution with a total of 21 × 10 6 tetrahedral elements is chosen in the branch model.The edge length of the elements varies between 10 and 160 m with the coarsest resolution located at the bottom and the edges of the model and the highest resolution in the Purbeck and Malm units of interest for geothermal exploration (Fig. 3).

Model calibration procedure and two-stage approach
The calibration of the root model with stress magnitude data is achieved by applying two Dirichlet boundary conditions, each on one of the perpendicular sidewalls of the model (Fig. 4, left row).A single S h min magnitude data record can be exactly modelled by a certain combination of two boundary conditions.More precisely an unlimited combination of two boundary conditions exist to achieve an exact fit of a single S h min magnitude calibration point.This unlimited number of combinations of displacement boundary conditions is a linear function of the E-W and N-S displacements and is displayed as a linear slope in Fig. 4a with displacement in N-S direction on the x axis and displacement in E-W direction on the y axis.Due to the assumed linear elastic model rheology, each combination of east-west and north-south displacement that lies on the slope leads to an exactly calibrated model (Fig. 4a).
If several S h min magnitudes are available for calibration, each of them can be exactly reproduced by an unlimited number of combinations of displacement boundary conditions.However, to achieve a calibration which works for all of the observed S h min magnitudes, a single "best-fit" slope is derived from the linear slopes for the individual calibration points using a linear regression (Fig. 4b) (Reiter and Heidbach, 2014).Each combination of displacement boundary conditions specified by this slope results in a best-fit model for all of the considered calibration points.
The same procedure is applied for the calibration of S H max magnitudes so that eventually a best-fit slope for both the S H max and S h min magnitude stress data records used for calibration are available (Fig. 4c).Displacement boundary conditions defined by the point where these two best-fit slopes intersect are used to compute the best-fit model that reproduces the S H max and S h min stress data records best (Fig. 4c) (Reiter and Heidbach, 2014).
Application of this calibration procedure is fast and simple since the best-fit boundary conditions can be found by combining two linear slopes based on the calibration data and the displacement boundary conditions.Therefore, to find the best-fit boundary conditions only three different models with arbitrary displacement boundary conditions are required (Fig. 5a).The modelled S H max and S h min magnitudes at the location of calibration points in each of the three models are compared to the actually observed data records (Fig. 5b, c).A linear regression with two unknown variables leads to the best-fit slopes for the combination of boundary conditions to model the S H max and S h min magnitude (Fig. 5d).At the in-tersection of the two slopes, the boundary conditions for the best-fit model can be derived (Fig. 5d).
It is assumed that the stress data records used for the calibration are the result of the far-field stress state and its interaction with structural features such as local density and/or strength contrasts represented within the root model.If the measurements were, e.g. the result of an unknown or unimplemented local active fault, the results of the calibration would not be reliable.Thus, in general, the data used for calibration should be representative for a large volume of the individual lithological layer.
Under this assumption the best-fit model simulates the stress state at discrete points in the entire model volume.Hence, information on the stress state is now also available in areas of the root model where previously no observables on the stress state were available.This means that in the branch model, which does not include any observed stress data records, simulated information on the stress state is also now available from the root model and can be used to calibrate the branch model (Fig. 5d, f).
Since the branch model is calibrated in the same way as the root model (but with a simulated stress state from the root model as calibration points instead of observed stress data records), a large number of potential calibration points with S h min and S H max magnitudes are available.The S H max and S h min magnitudes at each calibration point can be modelled individually in the branch model by combinations of boundary conditions, each described by a linear slope (Fig. 4a).For all S H max and S h min magnitudes a best-fit slope is derived, based on the individual linear slopes (Fig. 4b).Two best-fit slopes describe the combinations of boundary conditions which model the S H max and S h min magnitudes best.The intersection of these two best-fit slopes defines the boundary conditions, which are used to compute the best-fit branch model (Fig. 4c).This calibration procedure is performed analogously to that of the root model (Fig. 5e-h).
For a successful transfer of the stress state from the root to the branch model, it is critical that the stress state used for the calibration of the branch model is obtained at discrete points of the root model and not in its volume.Otherwise the stress state extracted from the root model is potentially biased due to interpolations from discrete points into the volume, which are performed by the visualization software.Since the large number of possible calibration points can be chosen arbitrarily, their locations need to be considered.We recommend using calibration points close to the border of the branch model but outside the zone prone to boundary effects.Calibration points from the root model in the centre of the branch model are a contradiction of the two-stage approach which aims at finding local stress changes due to high-resolution structural features that are only present in the branch model.Due to the lack of any other stress data in the branch model area, the calibration procedure imposes the root model's basic stress state on the branch model, which prevents such local stress perturbations.Hence, this necessary imposition should be re- duced to the boundaries of the branch model that are not used for interpretation.Furthermore, the calibration points should be evenly distributed along the branch model boundary and represented in all stratigraphic units to account for different material properties.Special attention needs to be paid to units which are either only present in the root or the branch model or have a significantly different geometry or rock properties in the two models.

Model results
In the following two sections we present the results of two model scenarios for the root model that fit equally well the observed stress data, but with different stress regimes (Fig. 6).For the branch model we present our results on one scenario that can be considered our best-fit model (Fig. 7).

Root model
The best-fit root model of the stress state at discrete points in the Greater Munich area is calibrated using S h min magnitudes from the two LOTs and one S H max magnitude described in detail in the stress data in Sect.3.2.The best-fit model has a good fit to the two S h min calibration data points and an almost perfect fit for the single S H max calibration point.Deviations between observed and modelled data are on average 0.4 MPa for the two S h min calibration points and 0.04 MPa for the single S H max calibration point.
Figure 8 shows the best-fit model results along the Sauerlach borehole profile along with the FIT data of Seithel et al. (2015).The black line shows the vertical stress magnitude with depth, which depends only on the chosen rock density.The blue line is the S h min magnitude, which is larger than all FIT values at all depth sections.The blue star represents the magnitude and depth of the S h min magnitude inferred from a FIT with leak-off.The red line is the S H max magnitude in the best-fit model while the dashed line represents S H max in another model scenario.The red star marks the depth and magnitude of S H max in the best-fit model.The shaded areas show the modelled magnitudes for model scenarios, which use S H max magnitudes between 92 and 118 MPa in a depth of 4 km below the Sauerlach site for calibration.This demonstrates that the single S H max magnitude derived in conjunction with the ambiguity of the stress regime opens up a wide range of model scenarios which all equally well fit the S h min data.Even though a compressional regime can be excluded by the available data in Sauerlach, no indication exists of whether S H max > S v or S H max < S v .That means that the prevalence of a normal faulting or a strike-slip stress regime is possible.To account for this variability, several different scenarios have been computed, two of which are compared in Fig. 6.Note that the only difference between these scenarios is the S H max magnitude value used for the root model calibration (Fig. 6a 96 MPa, Fig. 6b 112 MPa); the fit to the S h min data from the LOTs is equally good (Fig. 6).
In Fig. 6 we show a number of scalar stress values derived from the modelled 3-D stress tensor on cross sections and within stratigraphic units for the aforementioned two model scenarios.The figure shows that the values vary depending on the stratigraphic units horizontally and laterally.More importantly, the results from the two model scenarios which fit the model-independent calibration data equally well are quite different.The first row of Fig. 6 shows the variability of the stress regime using a continuous scale, the so-called regime stress ratio (RSR) from Simpson (1997).Close to the surface a strike-slip regime dominates with compressional components in some areas.With increasing depths this changes to a prevailing extensional regime.Moreover, some changes from strike-slip to extensional and back to strike-slip can be observed.They are not a smooth linear trend but are highly dependent on the lithology.
The second row of Fig. 6 displays the horizontal stress anisotropy as a stress magnitude ratio of S H max /S h min on a N-S and E-W cross section through the two model scenarios of the root model.It is clearly visible that the ratio varies significantly with depth and between the model scenarios.The southward-dipping Malm and Purbeck units have stress ratios of up to 0.15 higher than the basement layer and overlying sediments respectively.
The last row in Fig. 6 shows the differential stress in the middle of the Malm unit.Both model scenarios show higher differential stresses in the south where the Malm units are deeper than in the north.The largest N-S difference is 7 MPa in model scenario (a) in contrast to 12 MPa in model scenario (b), even though the relative pattern of the differential stresses in the Malm unit is quite similar in both model scenarios.This pattern highlights the main trend of an increasing differential stress towards the south.At the same time significant changes of the differential stress within less than 10 km of up to 1 MPa are predicted.

Branch model
In this section we show the model results of the branch model (Fig. 7) that uses the stress data derived from the root model scenario displayed in the right row of Fig. 6.In order to visualize the criticality of the reservoir, we use two scalar values which are computed from the modelled 3-D stress state.The first one is the fracture potential (FP) of intact rock volume (Connolly and Cosgrove, 1999).It is computed as and an alternative scenario that fits the S h min values equally well, but is calibrated against a lower S H max value (a).The overall difference that results from the different S H max values used for the calibration is expressed in the continuous scale of the regime stress ratio (RSR), which is between 0.5 (normal faulting regime), 1.5 (strike-slip), and 2.5 (thrust faulting regime) (Simpson, 1997) in the model volume.The horizontal stress anisotropy expressed in the ratio of S H max /S h min is shown on two cross sections which intersect below Munich.The differential stress (difference between the maximum and minimum principal stress, lowermost panel) is mapped on a surface which is vertically centred in the Malm α − γ units.(Morris et al., 1996).It is a measure of the criticality of faults which is illustrated as a scalar value for the distance to failure derived from the stress tensor with values between 0 (safe) and 1 (failure).Slip tendency is computed for faults or fault segments of a certain orientation and is defined as with the maximum shear stress τ max , the normal stress σ n , the friction angle = arctan(µ), and the friction coefficient µ.The application of these two values is shown in the branch model with generic faults in Fig. 7.
The high dependence of slip tendency on the orientation, friction, and cohesion of the fault is displayed in Fig. 7.A high variability of slip tendency between 0.05 and 0.3 is observed on the generic faults.This variability is induced by the 3-D stress tensor and the curved fault surfaces.Furthermore, due to differently assumed friction and cohesion of the rocks, the Malm δ -Purbeck units have a clearly smaller value of slip tendency compared to the Chattian units in the hanging wall and the Malm α − γ in the footwall.The fracture potential in the basement generally lies between 0.1 and 0.2, which is quite low, hence it requires high pressure for hy- draulic fracturing operations to enhance the permeability of the fracture network.
Information provided by the branch model is used in an early pre-drilling stage of a project to assess whether the initial conditions of the reservoir and its criticality allow safe production; i.e. both slip tendency and fracture potential have low values as in Fig. 7. Before the drilling of the borehole begins the planning of the drill paths can be optimized.Especially if intersections with faults are required, their paths can be monitored and adapted in a way that they circumnavigate fault segments which have a higher value of slip tendency, meaning that this fault segment is more favourably oriented for a potential failure compared to other fault segments.In Fig. 7 areas with cool colours are preferred for intersections of boreholes with faults compared to areas with hot colours.In Fig. 7 the Malm δ -Purbeck unit is mostly blueish coloured which indicates the lowest slip tendency values.Hence these are the best units for the intersection of boreholes with faults.An intersection with the northernmost fault in the red areas should be avoided.Lines show the results of the best-fit root model: blue for the S h min magnitude, black for the vertical stress S v , and red for the S H max magnitude.The blue dots are formation integrity tests (FITs), which are a lower boundary for the magnitude of S h min and not used for calibration, the blue star represents the suspected LOT, the red star shows the S H max magnitude of 112 MPa used for calibration (Seithel et al., 2015).Shaded areas in the same colour around the lines show the range of model scenarios that fit equally well to the model-independent constraints.The dotted red line shows the S H max magnitude for the model scenario in Fig. 6a.

Reliability of the model results
One of the key points in geomechanical modelling is the reliability of the model results in terms of the predicted processes and the presented multi-stage simulation of the in situ stress field.As already mentioned in the result Sect. 5 the calibration procedure introduces uncertainties due to the low number of data points as well as their relatively large uncertainties.Further uncertainties are introduced by the model input, e.g.calibration data, rock properties, and structure.Hence, the reliability of the model depends on the uncertainties of the input data used for the model.To quantify the model's reliability we use the already presented scalar value slip tendency (Morris et al., 1996), for which variability is introduced by the uncertainties in different input data.
We compute the slip tendency for model scenarios which use the extreme values of the input parameters range of uncertainties.The model's linear elastic behaviour allows the individual quantification of the impact of different model pawww.solid-earth.net/7/1365/2016/Solid Earth, 7, 1365-1382, 2016 rameter uncertainties on the model's reliability.Therefore we compute several model scenarios in which sequentially only a single parameter is changed to an extreme value.This enables us to derive the individual impact of different parameters and quantify the most important ones.The results of the slip tendency for each model scenario are subsequently compared to the best-fit slip tendency values from the best-fit model (Table 2).The variations of slip tendency introduced by the different independent parameters are added together, which leads to an expected maximum variability in slip tendency of ±0.57.
The two main sources for the variability of slip tendency can be identified as the model-independent data for the S H max magnitude used for the model calibration and the rock properties density, Young's modulus, and Poisson ratio (Table 2).A high variability of slip tendency of 0.14 is introduced by uncertainties in the S H max magnitude.Since the S H max magnitude is derived under several assumptions, a wide range of possible S H max magnitudes is used for the calibration of the slip tendency model scenarios.Due to the fact that only limited knowledge and measurements of the rock properties are available a wide range of values are possible and they introduce a high variability of 0.18 in slip tendency.
Slip tendency proves to be quite robust (±0.01) to the small uncertainties in the S h min magnitude under the assumption that the available data used for calibration is a valid proxy for the entire model (Table 2).Likewise only small variations in slip tendency are introduced by changes of ±10 • in the fault strike (±0.02) and dip (±0.03).The cohesion and friction angle act as more sensitive parameters (each ±0.07).Finally the two-stage calibration procedure itself introduces some moderate deviations (±0.05) with a large number of calibration points and their individual locations used in the branch model.

Discussion
The objective of this work was to demonstrate the multistage approach for a high-resolution 3-D geomechanicalnumerical modelling workflow assessing the criticality in reservoirs.In contrast to a single model, which includes both stress data records for calibration and high-resolution representation of a local reservoir structure, we use two models of different sizes.The regional-scale root model is calibrated on stress data records and provides the stress state for the calibration of the reservoir-scale branch model.This approach provides a cost-efficient, quick, and reliable state-of-the-art calibration of geomechanical-numerical models of the contemporary 3-D in situ stress field across scales.It is used to assess the criticality of reservoirs which can be quantified by scalar values such as slip tendency.If detailed information on the fracture behaviour of the rock are known, more elaborate fracture criteria than Mohr-Coulomb (e.g.Sulem, 2007;Zang and Stephansson, 2010) can be applied to analyse the model results.Furthermore, the approach provides the initial stress state for local application such as in THM models.

Workflow
A single model with the same functionality as the two models in the multi-stage approach needs to account for the required high resolution in the reservoir area and the large model extent to include data for calibration.These two requirements are not contradictory per se but prolong the process of mesh generation, e.g. by needing to harmonize a regionalscale low-resolution and local-scale high-resolution structural model in the area of overlap.Furthermore, the manageability of the model (e.g.logical size) and the available time for computation (number of elements) in most instances requires a variable resolution which is refined only in the target area.Such a change in element size in a single model is possible but the mesh generation is cumbersome and needs a high number of elements.For a THM simulation of production and (re)injection, incrementation over time significantly increases the computation time for each single element.Furthermore, in a single large model, only a very small area is of interest, hence large areas are simulated to no purpose while at the same time the logical size, computation time, and effort are increased.If a multi-stage approach with two models is applied, each model has its own fixed resolution with no required variation in element size (Fig. 2c).This significantly speeds up and simplifies the process of model generation since neither the structural models need to be harmonized nor a large difference in elements size needs to be implemented.Regarding the same resolution in the target area, the time required for computation decreases, but not as much as the logical size of the models, which improves the model's manageability.A geomechanical root model can also provide the stress state for a THM branch model which helps to save computation time by focussing the time-consuming THM simulation on the actual area of interest.Calibration data records for the additionally required scalar values on the pore pressure or temperature are provided in the literature or by dedicated models, e.g.Przybycin et al. (2015).
In addition the application of two models opens further possibilities for improved and safer exploration and drilling.Structural features and stress magnitude measurements recorded during advanced exploration or even initial drillings can be implemented into the model workflow due to the simplified mesh regeneration.Even a change in the target area within the root model can be more easily implemented in the workflow since only a new branch model is required.The calibration of the root model can be updated with new stress data records as soon as they become available.Finally, a large calibrated root model may include several target areas and can be reused and applied for more than one project.

Calibration
The two models in the presented two-stage approach are calibrated with different Dirichlet boundary conditions applied to an initial stress state.This approach follows the modelling procedure using isotropic elastic materials described by e.g.Reiter and Heidbach (2014), Hergert et al. (2015), and Gunzburger and Magnenet (2014).Almost identical results can be achieved by the application of according orthotropic elastic material and gravity loading only (Cornet and Magnenet, 2016).For deep lithosphere and asthenosphere models elasto-plastic materials with the application of gravity but no further boundary conditions can be applied and yield similar results (Maury et al., 2014).
Our root model is calibrated with data records which display the stress state as a result of the geologic history and tectonic evolution.In the presented region the stress field is very homogeneous but in other regions significant local lateral variations exist and need to be accounted for.This can be accomplished for example by lateral variations of the material properties or faults.It is crucial to ensure that the data used for the calibration is representative for the regional material and geometry in the root model.
The branch model, however, is calibrated on the stress state simulated in the root model.Both calibration procedures are not limited in the number of calibration points and a weighting of the calibration points according to reliability can be easily realized.An extension of the two-stage approach to include three (or even more) models of different sizes is possible.Furthermore, the calibration procedure allows running several alternating models with different calibration data or differently weighted calibration points as well as variations in rock properties to quantify model-specific variations.This ability was used to quantify the reliability of the model's results.It is also useful for future attempts at statistically determining uncertainties in the model's results.
Even without any additional computations, a first-order assessment of the impact of individual data records on the model calibration can be made by assessing changes in the boundary conditions.Therefore the best-fit boundary conditions derived with and without certain data records are compared.Such a data record could be a newly performed hydraulic fracturing experiment which provides an additional S h min magnitude data record.The variation of the derived boundary conditions induced by such a new data record provides a first idea of the variation of the stress state.Although, this feature cannot be used as a replacement for computations it helps to identify whether the newly included calibration point yields a significantly different stress state which requires a reassessment of the situation or if the changes are minor and the exploration can be continued as planned.
The models showed in this work do not include any implicit faults and no strain partitioning is assumed.The calibration of a model including faults and fault-specific behaviour, e.g.strain weakening or hardening or long-term relaxation of the gauge material, is possible as well if sufficient information on the fault properties are available.However, due to the non-linearities introduced by active faults the calibration process requires a regression analysis of a higher degree, hence several more test scenarios.This is beyond the scope of this work.

Model independent reliability
Apparently the model's reliability is mainly affected by the lack and high uncertainty of S H max magnitude data.The large influence of the S H max magnitude is shown by two different models for viable S H max magnitudes in Fig. 6.A feasible method to narrow down the S H max variability is to enhance the knowledge of the Andersonian stress regime, e.g. by gathering information on earthquake focal mechanism data (if available) or the crack orientation induced by leak-off tests or hydraulic fracturing (Haimson and Fairhurst, 1969;Hubbert and Willis, 1972;Zang and Stephansson, 2010).Such information is most likely available in the model area but not publicly accessible.Furthermore, an array of many expensive deep overcoring measurements (several per borehole) could provide valuable information on the stress state and S H max in particular (Hast, 1969;Sjöberg et al., 2003).
The uncertainties related to the material properties are another large factor that limits the model's reliability.This can be mitigated at least partly by using data from extensive databases (e.g.Bär et al., 2015;Lama and Vutukuri, 1978;Koch, 2009) or by converting seismic velocities which are founded on empirical relations (Mavko et al., 2009).Finally, averaging mean values from several laboratory tests of rock samples from the area and lithologic formations of interest are the safest but most expensive ways to retrieve reliable information of rock properties.
The uncertainties in the strike and dip of faults have a comparably small share in the reliability of the model while be-ing challenging to mitigate due to the general uncertainties in the interpretation of 3-D subsurface structures.The fault parameters cohesion and friction angle which are even more difficult to determine compared to the orientation reduce the model's reliability to a slightly higher degree compared to strike and dip.Increasing the model's reliability through a better understanding of these parameters is possible but requires a detailed understanding of the great variability of in situ fault zone behaviour and extent at depth.
Statistical methods to quantify uncertainties in the subsurface geometry exist for purely static structural models (Wellmann, 2013;Wellmann and Regenauer-Lieb, 2012).However, the computation time required to extend this approach to a 3-D geomechanical-numerical modelling approach and the ensuing analysis is beyond the scope of this work.A further investigation should be conducted as a sequel to the work by Bond et al. (2015) in a generic approach including geomechanical-numerical modelling.

Model dependent reliability
This model focusses on the stress tensor in the uppermost part of the crust and its extent is accordingly chosen.Deepseated processes at depths larger than 9 km are, therefore, not represented in the model.However, as shown by Maury et al. (2014), the lateral variations in the differential stress in the depths are small compared to variations introduced by the uncertain material properties and magnitude of S H max in our model.Furthermore the influence of deep structures such as the Moho geometry is minor, as shown by Reiter and Heidbach (2014) or Hergert et al. (2011).
The model does not include any faults.The inclusion of faults makes sense in situations where detailed information on fault geometry, extent, and parameters are available and a significant impact of the faults on the regional stress field or (re)activation is expected.However, in this example, the available stress data suggests that no faults with a major impact are located neither within the root model nor the branch model area.Therefore the variations introduced by omitting faults is assumed to be small.
Variations of the model results are also introduced by the multi-stage calibration approach itself and cannot be mitigated due to both models 3-D stress state with lateral and vertical variations.The model's calibration, however, depends on the variations of only two independent boundary conditions.Additionally, small variations may be introduced by the model assumptions.However, these variations can be disregarded in the light of the major reasons for variations due to the small amount of stress magnitude data and rock properties.Table 2 clearly shows that any further advances in modelling are not efficient as long as the amount and quality of input data (S H max , rock properties) is not increased.

Application in geoengineering
Hydrocarbon reservoirs are currently exploited on a minor level in the Alpine Foreland (Lemcke, 1988;Sachsenhofer et al., 2006) and some of the former reservoirs are used for oil and gas storage (Sedlacek, 2009).However, hydrothermal reservoirs of economic interest for district heating or power generation are available (Lemcke, 1988;Bachmann et al., 1982;Fritzer et al., 2012).These reservoirs are situated in highly karstified limestones of the Late Jurassic which are locally referred to as Malm formations (Lemcke, 1988).As of 2016 those deep reservoirs have already been exploited by 21 municipal geothermal power plants and district heating projects of which Aschheim, Dürrnhaar, Erding, Freiham, Garching, Holzkirchen, Ismaning, Oberhaching, Poing, Riem, Sauerlach, and Unterschleißheim are in the root model area (Bundesverband Geothermie, 2016).Borehole data from these projects could be easily implemented in the calibration of the root model and would increase its reliability if they became publicly available.
Within the root model perimeter, several geothermal projects are currently at the planning stage, namely Bernried, Gräfelfing/Planegg, Königsdorf, Markt Schwaben, Puchheim/Germering, Raststätte Höhenrain, Starnberg, Weilheim/Wielenbach, and Wolfratshausen (Bundesverband Geothermie, 2016).In addition the municipal energy supplier of Munich (SWM) plans to install an extensive geothermally driven district heating grid for the entire city (Stadtwerke München GmbH, 2012).Therefore, a 3-D seismic survey was conducted in the entire southern part of Munich in winter 2015/16 (Bundesverband Geothermie, 2015).The presented root model provides data for a first-order assessment of the in situ stress state at the exact locations of these planned geothermal projects.Furthermore, it provides calibration data for local-/reservoir-scale models based on high-resolution 3-D seismic surveys which simulate the stress state, its criticality, and the possibility of subsidence due to the production and reinjection of fluid and heat.
Furthermore, the two-stage approach could be extended to a three-stage approach which incorporates a global model of the entire Bavarian Molasse Basin.More data for calibration, as well as more potential applications, might be available in such an enhanced area.Thereby the regional or global root model could be established as a community model which provides the stress state for further applications and/or local models for planned projects.

Conclusions
In this work we present a multi-stage 3-D geomechanicalnumerical modelling approach, which provides a costefficient, reliable, and fast way to generate and evaluate the criticality of the stress state in a small target area where, in general, no stress data for model calibration are available.The approach uses a large-scale root model which is calibrated on available stress data and a small-scale branch model which is calibrated on the root model.We exemplify this in a two-stage approach in the German Molasse Basin around the municipality of Munich.The discussion of reliability of the model results clearly shows (1) that variations are large and (2) that they are mainly introduced by the uncertain material properties and missing S H max magnitude data.At this stage, the model's quality depends on the amount and quality of available input data and not on the modelling technique itself.Any further improvements in the model's resolution and applied techniques will not lead to an increase in reliability.This can only be achieved by more high-quality data for calibration.

Data availability
The stress orientation data used for model set-up and calibration is available from Reiter et al. (2016) and Heidbach et al. (2016).

Figure 1 .
Figure 1.Stress map of the Bavarian Molasse with 172 A-C quality data records based on the World Stress Map database release 2008(Heidbach et al., 2008(Heidbach et al., , 2010) ) and additional data fromReiter et al. (2015) andHeidbach and Reinecker (2013).Lines show the S H max orientation with line length proportional to WSM data quality(Heidbach et al., 2010).Colour coding of the data shows the stress regime with red for normal faulting (NF), green for strike-slip (SS), blue for thrust faulting (TF), and black if the regime is unknown (U).The star marks the location of the Sauerlach project where information on the S h min magnitude is available(Seithel et al., 2015).The orange box shows the lateral boundaries of the 3-D geomechanical-numerical model area (70 × 70 km 2 ) and the small black box indicates the typical dimensions of a reservoir model (10 × 10 km 2 ).The cross section A-A' (based onPrzybycin, 2015) shows a 1 : 2.5 exaggeration of the area with red lines being the borehole sections and stress indicators within the model area.

Figure 2 .
Figure2.Different types of modelling approaches.(a) A refined mesh in the area of interest is expensive and inefficient due to the large number of elements required for the discretization of the gradient in resolution.Furthermore, it requires a complete remesh and re-evaluation in case of any change in the geometry or input data.(b) A local model nested within a regional model matches the variables at the boundary.A complete remesh and re-evaluation is required in case of geometry or input data changes.(c) A multi-stage approach has the easiest mesh generation since the differently sized models are generated independently.Furthermore, several reservoir models can be based on the same regional model.
set-upBoth the regional-scale root model and the reservoir-scale branch model are based on the same modelling assumptions.Assuming that accelerations other than gravity can be neglected, the models solve the partial differential equation www.solid-earth.net/7/1365/2016/Solid Earth, 7, 1365-1382, 2016

Figure 3 .
Figure 3.The root and branch model discretized with 10 6 hexahedral and 21 × 10 6 tetrahedral elements respectively.Please note that to improve visibility the discretization of the branch model is only displayed within the magnified inset.Both magnified regions show the Malm α − γ (turquoise) and Malm δ − ζ and Purbeck (purple) units, which are the predominant target units for geothermal exploration in the Bavarian Molasse Basin.

Figure 4 .
Figure 4. Left: exemplified schematic models with the data records used for calibration (star: S h min magnitude, circle: S H max magnitude).Right: linear slopes that display the magnitudes of possible combinations of displacement boundary conditions applied normal to the E-W (y axis) and N-S (x axis) sidewalls of the model.For each data record an individual slope defines the possible combinations of boundary conditions to fit the model to this calibration data record.(a) A single S h min magnitude can be calibrated by an unlimited number of combinations of boundary conditions which are on a linear slope.(b) Several S h min magnitudes usually cannot be calibrated to an exact fit.However, a linear regression of all the linear slopes derived for the calibration of each individual data record provides a best-fit slope.This slope defines combinations of bestfit boundary conditions that fit the data records used for calibration equally well.(c) Several S h min and S H max magnitude data records used for calibration results for each S h min and S H max in a linear slope of combinations of best-fit boundary conditions.At the intersection of these two slopes the best-fit boundary conditions (indicated by a star) are found for the calibration of S H max and S h min together.

Figure 5 .
Figure 5.The calibration workflow for the root and branch model.(a) Three models with different Dirichlet boundary conditions provide stress data comparison values for a calibration with (b) observed magnitude stress data.The deviation of the modelled from the observed stress state of each of the three scenarios (c) is used in a linear regression to derive the boundary conditions to compute the best-fit root model (d).(e) Three different branch models provide stress data comparison values for a calibration with magnitude data from the best-fit root model (f).The deviation of the modelled stress state to that provided from the root model for each of the three scenarios (g) is used in a linear regression to derive the boundary conditions required to compute the best-fit branch model (h).

Figure 6 .
Figure6.Results of the best-fit root model (b) and an alternative scenario that fits the S h min values equally well, but is calibrated against a lower S H max value (a).The overall difference that results from the different S H max values used for the calibration is expressed in the continuous scale of the regime stress ratio (RSR), which is between 0.5 (normal faulting regime), 1.5 (strike-slip), and 2.5 (thrust faulting regime)(Simpson, 1997) in the model volume.The horizontal stress anisotropy expressed in the ratio of S H max /S h min is shown on two cross sections which intersect below Munich.The differential stress (difference between the maximum and minimum principal stress, lowermost panel) is mapped on a surface which is vertically centred in the Malm α − γ units.

Figure 7 .
Figure 7.The generic branch model results are shown by means of slip tendency (ST) values (Morris et al., 1996) mapped on generic faults in the Chattian, Purbeck, and Malm units and by means of the fracture potential (FP) (Connolly and Cosgrove, 1999) displayed for the model volume of the basement.Both values vary from zero to one indicating low and high criticality.Note, that the colour map of these values is non-linear.The results clearly indicate that the generic faults are far away from failure with the largest value of ST of 0.3.The low FP values (max.0.38) give an estimate on how much fluid pressure would be needed to fracture the intact rock in a stimulation experiment to enhance the permeability.

Figure 8 .
Figure 8. Stratigraphy and model result of the root model along the borehole of the geothermal project in Sauerlach.Lines show the results of the best-fit root model: blue for the S h min magnitude, black for the vertical stress S v , and red for the S H max magnitude.The blue dots are formation integrity tests (FITs), which are a lower boundary for the magnitude of S h min and not used for calibration, the blue star represents the suspected LOT, the red star shows the S H max magnitude of 112 MPa used for calibration(Seithel et al., 2015).Shaded areas in the same colour around the lines show the range of model scenarios that fit equally well to the model-independent constraints.The dotted red line shows the S H max magnitude for the model scenario in Fig.6a.

Table 1 .
The stratigraphic units, their discretization, and according rock properties, which are present in the root and branch models.Units which are only preserved in parts of the root model are marked with *.

Table 2 .
The expected maximum variations in slip tendency (ST) introduced by the uncertainties of the model parameters.This comparison is made at 40 locations in the Malm α−γ and Purbeck target units and an arithmetic mean is computed for each model parameter.