Stress rotation – impact and interaction of rock stiffness and faults

It has been assumed that the orientation of the maximum horizontal compressive stress (S_Hmax) in the upper crust is governed on a regional scale by the same forces that drive plate motion. However, several regions are identified where stress orientation deviates from the expected orientation due to plate boundary forces (first-order stress sources), or the plate wide pattern. In some of these regions, a gradual rotation of the S_Hmax orientation has been observed.

Several second- and third-order stress sources have been identified in the past, which may explain stress rotation in the upper crust. For example, lateral heterogeneities in the crust, such as density and petrophysical properties, and discontinuities, such as faults, are identified as potential candidates to cause lateral stress rotations. To investigate several of these candidates, generic geomechanical numerical models are set up with up to five different units, oriented by an angle of 60^∘ to the direction of shortening. These units have variable (elastic) material properties, such as Young's modulus, Poisson's ratio and density. In addition, the units can be separated by contact surfaces that allow them to slide along these vertical faults, depending on a chosen coefficient of friction.

The model results indicate that a density contrast or the variation of Poisson's ratio alone hardly rotates the horizontal stress (≦17^∘). Conversely, a contrast of Young's modulus allows significant stress rotations of up to 78^∘, even beyond the vicinity of the material transition (>10 km). Stress rotation clearly decreases for the same stiffness contrast, when the units are separated by low-friction discontinuities (only 19^∘ in contrast to 78^∘). Low-friction discontinuities in homogeneous models do not change the stress pattern at all away from the fault (>10 km); the stress pattern is nearly identical to a model without any active faults. This indicates that material contrasts are capable of producing significant stress rotation for larger areas in the crust. Active faults that separate such material contrasts have the opposite effect – they tend to compensate for stress rotations.

1 Introduction

Knowledge of the stress tensor state in the Earth's upper crust is important for a better understanding of the endogenous dynamics, seismic hazard or exploitation of the underground. Therefore, several methods have been developed to estimate the stress tensor orientation and the stress magnitudes. Stress orientation data are compiled globally in the World Stress Map database (Zoback et al., 1989; Zoback, 1992; Sperner et al., 2003; Heidbach et al., 2010, 2018). Based on such data compilations, it was assumed that patterns of stress orientation on a regional scale are more or less uniform within tectonic plates (Richardson et al., 1979; Klein and Barr, 1986; Müller et al., 1992; Coblentz and Richardson, 1995).

The plate-wide pattern is overprinted on a regional scale by the contemporary collisional systems. Recent examples in Europe are the Alps (Reinecker et al., 2010), the Apennines (Pierdominici and Heidbach, 2012) or the Carpathian Mountains (Bada et al., 1998; Müller et al., 2010). Closely related to that is the variability of crustal thickness, density and topography (Artyushkov, 1973; Humphreys and Coblentz, 2007; Ghosh et al., 2009; Naliboff et al., 2012). It was suggested that remnant stresses due to old plate tectonic events are able to overprint stress orientation on a regional scale (e.g. Eisbacher and Bielenstein, 1971; Tullis, 1977; Richardson et al., 1979). Such old basement structures also present geomechanical inhomogeneities and discontinuities, which have the potential to perturb the stress pattern. However, pre-Cenozoic orogens (or “old” suture zones), often covered and hidden by (thick) sediments, were rarely indicated as causes of significant stress rotation. In many cases it is the opposite: old orogens have apparently no impact on the present-day crustal stress pattern, e.g. the Appalachian Mountains (Plumb and Cox, 1987; Evans et al., 1989) or Fennoscandia (Gregersen, 1992). Deviations from the assumed uniform plate-wide stress pattern (here called stress rotations) are observed recently in several regions, such as in Australia, Germany and North America (Reiter et al., 2015; Heidbach et al., 2018; Lund Snee and Zoback, 2018, 2020). However, these effects can only be partly explained by the topography or lithospheric structures.

The complex stress pattern in central–western Europe was a subject of several numerical investigations in recent decades (Grünthal and Stromeyer, 1986, 1992, 1994; Gölke and Coblentz, 1996; Goes et al., 2000; Marotta et al., 2002; Kaiser et al., 2005; Jarosiński et al., 2006). Apart from a recent 3-D model (Ahlers et al., 2020), the previous models were limited to 2-D. These 2-D models were able to reproduce some of the observed stress patterns by considering variable lateral elastic material properties or discontinuities.

However, 2-D models have some limitations: they have to integrate topography, crustal thickness and stiffness to one property, and they potentially overestimate the horizontal stress magnitude (van Wees et al., 2003; Ghosh et al., 2006). Furthermore, none of these previous studies investigated the impact of the influencing factors separately.

In this work, a series of large-scale 3-D generic geomechanical models is used to determine which properties can cause significant stress rotations at a distance (>10 km) from material transitions or discontinuities. The model geometry is inspired by the crustal structure and the stress pattern in the German Central Uplands, where the S_Hmax orientation is 120 to 160^∘. This is in contrast to a N–S orientation (∼0^∘) of S_Hmax to the north and to the south of the uplands (Fig. 1, Reiter et al., 2015). The basement structures there are striking 45 to 60^∘, which is almost perpendicularly to the observed S_Hmax orientation. The influence of the structures on the stress field will be tested with a generic variation of Young's modulus, Poisson's ratio, the density and vertical low-friction discontinuities, which separate the crustal blocks. Each property is tested separately first, to avoid interdependencies; possible interactions are tested afterwards.

2 Stress rotation in the upper crust

2.1 Concept of stress rotation

This study focuses on stress rotations that occur horizontally, i.e. in the map view. A vertically uniform stress field is assumed, which is consistent with previous studies (Zoback et al., 1989; Zoback, 1992; Heidbach et al., 2018). Stress rotations with depth are occasionally observed within deep wells (Zakharova and Goldberg, 2014; Schoenball and Davatzes, 2017), due to evaporites (e.g. Roth and Fleckenstein, 2001; Röckel and Lempp, 2003; Cornet and Röckel, 2012) or man-made activities in the underground (e.g. Martínez-Garzón et al., 2013; Ziegler et al., 2017; Müller et al., 2018).

On a map view, several potential sources of stress can superpose on another and the resulting stress at a certain point comprises the sum of all stress sources from those plate-wide to very local stress sources. Differences between the resulting stress orientation and the regional stress source can be described by the angular deviation (Sonder, 1990), which can be substantial and can lead to a change of the stress regime (Sonder, 1990; Zoback, 1992; Jaeger et al., 2007). The stress regime (Anderson, 1905, 1951) is defined by the relative stress magnitudes, which are a normal faulting regime ($S_{\mathrm{V}}>S_{\text{Hmax}}>S_{\text{hmin}}$), strike slip regime ($S_{\text{Hmax}}>S_{\mathrm{V}}>S_{\text{hmin}}$) and thrust faulting regime ($S_{\text{Hmax}}>S_{\text{hmin}}>S_{\mathrm{V}}$), where S_V is the vertical stress and S_hmin and S_Hmax are the minimum- and the maximum horizontal stress, respectively. The difference between the largest and smallest principal stress is the differential stress (${\mathit{\sigma }}_{\mathrm{D}}={\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{3}}$), while the deviatoric stress is the difference between the stress state and the mean stress ($\mathit{\delta }{\mathit{\sigma }}_{ij}={\mathit{\sigma }}_{ij}-{\mathit{\sigma }}_{\mathrm{m}}$;  Engelder, 1994)

