Impact of faults on the remote stress state

The impact of faults on the contemporary stress field in the upper crust has been discussed in various studies. Data and models clearly show that there is an effect, but so far, a systematic study quantifying the impact as a function of distance from the fault is lacking. In the absence of data, here we use a series of generic 3-D models to investigate which component of the stress tensor is affected at which distance from the fault. Our study concentrates on the far field, located hundreds of metres from the fault zone. The models assess various techniques to represent faults, different material properties, different boundary conditions, variable orientation, and the fault's size. The study findings indicate that most of the factors tested do not have an influence on either the stress tensor orientation or principal stress magnitudes in the far field beyond 1000 m from the fault. Only in the case of oblique faults with a low static friction coefficient of μ=0.1 can noteworthy stress perturbations be seen up to 2000 m from the fault. However, the changes that we detected are generally small and of the order of lateral stress variability due to rock property variability. Furthermore, only in the first hundreds of metres to the fault are variations large enough to be theoretically detected by borehole-based stress data when considering their inherent uncertainties. This finding agrees with robust stress magnitude measurements and stress orientation data. Thus, in areas where high-quality and high-resolution data show gradual and continuous stress tensor rotations of >20^∘ observed over lateral spatial scales of 10 km or more, we infer that these rotations cannot be attributed to faults. We hypothesize that most stress orientation changes attributed to faults may originate from different sources such as density and strength contrasts.

1 Introduction

The crustal stress field is a key driver of geodynamic processes such as the earthquake cycle (Heidbach and Ben-Avraham, 2007; Wang et al., 2015; Hardebeck and Okada, 2018; Brodsky et al., 2020) and is of great importance for the safe exploitation of geo-reservoirs and storage of energy or waste in the subsurface (Fuchs and Müller, 2001; Zoback, 2010; Smart et al., 2014). In this context the interaction between the stress field in the Earth's upper crust and pre-existing faults is a crucial issue (Yale, 2003; Schoenball and Davatzes, 2017; Blöcher et al., 2018; Kruszewski et al., 2022; Röckel et al., 2022; Li et al., 2023).

For practical applications it is important to understand and to quantify on which spatial scale the fault changes the stress state. Exemplified with the site selection process for a deep geological repository for high-level radioactive waste, the interest is to know the distance to a fault at which no significant changes in the stress components occur in order to build the repository in a rock volume with homogeneous stress field conditions. In contrast to this, deep geothermal exploration targets faults or fault networks since they provide higher permeabilities compared to the rock matrix. Thus, the changes in the stresses in the near field of the fault and in its core or fracture network are of key interest to assess its dilation tendency (Moeck and Backers, 2011; Seithel et al., 2019; Ferrill et al., 2020). Stress perturbations are also significant for evaluating secondary fracturing near faults and its associated permeability, which encompasses joint orientation, secondary faulting, and bed-parallel slip (e.g. Kattenhorn et al., 2000; Maerten et al., 2002; Delogkos et al., 2022).

One of the key questions is on what spatial scale faults change the stress field and how to quantify which stress components are affected. The only component of the 3-D stress tensor that is systematically compiled is the orientation of maximum horizontal stress (S_Hmax, Heidbach et al., 2004, 2018). Areas with high data density revealed that the S_Hmax orientation can rotate significantly on scales from tens to hundreds of kilometres (Tingay et al., 2006; Heidbach et al., 2007; Rajabi et al., 2017b; Heidbach et al., 2018; Lund Snee and Zoback, 2020). The cause of this spatial variability has been investigated with generic geomechanical–numerical and analytical modelling (e.g. Sonder, 1990; Reiter, 2021). These studies show that stiffness, strength, and density contrasts are certainly key drivers of spatial distributed changes in the S_Hmax orientation.

Besides these findings, it was also hypothesized that active faults can cause rotations or magnitude variations as well (Dart and Swolfs, 1992; Yale, 2003; Faulkner et al., 2006; Konstantinovskaia et al., 2012; Schoenball et al., 2018; Li et al., 2023). This is confirmed on the borehole scale since logging data show stress rotations on the metre scale by means of abrupt changes in the orientation of borehole breakouts and drilling-induced tensile fractures (e.g. Barton and Zoback, 1994; Zoback et al., 2011; Rajabi et al., 2017c). It clearly showed that stress rotations on scales of 1 m to several hundred metres indeed occur due to faults and that the amount of rotation changes with distance to the fault core (Hickman and Zoback, 2004). Significant variation of stress magnitudes in the vicinity of faults has been reported for China and Scandinavia (Stephansson and Ångman, 1986; Li et al., 2023), but from these studies it is not clear which stress tensor component is affected as a function of distance to the fault. Furthermore, the mix of different methods that are used to estimate the stress parameter from very shallow locations near the surface and the lack of a rigorous uncertainty assessment make it difficult to assess whether the observed changes are significant and if they can be exclusively attributed to the nearby fault.

The only methods to test this are generic models, using geomechanical–numerical methods. There are several technical methods available to represent faults or fault zones numerically; for a method overview see Henk (2020). When using the continuum method, a fault is represented by selected elements with different behaviour, e.g. a lower Young's modulus (e.g. Cappa and Rutqvist, 2011), plastic behaviour (e.g. Mohr–Coulomb), or viscous behaviour. In contrast to that, using the discontinuous method, the fault is represented by contact elements (e.g. Buchmann and Connolly, 2007; Hergert et al., 2015) which allow offset along these structures (Fig. 1, Table 1). The finite-element method (FEM) is often used for such studies. Another discontinuous method, where the geometry is divided into several individual elements (circles or spheres, etc.), is the discrete-element method (DEM, e.g. Cundall and Hart, 1992; Yoon et al., 2014), which will not be used here. Physical models, using a photo-elastic material (e.g. de Joussineau et al., 2003), are also an option.

Figure 1 The general structure of a fault is described by the fault core, the damage zone, and the host rock (e.g. Caine et al., 1996; Faulkner et al., 2003). The purpose of this study is not to explore the effect of a fault on the stress state in the near field. This would include the fault core, the damage zone, and the neighbouring host rock. The study is focused on the far-field stress state, which is located several tens or hundreds of metres away from the fault and can extend up to a few kilometres at most. Numerical models typically employ one or a combination of two principle technical fault representations. Contact surfaces are a discontinuity within the mesh, where relative offset of the mesh is allowed, mainly depending on the friction. The second method uses a continuous mesh with elements having a lower stiffness or a failure criterion which results in a distributed deformation within the defined fault representation elements.

Tommasi et al. (1995)Buchmann and Connolly (2007)Xing et al. (2007)Hergert et al. (2011)Reiter and Heidbach (2014)Pereira et al. (2014)Hergert et al. (2015)Franceschini et al. (2016)Zhang et al. (2016)Meier et al. (2017)Schuite et al. (2017)Treffeisen and Henk (2020b)Reiter (2021)Homberg et al. (1997)Sánchez D et al. (1999)Maerten et al. (2002)McLellan et al. (2004)Camac and Hunt (2009)Cappa (2009)Cappa and Rutqvist (2011)

Table 1 List of some studies exemplifying the use of either a continuous or discontinuous mesh for fault representation (Fig. 1). Discontinuities are represented by contact elements or comparable methods (Contact). Another method of modelling faults is utilizing a continuous mesh that possesses a material definition slightly or significantly weaker (elastic, plastic, or viscous). These are 2-D elements within a 2-D mesh or 3-D elements in a 3-D mesh (Volume). Many models apply the finite-element method (FEM), while others use the finite-difference method (FDM), finite-volume method (FVM), or discrete-element method (DEM). The list does not claim to be complete.

