An Analytical Framework for Stress Shadow Analysis During Hydraulic Fracturing-Applied to the Bakken Formation, Saskatchewan, Canada

This paper presents selected results of a broader research project pertaining to the hydraulic fracturing of oil 10 reservoirs hosted in the siltstones and fine grained sandstones of the Bakken Formation in southeast Saskatchewan, Canada. The Bakken Formation contains significant volumes of hydrocarbon, but large-scale hydraulic fracturing is required to achieve economic production rates. The performance of hydraulic fractures is strongly dependent on fracture attributes such as length and width, which in turn are dependent on in-situ stresses. This paper reviews methods for estimating changes to the in-situ stress field (stress shadow) resulting from mechanical effects 15 (fracture opening), poro-elastic effects, and thermo-elastic effects associated with fluid injection for hydraulic fracturing. The application of this method is illustrated for a multi-stage hydraulic fracturing operation, to predict principal horizontal stress magnitudes and orientations at each stage. A methodology is also presented for using stress shadow models to assess the potential for inducing shear failure on natural fractures. The results obtained in this work suggest that thermo and poro-elastic stresses are negligible for hydraulic fracturing in the 20 Bakken Formation of southeast Saskatchewan, hence a mechanical stress shadow formulation is used for analyzing multistage hydraulic fracture treatments. This formulation (and a simplified version of the formulation) predicts an increase in instantaneous shut-in pressure (ISIP) that is consistent with field observations (i.e., ISIP increasing from roughly 21.6 MPa to values slightly greater than 26 MPa) for a 30-stage fracture treatment. The size of predicted zones of shear failure on natural fractures are comparable with the event clouds observed in microseismic monitoring when assumed values of 115°/65° are 25 used for natural fracture strike/dip; however, more data on natural fracture attributes and more microseismic monitoring data for the area are required before rigorous assessment of the model is possible. https://doi.org/10.5194/se-2021-1 Preprint. Discussion started: 13 January 2021 c © Author(s) 2021. CC BY 4.0 License.


