Injection Modeling and Shear Failure Predictions in Tight Gas Sands — A Coupled Geomechanical Simulation Approach

Abstract

This work presents theory for modeling of fracture propagation within reservoir simulator, history matching of field injection pressure using uncoupled and fully coupled geomechanical injection models, and sensitivity study of various parameters such as permeability enhancement/reduction functions, limiting length of fracture propagation, stress factor, and Biot’s constant. Two wells completed in tight gas sands in Western Canadian sedimentary basin were studied. The wells were fractured with different techniques (i.e., X-link gelled water fracs (Well A) and un-gelled slick water fracs (Well B)) and were both successfully matched with coupled geomechanical model.

1. Introduction

Fracture propagation modeling is an important part of reservoir geomechanics and must be considered in injection modeling of wells. Classical modeling of fracture geometry is well established and documented in literature of rock mechanics and stimulation [3]. Direct coupling of fracture propagation (fracture dynamics) and fluid flow is computationally very expensive [4, 5]. The modeling of “complex” fracturing [6] is also expensive. However, proper representation of dynamic propagation in which the fracture is directly coupled into a reservoir simulator is important for many applications. It has been shown that the some degree of coupled treatment of fracture mechanics, reservoir modeling and geomechanics is important for better understanding of the unconventional fracturing applications as well as for tight gas fracturing treatments such as waterfracs [7, 8]. While the fully coupled approach [5, 12] is not yet feasible for practical work, we have developed a simplified method which will be described briefly here. Such coupling, described next, can be achieved in a simplified fashion using only a reservoir simulator, which makes the modeling computationally efficient. The method has been developed and refined over two decades and proved to be successful in modeling a large variety of injection processes.

2. Theory of fracture propagation modelling

The method is based on modifications of transmissibility in the reservoir flow model. For injection scenarios transmissibilities are modified dynamically. To model dynamic fracturing process, a transmissibility multiplier function is assigned to a line (or plane) of grid blocks assumed for fracture propagation extending from the well. The multiplier function is a table that can be derived from simple 2-D analytical fracture models which approximate the actual fracture. In an uncoupled modeling, transmissibility multipliers are a function of fracture injection pressure, while in a coupled (geomechanical) system they are a function of a minimum effective stress. The multipliers are calculated based on the estimation of a 2D crack opening in a cross-section by Equation (1) [9], then calculating the fracture permeability by Equation (2) as a function of the net pressure in the fracture, and finally calculating the transmissibility multiplier (T_r) on the reservoir transmissibility of the block containing the fracture as described in Equation (3).

$$w_f=\mathrm{ }\frac{{∆P}_f{H}_f}{\stackrel{-}{E}}=\mathrm{ }\frac{4\mathrm{ }{∆P}_f{H}_f\left(1-{\upsilon }^2\right)}{E}=\mathrm{ }\frac{4\mathrm{ }\left(P_f-P_{foc}\right)H_f\left(1-{\upsilon }^2\right)}{E}\mathrm{ }$$E1

$$K_f=\mathrm{ }{R}_{fa}\frac{w_f^2}{12}\mathrm{ }$$E2

$$T_r=\mathrm{ }\frac{K_m{A}_m+\mathrm{ }{K}_f{A}_f}{K_m{A}_m}=1+\mathrm{ }{R}_{fa}\frac{w_f^3}{12\mathrm{ }{K}_{m}w}=f\left(P_f-P_{foc}\right)$$E3

H_f in Equation (1) is the estimate of fracture half-height based on the 2-D Perkins-Kern geometry assumption of vertical fracture with smooth closure at the top and bottom [10]. Permeability reduction factor R_fa in Equation (2) is a correction factor that accounts for deviations from Poiseuille law such as roughness, tortuosity and non-Darcy flow. It can be several orders of magnitude less than 1. Fracture opening or closing pressure (P_foc) is normally taken as the initial minimum horizontal total stress acting perpendicular to fracture face. Fracture opening or closing pressure may be actually greater (or lower) than the initial minimum horizontal total stress due to poroelastic and thermo elastic effects. The method allows one to create tables of dynamic transmissibility multipliers (as a function of pressure or stress) which are then used in a conventional reservoir simulator to propagate the high permeability (in the fracture plane) in time. As such, the method does not directly solve any fracture mechanics equations; although an estimate of fracture width can be obtained from Equation (3) from the known T_r at any grid cell in fracture plane and in time. Typically a single function is used, but variations of confining stress can be modeled by the use of multiple functions with different P_foc values. Finally, the created fracture volume can be accounted for by a similar function applied to porosity.