Grünthal and Stromeyer (1986)Bell and Lloyd (1989)Bell and McCallum (1990)Sonder (1990)Grünthal and Stromeyer (1992)Grünthal and Stromeyer (1994)Spann et al. (1994)Zhang et al. (1994)Gölke and Coblentz (1996)Homberg et al. (1997)Mantovani et al. (2000)Marotta et al. (2002)Yale (2003)Jarosiński et al. (2006)Mazzotti and Townend (2010)

Table 1 Comparison of selected previous observations or models on the subject of stress rotation in the context of faults, elastic material properties, density or topography variation. The characters “X” and “V” indicate whether the property is included or varied; “(X)” means that the subject is included indirectly. The characters “ < ” and “ > ” indicate that significant rotation occurs near (<10 km) or at greater distance (>10 km) from the fault or material transition.

Stress rotation within this study means an angular deviation of the S_Hmax orientation from the large-scale stress pattern. In the following subsections, previous observations and models on the respective causes are reviewed and also summarized in Table 1.

2.2 Density contrast and topography

Variability of density within the crust or lithosphere has a significant impact on the stress state (Frank, 1972; Artyushkov, 1973; Fleitout and Froidevaux, 1982; Humphreys and Coblentz, 2007; Ghosh et al., 2009; Naliboff et al., 2012). Assameur and Mareschal (1995) showed that local stress increases due to topography and crustal inhomogeneities are in the order of tens of MPa, which is on the order of stresses resulting from the plate boundary forces.

Gravitational forces are also derived by surface topography (Zoback, 1992; Miller and Dunne, 1996). Within mountains, S_Hmax is oriented parallel to the ridge and perpendicular to the ridge at the base of the mountain chain. Along passive continental margins, effects similar to those due to topography can be observed (Bott and Dean, 1972; Stein et al., 1989; Bell, 1996; Yassir and Zerwer, 1997; Tingay et al., 2005; King et al., 2012).

Sonder (1990) investigated the interaction of different regional deviatoric stress regimes (δσ_ij) with stresses arising from buoyancy forces (σ_G) and observed a rotation of S_Hmax of up to 90^∘. According to that, S_Hmax rotates toward the normal trend of the density anomaly. If regional stresses are large, compared to stresses driven by a density anomaly ($\mathit{\delta }{\mathit{\sigma }}_{ij}/{\mathit{\sigma }}_G\gg \mathrm{1}$), the influence of a density anomaly is small and vice versa: if the regional stress is small compared to the stress driven by the density anomaly ($\mathit{\delta }{\mathit{\sigma }}_{ij}/{\mathit{\sigma }}_G\ll \mathrm{1}$), the impact of a density anomaly on the resulting stress field is large. In the case that both stress sources are on a similar level ($\mathit{\delta }{\mathit{\sigma }}_{ij}/{\mathit{\sigma }}_G\approx \mathrm{1}$), small changes of one of the stress sources are able to change the stress regime, and thus potentially the stress orientation.

2.3 Stiffness contrast

Mechanical stiffness describes the material behaviour under the influence of stress and strain. The focus here is on linear elastic material properties, characterized by Young's modulus and Poisson's ratio. Stress refraction between two elastic media can be calculated, but only at the interface of the two media, based on the known stress state on one side of the interface and Young's modulus on both sides (Spann et al., 1994). Stress rotation due to stiffness contrast is for example reported for the Peace River Arch in Alberta, Canada (Fordjor et al., 1983; Bell and Lloyd, 1989; Adams and Bell, 1991). Potential stress rotation is supported by several numerical studies (e.g. Bell and Lloyd, 1989; Grünthal and Stromeyer, 1992; Spann et al., 1994; Zhang et al., 1994; Tommasi et al., 1995; Mantovani et al., 2000; Marotta et al., 2002).

2.4 Discontinuities

Discontinuities are planar structures within or between rock units, where the shear strength is (significantly) lower than that of the surrounding rock. Genetically, discontinuities can be classified into bedding, schistosity, joints and fault planes. In the context of this study the term discontinuity refers to fault planes or fault zones. Similar to the Earth surface, (nearly) frictionless faults without cohesion act like a free surface in terms of continuum mechanics (Bell et al., 1992; Bell, 1996; Jaeger et al., 2007). One of the three principal stresses must be oriented perpendicular to the frictionless fault; the two remaining ones are parallel to the discontinuity. For this reason, the stress tensor rotates near a frictionless fault, depending on its orientation. Significant stress rotation in the context of faults is reported (Bell and McCallum, 1990; Adams and Bell, 1991; Yale, 2003; Mazzotti and Townend, 2010). However, Yale (2003) assumes that stress rotation occurs only within several kilometres from the fault. Large differential stress leads to a more stable stress pattern (Laubach et al., 1992; Yale, 2003), whereas low differential stresses allow a switch of the stress regime caused by faults. The impact of faults on stress rotation has been investigated analytically (Saucier et al., 1992) and by numerical models (e.g. Zhang et al., 1994; Tommasi et al., 1995; Homberg et al., 1997).

3 Regional setting

3.1 Stress orientation in central Europe

Crustal stress data from Europe have been collected since the 1960s (e.g. Hast, 1969, 1973, 1974; Greiner, 1975; Ranalli and Chandler, 1975; Greiner and Illies, 1977; Froidevaux et al., 1980; Kohlbeck et al., 1980), later as part of the World Stress Map database from Zoback et al. (1989) and more recently by Heidbach et al. (2018).

S_Hmax orientation in western Europe is 145${}^{\circ }\pm \mathrm{26}$^∘ and rotates clockwise by about 17^∘ (Müller et al., 1992) to the direction of absolute plate motion from Minster and Jordan (1978). This is in agreement with Zoback et al. (1989), who obtained a better fit for relative plate motion between Africa and Europe than for absolute plate motion. As the major causes of the observed stress pattern in western and central Europe, the ridge push of the Mid-Atlantic ridge and the collisional forces along the southern plate margins are identified (Richardson et al., 1979; Grünthal and Stromeyer, 1986; Klein and Barr, 1986; Zoback et al., 1989; Grünthal and Stromeyer, 1992; Müller et al., 1992; Zoback, 1992; Gölke and Coblentz, 1996; Goes et al., 2000).

A fan-like stress pattern has been observed in the western Alps and Jura mountains, where S_Hmax in front of the mountain chain is perpendicular to the strike of the orogen (Fig. 1). Müller et al. (1992) assume that these structures only locally overprint the general stress pattern. However, in light of the recently available data, it is assumed that the S_Hmax orientation is rather controlled by gravitational potential energy of the alpine topography than by plate boundary forces (Grünthal and Stromeyer, 1992; Reinecker et al., 2010).

The stress pattern in western and central Europe has been the subject of several modelling attempts in the last three decades (Grünthal and Stromeyer, 1986, 1992, 1994; Gölke and Coblentz, 1996; Goes et al., 2000; Marotta et al., 2002; Kaiser et al., 2005; Jarosiński et al., 2006). In particular, these previous studies investigated the impact of a lateral stiffness contrast in the crust (Grünthal and Stromeyer, 1986, 1992, 1994; Jarosiński et al., 2006; Kaiser et al., 2005; Marotta et al., 2002), the elastic thickness of the lithosphere (Jarosiński et al., 2006), the stiffness contrast of the mantle (Goes et al., 2000), a lateral density contrast or topographic effects (Gölke and Coblentz, 1996; Jarosiński et al., 2006), the post-glacial rebound in Scandinavia (Kaiser et al., 2005), and activity on faults (Kaiser et al., 2005; Jarosiński et al., 2006).

