Modeling of the in situ state of stress in elastic layered rock subject to stress and strain-driven tectonic forces

In this study we describe and compare eight different strategies to predict the depth variation of stress within a layered rock formation. This reveals the inherent uncertainties in stress prediction from elastic properties and stress 10 measurements, as well as the geologic implications of the different models. The predictive strategies are based on well log data and in some cases on in situ stress measurements, combined with the weight of the overburden rock, the pore pressure, the depth variation in rock properties, and tectonic effects. We contrast and compare stresses predicted purely using theoretical models with those constrained by in situ measurements. We also explore the role of the applied boundary conditions mimicking two fundamental models of tectonic effects, namely the stress or strain-driven models. In both 15 models layer to layer tectonic stress variations are added to initial predictions due to vertical variation in rock elasticity, consistent with natural observations, yet describing very different controlling mechanisms. Layer to layer stress variations are caused by either local elastic strain accommodation for the strain-driven model, or stress transfers for the stress-driven model. As a consequence, stress predictions can depend strongly on the implemented prediction philosophy and the underlying implicit and explicit assumptions, even for media with identical elastic parameters and stress measurements. 20 This implies that stress predictions have large uncertainties, even if local measurements and boundary conditions are honored.

Information on the in situ stresses, like the magnitude of the minimum principal stress or the orientation of maximum horizontal stress, may be assessed using different techniques such as extended leak-off tests, hydraulic fracturing treatments, borehole breakouts, earthquake source mechanisms, geological indicators, etc. (Terzaghi, 1962;Hast, 1967;Zoback, 2010;Zoback et al., 1985;Amadei and Stephansson, 1997;Zang and Stephansson, 2010;Schmitt et Published by Copernicus Publications on behalf of the European Geosciences Union. al., 2012;Reiter et al., 2014;Reiter and Heidbach, 2014).Maximum horizontal stress cannot be measured directly.In rocks composed of layers with different elastic properties, like in sedimentary or volcanic areas, layer-to-layer variations in horizontal stresses may arise.This phenomenon occurs in several lithology types such as coal, mudstone, siltstone, sandstone, limestone, shales, lava flows, intrusions, and pyroclastics (Haimson and Rummel, 1982;Warpinski et al., 1985;McLellan, 1987;Evans et al., 1989;Warpinski and Teufel, 1991;Cornet and Burlet, 1992;Gunzburger and Cornet, 2007;Cornet and Röckel, 2012).Such stress variations are well documented, but local information on the stress field is often sparse and incomplete; a continuously sampled stress profile is rarely available because extended leak-off tests often concentrate only on the target formations and are rarely published.In this case, analytic formulations (McGarr, 1988;Gunzburger and Cornet, 2007;Jaeger et al., 2009), or numerical models (Teufel and Clark, 1984;Gudmundsson, 2006;Roche and Van der Baan, 2015), have to be used to assess the stress variations.
This paper focuses on the 1-D depth variation in stress magnitude, and more specifically on the local variation in minimum horizontal stress in an extensional regime.
A common methodology to predict stresses in the crust beyond individual measurements consists of applying a stress update to a pre-chosen initial stress model.The stress update can represent processes such as tectonic effects, uplift or subsidence (Haxby and Turcotte, 1976;McGarr, 1988;Gunzburger and Cornet, 2007), temperature change (Voight and St. Pierre, 1974;Haxby and Turcotte, 1976;McGarr, 1988;Blanton and Olson, 1999), and viscous behavior (Gunzburger and Cornet, 2007;Cornet and Röckel, 2012).In most of the cases, the stress update depends on the local rock behavior and can be assessed using the elastic properties of the rock.Different models can be used as an initial state of stress.The simplest model, called the lithostatic model, assumes an isotropic state of stress equal to the vertical stress, as calculated from the weight of the overburden rock (Jaeger et al., 2009).Otherwise, in the uniaxial strain model, stress depends on the Poisson's ratio (Warpinski et al., 1985;McLellan, 1987;Savage et al., 1992;Addis et al., 1996;Jaeger et al., 2009).
An alternative prediction strategy is to assume that the state of stress in the crust is close to the maximum strength that rocks can support at the large scale (i.e., scales of a kilometer or more).The critical stress model then becomes applicable.In this model, depending on the frictional strength of pre-existing fractures and faults, the magnitude of the minimum and maximum principal stresses can be assessed if the magnitude of one of these is known (Zoback, 2010).Since the earth may not be critically stressed everywhere, this model can provide bounds for the minimum and maximum principal stresses (Brace and Kohlstedt, 1980;Townend and Zoback, 2000;Zoback, 2010;Konstantinovskaya at al., 2012;Meixner et al., 2014).The rationale for this ap-proach is that faults and fractures are the weakest structures in rocks.Therefore, in theory, the stress state cannot exceed the stress-state-inducing slip on an optimally oriented fault.Any excess stress will lead to rupture, resetting the stresses to below the critical stress.
The use of the critical model to obtain bounds for the upper and lower stresses was advocated by Brace and Kohlstedt (1980) since stress prediction from the elastic properties of rock is prone to large uncertainty, in particular for areas where the crust has a long and complex tectonic history, possibly resulting in nonelastic behavior.However, stress predictions based on elasticity are commonly used, notably because elastic properties are easier to obtain from well logs than strength properties and in situ stress measurements.In addition, numerous applications in seismotectonics, volcanotectonics, structural geology, and hydrocarbon exploitation require more detailed stress information than upper and lower bounds.
The objective of this work is to analyze different methodologies for 1-D stress prediction.We explore the roles of two different boundary conditions.They correspond either to a regional strain or stress perturbation, referred to as the strainand stress-driven models, respectively.In the first case, horizontal strain results from a displacement caused by plate tectonics (Savage et al., 1992;Blanton and Olson, 1999;Beaudoin et al., 2011;Song and Hareland, 2012;Reiter and Heidbach, 2014).In the second case, regional tectonic forces impose external stress instead of strain (Teufel and Clark, 1984;Bourne, 2003;Roche et al., 2013).Then, we assess the inherent uncertainty in stress predictions due to lack of information and/or the use of different predictive strategies.Eight different predictive strategies (Table 1) are carried out in order to explore the role of (1) the initial state of stress (uniaxial or lithostatic), (2) the implemented boundary conditions (strain-driven or stress-driven), and (3) the benchmarking and updating procedure (based entirely on the conceptual critical stress model or available in situ local stress measurements) in final predicted stress profiles.For reference, the predictive strategy using an initial uniaxial strain model, in situ stress measurements, and the strain-driven model is the model most commonly used in the literature (Savage et al., 1992;Blanton and Olson, 1999;Beaudoin et al., 2011;Song and Hareland, 2012).For illustration, all modeling strategies are applied to a field study case, where the resulting stress predictions are compared and contrasted.