The impact of faults has also been investigated by several authors using forward models. These studies (e.g. Table 1) either focus on how to technically implement faults into geomechanical–numerical models (Prévost and Sukumar, 2016; Treffeisen and Henk, 2020b) or on specific geological settings (Chéry et al., 2004; Fitzenz and Miller, 2001; Hergert and Heidbach, 2011; Meier et al., 2017; Yoon et al., 2017). As an example, Fig. 2 plots stress components along a horizontal line at sea level within a model from northern Switzerland (Hergert et al., 2015). The magnitudes of the stress tensor vary significantly close to the faults. However, resulting stress changes are affected by other factors too, such as topography or variable material properties.

Figure 2 Plot of stress components along a NNW–SSE profile at approximately 400 m depth (sea level) within the Nördlich Lägern model (Hergert et al., 2015). The largest and smallest horizontal principal stress (S_Hmax and S_hmin) as well as the vertical stress (S_V) are shown. Additionally shown is the topography from the model. The locations of the implemented Siglistdorf and Stadel–Irchel faults are indicated by vertical black lines. Stress magnitude changes are significant next to the faults, but stresses are also variable due to a variable topography, rock stiffness, or other factors. The significant variation in S_Hmax is attributed to material changes, as the stratigraphic boundaries dip slightly towards the south.

Previous studies show that faults certainly have an impact, but a systematic approach is still missing. They do not provide a quantification of which component of the stress tensor is affected by the stress changes near the fault. In this paper we systematically investigate the change in individual stress tensor components with distance to the fault. In particular, we determine the changes in the maximum and minimum horizontal stress (S_Hmax, S_hmin), the vertical stress (S_V), and the von Mises stress as well as the orientation of the stress tensor by means of the S_Hmax azimuth in different settings regarding fault and rock properties, stress regime, and fault structure. Again, our focus is the far-field perspective, i.e. at distances far beyond 100 m from the fault core (Fig. 1). Thus, this work does not aim to answer the question of to what extent the stress tensor components are affected in the near field.

2 Model set-up

2.1 Model concept

We set up generic 3-D models with model dimensions, rock properties, and an initial stress state that are like the ones from a 3-D geomechanical–numerical model of a potential siting area for a high-level radioactive waste disposal site in northern Switzerland, presented by Hergert et al. (2015). For implementation in the model, faults are represented by contact elements, which allow an offset, or 3-D elements, which are elastically or plastically weaker than the surrounding rocks (Fig. 1).

2.2 Partial differential equation and solution scheme

The two key components of a static stress state are a result of volume forces due to gravity and surface forces from plate tectonics. Neglecting acceleration, the resulting partial differential equation is the equilibrium of forces. For the upper crust assuming linear isotropic elasticity is a good approximation to describe the stress–strain relation (e.g. Tesauro et al., 2012). Thus, for simplicity the three key model parameters in our study are density (ρ), Young's modulus (E), and Poisson's ratio (ν). Additionally, the Mohr–Coulomb criterion, using the friction (μ) and the cohesion (C), will be used for some models. As we introduce a fault in our model with different techniques, we solve the problem numerically using the FEM.

2.3 Geometry and material properties

The reference model has an extent of 10 km in each horizontal direction and 3 km in the vertical direction (Fig. 3). The model is intersected in its entirety in the centre by a 60^∘ inclined fault, represented by cohesionless contact elements with a friction coefficient of μ=0.4 (friction angle ϕ=21.8^∘). The main shortening direction is perpendicular to the strike of the fault. Homogeneous linear elastic and isotropic material properties are assigned to the reference model, having a Young's modulus of E=15 GPa, a Poisson's ratio of ν=0.27 and a density of ρ=2550 kg m⁻³. The FE mesh for the reference model has a resolution of 50 m in the X and Z direction and 500 m in the Y direction. The mesh was created with HyperMesh 2017 and 2019, respectively; the solver used is Abaqus 6.14.1.

Figure 3 The model extent ($\mathrm{10}\times \mathrm{10}\times \mathrm{3}$ km³, in green) with the fault (blue plane) inclined by 60^∘ (dip angle); in red is the visualization path at a depth of 660 m along which the stress magnitudes are presented for the majority of the figures. The vertical blue line indicates the location of a virtual vertical borehole (Fig. 4). The displacement boundary conditions in purple represent 10 m shortening ($\mathit{ϵ}=-\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}$) in the X direction (perpendicular to the strike direction of the fault), which governs the S_Hmax magnitude, and 2 m of dilation ($\mathit{ϵ}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$) in the Y direction (parallel to the strike direction of the fault), which drives the S_hmin magnitude.

2.4 Model scenarios

The scope of the study is to investigate factors that affect the stress state in the broader vicinity of faults. These include the element resolution (pre-tests), the representation of the fault by contact elements with a variable friction coefficient, representation of the fault by elastic weaker elements or by elements with elasto-plastic rheology, the inclination of the fault, the strike direction relative to the shortening direction, the variation of the rock stiffness (Young's modulus), and the size of the fault and model itself. In order to allow good readability of the study, specific variations of the model are always briefly explained before presenting the modelling results.

2.5 Initial stress state and boundary conditions

We implement an initial stress state of the model that is in equilibrium with the gravitational forces without resulting in any significant displacement along the fault and the model geometry. We follow the technical procedure as explained in Hergert et al. (2015). In a second step we apply along the model lateral boundaries displacement boundary conditions that result in tectonic stresses throughout the model volume. The main shortening of the reference model is perpendicular to the fault (X direction) of the order of −10 m ($\mathit{ϵ}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}$), which then corresponds to the S_Hmax orientation. Parallel to the fault strike (Y direction), the model undergoes a slight dilation of 2 m ($\mathit{ϵ}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$), which is then the orientation of S_hmin (Fig. 3). The stress magnitudes resulting from the boundary conditions are shown in Fig. 4 along a vertical synthetic well path in the centre of the model. This stress state is in general agreement with stress magnitude data that were derived from a measurement campaign in northern Switzerland using >150 micro-hydraulic fracturing and sleeve reopening tests (Desroches et al., 2021).

Figure 4 Virtual vertical well path in the centre of the reference model. Shown are the resulting stress components, which are S_Hmax, S_hmin, S_V, and the von Mises stress. The fault with a friction coefficient of μ=0.4 is traversed at a depth of −660 m. When crossing the fault from the hanging wall to the footwall block, there is a sudden increase in S_V for S_Hmax and for S_hmin a little less.

2.6 Stress definition and visualization

The 3-D stress state of the Earth's crust is described by a second-rank tensor (σ, Jaeger et al., 2011) with nine components, but due to its symmetry only six components are independent from each other. As is common in geoscience, compressive stress magnitudes are positive and tensile stresses are negative. The stress state can also be described with the magnitudes and orientations of the three principal stresses. These principal stresses are named from the largest to the smallest as ${\mathit{\sigma }}_{\mathrm{1}}>{\mathit{\sigma }}_{\mathrm{2}}>{\mathit{\sigma }}_{\mathrm{3}}$.

Furthermore, in our model the vertical stress (S_V) is a principal stress (Eq. 1). As a result, the two other principal stresses are in the horizontal plane and are labelled as the minimum and maximum horizontal stresses (S_hmin and S_Hmax).

$$\begin{array}{}\text{(1)}& S_{\mathrm{V}}=\underset{\mathrm{0}}{\overset{z}{\int }}\mathit{\rho }gz\end{array}$$

The relative ratio of these three principal stresses defines the stress regime (Anderson, 1905, 1951):