Obviously, the method is approximate in several respects. The fracture opening is computed from PK geometry for a fixed height while in reality the height varies along the distance from the well. The fracture opening also varies in the vertical direction. In spite of these approximations, the method provides a very realistic approximation to the results obtained by models based on fracture mechanics [4], and is capable of history matching complex injection sequences [13].

3. The field study

The wells studied are located in the tight gas sands in Western Canadian sedimentary basin. They were fractured with different techniques – Well A using cross-linked gelled water fracs, and Well B with un-gelled water fracs (slick water fracs). In both cases, the entire well was fractured through an open hole as opposed to multi-stage fracturing (which is the more common technique). However, microseismic monitoring and other techniques have shown that a number of fractures were created, with fairly regular intervals.

Field bottomhole injection pressure (BHIP) for wells A and B is given in Figure 1 and 2 respectively. Well B which is deeper than Well A has a higher bottomhole injection pressure than well A. Breakdown pressure, maximum pressure required to initiate fracture in formation, is also higher for well B. Complete set of reservoir, geomechanical and stimulation data has been provided in references [1, 2]. For modeling it was assumed that 15 fractures were created along the well with spacing of 50 m, as it was indicated by microseismic monitoring.

The simulations were carried out both in an uncoupled mode and coupled mode (solving both fluid flow and geomechanics). In both we employed the technique described above for fracture modeling, and also the pressure or stress dependent matrix permeability changes (which model the permeability enhancement in the SRV). The results of uncoupled modeling are not presented in detail, but in this case it is shown that coupled modeling is necessary to obtain history match (see [1, 2] for details).

Figure 1.

BHIP, Wellhead pressure and Proppant concentration - Well A

Figure 2.

BHIP, Wellhead pressure and Proppant concentration - Well B

4. History matching of field injection pressure — Uncoupled geomechanical injection models

The concept used for approximation of geomechanical effects in an uncoupled model and equations developed for production modeling [1, 2] can be used also for injection modeling. For injection cases change in pressure is always positive and consequently the effective mean stress is always larger when poroelastic effects are considered. For uncoupled modeling it can be assumed that the stress changes have stabilized and fracture opening or closing pressure is equal to the minimum (adjusted) horizontal total stress during the treatment.

In uncoupled history matching, any changes in stresses due to poroelastic effects must be incorporated manually in reservoir simulator for permeability multipliers and transmissibility calculation. Modified stresses are used in uncoupled model assuming that the hydraulic fracture increases in-situ stresses near wellbore and around fractures due to poroelasticity and permeability enhancement in that region must be predicted by using these modified stresses. Fracture height of 50 ft, Poisson’s ratio of 0.125 and Elastic modulus of 7.99 E6 psia is used for fracture transmissibility calculation. The method of SRV permeability multiplier calculation is explained in references [1, 2]; the strength of the nonlinearity is given by “stress factor” S. Value of S=6.0 was used for all simulation runs. Injection model was setup in reservoir simulator and run for given injection period by using the actual injection rate, subject to maximum bottomhole injection pressure of 10,000 psi.

4.1. Effects of fracture permeability reduction factor (R_fa)

Fracture permeability given by Equation (2) becomes an intrinsic permeability of smooth open fracture if the reduction factor is taken as 1. Fracture permeability in reality is much smaller than this value due to tortuosity, asperities interlocking, rock chipping at fracture face, unequal and rough surface of rock faces, and fracture degradation. Three injection cases were run to match field injection pressure with simulation injection pressure by varying only the fracture permeability reduction factor.

