Mechanical controls on horizontal stresses and fracture behaviour in layered rocks: A numerical sensitivity analysis

In layered materials, the deformation style, orientation, confinement, and 3D connectivity of natural fractures is generally impacted by changes in sedimentary facies and alternations in mechanical properties. In this study we address this effect and perform a numerical sensitivity analysis. Mechanical properties, confining pressures, and interfacial frictions are varied for a three-layered model, to investigate and quantify the relation between contrasting material properties, the principal horizontal stresses and fracture behaviour (i.e. deformation style and orientation). Firstly, the results show that tensile stresses develop in the stiffer layers due to the contrasting elastic parameters. The magnitude of these stresses is dependent on the ratio between the elastic parameters of stiffer and softer layers (i.e. Estiff=Esoft and νstiff=νsoft). There are no horizontal tensile stresses, when applying a compressive horizontal confining pressure (approx. 1/5 of the applied vertical stress). Implementing an interfacial friction lower than 0.2 will result in decoupling of the layers, resulting in slip on the layer boundaries and no tensile stresses within the stiffer layers. Further, the acquired numerical results are in good agreement with previously conducted laboratory work. Finally, we discuss whether the presented results can be used for better relating contrasting mechanical properties to potential fracture deformation styles and orientations in layered outcrops or subsurface reservoirs.


Introduction
As the available high permeable reservoirs for petroleum and geothermal exploitation are declining, exploration companies generally turn to less permeable, unconventional reservoirs. In these normally tight rocks, open natural fractures have a strong control on the effective permeability and fluid flow, and can therefore significantly enhance production (Gale et al., 2010;Olson et al., 2007;Toublanc et al., 2005). On the other hand, these fractures can result in unpredictable or even unwanted flow behaviour such as fluid flow channelling (Jolley et al., 2007;Toublanc et al., 2005;Wang et al., 2017). Therefore, multiple studies have argued that understanding the extent and 3D connectivity of naturally fractured systems is key in optimizing production from low permeable and structurally complex reservoirs (Jolley et al., 2007;Maerten et al., 2006;Odling et al., 1999;Toublanc et al., 2005;Wang et al., 2017).
Rocks are generally deposited in alternating lithological units with significant variation in mechanical properties. These properties generally define the tensile and compressive rock strength Smart et al., 2014), indicating under which stress state a rock fails (Roche et al., 2013). Field observations have shown that geological layering can also influence the mode, orientation, and abutment of natural fractures Ferrill and Morris, 2003;Roche et al., 2013). This implies that the natural fracture connectivity and geometry is at least partly dependent on geological layering and the coinciding alternating mechanical properties .
In field geology, contrasting or changes in mechanical properties are generally referred to as the Mechanical Stratigraphy (MS) (i.e., mechanical properties, layer thickness, and the nature of the layer interface) (Laubach et al., 2009). The effect of MS on fracture deformation style, orientation and confinement is widely observed in field examples. For instance, previous field studies on fractures in layered rocks have shown that the MS generally controls the fracture frequency and confinement in different layers (Cooke and Underwood, 2001;Hooker et al., 2013;Rijken and Cooke, 2001). In addition, other field observations show that MS also affects the fracture deformation style and orientation (Brenner and Gudmundsson, 2004;Ferrill et al., 2017Ferrill et al., , 2014Ferrill and Morris, 2003;Larsen and Gudmundsson, 2010;Smart et al., 2014). However, although the link between fracture behaviour and MS is clear, it should be noted that the MS can vary with diagenesis, and this may imply that the observed fracture stratigraphy may not always be representative of the MS which was present at the onset of fracturing (Laubach et al., 2009). This should be taken in to account, else erroneous interpretations can ensue (Laubach et al., 2009).
The change in fracture behaviour as a result of a mechanical stratigraphy is also addressed by multiple laboratory and numerical studies (Bourne, 2003;Douma et al., 2019;Guo et al., 2017;Sch€ opfer et al., 2011;Teufel and Clark, 1984;Warpinski et al., 1981). These studies showed that the change in fracture behaviour could be attributed to local horizontal principal stress changes which occur in response to alternating mechanical properties, since numerical models of the same experiments showed that horizontal tensile stresses are present in the stiffer layers (high E, low ν), whereas horizontal stresses are compressive in the softer layers (low E, high ν) (Bourne, 2003;Douma et al., 2019;Guo et al., 2017;Teufel and Clark, 1984).
Apart from contrasting mechanical parameters, laboratory and numerical studies also showed that the interfacial friction and applied confining pressures have a high impact on the horizontal stress distribution and fracture behaviour. For instance, studies addressing the impact of the friction between the layers showed that a low interfacial friction property resulted decoupling of the modelled layers, which minimized the effect of contrasting mechanical properties, resulting in no horizontal stress changes within the different layers (Bourne, 2003;Guo et al., 2017). Studies addressing the impact of confining pressures showed that the application of horizontal confining pressures resulted in elimination of the observed tensile stress and fracture containment in the weaker layers Douma et al., 2019;Kettermann and Urai, 2015;Nguyen et al., 2011;Ramsey and Chester, 2004).
The above stated field, laboratory, and numerical studies highlight that the relation between the mechanical stratigraphy, observed fracture behaviour and horizontal stress changes is well documented. However, these studies have yet not addressed the range of conditions at which this different stress and fracture behaviour occurs.
Therefore, in this study, we perform a numerical sensitivity analysis, which assesses systematically the impact of changing the contrast in elastic properties (i.e. E and ν), interfacial friction (μ) and confining pressures (P conf ) on the horizontal stress magnitude, for a wide range of material parameters and confining pressures. Our 3D three-layer modelling design and applied boundary conditions are based on a previously conducted laboratory study (Douma et al., 2019). For the sensitivity analysis, we first test the impact of contrasting mechanical properties for the different layers, by changing the Young's modulus and Poisson's ratio of each layer. Secondly, we test the impact of the interfacial friction by changing the friction parameter between the layer interfaces. Thirdly, the impact of the applied confining pressure is tested. Further, we discuss the numerical results and investigate the relation between the fracture behaviour and the implemented mechanical contrast. Finally, we discuss whether our results can be used for predicting fracture behaviour in the subsurface, or whether fracture observation made on outcrops can be used for assessing the contrast in elastic properties. Throughout this study, we define stiff layers as having a high Young's modulus (E) and low Poisson's ratio (ν), and soft layers as having the opposite mechanical properties. Further, fracture deformation style, orientation and confinement will be referred to as fracture behaviour.

