#### Abstract

After multistage fracturing, the flowback of fracturing fluid will cause two-phase flow through hydraulic fractures in shale gas reservoirs. With the consideration of two-phase flow and desorbed gas transient diffusion in shale gas reservoirs, a two-phase transient flow model of multistage fractured horizontal well in shale gas reservoirs was created. Accurate solution to this flow model is obtained by the use of source function theory, Laplace transform, three-dimensional eigenvalue method, and orthogonal transformation. According to the model’s solution, the bilogarithmic type curves of the two-phase model are illustrated, and the production decline performance under the effects of hydraulic fractures and shale gas reservoir properties are discussed. The result obtained in this paper has important significance to understand pressure response characteristics and production decline law of two-phase flow in shale gas reservoirs. Moreover, it provides the theoretical basis for exploiting this reservoir efficiently.

#### 1. Introduction

In order to obtain industrial gas flow, horizontal drilling and multistage hydraulic fracturing techniques have been widely used in SGR (shale gas reservoirs) development. Compared to the conventional reservoirs, SGR have their own uniqueness; for example, the shale is a self-sourcing reservoir in which the gas is stored either by compression or by adsorption on the surface of the solid material. Transport of gas in SGR is a complex multiscale flow, which is from macroscale to the molecular scale. Numerous researchers have studied migration mechanisms in SGR, but the two-phase flow caused by multistage hydraulic fracturing is usually ignored for analyzing and forecasting a shale gas well.

According to Hill and Nelson [1], 20% to 85% of the total gas storage in SGR may be subject to sorption. Montgomery et al. [2] pointed out that the absorbed gas in SGR is mainly stored in pores, natural fractures, and matrix surface and desorption rate would change with the difference of desorption time, desorption pressure, and organic matter volume. Freeman et al. [3] and Cheng [4] employed a numerical simulator to study the flow regimes for a multistage hydraulic fractured horizontal well in SGR and the effect of gas desorption was briefly discussed. The studies of hydrofracture are mainly concentrated on the fractures’ characteristics and fluid flow mechanism. Xiao et al. [5] presented an analytical model for effective thermal conductivity of nanofluids by considering the effect of Brownian motion of nanoparticles. Tan et al. [6] presented a dual fractal reservoir transient flow model created by embedding a fracture system simulated by a tree-shaped fractal network into a matrix system simulated by fractal porous media and discussed the influence of different fractal factors on pressure transient responses. Cai et al. [7–9] analyzed the capillary-driven flow in gas-saturated porous media based on fractal theory and modified the classical Lucas-Washburn equation. Later, the hydraulic conductivity and effective porosity in porous media have been investigated by many researchers, such as Tan et al. [10]. Freeman et al. [11] presented a complicated numerical model to study the adsorption and desorption in SGR and described the change of gas composition during wells production. After that, some special fractal reservoir models have been presented by Tan et al. [12, 13]; in these papers they discussed dual fractal flow characteristic and liquid carrying phenomenon.

Some scholars employed new methods to analyze production decline of SGR. Medeiros et al. [14] discussed the analysis of production data from hydraulically fractured horizontal wells in shale reservoirs. Ozkan et al. [15] presented a discussion of fractured horizontal-well performance in conventional and unconventional reservoirs and provided interpretation of the different production performance of fracturing horizontal wells in both types of formations. After that, scholars presented rate prediction models [16–18] in SGR, but the proposed models may have some limitations when applied to analyze production data in a shale well. Cho et al. [19] presented the results of an experimental study of pressure-dependent natural-fracture permeability in SGR and found that the effect of pressure-dependent natural-fracture permeability on shale gas well production was a function of the permeability of the matrix system. J. H. Lee and K. S. Lee [20] studied the multiscale heterogeneity in the existence of various permeability magnitudes and spatial distributions of permeability. Du et al. [21] presented a new experimental method explaining the flooding tests in the low-permeability reservoirs. Alkouh and Wattenbarger [22] presented a gas-water model to solve the challenge of simulating a natural fracture with trapped water without imbibitions in SGR. In the same year, Alkouh et al. [23] added post-frac water flowback and long-term water production to the analysis and showed the benefits of a new method for combining water flowback and long-term water production data in shale gas analysis. Adefidipe et al. [24] presented a new model to describe the two-phase flow in shale gas wells and analyzed the obvious gas phase coning phenomenon. Xie et al. [25] presented a transient flow model which considered two-phase flow after hydrofracture in shale gas reservoirs to analyze the transient pressure response.

