Abstract

After a shale gas reservoir is fractured, hydraulic fractures interweave with natural fractures to form a fracture network. Numerical simulation based on the continuous fracture network model is a relatively economical and convenient method to predict fracture network morphology and size in the field application. However, some important factors, such as fracture height variation and filtration loss, have not been considered in the past continuous fracture network models. Therefore, this paper is aimed at establishing a novel continuous fracture network model to improve simulation accuracy. Firstly, this paper established a method to judge whether natural fractures develop or not. Then, a novel continuous fracture network model considering fracture height variation and asymmetry, filtration loss, fluid flow, and other key factors was established, and the forward algorithm and inverse algorithm of the model were proposed. At last, this model was applied in a field case to verify accuracy, and the average accuracy is more than 90%. Compared with the traditional Meyer software, the average error of prediction was reduced by 7.86%.

1. Introduction

Hydraulic fracturing has been used commonly in the development of shale gas reservoirs to initiate several dominant hydraulic fractures from the wellbore [14]. The technique is aimed at opening the widely existing natural fractures in the shale formation and generating the largest possible hydraulic fracture network for transporting hydrocarbon to the wellbore [59]. Thus, accurate knowledge of the fracture initiation and propagation of hydraulic fractures is significant for optimizing the fracturing treatment design and, ultimately, improving well productivity [1012].

Many key parameters affect the complexity of a fracture network, including the existence of natural fractures, anisotropy of stress, and heterogeneity of rock properties [13, 14]. Moreover, simulating hydraulic fracture propagation in a formation with preexisting natural fractures is very complex and requires proper consideration of key physical elements, including fracture propagation [15], fluid flow in the fracture networks [16], fracture height growth [17], interaction between hydraulic and natural fractures [18], and proppant transport in the fracture networks [19]. Microseismic and other geophysical methods are the main techniques for monitoring a fracture network in a shale reservoir [20]. However, the uncertainty of microseismic events and the high cost of fracture monitoring make the implementation of these techniques difficult for multiple wells [21].

Numerical simulation is a relatively economical and convenient method to predict fracture morphology and size [22]. Different modeling approaches have been developed recently to simulate fracture networks, mainly divided into the discrete fracture network model [23, 24] and the continuous fracture network model. The main difference between these two models is that the discrete model fully considers the influence of natural fractures, while the continuous model uses the orthogonal fracture model to approximately describe the natural fracture network. In theory, the discrete model is closer to reality. However, in field applications, natural fracture diagnosis is limited by technology and high cost. Therefore, both methods have their own research value. This paper is aimed at developing and applying a new continuous fracture network model.

The important breakthroughs in the development of fracture network models are summarized below: In 2008, Olson presented a continuous fracture network model capable of predicting hydraulic fracture propagation and interaction with preexisting natural fractures [25]. However, the model is based on fracture mechanics alone and does not include fluid flow or proppant transport. From 2009 to 2010, Xu et al. proposed a wire-mesh model to estimate fracture network dimensions and proppant placement in a network [26, 27]. This model has the advantage of fast computation speed due to being a semianalytical solution. But its fracture network pattern (i.e., fracture spacing) is relatively simple in simulating the geometric shape of fractures. In addition, the wire-mesh model neglects changes in fracture height and filtration of fracturing fluid. After 2011, Meyer established a relatively complete continuous fracture network model and developed commercial software [28], but its core problem was that the algorithm of model was not fully published, and this model did not take into account asymmetry of fracture height. In recent years, the discrete fracture network (DFN) model proposed by Meyer has been improved and extended by many scholars [29, 30]. For example, in 2017, Nejadi et al. proposed history matching and uncertainty quantification of DFN models in fractured reservoirs [31]. In 2020, Yao et al. discussed the role of natural fracture characteristics on the performance of hydraulic fracturing for deep energy extraction using DFN [32]. In 2021, Hyman and Dentz discussed the transport upscaling under flow heterogeneity and matrix-diffusion in three-dimensional DFN [33]. The development of the fracture network model is summarized in Table 1.