Model design and material parameters
We used the Abaqus Finite Element Tool (SIMULIA Abaqus FEA®) for the 3D models. The modelling set-up and boundary conditions were largely based on the laboratory design presented by Douma et al. (2019) (Fig. 1a). In our numerical experiments, we implemented three different layers, with an adjustable width and height (Fig. 1a). We assumed linear elastic isotropic behaviour for each layer, so that the materials are defined by a constant and isotropic Young's modulus (E) and Poisson's ratio (ν), following Hooke's law for 3D stress. The elastic properties (E; ν) can be adjusted for each layer. Throughout this study, we used identical properties for the top and bottom layers (Table 1). Further, we used the geological sign convention, making tensile stresses or extensional strains negative and compressive stresses or contractional strains positive.
The top and bottom layer were fixed in the lateral directions (x, y) to prevent any lateral rotation or translation (Fig. 1a). Vertical translation and rotation were prevented by fixing the bottom layer in the z-direction (Fig. 1a). The four corner lines of each layer were fixed in the x or y direction, to ensure that the deformation remains isotropic and no lateral rotation of the centre layer occurs (Fig. 1a). Interfacial friction between the different layers follows Coulomb friction τ s ¼ μσ n , where: τ s is static shear stress, μ the friction coefficient, and σ n the stress acting normal to the interface, which in this case is the applied vertical stress (Abaqus 6.14, User Manual). Vertical and confining stresses (σ v ; P conf ) were applied to the top and side boundaries, respectively (Fig. 1a). Fig. 1b shows how the imposed boundary conditions prevent any lateral translation and rotation of the model and resulting in uniform lateral xdisplacement (note: deformation is exaggerated 40 times).
To compare the horizontal stresses for each layer, two sampling location were chosen. These locations were chosen so that they best account for the imposed boundary conditions and therefore best describe the horizontal stress behaviour of each layer (Fig. 1c). For the inner layer, the sampling cube was taken at the centre of the layer. To account for the boundary conditions assigned at the top and bottom of the model (i.e. no lateral movement (x, y ¼ 0)), the sampling cube for the outer layers was taken closer to the bounding interface ( Fig. 1a and c). For this study, the sampling cube covers multiple elements (n ¼ 144) and the stresses calculated for each element within the cube are averaged to get one representative number for the respective layer.

