A Hybrid Model for Simulating Fracturing Fluid Flowback in Tight Sandstone Gas Wells considering a Three-Dimensional Discrete Fracture

Two-phase (gas + water) flow is quite common in tight sandstone gas reservoirs during flowback and early-time production periods. However, many analytical models are restricted to single-phase flow problems and three-dimensional fracture characteristics are seldom considered. Numerical simulations are good choices for this problem, but it is time consuming in gridding and simulating. This paper presents a comprehensive hybrid model to characterize two-phase flow behaviour and predict the production performance of a fractured tight gas well with a three-dimensional discrete fracture. In this approach, the hydraulic fracture is discretized into several panels and the transient flow equation is solved by the finite difference method numerically. A three-dimensional volumetric source function and superposition principle are deployed to capture the flow behaviour in the reservoir analytically. The transient responses are obtained by coupling the flow in the reservoir and threedimensional discrete fracture dynamically. The accuracy and practicability of the proposed model are validated by the numerical simulation result. The results indicate that the proposed model is highly efficient and precise in simulating the gas/water two-phase flow and evaluating the early-time production performance of fractured tight sandstone gas wells considering a three-dimensional discrete fracture. The results also show that the gas production rate will be overestimated without considering the two-phase flow in the hydraulic fracture. In addition, the influences of fracture permeability, fracture half-length, and matrix permeability on production performance are significant. The gas production rate will be higher with larger fracture permeability at the early production period, but the production curves will merge after fracturing fluid flows back. A larger fracture half-length and matrix permeability can enhance the gas production rate.


Introduction
Hydraulic fracturing technology can achieve the extraction of hydrocarbon from unconventional reservoirs at a low cost and improve the well performance significantly, which has been attracting considerable attention in recent years [1][2][3][4][5][6][7]. To generate highly conductive hydraulic fractures to increase well productivity, almost thousands of fracturing fluids are pumped into the reservoir [8]. The two-phase (gas + water) flow in the hydraulic fracture is quite common during the flowback and early-time production periods [9][10][11][12]. A common challenge for tight sandstone gas reser-voir evaluation is the characterizations of the threedimensional fracture. The productivity prediction of tight gas wells is greatly affected by the characteristics of the gas-water two-phase flow and fracture seepage parameters. Much research has been dedicated to theoretical models of fractured tight gas wells [13][14][15][16], but most are only applicable to single-phase flow problems. Moreover, the threedimensional fracture characteristics is seldom considered in these models.
Recently, flow behaviour analyses of flowback and earlytime production from hydraulically fractured tight gas wells have been well known and researched. There is an array of analytical, semianalytical, and numerical simulation forecasting methods that can be used to analyze flowback data and make forecast for tight gas wells. Analytical models are based on the assumption that hydrocarbon does not breakthrough from fractured tight gas wells [17][18][19][20]. However, these methods are limited to the single-phase flow of fracturing fluid in the hydraulic fracture and are generally inappropriate for analyzing flowback data after hydrocarbon breakthrough. Numerical simulation can rigorously consider complex flow physics, such as multiphase flow, reservoir, PVT properties, and fracture characteristics. It has gained much attention and been employed to investigate the twophase flow behaviour and production performance of fractured tight gas wells during flowback [21][22][23][24]. However, time-intensive computation performance and the difficulty of gridding make this method less attractive. In recent years, several semianalytical models have been developed to simulate two-phase flowback characteristics in tight oil and gas reservoirs [25][26][27]. This approach can capture the twophase flow of gas and water during the entire flowback period and have much higher computational efficiency than numerical simulation methods, which has attracted significant attention. In their model, the details of the transient two-phase flow in the fracture are incorporated. However, to simplify the mathematical model, the fluid flow in fractures is assumed as a 1D flow and the flux variation along the vertical direction in the fracture is ignored. Therefore, the semianalytical models cannot accurately capture the three-dimensional fracture characteristics, resulting in significant errors in production performance analysis.
This work aims at providing a hybrid model for simulating fracturing fluid flowback in tight sandstone gas wells considering a three-dimensional discrete fracture. We discretize the fracture into several panels, and each panel is treated as a plane source. The orthogonal transformation method and Green's function method are employed to obtain the solution in the matrix system analytically. The finite-difference method is used to characterize the gas/water two-phase flow in the hydraulic fracture numerically. The solution of the model is semianalytically solved by coupling the transient flow in the fracture and matrix with the continuity of pressure and flux on the fracture surfaces. Based on the proposed model, the precision and efficiency of the proposed model are validated by using the commercial numerical simulator Eclipse. In addition, the influence of the two-phase flow in the hydraulic fracture on transient rate responses during flowback is analyzed. After that, the effects of several critical parameters on production performance are analyzed. The novelty of the proposed model is in the ability to semianalytically model the fracturing fluid flowback of fractured wells with three-dimensionally discrete fracture. The main advantage of the hybrid solution over more conventional numerical simulation is the reduced model setup and runtime without the loss of essential physics.