Stiffness variation in the lithosphere, e.g. in the Teisseyre–Tornquist Zone (TTZ) or the Bohemian Massif (BM), has been identified as a potential cause for the observed stress rotation in Central Europe (Grünthal and Stromeyer, 1986, 1992, 1994; Gölke and Coblentz, 1996; Reinecker and Lenhardt, 1999; Goes et al., 2000; Marotta et al., 2002; Kaiser et al., 2005). One example is the fan-shaped stress pattern in the North German Basin (NGB), with a rotation of S_Hmax from north-west in the western part to north-east in the eastern part of the basin as a product of the TTZ, which is the boundary between the Phanerozoic Europe (Avalonia) and the much stiffer Precambrian Eastern European Craton (Baltica).

Jarosiński et al. (2006) came to the conclusion that active tectonic zones and topography have major effects, whereas the stiffness contrast leads only to minor effects. Lateral variation of density does not have a significant impact on the stress pattern (Gölke and Coblentz, 1996); it causes only local effects. Finally, low differential stress allows significant stress rotation (Sonder, 1990; Grünthal and Stromeyer, 1992).

3.2 Basement structures in Germany

In large part, Germany consists of Variscan basement units, either exposed or covered by Post-Paleozoic basin sediments. The Variscan orogen is a product of the late-Paleozoic collision of the plates Gondwana and Avalonia (Laurussia) in late Devonian to early Carboniferous time, which lead to closure of the Rheic Ocean (Matte, 1986), and finally the formation of the super-continent Pangaea. Despite the fact that the European Variscides are well investigated in the last century and decades (e.g. Franke, 2000, 2006; Kroner et al., 2007; Kroner and Romer, 2013), it is for example still a matter of debate whether several microplates have been amalgamated in between or not.

Figure 1 Stress orientation in the German Central Uplands with the basement structural elements (separated by black lines), political boundaries (red) and major rivers (blue). Bars represent orientation of maximum horizontal compressional stress (S_Hmax); line length is proportional to quality. Colours indicate stress regimes, with red for normal faulting (NF), green for strike–slip faulting (SS), blue for thrust faulting (TF) and black for unknown regime (U). The Variscan basement structures introduced by Kossmat (1927) are visualized; the regional segmentation is as follows: BM = Bohemian Massif, MGCH = Mid-German Crystalline High, MZ = Moldanubian Zone, NPZ = Northern Pyllite Zone, RHZ = Rheno-Hercynian Zone, STZ = Saxo-Thuringian Zone, and VDF = Variscan Deformation Front. Other structures are as follows: MB = Molasses Basin; NGB = North German Basin, TS = Thor Suture, TTZ = Teisseyre–Tornquist Zone (redrawn after Franke, 2014; Grad et al., 2016).

Kossmat (1927) published the structural zonation of the European Variscides, which is still widely used (Fig. 1). The parts to the north-west of the Rheic Suture Zone are the Rheno-Hercynian Zone (RHZ) with the sub-unit of the Northern Phylite Zone (NPZ), both of Laurussian origin. South-east of the suture zone are the Mid-German Crystalline High (MGCH), the Saxo-Thuringian Zone (STZ) and the Moldanubian Zone (MZ); all except the MGCH were exclusively part of Gondwana.

The RHZ is exposed in the Rhenish Massif, in the Harz mountains and in the Flechtingen Hills. Dominant are Devonian to lower Carboniferous clastic shelf sediments (Franke, 2000; Franke and Dulce, 2017). These low metamorphic slates, sandstones, greywacke and quartzite are supplemented with continental and oceanic volcanic rocks, reef limestones and a few older gneisses. Further to the north of the RHZ are the sub-Variscan foreland deposits, consisting of clastic sediments and coal seams.

The NPZ is uncovered at the southern edge of the low mountain ranges Hunsrück, Taunus and eastern Harz. Petrologically it is probably the greenschist facies equivalent (Oncken et al., 1995) of the Rheno-Hercynian shelf sequence (Klügel et al., 1994), consisting of meta-sediments and within-plate metavolcanic rocks (Franke, 2000).

The MGCH is open in the Palatinate Forest, Odenwald, Spessart, Kyffhäuser, Ruhla Crystalline (Thuringian Forest) and Flechtingen Hills. It has been interpreted previously as a magmatic arc of the Saxo-Thuringian Zone. But Oncken (1997) assumes that the MGCH is composed of both Saxo-Thuringian and Rheno-Hercynian rocks. Composition and metamorphic grade vary considerably along the strike of the MGCH (Franke, 2000). It consists of late-Paleozoic sediments, meta-sediments, volcanic rocks, granitoides, gabbros, amphibolite and gneisses.

The Saxo-Thuringian Zone (STZ) is exposed in the Thuringian-Vogtlandian Slate Mountains, Fichtel Mountains, Ore Mountains, Saxonian Granulite Massif, Elbe Valley Slate Mountains and the Lausitz. It consists of Campro-Ordovician mafic and felsic magmatic rocks and late Ordovician to early Carboniferous marine and terrestrial sediments (Franke, 2000; Linnemann, 2004). These rocks underwent metamorphic overprint up to the early Carboniferous with different metamorphism stages up to eclogite- or granulite facies. These units are interspersed by late- or post-orogenic granites.

The MZ is exposed in the Bohemian Massif, the Bavarian Forest, the Münchberg Gneiss Massif, the Black Forest and the Vosges. They consist of mostly high-grade metamorphic crystalline rocks (gneisses, granulite, migmatite) and Variscan granites (Franke, 2000).

4 Model set-up

4.1 Model geometry

The chosen model geometry is inspired by the geometrical situation in the German Central Uplands (Fig. 1), but the overall intention is a generic model. To make it easy to understand, compass directions are used for the model description. The model geometry has a north–south extent of 400 and 300 km in the east–west direction, with a thickness of 30 km (Fig. 2). In the centre of the model, three diagonal units each with a width of 50 km are oriented 60^∘ from the north. The unit boundaries are vertically incident. A model variant is generated in which the unit boundaries allow free sliding, depending on a chosen friction coefficient. For each of the three central units, different material properties can be applied. The northernmost and southernmost block has always the same (reference) material properties, except for the realistic rock property scenario.

Figure 2 Reference model with the applied boundary conditions, used for all models, in map view and from the south. The model has a lateral extent of 300 km×400 km and a thickness of 30 km. It consists of five interconnected units, which have the same material properties. Blue visualizes the reference material (Table 2). The boundary conditions ban motion in the x direction on the western side, in the y direction on the northern side and in the z direction at the model base. A push of 400 m from the south and a pull of 60 m to the east is applied. The resulting S_Hmax orientation (north–south) at a depth of 1000 m is illustrated by the black bars. The red point (and line) indicates the location of the virtual well (Figs. 4 and 5). The four diagonal boundaries can be used as vertical faults with a chosen friction coefficient.

Figure 3 Selection of common elastic rock properties (Young's modulus and Poisson's ratio) and density (Turcotte et al., 2014). Coloured vertical bars indicate applied material properties; see Table 2.

4.2 Solution of the equilibrium of forces

The stress orientations in the models are investigated using the finite-element method (FEM). The usage of 3-D FEM models to investigate the stress state in the crust is a well-established technique (e.g. van Wees et al., 2003; Buchmann and Connolly, 2007; Hergert and Heidbach, 2011; Reiter and Heidbach, 2014; Hergert et al., 2015). The major reason that complex 2-D or 3-D models can be computed is the opportunity to use unstructured meshes.

The method in general computes the equilibrium of stresses arising from boundary forces (via displacement boundary conditions) and body forces (gravity) acting on the rock whose mechanical behaviour is characterized by a constitutive law and associated material parameters. The equilibrium of forces is represented by partial differential equations, which are solved numerically.

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle \frac{\mathit{\delta }{\mathit{\sigma }}_{ij}}{\mathit{\delta }{x}_j}}+\mathit{\rho }{x}_i=\mathrm{0},\end{array}$$