Normal faulting stress regime	NF	$S_{\mathrm{V}}>S_{\mathrm{Hmax}}>S_{\mathrm{hmin}}$
Strike-slip stress regime	SS	$S_{\mathrm{Hmax}}>S_{\mathrm{V}}>S_{\mathrm{hmin}}$
Thrust faulting stress regime	TF	$S_{\mathrm{Hmax}}>S_{\mathrm{hmin}}>S_{\mathrm{V}}$.

Additionally, we use the differential stress (σ_D) and its 3-D equivalent, the von Mises stress (σ_vM; Mises, 1913), to visualize the stress state (Eqs. 2 and 3).

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathrm{D}}={\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{3}}\text{(3)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathrm{vM}}=\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{2}}\right)}^{\mathrm{2}}+{\left({\mathit{\sigma }}_{\mathrm{2}}-{\mathit{\sigma }}_{\mathrm{3}}\right)}^{\mathrm{2}}+{\left({\mathit{\sigma }}_{\mathrm{3}}-{\mathit{\sigma }}_{\mathrm{1}}\right)}^{\mathrm{2}}}\end{array}$$

The model results are presented here in the same way whenever possible. Both the stress components (S_Hmax, S_hmin, S_V) and the von Mises stress are used to visualize the influence of a fault on the stress state. The results of the models are plotted along a horizontal path at a depth of −660 m (Figs. 3 and 5). This path is always parallel to the main shortening direction (X), except for the models with a variable fault strike. The visualization extends from the footwall block at −3000 m through the fault at 0 to +3000 m in the hanging wall block.

Figure 5 Stress magnitude visualization of the reference model from −3000 m (footwall – left) to +3000 m (hanging wall – right) for a constant depth of −660 m. Linear elastic material properties and a friction coefficient of μ=0.4 are used for the fault at 0 m, represented by the vertical black line. The dashed lines in comparison represent results of a similar model without a fault.

2.7 Pre-test: mesh resolution

The impact of the mesh resolution and sufficiency is investigated by varying the mesh size using elastic material properties only, like the reference model. A mesh resolution of 1000, 500, 250, and 100 m in all directions is tested; a finer resolution has been used with an element size of 50 m in the main shortening and depth direction (X and Z), for which the resolution parallel to the fault (Y) is 500 m.

The model with the coarsest resolution (1000 m) provides stress magnitudes that deviate significantly from the other models (red line in Fig. 6). Even for the model with a resolution of 500 m (magenta line in Fig. 6), the deviation from the higher-resolution models at a distance greater than 1000 m is clearly visible. All finer-resolution models (≤250 m) have only small differences close to the fault (<1000 m; Fig. 6). This shows that all models with a resolution of 250 m and finer have a sufficient mesh resolution. A finer mesh is only useful if the stress changes close to the fault are of interest, which is not the case in this study.

Figure 6 The impact of the mesh resolution is compared. The coefficient of friction of the fault is μ=0.4 for all models. The coarse-resolution models, with 1000 m (red) as well as 500 m (magenta), show significant deviations from the reference model with a resolution of 50 m (black), while the models with a resolution of 100 and 250 m (green and blue) show only slight deviations close to the fault.

3 Results

3.1 Reference model

Within the reference model, the fault is represented by a contact surface (μ=0.4, C=0). As a result, the components of the reduced stress tensor increase in the footwall close to the fault and decrease in the hanging wall (Fig. 5). S_V and S_hmin rise to a similar level (+3 and +1 MPa) within the footwall block near the fault. An opposite behaviour is observed for the von Mises stress. S_Hmax, however, increases only slightly close to the fault (<1 MPa), which is the reason for the decrease in the von Mises stress near the fault. Corresponding to these changes, the stress magnitudes decrease next to the fault within the hanging wall block; the largest amount is for S_V, resulting in a slight increase in the von Mises stress. Significant stress changes of more than 1 MPa occur within a distance of 1000 m from the fault. The S_Hmax orientation is not affected by the fault.

The results of all other models presented subsequently are displayed on a horizontal path at the same depth. For the reference model, the stress variation around the fault for different depth ranges is also shown in Fig. 7. It remains unchanged that stress variations > 1 MPa are limited to a distance of about 1000 m from the fault. Relatively large variations can be seen at shallow depths (blue, −200 m) in contrast to greater depths (red, −2800 m). The general patterns of stress variation are similar, except for the vertical stress component. S_V is smaller in the footwall block close to the fault and larger in the hanging wall block at a depth of −2800 m, in contrast to observations at shallower depths (<2000 m). The reason is that S_V becomes σ₁ (normal faulting regime) for a depth greater than 2000 m, while at shallower depths a transition from a thrust faulting to a strike-slip regime occurs (Fig. 4).

Figure 7 Variations of the stress components (S_Hmax, S_hmin, S_V, and von Mises stress) at different depth levels are shown with respect to the distance of the fault. Stress magnitude changes are visualized along a vertical line at depths of −200, −660 (reference depth, used by the other figures), −1000, −1400, −1400, −2000, and −2800 m.

3.2 Friction coefficient

In geomechanics and seismology faults are usually parameterized using the friction coefficient and the cohesion (e.g. Morris et al., 1996; Di Toro et al., 2011; Röckel et al., 2022). Commonly, a friction coefficient between 0.6 and 0.85 is assumed (Byerlee, 1978) but examples exist of significantly smaller friction coefficients (Di Toro et al., 2011). However, to investigate the influence of the frictional properties of a fault based on the reference model, the friction coefficient is varied from very low (μ=0.1) to very large (μ>1).

Using a very large friction coefficient (μ>1), there is no visible influence by the fault on the stress magnitudes (Fig. 8), and the stress magnitudes are identical to a continuous mesh without a contact surface (dashed line in Fig. 5). In contrast, for a low-friction case (μ=0.1), stress variation is significant near the fault. The general pattern is similar as for the reference model, but the increase (footwall: +8 MPa) and decrease (hanging wall: −5 MPa) in S_V are much larger. Similar but not large S_hmin changes of +2 MPa are observed for the footwall and −2 MPa for the hanging wall block. The drop of S_Hmax in the hanging wall block is significant (−4 MPa), whereas the increase in the footwall block next to the fault is negligible. However, a S_Hmax decrease of about −1 MPa is visible in both the footwall and hanging wall block, even between 1000 and 3000 m away from the fault. This is a result of stress dissipation due to larger fault offset in the case of low friction. Variation of the von Mises stress is mainly driven by the variation of S_V. It is mostly σ₃, determined by the fact that S_V becomes significantly larger than S_hmin in the footwall block about 500 m next to the fault for the models with low-friction contact definition.

Figure 8 Impact of a variable friction coefficient on the stress state. Plotted are the S_Hmax, S_hmin, and S_V as well as the von Mises stress. The graph with the friction angle of μ=0.4 is the reference model (Fig. 5). Large stress variations near the fault are a result of low friction.

Overall comparison of the models with a different friction in Fig. 8 shows that the stress perturbations gradually decrease with an increase in the friction coefficient. A stress variation of >1 MPa is limited to a distance of ≈1 km, except for S_Hmax in the hanging wall block. None of the variations result in a visible change in the S_Hmax orientation, and it is always parallel to the maximum displacement (X direction).

3.3 3-D fault representation by elastic weak elements

The representation of a fault by a 2-D plane is not realistic for the immediate vicinity of the fault where a zone of damaged rock is expected. A more realistic approach seems to be the representation by a layer of elements with an elastic rheology of reduced stiffness (Fig. 9). This simulates the damage zone around the fault core (e.g. Faulkner et al., 2006).

Figure 9 Sketch visualizing the representation of the fault zone by elastically weak 3-D elements with a thickness of (a) 10 m made from three elements or (b) 30 m made of nine elements. The elements outside this fault zone are not visualized.

