Numerical Simulations of Fracture Propagation in Jointed Shale Reservoirs under CO 2 Fracturing

Water-based hydraulic fracturing for the exploitation of shale gas reservoirs may be limited by two main factors: (1) water pollution and chemical pollution after the injection process and (2) permeability decrease due to clay mineral swelling upon contact with the injection water. Besides, shale rock nearly always contains fractures and fissures due to geological processes such as deposition and folding. Based on the above, a damage-based coupled model of rock deformation and gas flow is used to simulate the fracturing process in jointed shale wells with CO2 fracturing. We validate our model by comparing the simulation results with theoretical solutions. The research results show that the continuous main fractures are formed along the direction of the maximum principal stress, whilst hydraulic fractures tend to propagate along the preexisting joints due to the lower strength of the joints. The main failure type is tensile damage destruction among these specimens. The preexisting joints can aggravate the damage of the numerical specimens; the seepage areas of the layered jointed sample, vertical jointed sample, and orthogonal jointed sample are increased by 32.5%, 29.16%, and 35.05%, respectively, at time t = 39 s compared with the intact sample. The preexisting horizontal joints or vertical joints promote the propagation of hydraulic fractures in the horizontal direction or vertical direction but restrain the expansion of hydraulic fractures in the vertical or horizontal direction.


Introduction
In recent years, the application of hydraulic fracturing has significantly increased the production of oil and natural gas.It is estimated that hydraulic fracturing has improved recoverable reserves of oil by at least 30% and of natural gas by 90% [1].In addition, more than 60% of oil and gas wells need to be fractured first, especially for unconventional gas resource deposited in deep underground shale layers with extremely small permeability (usually less than 1 mD) [2,3].Hydraulic fracturing may increase permeability of shale by three to five orders of magnitude [4].Thus, the production of the fractured gas wells increases dramatically.In all, hydraulic fracturing introduced in oil and gas fields has revolutionized gas production around the world [5,6].
However, the water-based fracturing technology has some limitations and environmental concerns.First, hydraulic fracturing consumes large amounts of water.According to a report [7], the water use for hydraulic fracturing accounts for 9% of total freshwater consumption in Texas.The large consumption of water will restrict the oil and gas reservoir developments in a water-deficient area [8,9].Second, due to the large amount of water and chemical reagents used in the hydraulic fracturing, it may cause potential water pollution and chemical pollution if the treatment of flowback fluids with chemical reagents is insufficient [10][11][12].Third, for the reservoirs containing clay minerals, the permeability may decrease after water injection, thereby decreasing the production of gas reservoirs [13][14][15].The main reason is that, when the hydration minerals meet injection water, clay minerals can swell and result in blockage of seepage channels [16].All of these disadvantages of hydraulic fracturing promote the study and development of waterless fracturing [17].
Several waterless fracturing technologies have been introduced in oil and gas industries over the past few decades, including oil-based fracturing, N 2 fracturing, and CO 2 fracturing [18].Oil-based fracturing was first used in Colorado, Texas, and Kansas in late 1940s [19].Compared with hydraulic fracturing, it could be conducted in frozen areas.However, oil-based fracturing is expensive and may impair the effective permeability of wells [20].N 2 fracturing and CO 2 fracturing are the two most popular fracturing methods because they are more economical and efficient compared with hydraulic fracturing [21].According to engineering production data [22], the production of the reservoirs stimulated by CO 2 is 1.9 times as much as the production of those stimulated by N 2 .The laboratory experiment results indicate that the gas with lower viscosity, higher diffusivity, and lower surface tension can penetrate into smaller pore space to create more complex fracture networks compared with hydraulic fracturing [23,24].In addition, the fracture surface created by gas fracturing has a larger roughness and complexity, resulting in a greater increase in permeability [4].
In addition, shale reservoirs always contain kinds of joints caused by the geological deposition and folding [25][26][27].The existing joints have a significant influence on the initiation and propagation of the induced hydraulic fractures [28][29][30].The fracture networks of jointed reservoirs may be very complex due to the reopening of the existing joints, the expansion of hydraulic fractures, and the intersection between joints and hydraulic fractures [31].Nitrogen fracturing experiments were conducted on shale samples vertical and parallel to the bedding plane; the results indicated that a relative complex fracture surface is formed in the shale sample vertical to the bedding plane [32].He et al. [33] performed hydraulic fracturing on shale with bedding planes; the results showed that the bedding planes in shale formation have a significant influence on the propagation of hydraulic fractures.However, the mechanism of the fracture initiation and propagation in kinds of jointed reservoirs is not well investigated.It is important to learn the distribution of fracture networks for the successful design of stimulation in jointed reservoirs.
To this end, the numerical tools COMSOL and MATLAB are used to simulate the hydraulic fracture propagation driven by injection fluids in several jointed reservoirs.The distribution of fracture networks and the development of horizontal and vertical fracture radii are studied in this work.