Tectonic driving forces and stress corrections
For models depending on elastic rock behavior, prediction of the stress state is based on an initial stress model that is then modified by adding stress corrections (increments or decre- a and b indicate the regional tectonic strain perturbation and the regional tectonic stress perturbation, respectively. ments) to the horizontal stress components, that is where σ ' hl and σ ' Hl are the predictions for the effective horizontal stresses, σ hli and σ Hli are the local initial stress predictions, and σ hl and σ Hl are the local tectonic stress corrections.The subscripts "h" and "H" indicate the minimum and maximum horizontal stress, respectively.The initial stress predictions σ hli and σ Hli can originate from a uniaxial strain model, Eq. (A2) in Appendix A (Voight and St. Pierre, 1974;Haxby and Turcotte, 1976;Savage, 1992;Blanton and Olson, 1999) or a lithostatic model, Eq. (A1) (McGarr, 1988;Roche et al., 2013;Roche and Van der Baan, 2015).The initial model represents our best guess for the most representative baseline stress state.The stress corrections reflect the tectonic activity.
Various approaches exist to estimate the desired stress corrections σ hl and σ Hl .For instance, knowledge of the true magnitudes of the in situ local stresses σ hl and σ Hl permits their direct computation using where σ Hlm and σ hlm are local, independent measurements of both horizontal stresses.Alternatively, one could assume that the magnitude of the tectonic stress corrections σ hl and σ Hl increases until the initial stresses σ hli and σ Hli reach the critical stress model in the target formation (see Appendix A for background on the critical stress model).In this case, the stress corrections for the different stress regimes are obtained by combining Eq. (1) and Eq.(A3), yielding for a normal faulting regime, for a thrust-fault regime, and for a strike-slip regime, (5) The local critical-state stress corrections σ hl or σ Hl for the normal and thrust-fault regimes, Eqs. ( 3) and ( 4), respectively, may be obtained without any measurements.The local stress corrections applied to the intermediate principal stress, σ hl or σ Hl for the thrust and normal fault regimes, respectively, cannot be calculated in this case because the critical stress is based on the magnitude of the maximum and minimum principal stresses.For the strike-slip regime, both stress corrections σ hl and σ Hl are interdependent in Eq. ( 5).They can be calculated using Eq. ( 2) only if one is known from measurements.
It is important to distinguish here between the regional stress or strain corrections (boundary conditions) and the resulting local stress and strain variations (internal conditions).The local stress corrections, calculated with Eqs. ( 2), (3), (4), or (5), are defined for a specific layer.It is not possible to apply them to all other layers because in layered rocks stresses are prone to change as a function of depth since the stress corrections depend on the elastic rock properties.To calculate the local stress corrections for each layer, one needs to find the regional (external) tectonic perturbations applied to the rock formation that lead to the actual in situ local stresses in all layers.For convenience, we will refer to the regional stress or strain perturbations as the regional perturbations and the resulting local stress changes as the local stress corrections.
The regional perturbations are assessed according to strain-and stress-driven models.These two fundamental models are detailed in the next subsections.b) Strain-driven model with imposed regional strain perturbation that is positive (shortening), ε Hr , or negative (lengthening) ε hr .(c)-(f) Imposed regional stress perturbation that is positive (i.e., compressional), σ Hr , or negative (i.e., extensional), σ hr .(c)-(d) Non-coupled stress-driven model; (e)-(f) fully coupled stress-driven model.Each subfigure shows the boundary conditions (next to the bold horizontal arrows), a conceptual spring diagram, the resulting local stress and strain corrections (underneath the springs), and the resulting final local stresses in each layer.Subscripts "h" and "H" are the minimum and maximum principal stresses in the direction of the extension and compression, respectively.Subscripts c and s are the compliant and stiff layers.σ hls , σ Hls , σ hlc , and σ Hlc are the local stress corrections; ε hs , ε Hs , ε hc , and ε Hc are the local strain corrections resulting from the regional tectonic perturbations, σ hr , σ Hr , ε hr , and ε Hr .Absolute values of local stress and strain corrections are considered for the comparisons.Local stress equalization happens in the non-coupled stress-driven models (c, d).Both strain-driven and coupled stress-driven models can lead to similar local stress profiles, but with very different boundary conditions (a, b, e, f).

Strain-driven model
For the strain-driven model a biaxial strain model is appropriate in order to provide the magnitude of the local stress corrections since no variation in the overburden stress is assumed (Jaeger et al., 2009).In an isotropic elastic rock, the local stress corrections σ hl and σ Hl are then linked to the tectonic regional strain perturbations, ε hr and ε Hr , by (Savage, 1992;Economides and Nolte, 2000) and reciprocally (Jaeger et al., 2009) The local stress corrections σ hl and σ Hl depend on the vertical variation in Young's modulus E and Poisson's ratio υ.Horizontal stress corrections increase with increasing Young's modulus for regional shortening (Fig. 1a) and decrease for regional lengthening (Fig. 1b).Both horizontal regional strain perturbations ε Hr and ε hr are needed in order to assess the vertical variation in horizontal stresses.
Assuming an average Poisson's ratio of 0.3 or less, the stress correction in the perpendicular direction is merely 30 % of the one in the direction of the strain perturbation.Consequently, in cases where the regional strain correction occurs mainly in one direction, the second term in Eq. ( 6) is sometimes neglected, producing For instance, Eq. ( 8) may be used for a normal faulting regime where regional strain corresponds mainly to a lengthening in one horizontal direction.Equation ( 9) fits a thrustfault regime with a single direction of shortening.The regional strain perturbations, ε hr or ε Hr , may be locally derived from the horizontal stress corrections σ hl or σ Hl , respectively (Eq. 2), as obtained from in situ stress measurements.The problem becomes more complex if the second lateral regional strain perturbation is not negligible.Figure 1a and b conceptually illustrate the behavior of the strain-driven model.For such a model, the local strain corrections ε Hc and ε Hs , or ε hc and ε hs , are equal to the regional strain perturbations ε Hr and ε hr in each layer when a bilayer elastic medium is subjected to a compressive stress regime (a shortening in the direction of the maximum principal stress, i.e., positive ε Hr in Fig. 1a) or a normal faulting regime (a lengthening in the direction of the minimum principal stress, i.e., negative ε hr in Fig. 1b) since the two parts of the moving wall are vertically aligned.
Conversely, the local stress corrections σ Hl and σ hl change in each layer to account for uniform local strain perturbations due to different Young's moduli and Poisson's ratios.Consequently, the final stress profile shows layer-tolayer stress variations predominantly in the direction of the regional strain perturbations (i.e., σ H and σ h in Fig. 1a and b, respectively), and to a lesser extent in the orthogonal direction (i.e., σ h and σ H in Fig. 1a and b, respectively).
Figure 2a shows the final local stresses, σ 3 , and the local stress corrections, σ hl , calculated analytically with Eqs. ( 8) and ( 9), in various bilayer sections subjected to a regional extension.The increase in lithostatic stress with depth due to the weight of the overburden rock has been ignored in order to better highlight the effect of the layering.The local stress corrections σ hl and predicted stresses σ 3 depend solely on the local Young's moduli and Poisson's ratios of the individual layers; they are independent of the rock behavior or properties of the surrounding layers since no stress transfer or interaction occurs between layers.For instance, predicted stresses σ 3 are equal in the layers with a Young's modulus of 50 GPa, whether it is compliant or stiff compared with surrounding layers; predicted stress corrections are independent of the stiffness contrasts between layers.

Non-coupled stress-driven model
In the stress-driven model, the tectonic perturbations are driven by constant tectonic stresses, σ hr and σ Hr , applied to the rocks (Teufel and Clark, 1984;Mandl, 2000).In layered rocks, each layer then deforms independently if the layers are not coupled.The regional tectonic stresses induce local strain perturbations.For instance, the local strain correction for the non-coupled layers ε hnc in the direction of the minimum principal stress can be calculated using following Eq.( 7).The local strain corrections are discontinuous across layer boundaries due to the variation of the elastic properties, and slip may occur at the interfaces between layers.The local stress corrections, σ hl and σ Hl , are constant and equal to the regional stress perturbations, σ hr and σ Hr , in each layer, creating a homogeneous stress profile.This model is illustrated in Fig. 1a and d, in which the bilayer elastic medium from Fig. 1a and b are subjected to a regional stress that is either positive (compressional) in the direction of the maximum principal stress (i.e., σ Hr in Fig. 1c) or negative (extensional) in the direction of the minimum principal stress (i.e., σ hr in Fig. 1d).In both cases, the local stress and stress corrections σ h and σ hl in Fig. 1d, and σ H and σ Hl in Fig. 1c are constant and the resulting final stress profiles are independent of depth and the actual rock properties.

