Production Decline Analysis and Hydraulic Fracture Network Interpretation Method for Shale Gas with Consideration of Fracturing Fluid Flowback Data

After multistage hydraulic fracturing of shale gas reservoir, a complex fracture network is formed near the horizontal wellbore. In postfracturing flowback and early-time production period, gas and water two-phase flow usually occurs in the hydraulic fracture due to the retention of a large amount of fracturing fluid in the fracture. In order to accurately interpret the key parameters of hydraulic fracture network, it is necessary to establish a production decline analysis method considering fracturing fluid flowback in shale gas reservoirs. On this basis, an uncertain fracture network model was established by integrating geological, fracturing treatment, flowback, and early-time production data. By identifying typical flow-regimes and correcting the fracture network model with history matching, a set of production decline analysis and fracture network interpretation method with consideration of fracturing fluid flowback in shale gas reservoir was formed. Derived from the case analysis of a typical fractured horizontal well in shale gas reservoirs, the interpretation results show that the total length of hydraulic fractures is 4887.6m, the average half-length of hydraulic fracture in each stage is 93.4m, the average fracture conductivity is 69.7mD·m, the stimulated reservoir volume (SRV) is 418 × 104 m3, and the permeability of SRV is 5:2 × 10−4 mD. Compared with the interpretation results from microseismic monitoring data, the effective hydraulic fracture length obtained by integrated fracture network interpretation method proposed in this paper is 59% of that obtained from the microseismic monitoring data, and the effective SRV is 83% of that from the microseismic monitoring data. The results show that the fracture length is smaller and the fracture conductivity is larger without considering the influence of fracturing fluid.


Introduction
Due to the extremely low permeability of shale gas reservoirs, large-scale hydraulic fracturing treatment is often needed to achieve economic productivity. After multistage of fracturing, a complex fracture network is created near the horizontal wellbore. The analysis of flowback after hydraulic fracturing and production decline characteristics at the early stage of production provides a new idea for the evaluation of post-fracturing effect [1][2][3].
Since most shale gas wells exhibit long-term linear/bilinear flow characteristics in the early stage, the parameter inversion of production decline analysis mainly analyzes dynamic data through the approximate solution of linear/bilinear flow stage, establishes a suitable mathematical model, and draws pressure/production and its derivative curve, using local approximate solutions to obtain interpretation of key fracture parameters such as fracture half-length and conductivity.
The inversion of key parameters of complex fracture network in shale gas reservoir is a two-phase seepage problem. On the one hand, after hydraulic fracturing in many shale gas wells, a large amount of fracturing fluid is retained in the reservoir, causing reservoir pollution, and long-term gas-water two-phase flow characteristics appear. On the other hand, with the gradual maturity of fracturing horizontal well technology, shale condensate gas, volatile oil, and black oil resources in the Eagle Ford Basin in North America have also been industrialized [4].
Shale gas is produced by large-scale hydraulic fracturing, and tens of thousands of cubic meters of fracturing fluid are usually injected into the reservoir to obtain a complex fracture network. However, only 25-60% of the fracturing fluid is produced during the flowback and production process, causing a large amount of fracturing fluid to stay in the reservoir. Therefore, the early dynamic analysis of shale gas wells needs to consider both the gas-water two-phase seepage of the reservoir and the reservoir damage caused by fracturing fluid retention.
At present, some models used for fracturing fluid flowback data analysis consider the gas-water two-phase seepage problem in shale reservoirs. Clarkson and Williams-Kovacs [5] only considered the fluid in the fracture, ignoring the supply of the matrix to the fracture, and analyzed the early flowback dynamics. They found that the use of this model can predict the production dynamics in the early months, but not long-term. Dynamic analysis. Williams-Kovacs and Clarkson [6] considered the gas supply of the matrix to the fracture system and extended the flowback data analysis model of Clarkson and Williams-Kovacs [5].
For the two-phase seepage of oil and gas in shale reservoirs, some scholars have used the slab fracture model to propose methods for predicting the production of twophase oil and gas. Behmanesh et al. [7] treat saturation as a function of pressure, linearize the model by defining a two-phase pseudopressure function, and obtain a semianalytical solution for productivity. Due to the difficulty in solving the relationship between saturation and pressure, Zhang and Ayala [8] and Tabatabaie and Pooladi-Darvish [9] used the Boltzmann transformation method to obtain the self-model solution of the oil and gas two-phase seepage model. However, these models are only applicable to the situation before the pressure reaches the boundary. When the pressure reaches the boundary, the saturation and pressure no longer show a fixed relationship.
In summary, although many scholars have proposed complex fracture network inversion models for shale gas reservoirs, the existing gas-water two-phase dynamic analysis methods are mainly limited to early flowback simulations and cannot be used for long-term dynamic analysis. Therefore, considering the production dynamic analysis of gas-water two-phase flow under complex fracture network conditions has become a technical difficulty.
In this paper, based on the complex fracture network formed after hydraulic fracturing of shale gas reservoirs [10,11], a gas-water two-phase flow model was established considering complex fracture network conditions and integrates geological-fracturing-production data. Combining the analysis results of the flow stage and the automatic history matching correction fracture network model, a set of shale gas reservoir production decline analysis and fracture network inversion methods considering fracturing fluid flowback are proposed.