Herein, the fault zone has a width of 10 m represented by three elements normal to the fault (Fig. 9a). A Young's modulus of E=5, 1, and 0.25 GPa is tested, while the stiffer surroundings have E=15 GPa. The element resolution outside the fault area is 50 m in the X and Z direction and 500 m in the Y direction.

The stress magnitudes along the profile (Fig. 10) do not show a significant stress variability in the vicinity of the fault resulting from three of the less stiff elements. Stress changes are restricted to a very narrow domain, which are not visible; they are visually hidden behind the fault line. S_Hmax decreases depending on the decreasing stiffness. For the model with E=250 MPa fault representation, S_Hmax is always around 1 MPa lower because of stress dissipation by the low-stiffness fault domain. Therefore, the von Mises stress drops by the same amount. Stress dissipation also effects S_hmin, but by a much lower amount; for S_V no effect is visible.

Figure 10 Fault representation by a 10 m thin layer of three weak elements. The fault elements have a lower Young's modulus (E=5, 1, and 0.25 GPa) in contrast to the area outside this region (E=15 GPa). Shown in black is the reference model and vertically the implemented fault zone. Stress changes are narrowly limited to the area of the fault so that they are hidden by the visualization of the fault zone.

Another model version has a thicker fault of 30 m, represented by nine elements normal to the fault (Fig. 9b). Like the 10 m models, S_Hmax drops, especially for the model with the least stiff fault domain (E=250 MPa), by around −3 MPa (Fig. 11), again an effect of the stress dissipation. S_hmin decreases by almost 1 MPa, whereas S_V is stable. Near the fault, S_Hmax, S_hmin, and S_V decrease significantly, limited to a region narrow to the fault (<100 m). The von Mises stress variation is mainly driven by the reduction of S_Hmax because of the less stiff fault parts. There is no change in the S_Hmax orientation to observe and it remains parallel to the X direction.

Figure 11 Fault representation by 30 m (nine elements) of elastic weak elements having a lower Young's modulus (E=5, 1, and 0.250 GPa) compared to the area outside this region with E=15 GPa. In colours are the model results with the less stiff 3-D fault representation. Shown in black is the reference model using contact surfaces and vertically the implemented fault zone at 0 m.

3.4 3-D fault representation by elements with elasto-plastic rheology

As purely elastic elements do not allow failure, they cannot dissipate stresses such as a contact surface is able to do. To accommodate both the ability to dissipate stresses and the representation of a damage zone, elements with elasto-plastic rheology within the fault zone are now used. Out of a continuous mesh, elements close to the fault location were selected in a staircase-like manner, which have a specific plastic yield criterion. These fault elements laterally have a range of one (Fig. 12a) to eight elements (Fig. 12b). These elements have a friction angle of ϕ=30^∘ (friction coefficient μ=0.58) and a low cohesion of C=0.1 kPa. The dilation angle used is ψ=25^∘. In contrast to that, the non-fault elements have a much larger cohesion (C=500 kPa) but the same friction and dilation angle. The element resolution in the vicinity of the fault is 100 m in the X and Z direction and 500 m in the Y direction. The elastic material properties are the same as used by the reference model (E=15 GPa, ν=0.27, and ρ=2550 kg m⁻³).

Figure 12 Sketch showing the fault representation by selected elements out of the mesh, which plastify as a result of friction and a low cohesion of C=0.1 kPa. Elements outside this region (white area, mesh not shown) have a cohesion of C=500 kPa. A friction angle ϕ=30^∘ (friction coefficient μ=0.58) is used for all elements for the first test. Different numbers of lateral elements representing the fault are tested, ranging from one (a) to eight (b) lateral elements. As the element size is 100 m near the fault, the total width of the stairstep-like fault ranges from 100 to 800 m.

The representation by means of staircase-like elements with elasto-plastic properties (Fig. 13) shows that the impact on the stress components is nearly independent from the number of laterally used elements that allow plastification. S_Hmax, S_hmin, and S_V slightly increase in the footwall block near the fault domain and slightly decrease in the hanging wall block, again near the fault domain. The overall variation of S_Hmax, S_hmin, S_V, and the von Mises is <1 MPa. Stress magnitudes do not show any discontinuous behaviour at the fault zone, as the reference model do. Stress variations are restricted to a zone of about 1000 m next to the fault domain. Again, the S_Hmax orientation is not disturbed as a result of the fault.

Figure 13 Fault representation with staircase-like elements with elasto-plastic rheology (Fig. 12) that are allowed to deform non-elastically. Shown in black is the reference model with the implemented fault; in colours are the models with a continuous mesh with one (magenta) to eight lateral elements (dark blue). These elements have a low cohesion of C=0.1 kPa and a friction angle of ϕ=30^∘ (friction coefficient μ=0.58). The maximum width of eight elements is visualized by the dashed blue lines.

Figure 14 Fault representation with four staircase-like elements with elasto-plastic rheology (Fig. 12) that are allowed to deform non-elastically. Shown in black is the reference model and the fault centre (vertical); in colours are the models with a friction angle of ϕ=30, 25, 20, and 15^∘. The ϕ=30^∘ model is the same as the four-element model in Fig. 13. The width of four elements is visualized by the dashed blue lines.

The model having four weak elements laterally is used again to investigate the impact of the friction. Friction angles of ϕ=30, 25, 20, and 15^∘ are applied. The ϕ=30^∘ model has already been used for the variation of the number of lateral elements (Fig. 13: four elements). Modelling results in Fig. 14 show that a decreasing friction angle increases the stress variation near the fault. S_Hmax, S_hmin, and S_V increase in the footwall block near the fault, while a slight decrease can be seen in the hanging wall block. However, swing-in effects can be observed on both sides of the fault. The largest-magnitude changes are about +4.5 MPa for S_V, +5 MPa for S_Hmax, and +2.5 MPa for S_hmin. At a distance of >1400 m to the fault centre, the variation of the stresses is <1 MPa. The S_Hmax orientation remains unaffected.

3.5 Variation of the fault dip angle

To study the impact of the fault dip angle, several models with different fault inclination are prepared. These models have a dip angle of 30, 40, 50, 70, and 80^∘, in contrast to the reference model (60^∘, Figs. 3 and 5). Elastic material properties are the same as used in the reference model: E=15 GPa, ν=0.27, and ρ=2550 kg m⁻³ with fault representation by contact elements μ=0.4 and C=0.

In Fig. 15, it can be seen that the stress perturbation pattern is similar compared to the reference model. With increasing dip angle from 60^∘ (reference model) to 70 and 80^∘, the stress perturbation slightly decreases. The reduction is most significantly visible for the S_V magnitude in the footwall block next to the fault. Stress magnitudes at a distance from the fault increase slightly for the large-dip-angle models as the stress dissipation by the fault decreases.

Figure 15 Influence of the dip angle of the fault on the stress components S_V, S_Hmax, and S_hmin, as well as the von Mises stress. Shown are the models with a fault dip angle of 60^∘ (reference model), 70, and 80^∘. By increasing the dip angle, the magnitude of stress perturbation decreases.

Figure 16 Influence of the fault dip angle on the stress components. A range of fault dip angles is presented, including 60 (reference model), 50, 40, and 30^∘. By reducing the dip angle, the stress magnitude changes and the distance of the lateral stress perturbation increases. The most pronounced stress perturbation is seen for shallow dipping faults (30^∘ in red).

