Numerical Simulation of Gas-Water Two-Phase Flow in Deep Shale Gas Reservoir Development Based on Mixed Fracture Modeling

In this paper, the development of a three-dimensional, two-phase fluid flow model (Modified Embedded Discrete Fracture Model) to study flow performances of a fractured horizontal well in deep-marine shale gas is presented. Deep-marine shale gas resources account for nearly 80% in China, which is the decisive resource basis for large-scale shale gas production. The dynamic characteristics of deep shale gas reservoirs are quite different and more complex. This paper uses the embedded discrete fracture model to simulate artificial fractures (main fractures and secondary fractures) and the dual-media model to simulate the mixed fractured media of natural fractures and considers the flow characteristics of partitions (artificial fractures, natural fractures, and matrix). Gas desorption is considered in the matrix. Different degrees of stress sensitivity are considered for natural and artificial fractures. Aiming at accurately simulating the whole production history of horizontal well fracturing, especially the dynamic changes of postfracturing flowback, a postfracturing fluid initialization method based on fracturing construction parameters (fracturing fluid volume and pump stop pressure) is established. The flow of gas and water in the early stage after fracturing is simulated, and the regional phase permeability and capillary force curves are introduced to simulate the process of flowback and production of horizontal wells after fracturing. The influence of early fracture closure on the gas-water flow is characterized by stress sensitivity. A deep shale gas reservoir of Sinopec was selected for the case study. The simulation results show it necessary to consider the effects of fractures and stress sensitivity in the matrix when considering the dynamic change of production during the flowback and production stages. The findings of this study can help for better understanding of the fracture distribution characteristics of shale gas, shale gas production principle, and well EUR prediction, which provide a theoretical basis for the effective development of shale gas horizontal well groups.