where δσ_ij is the variation of total stress, δx_j is the spatial change and ρx_j represents the weight of the rock section (ρ=density). Linear elastic material behaviour expressed by Hooke's law is assumed. Two material properties, Young's modulus (E) and Poisson's ratio (ν) are essential. The stress state in this study will be calculated based on defined displacement boundary conditions (Fig. 2).

The lateral resolution of the model is about 3 km, consisting primarily of hexahedrons and some wedge elements (degenerated hexahedrons). Resolution into depth ranges from 0.44 km near the surface to about 3.4 km at the base of the model. In total, about 166 000 elements were used. The model version with contact surfaces uses 1725 contact elements along each contact surface. Model discretization was performed with HyperMesh^® v.2019. The equilibrium of forces (body forces and boundary condition) is computed numerically using the Abaqus^®/Standard v.6.14-1 finite-element software.

4.3 Mechanical properties

The main subject of this study is to investigate the impact of the variation of elastic rock properties, density and friction along faults on stress orientation in the upper crust in the given geometrical setting outlined in the previous sections (Fig. 2). To do this, each parameter is tested individually. Figure 3 visualizes the range of density (ρ), Young's modulus (E) and Poisson's ratio (ν) of representative rocks, taken from a textbook (Turcotte et al., 2014).

The reference material for this investigation has a density of $\mathit{\rho }=\mathrm{2.7}\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$, a Poisson's ratio of ν=0.25 and a Young's modulus of E=50 GPa. Such a material could represent for example granite or limestone. Based on this reference material, a lower and higher material value is always defined (Table 2), which is within the range of common rock properties (Fig. 3). The material with a low density ($\mathit{\rho }=\mathrm{2.2}\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$) may represent sediments (sandstone, limestone, shale etc.), whereas the high-density material ($\mathit{\rho }=\mathrm{3.2}\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$) could represent a rock from the lower crust or the upper mantle. A low Poisson's ratio (ν=0.15) may represent sediments (sandstone or shale), and a high Poisson's ratio (ν=0.35) could represent ultramafic rocks. Soft material with a low Young's modulus (E=10 GPa) may represent sediments, pre-damaged rock or weathered rock. Again ultramafic rock is an example of a stiff rock, having a large Young's modulus (E=100 GPa).

Table 2 Young's modulus, Poisson's ratio and densities used in the models. Bold numbers indicate the properties used, which differ from those of the reference material.

Figure 4 The stress ratio k (Eq. 3) is plotted vs. depth. Stress in the reference model is marked with the bold green line. Additionally, several data and defined stress ratios from the literature are visualized for comparison (Heim, 1878; Herget, 1973; Brown and Hoek, 1978; Lindner and Halpern, 1978; McCutchen, 1982; Sheorey, 1994; Brudy et al., 1997; Hickman and Zoback, 2004).

Laboratory rock experiments in the past delivered friction coefficients of about μ=0.6 to 0.85 (Byerlee, 1978). However, recent investigations using realistic slip rates for earthquakes decreased estimated friction coefficients by 1 order of magnitude down to μ<0.1 (Di Toro et al., 2011). Faults are represented by cohesionless contact surfaces in the models. The used friction coefficients are 0.1, 0.2, 0.4, 0.6, 0.8 and 1.0, which covers both slow and fast slip rates.

4.4 Initial stress state

The present-day stress state in the crust is a complex product of several stress sources from the past to the present. In order to model the stress state an initial stress state is defined, which is in equilibrium with the body forces (gravity) and which subsequently undergoes lateral straining to account for tectonic stress. Sheorey (1994) provided a simple semi-empirical function (Eq. 2) for the stress ratio k (Eq. 3), where E is Young's modulus and z is the depth in kilometres.

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }k=\mathrm{0.25}+\mathrm{7}E\left(\mathrm{0.001}+{\displaystyle \frac{\mathrm{1}}{z}}\right)\text{(3)}& {\displaystyle }& {\displaystyle }k={\displaystyle \frac{S_{\text{Hmean}}}{S_{\mathrm{V}}}}={\displaystyle \frac{S_{\text{Hmax}}+S_{\text{hmin}}}{\mathrm{2}{S}_{\mathrm{V}}}}\end{array}$$

Sheorey's equation (Eq. 2) is a reliable stress ratio vs. depth estimation, when compared to real-world data (Fig. 4). The model is pre-stressed with zero horizontal strain boundary conditions. The pre-stressing method used here has so far been used several times (Buchmann and Connolly, 2007; Hergert and Heidbach, 2011; Reiter and Heidbach, 2014). The model is allowed to compact several times under application of the body forces (gravity) using a Poisson's ratio of ν=0.396 during that procedure only. During the pre-stressing procedure, models with contact surfaces have a very large friction coefficient (μ=10) to prevent slip. At a virtual well in the centre of the model (x=150 km; y=200 km) stress was extracted from the model and compared to the stress magnitude data, which are visualized in Figs. 4 and 5, showing a good fit to stress–depth distribution assumptions (Heim, 1878; Herget, 1973; Brown and Hoek, 1978; McCutchen, 1982; Sheorey, 1994) and measured magnitude ratios also (Brown and Hoek, 1978; Lindner and Halpern, 1978; Brudy et al., 1997; Hickman and Zoback, 2004).

4.5 Boundary conditions

The overall S_Hmax orientation on a virtual profile along longitude 11^∘ (Fig. 1) displays a north–south orientation in the North German Basin (NGB) and in the Molasse Basin (MB) north of the Alps, except the Variscan basement units in between. Correspondingly, a north–south orientation of S_Hmax is intended for the reference model. In order to generate a meaningful stress state in the model, appropriate boundary conditions are required, which are technically applied by a defined lateral displacement. Results from a virtual well in the model centre are compared with data from deep wells. An extension of 60 m (${\mathit{ϵ}}_x=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$) in the east–west direction and a shortening of 400 m (${\mathit{ϵ}}_y=-\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}$) in the north–south direction (Fig. 2) provide a good fit of the reference model to stress magnitudes from selected deep wells (Fig. 5,  Brudy et al., 1997; Hickman and Zoback, 2004; Lund and Zoback, 1999). By fitting the data, the focus was more on the observed S_hmin magnitudes and to a lesser extent on the S_Hmax magnitudes. The latter are less reliable, as they are usually not measured; they are calculated on the basis of several assumptions. The determined boundary conditions are used for all models.

Figure 5 The stress magnitudes are plotted as a function of depth. The stress components from the virtual well in the model are illustrated by the coloured lines. The location of the virtual well and boundary conditions used are shown in Fig. 2. Due to the applied initial stress conditions, the stress regime changes from thrust faulting at a depth of 400 m to strike slip faulting, and finally to a normal faulting regime at a depth greater than 5500 m. Published stress magnitude data are shown for comparison (Brudy et al., 1997; Lund and Zoback, 1999; Hickman and Zoback, 2004).

Figure 6 Influence of density on the stress orientation. Black bars represent the orientation of S_Hmax at a depth of 1000 m. Colours indicate the material properties used. The medium blue area uses the reference material properties ($\mathit{\rho }=\mathrm{2.7}\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$), the light blue material uses a lower density (g: $\mathit{\rho }=\mathrm{2.2}\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$), the dark blue a larger density (G: $\mathit{\rho }=\mathrm{3.2}\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$).

4.6 Generic model scenario's