Coupled stress-driven model
If the layers are coupled together, slip cannot occur at the interface between the layers, and the strain has to be continuous throughout the medium.The regional stresses, σ hr and σ Hr , and the regional strain perturbations, ε hr or ε Hr , are constant across layer boundaries.The contrast in elastic properties then produces stress transfer between layers.Local stress corrections for individual layers thus depend on the elastic properties of all layers, contrary to the strain-driven model (Teufel and Clark, 1984;Mandl, 2000;Bourne, 2003).For a compressive regional stress perturbation σ Hr , a compliant layer is restrained from further shortening by any surrounding stiff layers, and the local stress correction σ Hl incorporates an additional layer-parallel extensional stress (Fig. 1e).Hence, the local stress σ H decreases in this layer and, in return, it increases in the stiff layer because of an additional layer-parallel compressive stress correction imposed by the shortening of the compliant layer.The opposite situation is encountered for an extensional regional stress perturbation, σ hr , due to the elongations of the stiff and compliant layers, requiring different local stresses to overcome the different Young's moduli (Fig. 1f).Notice that the coupled stress-driven model can lead to highly similar local stress profiles as for the strain-driven model (Fig. 1a and b), although the mechanisms are very different.
Numerical modeling is commonly used to compute the resulting stress profile for more complex media (Teufel and Clark, 1984;Bourne, 2003;Roche et al., 2013;Roche and Van der Baan, 2015).Appendix B contains the analytic formulation for this model for the simple case of two layers submitted to a tectonic force in a single direction.
Figure 2b and c show the profiles of the final local stresses, σ 3 , and the local stress corrections, σ hl , calculated with numerical modeling and the analytic solutions (Eqs.B11 and B12 in Appendix B), for the same media and initial lithostatic stress model as in Fig. 2a.The media are submitted to a negative (extensional) 10 MPa regional stress perturbation, σ hr .
For the numerical solution, the modeling is performed according to the same methodology as briefly described in Sect.3.4 and in detail in Roche and Van der Baan (2015).The local stress corrections σ hl reach extrema at the layer interface and tend to a constant value at the top and bottom, equal to the initial lithostatic stress minus the 10 MPa regional stress correction.The analytic formulation reproduces only the extreme values of the local stress corrections but not the local stress variation within a layer.For periodic media composed of two alternating layers, the local stresses are influenced by both surrounding layers.Consequently, the stress profile becomes almost constant within a layer (Roche et al., 2013), converging to the analytic solutions.
In the coupled stress-driven model, the magnitude of the local stress correction only depends on the contrast in Young's moduli, not their absolute values.The layer-to-layer local stress variation increases with increasing contrast in Young's moduli.For instance, in Fig. 2b, the local stress correction and the layer-to-layer stress variations are nearly constant for a ratio of 2 for different combinations of Young's moduli.They are significantly higher for a ratio of 5, despite similarity in Young's moduli in the different configurations.This mechanism is thus different from the straindriven model, where the local stress corrections depend only on the absolute values of the individual Young's moduli, and not their ratios (Fig. 2a).The stress transfer effect also depends on the thickness of the individual layers (see Eqs. B8 and B9), as confirmed by numerical studies (Roche et al., 2013), contrary to the strain-driven model.

Strategy
The four steps of calculation for the eight predictive strategies are presented in a synthesis view in Fig. 3 and are detailed in the next subsections (Table 1).The first step corresponds to the definition of the initial state of stress used as a base for the stress prediction.The second step is the definition of the calibration stresses.The third step incorporates the regional stress perturbations by combining the initial stresses and the calibrations (stress measurements or critical stress model).In the final step, the initial stresses are updated by including the tectonic effects, using either the stress-driven or strain-driven model.

Step 1: initial stresses
The stress predictions are based on two initial models with both horizontal stresses equal to either the horizontal stress of the uniaxial strain model or the lithostatic stress (Fig. 3).The average density is used to calculate the lithostatic stresses σ l , Eq. (A1).We use the lithology-dependent Poisson's ratios and pore pressures to calculate the uniaxial stresses σ u , Eq. (A2).For simplicity, we assume a Biot pore-pressure coefficient equal to 1, and the pore-pressure profile is calculated using a constant pore-pressure gradient.By comparing the stress predictions obtained with these two initial models, we explore the possible range of local stress corrections because the stress predictions obtained assuming an initial lithostatic and uniaxial state of stress are likely the upper and lower limits of the tectonic perturbations, respectively.

Step 2: preparation of the calibration stresses
To assess the magnitude of the tectonic effects, we compare the initial stress profile with a calibration stress model that is either the critical stress model or locally measured stresses.The choices for initial and calibration stresses produce four base predictions.These different choices for the calibration stress profiles allow us to compare stress predictions based entirely on conceptual models with those using in situ stress measurements.
In the case of observed data, the local stress corrections σ hl and σ Hl are directly calculated using Eq. ( 2).For the critical model, we must first assume a specific stress regime.For a normal faulting regime and a thrust faulting regime, the calculation of the critical horizontal stress is straightfor-ward since the vertical stress becomes the maximum principal stress or the minimum principal stress, respectively.For the strike-slip faulting regime, one needs at least one more data point to calibrate the stress in one horizontal direction at a specific depth.In theory, this can be done with in situ stress measurements of the minimum or the maximum principal stress, but in practice measurements of the maximum horizontal stress are subject to large uncertainty.In this study we assume a normal faulting regime.The critical horizontal stress σ 3c is calculated using the pore-pressure profile, an average coefficient of friction of 0.6, and the maximum principal critical stress σ 1c equal to the vertical stress σ v , Eq. (A3).

Step 3: computation of the tectonic perturbations (boundary conditions)
In order to compute the regional stress or strain perturbation, we compare the initial stress profiles to the calibration stress profiles.This comparison can be done in the directions of both horizontal principal stresses if the relevant calibration stresses are available.No stress measurements are available for the maximum horizontal principal stress in the following case study.In order to simplify the problem, avoid extra assumptions, and be able to compare stress predictions obtained with both calibration stresses independently, we postulate that only the minimum regional stress σ hr or strain ε hr perturbation exists and that the maximum regional stress σ Hr or strain ε Hr perturbation is negligible.
To get the regional stress perturbation σ hr , we calculate the average of the differences between the reference and initial stresses.The individual stress differences are given by Eq. ( 2) in case of stress measurements or Eq. ( 3) for the critical stress model.We use the average since the regional stress perturbation is assumed constant across all layers (see Sect. 2 and Fig. 1).For the regional strain perturbation ε hr , we convert all the stress differences into individual strain differences and subsequently compute their average, Eqs.(B4) and (B5).The stress-to-strain conversion thus assumes that the regional strain perturbation ε hr equals the average strain contribution of non-coupled layers.

Step 4: updating the initial stress to obtain the final stress profiles
The final step consists of modifying the initial stress profile obtained in step 1 in each layer as a function of the rock properties, depending on the regional tectonic stress or strain perturbation obtained in step 3, in order to get the final stress profile that takes tectonic effects into account.We focus on the minimum principal stresses, and the calculation of the intermediate horizontal principal stress is disregarded because the maximum and minimum principal stresses are more critical to assessing fracturing than the intermediate stress.Likewise, according to our assumption, i.e., a normal faulting regime and no tectonic perturbation in the direction of the intermediate principal stress, the magnitude of the layer-tolayer stress variation is potentially greater in the direction of the minimum principal stress.We explore both fundamental tectonic models, namely the coupled stress-driven and the strain-driven models, as applied to the four base predictions developed so far.We therefore obtain a total of eight different predictive strategies, leading potentially to eight different predictive stress profiles (Fig. 3).We use the regional stress σ hr and strain ε hr perturbations to obtain the final stress profiles, assuming respectively a stress-or strain-driven model (Table 1).For simplification, we assume that the regional stress σ hr and strain ε hr are constant along the layers.More complex models can be expected, for instance with increasing regional stress σ hr and strain ε hr .
For the strain-driven model, the local stress corrections σ hl are calculated analytically using the vertical variation of Poisson's ratio υ and Young's E with Eq. ( 8).Then, the depth-dependent local stress σ hl is calculated with Eq. ( 1), the pore-pressure P p , and the appropriate initial stresses σ hli .
For the stress-driven model, the local stress corrections σ hl and depth-dependent local effective stress σ ' hl are modeled using a discrete-element method (Cundall, 1988) assuming perfectly coupled layers.We use the same methodology as described in Roche et al. (2013) and Roche and Van der Baan (2015).In summary, we build a perfectly coupled 3-D elastic layered section that mimics the properties of the natural case (element size is 15 m).All the vertical changes in Young's moduli and Poisson's ratios as depicted in Fig. 4 are represented in the model.An initial state of stress corresponding to the chosen predictive base strategy (i.e., lithostatic stress σ l or uniaxial stress σ u ) is set within each element of the model.A boundary condition equal to the in situ stress plus constant regional tectonic stress perturbation σ hr is applied to the walls of the model in the direction of the minimum principal stress.Stress transfer across layers is properly accounted for once the model has reached equilibrium.
Finally, the critical stress model can be used to obtain upper and lower bounds for all predicted stresses (Brace and Kohlstedt, 1980).In the case study below, we assume a normal faulting regime such that only a lower stress bound is required since the maximum stresses are determined by the overburden weight.We will also assume a uniform coefficient of friction of 0.6.