The Model and Method
2.1. The Physical Model. In shale gas reservoir, the horizontal wells form complex fracture network after stimulated reservoir volume. The flowback and early-time production often shows two-phase flow. [12] In this paper, the discrete fracture network model is used to characterize the artificial fractures after fracturing, considering the water phase flow (fracturing fluid) in the discrete fracture network in the initial state and the matrix is mainly single-phase gas flow to establish fracturing flowback in shale gas reservoirs. The physical model of the gas-water two-phase flow in the process is shown in Figure 1.
As shown in Figure 2 above, this model considers the fractures after volume fracturing reformation as plane fractures with fully extended double wings, and the interconnected plane fractures form a complex fracture network. The discrete fracture network model is used to simulate the complex fracture network system after shale gas volume fracturing. The fractures are divided into several microelement segments, and the gas-water two-phase flow in any fracture system is considered, without considering the capillary force and gravity. The flowchart is shown in Figure 3.

The Gas-Water Two-Phase Equation of Discrete
Fracture System. In the process of fracturing fluid flowback in shale gas reservoirs, the fluid from the reservoir flows first from the matrix into the fracture and then into the wellbore. [13] The gas-water two-phase flow equation in the discrete fracture system is established considering the flowback of fracturing fluid and the gas-water two-phase flow in the fractures during the early production process and the singlephase flow of shale gas in the matrix, ignoring the pressure drop in the wellbore.
During fracturing fluid flowback and early-time production, the partial differential equation of gas phase transient flow in fully penetrated hydraulic fractures can be given by the following equation: ∂ ∂x where k f represents the absolute permeability of fracture, k rg represents gas phase relative permeability, φ f represents the fracture porosity, μ g represents gas viscosity, B g represents the formation volume factor of gas, s g represents the gas saturation, and p f represents the fracture pressure of the system. The right-hand cumulative term in Equation (1) can be written as 1 ∂t By substituting the definition of fracture porosity and gas volume coefficient into the above equation, it can be obtained:

Geofluids
∂ ∂x It can be seen from Equation (3) that the cumulative term in the diffusion equation is affected by gas-water two-phase fluid flow, fracture propagation, and gas saturation variation during fracturing fluid flowback and early production. In this paper, the following equation is introduced to define the compressibility of fluid.
In Equation (4), the compressibility coefficient of the fluid depends on the fluid saturation and the pressure in the fracture. Consider that the permeability in the fracture system is large and the pressure gradient is usually very small [14,15]. This paper proposes the average saturation of the fluid and the compressibility coefficient under the average pressure to simplify the equation, and the results are shown as follows: Equation (5) shows that the flow equation under average pressure and average saturation is similar to the flow equation for single-phase fluid. For the convenience of solving, the equation still needs to be linearized.

Material Balance Equation.
In this paper, the material balance method s [16][17][18][19] is used to solve the gas-water twophase flow equation in the discrete fracture system. In order to linearize the equation, a modified pseudotime term is introduced in this paper, as shown in the following equation: Above, ðc * mt Þ i represents the total compressibility of gas phase under the condition of the initial state, and c * mt represents the fracture system under the condition of average water saturation and pressure coefficient of the gas phase.
The modified pseudotime term (6) was substituted into the gas-water two-phase equation and further simplified as Based on the research results of Mattar and Anderson [15,16], a linear analysis formula considering the linear flow in the gas phase can be obtained, assuming that the fracture ends are closed and the matrix continuously supplies gas to the fracture and considering the variable working system conditions.

