A Novel Continuous Fracture Network Model: Formation Mechanism, Numerical Simulation, and Field Application

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 ﬁ eld application. However, some important factors, such as fracture height variation and ﬁ ltration 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, ﬁ ltration loss, ﬂ uid ﬂ ow, 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 ﬁ eld 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%.


Introduction
Hydraulic fracturing has been used commonly in the development of shale gas reservoirs to initiate several dominant hydraulic fractures from the wellbore [1][2][3][4]. 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 [5][6][7][8][9]. 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 [10][11][12].
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 threedimensional 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 net-work 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.

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.

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 (σ H ) and minimum horizontal principal stress (σ h ) of the far-field, and its angle from the maximum horizontal principal stress is θð0 < θ < ðπ/2ÞÞ.
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: 2 Geofluids 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 θ = ðπ/2Þ, 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 (τ e ) 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 θ = ðπ/2Þ, the net pressure is the largest, and p max is as follows: The natural fracture surface is unbounded (τ 0 = 0). 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.  [26,27] 2009 CFN Semianalytical solution model Meyer et al. [28] 2011 CFN Relatively complete model and software Nejadi et al. [31] 2017 DFN History matching and uncertainty quantification Yao et al. [32] 2020 DFN Role of natural fracture characteristics Hyman and Dentz [33] 2021 DFN Transport upscaling under flow heterogeneity and matrix-diffusion

Geofluids
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.

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 L f and the minor semiaxis w. At infinity, the major semiaxis is acted on by pressure stress (σ H ) and the minor semiaxis (σ h ); there is uniform pressure action in the fracture.
Boundary conditions: At y = 0, jxj < l f , At ffiffiffiffiffiffiffiffiffiffiffiffiffi Expressions of σ x , σ y , and τ xy can be obtained according to the solutions for plane problems of elastic mechanics. Stress distribution on the fracture surface is expressed as follows: where m = ðl f − wÞ/ðl f + wÞ.
Because l f ≫ w, m ≈ 1, and putting into Equation (8), we obtain the following: According to the elastic failure criterion, setting σ θ = −S t , we obtain the following: Therefore, when the value of net pressure in a fracture is larger than ðσ H − σ h Þ + S t , 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.

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  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 x axis was the direction of the maximum horizontal principal stress, and the y axis was the direction of the minimum horizontal principal stress. Suppose that the half-length of elliptical network fracture on x axis and y axis is L x and L y , respectively, and the fracture spacing on x axis and y axis is d x and d y , respectively, then the number of fractures parallel to x axis and y axis is 2n y + 1 and 2n x + 1, where According to the elliptic equation, the half-length of the j th fracture parallel to the x axis is The half-length of the ith fracture parallel to the y 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.

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 V f x and V f y are the volume of fractures parallel to x axis and y axis, respectively. V lx and V ly are the filtration volume of fracturing fluid for fractures parallel to x axis and y axis, respectively. According to the morphology and composition of the fracture network, the cross-section of fractures parallel to the x axis is approximately elliptical, and the volume of fractures is where w xj is the fracture width of the jth fracture parallel to the x axis. h xj is the fracture height of the jth fracture parallel to the x axis. Similarly, the volume of the fracture parallel to the y axis is where w yi is the fracture width of the ith fracture parallel to the y axis. h yi is the fracture height of the ith fracture parallel to the y axis. The volume of the fracture is equal to the sum of the volume of fractures parallel to the x axis and the volume of fractures parallel to the y axis, that is: According to the morphology and composition of the fracture network, the fracturing fluid filtration volume of fractures parallel to the x axis can also be obtained by the following equation: The fracturing fluid filtration volume of the fracture parallel to the y axis is In the jth fracture parallel to the x axis, considering the variation of crustal stress in the longitudinal direction, when the fracture penetrates the upper and lower layers (h xj ≥ H p ), the net pressure distribution in the fracture is

Geofluids
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 (h xj < H p ) or the longitudinal variation of crustal stress is ignored, the fracture width of the jth fracture parallel to the x axis is Similarly, when the fracture penetrates the upper and lower layers (h yi ≥ H p ) or the longitudinal variation of crustal stress is ignored, the fracture width of the ith fracture parallel to the y axis is where When the fracture does not penetrate the upper and lower layers (h yi < H p ) or the longitudinal variation of crustal stress is ignored, the fracture width of the ith fracture parallel to the y axis is

