3-D geomechanical–numerical model of the contemporary crustal stress state in the Alberta Basin (Canada)

In the context of examining the potential usage of safe and sustainable geothermal energy in the Alberta Basin, whether in deep sediments or crystalline rock, the understanding of the in situ stress state is crucial. It is a key challenge to estimate the 3-D stress state at an arbitrarily chosen point in the crust, based on sparsely distributed in situ stress data. To address this challenge, we present a large-scale 3-D geomechanical–numerical model (700 km×1200 km×80 km) from a large portion of the Alberta Basin, to provide a 3-D continuous quantification of the contemporary stress orientations and stress magnitudes. To calibrate the model, we use a large database of in situ stress orientation (321 SHmax) as well as stress magnitude data (981 SV, 1720 Shmin and 2 (+11) SHmax) from the Alberta Basin. To find the best-fit model, we vary the material properties and primarily the displacement boundary conditions of the model. This study focusses in detail on the statistical calibration procedure, because of the large amount of available data, the diversity of data types, and the importance of the order of data tests. The best-fit model provides the total 3-D stress tensor for nearly the whole Alberta Basin, and allows estimation of stress orientation and stress magnitudes in advance of any well. First-order implications for the well design and configuration of enhanced geothermal systems are revealed. Systematic deviations of the modelled stress from the in situ data are found for stress orientations in the Peace River and the Bow Island Arch as well as for leak-off test magnitudes. 1 Motivation The estimation of the in situ stress state in the upper crust, in addition to the understanding of earthquake cycles and plate tectonics, is crucial for exploration and production of energy resources. These include geothermal energy, hydrocarbons, CO2 sequestration (carbon capture storage – CCS) and geotechnical subsurface constructions such as mines, tunnels, interim storage sites for natural gas and nuclear waste deposits. Reliable estimates of orientation and magnitude of the crustal stress field are desired before drilling. This is important in terms of well stability (e.g. Bell and McLellan, 1995; Peska and Zoback, 1995), but is also related to the well configuration of several corresponding wells (e.g. Bell and McLellan, 1995) in the case of reservoir stimulation by hydraulic fracturing. This is important in geothermal reservoirs (enhanced geothermal systems – EGS) (Legarth et al., 2005; Wessling et al., 2009) and issues of inadequate understanding of the spatial stress pattern (e.g. Brown, 2009; Duchane and Brown, 2002). This is also true for hydrocarbon reservoirs or the evaluation of nuclear waste repositories (e.g. Fuchs and Müller, 2001; Gunzburger and Magnenet, 2014; Heidbach et al., 2013). The stress tensor and its components are not to be measured directly, but there are several stress indicators which allow estimation of several components of the stress tensor (e.g. Ljunggren et al., 2003; Schmitt et al., 2012; Zang and Stephansson, 2010). The following components of the stress tensor are potentially available: the azimuth of the maximum (or minimum) horizontal stress (SHmax), the vertical stress magnitude (SV), as well as the magnitudes of the maximum and minimum horizontal stress (Shmin and SHmax). However, Published by Copernicus Publications on behalf of the European Geosciences Union. 1124 K. Reiter and O. Heidbach: Calibration of a 3-D crustal stress model – Alberta Basin reliable estimation of the SHmax magnitude remains difficult, as only shallow in situ stress estimations are available, or numerous assumptions have to be made that impose high uncertainties. Furthermore, stress information is sparse, and extraor inter-polation of a few data records to the area or depth of interest is necessary. However, stress estimation via interpolation techniques becomes particularly questionable in the case of structural inhomogeneities like faults, detachments (Bell and McLellan, 1995; Röckel and Lempp, 2003; Roth and Fleckenstein, 2001; Yassir and Bell, 1994), or varying material properties (Roche et al., 2013; Warpinski, 1989). Furthermore, drilling down to a geothermal reservoir requires reaching greater depths, as available measurements are delivered in the context of hydrocarbon production. For example, in Alberta, deeper parts of the basin or the upper basement are the target depths for EGS (e.g. Hofmann et al., 2014; Majorowicz and Grasby, 2010a, b; Weides et al., 2013, 2014). Therefore, estimation of the stress state, especially at greater depths, is a challenge prior to drilling. An alternative approach to estimating the 3-D stress state is geomechanical–numerical modelling. This method has the advantage of incorporating structural and material inhomogeneities that impose local to regional changes on the stress field. There are several studies on tectonic plate scale stress orientation patterns in 2-D (e.g. Coblentz and Richardson, 1996; Dyksterhuis et al., 2005; Humphreys and Coblentz, 2007; Jarosinski et al., 2006), large-scale (regional) models in 3-D (Buchmann and Connolly, 2007; Hergert and Heidbach, 2011; Parsons, 2006), as well as local (reservoir-scale) 3-D models (e.g. Fischer and Henk, 2013; Heidbach et al., 2013; Henk, 2005; Orlic and Wassing, 2012; Van Wees et al., 2003). Modelling of the contemporary stress field mainly depends on the structural model, the material properties, the initial stress state and the applied displacement boundary conditions. However, the reliability of such models depends strongly on the model calibration towards in situ stress data. Usually, there are little in situ stress data available for model calibration in published studies (e.g. Buchmann and Connolly, 2007; Fischer and Henk, 2013; Heidbach et al., 2013; Hergert and Heidbach, 2011), which rules out any statistical validation. The Alberta Basin is a study area with well-understood structures and material properties, and a large collection of in situ stress data. We use this information to build a 3-D geomechanical–numerical model of the Alberta Basin and surroundings with an extent of 1200 km×700 km down to a depth of 80 km. The goal is to get the full tensor of the contemporary undisturbed stress state, called stress model in the following. These are 981 SV magnitude data, 321 SHmax azimuth data, 1720 Shmin magnitudes, and 2 measured (overcoring) and 11 calculated SHmax magnitudes within the model region. There is no other basin with a comparable range of available in situ stress data (Bell and Grasby, 2012). The availability of very good stress data allows for the calibration of the stress model vs. a never-reached diversity, and a number of in situ stress indicators. The model calibration will be done in three consecutive steps: (1) density of basin infill, using SV magnitude data, (2) orientation of displacement boundary conditions using SHmax azimuth data, and (3) magnitudes of displacement boundary conditions (strain) using Shmin and SHmax magnitude data. As linear elastic rheology is used for the model, the linear dependency between the two applied strain magnitudes (push and pull) along the outer edges of the model is calculated. This allows, via planar regressions the calculation of the optimal strain magnitudes, providing the best-fit model. The application of the model would be for exploitation of hydrocarbons and more for exploration and design of a geothermal plant in the Alberta Basin. Additionally it may be used in crystalline rocks, mainly in case of necessary hydraulic stimulation. Mistaken investments e.g. parts of the Fenton Hill project (Brown, 2009; Duchane and Brown, 2002) could potentially be avoided with a better previous understanding of the 3-D in situ stress state.


Abstract.
In the context of examining the potential usage of safe and sustainable geothermal energy in the Alberta Basin, whether in deep sediments or crystalline rock, the understanding of the in situ stress state is crucial. It is a key challenge to estimate the 3-D stress state at an arbitrarily chosen point in the crust, based on sparsely distributed in situ stress data.
To address this challenge, we present a large-scale 3-D geomechanical-numerical model (700 km×1200 km×80 km) from a large portion of the Alberta Basin, to provide a 3-D continuous quantification of the contemporary stress orientations and stress magnitudes. To calibrate the model, we use a large database of in situ stress orientation (321 S Hmax ) as well as stress magnitude data (981 S V , 1720 S hmin and 2 (+11) S Hmax ) from the Alberta Basin. To find the best-fit model, we vary the material properties and primarily the displacement boundary conditions of the model. This study focusses in detail on the statistical calibration procedure, because of the large amount of available data, the diversity of data types, and the importance of the order of data tests.
The best-fit model provides the total 3-D stress tensor for nearly the whole Alberta Basin, and allows estimation of stress orientation and stress magnitudes in advance of any well. First-order implications for the well design and configuration of enhanced geothermal systems are revealed. Systematic deviations of the modelled stress from the in situ data are found for stress orientations in the Peace River and the Bow Island Arch as well as for leak-off test magnitudes.