Methodology
2.1. Model Assumption. The scheme of a fractured tight sandstone gas well is shown in Figure 1(a). Figure 1(b) is the plan view of the well in the xy plane. The reservoir is box shaped, and the external boundaries of the top, bottom, and side are all assumed to be closed. The fracture is vertical, in the rectangle shape, and fully penetrate the formation. The fracture is initially filled with fracturing fluid after fracture stimulation and contributes to the fracture fluid recovery during flowback. The fluid flow in the fracture is considered as a 2D flow, and the gas enters the wellbore only through the fracture. This work assumes that water is only mobile in the fracture system to simplify the model. This physical problem is similar to what was assumed by Clarkson and Williams-Kovacs [11] and Jia et al. [28] to describe two-phase flowback behaviour in unconventional reservoirs. Some simplifying physical model assumptions for the derivation of the governing equation are listed as (1) the reservoir is assumed to be isotropic, (2) Darcy's flow is considered both in the fracture and formation, and (3) the gas and water are produced with a constant bottom hole pressure.

Mathematical Model.
As illustrated in the physical model, the hydraulic fracture system is the flow pathway directly connected with the wellbore. The fluids in the fracture system first flow into the wellbore, followed by the matrix-fracture flow. The fluid flow from the reservoir through the fracture and into the wellbore is composed of the two systems. This section will describe the mathematical model for the fluid flow in the hydraulic fracture system and in the matrix system.

Fluid Flow in the Hydraulic
Fracture. Compared with the length and height of the fracture, the width of the fracture is minimal. So, the fluid flow in the three-dimensional fracture can be simplified as a two-dimensional flow and the mass conservation equations for the gas and water phases are given as where u ξ and u η are the darcy flow rates, ξ and η represent the direction of the fracture, B g and B w are gas and water volume coefficients, respectively, ρ g and ρ w are gas and water densities, respectively, S g and S w are gas and water saturations, respectively, ϕ F is fracture porosity, and t is time. Based on the Darcy law, the flow equation in the hydraulic fracture is written as where μ g and μ w are gas and water viscosities, respectively, k rg and k rw are relative permeabilities of the gas and water phases, respectively, p is the reservoir pressure, and β is the transmissibility conversion factor, where β = 5:6146 for the field unit system. The saturation relation equation is expressed by Substituting equations (3) through (6) into (1) and (2), the governing equations in the three-dimensional fracture for the gas and water phases can be obtained.
For the gas phase, the governing equation can be described as For the water phase, For the initial condition, The well is assumed to produce at constant pressure, and the well production equation can be written as follows [29]: r eq = 0:14 where k F is fracture permeability, p F is fracture pressure, p i is initial reservoir pressure, ẽ q SCFg is a source term representing the gas flux entering the fracture from the matrix at a point on the fracture surface per unit volume at standard conditions, ẽ q SCWg and ẽ q SCWw represent the gas and water production rate per unit volume at standard conditions, respectively, r W is the wellbore radius, w F is the width of the fracture, Δl F is the length of the fracture panel, and Δ h F is the height of the fracture panel.