Study area
The field example is located northeast of Fort Nelson, British Columbia, in the Western Canada Sedimentary Basin (Mossop and Shesten, 1994).The hydrocarbon reservoir (gas) is stimulated by hydraulic fracturing treatments.The studied section includes the Evie, Muskwa, and Otter Park low-porosity shale members of the Horn River Formation (Fig. 4a).They lie above the carbonatic Keg River Formation and under the shale of the Fort Simpson Formation (Curtis et al., 2010;Chalmers et al., 2012).The shales are overpressurized with a pore-pressure gradient of 11 to 16 MPa km −1 (Hurd and Zoback, 2012).
The rock densities and the compressional and shear wave velocities are obtained from well data (Fig. 4b).The depth dependence of the dynamic elastic parameters, Young's modulus E d , and Poisson's ratio υ d are calculated from the density ρ and the compressional and shear wave velocities, V p and V s (Fig. 4d, e).The methodology is provided in Roche and Van der Baan (2015).
For practical purposes, in rocks mainly composed of shale like in this study, there is no difference between the dynamic and static Poisson's ratio (Mullen et al., 2007).The static Young's modulus E s is obtained from the dynamic moduli E D measured in the well logs using (Mullen et al., 2007) where ϕ is the total porosity of the rocks.For the porosity of the rocks, we use an average value equal to 2 %, although it varies slightly between formations.The shale rocks of the Horn River and Fort Simpson formations are anisotropic with transverse isotropic properties (Khan et al., 2011).This anisotropic effect is disregarded here in order to simplify calculations.Also, viscoelasticity likely occurs in some layers, but this behavior is not considered.In practice, viscoelasticity implies that the long-term asymptotic solution can be evaluated through a linear model (see discussion in Roche et al., 2013).In this case the elastic constants evaluated in this study are likely overestimated, especially in the shale layers.

Stress regime
The Western Canada Sedimentary Basin is stressed tectonically because of the Rocky Mountains located to the southwest, which results in a difference in the horizontal stresses (Beaudoin et al., 2011;Reiter et al., 2014).The orientation of the maximum horizontal stress is well defined all over the basin and trends NW-SW (Bell and Gough, 1979;Bell and Bachu, 2003;Reiter et al., 2014), which is normal for the Rocky Mountain trench and folding.However, the tectonic regime is not well delineated by stress measurements (Reiter et al., 2014).Studies highlight spatial variation in stress regime, with thrust faulting in the foothills, strike slip within the basin, and a normal faulting regime in the eastern part of the basin (Bell and Gough, 1979;Bell and Bachu, 2003;Reiter and Heidbach, 2014).Likewise, depth variations in stress regime are described, with thrust faulting at shallow depth (i.e., < 350-600 m), strike slip at intermediate depth 500-2500 m, and normal faulting at greater depths > 2500m (Fordjor et al., 1983;McLellan, 1987;Reiter and Heidbach, 2014).In most of the measurements performed in the basin, the minimum prin- cipal stresses are lower than the vertical stresses calculated from the weight of the overburden rock (Bell and Grasby, 2012).This precludes a thrust-fault regime.In the whole basin, the minimum principal stress gradient ranges from 12 to 27 MPa km −1 , with an average of 18 MPa km −1 (Bell and Grasby, 2012).The SHmax / Shmin ratio is about 1.3-1.6 (Fordjor et al., 1983).
More locally in the studied area, the bottomhole instantaneous shut-in pressures do not show consistent variation from the toe to the heel, with minimum principal stresses equal to 26 ± 7 and 38 ± 8 MPa in the Muskwa and the Evie members, respectively (Fig. 4).This corresponds to a minimum principal stress gradient equal to 13 MPa km −1 on average, which is low compared to values measured at a nearby site (i.e., ≈ 21 MPa km −1 ; Hurd and Zoback, 2012), or in the Western Canada Sedimentary Basin (Bell and Grasby, 2012).However, it is consistent with the stress field map provided by Bell and Grasby (2012).
The minimum principal stress is thus lower than the vertical stress.In addition, the ratio of the minimum principal stress and the vertical stress is very close to the ratio for a critical state of stress.Thus, it is very unlikely that the maximum horizontal stress is higher than the vertical stress because it would involve an overcritical state of stress.Our data therefore indicate that a normal or transtensional stress regime is more likely than a strike-slip regime.This may be due to re-gional strike-slip faults that induce significant stress perturbations responsible for local normal faulting.

Initial stresses and calibration stresses (steps 1 and
2) The uniaxial stress σ u and the lithostatic stress σ l are computed using the methodology described in Sect.3.1.Both stresses increase with depth due to the weight of the overburden rock (Fig. 5a).The increase is lower for the uniaxial stress σ u than for the lithostatic stress σ l .The lithostatic stress σ l increases linearly with depth, whereas layer-to-layer variations in uniaxial stress σ u occur due to the vertical variation in the Poisson's ratio υ.
For the calibration stresses (see Sect. 3.2), the minimum horizontal critical stress σ 3c increases linearly with depth due to the weight of the overburden rock (Fig. 5a).The increase is lower than for the lithostatic stress.The three values of minimum principal stress, as measured in the Muskwa and the Evie members, are indicated in Fig. 5.They are set to be equal to the average of the instantaneous shut-in pressures for the different stages along the well.The resulting value is generally considered to be a proxy for the minimum principal stress (Zoback, 2010).
www.solid-earth.net/8/479/2017/Solid Earth, 8, 479-498, 2017 (b-g) Stress predictions for the various predictive strategies, including the local stress corrections σ hl and the predicted minimum principal stress σ hl .Circled numbers are labels of the predictive strategies as presented in Table 1.The relevant initial stress (i.e., uniaxial σ u or lithostatic σ l ) and calibration stresses (i.e., in situ stress measurements and the critical stress σ 3c ) are also represented for each predictive strategy.The diamonds represent the average values of the minimum principal stress (set equal to the instantaneous shut-in pressure values).(e) and (i) Show simplified stratigraphic column of the studied area with the same legend as for Fig. 4. Cases 3 and 7 in Table 2 are not represented because the tectonic effects are negligible, and therefore the predicted minimum principal stress σ hl is equal to the calibration stress, σ u .The Muskwa (M.), Otter Park (O.P.) and Evie (E.) members together form the Horn River Formation, surrounded by the Keg River (K.R.), Fort Simpson (F.S.), and Red Knife (R.K) formations.