Motivation
The estimation of the in situ stress state in the upper crust, in addition to the understanding of earthquake cycles and plate tectonics, is crucial for exploration and production of energy resources. These include geothermal energy, hydrocarbons, CO 2 sequestration (carbon capture storage -CCS) and geotechnical subsurface constructions such as mines, tunnels, interim storage sites for natural gas and nuclear waste deposits. Reliable estimates of orientation and magnitude of the crustal stress field are desired before drilling. This is important in terms of well stability (e.g. Bell and McLellan, 1995;Peska and Zoback, 1995), but is also related to the well configuration of several corresponding wells (e.g. Bell and McLellan, 1995) in the case of reservoir stimulation by hydraulic fracturing. This is important in geothermal reservoirs (enhanced geothermal systems -EGS) (Legarth et al., 2005;Wessling et al., 2009) and issues of inadequate understanding of the spatial stress pattern (e.g. Brown, 2009;Duchane and Brown, 2002). This is also true for hydrocarbon reservoirs or the evaluation of nuclear waste repositories (e.g. Fuchs and Müller, 2001;Gunzburger and Magnenet, 2014;Heidbach et al., 2013).
The stress tensor and its components are not to be measured directly, but there are several stress indicators which allow estimation of several components of the stress tensor (e.g. Ljunggren et al., 2003;Schmitt et al., 2012;Zang and Stephansson, 2010). The following components of the stress tensor are potentially available: the azimuth of the maximum (or minimum) horizontal stress (S Hmax ), the vertical stress magnitude (S V ), as well as the magnitudes of the maximum and minimum horizontal stress (S hmin and S Hmax ). However, reliable estimation of the S Hmax magnitude remains difficult, as only shallow in situ stress estimations are available, or numerous assumptions have to be made that impose high uncertainties. Furthermore, stress information is sparse, and extraor inter-polation of a few data records to the area or depth of interest is necessary.
However, stress estimation via interpolation techniques becomes particularly questionable in the case of structural inhomogeneities like faults, detachments (Bell and McLellan, 1995;Röckel and Lempp, 2003;Roth and Fleckenstein, 2001;Yassir and Bell, 1994), or varying material properties (Roche et al., 2013;Warpinski, 1989). Furthermore, drilling down to a geothermal reservoir requires reaching greater depths, as available measurements are delivered in the context of hydrocarbon production. For example, in Alberta, deeper parts of the basin or the upper basement are the target depths for EGS (e.g. Hofmann et al., 2014;Majorowicz and Grasby, 2010a, b;Weides et al., 2013. Therefore, estimation of the stress state, especially at greater depths, is a challenge prior to drilling. An alternative approach to estimating the 3-D stress state is geomechanical-numerical modelling. This method has the advantage of incorporating structural and material inhomogeneities that impose local to regional changes on the stress field. There are several studies on tectonic plate scale stress orientation patterns in 2-D (e.g. Coblentz and Richardson, 1996;Dyksterhuis et al., 2005;Humphreys and Coblentz, 2007;Jarosinski et al., 2006), large-scale (regional) models in 3-D (Buchmann and Connolly, 2007;Hergert and Heidbach, 2011;Parsons, 2006), as well as local (reservoir-scale) 3-D models (e.g. Fischer and Henk, 2013;Heidbach et al., 2013;Henk, 2005;Orlic and Wassing, 2012;Van Wees et al., 2003). Modelling of the contemporary stress field mainly depends on the structural model, the material properties, the initial stress state and the applied displacement boundary conditions. However, the reliability of such models depends strongly on the model calibration towards in situ stress data. Usually, there are little in situ stress data available for model calibration in published studies (e.g. Buchmann and Connolly, 2007;Fischer and Henk, 2013;Heidbach et al., 2013;Hergert and Heidbach, 2011), which rules out any statistical validation.
The Alberta Basin is a study area with well-understood structures and material properties, and a large collection of in situ stress data. We use this information to build a 3-D geomechanical-numerical model of the Alberta Basin and surroundings with an extent of 1200 km×700 km down to a depth of 80 km. The goal is to get the full tensor of the contemporary undisturbed stress state, called stress model in the following. These are 981 S V magnitude data, 321 S Hmax azimuth data, 1720 S hmin magnitudes, and 2 measured (overcoring) and 11 calculated S Hmax magnitudes within the model region. There is no other basin with a comparable range of available in situ stress data (Bell and Grasby, 2012). The availability of very good stress data allows for the cali-bration of the stress model vs. a never-reached diversity, and a number of in situ stress indicators.
The model calibration will be done in three consecutive steps: (1) density of basin infill, using S V magnitude data, (2) orientation of displacement boundary conditions using S Hmax azimuth data, and (3) magnitudes of displacement boundary conditions (strain) using S hmin and S Hmax magnitude data. As linear elastic rheology is used for the model, the linear dependency between the two applied strain magnitudes (push and pull) along the outer edges of the model is calculated. This allows, via planar regressions the calculation of the optimal strain magnitudes, providing the best-fit model. The application of the model would be for exploitation of hydrocarbons and more for exploration and design of a geothermal plant in the Alberta Basin. Additionally it may be used in crystalline rocks, mainly in case of necessary hydraulic stimulation. Mistaken investments e.g. parts of the Fenton Hill project (Brown, 2009;Duchane and Brown, 2002) could potentially be avoided with a better previous understanding of the 3-D in situ stress state.

Model assumptions
The compilation of stress data in North America by Adams (1987, Adams and Bell (1991), Bell et al. (1994), Fordjor et al. (1983), Gough et al. (1983), Sbar and Sykes (1973), Zoback and Zoback (1980, 1981 and recently by Reiter et al. (2014) resolved that the pattern of S Hmax orientations is largely uniform over thousands of kilometres. An assumption was that the same forces driving plate tectonics are the major control on the stress field, which is confirmed in first order (e.g. Richardson, 1992;Zoback et al., 1989;Zoback, 1992).
The stress pattern is driven and altered by several stress sources; they are discriminated depending on the scales into first-order (> 500 km), second-order (100-500 km) and third-order stress sources (< 100 km) (Heidbach et al., 2007Müller et al., 1997;Tingay et al., 2005;Zoback, 1992;Zoback and Mooney, 2003). First-order stress sources as the main driving forces are summarized as plate boundary forces, which are ridge push, slab pull, and trench suction, gravity and basal drag by mantle convection. Second-order stress sources are lithospheric flexure, localized lateral density, stiffness and strength contrasts, topography, large fault zones, and lateral contrasts of heat production. Third-order stress sources are local density, stiffness or strength contrasts, basin geometry, basal detachment, incised valleys, and anthropogenic stress changes.
Under the parameters to reproduce the crustal stress field of the Alberta Basin, the model has to be large enough to portray the first-and second-order stress sources. Whereas the first-order stress sources, which control plate tectonic motion, are represented by the displacement boundary conditions, the second-and third-order stress sources are represented by the model geometry. This is possible when structures are known and convertible to the model. As inhomogeneous topography and mass distribution within the lithosphere have a major impact on the stress orientation (Camelbeeck et al., 2013;Ghosh et al., 2009;Humphreys and Coblentz, 2007;Naliboff et al., 2012), it is crucial to incorporate the major structural units into the crust and upper mantle in the model geometry.
Linear elastic material properties are an accurate approximation as long as the strain ( ) is small enough, and no failure occurs. This might be assumed for the Alberta Basin, which is seismically relatively quiescent. Documented earthquakes are usually restricted to the Rocky Mountains, foothills, and suspected man-made clustered events (e.g. Baranova et al., 1999;Schultz et al., 2014).
Our focus is the Alberta Basin and the uppermost basement below the basin due to our key interest in the investigation of the potential usability of deep geothermal energy (EGS). Furthermore, the calibration data are from the sediments up to a depth of about 5 km when no deeper stress indicators are available. Exceptions to this are three S Hmax azimuths, derived from focal mechanism solutions. Viscous rock deformation, acceleration, changing pore pressures as well as other thermal effects influence the stress state, but as we strive to model the contemporary undisturbed total stress tensor (stress model), these processes are assumed to have only marginal influence, and are disregarded in the model. We assume that the following model assumptions are sufficient to estimate the contemporary 3-D stress state within the Alberta Basin and upper basement: -large-scale geometry of the model down to the upper mantle is crucial (main structural units only).
-linear elastic material properties are used.
-gravity as the body force.
-the lateral displacement boundary conditions of the model are a parameterization of ongoing plate tectonic motion, effects of lateral density contrasts (gravitational potential energy) of outside of the model, and remnant stresses from terminated tectonic processes.
Due to the complex structures and inhomogeneous materials properties, an analytical solution cannot be estimated. Thus, we use for the discrete solution the finite element method (FEM), as it allows us to use unstructured meshes, for a good representation of the 3-D model structure. With these assumptions, the model is described with the partial differential equations of the equilibrium of forces: Model-independent data: Figure 1. Sketch of the general workflow. The geomechanical model is prepared based on the model geometry, the material properties, the variable displacement boundary conditions and the initial stress state. The numerically modelled total stress tensor is calibrated on model-independent in situ stress data until the model fits the calibration data.
where ∂σ ij is the change in total stress, ∂x the change in length, and ρx represents the weight of the rock section (ρ = density). After model design (structural model) and definition of material properties (Poisson's ratio, Young's modulus and density), the partial differential equation (Eq. 1) can be calculated within given displacement boundary conditions. The latter will be varied to find the best-fit model. This is the stress model together with material properties and displacement boundary conditions, which deliver the best-fit for all in situ stress data.

General workflow of model calibration
Generally, a model has to be calibrated before application or interpretation. The general concept, independent of the technical context, is to test the model's outcome vs. in situ data. Such data are called model-independent data, in contrast to model-dependent data, which are used to generate the model.
In this study, the lithological and tectonic structures, the rheology, the body force, the initial stress state, and the displacement boundary conditions, are the model-dependent data (Fig. 1). Based on these data, the structural model is defined, which is discretized to a (unstructured) mesh and assembled together with the material properties, body forces, the boundary conditions, and the initial stress state. Available in situ stress data are the model-independent data. These are the S V magnitudes, the S Hmax azimuth data, and S hmin and S Hmax magnitudes. The modelled stress tensor is tested against the in situ data. When one data set is tested successfully, the next data set is used in the next calibration step. Otherwise, the material properties or boundary conditions are optimized as long as the test is successful.
First, the stress model is tested vs. in situ S V magnitudes, to conclude estimation of density (material properties and body force). In the second step, the S Hmax orientation is tested to determine the orientation of applied displacement boundary conditions. In the final step, S hmin and S Hmax magnitudes are used to calibrate the applied magnitudes of the displacement boundary conditions. When all model-independent data sets are tested successfully, the bestfit model is found and is a subject of further use (interpretation and application).
The model-dependent data, construction and compilation process of the geomechanical model are described in Sect. 3, whereas the model-independent data are introduced in Sect. 4. The calibration procedure is presented in detail in Sect. 5. Finally, the discussion can be found in Sect. 6.

Tectonic and sedimentary history of the Alberta Basin
The Alberta Basin (Fig. 2) occupies a large portion of the much larger Western Canada Sedimentary Basin (WCSB). Starting from the northeast, clockwise it is bounded by the Canadian Shield, the Bow Island Arch, the Rocky Mountains and the Tathlina High in the north. The crystalline basement of the WCSB and, implicitly, of the superposed Alberta Basin, is the North American craton exposed by erosion to the northeast as the Canadian Shield Flowers et al., 2012;Hoffman, 1989;Ross et al., 1994Ross et al., , 2000. The main structural units of the Alberta basement are the Buffalo Head Terrane (Aulbach et al., 2004), the Taltson Magmatic Zone (e.g. Chacko et al., 2000), the Hearne Province  and the Trans-Hudson Orogen (e.g. Corrigan et al., 2005;Németh et al., 2005) and other smaller units, which welded together between 1.8 and 2.0 Ga. There are two important lineaments, the Snowbird Tectonic Zone (STZ - Ross et al., 2000, and references therein) and the Great Slave Lake Shear Zone (GLS - Sami and James, 1993), and their continuation (Hay River fault zone). Sediments were deposited in the basin, interrupted by a few discontinuities during the whole Phanerozoic (Mossop and Shetsen, 1994a, Chapter 6-26). Mainly shelf sediments deposited onto the craton as recently as the Upper Jurassic. At that time, sedimentation character changed, and the basin developed to a rapidly subsiding fore-deep trough (Poulton et al., 1994). Mature sediments were previously derived from the northeast, and changed to less mature sedi- ments derived from the west. The change to terrestrial deposits in Early Cretaceous (Smith, 1994) coincides with first the Ominicean Orogeny, and later the Lamariden Orogeny (Porter et al., 1982;Price, 1981;Wright et al., 1994). Jurassic to Palaeocene strata mainly deposited in the western part of the Alberta Basin and have been incorporated in the Rocky Mountains fold-and-thrust belt (foothills and front ranges - Fig. 3). This is bound farther west (main ranges) in British Columbia by the Rocky Mountain Trench. The final shape of the Alberta foreland basin developed by downward flexing of the Canadian Shield due to lithospheric loading and isostatic flexure in a retro-arc setting (Leckie and Smith, 1992;English and Johnston, 2004), together with the sediments derived from the developing Canadian Cordillera (Gabrielse and Yorath, 1989). The Alberta Basin consists of a nearly undeformed sedimentary wedge (Fig. 3) that increases in thickness from zero at the Canadian Shield to approximately 5500 m near the fold-and-thrust belt. The overall wedge shape in the Alberta Basin, perpendicular to the Rocky Mountains, is quite homogeneous from northwest to southeast.
Only the Peace River Arch close to the Rocky Mountains is striking within the homogeneous wedge, which is indicated by several geophysical investigations. There are several explanations: an elevated Precambrian basement (Bell and Babcock, 1986;Bell, 1996b;Bell and Grasby, 2012;Halchuk and Mereu, 1990), the occurrence of mafic sills which intruded in the upper crust of the Peace River Arch (Eaton et al., 1999), and/or lateral heterogeneities (transfer zone or local rheological properties in Bell and McCallum, 1990), a softer inclusion (Dusseault and Yassir, 1994), or crustal thinning caused by extension (Bouzidi et al., 2002).

Model geometry
The structural-geological model is of major importance, and is prepared with the GOCAD ® 3-D geomodelling system. Faults and lithological boundaries are defined as discrete triangulated surfaces. These are built based on points (stratigraphic borehole data and seismic data), curves (seismic and interpreted cross sections, lineaments from the geological map) or point clouds that describe surfaces (DEM). During surface generation, data are honoured as soft or hard constrained, depending on their quality (e.g. Ross et al., 2004). The roughnesses of the surfaces are polished with the discrete smoothing interpolation (DSI) algorithm (Mallet, 1992(Mallet, , 2002. The model box of Alberta, indicated in Fig. 4, is oriented parallel or perpendicular to the observed basin structure (Fig. 2), the orientation of S Hmax (e.g. Bell et al., 1994;Reiter et al., 2014, see rose diagram in Fig. 4), the wedge shape of the Alberta Basin (Fig. 3), the thermally defined Cordillera-Craton boundary (Hyndman et al., 2009) and the overall plate motion of the North American Craton, measured by GPS (Henton et al., 2006;Mazzotti et al., 2011).
The model has a southwest-to-northeast striking extent of 700 km, 1200 km in the northwest-to-southeast direction (Fig. 4), and 80 km in depth. For the definition of the model geometry, it was necessary to choose the geometrically relevant structures, strength contrasts or density variations. These can potentially affect the stress field, while considering limitations of the possible number of finite elements. The main structural units are the mantle, the crustal basement, the sedimentary basin, the foothills, the Rocky Mountains and the Elk Point evaporates within the basin, due to their Implemented are the main structural features (red lines), which are the front of the Rocky Mountains and the foothills, respectively, as well as the Snowbird Tectonic Zone and the Great Slave Lake Shear Zone. For comparison, see the tectonic map (Fig. 2). Push and pull along model sides and the allowed lateral motion are indicated by blue arrows and circles, respectively. The mean orientation of S Hmax is indicated by a rose diagram; note that stress orientation is parallel and orthogonal, respectively, to the model box. The trace of the cross section in Fig. 3 is indicated by the green line. potential to detach the stresses of the supra-salt units from the sub-salt units. Furthermore, the Snowbird tectonic zone and the Great Slave shear zone are incorporated, and cut the basement and the sediments.
The basement top ( Fig. 6) is defined as the boundary between the basement (southwestern continuation of the Canadian Shield) and the sedimentary basin. It is constructed by DSI, based on 7257 well data available from the Alberta Geological Survey (AGS).
The strata overlying the Precambrian basement are subdivided into the Rocky Mountains, the foreland fold-and-trust belt (foothills) and the sediments within the basin (Fig. 7). The first set consists of allochthonous Palaeozoic strata, whereas the foothills consists of the same stack of sediments from the entire Phanerozoic, deposited in the Western Canadian Sedimentary Basin. For a definition of the boundary between these parts, geological maps and interpreted cross sections were used (Mossop and Shetsen, 1994b;Wright et al., 1994).
Salt deposits within the sedimentary column have the possibility to geomechanically detach the stresses in the upper rock units from the long wave-length stresses at depth (Bell, 1993;Roth and Fleckenstein, 2001;Röckel and Lempp, 2003;Tingay et al., 2005). The Devonian Elk Point Group contains several salt deposits (Grobe, 2000;Meijer Drees, 1994). These up to 380 m thick deposits (Figs. 3 and 7) have been recognized within five stratigraphic formations. They are ordered from oldest to youngest: Lower Lotsberg, Upper Lotsberg, Cold Lake, Prairie Evaporate and Hubbard Evaporate salts. The Elk Point evaporates are separated from the other basin sediments (Mossop and Shetsen, 1994a, Chapter 8-26 ) as an independent unit; evaporate strata with a thickness of ≥ 100 m are used based on data from Grobe (2000). All these interfaces are also generated with the DSI algorithm. Finally, the model box is completed with the digital elevation model (DEM) from the USGS (2008).

Model discretization into finite elements
Our key goal is to model the contemporary 3-D stress state within the basin and in the upper part of the basement. To reproduce the thin rock salt layer within the basin, it was necessary to have a minimum element amount of six elements within the basin in the z direction. This results in a vertical resolution of about 200 to 800 m for the upper model parts (Fig. 8). In the x and y directions, the resolution within the basin and the upper basement is about 5000 m. The element thickness increases with depth within the basement. In deeper parts (−25 to −80 km), the resolution is approximately 20 km in all directions (Fig. 5)  functions. The partial differential equation of the equilibrium of forces is solved numerically using the Abaqus ® /Standard v.6.11 finite element software.

Rock properties
To calculate the stresses, Young's modulus (E) and Poisson's ratio (ν) are the essential geomechanical material properties. The body forces of the rock units are represented by the density (ρ). Mantle density below Alberta ranges from 3346 to 3366 kg m −3 according to White et al. (2005). For this model, a density of 3350 kg m −3 for the mantle is used ( Table 1). The density of the Canadian Shield ranges from 2640 to 2830 kg m −3 , with this model using a value of 2800 kg m −3 . Young's modulus and Poisson's ratio of the basement are calculated based on the V p and V s data from northern Alberta (Dalton et al., 2011). The dynamic Young's modulus and Poisson's ratio (0.21-0.22) are calculated according to Mavko et al. (2009). Based on the dynamic Young's modulus, the static Young's modulus is calculated according to King (1983) and Wang et al. (2000), with a range of 1.02×10 10 to 8.56×10 10 Pa. In the model, 0.21 and 7.0×10 10 are used as Poisson's ratio and Young's modulus respectively for the basement, which is in agreement with data from Turcotte and Schubert (2002). Most Phanerozoic sediments overlying the basement, including the foothills and the Rocky Mountains, are mainly clastic sediments (e.g. sandstone or shale) and limestone, with the exception of the separated evaporates. These material properties are estimated based on Fossen (2010), Okrusch and Matthes (2005) and Turcotte and Schubert (2002); see Table 1.

Initial stress state
Deformation of the model due to gravity-driven subsidence is not desired. Therefore, an initial stress state of the model is derived, which is in equilibrium with the body forces (gravity). For the initial stress state, uniaxial strain conditions (Eq. 2) or lithostatic stress conditions for greater depths (Heim, 1878, Eq. 3) are often assumed.
Using uniaxial strain conditions (k = 1/3 when ν is 0.25, Eq. 2) or lithostatic conditions (k = 1, Eq. 3), the stress ratio k (Eq. 4) is constant for both, when plotting vs. depth, but when k is plotted vs. depth, based on in situ data, the discrepancy is obvious (e.g. Brown and Hoek, 1978;Gay, 1975, Fig. 9a). Visible are increasing k values close to the surface. Thus, assuming lithostatic or uniaxial conditions is apparently insufficient for appropriate initial stress conditions. Sheorey (1994) provides a simple spherical earth model for tectonically calm regions with no significant lateral density and strength contrasts. In this model, k is a function (Eq. 5) of Young's modulus (E in GPa) and depth (z in m). This was confirmed by later published in situ stress magnitudes from the KTB borehole (Brudy et al., 1997, see Fig. 9a).
When the model is embedded in an extended model with inclined edges, it is possible to find a fit of the k values vs. depth from the model. This is in comparison to a synthetic depth distribution, based on the Sheorey equation (Eq. 5). This technique has so far been used only occasionally (Buchmann and Connolly, 2007;Hergert and Heidbach, 2011). By generating the initial stress model, settlement due to the gravitational load occurs. Using the initial stress condition in the stress model, settlement in the model (< 1 m) can be neglected in relation to the model size.
For comparison, the calculated k ratio, Eq. (5) from Sheorey (1994), and the initial k ratio from the model of Alberta, are plotted vs. depth. Material properties are adjusted for the initial model only until good agreement is obtained. Exemplarily, two of them are plotted for illustration ( Fig. 9b and c). From a purely technical point of view, the initial stress conditions were determined after calibration of the used sediment density. Henton et al. (2006) and Mazzotti et al. (2011) showed that surface strain measured by GPS indicates strain rates are below the measurement error within Alberta and the Rocky Mountains. More to the west, in the Intramontane Belt, the values are also very low, yet in the coastal cordilleras, rates of about 10-15 mm yr −1 in the northeasterly direction with respect to stable North America are observed. The North American Eulerian rotation pole is located southwest of Ecuador, resulting in an anticlockwise rotation of about 20 mm yr −1 in the southwesterly direction in Alberta (Henton et al., 2006). Flesch et al. (2007) found that (deviatoric) stresses associated with the accommodation of relative plate motion are of the same order of magnitude as buoyancy forces (gravitational potential energy -GPE). The orientation of observed North American rotation, shortening in the Canadian Cordillera (Henton et al., 2006;Mazzotti et al., 2011), and GPE gradient orientation (Flesch et al., 2007) correspond to the observed average S Hmax azimuths in Alberta (see the rose diagram in Fig. 4).

Boundary conditions
As the model edges are parallel and perpendicular, respectively, to the observed plate motion, GPE and horizontal stress azimuth, displacement at the model boundaries will be applied orthogonally to the side walls of the model box. Horizontal and vertical motion is allowed along the side walls (Fig. 4). The applied amount and orientation of push (towards the northeast) and pull (towards the southeast) along the model will be tested during the calibration phase of the model. The bottom of the model is fixed in the z direction; lateral motion within the extent of the model box is allowed.

In situ stress
This section presents a short introduction to the terminology used for the stress data during the model calibration procedure.

Orientation and magnitudes of stresses in sedimentary basins
The 3-D stress in rock (σ ) is described with a second-order tensor. By choosing an principal coordinate system, the stress tensor (σ ij ) can be expressed with the three principal stresses: These act normally to the principal planes and are the following: σ 1 > σ 2 > σ 3 , in the order of magnitude. As the earth's surface is a free surface, and sedimentary basins are roughly flat at the top, it is often assumed that the vertical stress (S V ) is a principal stress. With this assumption, the minimum horizontal stress (S hmin ) and the maximum horizontal stress (S Hmax ) (e.g. Jaeger et al., 2009;McGarr and Gay, 1978;Schmitt et al., 2012) are also principal stresses that are orthogonal to each other (Fig. 10). Their relative magnitudes determine the stress regime (Anderson, 1951, cited in Kanamori andBrodsky, 2004): -normal faulting: S V > S Hmax > S hmin -strike slip: S Hmax > S V > S hmin -reverse faulting: S Hmax > S hmin > S V .

Contemporary stress field in the Alberta Basin
The present-day stress field in Alberta has been the subject of several studies. It started with Bell and Gough (1979) recognizing in the Alberta Basin that borehole breakouts are an indicator of crustal stress orientation (Fig. 10). They found that the S Hmax azimuth is uniformly oriented southwest to northeast in substantial parts of the Alberta Basin (Fig. 11). This observed orientation is perpendicular to the Rocky Mountain trench, which was confirmed by Adams and Bell (1991), Bell and Gough (1981), Bell et al. (1994), Fordjor et al. (1983) and recently by Reiter et al. (2014). Orientation data are derived from a large variety of rock types, depths, and different indicators. These are borehole breakouts at a depth range of 113-5485 m (e.g. Bell et al., 1994), geological indicators (Bell, 1985), and drilling-induced tensile fractures , and seismological studies in the Canadian Cordillera (Ristau et al., 2007) confirmed the overall orientation pattern . Only an anticlockwise rotation of about 10-20 • is observed in northern Alberta over the Peace River Arch. The same homogeneous stress orientation is observed over wide areas of the North American plate (Bell and Gough, 1979;Adams and Bell, 1991;Fordjor et al., 1983;Gough et al., 1983;Reiter et al., 2014;Sbar and Sykes, 1973;Zoback and Zoback, 1980), which indicates that southwestto-northeast stress orientation is present over the whole lithosphere rather than sediments only . This implies also that the sediments are attached to the basement (Bell, 1996b). The S Hmax orientation is at a right angle to the Rocky Mountains fold axis. Therefore, the stress field responsible for thrust faulting in Mesozoic time is still present (Bell and Gough, 1979). The driving force of the observed stress pattern is plate tectonics, either by drag resistance of the lithosphere sliding over asthenosphere (Bell and   lan, 1995; Zoback and Zoback, 1980) or mantle convection propelling the lithosphere (Bell and Gough, 1979;Fordjor et al., 1983;Gough, 1984). The depth gradients of S V and S hmin increase from the basin centre towards the foothills and the Rocky Mountains (Baranova et al., 1999;Bell, 1996b;Bell and Bachu, 2004;Bell and Grasby, 2012). This trend coincides with higher organic maturity (England and Bustin, 1986;Nurkowski, 1984) and larger compaction (Bell and Bachu, 2004) Figure 11. Crustal stress map of Alberta. Lines represent orientations of maximum horizontal compressional stress S Hmax ; line length is proportional to the data quality. Colours indicate stress regimes, with green for strike-slip faulting (SS), blue for thrust faulting (TF), and black for unknown regime (U). In summary, there are 321 S Hmax azimuth data available within the modelled region. Data are from the latest update of the Canadian stress map (Reiter et al., 2014). direction, which is related to the depth of present and past burial. The maximum erosion of basin sediments is by about 1400 m (Woodland and , uplift occurring since the mid-Cenozoic time, mainly in the foothills (Bell and McLellan, 1995).
The stress regime in the basin sediments changes from thrust faulting in the foothills to strike slip within the basin, up to a normal faulting regime further east in Saskatchewan (Bell and Gough, 1979;Bell et al., 1994;Bell and McLellan, 1995;Bell and Bachu, 2003;Bell and Grasby, 2012;Woodland and Bell, 1989). A similar change from surface to depth is observed: from thrust faulting at < 350-600 m in depth, strike slip in a depth range of about 500-2500 m, down to normal faulting at greater depths > 2500 m (Bell and Babcock, 1986;Fordjor et al., 1983;Jenkins and Kirkpatrick, 1979). There is also a varying S hmin gradient discussed (Bachu et al., 2008;Bell and Grasby, 2012;Hawkes et al., 2005), but this is may be due to different measurement methods (Bell et al., 1994) or man-made stress changes. The S Hmax /S hmin ratio in the Alberta Basin is about 1.3-1.6 .
Man-made stress perturbation due to hydrocarbon production or acid gas injections (e.g. Bachu et al., 2008;Bell and Grasby, 2012;Woodland and Bell, 1989) reduces or in-creases reservoir fluid pressure respectively, but has likely only local effects (e.g. Altmann et al., 2010). Furthermore, Baranova et al. (1999) found a strong correlation between rates of gas production and the number of seismic events, which is reasonable because production lead to decrease of S V and increase of S Hmax -consequentially increasing differential stresses. The stress change due to the gas extraction point to a regime, which favours thrust faulting (Baranova et al., 1999). Hydraulic fractures applied for hydrocarbon industry or for enhanced geothermal systems deeper than 350 m will open parallel to southwest-to northeastoriented S Hmax orientations, except that in the Peace River Arch, they will tend to south-soutwest to north-northeast (Bell et al., 1994;Bell and Grasby, 2012). Close to the Rocky Mountain foothills, northwest-to southeast-oriented hydraulic fractures are possible parallel to the thrust planes and the fold axes (Bell and Babcock, 1986). However, horizontal wells e.g. for EGS should be designed parallel to the S hmin orientation (Bell and Grasby, 2012).

Vertical stress (S V )
The vertical stress (S V ) is the overburden load, which is estimated using density logs (e.g. Gardner and Dumanoir, 1980) in a well: For the Alberta model region 981 S V magnitude data sets are available (provided by the AGS), these are indicated by black points in Fig. 12. S V magnitude data vary only slightly, even in greater depths, the lateral variation is less than 5 MPa.

Orientation of maximum horizontal stress (S Hmax )
The orientation of S Hmax is indicated by borehole breakouts, focal mechanisms, hydraulic fracturing, overcoring, and drilling-induced fractured and geological indicators (for an overview, see Bell, 1996a;Ljunggren et al., 2003;Schmitt et al., 2012;Zang and Stephansson, 2010;. 321 S Hmax azimuth data sets are available for the modelled region in Alberta; these are indicated in Fig. 11, based on the latest update of the Canadian stress database (Reiter et al., 2014).

Magnitude of minimum horizontal stress (S hmin )
The S hmin magnitudes are measured by hydraulic fracturing or the similar leak-off test. During hydraulic fracturing (Bell, 1996a;Haimson and Cornet, 2003;Hubbert and Willis, 1957; and leak-off tests (e.g. Li et al., 2009;White et al., 2002;Zhou, 1997), the down-hole pressure is increased up to pressure loss due to fluid leakage in the  Figure 12. Depth plot of the in situ stress magnitudes. These are 981 S V magnitude data (black), 1720 S hmin magnitudes (several colours), and 2 measured as well as 11 calculated S Hmax magnitudes (highlighted red points). S hmin data are colour coded depending on the test type, which is taken over from the original database. rock mass. This happens, when the hydraulic fracture splits apart the surrounding rock perpendicular to the least principal stress (σ 3 ), usually assumed to be S hmin in sedimentary basins, and therefore the fracture opens in S Hmax orientation (Fig. 10). The highest pressure is the fracture breakdown pressure (FBP, Haimson and Fairhurst, 1969), which is S hmin + rock resistance up to failure. When the pressure at which the fracture closes or re-opens is less than S V , it is assumed that S hmin is measured (Haimson and Fairhurst, 1969). The mini-frac test (e.g. McLellan, 1987;Woodland and Bell, 1989) and the micro-frac test ) as hydro-fracturing methods estimate the closure pressure by opening and closing the fracture several times, but differs by the injected fluid volume. The term "leak-off tests" is used variably, and can be distinguished by their aim into formation integrity tests (FIT), "classic" leak-off tests (LOT) and extended leak-off tests (XLOT) (White et al., 2002). The general method is similar, but differs in pumping cycles and the point at which the pumping is stopped. Usually, leak-off tests (LOT) are meant, which provide the upper limit of S hmin and measure the fracture closure pressure (FCP) or the instantaneous shutin pressure (ISIP) (White et al., 2002). Extended leak-off tests (XLOT) allow measurement of the fracture re-opening pressure, like the original hydro-fracture tests.
For the model region, 1720 S hmin magnitudes data are available, provided by the AGS; see Fig. 12. These are 784 leak-off magnitude data and 936 magnitude data from hydraulic fracturing. The different hydraulic fracturing methods are 14 Micro-frac, 91 Mini-frac, 250 Hydro-Frac-AIP, and 581 Hydro-Frac-FBP data. The data scatter strongly, independently from the test method or lithology. Further detailed information about the measurements are not available; that would allow whether the data represents the undisturbed stress state or not. The scatter either reflects the spatial anisotropy of the in situ stress or that the data set is noisy, i.e. a mix of in situ stress information and data from areas with a disturbed stress field.

Magnitude of maximum horizontal stress (S Hmax )
The magnitude of S Hmax is measured via the overcoring method (McGarr and Gay, 1978;Obert, 1962), which isolates a rock cylinder from the surrounding rock and measures the elastic relaxation of the rock cylinder. This is equivalent to the stress magnitude as well as the stress orientation, before removal of the surrounding rock. The drawbacks are the small quantity of inspected rock mass and the fact that the application is usually close to the surface. Furthermore, there are several methods used to calculate S Hmax , based on S hmin magnitudes and known rock properties (e.g. Schmitt et al., 2012). For the model region, 11 calculated data (Bell et al., 1994) and 2 shallow measured data (overcoring from Kaiser et al., 1982) are available (see Fig. 12).

General comparison technique
The 3-D geomechanical-numerical model (with the initial stress state) from Alberta will be calibrated in the following according to the work flow scheme (Fig. 1). Each type of in situ stress data will be used step by step to calibrate the model. We first use the S V data to calibrate the density (technically, the initial stress state is found after this step); then, we use the S Hmax azimuth data to calibrate the orientation of applied displacement boundary conditions. Finally, the S hmin and S Hmax magnitudes are used to calibrate the magnitude of applied displacement boundary conditions, i.e. push and pull at the edges of the model box.
In each step, the modelled stress tensor is interpolated via inverse distance interpolation onto each point, where in situ stress data are available. The difference ( S) between measured stress (S Measured ) and the modelled stress (S Model ) is calculated in the following way: deviation is normalized by the modelled stress value: To evaluate the differences between each in situ data set and the model as a whole, the median of S ( S) is calculated as a single value for each model. In the case of the best-fit model, the S shall be close to zero. To estimate the influence of outliers and the variation of the data, the mean ( S) and the standard deviation (SD) is also calculated. The linear correlation between the in situ data and the model data is represented by the Pearson product-moment correlation coefficient (r), where r = 1 indicates total positive correlation, r = −1 total negative correlation and r = 0 no correlation. The data sets are contaminated with unlikely in situ data; such data are often sorted out (e.g. Bell and Bachu, 2004;Bell and Grasby, 2012) for interpolation. As statistical tests are used in this study, data weed out is not required.

Calibration of material density on S V data
The density of the sedimentary basin is calibrated based on S V magnitudes (n = 981, Fig. 12); all other material properties are defined in Sect. 3.2. The overall density of the modelled basin sediments is tested with 2200, 2300 and 2400 kg m −3 . According to Eq. (9), the difference between the measured and modelled S V ( S V ), as well as the normalized S V (n S V ), are calculated: The model with a density of 2200 kg m −3 has a S V which is close to zero (−0.09 MPa) and a S V of 0.28 MPa, which is also close to zero (Fig. 13). A standard deviation of 5.58 MPa as well as the Gaussian distribution in the normalized histogram (Fig. 14a) indicate that there is no large data drift. The correlation coefficient of r = 0.935 indicates a good fit.

Calibration of the orientation of displacement boundary conditions based on S Hmax azimuth data
The S Hmax orientation is controlled by the applied displacement boundary conditions. Within the model region, 321 S Hmax orientation data are available. They are displayed together with the data aside of the model in the stress map from Alberta (Fig. 11). The observed stress pattern is quite homogeneous.
The displacement boundary conditions act orthogonally to the model margins, in a horizontal direction. Whereas shortening is applied in a northeasterly direction to the model, extension is applied in a southeasterly direction (Fig. 4). According to Eq. (9), the difference between the measured and The histogram of the S Hmax azimuths (Fig. 14b) displays a main cluster around zero with a S Hmax azimuth of −3.44 • . The main cluster ranges from −40 • to 40 • ; a second (smaller) cluster ranges from −50 • over ±90 • to 60 • , with a slight peak around ±90 • . This is exactly orthogonal to the main cluster and explains the large SD of 26.9 • ; the mean ( S Hmax = 3.42 • ) is similar to the median. The best-fit orientation is found for a large range of push and pull magnitudes. Therefore, different oriented boundary conditions are not further tested. The normalized S Hmax magnitude histogram for the 13 available data displays two negative outliers (the two shallow (overcoring) data). However, the 11 calculated data are arranged around zero.

Calibration of the magnitude of displacement boundary conditions by S hmin and S Hmax magnitude data
The S hmin (n = 1720) and S Hmax (n = 13 : 2 measured and 11 calculated) magnitude data (Fig. 12) are used to calibrate the magnitude of applied push and pull along the model edges (Fig. 4). The aim is a model, which mimics the S hmin and S Hmax in situ magnitude data quite well. Several scenarios with different amount of push and pull are calculated, to estimate the range of push and pull, close to the best-fit model. In the following we focus in only four scenarios with a different amount of push and pull (Table 2). According to Eqs. (9) and (10) the difference between the measured and modelled S hmin and S Hmax magnitudes are calculated as well as the normalized difference.
The calculated S hmin and S Hmax of four model runs (Table 2) are plotted in the push vs. pull diagram ( Fig. 15a and b). To highlight the linear dependency between push and pull in an elastic model, colour coded isolines are plotted. Each model along the light blue line (Fig. 15a) would derive a model, which fits well with the in situ S hmin data. The same stands for the light blue line in Figs. 15b and S Hmax data. As the determination of the best-fit model is intended, the intersection of both light blue lines from Fig. 15a and b would derive such a model. This is done with a bivariate linear regression based on the spatial distribution of the S hmin and S Hmax (Fig. 15c). This method Table 2. Overview of major push-and-pull experiments. The orientation of the displacement boundary condition is indicated in Fig. 4. Four test scenarios with different push and pull magnitudes are displayed, same as in Fig. 15a  provides the following equations, which describes the zero isoline (light blue line) as a linear function. These are for median S hmin = 0: and for the S Hmax zero isoline: the best-fit model has a push from southwest of 86.24 m and a pull in the southeasterly direction of 194.52 m ( Table 2,  the mean values ( S hmin =0.03 and S Hmax = −2.76).
The distribution of the normalized S hmin (Fig. 14c) displays a positive screwed distribution. The correlation coefficient of the in situ S hmin magnitude data and the modelled S hmin is r = 0.835. The normalized best-fit of S Hmax (Fig. 14d) displays two outliers; these are the only two measured S Hmax magnitudes, measured at a depth of 152 m.

Workflow and calibration
The general workflow of model calibration (Fig. 1) is similar to other studies on numerical stress field modelling (e.g. Buchmann and Connolly, 2007;Fischer and Henk, 2013;Heidbach et al., 2013;Hergert and Heidbach, 2011). However, in contrast to former studies, the number of in situ stress data from the Alberta Basin allows a statistical comparison with the model results.

S V calibration
Given the large model size (1200 km × 700 km × 80 km), the Alberta Basin infill is considered to be one material type only, except for the Elk Point evaporates. The available data set of 981 S V magnitude data points (Fig. 12) could be used via linear regression to calculate the overall density of the sedimentary basin, but to incorporate the (minor) lateral effects of topography, the overall density is determined in the calibration process. Plotting the distribution of the normalized deviation (n S V ) in histogram Fig. 14a demonstrates a Gaussian distribution, which implies there is no process affecting data drift. As data spreading does not depend on the vertical depth, a slightly higher lithological resolution with linearly increasing density into depth would most likely not deliver a much better data fit. This could be solved by incorporation of all stratigraphic units, which would go far beyond the goals of this study. The spatial plot of S V (Fig. 16) shows that in situ S V magnitudes are slightly higher close to the foothills. This is expected from former studies, showing S V increases in the southwesterly direction (e.g. Bell and Bachu, 2004;Bell and Grasby, 2012).

Calibration of S Hmax orientation
The 321 data records of the S Hmax azimuth (Fig. 11) are used to test the orientation of the applied displacement boundary conditions. As long as a certain push to the northeast and a pull in the southeasterly direction are applied orthogonally to the model box (Fig. 4), a good fit of the stress orientation ( Fig. 14b and 17) is achieved. No variation of the boundary conditions (orientation of push and pull) was necessary, due to the appropriately chosen model orientation.
The histogram of the S Hmax azimuth (Fig. 14b) displays two data clusters. The larger cluster displays a normal distribution around zero, which is confirmed by a S Hmax azimuth of −3.44 • . A second data cluster is distributed around ±90 • , which explains the high SD of 26.9 • . The second cluster with a deviation of around ±90 • is the orientation of S hmin indicating that some of the used in situ S Hmax azimuth data, are likely miss-interpreted S hmin orientations. Such incorrect interpretations are sometimes observed in borehole breakout data (Brudy and Kjørholt, 2001;Barton and Moos, 2010); in such cases drilling induced tensile fractures, originated during drilling, are misinterpreted as borehole breakouts. Other reasons for orientation of data with right angles to the major population (S hmin ) are mud cake padding along caved zones and the collapse of pre-existing open fractures, again parallel to S Hmax (Bell and Babcock, 1986;Bell and Grasby, 2012). Therefore, this second (smaller) cluster around ±90 • rather confirms then disproves the good data fit by the model.
The alternative explanation would be a horizontal stress state close to isotropic. This would allow large stress rotation due to small local stress sources (Heidbach et al., 2007). However, from the provided data, this explanation can be ruled out, and the previously stated explanation is much more likely.
There are two areas in the modelled region, where a systematic difference of the S Hmax azimuth between the in situ data and the model is visible (Fig. 17). These are the western Peace River Arch (56 • N, 118 • W) and the Sweetgrass Arch in the very south, close to the southern Bow Island Arch (50 • N, 114 • W). In the Peace River Arch, in situ S Hmax is rotated by about 20-30 • anticlockwise (Bell et al., 1994;Bell and Grasby, 2012). The causes of this rotation, (see discussion in Bell and Babcock, 1986;Bell and McCallum, 1990;Bell, 1996b;Bell and Grasby, 2012;Bouzidi et al., 2002;Dusseault and Yassir, 1994;Eaton et al., 1999;Halchuk and Mereu, 1990), are not well represented by the model. Bell and Gough (1979) suggested that the S Hmax orientation is orthogonal to the topography of the Rocky Mountains, but a comparison close to the topography in the very south of Alberta displays a good fit between the in situ data and the model. They appear more influenced by the overall orientation than by the topography, which is also found by Reiter et al. (2014).
In contrast, the clockwise rotation of about 25 • with respect to the regional trend is obvious close to the Bow Island Arch. Likely this systematic rotation is caused by structural features along the Bow Island Arch which are not incorporated into the model, then by the Rocky Mountain topography.
The Bow Island Arch separates the Alberta Basin and the Walliston Basin. It is a northeastward plunging Precambrian basement feature, which was activated during the Laramiden orogeny and which may be associated with intrusions, similar to Eocene intrusions, about 200 km to the south in Montana (Podruski, 1988). The systematic S Hmax rotation in that region is most likely affected by these basement features along the Bow Island Arch.

S hmin calibration
The largest number of stress data are the S hmin magnitudes (n = 1720, Fig. 12). These, with a few S Hmax magnitudes (n = 13 − 2 measure and 11 calculated, Fig. 12), are used to find the best-fit magnitudes of the utilized boundary conditions.
A large number of models were tested, but only four of these are shown here (Table 2 and Fig. 15a and b). Based on these four test scenarios with different strain magnitudes, the best-fit model is determined via bivariate regression. Calculated is the intersection of zero-isolines of S Hmax and S Hmin (Fig. 15c) based on a plot of push vs. pull (Table 2, Fig. 15a and b). This is possible as linear elastic rheology is used in the model.
An evaluation as to whether the measured S hmin magnitudes really represent S hmin or only σ 3 , because hydraulic fracturing tests provide the information on the smallest principal stress. The correlation between the in situ S hmin magnitudes vs. modelled σ 3 (r = 0.837) is negligible higher then vs. modelled S hmin (r = 0.835). This indicates that in situ S hmin measurements in the Alberta Basin likely represent the magnitude of S hmin and confirms the assumptions for sedimentary basins, being that S hmin is the smallest principal stress (e.g. Jaeger et al., 2009;McGarr and Gay, 1978;Schmitt et al., 2012).
The spatial distribution of the S hmin (Fig. 18) indicates that larger differences between the in situ and modelled magnitudes mainly occur in regions with clustered data. Slightly higher in situ magnitudes are observed in the region 56 • N, 121 • W, in contrast to the slightly lower in situ S hmin magnitudes in region 55 • N, 119 • W.
The n S hmin histogram in Fig. 14c displays a positive screwed distribution. This indicates that the model underestimates a larger portion of the S hmin magnitudes, in contrast to the in situ data.
To examine deviation reasons, Fig. 19 plots n S hmin depending on depth with the measuring method is indicated. In situ S hmin LOT data provide rather positive n S hmin values in shallow depths (< 500 m) and negative values in depth > 500 m relative to the modelled S hmin magnitudes. This implies that the model derives smaller magnitudes compared to shallow LOT magnitudes and larger ones with respect to deeper LOT magnitudes. In contrast hydraulic fracturing data did not indicate systematic deviations.

Deviation of Leak-off test (LOT) data vs. stress model
There are several reasons for the discrepancy in S hmin in the in situ LOT data vs. stress model; they are -the thrust faulting regime (S hmin > σ 3 ), -systematic measurement errors, -systematic model errors, and -disturbed in situ measurements.

Thrust faulting regime (S hmin > σ 3 )
Hydraulic fracturing (HF) and leak-off tests (LOT) measure the smallest principal stress (σ 3 ), which is expected as S hmin in sedimentary basins with normal faulting or strike slip stress regime. In a thrust faulting stress regime, in situ data would underestimate S hmin , as S V is measured. However, the thrust faulting stress regime is expected in the Rocky Mountains, as well as in and close by the foothills (Bell and Gough, 1979;Bell et al., 1994;Bell and McLellan, 1995;Bell and Bachu, 2003;Bell and Grasby, 2012;Woodland and Bell, 1989) and in shallow depths (up to 600 m) (Bell and Babcock, 1986;Fordjor et al., 1983;Jenkins and Kirkpatrick, 1979). Measuring S V = σ 3 at shallow depths (< 500 m) instead of S hmin would indicate a thrust faulting regime for some regions, but shallow LOT magnitudes are systematically larger than expected by the stress model, which excludes this attempt at explanation. *

Systematic measurement errors
Alternatively, an overestimation of S hmin by LOT could be explained when the formation breakdown pressure (FBP) or leak-off pressure (LOP) (White et al., 2002) is measured. Additionally, LOT performed at shallow depths (< 300 m) are less reliable, because the tensile strength of the rock plays a more significant role for the measured pressure (Bachu et al., 2008). These reasons would explain larger S hmin magnitudes from LOT, compared to the model at shallow depths (< 500 m). LOT are also used as formation integrity tests (FIT) to determine whether the well bore can sustain the stresses expected during drilling and production, and then to determine . The distribution of the normalized S hmin vs. depth, the measurement method is colour coded. Results close to the surface up to about −500 m (indicated by greyish haze) have to be interpreted with care, as interpolation of the integration points to the nodes of the finite elements at the surface is problematic. Note that shallow (< −500 m) leak-off tests (LOT) deliver systematic higher magnitudes than the stress model. In contrast to that are deeper (> −500 m) LOT data, which have systematically smaller magnitudes as the model. Hydraulic fracturing (HF) data are unconcerned from such systematic drift.
stress magnitudes (e.g. White et al., 2002). Such FIT data derives smaller magnitudes then the formation S hmin magnitude. Furthermore, poor cement seal between the well bore and the casing close to the LOT can reduce measured magnitude (Edwards et al., 1998). Both reasons could explain smaller S hmin measurements in greater depth (> 500 m), derived from LOT magnitudes. *

Systematic model errors
In situ stresses are affected by the lithology at the locality (Roche et al., 2013;Warpinski, 1989). This is observed in the Alberta Basin, where sandstone exhibits lower in situ S hmin magnitudes than shale (Bell and Grasby, 2012;Kry and Gronseth, 1983). As the modelled basin has only one material property, likely the evaporate layer, the modelled stresses represent stress conditions from rocks where in situ stress data are derived. These are mainly sandstone and limestone, whereas leak-off tests are usually conducted in shale (Bell and Grasby, 2012). Therefore, the drift could be explained by shallow shale and deeper sandstone or limestone. Furthermore, extrapolation drift close to a free surface of an FEM model are sometimes observed. *

Disturbed in situ measurements
Man-made stress changes perturb the juvenile in situ stress due to production (Bell and Grasby, 2012) and injection of fluids (Bachu et al., 2008), besides other mining activity. This is obvious where induced seismicity has been reported (e.g. Baranova et al., 1999;Schultz et al., 2014). Such effects are not restricted to the reservoir alone; the country rock is affected too, in various styles, depending on the relative position (e.g. Segall, 1989).

S Hmax calibration
The quantity and quality of the S Hmax magnitude data are rather poor compared to the S hmin magnitude data, but are very helpful for constraining the best-fit model. Otherwise, only a linear best-fit function could be estimated. The two outliers (Fig. 14d) are the only two measured S Hmax magnitude data from Kaiser et al. (1982) at a depth of 152 m in clay shale. As the measured in situ data are more reasonable then the modelled magnitudes, the reasons for the large deviation are most likely extrapolation problems close to the model surface, as discussed in the previous chapter.

Reliability of the predicted 3-D stress field
The calibration procedure presented in Sect. 5 ensures that the resulting stress field is from a statistically point of view the best-fit model. Nevertheless, there are three main issues that limit the reliability of the model output: (1) The size of the model limits the local resolution. Local structural features such as small faults, material inhomogeneities or stratigraphy are not represented in the model. Thus the predicted stress field in particular close to local structures cannot resolve potential local stress field variations. However, the objective of the large-scale model is to provide boundary and/or initial conditions of regional to local stress field model that address the local stress field variations (e.g. Heidbach et al., 2013). (2) The predicted S Hmax orientations as well as the S V magnitudes are quite reliable given that density is quite well known (on the level of the resolution of the model) and due to the fact that the wave-length of the S Hmax orientation is high and only minor affected by local features (Reiter et al., 2014). Also, the reliability of the predicted S hmin magnitudes is high, due to the large number of available data in the model volume. However, to what extend the S hmin magnitudes that we use for the calibration of the undisturbed in situ stress are unaffected by production is not clear. Nevertheless, assuming that the majority of these data reflect the in situ stress, the outliers of S hmin magnitude data that are affected by production are weeded out by the statistic tests. (3) The least reliable stress component of the predicted stress field is the S Hmax magnitude. Only 13 S Hmax magnitude data were available for the model calibration and these data have much higher inherent uncertainties in comparison with the S hmin magnitude data. Given this, the statistical fit to the S Hmax magnitude presented in Figs. 14 and 15 is not well constrained.
In particular in areas where S Hmax is close to the S V magnitude the stress regime is not resolved and might change from strike-slip to normal faulting. This is critical for e.g. borehole stability or for the assessment of potential fault reactivation (Fuchs and Müller, 2001;Moeck and Backers, 2011). This is not a problem of the presented model in particular, but is a key challenge in general as this stress tensor component is hardly constrained from model-independent data.

Impact of fault activation
The Great Slave Lake Shear Zone (GLS) and the Snowbird Tectonic Zone (STZ) are incorporated within the basement and the basin as vertical contact surfaces. The contact between the Alberta Basin and the foothills (foothill front), as well as the contact of the foothills to the Rocky Mountains (Rock Mountain front), are defined in the model as contact surfaces too. During the model calibration, all these contact surfaces are handled as locked faults with a high friction coefficient. To test the impact of fault re-activation on the stress field, we use in a model variant friction coefficients of 0.3 for STZ and GLS within the basement. For the activation of the basement tectonic zones, the found correlation coefficient for S hmin has been lowered only slightly: STZ alone (r = 0.808), GLS alone (r = 0.828) and STZ together with GLS (r = 0.801), compared to the best-fit model (r = 0.835). When the friction is lowered at the foothill front (r = 0.836) and the Rocky Mountains front (r = 0.835) alone, the correlation coefficient did not change. Only when both, the foothill and the Rocky Mountains front are active, the correlation declines (r = 0.701).
The S Hmax orientation changes slightly (up to 2 • ) for all the fault activation. The exception is the Rocky Mountains front, where the S Hmax orientation is equal to the best-fit model. This is expected, as only a few S Hmax indicators are derived close to the Rocky Mountains front.

Impact of Moho depth variation
To test the influence of the Moho topography (Fig. 5) on the stress state within the Alberta Basin, the best-fit model is modified. The Moho depth is uniform (z = −50 km) over the entire model region. The results of S hmin magnitudes show that this model fits all data, similar to the best-fit model (r = 0.835 for both model runs). S Hmax orientation did not change between the models. Probably, stress magnitudes and orientations are only slightly influenced by the Moho topography in this region.

Model application for deep geothermal reservoirs
To generate electricity, water with temperatures of 120-150 • C is needed. This requires well depths of 4000-6000 m in Alberta (Majorowicz and Grasby, 2010a, b). However, stimulation is required to enhance permeabilities (enhanced geothermal systems -EGS) at such depths. Furthermore, less hot water has potential as a domestic heat source.
Three major issues for application of EGS are related to the crustal stresses. These are (1) the orientation of S Hmax , as induced fractures open parallel to S Hmax . This is important for configuration of injection and production wells. (2) The tectonic stress regime determines whether fractures open horizontally (thrust faulting) or vertically (strike slip and normal faulting regime). Furthermore, well stability is a major issue for deep and of course very expensive wells; thus, (3) stress orientation and magnitude as well as differential stresses are important parameters for safe drilling.
S Hmax orientations have been well understood in the Alberta Basin for decades (e.g. Bell and Gough, 1979;Bell et al., 1994;Fordjor et al., 1983;Reiter et al., 2014), mainly homogeneous in southwest-to-northeast orientation, with the exception of the Peace River Arch, and close to the Bow Island Arch.
Horizontal wells, oriented parallel to S hmin (southeast to northwest) (Bell and Grasby, 2012) with multiple fractures, open in the S Hmax direction. This creates several fluid propagation paths, which has the potential to provide cost-efficient energy (Hofmann et al., 2014).
We chose three locations in Alberta that have been identified as possible geothermal sites Pathak et al., 2013), and show S V , S hmin and S Hmax along the virtual well path (Fig. 20). They are (1) the Edson-Hinton region, with sediment thicknesses of 4000-6000 m and potential temperatures of 100-150 • C, (2) the village of Leduc, 30 km south of Edmonton, with the potential to heat houses , and (3) the Hunt well site (e.g. Majorowicz et al., 2012) close to the town of Fort McMurray, where heat is needed for industrial application (Pathak et al., 2013). The virtual wells in Fig. 20 show that thrust faulting regimes occur close to the Rocky Mountains and foothills at shallow depths. A strike-slip regime is common within the basin from the surface to about 1500 to 3000 m in depth and along the foothills from about 1000 to 4000 m in depth. At greater depths, from about 4500 m in depth for the foothills, 3000 m in depth for Edmonton, and 2000 m in depth for Fort McMurray, a normal faulting regime is expected. This confirms Bell and Gough (1979), Bell et al. (1994), Bell and McLellan (1995), Bell and Babcock (1986), Bell and Bachu (2003), Bell and Grasby (2012), Woodland and Bell (1989), Fordjor et al. (1983) and Jenkins and Kirkpatrick (1979). Therefore, opening of induced fractures horizontally can be expected only close to the foothills, at depths less than 1000 m.

Conclusions
A large data set of stress orientation and stress magnitude data is used to calibrate a 3-D geomechanical-numerical model of the Alberta Basin, which provides a good firstorder estimation of the contemporary stress tensor. During calibration procedure, the density of the sediments, the orientation of the displacement boundary conditions, and the magnitude of applied shortening of the model along the model boundaries, are calibrated. As linear elastic material properties are used, the magnitude of applied displacement boundary conditions for the best-fit model can be determined by bivariate linear regression. This is based on only three (or more) models with variable boundary conditions. The stochastic verified calibration allows the evaluation of measurement outlier and systematic uncertainties. Variations of the best-fit model suggest that main faults have only local effects on the stresses and that the Moho topography has only a negligible impact on the model results. A systematic drift of S hmin magnitudes from leak-off tests against the stress model is obvious, but may be affected by multiple reasons.
The best-fit model applies for potential EGS reservoir horizontal wells, oriented northwest to southeast. A virtual well path or cross estimation of the full contemporary stress tensors can be provided by the model in advance of any drilling. The model has the potential to derive boundary conditions for local or reservoir models (e.g. Reiter et al., 2013), where petrological and tectonic inhomogeneities could be respected in more detail.