Performed numerical experiments
In order to test the impact of contrasting mechanical properties, we changed the ratio in elastic parameters between the inner and outer layer, for 32 unconfined pressure tests (Table 1). The ratio between the elastic parameters (i.e., E inner =E outer and ν inner =ν outer ) was changed by increasing or decreasing the Young's modulus and Poisson's ratio for the inner layer. The elastic parameters of the outer layers were kept constant (Table 1). Further, the different layers were assumed to be mechanically coupled (i.e. μ ¼ 1.0), which is a feasible assumption for layers experiencing high normal stresses. The impact of the interfacial friction was tested by changing the friction coefficient (μ) whilst keeping the mechanical properties constant. The impact of the confining pressure was tested by repeating the elastic parameter tests whilst applying confining pressure (P conf ) of 10 MPa (Table 1). For all models, the applied vertical stress was set at 50 MPa. To account for local stress deviations caused by the model aspect ratio (height/width), the width and height of the models were set equal for the numerical sensitivity analysis (Table 1). It should be noted that the model validation uses different dimensions because it was based on a previously conducted laboratory experiment (Table 1 and Fig. 2).
The impact of internal pore fluid pressures will not be addressed in this sensitivity analysis. However, the impact potential impact of implementing a pore fluid pressure on the horizontal stress distribution and fracture behaviour will be discussed in section 4.2.

Model validation
We validated the model set up by making a numerical realization of Q.D. Boersma et al. an experiment performed in the laboratory study performed by Douma et al. (2019), using identical boundary conditions and mechanical properties. For this experiment, the inner layer was relatively soft, whereas the outer layers were relatively stiff (i.e. E inner = E outer < 1, ν inner = ν outer > 1) (Fig. 2a). Further, no layer parallel shear was observed within the laboratory tests, which implies that the inner and outer layers were mechanically coupled. The results of the laboratory experiment indicate that mode II fracturing occurs in the inner layer, whereas the outer layers show mode I fracturing (Fig. 2a). Our numerical realization of the laboratory experiment, highlights that the differences in fracture behaviour in each layer can be explained by the differences in the modelled stress field, which are compressive within the softer inner layer and tensile at the layer interfaces of the top and bottom layers, respectively ( Fig. 2b and c).
The laboratory experiment (Fig. 2a) also highlights that the mode I fractures in the outer layers stopped propagating in close proximity to the lower and upper bounds, respectively. This is consistent with our numerical results, which show compressive stresses at the upper and lower bound of the outer layers as a result of the implemented boundary conditions (x, y and z ¼ 0 at the lower and upper bound) (Figs. 1a and 2b-c).
The horizontal stress differences between the three layers can be explained by the differences in horizontal deformation, which is depicted by the horizontal strain field (Fig. 2d). Following from Hooke's law, the inner layer shows more horizontal strain with respect to the two outer layers (i.e. higher Poisson's ratio (ν inner ) and lower Young's modulus (E inner )). Since the inner and outer layers are numerically coupled (μ ¼ 1.0), the strain difference between the three layers results   Boersma et al. in an extensional pull on the outer layers, resulting in tensile stresses (Fig. 2).

