Wellbore stability model based on iterative coupling method in water alternating gas injection

Ensuring wellbore integrity is the most important factor in injection well design. The water alternating gas (WAG) injection is increasingly applied globally as the effective enhanced oil recovery (EOR) method in oil wells. High injection pressure or low injection temperature could lead to compressive wellbore failure. The rock stress around the wellbore is a function of the wellbore fluid flow and it should be precisely determined to avoid the wellbore failure. The purpose of this study is to propose a method to ensure the stability of the wellbore for the WAG process using iterative coupling method. The parameters of pressures, temperature, saturations and stresses are obtained for the multiphase flow condition using mathematical modeling. In this study, finite difference method is used to solve pressure, temperature and saturations; and finite volume method is acquired to solve the rock stresses. Iterative coupling method is employed to improve the accuracy of the results. This study introduces improved iterative coupling method between flow and stress models to reduce the processing time of obtaining corrected stress and failure results. Wellbore stability model is developed to determine the maximum pressure values, which lead to wellbore failure in WAG injection process for some different boundary conditions.

Abstract Ensuring wellbore integrity is the most important factor in injection well design. The water alternating gas (WAG) injection is increasingly applied globally as the effective enhanced oil recovery (EOR) method in oil wells. High injection pressure or low injection temperature could lead to compressive wellbore failure. The rock stress around the wellbore is a function of the wellbore fluid flow and it should be precisely determined to avoid the wellbore failure. The purpose of this study is to propose a method to ensure the stability of the wellbore for the WAG process using iterative coupling method. The parameters of pressures, temperature, saturations and stresses are obtained for the multiphase flow condition using mathematical modeling. In this study, finite difference method is used to solve pressure, temperature and saturations; and finite volume method is acquired to solve the rock stresses. Iterative coupling method is employed to improve the accuracy of the results. This study introduces improved iterative coupling method between flow and stress models to reduce the processing time of obtaining corrected stress and failure results. Wellbore stability model is developed to determine the maximum pressure values, which lead to wellbore failure in WAG injection process for some different boundary conditions.

List of symbols G
Acceleration of gravity (vector) (m/s 2 ) I Identity tensor (dimensionless) K Intrinsic permeability tensor (m 2 ) Krw Relative permeability for flow in phase w -l, g (dimensionless) qh Flux density for energy over all phases (J/ (m 2 s)) Sw Saturation of phase w -g, l (dimensionless) t Time (s) T Absolute temperature (K) u Velocity (m/s) z Elevation (m) a T Linear thermal expansion a Biot's constant for a porous media (dimensionless) b Turbulency factor e Total strain tensor (dimensionless) u, u f Porosity in general and porosity (dimensionless) lw Dynamic fluid viscosity of fluid phase w -l, g (Pa s) q l , q g Liquid and gas density (kg/m 3 ) r Macroscopic total stress tensor (tension positive) (MPa)