Governing Equations
In the numerical simulation, CO 2 is injected into the borehole.Then, the rock mass begins to fracture with the increasing injection pressure.The process of CO 2 fracturing involves solid deformation and fluid seepage.In this part, a series of governing equations are set up for solid mechanic field and flow field.Besides, damage equations are introduced to describe the destruction of the calculation elements.

Rock Deformation and Damage Evolution Equations.
In this work, shale rock is assumed as an elastic continuum material, whose constitutive relation satisfies with the physical equation of elasticity.It should be noted that the influence of pore pressure on stress distribution is also considered in the equation.Thus, the modified physical equation can be induced as where σ ij is the total stress tensor, ε ij is the total strain tensor, G = E/2 1 + ν is the shear modulus of rock, E is the elastic modulus of the rock, ν is Poisson's ratio of the rock, ε v = ε 11 + ε 22 + ε 33 is the volumetric strain, δ ij is the Kronecker delta, α is the Biot coefficient, and p is the pore pressure.
The relationship between strain and displacement is expressed by a geometric equation as follows: where u i and u j are the components of displacement in i and j directions, respectively.Substituting the modified physical equation ( 1) and the geometric equation ( 2) into the equilibrium equation, then the modified Navier-type equation is induced as where f i is the component of the net body force.
Since the initiation and propagation of hydraulic fractures are studied in this work, a damage model is introduced to characterize the damage condition during the injection process.The damage model is used to determine whether shale damage occurs after every calculation step.For the calculation element, when the stress state meets the maximum tensile stress criterion or the Mohr-Coulomb criterion, the tensile crack or shear crack occurs.It should be noted that the tensile crack is first generated, because the compressive strength is ten times greater than the tensile strength.Equations (4) and ( 5) are the maximum tensile stress criterion and the Mohr-Coulomb criterion, respectively: where σ 1 and σ 3 are the first and third principle stresses; f t0 and f c0 are the tensile strength and compressive strength of rock, respectively; and φ is the internal friction angle.
When elements start to be damaged, the elastic modulus reduces correspondingly according to damage theory [34].The evolution of elastic modulus is defined as where E 0 is the initial elastic modulus of rock and D is the damage variable and is calculated as [35][36][37] where ε 1 and ε 3 are the first and third principal strains and ε t and ε c are the tensile strain and compressive strain, respectively.
2.2.Gas Flow Equation.Gas flow equations are defined to describe the injection gas flow in this part.The gas continuity equation during gas transportation is defined as where m is the gas mass per volume of rock, ρ g is the density of the injection gas, q g is the seepage velocity of the gas, Q m is the source origin, and t is the time variable.On the basis of Darcy's law, the seepage velocity of gas is shown as where k is the permeability of the rock and μ is the dynamic viscosity coefficient.
Assuming that shale rock is saturated by CO 2 after the injection, the gas content per volume of rock can be defined as m = ρ g ϕ, and ϕ is the porosity of the rock.Injected CO 2 gas enters the supercritical state when the pressure exceeds 7.56 MPa at the temperature 76.8 degrees Celsius [38].When the injection CO 2 is transformed from the gaseous state to the supercritical state, density and viscosity change dramatically under the different pressure.The evolution of density and viscosity of CO 2 varying with pressure is shown in Figure 1.The relationship between density, viscosity, and pressure can be described by interpolating function in the model calculation.Thus, the first item in equation 7 is induced to equation (10) [39] as The density of CO 2 changes rapidly with pressure whilst the porosity of shale changes slightly during hydraulic fracturing.Thus, ρ g ∂ϕ/∂t is ignored since ρ g ∂ϕ/∂t is much smaller than ϕ ∂ρ g /∂t , where c = 1/ρ g ∂ρ g /∂p is the compressibility coefficient of CO 2 and c can be calculated from Figure 1.
Putting equations ( 9) and ( 10) into equation ( 8), the gas continuity equation can be induced into equation (11) as Considering the influence of solid mechanics on the evolution of porosity, the dynamic evolution of porosity can be described by [40] where ϕ 0 is the initial porosity of shale rock, ε v0 is the initial volumetric strain, p 0 is the initial pore pressure, and K s is the bulk modulus of rock grains.