Local stresses as a function of contrasting elastic parameters
The modelling results ( Fig. 2) show that the horizontal stresses are dependent on the difference between the mechanical properties implemented within the inner and outer layers. This implies that changing these properties will result in different horizontal stress distributions. In this study, the mechanical contrast characterised by the ratio between the Young's modulus (YM) of the inner and outer, and the ratio between the Poisson's ratio (PR) of the inner and outer layers, respectively (i.e. E inner =E outer and ν inner =ν outer ). For example, the results depicted by Fig. 2 had a low ratio in YM ratio and a PR ratio of approximately 1, and this implied that the inner layer is softer with respect to the outer layers.
The numerical experiments testing the relation between a contrasting YM and the horizontal stresses for the indicate that a high YM ratio between the inner and outer layer (E inner =E outer > 1:0), results in tensile stresses within the inner layer and compressive stresses within the outer layers (Fig. 3a). A low ratio in YM (E inner =E outer < 1:0), results in tensile stresses in the outer layers and more compressive stresses in the inner layer (Fig. 3a). Further, the results show that the magnitude of the difference in the implemented YM defines the difference in modelled horizontal stresses, with a higher YM difference resulting in highest stress difference (Fig. 3a).
Changing the Poisson's Ratio (PR) of the inner and outer layers gives opposite results, with respect to changing the YM ratio ( Fig. 3a and b). A low ratio between the PR (i.e. ν inner =ν outer < 1:0) results in tensile stresses in the inner layer and compressive stresses in the outer layers ( Fig. 3b), whereas a high ratio (ν inner =ν outer > 1:0) gives compressive stress in the inner layer and tensile stress in the outer layer (Fig. 3b). Similarly, to was observed with changing the YM ratio, the magnitude in the stress difference is defined by the difference in the PR of the inner and outer layers (Fig. 3b).
By repeating this experiment for different ratios of the YM and PR (see Table 1), we can depict the relation between contrasting elastic properties and the horizontal stress field, for wide range of elastic properties (i.e., different rock layers). The horizontal stress field between the different experiments is linearly interpolated. Fig. 3c and d shows the resulting stress maps for the inner and outer layers, respectively, and highlight the elastic control on the modelled horizontal stress field in each layer. For example, these plots show that scenarios which have a low ratio in YM and high ratio in PR will results in compressive stresses in the inner layer (i.e. σ h � 0:2σ v ) and tensile stresses in the outer layers (i.e. σ h � À 0:2σ v ) ( Fig. 3a and b). Alternatively, implementing opposite YM and PR ratios (high YM ratio and low PR ratio), results in the opposite horizontal stress effect (Fig. 3 c and  d). Finally, it should be noted that the stress deviations away from purely elastic behaviour (i.e. σ h ¼ 0:0 MPa at E inner =E outer ¼ 1 and ν inner = ν outer ¼ 1) are caused by the imposed boundary conditions and chosen sampling location for the outer layers ( Figs. 1 and 3).

Impact of friction between layers
The numerical experiments testing the impact of changing the interfacial friction indicate that a low friction parameter (μ ¼ 0:1)

Fig. 2.
Validation of the numerical model, by comparing an interpreted micro-CT scan with our numerical results. Applied vertical stresses are set at 50 MPa and no confining pressures were applied a) Highlighted fractures on a x-z cross section through a micro-CT scan result taken from Douma et al. (2019). The implemented Young's modulus and Poisson's ratio are listed above the interpreted micro-CT scan. b) Horizontal stress results of the numerical analysis. c) Horizontal stress in the x-direction depicted as a vector field. d) Horizontal strain field in the x-direction. It should be noted that the horizontal principal stress and strain vectors show a radial pattern (equal in all directions) as a result of shear stresses acting on the layer interfaces and the isotropic material properties. Therefore, we have chosen to depict the stress field in the x-direction on a x-z cross section. Due to the high vertical stress, rotation in the z-direction is minimal. (Fig. 4a), results in decoupling of the three layers, preventing horizontal stress transfer between the different layers. This implies that the stress difference between the modelled layers is low (Fig. 4a). In addition, no tensile stresses are observed within the inner layer. As expected, increasing the implemented friction parameter results in a high stress difference between the inner and outer layers (Fig. 4b), with the inner layer having significant tensile stresses (σ h ¼ À 0:2σ v ). Following from the Coulomb friction law, layer decoupling and the coinciding prevention of horizontal stress transfer is dependent on the interplay between the interfacial friction coefficient and the applied normal stress (in our models the applied vertical stress). This implies that for our modelling set-up, the effect of slip between the interfaces on the modelled horizontal stress field only becomes apparent for relatively low friction parameters (μ � 0:2) (Fig. 4c), which is lower than generally observed on rock interfaces (0:2 � μ � 0:85) (Byerlee, 1978;Sch€ opfer et al., 2011).
Therefore, these results imply that most of the layers will be coupled for relatively high layer-normal stresses (i.e. buried rocks).