Introduction
Rocks are a combination of different materials. However, rocks exhibit poroelastic response. The amount of the stress indexed by pore pressure depends on pore content. The study of stress in a two-phase medium and in void space is essential for well integrity in oil production. The study of temperature is also important in defining the stress. The theory of thermo-poroelasticity or porothermoelasticity is developed by combining the influence of thermal stress and the difference between solid and fluid expansion (Espinoza 1983;Fredrich et al. 2000;Zare 2012). Enhanced oil recovery refers to several processes conducted to increase oil production from a reservoir after primary and secondary recoveries, which are typically performed by injecting water or gas. The injected fluid may push the oil or interact with the reservoir rock oil system to prepare suitable conditions for recovery. Injecting water alternating gas (WAG) shows better sweep efficiency than injecting water or gas alone (Irawan and Bataee 2014;Chalaturnyk and Li 2004;Chin et al. 2002).
Thermo-poroelasticity describes the effect of the change in temperature and fluid flow on the stress in the borehole and reservoir. Injecting water results in changes in temperature, pore pressure, and stress in the reservoirs and affects the reservoir permeability and porosity. Most reservoir simulators undergo stress changes and rock deformations during production because of the considerable physical effect of the geomechanical aspects of reservoir behavior (Lewis et al. 1986;Geertsma 1973;Hansen et al. 1995). Freeman et al. (2009) studied the geomechanics of bitumen formations. They had used two different simulators and compared the results (Freeman et al. 2009). Du and Wong (2005) had developed the coupled geomechanical thermal flow simulator. They had used the finite element method to express the reservoir model (whereas in most of the studies the flow is modeled by the finite difference method). The finite element formulation is explained well, but the result was not comprehensively expressed (Du and Wong 2005 (Yin et al. 2009).
The simulation study had used the iterative coupling method and studied the result of stability with the different initial values of parameters (Joseph et al. 2011). Safari and Ghassemi (2011) analyzed the geomechanical aspect of huff and puff process. They had done the study for a fractured geothermal reservoir. Their model had shown a good agreement with the field measurements. They had analyzed different geomechanical and flow behavior of the fractures after some years (Safari and Ghassemi 2011).
Chiaramonte and Zoback had published some books on the subject of reservoir geomechanics and CO 2 sequestration simulation. They had done another CO 2 -EOR simulation project in a fractured reservoir (2012). They had investigated the mobility of CO 2 in a fractured field (Chiaramonte 2012). Tran et al. (2004) developed new iterative coupling method and had applied it in CMG reservoir simulator. They had also corrected the porosity formula for the method. They call it pseudo-coupling and their study had shown that the result of this model is like the fully coupling method, but with higher speed. One year later (2005), they compare the different methods of coupling and their results (Tran et al. 2004). Mendes et al. (2012) had done a study of coupling with heterogeneity. They get their special boundary conditions and solve the two-phase flow problem using Monte Carlo algorithm. They reach to the result of locally conservative numerical solution and impress that there is an obvious change in production resulted by heterogeneity (Mendes et al. 2012).
Some research tried to model the sand production around the wellbore. Bianco and Halleck analyzed the mechanisms of arch instability and sand production in twophase saturated poorly consolidated sandstones (Bianco and Halleck 2001). Wan and Wang starts to model the sand production within a continuum mechanics framework (Wan and Wang 2000). He continued his work on sand production and published a paper 4 years later on the topic of ''Analysis of Sand Production in Unconsolidated Oil Sand, Using a Coupled Erosional-Stress-Deformation Model'' (Wan and Wang 2004). Some years later he analyzed the borehole failure modes and the pore pressure effects on it (Papamichos 2010). Papamichos developed the relation between sand production rate under multiphase flow and water breakthrough (Papamichos 2010).
A flow simulation is based on time. This process determines initial conditions and goes through a time in a defined time step (Irawan and Bataee 2014a, b). Geomechanical calculations are not based on time, except for such phenomena as creep that are usually ignored (Lewis et al. 1989;Rutqvist et al. 2010). However, the deformation and pore volume changes feed back to the time-based flow results. The degree of frequency of this updating procedure (implicitness) significantly affects running speed and result accuracy (Pao et al. 2001;Rutqvist 2011;Edalatkhah 2010). Such frequency can be categorized as follows: Full coupling: Flow and geomechanics variables (pressure, temperature, stresses, and strains) are solved simultaneously. Full coupling provides accurate solutions. However, this approach requires the solution of a large matrix and processing time is long.
Iterative coupling: Flow and geomechanics variables are solved separately and in sequence. This method has accuracy close to that of full implicit coupling but with higher speed.
Explicit coupling: Pressures, saturations, and temperatures data are called from the flow simulator to the stress calculations. However, the change in porosity, permeability and hence, the corrected pressures is not calculated. This method is called one-way or explicit coupling. Explicit coupling is fast, and lots of the geomechanics simulators use this method. However, the accuracy of this method is low because the flow characteristics depend on geomechanics and it is ignored in this method.
Pseudo-coupling: Some correlations between porosity and stress are used in flow calculations to identify compaction and dilation. However, this method does not process geomechanics (e.g., stress field), and simple formulas are used in a reservoir simulator to calculate subsidence during the process. The running speed of this method is high (Bataee and Irawan 2014;Bataee and Kamyab 2010;Settari et al. 2001).

Methodology
In most coupling studies, the parameters of pressure, temperature, saturation, stress, and strain are integrated. The full coupling can be performed by several methods, such as finite difference method (FDM), finite element method, and finite volume method (FVM). Thus, a large matrix can be solved. In this study, stress and strain solved using FVM, which is a suitable method for large meshes and able to deal with mesh concentration in high-stress sensitive parts. The other parameters (i.e., pressure, temperature, and saturation) calculated by FDM.
The relation between the change in porosity and strain change was considered. The new values of porosity and permeability are integrated into the flow equation to obtain the corrected values of pressure, temperature, and saturation. The stress and strain are updated. This procedure continues until the convergence condition is obtained under a certain level of accuracy.
The overall flowchart of study is expressed in Fig. 1. As in the flowchart, there are 4 main steps in the study that should be done to model the stability of the wellbore.
Step 1-Flow model: In this step, pressures and saturations around the wellbore are calculated. Temperature values are calculated after the results of pressures and saturations.
Step 2-Stress model: Stress and strain parameters are obtained in this part, based on the mechanical properties of wellbore rock. The stress equation is the main equation that should be discretized. The change in pressure and temperature is ignored in this part and its effect processed in next step to find the corrected stress values.
Step 3-Iterative coupling between two models: The fluid flow in porous media alters the stresses, porosity, and permeability. This effect results in the change in pressure distribution, which requires using two-way coupling method. The coupling study conducted after obtaining the required parameters. The nodes in the FDM and FVM are in different positions. Thus, the values are obtained using bi-linear interpolation of the nearby nodes. The relation between the change in porosity and permeability change is considered.
Step 4-Wellbore Stability model: Based on the results of stress values, which are obtained Stability model is based on the result of corrected stress values (from step 3). Proper failure criterion should be used to calculate the failure index around the wellbore.
Step 1: flow and energy modeling FDM employed to solve the parameters, namely, pressure, temperature, and saturation. The continuity equation for the water flow expanded with following conditions: cylindrical Fig. 1 The overall flowchart of this study. The flow part and stress part programed separately. After that, the iteration implemented, regarding the change in porosity and permeability J Petrol Explor Prod Technol (2016) 6:755-775 757 model, considering gravity, considering turbulence effect attributed to high injection rate, incompressible water, and slightly compressible oil and compressible gas with the use of implicit pressure-explicit saturation (IMPES) method while ignoring tangential flow. The turbulence effect is only considered in the wellbore with the Muscat equation. The condition for the pressure and temperature is chosen for immiscible WAG (iWAG). The water is injected, and the gas is injected subsequently under an immiscible condition. The constants for the equation related to rock properties were chosen for the sandstone reservoir. The energy balance equation for the three phases can be expanded after the saturation (from the flow part) is obtained to calculate the temperature values. The conditions for solving the temperature are as follows: cylindrical, implicit method and convection and conduction with the tangential flow ignored. The flowchart of the study in this part is expressed in Fig. 2. The programs are based on the equations obtained. First, the initial and boundary conditions for pressure and saturation should be set. Subsequently, the water injection (i.e., IMPES method) should be set. This process considers the loop of pressure in the loop of saturation for each time step. Gas injection is considered after a certain period of water injection. A thermal study was conducted after flow study and acquisition of all pressure and saturation data for any time of injection. The pressure in this study is expressed for different radius, degree, height, and time.
If we want to deal with the wellbore flow, the equations should be defined first and then discretized for the FDM. Thereafter, a program that can be used to determine the pressure and saturation for each node should be developed. The energy balance equation that uses the saturation obtained in the previous part should be discretized. After obtaining the temperature, the heterogeneity study can be conducted with the random matrix for porosity and permeability distribution through special random functions. The boundary condition may be changed to find the values in different cases.

Pressures and saturations
The steps in this part include water injection, gas injection, water-oil system, and three-phase flow. The continuity equation for the water flow expanded for the finite difference method with the following conditions: cylindrical model, considering gravity, considering turbulence effect attributed to high injection rate, incompressible water, and slightly compressible oil and compressible gas with the use of implicit pressure-explicit saturation (IMPES) method while ignoring tangential flow.
These are boundary conditions for flow cases: Water flow After considering and expanding the mass conservation equation and momentum conservation in the form of Darcy's law, the following equation is obtained: The flow in the tangential direction (the second term) is eliminated in this study. The derivation of z over z is equal to 1. The equation should then be prepared for the finite difference study: The total compressibility used in the case of oil flow equation is: Gas flow By considering and expanding the gas density and compressibility in the form of Darcy's law, the following equation obtained: The equation should then be prepared as follows for the finite difference study: Two-phase immiscible water-oil To solve the two-phase problems, four equations require obtaining the four unknowns, namely, P w , P o , S w , and S o . The four equations are water flow, oil flow, capillary, and total saturation equations.
Water flow system Total saturation Capillary equation For the water and oil relative permeabilities the Pirson's correlation is used: in which S w * is obtained from the following equation: The water and oil flow equation should be prepared for the finite difference study: and for the oil phase: Three-phase immiscible water-oil-gas Six equations require to obtain the six unknowns, namely, P w , P o , P g , S w , S o , and S g , to solve the three-phase problems. These six equations are water flow, oil flow, gas flow, wo-capillary, OG-capillary, and total saturation equations. WO-capillary equation OG-capillary equation Total saturation Equations 4, 5, and 6 are the same as before, but the relative permeability for the three phases adopted from Wyllie's correlations is as follows: The water, oil and gas flow equation should be prepared for the finite difference study as follows: For the oil phase: and for the gas phase:

Temperature
The energy balance equation for the three phases may be expanded after the saturation (from the flow part) is obtained. The following are the conditions for solving the temperature: cylindrical, implicit method, convection, and conduction while ignoring tangential flow. The following is the boundary condition for the study: This is the energy balance equation: After expanding for the finite difference method this relation is obtained: Step 2: stress modeling FVM employed to determine the stress and strain parameters. The stress equation is the main equation that should be discretized. The change in pressure and temperature are ignored in this part.
The flowchart of the study is shown in Fig. 3. These procedures should be followed to determine the stress and strain parameters. First, the equation should be transformed into a weak form and must be solved for each node. Subsequently, the body shape of the study has to be meshed. The program developed to identify the strain and, subsequently, the stress for any tetrahedral meshed body. The input for the program includes the initial and boundary conditions, as well as the position of the nodes and their connectivity.
FVM is employed to develop the stress model in the wellbore. The procedure consists of transforming the equations to weak form, meshing the defined shape, and programming to obtain the values for each node.

Weak form of equations
The equations should be transformed into weak form for FVM. The pressure and temperature distributions remain constant. The main equation in this part is stress Eq. (30), where ''f'' is the body force and is assumed zero in this case. div The equation of the poroelasticity defines the relation between stress and strain.
In this equation, f is the strain of the fluid part and may be written as follows: The volumetric strain e vol is as follows: where the strain in the x-, y-, and z-directions are expressed as in Eqs. (34)(35)(36).
where the matrix is transformed into Eq. (40), where u, v, and w are the deformations toward the x-, y-, and zdirections.
The procedure to change the equation to the weak form begins, and the equations for one direction (x-direction) are written. The procedure for the others is the same. The factor is taken from the equations, and the derivation of u over x is transformed to F (F = du/dx). The procedure used to obtain N and H is the same.
After the transformation, the matrix transforms into Eq. (42).
Therefore, this new equation is obtained.
The equations should be multiplied by the unit volume.
After rearranging these formulas, Eqs. (45-47) obtained. Therefore, they should be solved for every single node with the summation equation for all the nearby elements for F, N, and H (Eqs. 48-50).
Tetrahedral meshed shapes The body shape meshed to find the values by FVM. The program in this study developed for only tetrahedral shapes. The advantage of this method is that any shape can be applied. The requirements are only meshed positions, connectivity, and Fig. 3 The procedure to obtain the stress and strain using FVM. It consists of the three parts; weak form of equations, shape of model and the program boundary condition. Therefore, the mesh node positions can be imported from any software to the program. Figure 4 shows the example of the pyramids attached to a single node. All properties should be solved with the values of all attached nodes. In this study, the values of all attached nodes should fed into Eqs. (48-50). Figure 5 shows the projection of the triangle in the x, y, and z planes. The values for Dsx, Dsy, and Dsz should be known for the main Eqs. (45-47).
The program reads the data of the node positions, connectivity, and initial and boundary conditions. Subsequently, the program calculates all the element volumes and surface projections. The program then uses the iteration method to solve the matrix of Eqs. (45-50) for all nodes.
The shape of the model should be sketched (Fig. 6) to obtain the stress and strain results of the wellbore.
The shape should then be meshed (Fig. 7). This meshing system can be provided by any mesh generator software that can prepare the nodes and connectivity for the program (Fig. 8). Table 1 presents the input data of this study. The boundary conditions of the stress and flow conditions are expressed as follows: The simple explanation about the case in this study is that there is a wellbore shape as in Fig. 6 (meshed as in Fig. 7 for FVM study). The input parameters are as in Table 1 and the boundary conditions are as in Eqs. 51-55. For the process of the flow around the wellbore, firstly, water injected into fully saturated oil medium. Then, gas injected to this OW two-phase medium. The stresses applied are horizontal maximum, minimum, and vertical stresses. Horizontal maximum in situ stress is in the direction of north-south.
Step 3: iterative coupling In this study, stress and strain are solved using FVM, and the other parameters (i.e., pressure, temperature, and saturation) calculated using FDM. The coupling study is conducted after obtaining the required parameters. The nodes in the FDM and FVM are in different positions. Thus, the values are obtained using bi-linear interpolation of the nearby nodes. The relation between the change in porosity and strain change and that between porosity and permeability was considered. The new values of porosity and permeability are integrated into the flow equation to obtain the corrected values of pressure, temperature, and saturation. The stress and strain are updated. This procedure continues until convergency obtained under a certain level of accuracy (0.01 Psi).
The coupling is the result of the change in pressure, temperature, saturation, stress, and strain after returning from each model. First, the pressure, saturation, and temperature should be determined. The stress and strain can be calculated thereafter. The porosity and permeability change with the change in pressure and strain (Eq. 56).
Permeability changes with the porosity and correlation for the change expressed in Eq. 57.
The new values of pressure, saturation, and temperature are calculated with the updated porosity and permeability. This iterative coupling procedure continues until it converges. The new values of pressure, saturation, and temperature are calculated with the updated porosity and permeability. This iterative coupling procedure continued until it converged.
Step 4: Wellbore stability model The proper failure criterion for this study should have the capability of calculation the stress for three-dimensional stress direction (capable of considering intermediate stress), and it should be suitable for wellbore rock model. Drucker-Prager is the most suitable rock failure criterion  among all failure criteria regarding the wellbore modeling. Therefore, failure index values will be calculated using calculated rock stresses around the wellbore. The Drucker-Prager failure criterion is expressed, using principal stresses, as in Eqs. 58-62. Equation 62 expresses the failure index and if it becomes a minus value, the failure happens.
where s oct ¼ 1 3 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ðr 2 À r 3 Þ 2 þ ðr 3 À r 1 Þ 2 þ ðr 1 À r 2 Þ 2 q ð59Þ Results of the program In this study, some different in situ stresses are applied to check the result of the failure index values in the wellbore. All the cases are investigated for the boundary condition equations of 51-55 and input parameters of Table 1. To find the results of the stress values, pressures and temperatures should be calculated. Temperature values are related to the saturations of each phase, therefore, saturations also should be defined for each node during the time of WAG injection. Calculation of these parameters is crucial in the first step of injection as they have a sharp change in the wellbore.
At the first stage, water injection implemented. The pressure increases as the injection continues. Figure 9 shows the graph of oil pressure versus time at 3.65 cm distance from wellbore wall. As the injection starts, the pressure increases sharply from 3600 (original reservoir pressure). It increases incrementally after 10 s of injection.
The oil pressure profile is expressed in Fig. 10. It is obtained after 4 s of water injection. The pressure is boosted regarding the injection and the pressure enhancement continues to go further in the wellbore as water injection continues.
The water saturation profile is depicted in Fig. 11. These results are obtained after 4 s of injection. The results show how oil is pushed by the water around the wellbore. After 4 s, the injected water breaks through at 8 cm distance of wellbore wall. It continued further as water injection was continued. It is obvious in the graph that after 8 cm, the  area in untouched and the saturation of the water is 0.22 which is the connate water saturation. Figure 12 shows the water saturation at the node of 3.65 cm far from wellbore wall. It shows that how saturation increases with time. After the bank of water reaches, there is a sharp increase of the water saturation. After about 15 s, the water saturation slightly increases until the point that no more oil could be removed from the pores (residual oil saturation).
Gas injection conducted after water injection. Figure 13 shows the results recorded after 12 s of gas injection into two-phase medium of oil-water around the wellbore. As it is shown in the figure, gas phase reaches to 14 cm during this 12 s of gas injection. The gas saturation of the untouched area is zero and only water and oil phases exist in that area.
Temperature distribution obtained after the saturations of the 3 phases had been determined. Figure 14 shows the temperature distribution in nearby wellbore during injection. The initial value of the reservoir temperature is 70°C (refer to the Table 1) and the injection temperature is 50°C. The temperature distribution around the wellbore is highly sensitive to the injection.
The stress values and the failure are studied in three different cases. The three cases show different failure conditions, based on in situ stress values and injection pressures. The first case shows the normal failure (as in fracturing). The second case expresses the compressive failure and the last case describes the shear failure (as in sand production phenomenon).

Case 1
The in situ stress values and injection pressure is expressed in Table 2. Radial stress around the wellbore is depicted in Fig. 15. The horizontal maximum and minimum stresses are 35 and 30 MPa, respectively (refer to the Table 2). In the 2D graph (left side of Fig. 15), the effect of horizontal maximum stress is obvious in the direction of east-west. Radial stress values are very important in normal and compaction failure. Tangential   stress values are depicted in Fig. 16. The highest tangential stress value is recorded at the wellbore and in the direction of minimum horizontal stress. Figure 17 depicts the vertical stress around the wellbore. Vertical stress is the weight of the overlying strata. The original value of the vertical stress, in this case, is 36 MPa and the distribution shows that it is almost the same around the wellbore, except in the region of the wellbore wall. It is important to note that horizontal maximum and minimum stresses have no effect on the vertical stress in each node and only the distance from the wellbore is important. The shear stress around the wellbore is depicted in Fig. 18. The direction of the maximum and minimum shear stresses is 45°with respect to the maximum and minimum horizontal stresses. The failure index is obtained according to the stress results and the procedure of step 4. Figure 19 shows the failure index. In this case, the failure is called ''normal failure''. The result of this type of the failure is fracturing. As seen in the input data from Table 2, the amount of injection pressure is very high. This high injection pressure induced high radial stress in wellbore wall (Fig. 15). Therefore, fracturing will happen in this case and it is shown (in Fig. 19) as blue color. The fracturing should be in the direction of the maximum horizontal in situ stress, and it is approved in the graph result of Fig. 19. Moreover, the direction of minimum in situ stress is in the most stable condition.
To sum it up, fracturing will happen if the injection pressure becomes very high. In this case, the wellbore wall cannot withstand the induced normal stress. The direction of the fracturing is the same as the direction of the maximum horizontal in situ stress.

Case 2
The in situ stress values and injection pressure are expressed in Table 3. The horizontal maximum and minimum stresses are 39 and 30 MPa, respectively (refer to the Table 3). Figure 20 shows the failure index. In this case, the failure is called, ''compaction failure''. The result of this type of the failure is crushing the wellbore in the direction of minimum horizontal in situ stress (opposite direction from the case of normal failure). The wellbore wall crashes because the values of maximum horizontal stress are very big with regard to the minimum in situ stress value. In this case, if the wellbore pressure is not sufficient, the compaction failure will happen (as seen in the input data from Table 2, the amount of injection pressure is low). It is shown (in Fig. 20) as the blue color and it is in the opposite direction from normal failure. Such a case will happen in drilling operations, or in the production operation if the wellbore pressure drops.  Case 3 The in situ stress values and injection pressure is expressed in Table 4. The horizontal maximum and minimum stresses are 35 and 30 MPa, respectively (refer to the Table 4). Figure 21 shows the failure index. In this case, the failure is called, ''shear failure''. The result of this type of the failure is removing circular layers from the wellbore (as sand production). The wellbore started to produce sand because the wellbore pressure values drastically reduced regarding the oil production.

Sensitivity analysis
The sensitivity of wellbore fracturing to different parameters is investigated in this part. The parameters of interest in this part are injection pressure, temperature and in situ stresses. The results show the effect of each parameter on wellbore stability. Table 5 is considered as the base case that shows the injection pressure which leads to the initiation of the failure. In this case, 3165.4 psi injection pressure will initiate the fracturing of the formation.
The effect of the change in horizontal maximum stress on fracturing is expressed in Table 6. It describes that in the case of constant horizontal minimum stress, fracturing would be easier if horizontal maximum stress increases.
The effect of the change in horizontal minimum stress on fracturing is showed in Table 7. It expressed that in the case of constant horizontal maximum stress the wellbore would be harder to break as horizontal minimum stress increases.

Validation of failure results
The proper equipment that can serve the polyaxial test is rare and the test is very costly. Therefore, there are limited studies in this case. There are two experimental studies that can be used as the reference for this study, because of the injection condition and core characteristics. The result of this study is not completely matched the experimental results; some different facts cause this difference which is explained in each case.  Case 1 Athavale implemented an experimental study on rock failure using polyaxial testing. Figures 22 and 23 show the rock sample after failure. Figure 24 shows the result of his study (Athavale 2007;Kwasniewski et al. 2012). The vertical, horizontal maximum, minimum stress values, and injection pressure are expressed in Table 8. The result of this study is compared to the result of the model in the same table (Table 8). Figure 25 shows the numerically simulated regarding the experiment. It is shown that the direction of fracturing is in the direction of horizontal maximum stress while the other direction is stable.

Case 2
The hydraulic fracturing results can be used as the source of data for the validation if all the stress amounts and parameters are known. The hydraulic fracturing field data are collected from some different papers and the result is expressed and compared to the model in this part (Table 9) (Kwasniewski et al. 2012;Rahman and Rahman 2013;Raaen et al. 2006). Figure 26 compares the results of the field data and the model of the study. The model predicted the results with the accuracy of 91 %.

Case 3
Al-Ajmi (2006) presented the collection of the octahedral shear stress data for different rock samples. The data are from polyaxial tests on Shirahama sandstone rock. Table 10 shows the result from the experiment for in situ stresses and results of octahedral shear stress (Kwasniewski et al. 2012;Al-Ajmi andZimmerman 2005, 2006; Al-Ajmi

Conclusions
To ensure the wellbore stability, stresses values should be obtained. The stress values interact with pressures, temperatures and saturations regarding the change in porosity     To sum up the results of flow around the wellbore, it is recorded that the wellbore pressures and saturations changed very fast. It is due to the small wellbore area and high injection pressure. After the gas injection, OW bank pushed out of the wellbore, however, some amounts of oil and water remained in the pores and need chemical treatment to be removed. Temperature values affect the wellbore stress; in the case of injection, wellbore cooling will happen and might cause stability problems. Temperature values are related to the saturation distribution around the wellbore for each phase. Therefore, these values calculated after flow study had completed.
Stress redistribution will happen around the wellbore after the injection. The stress value is a function of in situ stresses, pressures, and temperatures. Maximum values of radial stress are in the direction of horizontal maximum in situ stress. The values are important in wellbore failure because fracturing will happen in this direction. Maximum values of tangential stress are in the direction of horizontal minimum in situ stress. Vertical stress around the wellbore is not related to horizontal maximum and minimum in situ stresses. It is a function of wellbore radius, pressures, and temperatures. The direction of the maximum shear stress is 45°with respect to maximum horizontal in situ stress.
Three different cases are investigated to show the three different failure types. In the first case, high injection pressure leads to normal failure as fracturing; it started in the direction of maximum in situ stress. In the second case, compaction failure occurred which is caused by the high difference between in situ maximum and minimum pressure and lack of well pressure support. It started in the direction of minimum in situ stress. The third case investigated the shear failure as in sand production. The low wellbore pressure caused the layers of the sand separated from the wall; this type of the failure is common in production wells.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Fig. 27
Comparison between experimental and model results