Fluid Flow in the Tight Sandstone Gas
Reservoir. Based on the model assumption, for an isotropic-permeability reservoir, the governing differential equation for a single gas phase can be expressed as The initial condition is given by The external boundaries, including the top, bottom, and side, are assumed to be closed and given by for the side boundary, for the side boundary, 3 Lithosphere for the top boundary, and for the bottom boundary, where k M is the matrix permeability, p M is the matrix pressure, ϕ M is the matrix porosity, c tM is the total compressibility of the matrix system, x e , y e , and z e are dimensions of the reservoir, (x w , y w , z w ) is the center position of the plane source, andq is a point source influx from the matrix to the fracture.

Solution to the Mathematical Model
2.3.1. Solution of the Numerical Fracture Flow Model. In our work, the finite difference method is used to obtain the solution in the fracture system. The fracture is first discretized into small panels, as shown in Figure 2. Then, we apply the finite difference approximation to equations (7) and (8) at the time level n + 1 and the implicit finite-difference form of fracture panel ði, jÞ can be obtained. For the gas phase, For the water phase, where i = 1 ⋯ N x and j = 1 ⋯ N y , N F is the total fracture panels and N F = N x ⋅ N y , ψ i,j is a 2D array whose elements are the existing neighboring segments of fracture panel i, j, and T n+1 g l,ði,jÞ and T n+1 w l,ði,jÞ are the transmissibility of the gas and water phases between fracture panel ði, jÞ and its adjacent panels at time level n + 1 and given by: The geometric transmissibility between fracture panel ði, jÞ and its adjacent panels is given by where γ i,j is the geometric transmissibility between fracture panel ði, jÞ and the interface, which is expressed as The terms (C n+1 (18) and (19) are defined as where V b i,j is the bulk volume of fracture panel ði, jÞ, ww i, j are similar to the compressibility coefficient, but are functions of pressure and saturation, and the coefficients of ð1/B g i,j Þ′, ð1/B w i,j Þ′, and ϕ F i,j ′ are given by Figure 2: Discretization of the hydraulic fracture along both the horizontal axis and the vertical axis.
where B°w is the water volume coefficient at a reference pressure, ϕ°F is fracture porosity at a reference pressure, and c w and c F are the water compressibility and fracture compressibility, respectively. Applying the finite difference approximation to equation (10) at the time level n + 1, the implicit finite-difference form of production boundary condition can be expressed as where p n W is the wellbore pressure. The transmissibility of the gas and water phases between fracture panel l and the wellbore at the time level n + 1 is given by T n+1 Combing equation (30) with (18), the gas flow equation can be rewritten as Substituting equation (31) into (19), the water flow equation can be rearranged as: Equations (24) and (25) are the implicit finite-difference form of any fracture panels ði, jÞ at the time level n + 1. Applying the two equations to N F fracture panels, the fracture flow equations can be arranged into a matrix format yield where the matrix and vector in equation (34) are given by And the pseudotime is defined as where m is the pseudopressure and Z is the gas compressibility factor. To make the equations homogeneous, some dimensionless variables are defined and tabulated in Table 1. Then, using the dimensionless variables and applying the pseudopressure transformation to equations (12)-(17) results in for the initial condition, for the side boundary, for the side boundary, for the top boundary, and for the bottom boundary, where l R is the reference length, l R = Δl F /2, m i is the initial pseudopressure of the reservoir, m M is the pseudopressure of the matrix system, and T M is the temperature. In this section, the orthogonal transformation method is used to derive the analytical solution of the reservoir flow model. Applying the orthogonal transformation method to equations (39)-(44), one can obtain that Equation (45) is the initial value problem of ordinary differential equations, and the solution of the equation can be expressed bỹ where φ is a three-dimensional eigenfunction and λ is its eigenvalue. The 3D eigenfunction and its corresponding eigenvalue can be obtained by solving the three 1D eigenfunctions and eigenvalues in the x, y, and z directions, respectively.
The one-dimensional eigenfunctions and eigenvalues of equations (47)-(49) are given by Thus, the three-dimensional eigenfunction and the corresponding eigenvalue can be expressed as Dimensionless distance in the x direction x D = x/l R Dimensionless distance in the y direction y D = y/l R Dimensionless distance in the z direction z D = z/l R Dimensionless outer boundary in the x direction x eD = x e /l R Dimensionless outer boundary in the y direction y eD = y e /l R Dimensionless outer boundary in the z direction After that, utilizing the inverse transformation of equation (46), the pressure solution for a three-dimensional volumetric source can be obtained and given bỹ Combing the definitions of the parameters in Table 1, we can rewrite equation (52) as where m M is the current pseudopressure of the matrix system, m _ M is the initial pseudopressure of the matrix system, q sc is the withdrawal rate at standard condition, χ _ is the diffusivity and evaluated at the initial pressure, and α is the unit conversion factor. The fracture is discretized into several fracture panels, and each fracture panel can be regarded as a plane source, as shown in Figure 3. Based on the theory of point-source function, integrating equation (53) along with the directions of the fracture panel results in where x i,j , y i,j , z i,j are the center position of the fracture panel ði, jÞ and equation (54) is the pressure response at an arbitrary position ðx, y, zÞ at time t due to the contribution of a continuous plane source with a withdrawal rate of q sc . According to the superposition principle, for the system in Figure 3, the pressure drop at the point (x o , y o , z o ) caused by all fracture panels can be calculated and the equation is given by whereq k SCM i,j is the gas flow rate from the matrix to the fracture panel ði, jÞ per unit fracture area at time step k and SðÞ is the plane source function and expressed by To simplify the form of equation (55), let H n+1,k o,ði,jÞ denote Plane source

Substituting equation (57) into equation (55) yields
With the superposition principle, the pressure response of fracture panel ðI, JÞ caused by N F fracture panels at time level n + 1 can be calculated by Equation (59) can be arranged into a matrix format and yields In equation (60), the matrix H and vector h ! are defined as

Coupling Conditions for the Flow in the Two Systems.
The flow in the matrix system and the discrete fracture system should be coupled to obtain the solution. Based on the flux exchange of the two systems and pressure continuity conditions on the fracture surface, the gas pressure continuity condition can be expressed as And the gas flux continuity condition can be written as We can group the fracture equations and the matrix equations to develop a system of equations to characterize the two-phase flow performance of the three-dimensional fracture.
There are 3N F unknowns and 3N F independent equations in equation (65), and a unique solution is sure to be obtained. One should notice that the greatest challenge in solving equations (35) and (60) is to linearize equation (35). This is because the terms in matrix A and vector B are functions of pressure and saturation of the fracture panel, which are unknown parameters. This study uses the upstream weighting method for linearization in space and the simple iteration method for linearization in time. Then, the equations can be solved by a computer program. In the following, we analyzed the results based on this solution.

Results and Discussion
In this section, the proposed model is firstly validated using the numerical simulator Eclipse. Then, the effects of some critical fracture and formation parameters on rate transient behaviour during both the flowback period and the earlytime production period are investigated.

Model Verification.
To validate the model proposed in this paper, a fractured tight gas well model is built using the numerical simulator Eclipse and the schematic is shown in Figure 4. For comparison, the formation and fracture parameters are the same for both the Eclipse and the proposed model, which are listed in Table 2. The dimension of the whole reservoir is 600 × 800 × 10 m, and the grid sizes for x, y, and z directions are 20 m, 20 m, and 10 m, respectively. The dimension of the fracture is 200 × 15 × 0:01 m, and the grids near the hydraulic fracture face are gradually refined to capture the flow behaviour accurately. The fracture grids are initialized with two phases (gas + water), while the matrix is assigned 100% gas saturation. The gas PVT properties and the gas-water relative permeability curves in the hydraulic fracture system are shown in Figure 5.
The validation results are shown in Figure 6; it could be found that the result of the proposed model is in accordance with the numerical simulation result. The matching error is less than 2%. This shows the ability and precision of the proposed method in handling fracturing fluid flowback in tight sandstone gas wells considering three-dimensional discrete fracture. In addition, we also compared the computational speed of the two models. For the above simulations, the computation time for the numerical model is 48.2 s, while the computation time for the proposed model is 12.1 s. It is concluded that the proposed model is quite efficient in simulating the two-phase flow behaviour in three-dimensional discrete fractures.

Effect of the Two-Phase Flow in the Hydraulic Fracture.
The main difference between the proposed model and the traditional models is that the three-dimensional fracture characteristics and two-phase flow behaviour are incorporated in the proposed model. In this section, we simultaneously simulated the rate transients of both the proposed and conventional models using the same group of formation and well parameters in Table 2. The initial water saturation in the fracture is given 0.9 and 0 in the proposed and conventional models, respectively. The modelling results are shown in Figure 7. It is shown that two-phase flow in the fracture system significantly influences the gas production rate at the early time. However, there is little influence after 60 days of production. In addition, the early gas production rate is smaller if the water-phase flow is considered in the fracture system. This is possible because the gas production from the tight matrix acts as a displacement effect for the limited water in the fracture system, so the two-phase flow in the fracture system only affects the gas production within the early flowback period. Besides, the gas production rate  9 Lithosphere will be overestimated without considering the two-phase flow in the fracture.

Sensitivity Analysis.
For tight sandstone gas reservoirs, the fracturing parameters (fracture permeability and fracture half-length) and even matrix permeability play a significant role in well performance and long-term production forecasting. Based on the proposed model in this paper, the influences of fracture permeability, fracture half-length, and matrix permeability on rate transient behaviour during the flowback period are analyzed. Except for the parameters analyzed, other parameters are the same and are shown in Table 2. 3.3.1. The Effect of the Fracture Permeability. The effect of the fracture permeability is analyzed in this part. Reference fracture permeability is given as 100, 300, 500, and 1000 mD, and the modelling results are shown in Figure 8.
It should be noticed that gas and water exhibit higher rates with larger fracture permeability at the early production period. It is also observed that the fracture permeability significantly affects the gas rate during the early flowback period, while the curves will merge after 50 days of flowback. It is possible because the gas and water flow resistance in the fracture reduces with larger fracture permeability. In addition, the larger fracture permeability accelerates the water depletion in the fracture and the system will be dominated by the matrix earlier. So, the gas production rate becomes larger and the peak appears earlier as the fracture permeability increases.

The
Effect of the Half-Length of the Hydraulic Fracture. Figure 9 illustrates the effect of the half-length of the hydraulic fracture on gas and water production rates. The halflength of the hydraulic fracture is given 80, 100, 120, and 150 m for the case study in this section. We could find that the effect of the length of the hydraulic fracture on the gas and water production rates is significant throughout most of the flowback period. As the length of hydraulic fracture increases, water exhibits higher production rates early but declines faster later. Besides, the gas production rate increases rapidly and the production peak appears later during the early flowback period with a larger half-length. This is possible because the half-length of hydraulic fracture reflects the relative capacity of fracturing fluid stored in the fracture system; a bigger half-length is the response of relative abundant reserves in the fracture system. For the fracture with a large length, the fracturing fluid at the far end generates strong support for water saturation near the wellbore. That is the reason why the gas production peak appears later as the fracture half-length increases.
3.3.3. The Effect of the Matrix Permeability. Figure 10 shows the effect of the matrix permeability on gas and water production rates. Herein, the matrix permeability is given 0.01, 0.05, and 0.1 mD for the case study and noted that matrix  Water rate-proposed model Water rate-numerical simulation

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

Conflicts of Interest
The authors declare that they have no conflicts of interest.