A decrease in dip angle of the fault results in a significantly more pronounced increase in the stress perturbation near the fault (Fig. 16). This results in an increase in S_Hmax by >4 MPa in the footwall block and a decrease of about −4 MPa in the hanging wall block using a fault inclination of 30^∘. An increase in the S_V and S_hmin magnitudes in the footwall block and a decrease in the hanging wall block are clearly visible. The influence of the fault on the S_Hmax magnitude on both the footwall and hanging wall block is at a distance to the fault of about 1500 and 2000 m. However, the large distance is an effect of the small fault dip; the real distance is half of the values for the 30^∘ model. There is no perturbation of the S_Hmax orientation.

3.6 Variation fault strike angle

In addition to the influence of the dip, the influence of the S_Hmax orientation with respect to the fault strike is investigated. Thus, a strike angle of 90^∘ as the reference model is compared with other models where the fault is striking with an angle of 75, 60, 45, 30, and 15^∘. To geometrically allow such strike angles, the models are extended in the X direction from 10 to 20, 30, and 50 km for the models with a fault strike of 45, 30, and 15^∘, respectively. The resulting boundary conditions are adjusted to ensure comparability.

Results of the strike angle variation (Fig. 17) are shown perpendicular to the strike direction of the fault. The impact of the fault strike variation on the S_Hmax and S_hmin magnitude is minimal. Clear deviations are only observed for S_V in the footwall block, where S_V is smaller compared to the reference model. As a result, the von Mises stress is also less variable in the footwall block next to the fault. The variation of the S_Hmax orientation varies with distance to the fault but does not exceed 1.5^∘, which is significantly smaller than the uncertainties of orientation data records (Heidbach et al., 2018). Therefore, no visualization of that is shown.

Figure 17 Stress components are shown for the models with a variation of the strike angle relative to the orientation of the maximum shortening using a friction coefficient of μ=0.4. In contrast to the reference model with a fault strike angle of 90^∘, the varied models have a strike angle of 75, 60, 45, 30, and 15^∘. The stress components are plotted perpendicular to the strike of the fault.

Since the models with the fault strike variation and the friction coefficient of 0.4 only cause small S_Hmax rotations, the influence of a lower friction (μ=0.1) is also investigated. The plot of the stress magnitudes (Fig. 18) shows a visible variation of the magnitudes for the different orientations of the fault. The general pattern is similar to the reference model. For S_Hmax, significant variations in stress magnitude are observed between the models due to stress dissipation resulting from low friction at the fault. The largest magnitudes are for the reference model (90^∘) as well as the 15^∘ model. In contrast, the 45, 30, and 60^∘ models have the largest S_hmin magnitudes. As a result of the largest variation of the S_Hmax magnitudes, the lowest von Mises stresses are observed for 45, 30, and 60^∘ models.

Figure 18 Variation of the strike angle, with 90 (reference), 75, 60, 45, 30, and 15^∘ relative to the orientation of the direction of maximum shortening. A friction coefficient of μ=0.1 is used in contrast to the similar models with a friction coefficient of μ=0.4 (Fig. 17). The stress components are plotted perpendicular to the strike of the fault.

For the first time in the model series, a significant variation in the orientation of S_Hmax is clearly visible with a fault strike variation using a friction coefficient of μ=0.1 (Fig. 19). The deviation of the orientation reaches up to about 14^∘ in the hanging wall block for the model with a fault strike of 30^∘, closely followed by the 45^∘ model. The S_Hmax rotation for the 30, 45, and 15^∘ models is clockwise, parallel to the strike of the fault, while in the models with a strike of the fault of 60^∘ as well as 75^∘, S_Hmax orientation is counterclockwise, i.e. tends to be perpendicular to the orientation of the fault. In the footwall block, the rotation of S_Hmax is also visible but less so than in the hanging wall block, with a maximum of about 6^∘.

Figure 19 Variation of the strike angle (75, 60, 45, 30, and 15^∘) relative to the orientation of the maximum shortening direction using a friction coefficient of 0.1. Shown are the variations of the S_Hmax orientation compared to the reference model with a fault strike angle of 90^∘ and a constant S_Hmax orientation of 0^∘. The angular variation is plotted perpendicular to the strike of the fault.

3.7 Young's modulus

Since the elastic material properties have a significant influence on the deformation on the rock on both sides of the fault, the Young's modulus of the host rock is varied. In addition to the Young's modulus of the reference model (E=15 GPa), stiffnesses of 5, 20, 30, 40, 60, 80, and 100 GPa are tested. In order to keep the model comparable, the boundary conditions were adapted (Table 2) so that the far-field stress magnitudes of the different models were equal.

Table 2 Boundary conditions are chosen depending on the Young's modulus to generate equal far-field stress magnitudes for the different models. The boundary conditions for 15 GPa are the reference model settings.

The variation of the Young's modulus has a limited effect on S_Hmax in the footwall block (Fig. 20), where in the hanging wall block S_Hmax decreases by up to −4 MPa with increasing Young's modulus next to the fault. S_hmin increases slightly with the Young's modulus in the footwall block and decreases in the same way in the hanging wall block slightly by up to −2 MPa. The S_V magnitude shows the same pattern, but the stress deviation is much larger near the fault: up to +7 MPa in the footwall and −4.5 MPa in the hanging wall block. The von Mises stresses decrease with increasing Young's modulus in the footwall block next to the fault and increase in the hanging wall block next to the fault.

Figure 20 The influence of Young's modulus on the stress perturbation is investigated. The models have a Young's modulus of E=5, 15 (reference model), 20, 30, 40, 60, 80, and 100 GPa.

In general, the stress perturbation increases due to a larger Young's modulus as stress dissipates on the fault. The lateral influence of the fault on the stress components, producing a stress variation of more than 1 MPa, is limited to a range from −1000 to +1000 m next to the fault. Again, the S_Hmax orientation is always parallel to the direction of principal shortening.

3.8 Model size

It is obvious that the influence of the fault on the stress state also depends on the size of the fault surface or on the overall size of the model. For this purpose, the size of the active fault surface using the reference model geometry is reduced to 4000×1000 and 4000×3000 m² (Fig. 21). Also, the reference model with the full fault surface is doubled and quadrupled in size. The resulting models then have dimensions of $\mathrm{20}\times \mathrm{20}\times \mathrm{6}$ and $\mathrm{40}\times \mathrm{40}\times \mathrm{12}$ km³, respectively. The resulting mesh resolution is then 100 and 200 m in the X and Z directions, respectively, and 1 and 2 km in the strike direction (Y) of the fault, which is parallel to S_hmin. The boundary conditions were adjusted accordingly to generate a similar stress state.

Figure 21 Model sketch with a reduced fault surface area of 4×3 km² (light and dark blue areas together) and 4×1 km² size (dark blue area only). Everything else is the same as shown by Fig. 3.

The comparison of the results in Fig. 22 shows that as the size of the fault increases, the magnitude deviation near the fault increases. Thus, in the hanging wall block S_Hmax is reduced by almost −3 MPa, while S_V in the footwall block increases by more than +5 MPa for the model with a side length of 40 km. As a result, the von Mises stress in the footwall block decreases more significantly close to the fault. However, the increase in the fault surface area does not have a significant influence on the far-field stress pattern. Significant stress changes (>1 MPa) occur up to about 1000 m from the fault. No rotation of the S_Hmax orientation can be observed.

Figure 22 Influence of the fault size on the stress components S_Hmax, S_hmin, and S_V, as well as the von Mises stress. Models with a reduced fault surface area with a size of 4×3 and 4×1 km² (Fig. 21), as well as models like the reference model with a total size of $\mathrm{20}\times \mathrm{20}\times \mathrm{6}$ km³ (double size) and $\mathrm{40}\times \mathrm{40}\times \mathrm{12}$ km³ (quadruple size), are shown.