Stress and strain regional perturbations
By comparing the initial and calibration values of stress, we obtain the different regional tectonic perturbations ε hr and σ hr (Sect.3.3).The uniaxial stress σ u is higher than the in situ measurements (Fig. a).This corresponds to negative local stress adjustments that are equal to −14 and −2 MPa in the Muskwa and the Evie members, respectively.This involves extensional stress and strain regional tectonic pertur-bations σ hr and ε hr equal to −6 MPa and −0.17 mm m −1 (the minus sign indicates extension).These are the basis for predictive strategies 5 and 1 in Table 1.
The lithostatic stress σ l is also higher than the in situ measurements (Fig. 5a).Since the lithostatic stress σ l is higher than the uniaxial stress σ u , the resulting stress and strain regional perturbations σ hr and ε hr are higher for the predictive strategies based on the lithostatic stress and they equal −34 MPa and −0.94 mm m −1 .The lithostatic stress σ l is higher than the minimum horizontal critical stress σ 3c .The resulting stress and strain regional perturbations σ hr and ε hr are equal to −21 MPa and −0.75 mm m −1 , respectively.These are the basis for predictive strategies 6 and 2 in Table 1.
The uniaxial stress σ u is very close to or slightly higher than the critical stress σ 3c in some formations.The resulting stress and strain regional perturbations σ hr and ε hr are thus very small, i.e., 1 MPa and 0.006 mm m −1 , and compressive (predictive strategies 7 and 3 in Table 1).
For nearly all the predictive strategies, the regional tectonic perturbations σ hr and ε hr are extensional.The exception is the predictive strategy that is based on the uniaxial and the critical models, where the tectonic perturbation is compressive.However, this perturbation is so low that its effect is negligible.This case is not illustrated in Fig. 5.These results show that a normal faulting regime does not necessarily involve extensional tectonic perturbations.A compressive perturbation can also result in a normal faulting regime if it is applied to a low-magnitude initial state of stress, like in the case of a uniaxial state of stress.The effect of the various regional tectonic perturbations on the final predicted stresses is presented in the next subsection.

Local stress corrections and predicted stress profiles
Next we compute the local stress corrections σ hl and the final stress profiles (see Sect. 3.4).In this section we first present the results obtained for the strain-driven model, then for the stress-driven model.In the next subsection we address the observed differences between both models.The local stress corrections σ hl are negative for an extensional regional tectonic strain or stress perturbation ε hr and σ hr (Fig. 5).For a strain-driven model, the local stress corrections σ hl are positively correlated to the Young's modulus E. Their absolute values are therefore greater in the Keg River Formation, followed by the Horn River Formation and finally the Fort Simpson Formation (Fig. 5b, c, and d).The magnitude of the local stress corrections σ hl in a specific layer also depends on the magnitude of the regional strain perturbations ε hr .They are maximal for strategy 2 (Table 1) based on the lithostatic model and the in situ stress measurements.
The predicted minimum principal stress σ hl depends on the initial state of stress and the local stress corrections σ hl .Because the local stress corrections σ hl are negative, the predicted minimum stress σ hl decreases due to inclusion of the tectonic effects.For all prediction strategies, due to the variation in Young's modulus E, the decrease is most important in the Keg River Formation, followed by the Horn River Formation, and finally the Fort Simpson Formation (Fig. 5b,  c, and d).Likewise, the decrease is slightly lower in the Otter Park Member than in the surrounding Muskwa and Evie members.Consequently, in regard to the Muskwa and Evie target members, there is a negative change in the minimum horizontal stress at the lower interface with the Keg River Formation, and a positive change at the upper interface with the Fort Simpson Formation.The intermediate Otter Park Member exhibits a higher stress than the surrounding target formations (Fig. 5b, c, and d).
The predicted minimum stress σ hl also depends on the magnitude of the regional strain perturbations ε hr .The changes at formation interfaces are maximal for strategy 2 (Table 1) based on the lithostatic model and the in situ stress measurements.In this case, the layer-to-layer stress variation reaches up to 13 MPa between the Horn River and the Keg River formations, 10 MPa between the Horn River and Fort Simpson formations, and 3 MPa between the Otter Park member and the surrounding Muskwa and Evie members (Fig. 5d).
For the predictive strategies 1 and 3, based on the uniaxial strain model (Table 1), layer-to-layer stress variations are initially present before the tectonic perturbation due to the vertical variation in the Poisson's ratio (Fig. 5b).Depending on the rock properties, the local stress corrections σ hl may www.solid-earth.net/8/479/2017/Solid Earth, 8, 479-498, 2017 promote the initial layer-to-layer stress variations, e.g., between the Horn River Formation and the Fort Simpson Formation.In other instances the local stress corrections σ hl may inhibit them, e.g., between the Horn River Formation and the Keg River Formation.Depending on the magnitude of the local stress corrections σ hl , the layer-to-layer stress variations in this initial stress profile may be decreased or enhanced by the computed local stress corrections.Nevertheless, in most of the models, the local stress corrections σ hl are higher than the initial variations.This implies that the tectonic effects can outweigh or reverse any initial stress variations due to the Poisson's effect.
For the stress-driven model, similar local stress corrections σ hl and minimum horizontal stress σ hl trends and magnitudes are obtained as for the strain-driven model (compare Fig. 5b-d with Fig. 5f-h and Fig. 6).Similar to the straindriven model, for the Muskwa and Evie target members, the minimum horizontal stress decreases at the interface with the Keg River Formation and increases at the interface with the Fort Simpson Formation, and the Otter Park Member exhibits a higher stress (Fig. 5f-h).
The predicted magnitudes of the layer-to-layer stress variations are also very similar.The maximum values are obtained for prediction strategy 8 (Table 1), where layer stress variations reach 11 MPa between the Horn River and Keg River formations, 8 MPa between the Horn River and Fort Simpson formations, and 5 MPa between the Otter Park Member and the surrounding Muskwa or Evie members (light gray curves in Fig. 5c).In this case, however, predicted stress changes are due to stress transfer between layers and depend on the contrast in Young's moduli between layers and their thicknesses (see Appendix B).This causes important differences between the stresses predicted with the stress-or strain-driven models.This is explored in the next subsection.

Comparison of stresses predicted with the strainand stress-driven models
Almost no regional tectonic perturbations, σ hr and ε hr , occur for prediction strategies 3 and 7 (Table 1) based on the initial uniaxial strain model and the critical model.In this case, the final predicted stress profiles are visually similar for both the stress-and strain-driven models because they remain essentially equal to the uniaxial stress profile.
For the other strategies, the results obtained with the stress-and strain-driven models show similar trends (Fig. 6).For both models, the magnitude of the local stress corrections depends on the magnitude and polarity (i.e., compression or extension) of the tectonic perturbations ε hr and σ hr .Depending on the tectonic perturbation, the resulting local stress corrections σ hl may promote, maintain, or inverse the initial layer-to-layer variations in the stress uniaxial profile due to the depth dependence of the Poisson's ratio.However, we note two significant differences between the models.
First, for the strain-driven model, the local stress corrections σ hl are constant within a layer, whereas they are concentrated near the formation interfaces for the stress-driven model.As a consequence, the final local stress gradient for the strain-driven model remains equal to the gradient of the initial state of stress.In the stress-driven model the final predicted stress gradient changes within a layer, as well as between layers.The gradient may be higher or lower than the stress gradient of the initial state of stress.
Second, for the stress-driven model, the local stress corrections σ hl fluctuate around the regional perturbation σ hr , i.e., the local stress correction averaged over several layers equals the regional perturbation.This average value is independent of the layering.For instance, for strategy 6 (Table 1) displaying the largest regional perturbation, the average local stress corrections calculated over the Fort Simpson and Horn River formations equals −33, and it is −34 MPa over the Horn River Formation and the Keg River Formation.
Conversely, for the strain model, the average magnitude varies as a function of the elastic properties of the studied section.For instance, very different average local stress corrections are found for strategy 2 (Table 1), i.e., −23 and −40 MPa for the Fort Simpson-Horn River formations and the Horn River-Keg River formations, respectively.
As a consequence, a strain-driven model may exhibit lower stress corrections than a stress-driven model in the compliant layers, and higher local stress corrections occur for the strain-driven model in the stiffer lithologies.In the case of an extensional perturbation, like in the study case, the predicted minimum horizontal stress σ hl is thus higher in compliant layers, e.g., the Fort Simpson Formation, and it is lower in stiff ones, e.g., the Keg River Formation, for the strain models.The horizontal stresses become very similar in the Horn River Formation, which has a mean stiffness and where the tectonic perturbation has been originally calibrated (Fig. 5).