Impact of the confining pressure
As expected by the superposition principle (Jaeger et al., 2007), changing the applied confining pressure has a significant effect on the modelled horizontal stresses distribution. Our modelling results indicate that high confining pressures result in a complete elimination of the horizontal tensile stresses caused by contrasting material parameters ( Fig. 5a and b). Additionally, the numerical results also indicate that the applied confining pressures also reduces the strain difference between stiffer and softer layers. As was shown by Fig. 2, this difference in strain causes a tensile pull on the stiffer layer. Therefore, by mitigating this effect, confining pressures further prevent tensile stresses from forming within the stiffer layer (Fig. 5).
By applying a confining pressure of 10 MPa (i.e. P conf ¼ 0:2σ v ) to the models previously shown in Fig. 3, the effect of a horizontal confining pressure becomes especially apparent (Table 1) (Fig. 6). For all tested modelling scenarios, the applied confining pressures result in a complete mitigation of the tensile stresses caused by contrasting elastic properties within the inner layer (Fig. 6a). For the outer layers, tensile stresses are present when the contrast in elastic properties is extremely high (Fig. 6b). Therefore, these results indicate that the presence of tensile stresses in layered rocks are very sensitive to the presence of confining pressures (Figs. 5 and 6).  3. Results of the numerical sensitivity analysis assessing the effect of the ratio between the Young's moduli (E inner =E outer ) and ratio of the Poisson's ratio (ν inner = ν outer ), whilst applying no confining pressures and assuming that μ ¼ 1.0. Total number of models used for this analysis is 32. See Table 1 for the experimental boundary conditions. Normalization is done over the applied vertical stress. a) Normalized horizontal stress in the inner and outer layers as a function of the of the YM ratio. The ratio in PR is kept constant all models (ν inner =ν outer ¼ 1:0). b) Normalized horizontal stress in the inner and outer layers as a function of the of the PR ratio (E inner =E outer ¼ 1:0). c-d) Interpolated horizontal stresses as function of the ratios in YM and PR, for the inner layers (c) and outer layers (d).

Mechanical factors controlling fracture behaviour in layered materials
The effect that contrasting elastic properties, interfacial friction and confining pressures have on the observed fracture behaviour and horizontal stress distribution has been well documented and is addressed by multiple field (Brenner and Gudmundsson, 2004;Larsen and Gudmundsson, 2010;McGinnis et al., 2017), laboratory (Douma et al., 2019;Ramsey and Chester, 2004;Teufel and Clark, 1984) and numerical studies (Bourne, 2003;Douma et al., 2019;Guo et al., 2017). These studies showed that the change in fracture behaviour can be attributed to local stress changes, which are caused by differences in the horizontal deformation of the stiffer and softer layers (Bourne, 2003;Brenner and Gudmundsson, 2004;Douma et al., 2019;Ferrill et al., 2014;Guo et al., 2017;Teufel and Clark, 1984). Further, these studies showed how low interfacial friction and/or a high confining pressure resulted in layer decoupling and a complete removal of the effects caused by contrasting material properties (Bourne, 2003;Douma et al., 2019;Guo et al., 2017;Nguyen et al., 2011;Ramsey and Chester, 2004).
This study quantifies the effects caused by contrasting layer properties, interfacial friction and confining pressures on the horizontal stress distribution, for a wide range of elastic properties and applied boundary conditions. However, we did not relate these effects to a potential fracture behaviour. Therefore, for the purpose of relating our numerical results to potential fracture regimes, we propose that the stress results shown in Fig. 3 can be converted to potential modes of fracturing (i.e. mode I, mode II and hybrid fracturing). For this conversion, we assume classical fracture mechanics and Mohr-Coulomb failure (Fossen, 2010), so that the mode I fracture domain start when tensile stresses surpass À 0:075σ v . This 0.075 value comes from the average ratio between the Brittle Tensile Strength (BTS) and Uniaxial Compressive Strength (UCS), for most sedimentary rocks (Cai, 2010;Hoek and Brown, 1997). For the hybrid and mode II domain, we assume À 0:075σ v � σ h � 0:075σ v and σ h � 0:075σ v , respectively. The Fig. 4. Impact of the interfacial friction on the least principal horizontal stress. a) Horizontal stress field for a friction parameter (μ) of 0.1. b) Horizontal stress field for a rough interface (μ ¼ 1:0). c) Normalized horizontal stress vs the implemented friction parameter for the inner and outer layers, respectively. Normalization is done over the applied vertical stress.  5. Impact of the applied confining pressure. a) Horizontal stress distribution for a low confining pressure (P conf ¼ 0:1σ v ð5 MPaÞ). b) Horizontal stress distribution for a high confining pressure (P conf ¼ 0:4σ v ð20 MPaÞ). For these two models the applied YM and PR ratios were 2.5 and 0.33, respectively. interpreted fracture regimes as function of contrasting elastic properties for inner and outer layer are shown by Fig. 7a and b, respectively.
Using these converted plots, we can relate observations made in the laboratory to our numerical results (Fig. 7). Fig. 7c shows interpreted microCT scans of uniaxial compression tests with three layered samples having different mechanical properties (modified from Douma et al. (2019)). The first test (Sample A) shows a low elastic contrast between the inner and outer layers (i.e. E inner =E outer ¼ 1:044 and ν inner = ν outer ¼ 1:200), and shows hybrid/mode I fracturing for both layers (Fig. 7c). Fig. 7a and b also highlight that under these elastic ratios, hybrid fracturing is the dominant regime for both the inner and outer layers. The second and third tests (Samples E and F) show a high contrast in mechanical properties (i.e. the inner layers are stiff and outer layers are soft). Both interpretations of the experiments show dominantly mode II fracturing in the outer layers and mode I fracturing in the inner layer (Fig. 7a). Placing the two experiments in to the converted fracture regime plots ( Fig. 7a and b), indicates similar results. Here, for both tests, the inner layers fall within the mode I regime and the outer layers fall in the mode II regime. Finally, the fourth test (Sample D) has a softer inner layer and stiffer outer layers, and the interpretation shows mode I fracturing in the inner, and mode II fracturing in the outer layers, respectively (Fig. 7a). This is also highlighted by the fracture regime plots, with the outer layers plotting in mode I regime and the inner layer plotting in the mode II regime ( Fig. 7a and b).
These similarities imply that fractures can be used as clear indicators for estimating the ratios between contrasting material parameters. On the other hand, if the ratio between the elastic parameters is known, the modes of fracturing in the respective layers can be predicted. However, it should be noted that this predictability is only applicable for low effective confining pressures, since an applied confining pressure of 0:2σ v already results in a complete removal of the observed tensile stresses ( Fig. 6a and b).
Apart from indicating how contrasting material properties can control the fracture behaviour, our results also indicate how the applied confining pressure and interfacial friction can affect the modelled horizontal stresses Figs. 4, 5 and 6. As was shown here and by previous modelling and laboratory work (Douma et al., 2019;Guo et al., 2017;Kettermann and Urai, 2015;Nguyen et al., 2011;Ramsey and Chester, 2004), a high confining pressure and low interfacial friction can result in a significant reduction or even a complete removal of the horizontal tensile stresses Figs. 4, 5,6. Further, Douma et al. (2019) concluded that increasing the confining pressure resulted in a high probability of fracture confinement in the weaker layers. Therefore, by combining our results with these previous findings, we can derive expected fracture behaviour or abutment as a function of the contrasting layer properties, confining pressure and interfacial friction. This behaviour is summarized in Fig. 8.