Introduction
With the rapid development and maturity of horizontal well and fracturing technologies, large-scale industrial production has been achieved for the shale oil and gas reservoirs. Based on domestic researches for the unconventional reservoirs in China, it is widely recognized that the deep-marine shale gas reservoirs account for nearly 80% [1][2][3], which is very high and decisively assures the resource basis for the large-scale production of shale gas. Among these deep-marine reservoirs, the shale gas covering area and resource amount in the Jiaoshiba-2 reservoir is particularly prominent. Till 2019, production plans and implementations have been proposed for the Pingqiao and Jiangdong plays [4]. In southern Sichuan, a credible estimation of shale gas resource can reach 10,400 billion m 3 , in which the amount of the Longmaxi play accounts for over 50%, i.e., 5,220 billion m 3 . Moreover, in the Longmaxi play, the resource amount of the reservoir the depth of which is larger than 3,500 meters accounts for nearly 94.1% [5], i.e., 4,910 billion m 3 .
In recent years, SINOPEC has made remarkable progress in the production of deep shale gas reservoirs in southern Sichuan. In the multiwell test for the Longmaxi play, industrial-scale gas production rate has been achieved, which can reach several hundreds of thousand cubic meters per day [6]. In spite of these industrial practical progresses, there still remains a lack of understanding on the storage and transport mechanisms embedded in the deep shale gas reservoir, which has typical high temperature, high pressure, and high stress. The above features lead to many practical problems, such as the difficulties of reservoir stimulation, the high flowback rate in the early production stage, and the complex gas-water relationship, which limits the application of conventional gas reservoir simulators for the production of deep shale gas reservoirs.
To solve the above issues, many simulation techniques have been proposed in recent decades, for example, the EDFM (Embedded Discrete Fracture Model) [7][8][9][10]. The EDFM technique is first proposed by Lee et al. in order to improve the flexibility of fracture characterization. The concept of multisegment wells is utilized to handle the complex fracture networks based on Peaceman's representation for transmissibility [11][12][13][14] and to construct the conductivity between connected reservoir grids. Different from the DFN (discrete fracture network) method in which a complex unstructured mesh is constructed to conform both fracture geometry and matrix grids, the EDFM utilizes the transmissibility function to describe the fluid transfer between fractures and matrix, and instead of introducing a complex unstructured mesh, regular convectional structured grids can be utilized for the discretization. Recently, the EDFM model has been proven to be a high-resolution simulation technique which is similar to the DFN method in unconventional reservoirs [15][16][17][18]. With special treatments, both the high-resolution characteristics of matrix grids and local fracture features can be well maintained during simulation, particularly for the extremely low permeable reservoirs. Basically, in the case of extremely low permeability, there exists a very long period for the unsteady diffusion process [19], which leads to a larger inconsistence error of the matrix-fracture transfer equation. This problem is also reported in previous DPDP (dual porosity and dual permeability) model-based literatures, which is more computationally efficient than the EDFM model and usually utilized for simulation with a large amount of small-scale dense fractures [20][21][22].
The flowback characteristics of the fracturing fluid are an important component of stimulation features in shale reservoirs [23][24][25]. A large amount of literatures have shown the importance of the role of the capillary effect in the redistribution process of the water and gas phases. During the close-in stage, the amount of penetrated water is highly related to the capillary force and soaking time [26,27]. The geomechanical response due to the stress release during the flowback period leads to the closure of microfractures and separates these microfractures from the stimulated hydraulic fractures, which enhanced the residual water amount in the reservoir. The Multiscale Finite Volume (MsFV) method is used to simulate highly heterogeneous porous media [28]. Above all, the state-of-the-art research on the fluid storage state after hydraulic fracturing is still far from complete [29], and more investigations for determining the initial pressure of water and gas phases, as well as saturation features, are still needed during reservoir simulations.
Therefore, in this paper, with full consideration of the complex fracture networks formed during the multisegment/cluster fracturing in shale reservoirs, the characteristic transport mechanisms embedded in the shale matrix, and the gas-water two-phase dynamic flow features, we aimed to construct a gas-water two-phase reservoir simulation approach oriented to the deep shale reservoir applications based on the EDFM method. With the proposed framework, it is possible for us to provide effective plans for the industrial development of deep shale gas reservoirs.
First, the methodology is presented. Then, a case study is carried out to show the accuracy of the method and discuss the well performance and two-phase flow mechanism among the fracture and matrix in deep-marine shale gas reservoir. (1) The shale reservoir is composed of three parts, i.e., shale matrix, natural fractures, and hydraulic/artificial fractures (2) Natural fractures and shale matrix are both regarded as continuum with anisotropic permeability tensor, i.e., k fh ≠ k fv (3) The hydraulic/artificial fractures are represented based on the EDFM method (4) Gas in the fracture networks is assumed to be free shale gas, which obeys Darcy's law (5) In the shale matrix, there exist two storage patterns for shale gas, i.e., adsorbed gas and free gas (6) The desorption behavior of the adsorbed gas in shale matrix can be described by Langmuir's isothermal adsorption law (7) The reservoir is assumed to be already reaching an equilibrium state before exploitation, and the adsorbed gas and free gas are also in a dynamic equilibrium state (8) The gas-water two-phase flow is assumed to be isothermal, and the gravity and capillary effects are considered 2.2. Governing Equations. In the development of shale gas (dry gas) reservoirs based on hydraulic fracturing and horizontal wells, there exist multiscale flow behaviors in the system of the matrix, natural fracture, and hydraulic fracture.

Methodology
To accurately describe this system, gas-slip, adsorption-2 Lithosphere desorption, and stress sensitivity effects should be considered for describing the flow in the matrix. Besides, there exists a flux transfer among the matrix, natural fractures, and hydraulic fractures. For this three-dimensional gaswater two-phase system, the related governing equations are listed as follows.
2.2.1. Two-Phase Flow in Matrix. For the matrix system, we consider two important microscale mechanisms, i.e., gas adsorption-desorption and phase-slip effects, and individual relative permeability and capillary force models are adopted to describe the multiphase flow in the matrix. Then, the governing equation of the gas phase and water phase can be written as where V L is the Langmuir volume, m 3 /t; p L is the Langmuir pressure, MPa; p gm is the pressure of the gas phase in the matrix, MPa; k rg is the relative permeability of the gas phase in the matrix; q gmnf is the gas flux transfer between the matrix and natural fractures; and q gmhf is the gas flux transfer between the matrix and hydraulic fractures.
where p wm is the pressure of water phase in matrix, MPa; k rw is the relative permeability of the water phase in the matrix; q wmnf is the water flux transfer between the matrix and natural fractures; and q wmhf is the water flux transfer between the matrix and hydraulic fractures.