Until now, the studies of a multistage hydraulic fractured horizontal in SGR are mainly concentrated on single phase gas flow; the two-phase flow caused by flowback after hydrofracture in SGR does not attract much attention. In order to understand the effect of two-phase flow in SGR, a two-phase flow transient analysis model of multistage fractured horizontal well in SGR is presented, and the accurate solution to the flow model is obtained using source function theory [26, 27], Laplace transform, three-dimensional eigenvalue method, and orthogonal transformation. Type curves of pressure response and production decline are plotted and different flow regimes are identified. The influence of different factors on production decline performance is discussed.

#### 2. Physical Model

Figure 1 shows a horizontal well with hydraulic fractures in the SGR center. The physical model assumptions are as follows.(1)The SGR are dual mechanism and with sealed boundary.(2)Two-phase flow exists in fractures and horizontal well.(3)All fractures are assumed to be equally spaced, to penetrate the formation completely, and to be perpendicular to the horizontal wellbore. Both the horizontal wellbore and the fractures are assumed to be infinitely conductive.(4)Ignoring the impact of gravity and capillary forces.

#### 3. Mathematical Model

##### 3.1. The Two-Phase Flow Governing Equations in Fracture System

Governing flow equations of two-phase flow in fractures under three-dimensional Cartesian coordinates can be written as follows.

Gas flow equation in fractures is

Liquid flow equation in fractures is

The transient diffusion rate of the shale gas from the shale matrix can be expressed as follows:

The two-phase flow pseudopressure is employed and defined as follows:

With (1), (2), and (4), the two-phase flow governing equation can be written as follows:

Employing Taylor expansion to simplify (5), we can obtain the following equation:where

Employing source function into (6) yields

To simplify the mathematical model, dimensionless parameters are defined as follows.

The dimensionless two-phase pseudopressure is

The dimensionless wellbore storage coefficient is

The dimensionless time is

The dimensionless coordinate value of is

The dimensionless coordinate value of is

The dimensionless coordinate value of is

The dimensionless relative mass concentration is

The dimensionless radial distance in matrix system is

Using the dimensionless parameters, (8) can be written as follows:

The two-phase flow mathematical model in Laplace space is obtained by taking the Laplace transformation [28] of (17). The two-phase flow mathematical model is as follows:

With the condition of sealed boundary, the dimensionless form of the initial condition is

The dimensionless forms of the outer boundary conditions are

##### 3.2. The Transient Diffusion Rate Equations in Matrix System

With Fick’s second law and the assumption of spherical matrix blocks, the transient diffusion rate equation can be written as follows [29]:

Using dimensionless parameters and Laplace transform to rewrite (21),

For transient diffusion case, the center of the spherical matrix block can be treated as a no-flow boundary, so the inner boundary condition can be written as follows:

The concentration on the external surface of matrix block is in equilibrium with the pressure in the fracture. The outer boundary condition can be written as follows:

The general solution of (22) can be expressed as follows:

With the boundary condition and general solution expression, the accurate solution of the dimensionless mathematical model of matrix system for transient diffusion case can be obtained as follows:

With the definition of hyperbolic functions, (26) can be rewritten aswhere is the dimensionless equilibrium volumetric gas concentration. Employing Langmuir isotherm [30], there is

Equation (28) can be rewritten as

Inserting (29) into (27) and getting derivative with respect to , the solution of the matrix model in the Laplace domain is as follows:

##### 3.3. The Two-Phase Flow Mathematical Model for Multistage Fractured Horizontal Well in SGR

Inserting (31) into (18) results in

To simplify, (32) can be rewritten as

#### 4. Mathematical Model Solution

##### 4.1. Eigenvalue Method and Orthogonal Transformation

Employing three-dimensional eigenvalue method and orthogonal transformation to solve the accurate solution of the mathematical model,where is the expression of dimensionless two-phase pseudopressure in three-dimensional eigenvalue space, and we can transform the three-dimensional problem to three one-dimensional problems:

With (36), eigenvalue in direction can be obtained as , and eigenfunction in direction is , .

In the same way, direction’s eigenvalue can be obtained as , and eigenfunction in direction is , .

direction’s eigenvalue can be obtained as , and eigenfunction in direction is , .

Use the one-dimensional eigenvalues and eigenfunctions to obtain three-dimensional eigenvalue and eigenfunction as follows:

The eigenfunction in (38) is a complete orthonormal system in three-dimensional space; the completeness orthogonality equation is as follows:

Take the eigenvalues and eigenfunctions into (38) to rewrite the completeness orthogonality equation as follows:where

Employing to define orthogonal transformation, then (33) can be solved as follows:

With the definition of orthogonal inverse transformation and norm, there are equations as follows.

Orthogonal inverse transformation is as follows:

Norm is as follows:

The value of norm can be obtained as follows:

Taking (43) and (44) to transform (42), the point source expression can obtained as follows:

##### 4.2. Model Solution

Based on the assumptions, all the fractures fully penetrated the SGR and are infinitely conductive. According to the Laplace transform,

There are transform equations as follows:

The is the point source solution of the fractures, with the property of point source function to rewrite (48), and the result can be obtained as follows:where

Two-phase dimensionless characteristic production rate and other two parameters can be defined as follows:

Thus, (49) can be rewritten as follows:

Integrating (52) with respect to over the interval 0 to results inwhere

As we can see in Figure 2, every fracture is divided into equal segments in direction. Each discrete segment has a midpoint, () is the coordinate of one discrete segment, and and are the start point and the end point of one discrete segment.

Integrating (53) with respect to over the interval to , the result is as follows:

The dimensionless pseudopressure response for fractures can be expressed as follows:

Fluid flow in wellbore and fractures is under the assumption of infinite conductivity, so the pressure at each point within the fractures and the horizontal wellbore is identical to bottom hole pressure, and there is an expression as follows:

With the assumption of constant flow rate, there is an expression as follows:

With (58) and (56), we can establish a matrix equation group. Employing Duhamel’s theorem, the accurate solution incorporating the influence of skin factor can be obtained as follows:

The solution in Laplace space to the two-phase transient flow model at constant wellbore pressure production can be expressed as

#### 5. Analysis of Type Curves Characteristics

##### 5.1. The Type Curves of Pseudopressure Response and Production Decline

The type curves of a two-phase multifractured horizontal well in SGR can be plotted by Stehfest [31] numerical invention algorithm. Intuitively, the type curves graphically show the process and characteristics of fluid flow in reservoirs.

In the condition of sealed boundary, the two-phase pseudopressure-transient pressure response of multistage fractured horizontal well in SGR, which has six flow regimes, can be clearly shown in Figure 3. Regime 1 is the pure wellbore storage regime; the pseudopressure curve and its derivative curve have the same trend as conventional gas reservoir. Regime 2 is hump transition regime, which is caused by fluids difference in two-phase flow; the shape of the derivative curve looks like a concave. Regime 3 is early linear flow regime; the flow energy in this regime is provided by free gas flow. Formation pressure is still over desorption critical pressure, so the impact of desorbed gas cannot be seen in this period. Regime 4 is early radial flow regime, which behaves as a horizontal line in type curve. The clarity of this period primarily depends on the fracture length and spacing. When facture spacing is small enough and fracture length is long enough, this period may not be so obvious in type curve. Regime 5 is desorbed gas transient diffusion flow regime, and the derivative curve looks like an “arch,” which is caused by desorption. It is obviously influenced by absorption index. Regime 6 is late-time compound flow regime, and because of sealed boundary condition, it behaves like an up-warp in the type curve.

As we can see in Figure 4, the dimensionless production decline type curve is plotted by using the Stehfest numerical invention algorithm and Blasingame production decline analysis method. The dimensionless production decline type curve can be divided into six flow stages: Stage 1: wellbore storage which has the same characteristics as the conventional reservoir; Stage 2: formation linear flow period, which is characterized by a horizontal line on both dimensionless production integral and its derivative curves; Stage 3: interporosity flow period, which is characterized by a V shape in the derivative curve, and in this period free gas has exhausted, and desorbed gas from the shale matrix surface diffuses to factures by pressure difference; Stage 4: fracture linear flow period, which is characterized by a drop line in type curves; Stage 5: compound radial flow period, which is characterized by a drop line in dimensionless production integral curve and a rising line in dimensionless production derivative curve; Stage 6: boundary control flow period, which is characterized by a rapid drop line on both dimensionless production integral and its derivative curves for the reason of sealed boundary condition.

##### 5.2. Sensitivity Analysis of Production Decline

Figures 5–9 show the type curves of production decline analysis for a two-phase multifractured horizontal well in SGR. These figures are the typical response of two-phase multifractured horizontal well percolation in SGR at constant wellbore pressure production.

Figure 5 shows the production decline type curve characteristics affected by the initial saturation. It can be seen from the figure that the curves of the dimensionless production rate integral and its derivative curves go down with the increase of gas initial saturation.

Figure 6 shows the production decline type curve characteristics affected by the skin factor. It can be seen from the figure that the curves of the dimensionless production rate integral and its derivative curves go down with the increase of skin factor. This characteristic can be obviously seen in the early two stages: wellbore storage period and formation linear flow period.

Figure 7 shows the production decline type curve characteristics affected by the absorption index. It can be seen from the figure that the interporosity flow period may appear earlier with the increase of skin factor, and the fracture linear flow period will last even longer. The effort of absorption index cannot be observed in early product period and later product period.