Comparison of the stress predictions with a critical state of stress
On average, the predicted minimum principal stress σ h is close to the critical state of stress for all predictive strategies, except for strategy 6 (Fig. 5), implying generally consistent results.It is worth noting that the regulator has declared that hydraulic fracturing in this area has led to human-induced seismicity (BCOGC, 2012(BCOGC, , 2014)).Most of the larger magnitude events occurred within the Keg River Formation or below.The prediction of a critically stressed area is thus also reasonable.A few layers are predicted to be slightly infracritically stressed (i.e., σ 3c < σ h ) but more commonly, layers are supra-critically stressed (i.e., σ 3c > σ h ).This may indicate that the minimum principal stress σ h is locally underestimated.
A closer scrutiny of the various stress profiles in Fig. 5 shows that the strain-driven tectonic models (strategies 1-4, Table 1) predict that the largest supra-critical stresses occur in the Keg River Formation, whereas the stress-driven models (strategies 5-8, Table 1) also show supra-critical stresses above this formation.A critically stressed Keg River Formation is more likely given that the induced seismicity occurred within or below this formation.
Finally, predictions 2 and 6 (Table 1), based on an initial lithostatic stress and the in situ stress measurements, show the most supra-critical stresses independent of the chosen tectonic model.This can be explained by the fact that the in situ stress measurements are substantially lower than the critical stresses.Furthermore, using the lithostatic stress as the initial state induces higher tectonic effects than using any other initial stress state.As a consequence, stress predictions based on the lithostatic stress and in situ stress measurements creates strongly supra-critical predictions in this case.

Uncertainty in stress prediction
The various inputs and/or models give rise to inherent uncertainty in stress prediction.In the study case, the maximum differences in final stress obtained between all the stress predictions reach 20, 18, and 28 MPa in the Fort Simpson, Horn River, and Keg River formations, respectively (Fig. 5).For the different stress predictions that are based on a uniaxial stress, the maximum difference reaches 7, 6, and 8 MPa in the Fort Simpson, Horn River, and Keg River Formations.For an initial lithostatic stress, they reach 20, 15, and 28 MPa in the same formations.When using the critical stress as calibration stress, the maximum changes between the various stress predictions are 6, 9, and 16 MPa in the Fort Simpson, Horn River, and Keg River Formations.They are 15, 4, and 16 MPa in the same formations when using the in situ stress measurement as calibration stresses.Assuming a strain-driven model, they reach 8, 8, and 18, and they reach 14, 18, and 12 MPa assuming a stress driven model.Likewise, the range of predicted stress discontinuities between layers is 1-7 and 1-12.5 MPa between the Fort Simpson and Horn River formations and between the Horn River and Keg River formations, respectively.These uncertainties in stress prediction have fundamental implications on fracturing and containment capacity.

Discussion
In situ stress depends on several effects, including pore pressure, rock viscosity, rock anisotropy, preexisting fracturing, thermal effects, etc. (Cornet and Burlet, 1992;Addis et al., 1996;Khan et al., 2011;Meixner et al., 2014;Roche and Van der Baan, 2015).In this paper we focus on tectonic forces in combination with the vertical variation in rock properties while disregarding other potential effects.First, we discuss the choice of the initial and the calibration stress, then the stress-and the strain-driven model, and finally uncertainties in stress prediction.The choice of the initial regional tectonic regime and whether tectonic perturbations in one direction or in both horizontal directions are taken into account are also critical, but these are beyond the scope of the paper.

Initial state of stress
Stress predictions are based on an initial stress model that is modified to account for the tectonic effects.We used two models for the initial stress: the lithostatic and uniaxial strain models.In the literature most stress predictions are based on an initial uniaxial strain model (Voight and St. Pierre, 1974;Haxby and Turcotte, 1976;Savage, 1992;Blanton and Olson, 1999), and few are based on an initial lithostatic model (McGarr, 1988;Roche et al., 2014;Roche and Van der Baan, 2015).The choice of the initial stress model influences the final stress predictions for several reasons.
For a uniaxial strain model, initial layer-to-layer stress variations occur due to the Poisson's effect.These variations are added to those created by tectonic effects.In the study case, the tectonic effects dominate Poisson's effects.Still, the Poisson's effect may change the magnitude and direction of the final layer-to-layer stress variations if the tectonic effects are low, the variation in the Poisson's ratio is significant, and the pore pressure is low.
The magnitude of the regional tectonic perturbation depends on the chosen initial stress model.For instance, the initial stresses calculated with a uniaxial strain model are lower than those calculated with a lithostatic model.Hence, for an in situ stress measurement that is higher than the lithostatic stresses (i.e., compressive regime), the magnitude of the tectonic perturbations is greater for an initial uniaxial strain model than for a lithostatic model.As a consequence, stress predictions based on a uniaxial strain model tend to underestimate the layer-to-layer stress variations.For an extensional regime, the tectonic perturbation may be either greater or smaller for the lithostatic model, if the in situ stress measurement is closer to the lithostatic state of stress or to the uniaxial stress, respectively.
The choice of the initial model is potentially more critical when looking at its effect on the polarity of the regional tectonic perturbations.For instance, perturbations that are in compression or in extension may be obtained if the in situ stress measurements are comprised between the lithostatic and uniaxial stresses.In such a case, by assuming a tectonic perturbation in one direction, final stress predictions based on an initial uniaxial stress are likely to overestimate and underestimate the minimum stress in the stiff and compliant layers, respectively.Consequently, tensile and shear failures are more likely to appear inhibited in the stiff layer and promoted in the compliant layer.
This discussion raises the question: which model is more likely to dominate?The stress state predicted by the uniaxial strain model has rarely been observed in field data (Jaeger et al., 2009;McGarr, 1987McGarr, , 1988)).Also, it creates a bias toward highly extensional tectonic regimes and excludes thrust www.solid-earth.net/8/479/2017/Solid Earth, 8, 479-498, 2017 and strike-slip regimes (McGarr, 1987(McGarr, , 1988)).However, the lithostatic model also appears implausible.There is little experimental evidence to support this model because the timedependent failure mechanisms that could remove all the deviatoric stresses for indefinite periods of time do not seem to exist in the upper crust (Kirby and McCormick, 1984;Mc-Garr, 1988).Nevertheless, the lithostatic stress state appears to be the more realistic calibration stress (McGarr, 1988).

Calibration state of stress
The initial stresses influence the computed tectonic perturbations.Their magnitudes and polarities also depend on the calibration stresses.Two types of reference are used here: in situ stress measurements and critical stresses.In situ stress measurements are a better choice because they provide a direct measurement of the local stress.However, such data are often not available, have uncertainty, or are limited to specific depths.Thus, it is useful to be able to predict stress based only on models.This can be done using the critical state of stress.In the latter case, we assume that the stress is controlled by friction on preexisting cohesionless faults.Wrong estimation of this behavior will lead to a change in the critical stress.For instance, in clay-rich formations, a lower coefficient of friction may be expected.Likewise, it has been shown that the stress observed close to an active fault may imply a plastic behavior for the fault gouge that does not satisfy Coulomb failure criterion, but rather obeys a nonassociative plastic flow rule (Sulem, 2007).In this case, we may expect a lower critical stress than the one used in this paper.
If computation of the tectonic perturbation is based on stress measurements, then likely only a few points will determine their magnitudes.This is not the case if the critical stress model acts as a reference since then local differences are available for all layers.This points to a likely trade-off in the final predictions between paucity of information (lacking and uncertain measurements) versus potential for systematic errors due to an incorrectly assumed calibration state of critical stress.Likewise, uncertainty in the magnitude of the stress measurements is likely to be responsible for significant variation in the predicted stresses.
The magnitude of the maximum horizontal stress cannot be measured accurately or precisely (Schmitt et al., 2012).Hence, although we provide formulations involving the measurement of the maximum horizontal stress as a theoretical possibility, the use of in situ stress measurements as calibration stresses is only possible for a normal faulting regime since maximum horizontal stress measurements would be needed in case of a strike-slip or thrust faulting regime.
In this study, the difference between the calibration stresses, i.e., in situ measurement and critical stress, is relatively small compared to the difference between the initial stresses, i.e., lithostatic and uniaxial stresses.Hence, the change in tectonic perturbation, obtained for a similar initial stress but different calibration stresses, is only 25 % of the change obtained for similar calibration stress but different initial stresses.In our case, the stress predictions obtained using in situ measurements and critical stresses are thus similar.