Two-Phase Flow in Natural
Fractures. For the flow in natural fractures, we consider the stress sensibility, and individual relative permeability and capillary force models are adopted to describe the flow. The governing equations of gas and water phases are then written as where p gnf is the gas pressure in natural fractures, MPa; q gmnf is the gas flux transfer between the matrix and natural fractures; and q gnhf is the gas flux transfer between natural fractures and hydraulic fractures.
where p wnf is the gas pressure in natural fractures, MPa; q wmnf is the water flux transfer between the matrix and natural fractures; and q wnhf is the water flux transfer between natural fractures and hydraulic fractures.

Two-Phase Flow in Hydraulic
Fractures. Similar with the case in natural fractures, the governing equations describing the gas-water two-phase flow in hydraulic fractures can be written as where p ghf is the gas pressure in hydraulic fractures, MPa.
where p whf is the water pressure in hydraulic fractures, MPa. Note that the flux transfer between natural fractures and matrix is calculated by dual-media model, while flux transfer between hydraulic fractures and matrix/natural fractures is calculated based on the EDFM method.

Numerical Issues
2.3.1. Discretization Based on Finite Volume Method. In this work, the finite volume method is adopted to discretize the governing equations discussed in Section 2.2. Spatially, the reservoir is discretized into regular Cartesian grids (control volume), and Figure 1 shows the two-dimensional case. Particularly under this condition, the line connecting grid i and grid j is perpendicular to the public shared interface between two grids, and by utilizing Gauss' theorem, the volume integral of the divergence operator can be converted into the 3 Lithosphere surface integral of certain variables surrounding control volumes.
Based on the finite volume method, the governing equations of the gas phase can be discretized into a numerical scheme of pseudopressure, including the transmissibility calculation, the flux transfer equation in dual-media model, and the well flow equations. Because the compressibility of water is comparably small than that of the gas phase, there is no need to carry out functional transformation or introducing "pseudopressure" for the governing equations of water. After the discretization based on the finite volume method for this two-phase nonlinear system, the Newton-Raphson iteration technique is then utilized to solve the problem.
As mentioned before, in our model, we solve "pseudopressure" for the gas phase instead of real gas pressure. Particularly, the pseudopressure of the gas phase can be expressed as Note that the gas pressure is the primitive variable in our framework, and the pressure of the water phase can be obtained by subtracting the gas pressure by capillary pressure. By considering the gravity effect of the gas phase, the above equation can be also expressed through the potential function: From equations (7) and (8), it is implied that for the gradient of the gas pressure and pseudopressure, we have Also, for the gradient of the potential function, we have Based on the TPFA scheme, the transmissibility can be calculated by Note that the gas pressure and gravity term can be combined and expressed by the pseudopressure; then, a simpler formula for transmissibility calculation reads as For temporal discretization, the backward implicit Euler scheme is adopted. As a summary, the discretized equations can be expressed as follows.
Gas equation in fractures: Water equation in fractures: Note that there exist two types of flow in the matrix. One is the flow in the matrix outside of SRV regions (no fractures), and another one is the flow in SRV regions which is described by the dual-media model.
Gas equation in matrix: Water equation in matrix: Gas flux transfer rate between matrix and fractures can be calculated by  Then, the flux transfer between fractures and the matrix can be written as a boundary integral form: where ðξ, ηÞ indicates inner points in grids, ðx, yÞ indicates boundary points of grids, ∂Ω is the grid boundary illustrated as red lines in Figure 2, n is the unit normal vector of grid boundary, sðx, yÞ is the unit direction vector of grid boundary, and λðξ, ηÞ is the coefficient of the to-be-solved point. Namely, when the point locates at the grid boundary, the coefficient is 0.5. When the point locates outside the grid, the coefficient is 0. When the point locates inside the grid, the coefficient is 1. n f is the fracture element number in the current grid, and ∂Γ is the fracture geometry.
The function G L ðx, y ; ξ, ηÞ is the basic solution of the 2D Laplacian equation: In our work, fracture elements are considered as constant elements. The pressure of each element is represented by the value at the element centroid (p f 1,2,3,4 in Figure 2), and pressure on the grid boundary (p ed1,2,3,4 ) is obtained through interpolation with the current grid and neighbor grids. For example, p ed1 can be calculated by where k x0 and k x1 are permeability values of two neighbor grids in the X direction and dx 0 and dx 1 are grid sizes in the X direction. Then, the flux transfer of the oil phase between the matrix and fractures is The equations for water and gas are presented in the Appendix.
The flow inside fractures is considered as a onedimensional flow, of which discretization is similar with  5 Lithosphere the flow in the matrix. It is worth noting that we adopt a star-like transformation to calculate the conductivity terms for fracture grids, which is illustrated in Figure 3.
Then, the conductivity of arbitrary neighborhood fracture grids can be expressed as where where w f i is the width of the i-th fracture grid and l f i is the length of the i-th fracture grid.