Introduction 30
During fracturing operations, the changes in principal stress magnitudes and orientations induced in the reservoir during a given fracture stage may alter the conditions for the subsequent stages, giving rise to the so-called "stress shadow" effect (see Figure 1). Stress shadows can impact the effectiveness of fracturing operations by increasing the injection pressures required to create and propagate fractures, reducing fracture width, and potentially altering fracture trajectory (Nagel et al., 2013; The relatively simple stress shadow formulation developed in this project can be used easily on site to assess the effects of stage spacing on fracture pressure and (potentially) fracture orientation, improving fracture stimulation design.
The methodology developed in this work can also be extended to predict zones of possible rock failure (due to shearing) around a hydraulic fracture. If there is a large-scale discontinuity around the wellbore, failure on that discontinuity due to hydraulic 65 fracture can be assessed and that assessment can mitigate the potential for operational problems due to reactivation. Predicting the failure zone associated with slip on smaller natural fractures around a hydraulic fracture is also valuable for predicting stimulated reservoir volume. As such, the following task was additional undertaken in this work:  Using stress shadow modeling results to assess the potential for induced shear failure on natural fractures, as a means of predicting the dimensions of the microseismic cloud area around a hydraulic fracture. 70

Workflow for Using Analytical Stress Shadow Models in Tandem with a Numerical Fracture Simulator
This section presents equations to analyze and calculate the stress shadow existing at each stage of hydraulic fracturing due to the cumulative effects of the previous stages. The workflow used to implement these stress shown models in combination with a numerical fracture simulator is shown in Figure 3. In this work, hydraulic fracture dimensions and pressures were predicted 75 using the commercial fracture simulator GOHFER (Barree and associates, 2016).
The process begins by simulating the first fracture stage with GOHFER using the initial in-situ stress state as model input. The time interval simulated includes the injection period as well as the ensuing pressure falloff, up to the point in time when the second fracture stage is undertaken. At this point, key attributes of the fracture (length, average width, average pressure) are taken and used as inputs for the stress shadow equations. In turn, these equations provide estimates of stress changes induced 80 at the location of the second fracturing stage, which are superimposed on the original in-situ stress state to yield a revised stress state. This revised stress state is used as input for a GOHFER simulation of the second fracture stage. The workflow is then repeated for each subsequent fracture stage, with the stress shadow effects of all preceding stages accounted for.
GOHFER was chosen for this work as a simulator. The rationale for this choice was based on the fact that this modeling tool strikes a balance between being physically robust (in terms of physical processes directly associated with hydraulic fracture 85 behaviour) while being sufficiently practical and accessible to be widely used in the petroleum industry. As such, the results of this work would be relatable in form and content to industry personnel, and new insights and methods would have greater potential for being adopted in practice.
The following sections present the equations that have been identified to enable prediction of stress shadow are based on mechanical, thermal and poro-elastic mechanisms.    https://doi.org/10.5194/se-2021-1 Preprint. Discussion started: 13 January 2021 c Author(s) 2021. CC BY 4.0 License.

Mechanical Stress Shadow
For the calculation of mechanical stress shadow around a fracture, an improved version of the Pollard and Segall (1987) method was used (Figure 4). Ge and Ghassemi (2008) have also used this technique. The hydraulic fracture is assumed to be a 2-D 105 crack with an average propped width of . The equations from Pollard and Segall (1987) are given in equations 1 to 9, as follows. (1) (2) And ∆ = (∆ + ∆ ) (5) r=√ 1 2 , = 1 + 2 2 and = in-situ shear stress and shear stress on the crack surface (assumed zero because hydrostatic pressure of fracture fluid, inside of the crack, doesn't create shear stress on the crack surface) ∆ , ∆ and ∆ = induced stress in X, Y and Z directions.
Negative values of angles , 1 and 2 (shown in Figure 3) should be corrected by adding 180° (π radians). Values for L and are taken from Gohfer simulation output for the prior hydraulic fracture stages.

Thermo-elastic Stress Shadow
In Most studies in thermo-elastic stresses have been conducted in the field of geothermal production (Ge & Ghassemi, 2008). 115 In the petroleum industry, the thermo-elastic stress shadow is often considered to be insignificant due to the small difference in temperature between the injected fluid and in-situ rock. As a result, it is assumed that temperature is distributed uniformly and elliptically around the fracture and that the rock temperature increases outwards as a fracture propagates and injection https://doi.org/10.5194/se-2021-1 Preprint. Discussion started: 13 January 2021 c Author(s) 2021. CC BY 4.0 License.
continues. The cool region (i.e., rock within the ellipse of semi-axes a0 and b0 in Figure 5) is assumed to have the same temperature as the injected fluid, and the unaffected region (i.e., rock exterior to the ellipse of semi-axes a0 and b0) remains 120 constant at the in-situ reservoir temperature.
For the cool region, thermo-elastic stresses perpendicular ( 1 ) and parallel ( 2 ) to the fracture surface are calculated as follows (Perkins and Gonzales, 1985):

Poro-elastic Stress Shadow
For the flood region zone (i.e., the rock within the ellipse of semi-axes of a1 and b1 in Figure 5), poro-elastic stresses perpendicular (Δ 1 ) and parallel (Δ 2 ) to the fracture surface are calculated as follows (Perkins and Gonzales, 1985):

155
( For calculating the semi-axes of the flood region a1 and b1 within the ellipse, the following equations are used: 160 Where Vwt = total volume of the flooded region

F2 = an intermediate calculation parameter
And all other parameters are the same as those defined in Section 2.2.
Pore pressure distribution can be obtained by using equations 24 to 27 (Koning, 1985;Ge & Ghassemi, 2008;Perkins and 165 Gonzalez, 1985): Where Δ 1 = pore pressure increase at the elliptical boundary of the flood front , cw and cf = grain, water and formation compressibility, respectively = reservoir permeability = relative permeability to "water" (i.e., injected fracture fluid filatrate) at residual oil saturation 180 0 = viscosity of oil at reservoir temperature ℎ = viscosity of fracture fluid filtrate at reservoir temperature = viscosity of fracture fluid filtrate at cool temperature and = axial parameters for an elliptical coordinate system, as shown in Figure 6. 0 , 1 and 2 = the boundaries of cool region, flood region, and pressure front, respectively 185 For convenience of analysis and visualization, it is beneficial to use an elliptical coordinate system for this type of analysis. An elliptical coordinate system ( − ) is shown in Figure 6. Confocal ellipses and hyperbolae create two-dimensional orthogonal elliptical coordinates. Two foci are located at fixed positions of -Lf and Lf on the X-axis of the Cartesian coordinate system (X-Y). For any point around the fracture in an X-Y coordinate, the elliptical coordinate can be set as (Ge & Ghassemi, 190 2008

Analysis of Stress Rotation Due to Stress Shadows
The induced stress changes due to stress shadow effects will generally result in principal stresses that are rotated with respect to the in-situ principal stresses. For analyzing stress rotation around a fracture, principal stress magnitudes are calculated using equations 30 and 31. 200 https://doi.org/10.5194/se-2021-1 Preprint. Discussion started: 13 January 2021 c Author(s) 2021. CC BY 4.0 License.
Where 1 and 3 = principal stresses in the X-Y plane Principal planes are defined as the planes on which the principal stresses act and the shear stress is zero. Stress rotation angle is determined by using Equation 32.
Where = stress rotation angle around the vertical axis. 210

Evaluation of Shear Failure Potential on Natural Fractures
The shear failure potential of natural fractures around an induced hydraulic fracture was analyzed in this work, because slip on natural fractures can induce microseismic events which can be detected and located if appropriate monitoring technologies are used. Several failure criteria for shear failure for rock joints have been developed by various scholars including the  Bandis failure criterion (Barton & Choubey, 1977;Barton andBandis, 1980, Barton, 1976). Various investigators have offered criteria for shear failure of rock mass (Sheorey, 1997) such as the Bieniawski-Yudhbir criterion, the Ramamurthy criterion, and the Hoek-Brown criterion (Hoek & Brown, 1980). In petroleum geomechanics, the Mohr-Coulomb failure criterion is widely used due to its simplicity. The Mohr Coulomb failure criterion for natural fracture is given as follows:

= natural fracture friction angle
Failure along the natural fracture can be quantified in terms of a slip criterion which is defined as follows: 225 Where Slip = slip parameter (values ≥ 0 denote slip) ′ = effective normal stress | | = magnitude of shear stress, acting on natural fracture face A negative value of the parameters indicates that shear failure is not predicted for the fracture.
For analyzing this slip criterion, shear and normal stress on fracture plane are determined by using the following methodology 230 (Zoback, 2007): 1. Transforming stress from the X-Y-V coordinate to the geographic coordinate (N-E-D) (Figure 7). For this, the stress tensor around the hydraulic fracture in the X-Y-V coordinate system, denoted by S (Equation 35) is transformed to Sg, which denotes stress in the geographic coordinate system (equations 36 and 37) Because is vertical: 2. Transforming stress from the geographic coordinate (denoted by Sg) to an arbitrary coordinate system (denoted by 240 Sf), which is defined by the orientation of a natural fracture (Equation 38).
245 Str = strike (measure clockwise from north; see Figure 8 for definition of strike using the right-hand rule) dip = dip angle (measured downwards from the horizontal plane) Figure 8 shows the natural fracture in three dimensions. If shear failure occurs in the fracture plane the slip direction is denoted by the rake angle (i.e., the angle between the slip direction and a horizontal line contained within the fracture plane).
Finding the shear and normal stresses acting on the fracture plane as follows: 250 (hence potential slip direction) acts in the direction of the rake angle (see Figure 8 for definition), which is determined as follows: If (3,2) > 0 and (3,1) > 0 or (3,2) > 0 and (3,1) < 0; then rake = arctan ( If (3,2) < 0 and (3,1) > 0; then 255 rake = π -arctan ( If (3,2) < 0 and (3,1) < 0; then In this work, the values of , , and used as input for slip analysis were calculated based on in-situ stresses plus stress shadow components (as defined following Equation 29). As such, the effects of induced stresses on natural fractures were considered. subdivided informally by many authors (Lefever et al., 1991;Gorjian, 2019). The middle member serves as the oil reservoir; 275 multi-stage hydraulic fracturing is generally required to achieve economic production rates from this member. The authors previously completed an extensive study of the Bakken Formation in the Viewfield region of southeast 280 Saskatchewan, including geomechanical site characterization using data from 13 wells, numerical simulation of hydraulic fracture simulation using GOHFER, and comparison of the simulation results to field data. Interested readers are referred to Gorjian (2019) for details. In this paper, we focus on stress shadow analyses conducted on one well from the original study, referred to here as well S for simplicity (the actual well name is 191/14-15-007-07w2). Well S was chosen because of the availability of instantaneous shut-in pressures (ISIP's) from each stage of a multi-stage hydraulic fracture treatment, and its 285 proximity to another well (11/16-28-008-07w2) which had been used for rock mechanical properties analysis. For initial conditions (i.e., no stress shadow effects; key input parameters as summarized in Table 1), the numerical simulator predicted a maximum fracture height of 19.8 m, half-length of 81 m, and average propped width of 5.7 mm.  Figures 12 and 13 were generated by the author using Matlab code written to solve the equations presented in sections 2.1 and 2.4. Figure 11 shows the normal stresses perpendicular and parallel to the fracture at the end of treatment (immediately after shut-in), whereas Figure 12 shows the shear stress in the horizontal plane and the rotation angle of the principal stresses.
According to these results, if hydraulic fracturing of stage 3 began immediately after shut-in of stage 2, the minimum horizontal stress would be greater than the original in-situ value (22.7 MPa compared to 20.8 MPa), the maximum horizontal stress would 305 be less than the original value (25.9 MPa compared to 26.1 MPa), and the tips of the hydraulic fracture would be rotated 11° towards stage 2 (rather than parallel to the stage 1 fracture).
In general, it is most important to note that the stress shadow results in an increase in minimum horizontal stress, hence an increase in the pressure required to initiate and propagate a fracture, and in the tendency for the induced fracture to deviate towards the previous fracture stage. Furthermore, sensitivity analyses showed that these effects become more pronounced as 310 fracture stage spacing is decreased, and the fracture rotation tendency can become particularly pronounced if the anisotropy in horizontal stresses is small. One final note pertaining to the results shown here is that these results represent an upper bound because they were generated for a scenario with no lag time between stages 2 and 3. In reality, pressure dissipation will occur with time following stage 2, which would diminish the magnitude of the stress shadow. Eventually, the fracture width would become constant with time, once both faces are in firm contact with proppant. 315

Thermo and Poro-elastic Stress Shadows
Thermo and poro-elastic stress shadows were analyzed by writing code in Matlab to solve the equations presented in sections 2.2 and 2.3. It should be noted that all the key input parameters in this analysis is given in Table A.1. To the best of authors' knowledge, it is the first time that the thermal properties (Thermal conductivity, linear coefficient of thermal expansion, specific heat of mineral grains per unit volume) for Bakken Formation are measured in laboratory and reported. These shadows 320 were found to be insignificant for the Bakken Formation, largely due to its low permeability (hence the limited extent of pressure increase and cooling around fractures). To illustrate the use of these solutions, extreme cases were analyzed; the results are given in Appendix A (thermo-elastic) and Appendix B (poro-elastic). It is worth noting that it is not necessary to analyze the distribution of poro-and thermo-elastic stress shadows throughout the entire model domain, if seeking to minimize analytical effort while assessing the relative significance of the different stress shadow mechanisms. Rather, it would be 325 sufficient to calculate the change in minimum horizontal stress at a point corresponding to  = 90° and R = fracture stage spacing for each stress shadow mechanism (see Figure 4). Any stress change deemed negligible at this point should have a negligible influence on a fracture being initiated at this point. For example, in the case analyzed here, the horizontal stress change predicted at this point is 1.9 MPa for the mechanical stress shadow mechanism, 0 MPa for the thermo-elastic and 0 for the poro-elastic.

Stress Shadow Modeling for Multistage Hydraulic Fracturing
The lateral spacing and time lag between successive stages for well S are shown in figures 13 and 14, respectively.

Stress Shadow of Stage 2 on Stage 3
In Section 3.2.1, the stress shadow effect for stage 3 was analyzed for conditions existing immediately after shut-in for stage 2. However, in practice there is generally a time lag between performing the successive fracture stages, which should be 335 considered because during that time fluid can leak off from the fracture to the formation, thus decreasing pressure and the magnitude of the stress shadow effect. Based on a 73 minute time lag between stages 2 and 3 (see Figure 14), the calculated distributions of normal stresses are shown in Figure 15 and shear stresses and principal stress rotation in XY plane are shown in Figure 16.

Cumulative Stress Shadow of Effects of Stage 2 and 3 on Stage 4
By considering 219 minutes and 146 minutes lag times after stages 2 and 3 respectively, the predicted distribution of minimum horizontal stress and principal stress rotation angle at Y = 96 m (i.e., the location of stage 4) are shown in Figure 18. Figure  390 19 shows the minimum horizontal stress and maximum horizontal stress versus Y at X=0 (center of fracture). Notable observations based on Figures 18 and 19 are as follows:  The predicted principal stress rotation angle doubles from 6° to 12° when accounting for a cumulative stress shadow on stage 4 caused by stages 2 and 3.

Stress Shadow for the Final Fracture Stage ("Stage 30")
Calculations for stage 30 were undertaken in a manner consistent with those described for stage 4, except in this case the time lag and distance between stage 2 and 30 were used to calculate stage 2's stress shadow for stage 30. Then the time lag and distance between stages 3 and 30 were used to calculate stage 3's stress shadow for stage 30, and so on (up to stage 29). Upon 410 summing the stress changes resulting from stage 2 through 29, the stress conditions during stage 30 were obtained.  Figure 21. This figure shows a rotation angle of 35° at the tip of the fracture, which is greater than the 6° predicted for the effect of stage 2 on stage 3 ( Figure 16). Figure 22 shows horizontal stress magnitudes in the X 415 and the Y directions versus Y at X = 0 (center of fracture) immediately before stage 30. This figure also reveals that the vertical line at Y = 1294 m is the boundary of a stress reversal zone. Beyond this boundary (i.e., for Y > 1294 m), the horizontal stress state become nearly isotropic up to Y = 1341 m, beyond which normal stress in the X direction becomes greater than normal stress in the Y direction. This result is important because it suggests that fracture orientation could rotate by up to 90, or there might not be a strongly preferred growth direction at all. 420

Model Simplification
Given that the stress changes occurring along the well trajectory (at each fracture port) are the only components of stress change that can be used as input to Gohfer, calculation of the full stress field is not useful for fracture simulations. For such conditions, the following simplified equations can be used for stress shadow modeling. For X = 0 (for any point along the well trajectory), the stress shadow in the Y and X directions for stage n can be calculated by simplifying equations 1 and 2 as 435 follows: Given that horizontal stresses may be reversed within the stress reversal region, the minimum horizontal stress must be assessed as follows: The model may be further simplified by assuming constant values for fracture length (Li), average width ( ) and fracture height (hi) by assuming the propped values for these parameters (i.e., rather than modeling leakoff for each stage, taking the 445 final values for these parameters after pressures have dissipated and the fracture dimensions are controlled solely by the proppant contained within the fractures). Results obtained using this simplified approach seem promising, as illustrated in  properties that make them amenable to shear failure ("slip") during the injection of fracture fluid, dilation of these fractures may result in a zone of increased permeability (i.e., a stimulated zone). Furthermore, slip on natural fractures will generally result in the release of acoustic energy (microseismic events), which means there is a close association between the interpretation of microseismic monitoring data and natural fracture attributes (e.g., Zoback, 2007). 460 In the stress shadow analysis presented in Section 3, the effects of possible natural fractures were not considered. According to core logging conducted by the authors, natural fractures were observed in one well in the study area, possessing a strike of 155° and dip of 65° (Gorjian, 2019). shown has reached its full length. It is expected that the shear failure zone would have started near the well and grown progressively outwards as the hydraulic fracture propagated, hence resulting in a zone of slip (and associated microseismic events) extending from the well towards an outer bound delineated by the darkly shaded regions in Figure 24.

Assessment of Stress Shadow Effect Based on Field Data
To assess the validity of the stress shadow modeling workflow developed in this research (see Section 3.1.2), the instantaneous shut-in pressure (ISIP) values measured for stages 2 through 30 during field operations at well S were plotted and compared 480 against the model-predicted values of the minimum horizontal in-situ stress (Figure 25). Bottomhole ISIP may be regarded as an estimate of minimum in-situ stress measurement, especially in fracture treatments involving relatively small, planar fractures. The predicted trend of increasing minimum horizontal stress with increasing stage number (Figure 25

Microseismic Event Evaluation
To evaluate potential hydraulic fracturing-induced microseismic events in the Bakken Formation, shear failure analysis was conducted for an assumed (hypothetical) natural fractures surrounding a hydraulic fracture in the Bakken Formation for conditions existing at well L (191/05-16-009-10w2), which was the only well for which microseismic monitoring data were publically available. Mechanical stress shadow modeling (for a single fracture stage) was repeated for well L, in the same 500 manner as presented for well S in section 3.2.1. Similarly, slip analysis for well L was repeated, in the same manner as presented for well S in section 4. Given the similarity in conditions at wells S and L, and for the sake of brevity, the results are not presented here. Readers may refer to Gorjian (2019) for more details regarding well L. Table 2 shows a comparison of horizontal extents of shear failure zones with zones of microseismic events logged at well L (see Appendic C). The extents of the predicted zones of shear failure compare favourably with the field data for the scenario 505 with natural fractures striking at 115° with a dip of 65° and friction angles in the range of 25° to 30°. This dip is consistent with observations reported in the lone well that was logged and found to have natural fractures, however the strike is 40° less than the value of 155° that was measured. Given that the strike direction of 155° is based on a limited dataset observed at a different location, it is not unreasonable to expect that fractures striking at 115° might be present near well L. Clearly, more investigation into natural fracture attributes in the Bakken Formation and more public reporting of microseismic monitoring 510 results is required before more confident analyses can be undertaken. When additional data are acquired, the modeling workflow presented here will serve as an effective tool for predicting estimates of the shear failure zone (hence the stimulated reservoir volume).

Conclusion 520
The following is a list of conclusions based on this study:  A new toolbox, based on analytical solutions, was developed for predicting stress shadow effects in multi-stage hydraulic fracture, and was demonstrated to compare favourably against ISIP data collected from a well in the study area. The formulation used in this work only considered mechanical stress shadow effects, as thermo and poro-elastic stress shadows were found to be negligible in the Bakken Formation due to the small injection volumes and times 525 and low permeability. However, if the permeability was 2.5 md instead of 0.1 md, the poro-elastic stress shadow should be considered as well. Thermal stress shadow at the tip of fracture can increase the length of fracture; however, it can also decrease the length of fracture if it creates secondary thermo-elastic fractures perpendicular to the primary hydraulic fracture. Analysis showed that the thermo-elastic stress shadow was still too low to create secondary thermal fractures. The effect of thermo-elastic stress shadow at the fracture tip was not analyzed in this work. 530  Shear failure of natural fractures around a hydraulic fracture was analyzed and compared against microseismic monitoring data reported for one well in the study area. The results showed a favourable comparison for a scenario with natural fractures striking at 115° with a dip of 65° and friction angles in the 25° to 30° range. More investigation of natural fracture attributes in the Bakken Formation is required in order to implement this type of analysis in the future, in order to predict stimulated reservoir volumes. 535

Conflict of interest
The authors know of no conflicts of interest with this publication.