Characterization of Rock Heterogeneity.
Since rock is heterogeneous material with natural fissures and induced fissures [41], Weibull distribution function is introduced to represent the heterogeneity of shale rock in this part [42].
In this work, parameters, such as elastic modulus and strength, are assumed to conform to the Weibull distribution and produced in MATLAB.Weibull statistical distribution and the probability density function are defined as where x is the mechanical parameter of rock, x is the average value of x, and λ is the coefficient of heterogeneity. 3 Geofluids plane cross-section through a horizontal well, which is a 2D plane square with a length of 0.2 m and a borehole with a diameter of 0.04 m drilled in the center.The distribution of preexisting horizontal joints, vertical joints, and orthogonal joints is showed in Figure 2. The thickness of the joints is 1 mm among these models.For mechanical boundary settings, the load applied on the left boundary is equal to that on the upper boundary, the value is 5 MPa, whilst the right boundary and lower boundary are roller boundaries; these boundary settings are the same for the intact sample, layered jointed sample, vertical jointed sample, and orthogonal jointed sample.The gas is injected from the inner boundary with the injection rate of 0.0106 m 3 /s.As for seepage boundaries, there is no flow boundary on the outer boundaries.The parameters for the rock matrix and joints are listed in Table 1; it should be noted that the elastic modulus, tensile strength, and compressive strength of joints are only half of those of the rock matrix.

Numerical Implementation of the Model.
In the front parts, the coupling equations and model geometry are established.Due to the complexity of the coupling relationship between solid deformation and fluid flow, these equations are difficult to be directly calculated.Thus, the finite element method is adopted to solve coupling equations via COMSOL  Multiphysics and MATLAB.The distribution of stress and pore pressure is obtained to further discuss the condition of the numerical sample.The flow chart of numerical calculation procedures is showed in Figure 3.

Model Validation.
It is important to validate the effectiveness of the model before it is used to simulate the fracture propagation in the jointed reservoirs.The classical theoretical solution for forecasting the breakdown pressure [43] is shown as where p b is the breakdown pressure, σ t is the tensile strength of the rock, σ 1 and σ 3 are the maximum and minimum principal stresses, respectively, and p 0 is the initial pore pressure.
The numerical specimen used in this part is the same with the intact sample shown in Figure 2 (a).In this part, the initial pore pressure is 1 MPa, the average tensile strength and compressive strength of the specimen are 5.9 MPa and 59.4 MPa, respectively, and the average elastic modulus of the numerical specimen is 35.9GPa.The horizontal stress is fixed at 25 MPa, and the vertical stress is varying from 10 MPa to 25 MPa.A tectonic stress coefficient β σ x /σ y is defined in this part, and the breakdown pressure versus the tectonic stress coefficient is calculated in this part.
A comparison of the breakdown pressure with different tectonic stress coefficients by the theoretical solution and numerical solution is shown in Figure 4.The results show that the numerical simulation results agree well with results obtained from the theoretical equation, indicating that the proposed model is suitable to simulate the hydraulic fracture propagation under CO 2 fracturing.