Are the modelled tensile stresses high enough to cause tensile failure?
Using the numerical results, we were able to relate fracture behaviour to contrasting material parameters and other implemented boundary conditions (Figs. 7 and 8). However, equating the modelling stresses with measured BTS (Brittle Tensile Strength) magnitudes, shows that the modelled tensile stresses are lower than the implemented rock strength (i.e. σ h= BTS < 1.0) ( Table 2). This essentially implies that the created tensile stresses are insufficient to explain the observed mode I fractures in the stiffer layers (Fig. 7), and that additional processes are needed to describe the observed features.
Laboratory experiments performed by Douma et al. (2019) and numerical experiments performed by Guo et al. (2017) showed that in layered materials, fractures initiate within the weaker brittle layers, showing mode II behaviour. Further, these studies showed that these mode II fractures propagated in to the stiffer brittle layers showing mode I behaviour. This observation can in part explain why mode I failure occurs under σ h= BTS ratios lower than 1.0, since following from Linear Elastic Fracture Mechanics (LEFM) (Irwin, 1957) tensile stresses are localized in close proximity to crack tips. These additional local tensile stresses may result in σ h= BTS � 1:0 near mode II crack tips, resulting in either hybrid or mode I fracturing in the stiffer layers (Figs. 7 and 8). The presence of pore fluids can also explain why mode I failure occurs for relatively low σ h= BTS ratios. For instance, the presence of fluid overpressures within certain layers may result in the development of sufficient tensile stress for mode I failure to occur (Roche and van der Baan, 2015). Additionally, the presence of pore fluids within geological materials allows for sub critical crack growth to occur (Atkinson, 1984;Olson, 1993). Under these conditions, fractures can propagate under significantly lower tensile stresses than the implemented rock strength (Boersma et al., 2018;Ko and Kemeny, 2011;Nara and Kaneko, 2005;Nara et al., 2012;Olson, 2007).