Geofluids
The gas phase is as follows: Besides, It can be seen from Equation (8) that in the fracturing fluid flowback and early production stage, the gas-water two-phase flow in the fracture system presents a linear relationship as shown below: where N w represents the initial volume of water and s wi represents the initial water saturation and usually takes 1 to represent the amount of fracturing fluid produced.
In the early stage of fracturing fluid flowback and production, the unsteady flow equations of water and gas in the fracture can be expressed by the following equation: The average viscosity and average volume coefficient of gas can be calculated through the average pressure in the fracture. Combined with Equations (12) and (13), the relationship between the gas-water production ratio and the relative permeability (gas-water ratio) can be obtained as follows: In combination with Equations (11) to (14), the rela-tionship between average water saturation and average fracture pressure can be determined.

Integrated Interpretation Method for Key Parameters of Fracture Network
In order to reduce the uncertainty of the key parameters of fractured horizontal wells in shale gas reservoirs and make full use of existing data for fracture characterization, this paper integrates multiple data such as cores, microseismic monitoring data, hydraulic fracturing construction data, and production performance data, which formed a set of integrated inversion methods for key parameters of fracture network of fractured horizontal wells in shale gas reservoirs.
The specific process is shown in Figure 3. The method mainly includes the steps of random modeling of natural fractures, artificial fracture propagation simulation, gaswater two-phase straight-line analysis, and history matching.

3.1.
Step 1: Natural Fracture Stochastic Modeling. Mayerhofer et al. [17] and Gamboa et al. [18] showed that natural fractures would reopen under hydraulic fracturing, and the energy of the reopened natural fractures would trigger seismic wave events. Therefore, the microseismic data was utilized to determine the location of natural fractures, combining well-logging interpretation results. We can obtain the direction, length, and number of natural fractures and the frequency distribution of the length and direction of each fracture. And then, the length and orientation of each natural fractures could be generated randomly according to the fracture parameter distribution.

3.2.
Step 2: Hydraulic Fracturing Simulation. Based on natural fractures, a hydraulic fracture propagation model is established. In this paper, each stage of fracturing is divided into three clusters. It is assumed that a hydraulic fracture would be generated in each cluster after hydraulic fracturing, and the amount of fracturing fluid and proppant in each cluster of the same stage is the same [20]. In the matrix, as we all know, hydraulic fractures propagate along the maximum horizontal principal stress until they intersect with a natural fracture. If a hydraulic fracture intersects a natural fracture, it extends along that natural fracture. In the process of fracture propagation, it is assumed that the hydraulic fracture propagates at the same rate in the matrix on both sides of the horizontal well. When one side intersects the natural fracture, the other side stops propagating in the matrix. This is because the fracture propagates more easily along the natural fracture than in the matrix. In addition, the total length of hydraulic fracture is proportional to the volume of fracturing fluid [20][21][22]. Based on the above process, the generated seam net can be obtained.
Above, X F,i is the total length of Class I hydraulic fracture (m), β F is coefficient, and V F,i is the total amount of stage i fracturing fluid used (m 3 ).

Step 3: Linear Flow Analysis of Gas and Water
Production Data. Based on the established natural fracture and artificial fracture propagation models, the gas-water two-phase flow in the discrete fracture network was considered. Then, the analytical model of gas-water two-phase could be established. And finally, the typical flow stage and the key parameters (permeability of stimulated reservoir, volume of stimulated area, fracture permeability and fracture length) of the complex fractures were determined.

3.4.
Step 4: History Matching. Based on the discrete fracture network model, the linear analysis method proposed in this paper is utilized to adjust the key fracture network parameters, such as discrete fracture length, conductivity of discrete fracture, permeability of reformed zone, and volume of reformed zone, to match the production dynamic data (daily gas production and daily water production).

Field Application
In this article, a multistaged fracturing well in a shale gas window in China is analyzed. The reservoir depth is about 2,200 m, the effective thickness is about 15 m, the initial pressure of the reservoir is 28.5 MPa, the reservoir temperature is 358.1 K, and the initial gas saturation of the reservoir is 0.6. The horizontal well is hydraulically fractured into 26 sections, each stage fracturing 3 clusters, the production performance is shown in Figure 4. The logging interpretation results show that the matrix porosity is 0.064. Statistical imaging logging and core data to obtain azimuth distribution of natural fractures as the Table 1 shown. The maximum principal stress orientation is basically perpendicular to the horizontal well direction.