3.9 Strain variation

The effect of stress anisotropy is studied by defining variable lateral boundary conditions. The shortening perpendicular to the fault strike (X direction) is tested from 1, 2, 3, 4, 6, 8, 10 (reference model), 12, 14, 16, and 20 m ($\mathit{ϵ}=-\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$ to $-\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$), where the dilation to the fault (Y direction) remains identical to the reference model of −2 m ($\mathit{ϵ}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$). Everything else is identical to the reference model.

The different S_Hmax magnitudes directly result from the variable shortening applied to the model boundaries (Fig. 23). The overall pattern is like the reference model. The observed variation is low for low strain, where variation is larger for higher strain. S_Hmax is smaller for larger strain away from the fault and increases a bit next to the fault. In the footwall block, the pattern is clear: the closer to the fault, the smaller the S_Hmax.

Figure 23 Influence of a variable strain on the stress components is shown. The models have a shortening of 1, 2, 3, 4, 6, 8, 10 (reference), 12, 14, 16, and 20 m ($\mathit{ϵ}=-\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$ to $-\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$) perpendicular to the strike of the fault (X direction) and a constant dilation of 2 m ($\mathit{ϵ}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$) parallel to the fault (Y direction). To avoid information overload in the legend, only the 1, 10, and 20 m models are indicated there. As the different lateral strains along the model boundaries result in different stress magnitudes, only the relative stress changes (local stresses − far-field stress) are shown for S_Hmax, S_hmin, and the von Mises stress. The general pattern of stress variation is like the reference model, and the variation is smaller for less strain and larger for more strain. However, relative variations of the stress components are not bigger than about 1.5 MPa for S_Hmax and S_hmin, around 2 MPa for the von Mises stress, and about 2.3 MPa for S_V.

The variation of S_hmin is similar to S_Hmax: variation is small for less shortening and increases with increasing shortening of the model (Fig. 23). S_hmin increases in the footwall block next to the fault and is smaller next to the fault in the hanging wall block.

Larger variation can be seen for S_V, with an increase in the footwall bock and a decrease in the hanging wall block. The S_V magnitude variation in the footwall block increases from +0.4 MPa for 2 m of shortening to +2.3 MPa for 20 m of shortening. Nearly the similar amount of decrease happens in the hanging wall block.

The von Mises stress variation (Fig. 23) increases with the increase in shortening compared to the reference model. For the model with little strain (<4 m) the observed variation of the von Mises stress displays another pattern. The von Mises stress increases in the footwall block and decreases in the hanging wall block next to the fault. Again, major stress variations are limited to a distance of less than 1000 m next to the fault. The S_Hmax orientation is not affected for larger shortening perpendicular to the fault. For the models with a shortening of 1 and 2 m in the X direction, the stress magnitudes are horizontally isotropic (S_hmin=S_Hmax) and the S_Hmax orientation is not clearly defined.

4 Discussion

4.1 Model set-up and assumptions

The goal is to investigate the impact of faults on the far-field stress state (>100 m). The model design does not allow estimations of the stress state or stress perturbations close to a fault (<100 m). Investigating that, a much finer mesh resolution would be needed. It is also questionable whether and which methods of fault implementation are suitable for this purpose.

Like all generic models, those used here are a significant simplification of rock physics, geological structures, and the fault representation itself. Except for two scenarios, only linear elastic material properties are used to represent the rock volume. This neglects various rheological processes within the Earth’s crust. But Hooke's law seems to be a proper approximation for the major mechanical behaviour of rocks in the upper crust, as the elastic thickness of the crust (T_e) is usually much larger than the models used here (Burov and Diament, 1995; Hyndman et al., 2009; Tesauro et al., 2012). According to field investigations by Maerten et al. (2016), most brittle deformation can be explained using linear elastic material properties. Furthermore, the focus is not on stress changes during the co-seismic phase (e.g. Lin et al., 2013; Brodsky et al., 2020; Shi et al., 2020; Zhang and Ma, 2021) or deformation over several seismic cycles. The focus is on the quasi-static stress state in the inter-seismic phase.

The reference geometry is a normal faulting structure with a fault dip of 60^∘, but the applied boundary conditions result in a thrust to strike-slip faulting regime at the depth, where stresses are plotted, usually at −660 m. Even if most models use specific structures and specific stress regime conditions, other structural settings or faulting regimes are covered by some of the models or specific result presentations. These are the variation of the dip angle (Figs. 15 and 16), the variation of the strain (Fig. 23), and the variation of the depth for the reference model (Fig. 7). Therefore, results for all stress regimes and faulting structures are provided. However, the overall behaviour remains unchanged.

The specific objective was to investigate how faults can lead to stress rotations since this has been claimed to be the reason for observed stress rotations on scales of tens of kilometres. However, for most scenarios only stress magnitudes are shown here. This is of course due to the fact that many models do not show S_Hmax rotations. Visualizing the stress magnitudes gives a much broader insight into the effect of faults on the stress state. And, if the stress magnitudes change, stress rotation is possible, but if the magnitudes do not change, rotation can be ruled out. Therefore, the stress magnitude visualization used also acts as a proxy for potential stress rotation.

To allow good comparability of modelling results, constant boundary conditions have been used, with a few exceptions. The models with different strain have different stress magnitudes as a result. For models having a different extent or a variable Young's modulus, the boundary conditions were scaled accordingly to ensure comparability. The models with a lower Young's modulus in the fault zone and low-friction contact faults dissipate localized stresses, which has not been corrected, as the influence on the result is small.

4.2 Discontinuity approach: contact elements

Several of the model scenarios use contact elements to represent a fault within the model. This is the case for the reference model, the variation of the friction, the fault dip and fault strike angle, the Young’s modulus variation in the host rock, the model size, and the boundary conditions. The overall observation is an increase in the stress components (S_Hmax, S_hmin, and S_V) in the footwall block and a decrease within the hanging wall block, both next to the fault (Fig. 24a–d). In contrast, the von Mises stresses decrease in the footwall block and increase in the hanging wall block. This is the case as S_V varies more than the other stress components.

Figure 24 Summary illustration of the results from various presented models. Panel (a) shows the impact of the fault friction (μ=0.1, 0.2, 0.4, 0.6, 1.0, and >1.0) using contact elements (Fig. 8) and the influence of the fault size and model size (Fig. 22). Panel (b) displays the influence of a variable Young's modulus of the host rock on the stress state near the fault and far from the fault (Fig. 20). Panel (c) shows the impact of a variable fault dip (Figs. 15 and 16), where (d) illustrates the impact of a variable fault strike and additionally friction variation (μ=0.1, 0.2, 0.3, and 0.4) on the stress state resulting from a fault represented by a contact surface (Figs. 17 and 18). The impact of a fault representation by 3-D elements is shown, where (e) elastically weak elements have a different stiffness (Figs. 10 and 11) and (f) where the elements are allowed to plastify as a result of a variable low friction (Fig. 14) and a laterally variable number of elements (Fig. 13).

For these contact surfaces, no cohesion (C=0) is used, which is nevertheless a reasonable and conservative simplification in particular for pre-existing faults or fault zones, as granular material has a very low cohesion: C<1 kPa (Schellart, 2000). On the other hand, cohesion strengthening can increase the cohesion to C>1 MPa (van den Ende and Niemeijer, 2019), C=8 MPa (Muhuri et al., 2003), or C=35 MPa for very high temperatures (Tenthorey and Cox, 2006). According to Tenthorey and Cox (2006), cohesion will reach 3 MPa for a 100-year earthquake recurrence interval at a depth of about 2 km.