Figure 8 shows the production decline type curve characteristics affected by the different fracture stages. It can be seen from the figure that the curves of the dimensionless production rate integral and its derivative curves go down with the decrease of fracture stages. More fractures would make flow condition near the horizontal well better, but if there are relative large fracture stages, the increase rate of the type curves may not be so obvious.

Figure 9 shows the production decline type curve characteristics affected by the horizontal well lateral length. It can be seen from the figure that the effort of horizontal well lateral length can be obviously seen in the early two stages, especially when the horizontal well lateral length is relatively long. If the horizontal well lateral length is long enough, the dimensionless production rate derivative curve may go down later and be smooth in the early flow period.

#### 6. Conclusions

The transient flow model of two-phase multifractured horizontal well and desorbed gas transient diffusion in SGR is established and solved, the type curves of pressure response and production decline are illustrated, and the production decline performance under the effects of hydraulic fractures and SGR properties are discussed. The following conclusions were obtained.(1)The two-phase multifractured horizontal well model can be solved using the eigenvalue method and orthogonal transformation. The desorbed gas transient diffusion rate can be solved using Langmuir’s equilibrium sorption isotherm and Fick’s second law.(2)The type curves of production decline are dominated by initial saturation, skin factor, and fracture stages. The curves of the dimensionless production rate integral and its derivative curves go down with the increase of gas initial saturation. The skin factor has primarily effect on the early two stages; the larger the skin factor is, the lower the derivative curves are. More fracture stages can improve the flow condition near a horizontal well, and the number of stages indicate the value of fracture spacing; larger fracture spacing increases the duration of fluid flow around individual fractures, and the more the fracture stages are, the more the dimensionless production rate is.(3)The absorption index has a primary effect on the starting time of the interporosity flow period and the depth of the trough. The trough in the type curves is an indication of desorbed gas diffusion in SGR; the more the absorption index is, the later the decline appears, and, in addition, the trough becomes more significant with larger value of absorption index.(4)The horizontal well lateral length has a primary effect on the early two flow stages but does not affect the intermediate and later flow behavior; the smaller the horizontal well lateral length is, the faster the decline appears.

#### Nomenclature

: | Length of SGR (m) |

: | Width of SGR (m) |

: | Total mass concentration (kg/m^{3}) |

: | Wellbore storage coefficient (kg/MPa) |

: | System compressibility (kg/(m^{3}·MPa)) |

: | Diffusion coefficient (m^{2}/s) |

: | Characteristic function, dimensionless |

: | Geometric factor, dimensionless |

: | Thickness of SGR (m) |

: | Absolute permeability (μm^{2}) |

: | Gas permeability in fractures, dimensionless |

: | Liquid permeability in fractures, dimensionless |

: | Horizontal well lateral length (m) |

: | Half length of fractures (m) |

: | Initial pseudopressure (kg/(m^{3}·s)) |

: | Pseudopressure (kg/(m^{3}·s)) |

: | Langmuir constant (kg/(m^{3}·s)) |

: | Orthogonal factor in -axis direction, dimensionless |

: | Fracture stages, dimensionless |

: | Pressure in fractures (MPa) |

: | Langmuir pseudopressure (MPa) |

: | Standard state pressure (MPa) |

: | Production of gas flow (m^{3}/s) |

: | Two-phase characteristic production rate (m^{3}/s) |

: | Transient diffusion rate (kg/(m^{5}·s)) |

: | Total mass flow rate (kg/s) |

: | Production of water flow (m^{3}/s) |

: | Radial distance in matrix block (m) |

: | Gas constant (0.008314 MPa·m^{3}/(kmol·K)) |

: | Laplace transform parameter, dimensionless |

: | Skin factor, dimensionless |

: | Time (s) |

: | Standard state temperature (K) |

: | Langmuir volume (m^{3}) |

: | Space coordinates in Cartesian coordinates (m). |

*Greek Symbols*

: | Absorption index, dimensionless |

: | Orthogonal factor in -axis direction, dimensionless |

: | Orthogonal factor in -axis direction, dimensionless |

: | External radius of matrix block (m) |

: | Eigenvalue, dimensionless |

: | Gas viscosity (MPa·s) |

: | Liquid viscosity (MPa·s) |

: | Gas density (kg/m^{3}) |

: | Liquid density (kg/m^{3}) |

: | Porosity of fractures, dimensionless. |

*Abbreviations*

SGR: | Shale gas reservoir. |

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

This study was supported by National Science Fund for Distinguished Young Scholars of China (Grant no. 51125019), National Natural Science Foundation of China (Grant No. 51304165), and the 2014 Australia China Natural Gas Technology Partnership Fund Top Up Scholarships.