Workflow of Construction of Mixed Fracture Model.
As mentioned in previous sections, the embedded discrete fractures can be considered as arbitrary polygon plates which are embedded into regular matrix grids. Usually, structured grids are utilized to represent the matrix, such as Cartesian grids and corner-point grids. With these manipulations, the EDFM is capable of simulating complex fracture networks and maintaining geometric details of fracture distribution. Moreover, there is no need to carry out grid remeshing for matrix grids during simulation, which is more convenient for the geometric searching process. Based on the above discussions, the workflow of the mixture fracture model can be summarized as follows: (1) Carry out modeling of natural fractures in shale reservoirs and determine key parameters such as fracture distribution and density (2) Obtain the basic outlines of artificial fractures based on microseismic data or fracturing simulation data (3) Represent artificial fractures as embedded discrete fractures and calculate the conductivity between polygon elements and matrix grids as shown in Figure 4 (4) Build up a conductivity relationship between wellbores and artificial/natural fractures In the implementation of EDFM, the EDFM-related searching manipulation, discretization, and calculation of conductivity are all assembled into the reservoir simulator. As illustrated in Figure 5, first, fracture polygons are generated through conventional approaches, grids are then refined, and finally, connections between fractures and matrix are then appended into the workflow.

Field Application of Deep Shale Gas
Simulation and Analysis 4.1. Basic Information of Investigated Work Area. The investigated work area locates at the first segment of the Wufeng-Longma play in Southern Sichuan, of which the formation thickness is about 38 meters. The TOC of shale reservoir is larger than 2%, and the classification of organic matter is dominated by types I and II1. The porosity of this reservoir is in the range from 2.5% to 6.4%. The content of clay accounts for 37.5%, and the brittle matter accounts for 61~65%. Based on the geological interpretation, top and bottom formation surfaces, including the overall geological body, are constructed as totally 9 layers. The total area of investigated work area reaches 13.7 km 2 , and the total thickness in the vertical direction reaches 84 meters. Figure 6, in the geological model, the average thickness of the first 8 sublayers is 7.5 m, 3.  3.94%, and 3.85%. By utilizing the sequential Gaussian interpolation method, a spatially stochastic field describing porosity distribution is obtained which obeys Gaussian dis-tribution as illustrated in Figure 7. Particularly for this realization, the average porosity of the first 3 layers is 5.12%, and the average porosity of the first 5 layers is 5.1%.   The initial average water saturation of the shale reservoir is set as 39.48%, 20.97%, 40.16%, 49.01%, 46.55%, 48.39%, 50.93%, 51.15%, and 55.75%, for 9 sublayers. Similar with that of the porosity field, the sequential Gaussian method is again utilized to generate an initial water saturation field obeying Gaussian distribution, as illustrated in Figure 8   Based on the BHP detection, curvature distribution, and dynamic analysis, the connection map between wells through discrete fractures is then obtained as illustrated in Figures 11 and 12.

Characteristics of Fluid in Fractures.
After fracturing, the pump pressure is about 63~69 MPa, which implies that the BHP is about 100 MPa. Namely, the fluid pressure in fractures after which pumps stop working is about 100 MPa, and the fractures are full of water with 100% saturation.