Results and Discussion
4.1.Fracture Propagation in the Intact Sample.The hydraulic fracture initiation and propagation in an intact sample without joints are shown in Figure 5.To distinguish tensile cracks and shear cracks, the damage in tension is defined as negative whilst that in shear is defined as positive.It can be seen that the fractures first appear around the borehole and propagate gradually to the surrounding rock with the increase of injection gas pressure [44].Eventually, several main fractures are formed uniformly in the rock sample.As shown in Figure 5, hydraulic fractures are principally formed in the tensile mode in the intact sample during gas fracturing, and the number of damaged elements in the tensile mode accounts for 93.6% of the total number of damaged elements at time t = 37 s.
Figure 6 shows the development of the seepage area versus injection time.It should be noted that the cracks (damage area) are generated from the destruction of the elements.The  5 Geofluids seepage area keeps as 0 m 2 due to the low gas pressure from t = 0 s to t = 15 s.Then, the seepage area increases slowly from 3.27e -5 m 2 to 3.44e -3 m 2 from t = 16 s to t = 31 s with fracture initiation and gradual propagation.Besides, it increases dramatically from 3.44e -3 m 2 to 6.93e -3 m 2 from t = 31 s to t = 39 s, indicating that unstable fracture propagation occurs.The whole process of the evolution of the seepage area can be fitted by an exponential function as shown in Figure 6.
The development of horizontal and vertical fracture radii versus injection time in the intact sample is presented in Figure 7.It can be seen that hydraulic fractures initiate at time t = 16 s.As the injection gas pressure in the borehole further increases, the fracture radii in horizontal and vertical Step 18 Step 22 Step 29 Step 37 -0.8 0.8 -0.6 0.6 -0.4 0.4 -0.2 0.2 0   directions increase at a similar rate, reaching 6.3e -2 m and 6.49e -2 m at time t = 39 s, respectively.The slight difference of two curves shown in Figure 7 may result from the heterogeneity of mechanical parameters assigned to the sample.

4.2.
Fracture Propagation in the Layered Jointed Sample.
In the shale gas wells, the surrounding rock contains kinds of joints and fissures caused by the geological deposition and folding.The existence of these joints and fissures has a significant influence on the propagation of hydraulic fractures.Based on the above, hydraulic fracture propagation in kinds of joint samples will be discussed in the following parts.The fracture propagation in layered jointed samples is researched in this part.
The distribution and evolution of the fracture in the layered jointed sample are presented in Figure 8.At time t = 15 s, the sample begins to fracture.Hydraulic fractures first emerge around the drilling hole at time t = 18 s due to stress concentration in the inner boundary.With the increase of the injection time, hydraulic fractures propagate to the Step 18 Step 23 Step 30 Step 37  7 Geofluids surrounding rock.Finally, fractures propagate along the diagonal direction and horizontal direction forming several main fractures, as shown in Figure 8.The distribution of the fracture is different from that in the intact sample.This result from the fracture propagation favors the horizontal preexisting joints along which the least energy is dissipated due to their low strength and elastic modulus.

Horizontal
The seepage area begins to increase from time t = 14 s; then, it increases slowly from 7.92e -6 m 2 to 3.22e -3 m 2 from time t = 14 s to t = 29 s due to the low injection gas pressure at the initial stage.With the increase of injection gas pressure, the seepage area increases rapidly from 3.22e -3 m 2 to 9.18e -3 m 2 from time t = 29 s to t = 39 s.The seepage area increases by 32.5% in this case compared with that in the Step 18 Step 22 Step 29 Step 37 8 Geofluids intact sample at time t = 39 s.An exponential function is used to describe the development of the seepage area (Figure 9). Figure 10 is the development of horizontal and vertical fracture radii versus injection time in the layered jointed sample.Both the horizontal radius and vertical radius increase with the injection time.It can be seen that the horizontal radius is larger than the vertical radius during the whole fracturing process, due to the low strength of the joint in the horizontal direction.The horizontal radius and vertical radius are 8e -2 m and 5.55e -2 m at time t = 39 s, respectively.Compared with the intact sample (horizontal radius and vertical radius are 6.3e -2 m and 6.49e -2 m, respectively, at time t = 39 s), the preexisting joints promote the expansion of cracks in the horizontal direction but restrain the crack expansion in the vertical direction.