The model geometry consists of five units (Fig. 2). The northern- and southernmost blocks are always assigned the reference material properties (Table 2). In between there are three diagonal units in which material properties are varied. Along the vertical borders within the model, friction properties can be used. The lower (L) or higher (H) values of the material properties with respect to the reference material (B) will be varied in the following way: LLL, HHH, LBL, BLB, etc. When the model geometry mimics discontinuities using contact surfaces, all contacts have the same friction coefficient. In the figures showing the results the label “|” indicates contact. For example, HLH with four contacts is |H|L|H|.

The S_Hmax orientation is visualized at a depth of 1000 m below the surface using a pre-defined grid, where the lateral distance to the material transition or discontinuity is >12.5 km, as far-field effects are the main interest of this study. The variation of density, Poisson's ratio, Young's modulus and friction coefficient will be tested first. In addition, the variation of Young's modulus is tested in interaction with low-friction contacts.

Table 3 Material properties used for the scenario using realistic rock properties for Variscan basement units; properties are estimated based on Turcotte et al. (2014). MGCH = Mid-German Crystalline High, MZ = Moldanubian Zone, NPZ = Northern Pyllite Zone, RHZ = Rheno-Hercynian Zone, STZ = Saxo-Thuringian Zone.

4.7 Realistic rock property scenario

A reality-based rock property scenario, inspired by the structural zonation of the European Variscides according to Kossmat (1927), is tested. The RHZ and the NPZ are dominated by clastic shelf sediments with a low- or mid-metamorphic overprint, which is made of slate (RHZ) and phyllite (NPZ). This zone, the RHZ and the NPZ together, is the most flexible one and will have the lowest Young's modulus (Table 3). The MGCH consists of granitoids or gabbros and their metamorphic equivalents (gneiss, amphibolite), meta-sediments, and some volcanites. Therefore, this zone is a stiff unit. The Saxo-Thuringian Zone (STZ) is dominated by meta-sediments, mafic and felsic magmatites and their metamorphosed equivalents, and some high-grade metamorphic rocks (granulite, eklogite). Taking all the different rock types into account, the STZ is stiffer than the RHZ and softer than the MGCH. Mechanically, the MZ can be represented by high-grade metamorphic rocks (gneiss, granulite, migmatite) and granitoids and will be a stiff unit, similar to the MGCH. Therefore, the unit stiffnesses are different: they are from slightly deformable to rigid in the following order: $\text{RHZ}\approx \text{NPZ}<\text{STZ}<\text{MGCH}\approx \text{MZ}$. Material properties used are estimated based on typical rock values (Table 3). The same initial stress procedure, boundary condition and visualization procedure are applied as previously described.

Figure 7 Influence of Poisson's ratio on the stress orientation. Black bars represent the orientation of S_Hmax at a depth of 1000 m. Colours indicate the material properties used. The blue area uses the reference material properties (ν=0.25), the light purple area is characterized by a low Poisson's ratio (p: ν=0.15) and the dark purple one by a large Poisson's ratio (P: ν=0.35).

5 Results

5.1 Density influence

To identify the influence of a density variation, the reference density ($\mathit{\rho }=\mathrm{2.7}\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$) in blue is varied using a small density (g: $\mathit{\rho }=\mathrm{2.2}\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$), which is coloured in light blue, and a large density (G: $\mathit{\rho }=\mathrm{3.2}\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$), which is dark blue (Fig. 6).

The low-density anomaly (ggg) results in a slight counter-clockwise ($-\mathrm{6}{}^{\circ }$) rotation of the S_Hmax orientation in the reference material near the anomaly (Fig. 6). Within the low-density units near the reference material, nearly no rotation is observed ($-\mathrm{1}{}^{\circ }$), but S_Hmax orientation turns counter-clockwise ($-\mathrm{8}{}^{\circ }$) in the centre of the material anomaly. The angular variation of S_Hmax crossing the units is of the order of 7^∘. The high-density anomaly (GGG) results in a slightly clockwise rotation ($+\mathrm{7}{}^{\circ }$) in the reference material near the anomaly. In the high-density unit near the reference material, S_Hmax is minimally influenced ($+\mathrm{1}{}^{\circ }$) but rotates further clockwise ($+\mathrm{12}{}^{\circ }$) in the centre of the anomaly. Based on that, the variation across the units is about 11^∘. The models with mixed densities in the three units show a clockwise rotation ($+\mathrm{10}{}^{\circ }$) of S_Hmax within the lighter material next to the denser units. The high-density units show a counter-clockwise rotation ($-\mathrm{7}{}^{\circ }$) next to the low-density unit; therefore, the total variation of S_Hmax is 17^∘.

In general, S_Hmax tends to be oriented parallel to the anomaly in low-density units and perpendicular to the anomaly in large density units. In the centre of the low-density units (ggg), the stress orientation becomes perpendicular to the overall structure. In the centre of the high-density units (GGG) the opposite is true, and S_Hmax becomes parallel to the structure.

5.2 Influence of Poisson's ratio

The influence of Poisson's ratio on the stress rotation is tested by variation of the reference Poisson's ratio (ν=0.25) using a lower one (p: ν=0.15) in light purple and a larger one (P: ν=0.35) in dark purple (Fig. 7). The models with only a lower (ppp: −1.5^∘) and only a higher Poisson's ratio (PPP: +2.2^∘) show only little S_Hmax rotation (Fig. 7). Mixed models with largest Poisson's ratio variation (pPp and PpP) have some counter-clockwise rotation in the low Poisson's ratio units ($-\mathrm{3.0}{}^{\circ }$) and a clockwise rotation in the high Poisson's ratio units ($+\mathrm{4.2}{}^{\circ }$). Therefore, the total variance of S_Hmax is about 7.5^∘.

5.3 Impact of Young's modulus

Figure 8 Influence of Young's modulus variation on the stress orientation. Black bars represent the orientation of S_Hmax at a depth of 1000 m. Colours indicate the used Young's modulus; the blue area uses the reference material properties (B: E=50 GPa), the green material uses a low Young's modulus (e: E=10 GPa) and the red material has a large Young's modulus (E: E=100 GPa).

The impact of Young's modulus is investigated using the reference material (B: E=50 GPa) in contrast to a softer material (e: E=10 GPa) in green and a stiffer material (E: E=100 GPa) in red (Fig. 8). The models with the soft units (eee, eBe and BeB) exhibit a strong clockwise S_Hmax rotation (+56^∘) in the units with the reference material and a counter-clockwise rotation in the softer units (−22^∘) near the material transitions. For the models with three soft units (eee) the S_Hmax orientation decreases to −5^∘ in the centre of the units. This means that the S_Hmax variation within the soft units is considerable (17^∘). The resulting total variation is 78^∘. The models with the stiff units (EEE, EBE and BEB) exhibit a gentle counter-clockwise rotation in the units with the reference material (−5.5 to −7^∘) next to the stiff units. Within the stiff units, a significant clockwise rotation (+20 to +25^∘) is apparent next to the reference units. In the model with three stiff units (EEE), the S_Hmax orientation decreases to (+5^∘) in the centre. This is a considerable S_Hmax variation of 15^∘ within the stiff units. The total variation is 31^∘.

For the models with alternating soft and stiff material units (EeE and eEe), the soft units exhibit a counter-clockwise S_Hmax rotation (−19 to −22^∘), whereas the stiff units display a clockwise rotation (+53 to +56^∘). Consequently, the total variation between the soft and stiff units is 72 to 78^∘. The general observation is that next to the material transition, S_Hmax rotates perpendicular to the anomaly for the compliant units and parallel for the stiff units.