4.2. Effects of limiting length of fracture propagation

The initial runs produced a flat injection pressure while the field pressure is steadily increasing. A mechanism that would create larger pressure increase with time is required. One method is to restrict or confine fracture propagation in length (half length), which can be achieved by modifying transmissibility of grid blocks in fracture plane only within an assumed fracture half length. Possible justification is the scale-dependence of effective fracture toughness that was proposed theoretically and indicated by matching data [14]. Several simulation cases were run both for each well using different values of pre-determined maximum fracture half length. Reservoir parameters used in simulation runs are same as in base case.

Results of all uncoupled cases are not shown here except best matched case (see Figure 3 and 4) as history matching was not achieved. Only their effects are discussed here (detailed description of the results is provided in reference [2]). Results show that decreasing the R_fa factor pushes the injection pressure upward. Smaller reduction factor means smaller transmissbility mutlipliers and hence larger pressure drop down in the fracture. However, fracture propagation confinement did not improve the rising trend of injection pressure in uncoupled simulation.

5. History matching of field injection pressure — Coupled geomechanical injection models

In coupled geomechanical simulation, because stresses are continuously computed, P_foc in Equation (3) is updated at each time step and grid block in reservoir simulator by taking as an input effective stress from geomechanical part of the simulator. Therefore there is no need to modify stress data to correct for poroelastic effects. To run a fully coupled geomechanical simulation the original in-situ stress is used to calculate transmissibility and permeability multipliers, which are a function of effective stresses. Run times for coupled simulation are generally very large and consequently a detailed study for each parameter was not possible due to time constraints. The sensitivity study and calculations shown here are performed for well A. Only conclusions and end results are then applied to well B to get a history match.

Note that for fracture initiation (and propagation) minimum effective stress must be negative; in other words injection pressure should be higher than minimum horizontal total stress. Biot’s constant of 1.0 was initially used for effective stress calculation, but it was found that fracture initiation could not be achieved because the poroelastic stress component caused by injection pressure was too high and the total stress increased above the injection pressure limit (set at 10,000 psia). The smaller the Biot’s constant the slower is the increase in total stress and it is less difficult to fracture the rock. It was therefore concluded that Biot’s constant should be significantly less than 1.0.

5.1. Effects of limiting length of fracture propagation

Few simulation cases are run using different values of permeability reduction factor and confining length of fracture propagation. For this purpose, it was assumed that rock behaves as a perfectly elastic material which does not exhibit hysteresis during loading and unloading. A base case here (Case-1) was therefore set up allowing unlimited fracture propagation in y-direction and modifying Biot’s constant value in the geomechanical simulator to 0.65 (initial guess). Summary of history matching parameters are presented in Table 1. Simulation results for all these cases for well A and history matched case for Well B are presented in Figure 3 and 4 respectively.

Property – Well A	Case - 1	Case - 2	Case - 3	Case – 4
Biot’s constant (α)	0.65			0.75
Permeability reduction factor (Rfa )	0.00001		0.0000052	0.0000052
Fracture half length, ft	Not restricted	100	130	130

Table 1.

Parameters varied in coupled injection Cases 1 – 4 – Well A

Figure 3.

Comparison of simulation results and field BHIP - Well A

Figure 4.

Comparison of simulation results and field BHIP - Well B

The effect of fracture permeability reduction was discussed in detail in uncoupled simulation section; decreasing its value shifts pressure injection curve upward which can be observed in Figure 3. Although simulation results of Case – 4 of well A do not exactly match field injection pressure, it represents a reasonable history match. It is concluded that injection history match requires some mechanism to constrain fracture propagation at a late stage. This issue was not pursued further; however, the coupled cases show much improvement compared to the uncoupled simulations as shown in Figures 3 and 4.

5.1.1. Discussion on the late time history matching (Well A)

It is important to point out that history matching of field injection pressure after 190 minutes of injection for well A cannot be achieved through our simulation results (See Figure 3). It was observed in field treatment report (See Figure 1) that the injected proppant concentration was increased after 190 minutes to approximately three times of the overall average concentration. Our simulation study does not include coupling of fracture propagation simulation with proppant transport, modeling of fracture propagation based on downhole variable proppant concentration is not possible here and beyond the scope of this study. Fracture modeling in this work was performed based on total downhole amount of slurry injected. Late time history matching for well B is more acceptable.