Implications for subsurface reservoir studies
The impact that natural fractures can have on fluid flow and the effective permeability within tight reservoir rock is becoming increasingly recognized Gale et al., 2014;Smart et al., 2014), with the presence of fractures generally resulting in unpredictable or even unwanted fluid flow behaviour (Toublanc et al., 2005). Numerous outcrop, laboratory and numerical studies have shown that Fig. 6. a-b) Sample testing (ν inner =ν outer vs E inner =E outer ) results whilst applying 10 MPa confining pressure (i.e. P conf ¼ 0:2σ v ). See Fig. 3 for additional information.
Q.D. Boersma et al. the 3D geometry and 3D connectivity of these naturally fractured systems are in part controlled by layering and contrasting material parameters (i.e. Mechanical Stratigraphy) (Brenner and Gudmundsson, 2004;Cooke and Underwood, 2001;Douma et al., 2019;Ferrill et al., 2014Ferrill et al., , 2017Ferrill and Morris, 2003;Guo et al., 2017;Larsen and Gudmundsson, 2010;Rijken and Cooke, 2001;Smart et al., 2014;Teufel and Clark, 1984). Therefore, understanding how these parameters relate to potential fracture behaviour is key for making accurate predictions in layered rocks present in the subsurface (McGinnis et al., 2017).
The modelling results shown in this study quantified the relation between contrasting material properties, interfacial friction, confining pressures, the horizontal stress field and potential fracture regime, for layered materials (Figs. 2-8). We argue that these results can help to better predict fracture behaviour in layered subsurface reservoirs, and could for instance help in creating mechanically constrained 3D discrete fracture network descriptions..  . 7. a-b) Fracture regime as a function of contrasting mechanical properties, for the inner (a) and outer layer (b). These plots are created using the stress results depicted by Fig. 3b and c. For this figure, we assume P conf ¼ 0:0 and μ ¼ 1:0. See text for additional details. c) Fracture geometries in different layered experiments.
Interpretations and are based on Micro-CT scan results from Douma et al. (2019). Sample names were also taken from Douma et al. (2019). The ratios between Young's modulus and Poisson's ratio are highlighted on the figure. It should be noted that these ratios are based on the static Young's modulus and Poisson's ratio. The rock names, tensile and compressive rock strengths (Brittle Tensile Strength (BTS) and Unconfined Compressive Strength (UCS)) can be found in Table 2. To illustrate this argument, we have taken density and sonic velocity data of a layered reservoir rock from well K12-17, from the TNO geological survey well log database (NLOG) (data is freely available on: www.nlog.nl) (Fig. 9a-c). This well is located within the North Sea Basin (offshore The Netherlands), and targeted the layered Upper Slochteren Sandstone Formation (Permian).
By combining the well data with Figs. 3 and 7, the local horizontal stress field and fracture regime is estimated for the different mechanical units within the well bore. For this small modelling exercise, we assume that P conf ¼ 0.0 and that μ ¼ 1.0. We acknowledge that applying no confining pressures is a big assumption. However, multiple studies have shown and observed that high pore fluid pressures in porous rock can reduce the effective confining pressures, even under relatively deep conditions (Emery, 2016;Fisher et al., 2005;Gale et al., 2010;Imber et al., 2014;Warpinski et al., 2014).
The dynamic elastic properties ( Fig. 9d-g) are derived using the elastic wave equation for isotropic media, the density (ρ m ), P-wave velocity (V P ) and S-wave velocity (V S ) (Mavko et al., 2009). Here, the dynamic Young's modulus is calculated using equation (1): where the bulk and shear moduli (K dyn and u dyn ) can be expressed as: Poisson's ratio is calculated using equation (2): . 8. Impact of the tested parameters on the fracture behaviour and confinement. For this figure we assume that the inner layer is stiffer and stronger with respect to the outer layers. a) Impact of the assigned confining pressure. b) Impact of the implemented interfacial friction.  From these calculated elastic properties, 25 layers have been identified and for each of these 25 layers, an averaged YM and PR is calculated ( Fig. 9h and i) (See figure caption for additional details). Using these averaged elastic properties, the ratio between the YM and PR for each layer and the layer above is calculated (Fig. 9j and k). Finally, the calculated ratios in elastic properties are placed in Figs. 3 and 7, to infer to horizontal stresses and potential fracture behaviour (Fig. 9l).
The predicted horizontal stress and fracture behaviour highlights the impact of contrasting mechanical properties, with predicted fracture regimes mostly alternating between tensile hybrid fracturing (red) and compressive hybrid fracturing (blue) (Fig. 9l). The transition between the Upper Slochteren sandstone and Ameland claystone shows a sharp contrast in elastic properties and this transition therefore indicates relatively high tensile stresses and mode I fracturing (Fig. 9).
Although these results show a distinct impact of alternating layer properties (Fig. 9l), increasing the applied confining pressure results in a different estimated stress distribution (i.e. no tensile stresses and no mode I fracturing). Further, under subsurface conditions, these horizontal confining pressures are generally present (i.e. lithostatic and Poisson stresses), making the presence of tensile stresses in the subsurface highly unlikely (Fossen, 2010;Zoback, 2007).
While the above statements are true, alternating mechanical properties can also play role under normal confined conditions. For instance, a hydraulic stimulation case study done by Roche & van der Baan (2015) highlighted that as a result of contrasting layer properties significant horizontal stress differences occurred within the different layers. These authors showed that due to these stress differences, some layers where more prone to being hydraulically stimulated. Therefore, these and our results imply that properly accounting for changes in mechanical properties is key in identifying zones which either contain natural fractures or are prone to failure and hydraullic stimulation ( Fig. 9) (Roche and van der Baan, 2015).