From the above literature review, a number of remaining challenges need to be resolved. To properly simulate the propagation of continuous fractures developed during hydraulic fracturing and get more accurate simulation results, a new model is needed. Therefore, this work researches the mechanical mechanism associated with fracture network formation based on fracture mechanics and proposes a kind of extended model of a continuous fracture network. This study establishes a mathematical model of the fracture network, including equations of width, fracture length, and fracture height. The model solves a system of equations governing fracture deformation, height growth, fluid flow, and proppant transport so as to predict the size of a complicated fracture network more accurately.

The paper is organized as follows. Firstly, it provides mechanical condition requirements related to fracture network formation to determine whether natural fractures have developed or not. Secondly, it establishes a continuous fracture network extension model according to the propagation law of continuous fracture networking and considers the change in fracturing fluid filtration and fracture height. Finally, it uses field experimental data to verify the accuracy of the model.

2. Mechanical Mechanism of Fracture Network Formation Based on Fracture Mechanics

Warpinski et al. [28, 34, 35] pointed out that fracture networks are the result of interaction between hydraulic fractures and natural fractures; researchers have held that a continuous fracture network consists of one main fracture and several branch fractures, as shown in Figure 1. The key to fracture network formation lies in the fact that branch fractures develop around the main fracture [36, 37]. Some natural fractures develop, while some do not. This work studied the mechanical mechanism of reservoir formation with and without natural fracture development.

2.1. Mechanical Condition Requirements for Reservoir Formation with Natural Fracture Development

For reservoirs with natural fracture development, when the net pressure in the main fracture of a certain length increases after its formation, a natural fracture or weak plane opens, and then, a fracture network is formed. At present, the linear criterion put forward by Warpinski and Teufel [38, 39] is the most widely used fracture propagation criterion. The mechanical condition requirements for the formation of branch fractures in a fracture network can be analyzed based on the fracture propagation of a naturally fractured reservoir. Figure 2 shows a diagram of a reservoir fracture network with natural fracture development.

The closed natural fracture (or weak plane) is acted on by the maximum horizontal principal stress () and minimum horizontal principal stress () of the far-field, and its angle from the maximum horizontal principal stress is .

According to the 2D linear elastic theory, shear stress and normal stress can be expressed as follows [40]:

According to the Warpinski criterion widely applied to fracture propagation of naturally fractured reservoirs, when the fluid pressure of a natural fracture is larger than its normal stress, fracture extension will occur [41] as follows:

When two fractures intersect, the hydraulic fracture end is connected with the natural fracture, and fracturing fluid floods into the natural fracture. The pore pressure of the natural fracture’s near wall is as follows:

By substituting Equations (2) and (5) into Equation (3), the net pressure required for the tension fracture of the natural fracture is as follows:

When , the net pressure reaches the maximum:

According to the Warpinski criterion, when the pressure in the natural fracture is low, normal stress acting on the fracture surface is pressure stress, and the fracture is closed. When the shear stress acting on the natural fracture is high, the shear slip of the natural fracture will occur easily [42]:

Effective shear stress () acting on fracture surface is as follows:

The fracture criterion of a compression shear fracture is as follows:

By substituting Equations (1), (2), and (11) into Equation (9), we get critical pressure (pc) for the occurrence of shear slip and branch fracture formation:

According to Equation (12), when , the net pressure is the largest, and is as follows:

The natural fracture surface is unbounded (). Net pressure required for shear slip is horizontal principal stress difference ().

According to Equations (7) and (13), maximum net pressure is positively correlated with . Therefore, the reduction in initial horizontal stress difference can help expand the fracture network fully.

The calculation formula for the difference between two horizontal principal stresses of a reservoir can be derived based on the maximum horizontal principal stress formula [(36)]:

Therefore, according to , combined with the net pressure in the fracture propagation model, a judgment can be made on whether a fracture network can be formed.

2.2. Mechanical Condition Requirements for Reservoir Formation without Natural Fracture Development

For a reservoir formation without natural fracture development, branch fractures in the rock body are required for fracture network formation [43]. According to the mechanical model shown in Figure 3, there is an elliptical fracture in the infinite field with the major semiaxis and the minor semiaxis . At infinity, the major semiaxis is acted on by pressure stress () and the minor semiaxis (); there is uniform pressure action in the fracture.

Boundary conditions:

At ,

At ,