The friction coefficients used for the contact surfaces reach from 0.1 over 0.4 (reference model) to 1.0 and larger. In the past, it was assumed that the friction coefficient of faults is about 0.6 to 0.85 (Brace and Kohlstedt, 1980; Byerlee, 1978; Brudy et al., 1997). But the friction can be much smaller if clay minerals dominate (Byerlee, 1978; Lockner et al., 2011), in the case of dynamic offset (Di Toro et al., 2011; Boulton et al., 2017), or for high pore pressures (Blanpied et al., 1992; Byerlee, 1993).

Low friction is also expected for large fault (zones) or subduction zones (Bird and Xianghong Kong, 1994; Carena and Moder, 2009; Iaffaldano, 2012; Fulton et al., 2013; Carpenter et al., 2015; Houston, 2015). The friction coefficient is of the order of 0.08 for the 2011 Tohoku–Oki earthquake (Fulton et al., 2013), μ=0.12–0.25 or 0.05–0.2 for the San Andreas Fault (Bird and Xianghong Kong, 1994; Carena and Moder, 2009), and μ=0 to 0.1 for tremors in general (Houston, 2015). Iaffaldano (2012) assumes a friction coefficient of 0.01 to 0.07 for large-scale plate boundaries. However, the investigated range of friction covers this variation well, except for μ<0.1.

As a free surface, or a fault with a very low friction coefficient, is unable to build up shear stresses (Hafner, 1951), principal stresses will be parallel and perpendicular to the surface (Hudson and Cooling, 1988; Osokina, 1988; Rawnsley et al., 1992; Petit and Mattauer, 1995; Bell, 1996; Camac and Hunt, 2009). A classic example is the San Andreas Fault (Mount and Suppe, 1987), where the interpretation of borehole breakouts and drilling-induced tensile fractures from nearby borehole indicates S_Hmax orientations that are almost perpendicular to the fault (Zoback et al., 1987; Mount and Suppe, 1992). However, the distance of these boreholes is >1000 m from the fault core in most cases and it is thus questionable whether the derived S_Hmax orientations can be used as an observable for the fault strength. Hickman and Zoback (2004) show in their analysis of borehole breakouts and drilling-induced tensile failures of the SAFOD borehole through the San Andreas Fault that significant S_Hmax rotations can only be resolved in the near field of the fault.

4.3 Continuity approach: weak elements as fault zone

4.3.1 Young's modulus variation in the fault zone

Fault representation by elastic weak elements exhibits no significant stress variation pattern using three elements (Fig. 10) compared to the reference model using contact elements. Even if the number of elements representing the fault zone is increased to nine (Fig. 11), the stress pattern is hardly different. Only close to the fault can a stress drop be observed for S_Hmax, S_hmin, and S_V. The von Mises stress increases locally, as the S_Hmax decrease is lower than for S_hmin and S_V. Localized swing-in effects can be observed from the extent, most probably an artefact of the mesh resolution.

Fault zones are a 3-D structure consisting of the fault core and the damage zone embedded within the host rock (Chester and Logan, 1986; Caine et al., 1996; Faulkner et al., 2003, 2006). Previous work suggests that the Young's modulus of the host rock decreases towards the damage zone, where the Poisson's ratio increases in the same way (Casey, 1980; Faulkner et al., 2006; Isaacs et al., 2008). However, the variation of the Poisson's ratio is not tested here. The observed reduction of Young's modulus is from 55.4 GPa down to 16.2 GPa (Isaacs et al., 2008) or a reduction of about 6.5 GPa, e.g. from 66 to 59.5 GPa (Faulkner et al., 2006). The range from E=15 to 0.25 GPa investigated here covers a large material property range. According Treffeisen and Henk (2020a), the amount of Young's modulus contrast have a strong impact on the resulting stress perturbation. Overall, the fault representation by means of elastically soft elements did not provide a stress pattern as the contact surface method did. It is probable that representing a fault using only elastic weak elements is a method of stress dissipation rather than an accurate representation of low-friction faults.

4.3.2 Friction variation within the 3-D elements

Models having a 3-D representation of the fault with a laterally variable number of elements are allowed to fail according to the Mohr–Coulomb criterion. The resulting stress state by a friction angle of ϕ=30^∘ and a cohesion of C=0.1 kPa did not show much difference (Fig. 13) compared to a model without a fault representation. Magnitude changes are of the order of less than 1 MPa next to the fault zone. The models with lower friction (ϕ=25, 20, and 15^∘) display larger stress perturbation in the vicinity of the fault (Fig. 14). The magnitude of stress perturbation is larger for the model using a friction angle of 15^∘ compared to the reference model with contact surfaces. The overall pattern is complex and some of the trends are similar, but the stress magnitudes are not decoupled when crossing the fault zone. As previously discussed, low friction can be assumed for present-day fault activity. However, resulting stress patterns differ from the results using contact elements. The continuous finite-element mesh does not allow a mechanical decoupling. This may be different for other methods such as DEM where resulting behaviour depends on the number of elements and the friction (Hunt et al., 2004).

4.3.3 Cohesion variation within the 3-D elements

Usually, the key driver between intact rock and the fault using the Mohr–Coulomb failure criterion is not the friction coefficient, but the cohesion. Even from the modelling perspective, cohesion has the largest impact (Treffeisen and Henk, 2020a) on the stress state. Therefore, models with elements that have elasto-plastic rheology employ the same friction (ϕ=30^∘, or μ=0.58), but a very low cohesion C=0.1 kPa within the fault zone, in contrast to C=500 kPa outside this area. This is also the case for elements with elasto-plastic rheology, even when the number of parallel elements reaches eight.

4.4 Distance of stress disturbance to faults

4.4.1 Far field vs. near field

We have not specified the exact distance for the far field or near field, as such a distance depends on the orientation, properties, and size of the fault as well as on given stress field in the surrounding model volume. Figure 1 and the previous content suggest that the far field is beyond about 100 m to the fault for intact host rock. As the ratio of displacement to fault length is about 1:100 (Torabi and Berg, 2011), even for a fault with a length of 10 km, the fault offset can be up to 100 m. Depending on the faulting type, a limited correlation between fault displacement and thickness of a damage zone can be observed (Childs et al., 2009; Torabi and Berg, 2011). But the thickness of the damage zone is limited to a maximum of several hundred metres (Faulkner et al., 2010; Savage and Brodsky, 2011). However, for faults with wide damage zones the impact of such a zone on the host rock is unlikely to be greater than for narrow fault zones using the distance from the damage zone as a measure.

The impact of the different modelling approaches on the stress state differs. But a significant stress perturbation is spatially limited to a maximum distance of 1000–2000 m next to the fault. Figure 24 provides a visual overview of modelling results. This major assumption is supported by several authors using different approaches from a more map-view perspective (Petit and Mattauer, 1995; Su and Stephansson, 1999; Provost and Houston, 2001; Yale, 2003; Faulkner et al., 2006). Also, observations from wells support the idea that the stress perturbation is usually <200 m away from the fault (Stephansson and Ångman, 1986; Barton and Zoback, 1994; Brudy et al., 1997; Tamagawa and Pollard, 2008; Lin et al., 2010). A rotation of about 90^∘ within less than 200 m in the vicinity of a fault has been observed near the Taiwan Chelungpu fault (Lin et al., 2010) and at the Lansjärv well (Sweden, Bjarnason et al., 1989).