5.4 Influence of faults

Several models with the reference material properties separated by three discontinuities (|B|B|B|) with a friction coefficient (μ) from 0.1 to 1 are tested. The low-friction coefficient (μ=0.1) leads to a counter-clockwise S_Hmax rotation of only −3^∘ (Fig. 9). The maximum observed fault offset is about 16 m. By increasing the friction coefficient to μ=0.2, the S_Hmax rotation is −2^∘; for μ=0.4, S_Hmax rotation is only −1^∘. For larger friction coefficients, the S_Hmax rotation is below $-\mathrm{1}{}^{\circ }$. As the S_Hmax rotation is too small for a visual differentiation, only the μ=0.1 model is shown in Fig. 9.

Figure 9 Influence of low-friction faults on the far-field stress orientation. Black bars represent the orientation of S_Hmax at a depth of 1000 m. All areas have the properties of the reference material (Table 2). White lines indicate cohesionless discontinuities (vertical faults). The model using a friction coefficient of μ=0.1 along the three discontinuities is shown. The other models with a larger friction coefficient (up to 1 and larger) have similar results; they are waived out because of the visual similarity.

Figure 10 Influence of Young's modulus in interaction with low-friction faults on the far-field stress orientation. Black bars represent the orientation of S_Hmax at a depth of 1000 m. Colours indicate the material properties used. The blue area uses the reference material properties, the green material uses a low Young's modulus and the red material has a larger Young's modulus; see Table 2. White lines indicate cohesionless vertical discontinuities (faults) with a friction coefficient of μ=0.1.

5.5 Stiffness variation combined with low-friction faults

The interaction between a significant Young's modulus contrast and a cohesionless contact with a low-friction coefficient (μ=0.1) is tested along all four discontinuities. The model with three stiff units (|E|E|E|) provides only little counter-clockwise rotation ($-\mathrm{4}{}^{\circ }$) in the reference material near the material transition (Fig. 10). Similar clockwise rotation occurs in the stiff units ($+\mathrm{4}{}^{\circ }$) near the material transition and decreases to the centre of the units ($+\mathrm{1}{}^{\circ }$). The total S_Hmax variation is about 8^∘.

The model with the soft units and the low-friction discontinuities (|e|e|e|) shows larger rotations than the model with stiffer units. Clockwise rotation of +19^∘ occurs in the reference material and counter-clockwise rotation of −13^∘ in the soft units. This decreases towards the centre of the soft units (−9^∘). Overall rotation is about 32^∘.

The models with alternating stiffnesses and low-friction discontinuities (|E|e|E| and |e|E|e|) generate a counter-clockwise rotation of about −10 to −12^∘ in the soft units. Within the stiff units, the S_Hmax orientation is in the range of +2 to +7^∘. The total variation is up to 19^∘. The maximum observed fault offset is about 10 to 15 m.

5.6 Stress rotation for realistic material properties

The resulting S_Hmax orientation (Fig. 11a) of the model using realistic material properties (Table 3) indicates counter-clockwise rotation in the RHZ and NPZ and clockwise rotation within the MGCH and MZ units. The overall pattern of the simple model (Fig. 11a) shows only limited similarity with the observed and the mean S_Hmax orientation on a regular grid using a search radius of 150 km and a quality and distance weight (Fig. 11b and c). However, some similarities can be observed. For example, the simple model (Fig. 11a) shows a clockwise rotation from the NPZ to the MGCH and counter-clockwise from the MGCH to the STZ. In Figure 11b these areas show similar, but less pronounced, rotation of S_Hmax. The north-north-east S_Hmax orientation within the central part of the MGCH is similar between the model, the data and the mean S_Hmax orientation (Fig. 11a–c).

Figure 11 Comparison of orientations of the maximum horizontal stress (S_Hmax). The equivalent regions are the RHZ (Rheno-Hercynian Zone), NPZ (Northern Phyllite Zone), the MGCH (Mid-German Crystalline High), the STZ (Saxo-Thuringian Zone) and the MZ (Moldanubian Zone). (a) Model results, application of estimated material properties of the Variscan units (Table 3). Black bars represent the S_Hmax orientation at a depth of 1000 m. (b) Bars indicate the S_Hmax orientation data (Heidbach et al., 2018), and quality is indicated by shades of grey; see legend. (c) Mean S_Hmax orientation on a 150 km search radius with a distance and quality weight (n>3) using the tool stress2grid (Ziegler and Heidbach, 2017). Panels (b) and (c) have the same extent as Fig. 1.

6 Discussion

6.1 Model simplification

This study investigates the influence of elastic material properties, density and friction coefficient at vertical faults on the orientation of S_Hmax. The focus is not on stress rotation close to the material transition or discontinuity (<10 km), the priority is on the far-field effects (>10 km). Although the model is inspired by a particular region, the goal is to gain a better understanding on how the variable material properties affect the stress orientation. For this reason, the model geometry is very simple and some of the material properties used may have no proper natural equivalent.

Chosen properties are constant over a depth of 30 km, which is unlikely. Even for a given lithology, the properties can change with depth, as a result of the acting gravity and compaction, especially for sediments. Each lithological unit is at least partially affected by these changes. Linearly increasing rock properties with depth would account for this and be a more realistic representation. But this would not affect the resulting stress pattern, especially since a vertically uniform stress field is assumed (Zoback et al., 1989; Zoback, 1992; Heidbach et al., 2018), with a few exceptions.

The simple generic models neglect various rheological processes in the crust by applying linear-elastic material law. However, the overall geometry seems reasonable, as the brittle domain or elastic thickness of the lithosphere (Te), which is a measure of the integrated stiffness of the lithosphere, is of the order of 30 km and more in central Europe (Tesauro et al., 2012). The Moho depth in Germany or central Europe is also about 30 km (Aichroth et al., 1992; Grad and Tiira, 2009). Jarosiński et al. (2006) for example used a range of Te=30–100 km for their model of central Europe. Furthermore, results are represented and discussed mainly for a depth of 1000 m where elastic behaviour is certainly predominant.

The scenario models were tested with an additional very stiff mantle (Table 3) with a thickness of 30 km. This had no influence on the observed stress pattern at a depth of 1000 m. However, the models with the same geometry but a total thickness of only 10 km resulted in much lower stress rotation. Therefore, the elastic thickness of the lithosphere and the aspect ratio of thickness and width of the units are important constraints for the possible stress rotation. The depth at which the stress orientation is plotted is also important, as the stress rotation decreases with depth (Fig. 12), so that it disappears at about 10 km depth for the used configuration. As homogeneous material properties are used, smaller scaling of results seems to be reasonable, considering the aspect ratio.

Figure 12 North–south depth profiles displaying the S_Hmax orientation colour-coded for models with a variable Young's modulus. In the model without the discontinuities (eEe), S_Hmax is oriented around 40^∘ in the stiffer units next to the softer units near the Earth surface. A similar orientation can be observed in the soft units in the deepest parts. In contrast to that, in the model with the same material properties but low-friction faults (|e|E|e|), the S_Hmax orientation is nearly north–south for all units and depths. (Small coloured dots are artefacts.) The discontinuities with a low-friction coefficient counterbalance stress rotation due to the stiffness contrasts.

All models were loaded with the same displacement boundary conditions (Fig. 2). This results in slightly different stress magnitudes due to the variable material properties. Since these models have different mechanical properties depending on the unit, the question would arise, in which of the units identical stress magnitudes should be achieved? Even if each model were calibrated individually, this would not significantly change the results, as both the stress regime and stress orientation would remain nearly constant for slightly different boundary conditions. Therefore, constant boundary conditions are reasonable and applied to all scenarios.