5.2. Effects of stress factor ($\mathit{S}$)

Stress factor (S) defines shape of pressure/ effective stress dependent permeability curves and controls the permeability dependence on effective stress [1, 2]. The larger the value of S, the higher is the permeability dependence on stress. Permeability multipliers are applied in the whole reservoir except in fractured blocks. Increasing S value from 6 to 16 results in increase of permeability multipliers to several orders of magnitude but there is little difference in injection pressure except at late time when the pressures are lower with higher S value (see Figure5). We observed similar difference in results for injection pressures when the exercise was repeated for well B. It is therefore concluded that effect of leak off on injection pressures in low permeability formations is not considerable, although it affects fracture length and the match with microseismic (MS) data.

Figure 5.

Effects of Stress Factor (S) on BHIP - Well A

6. Failure predictions — Tensile and shear failure

The simulations of the injection process presented in this work showed that one must assume a substantial stress-dependent enhancement of permeability around the primary single plane fracture (SPF) to history match the injection pressures. Often it is postulated that the creation of this SRV is due to shear fracturing, i.e., creating shear failure. Coupled modeling provides us with the tool to investigate under what conditions shear fracturing occurs and what would be the extent of the SRV if it was caused purely by shear failure. This aspect is examined in the present section.

When tensile stress across a plane exceeds critical limit then tensile failure occurs. This critical limit is called the tensile strength or ultimate tensile strength (UTS). The tensile failure criterion is applied to determine the propagation of the main fracture (SPF) through grid blocks. In some rare instances, tensile failure can also occur in the reservoir around the SPF (e.g. due to thermal effects [11]).

When shear stresses along a plane in a specimen exceed shear strength of material, shear failure occurs. The shear strength of material / rock indirectly depends on the normal stress acting on the failure plane. There are different shear failure criterions available in literature such as Tresca, Mohr-Coulomb and Griffith. For this study Mohr-Coulomb criterion is used to predict failure mechanism during injection.

To investigate if tensile or shear failure will occur; time-history of pressure and stresses was extracted for specific grid blocks from a coupled simulation run. By plotting the Mohr circles in MATLAB^® we can make failure prediction of these grid blocks in graphical form. For this purpose, the history matched case, i.e., Case – 4 of well A was used. All the fractures behave the same way and pressure propagation is also approximately the same for all fractures. Therefore only one fracture is selected for this analysis which is fracture # 4 (4^th fracture from the line of symmetry) and conclusions drawn from this analysis will apply to all sets of fractures. Two grid blocks were selected and marked as shown in Figure 8, which represents cross section of the model in y-z plane of the fracture. The well is completed in x-direction and block 1 represents the perforation location.

6.1. Base case — High cohesion

The Mohr-Coulomb failure envelope was based on rock geomechanical data given in references [1, 2]. The base case used friction angle of 30⁰ and Uniaxial compressive strength of 321 Mpa (intact rock). Mohr – Coulomb circle progression is presented in Figure 6 for well A, where Circles 1 - 5 are for block 1 and circles 6 - 10 for block 2. Similar envelope can be constructed from coupled simulations output for well B. In Figure 6 there is no shear failure during injection because Mohr’s circles are much below the failure line. It is obvious that the dominant failure mechanism in these blocks is tensile because the SPF penetrated them. We also repeated the same exercise for all blocks/time steps and confirmed that no shear failure occurred.

Figure 6.

Mohr – Coulomb failure envelope for C_o= 46570psi – Well A

6.2. Case with low cohesion

More realistic case was run by reducing the uniaxial compressive strength by 10 times to a base value while keeping other parameters such as friction angle, elastic modulus and Poisson’s ratio the same as in the previous case. Mohr-Coulomb circles for this case are shown in Figure 7.

Figure 7.