Strain-driven models versus stress-driven models
For similar initial and calibration stresses, the stress predictions obtained with the strain-driven and the stress-driven models share a similar trend, consistent with in situ stress measurements, as a first approximation.For both models, stresses predicted in one specific layer depend on the tectonic perturbation and the elastic parameters.Nevertheless, significant differences occur in terms of stress magnitude and stress gradients (see Sect. 5.4).It is therefore fundamental to know which model dominates in nature, or if both models occur.
The strain-driven model has been widely used, notably as a standard method for oil and gas reservoir exploration (Thiercelin and Plumb, 1994;Blanton and Olson, 1999;Beaudoin et al., 2011;Song and Hareland, 2012).The stressdriven model may involve discontinuous strain across layer boundaries if the layers are not coupled together.Such a behavior may appear as a nonphysical consequence of the model, leading authors to disregard this model in favor of the strain-driven model (Blanton and Olson, 1999).This explains why the stress-driven model is scarcely used (Teufel and Clark, 1984;Bourne, 2003;Roche et al., 2013).However, if the layers are coupled together, both the regional stress and regional strain are continuous throughout the layering.Likewise, the possible occurrence of strain discontinuities does not appear decisive in disregarding the stress-driven model because strain decoupling likely exists in natural rocks, notably in the case of detachment faults (Cornet and Burlet, 1992;Meixner et al., 2014).Also, a recent study shows that fracturing depends on the contrast in elasticity between layers, as well as the layer thicknesses, rather than solely on the properties of the layer in which fracturing develops (Roche et al., 2014).This tends to support the stress-driven model rather than the strain-driven model.Lastly, bed-parallel faults occurring in tabular rocks have been described (Roche et al., 2012a, b).Such a structure may highlight strain discontinuities.This is also in accordance with the stress-driven model.Otherwise, for the strain-driven model, an additional mechanism must be involved to create these structures.