Expressions of , , and can be obtained according to the solutions for plane problems of elastic mechanics. Stress distribution on the fracture surface is expressed as follows: where .

Because , , and putting into Equation (8), we obtain the following:

According to the elastic failure criterion, setting , we obtain the following:

Therefore, when the value of net pressure in a fracture is larger than , the fracture in the rock body breaks and forms branch fractures, and then, a fracture network is formed.

To sum up, mechanical condition requirements for a fracture network in a naturally fractured reservoir lie in net pressure in the constructed fracture exceeding the horizontal principal stress difference of the reservoir. For a fracture network in a reservoir without natural fracture development, mechanical condition requirements lie in the fracture breaking in the rock body and the value of net pressure in the fracture being larger than the sum of horizontal principal stress difference and reservoir tensile strength.

3. Continuous Fracture Network Modeling

The physical model of network fractures is shown in Figure 4, which is composed of two clusters of orthogonal equally spaced vertical fractures. They are elliptical in the transverse direction, but the fracture height is not constant in the longitudinal direction, so they can be asymmetric in the upper and lower directions. The advantage of this physical model is that when the distribution of natural fractures is not clear, natural fractures and other factors are contained in the model, and it is a reasonable approximation of fractures in complex networks.

3.1. Building Fracture Network Modeling

In order to establish the mathematical model of the fracture network, the coordinate system was established. The origin of the coordinates was taken as the shaft injection point, the axis was the direction of the maximum horizontal principal stress, and the axis was the direction of the minimum horizontal principal stress. Suppose that the half-length of elliptical network fracture on axis and axis is and , respectively, and the fracture spacing on axis and axis is and , respectively, then the number of fractures parallel to axis and axis is and , where

According to the elliptic equation, the half-length of the th fracture parallel to the axis is

The half-length of the th fracture parallel to the axis is

In the process of network fracture propagation, the fluid flows from the center of the ellipse to the edge of the ellipse and presents laminar flow in the fracture. The mathematical model describing network fracture propagation is established by considering fracture fluid filtration and fracture height variation.

3.1.1. Volumetric Equilibrium Equation

According to the principle of mass conservation and ignoring fluid compressibility and initial filtration loss of fracturing fluid, the volume of fracturing injection is equal to the sum of the volume of network fractures and the filtration loss volume of fracturing fluid, that is where and are the volume of fractures parallel to axis and axis, respectively. and are the filtration volume of fracturing fluid for fractures parallel to axis and axis, respectively.

According to the morphology and composition of the fracture network, the cross-section of fractures parallel to the axis is approximately elliptical, and the volume of fractures is where is the fracture width of the th fracture parallel to the axis. is the fracture height of the th fracture parallel to the axis.

Similarly, the volume of the fracture parallel to the axis is where is the fracture width of the th fracture parallel to the axis. is the fracture height of the th fracture parallel to the axis.

The volume of the fracture is equal to the sum of the volume of fractures parallel to the axis and the volume of fractures parallel to the axis, that is:

According to the morphology and composition of the fracture network, the fracturing fluid filtration volume of fractures parallel to the axis can also be obtained by the following equation:

The fracturing fluid filtration volume of the fracture parallel to the axis is

3.1.2. Fracture Width Equation

In the th fracture parallel to the axis, considering the variation of crustal stress in the longitudinal direction, when the fracture penetrates the upper and lower layers (), the net pressure distribution in the fracture is where

According to the elastic deformation theory, the net pressure in the fracture is decomposed into the sum of even stress distribution and odd stress distribution, and then, the fracture width profile is calculated by the England and Green formulas. After deduction, the fracture width is

When the fracture does not penetrate the overlying and underlying strata () or the longitudinal variation of crustal stress is ignored, the fracture width of the th fracture parallel to the axis is

Similarly, when the fracture penetrates the upper and lower layers () or the longitudinal variation of crustal stress is ignored, the fracture width of the th fracture parallel to the y axis is where

When the fracture does not penetrate the upper and lower layers () or the longitudinal variation of crustal stress is ignored, the fracture width of the th fracture parallel to the axis is

3.1.3. Fracture Height Equation

According to the theory of linear elastic fracture mechanics, when the fracture penetrates the upper and lower layers (), the stress intensity factor at the upper and lower ends of the th fracture parallel to the axis is