Mohr – Coulomb failure envelope for C_o= 4657psi – Well A

Complete spatial map of the failure can be obtained by plotting “stress level” S_L, a feature offered in GeoSim® (Geomechanical Simulator) which represents the ratio of the size of the Mohr circle at any point to the circle at critical state, i.e., when the circle touches the failure line. Stress level therefore ranges between 0 and 1. When S_L < 1 there is no shear failure, and when shear failure is reached, S_L remains theoretically at 1. Stress level for fracture # 4 after 237 minutes of injection (end of injection) is shown in Figure 8 for the y-z cross section through the fracture plane. It should be noted that our simulations were not carried out using elasto-plastic modeling, but only linear elastic treatment, and therefore S_L can exceed 1. The modeling is not rigorous past shear failure, but it still provides useful picture of the possible extent of failure. In this case, failure is also predicted for planes adjacent to fracture plane.

From above simulation runs (which required increase of matrix permeability during injection), and knowledge of presence of micro cracks and heterogeneities in tight sands [2] as recorded by microseismic, it is concluded that the high original value of uniaxial compressive strength (which does not allow any shear events) is unlikely. Reducing the C₀ to account for weak planes and natural fractures then will predict possibility of shear fracturing and shear-generated SRV creation. However, the SRV based on shear failure is still very narrow and therefore one has to conclude that the majority of the matrix permeability enhancement should be contributed to matrix and micro-fractures. These results are preliminary and further work should be done using finer gridding and elasto-plastic modeling.

Figure 8.

Stress level after 237 mins of injection – Co= 4657 psi – YZ cross section – Well A

7. Conclusions

• The method used for modeling the fracture propagation is practical, and provides realistic representation of fracturing in reservoir models or coupled geomechanical models.

• Uncoupled modeling is not capable of history matching the injection pressures for the two wells studied.

• Coupled modeling achieves reasonable history match of both wells. The main factors that have been identified as important are the fracture permeability factor (R_fa) (which primarily shifts the pressure curve), the reservoir permeability dependence on stress and confining the length of fracture propagation (which causes to increase of pressure in later part of the job and thus improves the matches).

• Value of Biot’s constant controls the increase of effective stresses during pumping. For larger values of Biot’s constant it is very difficult to fracture the formation.

• Preliminary work on the modeling of shear failure region (SRV) shows that no shear events are detected when a high value of uniaxial compressive strength (UCS) of the rock is assumed, representative of intact rock. A narrow shear region is predicted when the UCS is lowered to represent media with pre-existing fractures or planes of weakness.

This work demonstrates the need for coupled geomechanical modeling in injection to capture poroelastic effects and stress alterations during stimulation.

Nomenclature

$A_f$= Fracture cross sectional area, ft²

$A_m$= Matrix block cross sectional area, ft²

BHIP = Bottomhole injection pressure, psi

$C_o=$Uniaxial compressive strength (UCS), psi

$E$= Elastic modulus, psi

$H_f$= Fracture half height, ft

$K_f$= Fracture permeability, mD

$K_m$= Matrix block permeability, mD

$L_f$= Fracture half length, ft

MS = Microseismic

$P_f$= Fluid (fluid) pressure, psi

$R_{fa}$= Permeability enhancement/reduction factor

$P_{foc}$= Fracture opening or closing pressure, psi

$S$= Stress factor

SPF = Single planer fracture

SRV = Stimulated reservoir volume, ft³

S_L = Stress level

$T_r$= Transmissibility multiplier

UTS = Ultimate tensile strength, psi

$W$= Grid block size in x-direction, ft

$W_f$= Fracture width, ft

$\alpha$= Biot’s constant

$\upsilon$= Poisson’s ratio

Acknowledgments

We would like to acknowledge help from Taurus Reservoir Solutions Ltd. for providing the TRS® reservoir and GeoSim® geomechanical simulator. We also wish to acknowledge the financial aid from the JIP consortium for Tight Gas Sands and Shale Gas Modeling at University of Calgary for supporting the research fund. We wish to thank Vikram Sen for helping during this project and Apache Canada for providing us data for this study.