Fracture Height Equation.
According to the theory of linear elastic fracture mechanics, when the fracture penetrates the upper and lower layers (h xj ≥ H p ), the stress intensity factor at the upper and lower ends of the jth fracture parallel to the x axis is According to the net pressure distribution in the fracture and let K Iu = K ICu and K Il = K ICl , the fracture height equation of the jth fracture parallel to the x axis is When the fracture does not penetrate the upper and lower layers (h xj < H p ) or the longitudinal variation of crustal stress is ignored, the fracture height equation of the jth fracture parallel to the x axis is Similarly, when the fracture penetrates the upper and lower layers (h yi ≥ H p ), the fracture height equation of the ith fracture parallel to the y axis is When the fracture does not penetrate the upper and lower layers (h yi < H p ) or the longitudinal variation of crustal stress is ignored, the fracture height equation of the ith fracture parallel to the y axis is

Fluid Flow Equation.
According to the mechanics through porous media, the fluid flow equation in the jth fracture parallel to x axis is where ξ is the axial ratio of the x axis to the y axis, ξ = L x /L y . η is the complete elliptic integrals of the second kind, η = Ð π/2 0 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 − ð1 − ξ 2 Þ sin 2 s p ds. k f xj is the permeability of the jth fracture parallel to the x axis.
Considering the flow in the fracture is in laminar, the permeability of the jth fracture parallel to the x 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 p 0 is the initial pressure. Expand the above equation by series, omit the infinitesimal of higher order, and get so get Thus, in the jth fracture parallel to the x axis, the fluid flow equation is where The initial conditions and boundary conditions are Similarly, in the ith fracture parallel to the y axis, the fluid flow equation is where The right end of the fluid flow equation is discretized as follows: where Δt is the time step and Δx is the length step in the x direction.
Input basic parameters, fracture spacing d x and d y . Take the initial value L x 0 and L y 0 . Let t = 0, n = 0 and given precision .
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 ith fracture parallel to the y axis, the difference equation of the fluid flow equation (3-108) can be expressed as follows: where Δy is the length step in the y direction.

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 E, Poisson's ratio ν and other rock parameters, fracturing fluid viscosity μ, injection rate q, injection time T, and other construction parameters); by solving the model, the fracture network parameters can be obtained (such as half-length L x and L y , fracture width, fracture height, and V f ).
Input basic parameters, fracture spacing L x and L y . Take the initial value d x 0 and d y 0 . Let t = 0, n = 0 and given precision .  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 (d x and d y ) and fracture half-length (L x and L y ), 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 (L x and L y ) of fracture network, the length, width, height, pressure, and V f of each fracture under the condition of known basic parameters and fracture spacing (d x and d y ).
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 L x and , J d1 can be calculated by the function. Theoretical analysis and practical calculation show that if J d1 = 0, then L x =L y . Therefore, for a fixed value L x , according to the condition J d1 = 0, it is easy to calculate the value of L y 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 L x , since the value of L y is determined by the condition J d1 = 0, J d2 is just a function of L x . Theoretical analysis and practical calculation show that if J d2 = 0, then the value of L x reaches volumetric balance. Therefore, according to the condition J d2 = 0, it is easy to calculate the value of L x by iterative method. At this time, the iteration formula is designed as follows:   10 Geofluids (3) The forward algorithm advances step by step by time step Δt. 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.

Inverse Algorithm for the Numerical Solution.
The forward algorithm for the numerical solution of the network fracture model requires the fracture spacing (d x and d y ) 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 (d x and d y ), the fracture length, the fracture width, the fracture height, the pressure in the fracture, and V f of each fracture through the algorithm under the condition that the basic parameters and the half-length (L x and L y ) 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:

Geofluids
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.

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 ðσ H − σ h Þ + S t . 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  Figure 9. The total injection time of well H03-11 is 126 minutes. At the early time stage, the fracture height increases rapidly. The latertime 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 L x , L y , and V f with different times; the results are shown in Figure 10. It can be seen from Figure 10 that the L x and L y increase rapidly at the early time stage and grow slowly at later-time. The growth of V f is Sshaped; 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 L x , L y , and V f . 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 L x and L y , but a greater effect on V f , only because viscosity has a greater effect on fracture height. (3) The injection rate has a great influence on V f , 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.

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 13 Geofluids 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

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