4.3.
Fracture Propagation in the Vertical Jointed Sample.The distribution and evolution of hydraulic fractures in the vertical jointed sample are presented in Figure 11.It can be seen that fractures are first generated around the borehole due to stress concentration.Then, fractures propagate along the diagonal direction and vertical direction with the increase of injection pressure, probably resulting from the equal load applied on the horizontal and vertical directions, and the preexisting joints in the vertical direction.Eventually, several main fractures are formed along the diagonal and vertical directions.Tensile damage is the main model of destruction, and the number of tensile cracks accounts for 93.9% of the total cracks.It can be seen from Figure 12 that the seepage area increases slowly from 3.96e -6 m 2 to 2.18e -3 m 2 with the injection time varying from t = 14 s to t = 26 s, owing to only a small amount of crack propagation at the low injection pressure stage.Then, the seepage area increases sharply from 2.18e -3 m 2 to 8.95e -3 m 2 with the injection time varying from t = 26 s to t = 39 s, due to the propagation of several main fractures at the same time in this period.Besides, the development of the seepage area is fitted by an exponential function, as shown in Figure 12.
The development of horizontal and vertical fracture radii versus injection time in the vertical jointed sample is presented in Figure 13.It can be seen from the curve that fractures begin to propagate at time t = 15 s.With the increase of injection gas pressure, both the horizontal and vertical radii increase gradually.It should be noted that the vertical radius is greater than the horizontal radius in the whole injection process.The horizontal and vertical radii are 5.1e -2 m and 8e -2 m, respectively, at time t = 39 s.Compared with the intact sample (horizontal radius and vertical radius are 6.3e -2 m and 6.49e -2 m, respectively, at time t = 39 s), the preexisting vertical joints promote the expansion of cracks in the vertical direction but restrain the crack expansion in the horizontal direction.

4.4.
Fracture Propagation in the Orthogonal Jointed Sample.Orthogonal joints are also common in naturally fractured reservoirs.The distribution and evolution of hydraulic fractures after gas fracturing in the orthogonal jointed sample Step 18 Step 23 Step 29 Step 9 Geofluids are shown in Figure 14.Fractures appear around the borehole at the initial stage (at time t = 18 s), and it propagates to surrounding rock along the horizontal and vertical directions at time t = 29 s.The preexisting orthogonal joints play a significant role in the propagation direction of the fracture: many generated fractures are extended in the direction of the preexisting natural fractures, as shown in Figure 14.During this period, fractures begin to expand gradually to form continuous cracks; the seepage area increases gradually from 1.98e -5 m 2 to 3.13e -3 m 2 .The main fractures are connected to form the crush area near the injection hole at time t = 37 s.As shown in Figure 15, the seepage area increases dramatically from 3.13e -3 m 2 to 9.36e -3 m 2 , and a complex fracture network is developed eventually.The seepage area increases 35.05% compared with that of the intact sample at time t = 39 s.Tensile damage is the primary damage mode in this situation with 93.02% of the damaged elements being in the tensile mode.
Figure 16 is the development of horizontal and vertical fracture radii with injection time in the orthogonal jointed sample.The fractures begin to propagate into the further surrounding rock at time t = 15 s.The horizontal radius and vertical radius are almost identical before time t = 26 s; then, the vertical radius becomes greater than the horizontal radius in the following process, which may result from the heterogeneous mechanical parameters used in this model.The horizontal radius and vertical radius are 7.3e -2 m and 8e -2 m at time t = 39 s, respectively.Compared with the intact sample (the horizontal radius and vertical radius are 6.3e -2 m and 6.49e -2 m, respectively, at time t = 39 s), the preexisting joints promote the expansion of fractures in the vertical direction.
Based on the numerical simulations conducted on the intact sample, layered jointed sample, vertical jointed sample, and orthogonal jointed sample, the numerical results can be concluded in Table 2.The numerical results indicate that the preexisting horizontal joints or vertical joints promote the propagation of hydraulic fractures in the horizontal direction or vertical direction but restrain the expansion of hydraulic fractures in the vertical or horizontal direction.And the preexisting orthogonal joints promote the propagation of hydraulic fractures both in the horizontal direction and vertical direction.The similar results can also be observed in the research of Wang et al. [28] as shown in Figure 17.In addition, the preexisting joints can aggravate the damage of the numerical specimens; the seepage areas of the layered jointed sample, vertical jointed sample, and orthogonal jointed sample are increased by 32.5%, 29.16%, and 35.05%, respectively, at time t = 39 s compared with the intact sample.