6.2 Stress rotation by density contrast

The lateral variation of the density is responsible for S_Hmax rotation in the range of 7 to 17^∘ (Fig. 6). In general, the S_Hmax rotates in the low-density units slightly toward parallel to the high-density unit (+10^∘), whereas S_Hmax rotates in the high-density units a little bit in the direction to the low-density units (−7^∘).

Taking a broad range of sediments into account (evaporites, shale, sandstone or limestone), they could have even a lower density than the lowest value used ($\mathit{\rho }=\mathrm{2.2}\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$). Most probably, models with a lower stiffness would result in larger stress rotation. However, sediments with a low stiffness could reach a thickness of several thousand metres, but not of the order of the model depth of 30 km or with such a low density due to increasing compaction with depth. Therefore, the impact of density variation on the stress orientation in nature will be much smaller, or on a very local scale. This agrees with the results of Gölke and Coblentz (1996), where a lateral density variation did not have a significant impact on the stress pattern; only local effects are observed.

This assumption seems to be a contradiction to the fact that the gravitational load is one of the main sources of stress in the Earth's crust. However, a density anomaly is a much smaller influencing variable on the stress state than density. According to Sonder (1990), the resulting stress rotations depend on the relative influence of regional stress sources as opposed to the density anomaly. Depending on the model scenarios used, the influence of the boundary conditions (regional stress sources) appears to be greater than that of the density anomaly. Therefore, the model results are probably not representative for regions with small horizontal differential stresses.

6.3 Stress rotation due to a variation of Poisson's ratio and Young's modulus

Model results suggest that the variation of Poisson's ratio can be responsible for a S_Hmax rotation of up to 7.5^∘ (Fig. 7). This is below the uncertainties of stress orientation estimations of about ±15^∘ and more (Heidbach et al., 2018). Therefore, the variation of Poisson's ratio can be neglected as a potential source of significant stress rotation.

In contrast to that, the lateral variation of Young's modulus can lead to significant S_Hmax rotation (Fig. 8). For the geometry and material parameters used, the relative rotations are up to 78^∘, which is not far from the maximal possible rotation of 90^∘. The largest rotation occurs in the units with a lower Young's modulus, for example the eee model has a total rotation of 78^∘, whereas the EEE model causes only a 31^∘ rotation. This is not surprising as Young's modulus is simply a measure of the stiffness. Therefore, the largest stress rotation due to stiffness contrast will happen in the soft units, not in the rigid ones. From this, it can be deduced that for units with smaller Young's modulus, the stress rotation is even greater.

S_Hmax will be oriented parallel to the structure for stiff units and perpendicular for soft units, which agrees with the literature (Bell, 1996; Zhang et al., 1994). The largest stress rotation occurs nearest to the material transition and decreases with distance to the material transition, similar to other models (Spann et al., 1994). Similar impacts of stiffness contrast have been described in previous studies (Grünthal and Stromeyer, 1992; Spann et al., 1994; Tommasi et al., 1995; Mantovani et al., 2000; Marotta et al., 2002). In contrast to that, Jarosiński et al. (2006) found that a stiffness contrast has only minor effects. But they did not test the stiffness contrast separately; they applied it only in combination with active faults in between the units. However, this agrees with the results of this study, as active faults balance stress rotation by stiffness contrast.

Within the units with a small Young's modulus, significant deformation is possible. For example, within the eEe model (Fig. 8), the soft units in green will be sinistrally deformed. The stiff unit in red cannot be deformed in the same way. But as the units are connected, the stiff unit is affected by the tangentially acting stress source. This leads to a S_Hmax orientation parallel to the structure, within the stiff units. As the soft one allows such deformation, S_Hmax will be oriented normally to the stiff unit.

At the interface between stiff and soft units differential stresses are greatest, as both units are differently deformable. This fits with the observation of concentrated intra-plate earthquakes around cratons (Mooney et al., 2012). On a smaller scale this has been observed for stiff sedimentary layers or rigid dykes, which attracts the occurrence of seismicity (Roberts and Schweitzer, 1999; Ziegler et al., 2015).

The observed radial stress pattern to the south of the Bohemian Massif (Reinecker and Lenhardt, 1999) agrees well with this study, where S_Hmax in the soft sediments of the Upper and Lower Austrian basin is perpendicular to the stiff crystalline Bohemian Massif. This is more ambiguously the case for the fan-shaped pattern in the western and northern part of the Alpine molasse basin (Grünthal and Stromeyer, 1992; Kastrup et al., 2004; Reinecker et al., 2010). The reason for this could be a lateral stiffness contrast of the rock, next to the topographic features of the mountain chain and the overall crustal structure. When comparing the stress rotation, it is important to consider the respective depth (see Fig. 12). For example, data in the north-western Alps originate from focal mechanisms, and in the foreland of the central Alps, the majority of data are from wells, which are more shallow (Reinecker et al., 2010).

Substantial stress rotations are not observed along major pre-Mesozoic boundaries and sutures in the eastern United States, like the Greenville front, a suture from Missouri to New York, or in the Appalachian Mountains (Zoback, 1992). Gregersen (1992) reports the same for Fennoscandia. In the case that these tectonic boundaries did not provide a significant stiffness transition, it is not a contradiction to this study. The mechanical contrast is important, not the relative ages.

6.4 Comparison of stress rotation due to elastic material properties

The rotation of S_Hmax perpendicular (counter-clockwise) to the structure can be observed most clearly in material with a lower Young's modulus next to a material transition, up to −22^∘. Rotation in the same direction, but with a lower amount, is observed in rocks with a larger density or a smaller Poisson's ratio. Within the units with a greater Young's modulus, S_Hmax rotates significantly parallel (clockwise) to the material transition, up to 56^∘. Similar rotation with a smaller magnitude can be observed in the low-density units or in the units with a larger Poisson's ratio. As rocks with a larger Young's modulus will usually have a larger density and vice versa (Fig. 3), real rocks will have less S_Hmax rotation as suggested by these generic models. But the aim of this study is to test and combine the possible range of variation, in order to identify the most important causes (Fig. 13).

Figure 13 Comparison of resulting maximum stress rotation, based on the geometry used (Fig. 2) and the varied properties (Table 2).

6.5 Failure criteria

As only elastic material properties are used, failure is not possible. To study the influence of this simplification, two models (EEE and eee) have been calculated using two different Coulomb failure criteria. The models are run first with a cohesion (C) of 30 MPa and a friction angle (FA) of 40^∘ and in addition with C=10 MPa and a FA of 30^∘.

For the EEE model with C=30 MPa and FA = 40^∘, no failure will be reached (yield criteria <1). For the eee model using that criteria and for both models (EEE and eee) using C=10 MPa and FA = 30^∘ failure occurs. Conditions of failure or close to failure (Yield criteria ∼1 or >1) occur only near the surface (a few kilometres) and close to the material transition (∼20–30 km). Around the material transition, (near) failure can be observed within the stiff units only. For the EEE model with C=10 MPa and FA = 30^∘, failure is more spaciously distributed near the surface.

In the case of failure, the S_Hmax orientation will be balanced, which means that S_Hmax rotates back in the north–south orientation, similar to the applied boundary conditions. In general, failure compensates for stress rotation in the same way as low-friction contact surfaces (faults). However, the stress orientation in the models that account for failure shows a similar stress pattern in the pre-failure phase to the models without failure. As a conclusion from this observation, the model results showing significant stress rotation are still valid for solid rocks.

6.6 Effect of faults on stress orientation