According to the net pressure distribution in the fracture and let and , the fracture height equation of the th fracture parallel to the axis is

When the fracture does not penetrate the upper and lower layers () or the longitudinal variation of crustal stress is ignored, the fracture height equation of the th fracture parallel to the axis is

Similarly, when the fracture penetrates the upper and lower layers (), the fracture height equation of the th fracture parallel to the axis is

When the fracture does not penetrate the upper and lower layers () or the longitudinal variation of crustal stress is ignored, the fracture height equation of the th fracture parallel to the axis is

3.1.4. Fluid Flow Equation

According to the mechanics through porous media, the fluid flow equation in the th fracture parallel to axis is where is the axial ratio of the axis to the axis, . is the complete elliptic integrals of the second kind, . is the permeability of the th fracture parallel to the axis.

Considering the flow in the fracture is in laminar, the permeability of the th fracture parallel to the axis is

Substitute into the fluid flow equation, it can be obtained:

Considering the change of porosity with pressure, it is represented by rock compression coefficient, that is

Integrate the above equation and the porosity can be obtained by: where is the initial pressure.

Expand the above equation by series, omit the infinitesimal of higher order, and get so get

Thus, in the th fracture parallel to the axis, the fluid flow equation is where

The initial conditions and boundary conditions are

Similarly, in the th fracture parallel to the axis, the fluid flow equation is where

The initial conditions and boundary conditions are

3.1.5. Coupling Conditions

Each fracture of the network (whether parallel to axis or axis) is relatively independent and interrelated, constituting a complex network fracture system. Each fracture is described by its own fluid flow equation, fracture width equation, and fracture height equation and is controlled by volume balance equation. At the same time, there is coupling condition; that is, the pressure at the origin of coordinates is equal, expressed as

3.2. Numerical Solution of the Fracture Network Model
3.2.1. Difference Method for Fluid Flow Equations

It is a nonlinear partial differential equation for the fluid flow equation of the th fracture parallel to the axis. We use the implicit scheme difference method to solve it, and the nonlinear problem is solved by the coefficient freezing method. Thus, the left end of the fluid flow equation is discretized as follows:

The right end of the fluid flow equation is discretized as follows: where is the time step and is the length step in the direction.

Thus, the difference equation of the fluid flow equation is

These are tridiagonal linear equations with dominant diagonal, which can be solved by the catch-up method combined with their initial conditions and boundary conditions.

Similarly, for the th fracture parallel to the axis, the difference equation of the fluid flow equation (3-108) can be expressed as follows: where is the length step in the direction.

3.2.2. Forward Algorithm for the Numerical Solution

The main purpose of establishing and solving the fracture network model is known basic parameters (such as elastic modulus , Poisson’s ratio and other rock parameters, fracturing fluid viscosity , injection rate , injection time , and other construction parameters); by solving the model, the fracture network parameters can be obtained (such as half-length and , fracture width, fracture height, and ).

It is not hard to see that each fracture of the network is described by its own fluid flow equation, fracture height equation, and fracture width equation, which is a complete mathematical model. In other words, under the condition of known basic parameters, fracture spacing ( and ) and fracture half-length ( and ), the fracture length, fracture width, and fracture height and pressure in each fracture can be calculated by using this model. Therefore, by using the volume balance equation and coupling conditions (the pressure is equal at the intersection), the algorithm is established to obtain some parameters of the fracture spacing and the fracture half-length.

The forward algorithm in this paper is to obtain the half-length ( and ) of fracture network, the length, width, height, pressure, and of each fracture under the condition of known basic parameters and fracture spacing ( and ).

The specific step of the forward algorithm is as follows: (1)The coupling conditions are transformed to construct the following function:

As you can see that for any and , can be calculated by the function. Theoretical analysis and practical calculation show that if , then =. Therefore, for a fixed value , according to the condition , it is easy to calculate the value of by iterative method. At this time, the iteration formula is designed as follows: where is the iteration factor which is set to 0.1 during trial calculation. (2)According to the volume equilibrium equation, the following function is constructed:

It is not difficult to see that for any value of , since the value of is determined by the condition , is just a function of . Theoretical analysis and practical calculation show that if , then the value of reaches volumetric balance. Therefore, according to the condition , it is easy to calculate the value of by iterative method. At this time, the iteration formula is designed as follows: (3)The forward algorithm advances step by step by time step . Because the implicit difference method is used to solve the fluid flow equation, the numerical calculation is unconditionally stable. However, the smaller the time step, the higher the accuracy of the calculation result

To sum up, the forward algorithm designed is shown in Figure 5.

3.2.3. Inverse Algorithm for the Numerical Solution

The forward algorithm for the numerical solution of the network fracture model requires the fracture spacing ( and ) to be input as known parameters. Usually, the forward algorithm is used to check and invert the microseismic data of nearby wells, borrowing the fracture spacing of adjacent wells.

In fact, parameter inversion is the inverse problem of parameter forward modeling, and the inversion algorithm is different with different parameters. The inversion algorithm in this paper is to obtain the fracture spacing ( and ), the fracture length, the fracture width, the fracture height, the pressure in the fracture, and of each fracture through the algorithm under the condition that the basic parameters and the half-length ( and ) of the fracture network are known.

The inversion algorithm of the fracture network numerical solution is similar to the forward algorithm, which is briefly described here. Take the coupling conditions of equal pressure at the origin of coordinates and construct the following function:

The following iterative formula is adopted:

At the same time, according to the volume equilibrium equation, the following function is constructed:

The following iterative formula is adopted:

Therefore, the process of the inversion algorithm is shown in Figure 6.

4. Results and Discussion

4.1. Application

Take the well H03-11 as an example. The proposed continuous fracture network mechanics and model from this research have been applied to this well. The main parameters used during the simulation process for well H03-11 are listed in Tables 2 and 3.

The formation of this well belongs to a naturally fractured developed reservoir, with the value of net pressure in the fracture larger than . According to the mechanical mechanism of fracture network formation proposed in this paper, well H03-11 can form a fracture network. Therefore, the continuous fracture network model established in this paper is appropriate for simulations of fracture propagation within well H03-11. Using the continuous fracture network model, the growth of fracture height was calculated. As shown in Figure 7, the heights of the upper and lower extensions are 19.53 m and 11.71 m, respectively. Therefore, the height of the fracture is 31.24 m.

and are the half-length of the fracture in the and directions, respectively. is the fracture height. The calculated results of the continuous fracture network are shown in Figure 8. is 335.42 m, is 105.72 m, and is 31.24 m.

4.2. Verification and Discussion

The height of a fracture network changes with time as shown in Figure 9. The total injection time of well H03-11 is 126 minutes. At the early time stage, the fracture height increases rapidly. The later-time fracture height grows slowly. In the production process, the height changes, so considering height variation, the simulation results are more accurate.

We discuss the changes of , , and with different times; the results are shown in Figure 10. It can be seen from Figure 10 that the and increase rapidly at the early time stage and grow slowly at later-time. The growth of is S-shaped; that is, the first growth is slow, the middle growth is fast, and the latter is slow.

We also discussed the influence of fracturing fluid viscosity, filtration coefficient, and injection rate on , , and . The results are shown in Figure 11. It can be seen from Figure 11 that (1) the higher the viscosity and injection rate are, the better the effect of reservoir simulation is. The larger the filtration coefficient is, the worse the effect of reservoir simulation is. (2) Viscosity has little effect on and , but a greater effect on , only because viscosity has a greater effect on fracture height. (3) The injection rate has a great influence on , indicating that the injection rate should be increased in field case under the premise of ensuring safety.

Microseismic monitoring technology is the main method for accurate fracture network design. In order to verify the accuracy of the continuous fracture network model in this paper, the fracture was monitored by microseismic monitoring technology; the results are shown in Figure 12. An industry-accepted technology, Meyer, was employed to predict the dimension of fractures with the same well H03-11 data. Calculated by the forward algorithm of continuous fracture network model, Meyer and the microseismic monitoring results are compared in Table 4.

As shown in Table 4, the prediction accuracy of the continuous fracture network is higher than that of Meyer. This is because the influence of fracture height and filtration coefficient is not taken into account in the Meyer. From the discussion in Figure 11, we know that these two factors have a great influence on the prediction results.