Implications for field geology
While our results can be used for predicting fracture behaviour, outcrop descriptions generally provide fracture network geometries rather than the underlying physical mechanisms (Fig. 10) (Brenner and Gudmundsson, 2004;Ferrill et al., 2014;Larsen and Gudmundsson, 2010). Unfortunately, quantifying these mechanisms from a known fracture network geometry is difficult since, multiple conditions can result in similar fracture network geometries. Furthermore, because it can vary with time due to diagenesis and exhumation, the observed MS derived from outcrops may not be representative for the MS which was present at the moment fracturing (Laubach et al., 2009). However, if these mechanisms are properly taken in to account, we propose that the numerical results presented in this study, can be used to qualitatively estimate ratios in mechanical properties and applied boundaries from natural outcrop data. For instance, outcrops showing alternating layers of stronger and weaker intervals ( Fig. 10a and b), generally show mode II fracturing in the weaker layers and mode I/hybrid fracturing in the more competent layers (Brenner and Gudmundsson, 2004;Ferrill et al., 2014). Following from Fig. 8, this behaviour requires relatively high ratios in mechanical parameters, low effective confining pressures and that the alternating layers are mechanically coupled. Alternatively, outcrops showing more homogenous layer alternations, generally depict mode II failure for all layers (Fig. 10c) (Larsen and Gudmundsson, 2010), which following from our models implies low ratios in elastic properties. Furthermore, the interpretations highlighted by Fig. 10c, show that fracturing generally abuts against layer/stylolite interfaces, which could result from low interfacial frictions and/or normal stresses (Fig. 8). Finally, it should be noted that although our results can help in estimating mechanical properties and applied boundary conditions from fracture data, a full quantification of this correlation requires a different modelling workflow, at which fracture propagation in layered materials can directly be modelled.

Conclusions
This study investigates and quantifies the relation between different mechanical parameters, boundary conditions and different horizontal stress and fracturing regimes by performing numerical sensitivity analysis using FEM. The modelling design and implemented boundary conditions are based on previous deformation experiments on layered rocks conducted by Douma et al. (2019).
Firstly, the numerical results highlight the relation between elastic parameters and the modelled horizontal stresses, for the multi-layered experiments. The results show that a high ratio between the Young's modulus (i.e. E stiff =E soft > 1) and low ratio between the Poisson's ratio (i. e. ν stiff =ν soft < 1), results in tensile stresses in the respective stiffer layers.
Furthermore, our results indicate that the magnitude of the modelled tensile stress is dependent on the difference between the implemented mechanical parameters.
Secondly, the results show that the confining pressure and interfacial friction have a significant impact on the horizontal stresses. Here, increasing the applied confining pressure results in a significant reduction in the observed tensile stresses, with tensile stresses being completely removed for confining pressures reaching approximately 0:2σ v . Our results also indicate that a low interfacial friction (i.e. μ � 0:2) will result in a complete decoupling of the modelled layers, resulting in no tensile stresses in the stiffer layers. Further, we show that the presented numerical results can be converted into plots that relate the elastic parameters to a certain fracture regime. These plots are in good agreement with previously conducted laboratory experiments  (Douma et al., 2019).
Finally, we argue that the results can be used to predict fracture behaviour in layered rock, assuming that the elastic parameters are known. Alternatively, these results can help in estimating the elastic ratios between different layers from observations of different modes of fracturing.