History Matching.
History matching is mainly accomplished through adjusting the conductivity of main fractures, fracture permeability of SRV regions, and shape factors of the SRV region. During the flowback stage, the BHP is assimilated given that water production rate is assumed to be constant. And during the production stage, the water production rate and BHP are assimilated assuming that gas production rate is constant. Because the production history of the platform is very complicated and well shutdown is very frequent, we carry out history matching by selecting the shutdown wells and assimilating one-week-period production data. Note that the history data during the production stage is obtained through the averaging method, while the BHP of shutdown wells is directly selected from the raw data.
The history matching curves are then obtained as illustrated in Figures 14-17. Figures 14 and 15 show the gas production rate and wellbore hole pressure of well 1-3HF. Figures 16 and 17 illustrate the BHP of well 1-1HF of two stages, i.e., before and after simulation-restart, respectively. Note that history matching curves are represented by solid lines, while actual observed production data are plotted as discrete dots.
Note that wells 1-1HF and 1-4HF are both utilized for production after several years of production of well 1HF. Figures 18 and 19 show the gas production rate and fluid production rate of well 1-1HF, respectively. Figures 20 and  21 show the gas production rate and fluid production rate of well 1-2HF, respectively. Figures 22 and 23 show the wellbore hole pressure of well 1-1 and well 1-2HF, respectively.

Analysis of Exploitation Status and Index Prediction of
Gas Wells. The well 1HF was first put into fracturing and            Figure 19: Fluid production rate (m 3 /day) of well 1-1HF.     11 Lithosphere appended into simulation. Finally, the five wells are put into production simultaneously. Most of the well tracelines go through sublayers 1~3. As illustrated in Figure 24, the traceline of well 1-3HF touches the top region of formation. Expect for well 1-3HF, most fractures locate at sublayers 1~4. Figure 25 illustrates the pressure distribution of sublayers 1~5 at 2020. Here, we summarize the exploitation status of each sublayer as follows, given a threshold value of pressure drop (10 MPa). The exploitation areas of sublayers 1~5 are 1.03 km 2 , 1.03 km 2 , 1.40 km 2 , 1.27 km 2 , and 0.76 km 2 , respectively. For a single well, the exploitation length is at the range from 80 m to 320 m, and the unexploitation distance between wells is at the range from 75 m to 280 m.
Through history matching, the distribution of fracture networks and the effective permeability of the SRV region can be obtained. Considering the case under conditions that the production rate is 55,000 m 3 /d, the minimum BHP is 7.5 MPa and 2.5 MPa, and total production time is 30 years; the EUR of the platform is then estimated as listed in Table 1.

Summary and Conclusions
In this paper, features of complex fracture networks and the gas-water two-phase flow are investigated during the exploitation process of deep shale gas reservoirs through horizontal wells and hydraulic fracturing techniques. We propose a simulation framework for describing gas-water twophase flow in deep shale gas reservoirs based on the mixed fracture method. The main conclusions are drawn as follows: (1) By considering the special gas-water two-phase flow mechanisms in shale gas reservoirs, a multiscale model is constructed for deep shale reservoirs by defining the pseudopressure function and derivation of governing equations, which suppresses the nonlinear characteristics of gas flow and improves computational efficiency for numerical solution (2) By combining multiple types of software such as Petrel, and a shale gas simulator, a standard investigation process is built up, i.e., prefracturing modeling, postfracturing modeling, dynamic analysis, and numerical simulation, which provides solid support for multiwell simulations (3) Through numerical simulation and history matching, it is found that under the condition of 20 m fracture height, the planar distribution range of the assimilated fracture can reach over 100 m. The reason why this value is smaller than the estimation from the microseismic data is that the assimilated fracture length actually reflects the "effective" length for gas-water flows (4) The united well-platform simulation is capable of estimating the exploitation status of whole formation and simulating the effects of random distributed fractures. Besides, the interference of fractures between wells may increase the EUR value for a particular period, while the total EUR is positively proportional to the extent of reservoir stimulation    (5) It is necessary to consider the effects of fractures and stress sensitivity in the matrix when considering the dynamic change of production during the flowback and production stages. Because the fluid pressure in fractures tends to decrease after flowback and gas production begins under the high stress in deep shale reservoirs, fractures are gradually closed. During the early stage of shale gas production, it appears that the production curve is typically dominated by linear flow, and the gas-oil ratio stays at a constant value. Therefore, for the prediction of production, it is an efficient approach to first use models for predicting the gas production rate and then predicting the oil production rate through the gas-oil ratio (GOR).
(6) The limitation of the study is that the transient flow during the early stage is effected by the grid size near the fracture, and a small grid size is recommended. However, large numbers of the grid will slow down the computation time; the LGR grid with the EDFM model will be investigated in the succeeding study