Conclusions
In this work, a series of numerical simulations are performed on jointed samples under the CO 2 -based hydraulic fracturing.A damage model is introduced to describe the initiation and propagation of microcrack in the fracturing process.The mechanical mechanism of the propagation of hydraulic fractures in several kinds of joint reservoirs is researched in this work.The numerical results of the intact sample, layered jointed sample, vertical jointed sample, and orthogonal jointed sample under CO 2 fracturing are shown in Table 2. Based on the results, the following conclusions are obtained.
It is shown that several continuous main fractures are formed uniformly in the intact sample.For the jointed  10 Geofluids samples, hydraulic fractures mainly propagate along the preexisting joints due to the lower strength and elastic modulus of joints compared with the rock matrix.Hydraulic fractures intersect with preexisting joints to form complex fracture networks.Besides, tensile damage is the main failure model in these numerical samples.The variations of the seepage area versus injection time among these conditions are similar.The seepage area increases smoothly at the beginning since hydraulic fractures propagate slowly under low injection pressure, and then, it increases rapidly with the propagation of a large amount of hydraulic fractures to the surrounding rock.The whole process can be described by an exponential function.The seepage areas of the layered jointed sample, vertical jointed sample, and orthogonal jointed sample are increased by 32.5%, 29.16%, and 35.05%, respectively, at time t = 39 s compared with the intact sample.For the complexity of the fracture networks, it can be described by the development of horizontal and vertical fracture radii.The preexisting horizontal joints or vertical joints promote the propagation of the fracture in the horizontal direction or vertical direction but restrain the fracture expansion in the vertical or horizontal direction.In addition, the hydraulic fracture propagation is promoted along the preexisting orthogonal joints to form the complex fracture networks.

Figure 2 :
Figure 2: Numerical samples with different kinds of joints: (a) intact sample, (b) sample with layered joints, (c) sample with vertical joints, and (d) sample with orthogonal joints.

Figure 3 :
Figure 3: Calculation procedures of the numerical model.

Figure 4 :
Figure 4: Comparison of the theoretical solution and numerical solution of breakdown pressure under different tectonic stress coefficients.

Figure 5 :
Figure 5: Distribution and evolution of the fracture in the intact sample.

Figure 6 :
Figure 6: Development of the seepage area versus injection time in the intact sample.

Figure 7 :
Figure 7: Development of horizontal and vertical fracture radii versus injection time in the intact sample.

Figure 8 :Figure 9 :Figure 10 :
Figure 8: Distribution and evolution of the fracture in the layered jointed sample.

Figure 11 :Figure 12 :Figure 13 :
Figure 11: Distribution and evolution of the fracture in the vertical jointed sample.

Figure 14 :
Figure 14: Distribution and evolution of the fracture in the orthogonal jointed sample.

Figure 16 : 2 )Figure 15 :
Figure 16: Development of horizontal and vertical fracture radii versus injection time in the orthogonal jointed sample.

Table 1 :
Mechanical parameters for the simulations.

Table 2 :
Numerical results of the intact sample, layered jointed sample, vertical jointed sample, and orthogonal jointed sample under CO 2 fracturing.