According to the model results, the influence of low-friction faults can be neglected concerning the orientation of the far-field stress pattern for homogeneous units. The low-friction faults (C=0, μ=0.1; FA 5^∘) lead to only 3^∘ S_Hmax rotation at a distance of about 12.5 km next to the fault zone. Observed stress rotation is lower than 1^∘ for a friction coefficient μ>0.4. This is not in contrast to the strong stress perturbation, observed in the vicinity of faults; as one of the three principal stresses must be oriented perpendicular to a fault, the two remaining ones are parallel to the discontinuity (Bell et al., 1992; Bell, 1996; Jaeger et al., 2007). Observations from meso-scale outcrops indicate stress perturbation within 2 km (Petit and Mattauer, 1995) or less than 1 km to a fault (Rispoli, 1981); larger stress perturbation can be observed at the termination of the fault (2–3 km). If S_Hmax is parallel next to the fault, it will rotate by 90^∘ at the termination of the fault (Rispoli, 1981; Osokina, 1988).

Yale (2003) suggests significant stress rotation as a product of active faults within a distance of several hundred metres for large differential stress provinces and several kilometres for regions with small differential stresses. This is supported by observed stress rotations near a fault within a range of a few hundred metres to a few kilometres (Brudy et al., 1997; Yassir and Zerwer, 1997). However, not all observed stress rotation agrees with the presented models, like observations offshore of eastern Canada (Bell and McCallum, 1990; Adams and Bell, 1991), where stress rotation occurs at a distance of about 10–15 km to a fault. Whether this is due to an inaccurate localization of the fault, a low Young's modulus or other causes cannot be clarified here.

Numerical models investigating stress rotations near a fault provide stress rotation between 20 and 60^∘, next to the fault, depending on the fault strike, the boundary conditions, and the friction or weakness of the fault. Near the termination of the fault, stress rotation increases to 50–90^∘ (Zhang et al., 1994; Tommasi et al., 1995; Homberg et al., 1997). However, rotation is observed by these models only within 2–3 elements, away from the discontinuity, which are anyway needed to distribute the deformation by such numerical models. Therefore, observed distances of rotation within these models are not considered here. To avoid this influence of an overly coarse mesh, the orientation of S_Hmax in this study is displayed at least four elements away from the contact surface.

6.7 Effect of faults combined with stiffness contrasts on stress orientation

The models with low-friction faults and a variable stiffness (Fig. 10) illustrate much lower stress rotation than the models without the faults (Fig. 8). The reason is that the soft units cannot transfer tangential shear stresses to the stiffer units. Therefore, each unit can be deformed independently from each other.

It seems to be that discontinuities play an important role in reducing stress rotations, produced by lateral Young's modulus variation (or other reasons). Regarding the used model geometry and materials, the S_Hmax variation is reduced for the soft models from 78 to 32^∘, for a comparison of eee and |e|e|e| in Figs. 8 and 10, using a friction coefficient of μ=0.1. Also for the mixed models, a reduction from 78 to 19^∘ is significant. Much lower is the rotation for the stiff model, with a reduction from 31 to 8^∘ (EEE in contrast to |E|E|E|).

The influence of low-friction faults in combination with variable mechanical properties was also investigated for the models with a density contrast. The observed effects are limited; therefore, presentation of the results is omitted.

6.8 Depth variation

The interaction of units with a variable Young's modulus and presence or non-presence of low-friction faults is visible in Fig. 12. While the model without active faults displays significant stress rotation, the same setting with low-friction faults shows only little rotation. The observed stress rotation within the model without faults strongly depends on the depth. In the soft units, S_Hmax rotates counter-clockwise near the surface (0 to −8 km). In contrast to that a clockwise rotation can be observed at greater depth (18–30 km). The stiffer units show a clockwise rotation in the upper part, about 0 to −8 km, which changes slightly to a counter-clockwise orientation in the deeper part.

This could be an indication that stress rotations due to stiffness contrasts in or near sedimentary basins can be significantly greater than in deeper material transitions, such as in a buried crystalline basement. This is all the more likely as sediments tend to be less stiff than crystalline rock.

6.9 Model using the Variscan zone rock properties

Figure 11 presents comparatively the results of S_Hmax orientation of the model using the chosen material properties, inspired by the Variscan units (a), the S_Hmax orientation data (b) and the averaged orientation of S_Hmax on a regular grid (c). A limitation of the mean stress orientation on a regular grid (Fig. 11c) is the calculation based on distance, and not depending on the specific unit. The similarity of the S_Hmax orientation is not very convincing, because the overall pattern cannot be reproduced. Some of the deviations from the trend are similar, some are not.

There are probably several reasons why the simple model is not able to reproduce the observed stress pattern in the German Variscides (Fig. 11). First of all, only one single elastic material composition represents each of the Variscan units. The model did not reproduce the complex and uncertain vertical variability of the deeper structures (Franke et al., 1990; Aichroth et al., 1992; Blundell et al., 1992), as no deep wells are present there. Only refraction seismic profiles from the 1980s (DEKORP) and their interpretations are available (Meissner and Bortfeld, 1990). Consequently, besides the uncertain structures, the material properties and the dip of the unit boundaries are also uncertain. More complex geometries and variable dip angle may result in different stress patterns to the ones obtained for vertical discontinuities. However, studying such variability is beyond the scope of this study.

Of course, it could also be possible that the units are decoupled. Thus, a scenario was calculated in which the units are separated by low-friction faults. But this did not provide a better fit of the stress orientation in contrast to the observation. Furthermore, there are no seismic or geodetic indications for such a decoupling of the units in that region.

Structures outside the Variscan zonation may also play a role. To the south, the stress pattern in the Molasses Basin is probably more governed by the structure of the Alpine chain (Reinecker et al., 2010) than older structures. The fan-shaped stress pattern in the eastern part of the North German Basin has been explained as an effect of the close boundary to the stiff Eastern European Craton along the north-west to south-east striking TTZ (Grünthal and Stromeyer, 1986, 1992, 1994; Gölke and Coblentz, 1996; Goes et al., 2000; Kaiser et al., 2005; Marotta et al., 2002). This interpretation agrees well with the results of the models, where S_Hmax becomes perpendicular in a soft unit (NGB) directed to a stiff region, like the East European Craton (e.g. model eee in Fig. 8).

7 Conclusions

The effect of varying elastic material parameters (Young's modulus and Poisson's ratio), density and low-friction discontinuities on the stress pattern in map view is investigated. Each property is tested separately to avoid interdependencies. This is performed with generic 3-D models using the finite-element method. Three units of variable material properties are included within the models, with boundary conditions determining the overall S_Hmax orientation. The variation of density and Poisson's ratio lead to small rotation (17 and 7.5^∘) of S_Hmax. In contrast, stiffness variation is able to produce significant stress rotation of 31 to 78^∘. Therefore, variation of Young's modulus in the upper crust is a potent explanation for observed stress rotation. Faults are represented in the models by cohesionless contact surfaces. The observed stress rotation in the far field due to low-friction faults (μ=0.1) is less than 3^∘. Implementation of low-friction discontinuities in models with a Young's modulus anomaly results in much smaller S_Hmax rotation, of the order of 8 to 32^∘. It follows that faults do not produce far-field stress rotation, but rather compensate for stress rotation that is an effect of the Young's modulus anomaly or other causes. Comparison of the model results with the observed stress orientation in the region that inspired the models provides only limited consistency. Nevertheless, the studies clearly show that fault systems are hardly the source of stress rotations on length scales of 100 km or larger. Furthermore, the study indicates that strength contrasts are promising candidates that have the potential to explain the slight stress pattern rotations in intraplate settings where topography and low-friction fault systems are missing.