Final remarks: uncertainty in stress predictions
Accurate prediction of the in situ stresses in the Earth is hampered by epistemic and aleatoric uncertainty.Epistemic (or systematic) uncertainty is caused by lack of knowledge and data; aleatoric (or statistical) uncertainty is caused by the inherent randomness of a phenomenon such as the rolling of dice.The variation in final predicted stresses stems from a range of different sources for uncertainty, including observational and interpolation-extrapolation uncertainty, parameter uncertainty, model uncertainty, and numerical-algorithmic uncertainty (Kennedy and O'Hagan, 2001).
Observational, interpolation, and extrapolation uncertainty arises, for instance, since the stress measurements used in step 2 are prone to observational error, thus influencing the final stress predictions.In addition, the stress measurements are often limited to specific layers and sites, thereby requiring interpolation and/or extrapolation to derive values appropriate for the areas and depth zones under consideration.This introduces uncertainty and spread in the final predicted stress values.
Next, parameter uncertainties are also important, for instance, in assumed elastic parameters (Young's modulus, Poisson's ratio) and friction coefficients.Some of these parameters are derived from well logs and are thus prone to observational error; yet even in the absence of such observational uncertainty, we lack knowledge of how to exactly convert measured parameters (e.g., dynamic moduli) to required modeling parameters (e.g., static moduli), the control of viscoelasticity, or the behavior for the faults, thus creating parameter uncertainty.Furthermore, each predictive strategy depends on different parameters, with their own uncertainties and possible biases, thus introducing further diversity in the final predicted stresses.
Moreover, we have model uncertainty because full knowledge of the actual driving forces and most appropriate geologic boundary conditions is lacking.For instance, both the stress-driven and strain-driven models are geologically plausible.Likewise, lithostatic, uniaxial, and critical states of stress are reasonable assumptions in various circumstances, yet the most appropriate one is rarely known.Also, assumptions implicit in the derivation of the provided equations create model uncertainty.Again, lack of knowledge (epistemic uncertainty) causes diversity in final predictions.Finally, there is numerical (or algorithmic) uncertainty, caused by numerical errors and approximations when solving for the stress state, e.g., using the discrete-element method in step 4 or in the implementation of any of the provided equations.We assume that these only play a very minor role in our case; for instance, convergence of solutions was closely monitored.Nonetheless, numerical uncertainty remains a possible source of uncertainty in all computer simulations.

Conclusions
Different strategies to predict the vertical variations in the in situ stresses lead to different answers.Such stress predictions take into account the weight of the overburden rock, the pore pressure, the variation of the rock properties, and the tectonic effects.In addition, they assume either stress-or strain-driven models, they assume an initial uniaxial or lithostatic model, and they use a critical model or in situ stress measurements.The different prediction strategies generally lead to similar trends in predicted stresses, yet differences appear both within a layer and between layers due to fundamentally different underlying mechanisms, assumptions, and governing parameters.The spread and diversity in final stress predictions is caused by mostly epistemic uncertainty, expressed as lack of knowledge, data, and/or observational errors in some measurements or variables.Nonetheless, the combined analysis of all eight stress predictions helps reveal the uncertainty, or conversely similarity, in all stress predictions.
V. Roche and M. van der Baan: Stress and strain-driven tectonic forces Appendix A: Lithostatic, uniaxial, and critical state of stresses In the lithostatic model, the principal stresses σ 1 , σ 2 , and σ 3 are equal to the lithostatic stress σ l calculated from the overburden stress as determined by the density ρ.Thus, where σ v is the overburden vertical stress, g is the acceleration due to gravity, and z is the depth.This model assumes that rocks cannot support differential stresses, such that horizontal stresses equalize to the overburden (vertical) stress over a long period of time due to inelastic deformations.This model may be used for low viscosity, plastic rocks, or for a relatively long period of time without modification in external boundary conditions.For instance, a stress field close to lithostatic has been described in shales (Warpinski et al., 1985).The uniaxial strain model assumes that no regional horizontal strains exist but horizontal stresses are imposed by Poisson's effect (that is, a horizontal force due to vertical loading; Savage et al., 1992;Jaeger et al., 2009).The overburden stress σ v creates a uniaxial horizontal stress σ u due to the lateral extension of the medium that is impeded by non-deformable walls (i.e., no lateral strains).In an elastic isotropic rock, we obtain (Engelder and Fischer, 1994;Addis et al., 1996;Blanton and Olson, 1999) where P p is the in situ pore pressure, υ is the Poisson's ratio, α is the Biot pore-pressure coefficient, and the apostrophe ' indicates an effective stress.Equation (A2) holds if no poroelastic coupling occurs.With a Poisson's ratio of 0.3, the ratio between the uniaxial and vertical stresses equals 0.43.This model explicitly excludes compressive regimes (Jaeger et al., 2009;McGarr, 1988).At best it shows qualitative agreement with stress measurements but rarely an accurate prediction, while in other cases, stress measurements are in contradiction to the model (Warpinski et al., 1985;McLellan, 1987;Whitehead et al., 1987;Ahmed et al., 1991;Plumb et al., 1991;Thiercelin and Plumb, 1994;Addis et al., 1996).However, its main advantage is that horizontal stresses are easy to compute since only knowledge of the density, pore pressure, and Poisson's ratios are required.
In the critical stress model, the horizontal stresses are assumed equal but the ratio between horizontal and vertical stress is set such that preexisting, optimally oriented faults and fractures are at the point of shear failure.The ratio of effective maximum to minimum critical principal stresses σ 1c / σ 3c then only depends on the frictional strength of the preexisting faults as their cohesion is set to zero (Zoback, 2010), that is, where µ is the coefficient of friction.For µ = 0.6, the ratio between the critical stresses equals 3.1.An equivalent ratio is obtained with a Poisson's ratio equal to 0.25 for the uniaxial strain model.The stresses calculated using the uniaxial strain model with Poisson's ratios lower than 0.25 are thus overcritical.

Appendix B: Coupled stress-driven model
In order to analytically describe such a behavior, we assume that only one regional tectonic perturbation σ hr exists applied in the direction of the minimum principal stress.In such a case, the regional strain continuity between the layers can be expressed with the following equation: where ε hrc and ε hrs are the regional strain perturbations in the compliant and stiff layers, respectively (Fig. 1f).The subscripts c and s refer to the compliant and stiff layers, respectively.The strain in each layer can be split into two components, namely a strain correction for the non-coupled case and a second correction due to coupling, producing ε hrs = ε hsnc + ε hst , and ε hrc = ε hcnc + ε hct , (B2) where ε hsnc and ε hcnc are the local strain corrections obtained for non-coupled layers and ε hst and ε hct are the local strain corrections, induced by the coupling.The latter lead to strain transfer, denoted by the subscript t.Combining Eqs.(B1) and (B2), we obtain The various strain corrections can be computed from Eq. ( 10) by assuming that a single strain perturbation in one horizontal direction occurs, with zero stress change and zero strain in the second horizontal direction, producing The equilibrium condition sets that there is zero net force across a plane that is normal to the layers (Holzhausen and Johnson, 1979;McGarr, 1988): where h s and h c are the stiff and compliant layer thicknesses, respectively.Using Eqs. ( B6) and (B7), the local stress transfers in the stiff and compliant layers become

Figure 1 .
Figure 1.Conceptual representation of the stress-and strain-driven models.(a)-(f) stiff layer on top with Young's modulus, E s , and compliant layer underneath with Young's modulus, E c , represented by thick and thin springs, respectively.Dashed wall is non-deformable (fixed).A single regional perturbation occurs in the direction of the compression (a, c, e) or extension (b, d, f).(a)-(b) Strain-driven model with imposed regional strain perturbation that is positive (shortening), ε Hr , or negative (lengthening) ε hr .(c)-(f) Imposed regional stress perturbation that is positive (i.e., compressional), σ Hr , or negative (i.e., extensional), σ hr .(c)-(d) Non-coupled stress-driven model; (e)-(f) fully coupled stress-driven model.Each subfigure shows the boundary conditions (next to the bold horizontal arrows), a conceptual spring diagram, the resulting local stress and strain corrections (underneath the springs), and the resulting final local stresses in each layer.Subscripts "h" and "H" are the minimum and maximum principal stresses in the direction of the extension and compression, respectively.Subscripts c and s are the compliant and stiff layers.σ hls , σ Hls , σ hlc , and σ Hlc are the local stress corrections; ε hs , ε Hs , ε hc , and ε Hc are the local strain corrections resulting from the regional tectonic perturbations, σ hr , σ Hr , ε hr , and ε Hr .Absolute values of local stress and strain corrections are considered for the comparisons.Local stress equalization happens in the non-coupled stress-driven models (c, d).Both strain-driven and coupled stress-driven models can lead to similar local stress profiles, but with very different boundary conditions (a, b, e, f).

Figure 2 .
Figure 2. Examples of stress predictions, σ 3 , and local stress corrections, σ hl , in the strain-(a) and stress-driven models (b, c).(ac) Bilayer media of a 200 m thick compliant layer (top) with Young's modulus E c and a 200 m thick stiff layer (bottom) with Young's modulus E s .Various combinations of Young's moduli are used as indicated in the figure.The Poisson's ratio is constant at 0.25.(a) A regional strain perturbation of 0.002 m m −1 extension is imposed.Stresses are calculated using the analytic solution described in Sect.2.2.(b, c) A regional stress perturbation of −10 MPa is imposed.Stresses are calculated using the numerical (b) and analytic solution (c) as described in Appendix B. The initial horizontal stress is the lithostatic stress (σ i = σ l ) calculated using an average 2400 m depth and an average density of 2600 kg m −3 .Local stress corrections ( σ hl ) and the minimum horizontal stress (σ 3 ) are the dotted and continuous colored lines for indicated combinations of Young's moduli.

Figure 3 .
Figure 3. Flow chart detailing the different prediction strategies.See Sect. 3 for details.Rock properties used in the different models are υ Poisson's ratio, µ friction, ρ density, and E Young's modulus.

Figure 4 .
Figure 4. Petrophysical and mechanical properties of the study case.(a) Lithological section: Fort Simpson, Muskwa, Otter Park, and Evie members are clay-rich rocks represented in various shades of gray.The Keg River Formation is carbonate.The Muskwa, Otter Park, and Evie members together form the Horn River Formation.The white stars represent the locations of stress measurements.(b, c) Depth dependence of the density and the P-and S-wave velocities are derived from the well data.(d, e) Variation of the Young's modulus and Poisson's ratio with depth are calculated from the density and velocities, respectively.In (d) the dynamic and static Young's moduli are indicated.

Figure 5 .
Figure 5. Calculated stresses for the various models and predictive strategies.(a) Uniaxial (σ u ), lithostatic (σ l ), and critical (σ 3c ) models.(b-g)Stress predictions for the various predictive strategies, including the local stress corrections σ hl and the predicted minimum principal stress σ hl .Circled numbers are labels of the predictive strategies as presented in Table1.The relevant initial stress (i.e., uniaxial σ u or lithostatic σ l ) and calibration stresses (i.e., in situ stress measurements and the critical stress σ 3c ) are also represented for each predictive strategy.The diamonds represent the average values of the minimum principal stress (set equal to the instantaneous shut-in pressure values).(e) and (i) Show simplified stratigraphic column of the studied area with the same legend as for Fig.4.Cases 3 and 7 in Table2are not represented because the tectonic effects are negligible, and therefore the predicted minimum principal stress σ hl is equal to the calibration stress, σ u .The Muskwa (M.), Otter Park (O.P.) and Evie (E.) members together form the Horn River Formation, surrounded by the Keg River (K.R.), Fort Simpson (F.S.), and Red Knife (R.K) formations.

Figure 6 .
Figure 6.Comparison of the stress-and strain-driven models.Stress predictions for the stress-driven and the strain-driven models for (a) predictive strategies 1 and 5, based on the uniaxial stress and in situ stress measurements; (b) predictive strategies 2 and 6, based on the lithostatic stress and in situ stress measurements; and (c) predictive strategies 4 and 8, based on the lithostatic stress and critical model.Also shown are the initial uniaxial stresses σ u (a), the initial lithostatic stresses σ l (a-c), and the critical stresses σ 3c (c).(d) Simplified stratigraphic column of the studied area with the same legend as for Figs. 4 and 5.
the compliant layer, where σ hct and σ hst are the local stress transfers in the compliant and stiff layers, respectively.E s and E c are the Young's moduli in the compliant and stiff layers, respectively, and υ s and υ c are the Poisson's ratios in the compliant and stiff layers, respectively.Using Eqs.
h s (A − B) σ hr h s B + h c A = σ hct and σ hst = − stress corrections then correspond to the sums of the local stress transfers due to coupling and the regional stress perturbations.That is, σ hls = σ hr + σ hst and σ hlc = σ hr + σ hct .(B10)Finally,we can calculate the local stress in each layer with the following equations:σ hlc = σ hi − P p + σ hr + σ hct ,(B11) σ hls = σ hi − P p + σ hr + σ hst .(B12) Equations (B11)-(B12) hold for periodic media composed of two layers with Young's moduli E and E c and thicknesses h s and h c , where all layers are coupled.www.solid-earth.net/8/479/2017/Solid Earth, 8, 479-498, 2017

Table 1 .
The eight predictive strategies.Predictive strategy 1 is commonly used in literature.
Hr and σ hl =