In order to further prove the prediction accuracy of the continuous fracture network model in the field case, we applied this model to a horizontal well TP12-06. The model in this paper is only used to predict the size of fractures, and the directions of fractures are calculated by the pressure field model. The results of microseismic monitoring of TP12-06 are shown in Figure 13(a). The comparison between the prediction results of microseismic monitoring and continuous fracture network model is shown in Figure 13(b) and Table 5. The average error of prediction of fracture length and fracture width is 3.64% and 5.70%, respectively.

5. Conclusion

This work contributes to a novel continuous fracture network model. The findings from this study allow the following conclusions to be drawn: (1)By comprehensively considering the development situations of natural fractures, this study provides mechanical condition requirements related to fracture network formation. The mechanical condition requirements for a fracture network in a naturally fractured reservoir lie in net pressure in a constructed fracture exceeding the horizontal principal stress difference of the reservoir. For a fracture network in a reservoir without natural fracture development, the condition requirements lie in the fracture breaking in the rock body and the value of net pressure in the fracture being larger than the sum of horizontal principal stress difference and reservoir tensile strength(2)A continuous fracture network model was established based on the mechanical mechanism for fracture network formation. A geometric fracture network model was constructed with the main fracture and branch fractures, mathematical equations were established to simulate geometric parameters of a shale gas continuous fracture network, and related calculation methods were derived(3)In this study, a numerical simulation program of a continuous fracture network was developed, overall consideration was made for an increase in fracture height and filtration of fracturing liquid, and a simulation study was performed on a continuous fracture network. Actual microseismic monitoring data verifies the accuracy of this model. The simulation results are of great significance for fracturing design optimization of shale gas fracture networks

6. Recommendations

In future research, the following aspects can be further explored: (1)The junction of fracture is assumed to be noninterference, but there is fracture interference in the actual situation, so this study needs to be added in the future(2)It is assumed that fractures are perpendicular to each other. However, for an oblique fracture, how to establish a model in accordance with the real situation is worthy of further discussion

Nomenclature

:Maximum horizontal principal stress, MPa
:Minimum horizontal principal stress, MPa
:Angle from the maximum horizontal principal stress, dimensionless
:Shear stress, MPa
:Normal stress, MPa
:Time, min
:Pressure in fracture, MPa
:Normal net pressure in fracture, MPa
:Permeability, 10−3 μm2
:Porosity, dimensionless
:Length of a natural fracture, m
:Coefficient of friction on the surface of a natural fracture, dimensionless
:Effective shear stress, MPa
:Tensile strength, Pa
:Half-length of elliptical network fracture on axis, m
:Half-length of elliptical network fracture on axis, m
:Fracture spacing on axis, m
:Fracture spacing on axis, m
:Volume of fractures parallel to axis, m3
:Volume of fractures parallel to axis, m3
:Filtration volume of fracturing fluid for fractures parallel to axis, m3
:Filtration volume of fracturing fluid for fractures parallel to axis, m3
:Fracture width of the th fracture parallel to the axis, m
:Fracture height of the th fracture parallel to the axis, m
:Fracture width of the th fracture parallel to the axis, m
:Fracture height of the th fracture parallel to the axis, m
:Fracture thickness, m
:Pressure of the th fracture parallel to the axis, MPa
:Pressure of the th fracture parallel to the axis, MPa
:Stress intensity factor of upper fracture, MPa·m0.5
:Stress intensity factor of lower fracture, MPa·m0.5
:Fracture toughness of upper rock, MPa·m0.5
:Fracture toughness of lower rock, MPa·m0.5
:Fracture toughness, MPa·m0.5
:Axial ratio of the axis to the axis, dimensionless
:Complete elliptic integrals of the second kind, dimensionless
:Rock compressibility, Pa-1
:Initial pressure, MPa
:Time step, min
:Length step in the direction, m
:Length step in the direction, m
:Iteration factor which is set to 0.1 during trial calculation, dimensionless.

Data Availability

The data used to support the findings of this study are included in the article.

Conflicts of Interest

The authors confirm that this article content has no conflict of interest.

Acknowledgments

This work was supported financially by the Science and Technology Cooperation Project of the CNPC-SWPU Innovation Alliance and Study on Production Law and Technology Policy of Coalbed Methane Reservoir (2017E-1405).