Only models with an oblique fault orientation relative to the maximum compression can achieve significant S_Hmax rotation. Models with low friction (μ=0.1, Fig. 19) in particular show S_Hmax rotation of up to 14^∘ next to the fault. However, at a distance of 1500 m the deviation is smaller than 5^∘, which is quite below the uncertainties of the stress orientation indicators. Only when the friction coefficient becomes unrealistically small for faults in the inter-seismic phase (<0.1) can larger rotations be observed by the models at distances of >1500 m away from the fault.

The relative stress state affects the spatial stress perturbation (Pollard and Segall, 1987). Therefore, Yale (2003) assumes that in the case of low differential stress, the spatial extent of stress perturbation is able to be observed up to several kilometres away from the fault. This in general fits the results of the models varying the lateral strain, where the stress magnitude variation near the faults increases with a larger differential stress. Some previous models show more spacious far-field stress perturbations (Tommasi et al., 1995; Sánchez D et al., 1999; Camac and Hunt, 2009; Maerten et al., 2002), which are most probably an artefact of a mesh resolution that is too coarse.

4.4.2 Vertical rotation of the stress tensor

Usage of the reduced stress tensor (S_Hmax, S_hmin, and S_V) is based on the assumption that S_V is a principal stress. However, near a weak and non-vertical fault, the principal stress orientation will be vertically distracted, as principal stresses are always parallel to oblique to a free surface. This leads to a variation of all reduced stress components, including the shown S_V magnitudes. In the case of a thrust faulting or strike-slip regime, S_V will be smaller in the hanging wall block and larger in the footwall block next to the fault (e.g. Fig. 5). The opposite can be seen for a normal faulting regime, e.g. stress plots at greater depth (Fig. 7 at −2800 m).

4.5 Magnitude of stress perturbation

A decrease in horizontal stresses near the faults in the hanging wall and an increase in the footwall are reported for the Forsmark DBT 1 well (Sweden, Stephansson and Ångman, 1986). Fewer borehole breakouts in the hanging wall block and more in the footwall block are observed from the KTB well (Germany, Barton and Zoback, 1994). A reduction of σ₃ by about 5 MPa has been observed within less than 10 m near a tunnel at the Grimsel test site (Switzerland, Krietsch et al., 2019). All these observations fit the results of the models having a fault representation by contact elements, where the horizontal stresses are smaller above the fault (Fig. 25), and the horizontal differential stress is smaller in the hanging wall block (Fig. 26). The latter would make the occurrence of borehole breakouts less likely in the hanging wall.

Figure 25 Stress magnitudes from a virtual well section for the depth range of −400 to −900 m of the reference model having contact surfaces (continuous line) and a model without a fault (dotted line).

Figure 26 The von Mises stress and the difference between the two horizontal stresses (S_Hmax−S_hmin) are shown for the reference model with contact surfaces (continuous line) and a model with a continuous mesh (dotted line).

In contrast to that, larger horizontal stresses above a fault have been observed for the Lansjärv well (Sweden, Bjarnason et al., 1989). The maximum horizontal stresses are observed about 100 m above the fault in the hanging wall block, which also points to other causes. One possible explanation is the lithological variation in that well, where several pegmatites and amphibolites in that depth range have been observed (Bjarnason et al., 1989), which eventually provide larger magnitudes as a result of a larger Young's modulus.

According to Su and Stephansson (1999) the magnitude variation is positively correlated with the stress ratio and negatively correlated with the friction. This can be clearly confirmed by this study (Figs. 8 and 23), where the stress variation near faults is largest for low-friction models and models with a larger strain variation. Observations indicate that stresses decrease near a fault after an earthquake (Zhou et al., 2012; Wang et al., 2015; Li et al., 2023). This can be confirmed by the models for the hanging wall, but not for the footwall block. Either the observations are from the hanging wall block only, or other factors, like the 3-D structure, are responsible, which are not represented by the models used here.

4.6 Other potential factors

All models analyse the variation of stress components and the orientation towards generic models with only one homogeneous fault. The extent to which the results can be applied to other scenarios remains questionable. There are some scenarios where we assume that other factors could have a greater influence on the stress state. These include extensive settings such as horst and graben structures, listric faults, or step-over zones. In such cases, whole blocks may be completely decoupled, either by faults or by any kind of decoupling horizon (salt, wet clay, or pore over-pressure). The stress state in such a block is then dominated by gravity only. One potential example of this is the Arches National Park in Utah, USA, where the joints are almost perpendicular to the normal faults and are constant over several hundred metres (Kattenhorn et al., 2000). Secondary faulting also provides a possible explanation for the complex stress pattern within the Viking Graben (North Sea, Maerten et al., 2002). According to Siler (2023), large stress perturbations can be caused by a fault step-over structure in hydrothermal systems over a distance of more than 1000 m in the Great Basin, western United States.

Faults or fault zones in nature are never planar structures, as assumed by the presented models. Roughness plays a role, but the roughness in the direction of previous slip is much less than in other directions (Power et al., 1987). The geometrical complexity is a result of non-planarity (bending, listric, bifurcation), combination or coalesce of faults (step-over or relay zones), or other factors (e.g. Roche et al., 2021). Fault zones can exist out of several single parallel faults, which would probably produce a more widely distributed area of stress perturbation. Pore pressure, especially above hydrostatics, has a significant impact on effective fault normal stresses (Blanpied et al., 1992; Byerlee, 1993). Despite the large number of models presented, such complex structures or properties have not been tested.

Stress changes near the fault tip (e.g. horsetail fault terminations) lead to a complex stress pattern (Segall and Pollard, 1980; Rispoli, 1981; Homberg et al., 1997; Siler, 2023). To model that, using only linear elastic material properties would result in unrealistic local stress peaks as elastic energy would not be dissipated by plastic deformation. Therefore, such structures are not considered here. However, it can be assumed that stress changes induced by fault tips are negligible at distances of a few kilometres from the fault (Segall and Pollard, 1980; Su and Stephansson, 1999).

5 Conclusions

The results of our study show that the static fault friction coefficient, rock strength, stiffness, and density contrast of the fault significantly affect the stress tensor beyond the fault core. However, the stress magnitudes and stress tensor orientation are not significantly changed beyond a distance of about 1000 m. S_Hmax rotation is only observable when the overall orientation of S_Hmax is oblique to the fault strike and the static friction coefficient is low (e.g. μ=0.1). From these findings we can conclude that many of the stress tensor rotations that are documented in recent publications based on high-density data sets (Heidbach et al., 2007; Pierdominici and Heidbach, 2012; Rajabi et al., 2016, 2017b; Lund Snee and Zoback, 2018, 2020) are probably not controlled by faults. Other factors probably play a greater role, like variable rock property (e.g. Reiter, 2021) or the superposition of plate boundary forces with different orientation and magnitude (Ferreira et al., 1998; Rajabi et al., 2017a). Specific fault settings could also play a roll, like decoupled graben blocks (Ferreira et al., 1998; Rajabi et al., 2017a), secondary faults in extensional settings (Maerten et al., 2002), or fault termination or transfer zones (Siler, 2023). However, it is doubtful that their far-field effect extends beyond 10 km.

Appendix A: Symbols

C	Cohesion
DEM	Discrete-element method
E	Young's modulus
FDM	Finite-difference method
FEM	Finite-element method
FVM	Finite-volume method
g	Gravitational acceleration
SHmax	Maximum horizontal stress
Shmin	Minimum horizontal stress
SV	Vertical stress
X, Y, Z	Coordinates (Cartesian)
z	Depth
ϵ	Strain
μ	Static friction coefficient
ν	Poisson's ratio
ρ	Density
σ	Stress tensor
σ1	Largest principal stress
σ2	Intermediate principal stress
σ3	Least principal stress
σD	Differential stress
σvM	von Mises stress
ϕ	Friction angle
ψ	Dilation angle