4.1.
Step 1: Reconstruction of Complex Fracture Network. Based on microseismic data from the fractured horizontal well in the shale gas reservoir, each microseismic data divides the fracture into two sections, assuming that the ratio of the lengths of the two sections obeys a uniform distribution. Combining the core and image logging interpretation results, a complex fracture network is randomly generated as shown in Figure 5.

4.2.
Step 2: Linear Analysis of Gas and Water Production Data. Based on the gas-water two-phase linear analysis method proposed above, the gas-water two-phase production dynamic data were combined with material balance time, normalized pseudopressure and normalized output to obtain the gas-phase flow characteristic curve, as shown in Figures 6, 7 and 8. It can be seen from the gas phase flow characteristic curve that the gas phase bilinear flow does not appear at the early stage of production, and the linear flow is obvious. After 2.5 months, most of the fracturing fluid has been discharged, and the fracture conductivity is high or the fracture half- 6 Geofluids length is large. After 20 months, interference flow between fractures appeared in the gas phase.
Combining the linear analysis method proposed in the previous article, the linear flow stage of the gas phase is fitted, and the fitted linear relationship (y = 0:0418x + 0:1025, R 2 = 0:816), based on the linear analysis method given above, by identifying the linear flow in the gas phase. At this stage, the average fracture half-length is 96.7 m, the fracture conductivity is 74.5 mD·m, the average width of the reformed area is 48.7 m, and the permeability of the reformed area is 6:01 × 10 −4 mD.

4.3.
Step 3: Production History Matching. According to the integration fracture network and the discrete fracture network model, the discrete fractures length, the conductivity of discrete fracture, and the permeability of the stimulated area were adjusted (as shown in Table 2) to fit the production performance data of the multistage fractured horizontal well in this shale gas reservoir, the basic grid size was set to 50 m × 50 m, and the grid was defined connected to the fracture to 25 m × 25 m. The fitting result is shown in Figures 9  and 10.
Combined with the fitting results of actual well production performance, the discrete fracture morphology, conductivity, stress sensitivity coefficient, and the scope and permeability of the stimulated zone were obtained. The closer it is to the wellbore, the higher the fracture conductivity and the smaller the stress sensitivity coefficient are, as shown in Figures 11 and 12.
The fitting results of the fracture network key parameters such as discrete fracture length, average fracture halflength, average fracture conductivity, and reconstructed fracture area were compared with microseismic, production decline analysis, and commercial software kappa, as shown in Table 3.

Geofluids
The integrated fracture network interpretation results show that the half-length of discrete fractures is 4887.6 m, the average half-length of fractures is 93.4 m, the average fracture conductivity is 69.7 mD·m, the reconstructed volume is 418 × 10 4 m 3 , the permeability is 5:2 × 10 −4 mD, and the porosity is 0.058. From the comparison results of key parameters of fracture network, the effective discrete fracture length obtained by the integrated fracture network interpretation method proposed in this paper is 59% of that of microseismic monitoring, and the effective reconstruction body is 83% of that of microseismic monitoring. Without considering the influence of fracturing fluid, smaller fracture length and larger fracture conductivity can be obtained.

Conclusions
(1) Based on the complex fracture network formed after hydraulic fracturing of shale gas reservoirs, a model considering gas-water two-phase flow in fractures is established absorbing geologicalfracturing-production data, and then, a comprehensive fracture network parameter inversion method for shale gas reservoir considering fracturing fluid flowback is developed by integrating flow stage analysis and automatic history fitting correction fracture network model (2) Taking a fractured horizontal well in a shale gas reservoir as an example, the discrete fracture halflength was obtained through the steps of random modeling of natural fractures, simulation of artificial fracture propagation, gas-water two-phase line analysis, and production history fitting. The results indicate that half-length of discrete fracture is 4887.6 m, the average fracture half-length is 93.4 m, the average fracture conductivity is 69.7 mD·m, the reconstructed volume is 418 × 10 4 m 3 , the permeability of the reconstructed area is 5:2 × 10 −4 mD, and the porosity of the reconstructed area is 0.058 (3) Compared with the microseismic interpretation results, the effective discrete fracture length obtained by the integrated fracture network interpretation method proposed in this paper is 59% of that obtained by the microseismic monitoring, and the effective SRV volume is 83% of that obtained by the microseismic monitoring. Without considering the influence of fracturing fluid, a smaller fracture length and larger fracture conductivity would be obtained

Data Availability
The data used to support the findings of this study are included within the article.