**1. Introduction**

Natural gas is an indispensable energy resource which plays an important role on the world's energy map. As global demand continues to increase, conventional sources of natural gas will not be able to meet the world's requirements. Consequently, the oil industry has shifted its research focus from conventional gas to unconventional natural gas, especially tight gas. Recently, alongside continuous developments in petroleum geology and the boost to oil and gas industry technologies, the exploitation and development of tight gas reservoirs has made considerable breakthrough and become the primary growth point

for oil and gas production [1–6]. Because of its strong heterogeneity, complicated microscopic pore structure, and poor connectivity of the e ffective sand bodies, the porous flow mechanism of tight gas reservoirs is obviously distinct from those of conventional reservoirs. Moreover, because of poor reservoir physical properties, small discharging radius, lack of or low natural productivity, high development di fficulty, and some other potential disadvantages, a MFHW with hydraulic stimulations is extensively used for developing tight gas reservoirs [7–16].

To date, many scholars are dedicated to the transient productivity of fractured wells in unconventional hydrocarbon reservoirs, considering the di fferent characteristics of nonlinear flows compared to multi-scale flows. The methods to calculate the transient productivity of fractured horizontal wells primarily include analytical and numerical solutions, among which the analytical solution includes the complex potential theory, conformal transformation, and the equivalent flow resistance method [14,17,18]. Research on the transient productivity of tight gas reservoirs is increasingly wide. Through an iterative algorithm, Zheng et al. [19] quantified the correlation between saturation and pressure of the infinitesimal coal. They used the Corey relative permeability model to describe relative gas/water permeability as a function of the pressure. By applying the inter-porosity flow equation based on a pseudo-pressure, Zhenzihao et al. [20] developed a density-based rate-transient analysis technique, using which they predicted the gas production rate and were able to assess initial gas content in the system. This methodology is also capable of accurately forecasting gas production rates by converting the response of the corresponding liquid (separating pressure-dependent e ffects), which was performed using dimensionless depletion-driven parameters. Feng et al. [21] also developed a model thoroughly analyzing MFHWs in heterogeneous gas reservoirs. For this purpose, they used source functions, Green's solution, and the element method to constrain a model in Laplace coordinates. However, by combining the pseudo-function method and material balance to describe homogeneous gas reservoirs, one cannot eliminate the nonlinearity of non-homogeneous gas reservoirs. To solve this issue, they calculated material balance by taking into account heterogeneous gas reservoirs from di fferent regions. Jie et al. [22] created a new analytical model, which considers fracture damage and complex gas transport mechanisms as well as the shale gas stimulated reservoir volume (SRV). Berawala et al. [23] presented a mathematical 1D + 1D model to investigate the main controlling factors during a continuum-flow regime in shale-gas production in the context of well-induced fractures. Wenxi et al. [24] simulated shale gas production using a multi-stage fractured horizontal well using another model established in their study, which was based on a trilinear flow and did not assume that the secondary fracture length was equal to half the distance between hydraulic fractures. Jinzhou et al. [9] calculated gas productivity for the multi-staged horizontal wells in fractured tight sandstone gas reservoirs using a semi-analytical mathematical approach. Their model considered how flow in the fracture network (natural fracture system) a ffected the outcome. This helped to create a production forecast model of complex hydraulics with simulated impacts arising from natural fractures. As a result, after using the Gaussian elimination method, they obtained a semi-analytical solution.

Most previous studies simplified the nonlinear flow mechanisms of tight reservoirs to some extent. It is di fficult to accurately and comprehensively reflect the actual flow state of fluid flowing into horizontal wells. To fill this gap, considering the seepage properties of the tight gas reservoirs and the steady-state seepage of fractured horizontal wells, we introduce the transient model of a matrix gas discharge radius. Considering the influence of various factors of tight reservoirs on the productivity, we also introduce the equivalent well diameter model. Under the same conditions, the pressure distribution equation of a single fracture is equivalent to the vertical well pressure distribution equation. Moreover, the transient productivity model describing a multi-stage fractured horizontal well is obtained. The actual data from a gas field in China were applied to verify the productivity model. It was demonstrated that two results obtained by the method proposed in this paper and dominant commercial software have an error of less than 1.62%, and the single parameter sensitivity analysis was also carried out. This paper provides a useful and feasible tool for reservoir engineers to predict, evaluate, and optimize the productivity of multi-fractured horizontal wells in tight gas reservoirs.

#### **2. Reservoir Characteristics Analysis and Physical Model Assumptions**

#### *2.1. Productivity Forecast Model for Individual Horizontal Wells*

As shown in Figure 1, for fractured horizontal wells, different production stages occur in different seepage areas within different flow media. Within the control area of individual wells, gas located in the fracture first flows to the wellbore. Next, the matrix near the wellbore area re-supplies the fracture. Therefore, the gas in the matrix flows toward the wellbore through the fracture. Then, the matrix far from the wellbore area supplies the matrix near the wellbore, and gas in the far wellbore area flows to the fracture through the matrix throat near the wellbore.

**Figure 1.** Schematics for each production period of the fractured horizontal well.

To achieve an accurate productivity prediction, based on the characteristics of different production stages of fractured horizontal wells in tight gas reservoirs, different flow media in each seepage area, and different seepage mechanisms, it is therefore necessary to establish corresponding productivity equations by considering the influence of factors such as: (a) threshold pressure gradient; (b) stress sensitivity effect; (c) slippage effect; and (d) high-speed non-Darcy effect. This comprehensive model is able to predict the productivity of the fractured horizontal wells accurately.

In the initial high-production stage, the artificial fractures are the primary flow medium, and hence the gas seepage occurs primarily in the artificial fractures. During this stage, the high-speed non-Darcy seepage cannot be ignored. Moreover, both stress-sensitive effects in artificial fractures and slippage effects play important roles. In the transition period, the gas seepage area is primarily in the proximity of the wellbore (SRV) and the flow medium is the matrix near the wellbore. In this stage, low-speed non-Darcy seepage plays a significant role, which is primarily affected by the threshold pressure gradient and the stress sensitive effect in the matrix. In the final stage of production, the seepage area is far from the wellbore (outside the SRV), and the flow medium is primarily the matrix far from the wellbore. In this stage, low-speed non-Darcy seepage plays a predominant role. However, the influence of the threshold pressure gradient, slippage effect, and stress sensitive effect in the matrix should be considered as well.

#### *2.2. Physical Model Assumption*

As shown in Figure 2, a fractured horizontal well with perforation completion is bored in a horizontal, infinite, homogeneous, and equal thickness tight gas reservoir. Some basic assumptions made for multi-stage fractured horizontal wells are as follows:

(1) This model is applicable for isothermal single-phase unstable flows and the influence of gravity is neglected.


**Figure 2.** The schematic of seepage field in fissures of fractured horizontal well.

#### **3. Mathematical Model**

#### *3.1. Mathematical Model for Nonlinear Flow Mechanisms*

Considering the threshold pressure gradient, stress sensitive effect, and slippage effect, the nonlinear flow mathematical model of tight gas reservoirs can be derived through the generalized form of Darcy's Law [25,26].

#### 3.1.1. Threshold Pressure Gradient Effect

Because of the differences of porosity structure characteristics, when gas receive the effect of surface molecules force, this will in turn contribute to the threshold pressure phenomena. Consequently, gas flow in low tight gas reservoirs obviously differs from that in conventional gas reservoirs with medium or high permeability [27]. The threshold pressure gradient, which is associated with non-Darcy flow in low-permeability reservoirs, is defined as the level of pressure gradient that must be attained to enable the fluid to flow.

$$w\_{\mathfrak{E}} = \frac{\mathbb{K}\_{\mathfrak{E}}}{\mu\_{\mathfrak{E}}} (\frac{dp}{d\overline{y}} - \lambda) \tag{1}$$

where *<sup>v</sup>*g is velocity of gas, m/s; *<sup>K</sup>*g is permeability measured by gas, mD; μg is gas viscosity, mPa·s; *p* is the pressure in the matrix, MPa; and λ is the threshold pressure gradient factor, MPa.

#### 3.1.2. Stress Sensitivity Effect

Pressure depletion easily causes the deformation of the pores across different scales. When affected by the failure of the propellant, the artificial fracture is deformed or even completely closed [28,29]. The relation of permeability stress sensitivity because of matrix pore throat deformation is as follows [30]:

$$K\_{\rm m} = K\_{\rm m0} e^{-\alpha\_{\rm F}(p\_i - p)} \tag{2}$$

where *K*m is permeability of matrix, mD; *K*m0 is the initial permeability of matrix, mD; α*F* is the coefficient of stress sensitivity, MPa−1; *pi* is the initial formation pressure, MPa; and *p* is current pressure, MPa.

#### 3.1.3. Gas Slippage Effect

Slippage is a phenomenon in which natural gas from a reservoir bypasses crude oil and water that is released from the capillary opening of porous oil formations and approaches to the mean free path of the natural gas. Gas slippage can be defined as the gas movement through liquid phase of the reservoir front. The gas slippage effect is an important factor affecting gas flow in compact porous media [31–33]. We used the Klinkenberg equation to describe its dynamic characteristics [34,35]:

$$K\_{\mathbb{R}} = \mathbb{K}(1 + b/\overline{p}) \tag{3}$$

where *<sup>K</sup>*g is the permeability measured by gas, mD; *b* is the slippage factor, MPa; *p* is the average gas reservoir pressure, MPa; and *K* is the permeability of a different medium, mD.

#### *3.2. Productivity Model for Porous Flow between Matrix and Fractures*

#### 3.2.1. Steady-state Productivity Model

There are two different porous flow zones: the matrix porous flow zone and the fracture porous flow zone. In the process of porous flow from matrix to fracture, the porous flow medium is the matrix. Each fracture can be simplified as a linear source. The two-dimensional non-Darcy elliptic flow model is assumed to occur in the formation; i.e., the conjugate isopiestic elliptic cylinder and hyperboloid streamline with the well as the center and fracture endpoint as focus is achieved [36]. When a well starts producing, an equal pressure elliptic cylinder is created in the formation (Figure 3), and the relationship between the rectangular coordinate and elliptical coordinate is defined as follows:

$$\mathbf{x} = a\cos\eta, \mathbf{y} = b\sin\eta \tag{4}$$

$$a = \mathbf{x}\_{\overline{F}} \cosh \xi, b = \sinh \xi \tag{5}$$

where *x* and *y* are Cartesian coordinates; *a* and *b* are major axis and minor axis, respectively; *xF* is the half length of fracture, m; and η and ξ are elliptical coordinates, m.

**Figure 3.** Flow field of the elliptic cylinder with a constant pressure.

By combining Equations (4) and (5), an isobaric elliptical equation and hyperbolic streamline equations can be expressed as Equation (6) and (7), and its flow field is illustrated in Figure 3:

$$\frac{x^2}{a^2} + \frac{y^2}{b^2} = 1\tag{6}$$

$$\frac{\mathbf{x}^2}{\left(\mathbf{x}\_F \cos \eta\right)^2} - \frac{y^2}{\left(\mathbf{x}\_F \sin \eta\right)^2} = 1\tag{7}$$

Combining Equations (4) and (5), the relationship between *p* and *y* under Cartesian coordinates can be transformed to the relationship between *p* and ξ under elliptical coordinates.

$$\frac{dp}{dy} = \frac{dp}{d\xi} \cdot \frac{d\xi}{dy} = \frac{\pi}{2x\_F \cosh \xi} \cdot \frac{dp}{d\xi} \tag{8}$$

The average mass flow rate of the cross section of the elliptical column in *y* direction is as follows:

$$q = \frac{45\pi a}{32} h \cdot v\_{\text{m}} \tag{9}$$

where *q* is the gas flow rate, m3/s; *h* is the reservoir thickness, m; and *v*m is the gas velocity, m/s.

By combining Equations (8)–(10), a porous flow equation considering stress sensitivity, threshold pressure gradient, and slippage effect can be expressed as follows:

$$\frac{\mathbb{K}\_{\text{m}0} \exp\left[-\alpha(p\_{\text{i}} - p)\right]}{\mu} \cdot (1 + b/\overline{p}) \cdot (\frac{\pi}{2\text{x}\_{\text{F}} \cosh \xi} \cdot \frac{\text{d}p}{\text{d}\xi} - \lambda) = \frac{ZT}{p} \frac{p\_{\text{sc}}}{\text{Z}\_{\text{sc}} T\_{\text{sc}}} \frac{32q\_{\text{sc}}}{45\pi \text{x}\_{\text{F}} \hbar \cosh \xi} \tag{10}$$

where *Z* is the gas compressibility factor, dimensionless; *T* is the reservoir temperature, *K*; *Zsc* is the gas compressibility factor under standard conditions, dimensionless; *Tsc* is the standard temperature, *K*; and *psc* is the standard pressure, MPa.

Assuming *f*(*p*) = *p* exp<sup>α</sup>(*p*−*pi*) μ*Z* , the pseudo-pressure function equation is obtained.

$$
\psi(p) = \int\_{p\_0}^p f(\delta)d\delta \tag{11}
$$

where δ is variable, MPa; *p*0 is the initial pressure, MPa; and *p* is the current pressure, MPa.

By substituting Equation (11) into Equation (10), a simplified porous flow equation is derived as follows:

$$
\psi(p\_i) - \psi(p) = \frac{64p\_{\text{sc}}T \int\_{\xi(0)}^{\xi(r\_\ell)} d\xi}{45\pi^2 K\_{\text{m0}} Z\_{\text{sc}} T\_{\text{sc}} (1 + b/\overline{p})h} \cdot q\_{\text{sc}} + \frac{2\mathbf{x}\_\text{F}}{\pi} \int\_{\xi(0)}^{\xi(r\_\ell)} \lambda f(\overline{p}) \cosh \xi d\xi \tag{12}
$$

where ψ(*pi*) is the initial pseudo-pressure, MPa<sup>2</sup>/(mPa·s); ξ(*re*) is variable upper limit, dimensionless; and ξ(0) is variable lower limit, dimensionless.

3.2.2. Transient Model of Discharge Radius in the Matrix

The discharge radius of the tight gas reservoir matrix has a non-transient effect and will expand as the pressure wave propagates. The entire area of gas flow at each instant actually involves the whole formation. Thus, the gas reservoir can be divided into two zones: the pressure-related zone and the un-flushed zone. The radius of the excitation zone (i.e., the radius of the matrix's deflation) increases with time. The motion law of the dynamic zone and the non-excited zone can be obtained using the material balance equation with boundary conditions.

From Equations (4) and (5), the relationship between ξ and *a* can be expressed as follows:

$$\mathcal{L} = \ln[(a + \sqrt{a^2 - \mathbf{x}\_F^2})/\mathbf{x}\_F] \tag{13}$$

Figure 4 is the schematic of elliptical flow area which illustrates the relationship between the major axis of ellipse *a* and the discharge radius of matrix *<sup>R</sup>*(*t*). Note that *a* and *R*(*t*) can be expressed as follows:

$$a = \mathbf{x}\_{\overline{F}} + R(t) \tag{14}$$

**Figure 4.** Schematic of the elliptical seepage area.

Note that the discharge radius of matrix *R*(*t*) is a transient value, which increases over time t, and the material balance equation is then used to determine *<sup>R</sup>*(*t*).

The original gas reserves in the region with radius *R*(*t*) can be calculated as follows:

$$M\_0 = \pi [R^2(t) - r\_w^2] h \phi \frac{\rho\_{dT}}{p\_{dT}} p\_a \tag{15}$$

where *M*0 is the initial reserve in place, kg; *rw* is the well radius, m; φ is the porosity, dimensionless; ρ*aT* is the gas density at standard condition, kg/m3; *paT* is the standard pressure, MPa; and *pa* is the real formation pressure, MPa.

The current reserves can be expressed using the average formation pressure *p*:

$$M\_{\rm f} = \pi [R^2(t) - r\_w^2] h \phi \frac{\rho\_{\rm dT}}{p\_{\rm dT}} \overline{p} \tag{16}$$

Moreover, the average formation pressure *p* can be calculated using the following equation:

$$\overline{p} = p\_i - \frac{p\_i^2 - p\_{wf}^2}{4p\_i \ln[R(t)/r\_w]} \tag{17}$$

where *pw f* is the bottom hole flow pressure, MPa.

Because the gas production rate is constant at *QaT*, the total mass of production gas is equal to ρ*aTQaTt* at time *t*. Therefore, the conservation law of matter can be expressed as follows:

$$M\_0 - M\_t = \rho\_{aT} Q\_{aT} t \tag{18}$$

Furthermore, the productivity equation of steady state seepage of tight gas reservoirs considering the threshold pressure gradient can be obtained as follows:

$$Q\_{dT} = \frac{\pi kh \left[p\_i^2 - p\_{wf}^2 - \lambda \overline{p} (R(t) - r\_{w\circ})\right]}{\mu p\_{d\,T} \ln[R(t)/r\_w]}\tag{19}$$

Combining Equations (16)–(19), the discharge radius equation, which considers the threshold pressure gradient, can be obtained as follows:

$$R^2(t) - r\_w^2 = \frac{4kp\_i}{\imath \phi} \cdot \frac{p\_i^2 - p\_{wf}^2}{\lambda \overline{p} [R(t) - r\_w]} t \tag{20}$$

Moreover, taking the stress sensitivity effect into consideration, the relationship between *R*(*t*) and *t* is expressed as:

$$R^2(t) - r\_w^2 = \frac{4k\_{\rm m0}e^{-\alpha(p\_i - \overline{p})}p\_i}{\nu \phi} \cdot \frac{(p\_i^2 - p\_{wf}^2) - \lambda \overline{p}[R(t) - r\_w]}{p\_i^2 - p\_{wf}^2}t \tag{21}$$

Placing *R*(*t*) = *e*ξ(*t*) and *rw* = 1 into Equation (21), the relationship between ξ(*t*) and *t* can thus be written as:

$$\epsilon^{2\mathcal{E}\_{\epsilon}(t)} - 1 = \frac{4k\_{\mathrm{m}0}e^{-\alpha(p\_{\mathrm{i}}-\overline{p})}p\_{\mathrm{i}}}{\mu\phi} \cdot \frac{(p\_{\mathrm{i}}^2 - p\_{wf}^2) - \lambda\overline{p}[e^{\mathcal{E}\_{\epsilon}(t)} - 1]}{p\_{\mathrm{i}}^2 - p\_{wf}^2}t \tag{22}$$

$$\overline{p} = p\_i - \frac{p\_i^2 - p\_{wf}^2}{4p\_i \xi\_c(t)}\tag{23}$$

Combining Equations (12), (22), and (23), an unstable-state productivity model from matrix to fracture can be obtained as follows:

$$
\psi(p\_i) - \psi(p) = \frac{64p\_{\text{sc}}T \cdot \int\_0^{\zeta\_c(t)} d\zeta}{45\pi^2 K\_{\text{m0}} Z\_{\text{sc}} T\_{\text{sc}} (1 + b/\overline{p}) h} \cdot q\_{\text{sc}} + \frac{2\mathbf{x}\_{\text{F}}}{\pi} \int\_0^{\zeta\_c(t)} \lambda f(\overline{p}) \cosh \xi d\zeta \tag{24}
$$

3.2.3. Transient Productivity Model for Fractured Horizontal Wells

The transient productivity model for fractured horizontal wells can be obtained by introducing the transient model of the discharge radius in the matrix. Assuming that the flow rate of gas in the gas reservoir-fracture is *qsc*1 (converted to ground standard condition), the pressure at the edge of the fracture is *p*1 and the corresponding pseudo-pressure is ψ(*p*1). Moreover, the pseudo-pressure distribution equation of transient seepage flow in gas reservoir-fracture is expressed as [37–39]:

$$
\psi(p\_i) - \psi(p\_1) = \frac{64p\_{\text{sc}}T \cdot \xi\_{\text{c}}(\text{t})}{45\pi^2 K\_{\text{m0}} Z\_{\text{sc}} T\_{\text{sc}} (1 + b/\overline{p})\overline{h}} q\_{\text{sc1}} + \frac{2\text{x}\_{\text{f}}}{\pi} \int\_0^{\xi\_{\text{c}}(t)} \lambda f(\overline{p}) \cosh \xi d\xi \tag{25}
$$

#### *3.3. Productivity Model for Porous Flows between the Fracture and Near the Wellbore Area*

In the fracture that is near wellbore seepage, the flowing medium is the fracture itself. Because of a high permeability of the fracture and the high velocity of gas flow in the fracture, Darcy's equation is not applicable, and rather the Forchheimer equation is required to describe the gas flow. Considering the effects of high velocity turbulence, stress sensitivity, and slippage, the mathematical model of gas seepage in rocks can be calculated using the following equations.

Because of high conductivity in artificial fractures, the gas flows quickly inside the fracture with characteristics of a turbulent flow, which conforms to the high-speed non-Darcy seepage principle [8]. The kinematic equation is then obtained as follows:

$$\frac{dp}{d\mathbf{x}} = \frac{\mu}{K\_F} \upsilon\_F + \beta\_\mathcal{\mathcal{S}} \rho\_\mathcal{\mathcal{S}} \upsilon\_F^2 \tag{26}$$

where *KF* is the permeability of fracture, mD; *vF* is the seepage velocity of gas in fractures, m/s; β*g* is the turbulence coefficient, m<sup>−</sup>1; and ρ*g* is the gas density, g/m3.

The stress-sensitive equation of fractures can be expressed as follows:

$$K\_F = K\_{F0}e^{-\alpha\_F(p\_i - p)}\tag{27}$$

where *KF* is the permeability of the fracture, mD; *KF*0 is the initial permeability of the fracture, mD; and α*F* is the coefficient of stress sensitivity, MPa−1.

The velocity of gas flows in the fracture can be obtained as follows:

$$w\_F = \frac{q\_2}{2w\_Fh} = \frac{Z}{P} \cdot \frac{p\_{sc}T}{Z\_{sc}T\_{sc}} \cdot \frac{q\_{sc2}}{2w\_Fh} \tag{28}$$

where *qsc*2 is the flow velocity of the linear flow zone in the fracture, m3/d; and *wF* is the width of fracture, m.

Assuming that the flow rate in linear flow zone is *qsc*2 (converted to standard surface conditions), the pressure at the interface between the linear flow zone and the radial flow zone (radius is h/2) is *p*2, and the pseudo-pressure is ψ(*p*2), the pseudo-pressure distribution equation of the gas reservoir near the wellbore unstable seepage flow is obtained as follows:

$$
\psi(p\_1) - \psi(p\_2) = \frac{p\_{\rm sc} T(\mathbf{x}\_{\rm F} - \mathbf{h}/2)}{2K\_{\rm F0} w\_{\rm r} \hbar Z\_{\rm sc} T\_{\rm sc}} \cdot q\_{\rm sc2} + \beta\_{\rm g} \frac{MTZ p\_{\rm sc}^2 (\mathbf{x}\_{\rm F} - \mathbf{h}/2) f(\overline{p})}{4R \overline{p} w\_{\rm F}^2 h^2 Z\_{\rm sc}^2 T\_{\rm sc}^2} \cdot q\_{\rm sc2}^2 \tag{29}
$$

where *qsc*2 is the flow rate in linear flow zone under standard surface conditions, m3/d.

β*g MTZp*<sup>2</sup>*sc*(*xF*−*h*/2)*f*(*p*) <sup>4</sup>*Rpw*<sup>2</sup>*Fh*2*Z*2*scT*2*sc* ·*q*<sup>2</sup>*sc*<sup>2</sup> is the additional pseudo-pressure drop caused by the high-speed turbulent velocity.

#### *3.4. Radial Porous Flow Model from Fracture to Wellbore*

When gas accumulates around the wellbore at the edge of fracture, the radial flow leads to an additional pressure drop near the wellbore. As shown in Figure 5, this phenomenon is called the radial concentration effect.

**Figure 5.** Schematic of porous flow from fracture to wellbore: (**a**) porous flow in one fracture of a vertical well; and (**b**) porous flow in one fracture of a horizontal well.

Assuming that the flow rate in radial flow zone is *qsc*3 (converted to standard surface conditions), the velocity of the radial flow can be calculated as follows:

$$w\_3 = \frac{q\_3}{2\pi r w\_\text{F}} = \frac{Z}{p} \cdot \frac{p\_{\text{sc}} T}{Z\_{\text{sc}} T\_{\text{sc}}} \cdot \frac{q\_{\text{sc3}}}{2\pi r w\_\text{F}}\tag{30}$$

where *v*3 is the velocity of radial flow in the fracture, m/s; *q*3 is the flow rate in the fracture under reservoir conditions, m3/d; and *qsc*3 is the flow rate in the fracture under surface condition, m3/d.

Combining Equations (25)–(30), the pseudo-pressure equation of fracture-wellbore transient flow is obtained.

$$
\psi(p\_2) - \psi(p\_{wf}) = \frac{p\_{\rm sc} T \ln(h/2r\_w)}{2\pi K\_{\rm F0} w\_{\rm F} Z\_{\rm sc} T\_{\rm sc}} \cdot q\_{\rm sc3} + \beta\_{\rm \%} \frac{MTZ(1/r\_w - 2/h)p\_{\rm sc}^2}{4\pi^2 R \overline{p} Z\_{\rm sc}^2 T\_{\rm sc}^2 w\_{\rm F}^2} \cdot q\_{\rm sc3}^2 \tag{31}
$$

#### *3.5. Productivity Model of a Single Fracture in Horizontal Wells*

The flow field of the fractured horizontal wells can be divided into two parts: (a) external seepage field (reservoir-fracture); and (b) internal seepage field (fracture-horizontal wellbore). According to the principle of hydropower similarity, the external seepage field will continuously supply gas to the internal seepage field. Therefore, the interfacial pressure and rate will be equal.

Considering the slippage effect, pressure sensitive effect, threshold pressure gradient, and high-speed non-Darcy turbulence effect, the productivity equation can be expressed by Combining Equations (25), (29), and (31).

$$\begin{split} \psi(p\_i) - \psi(p\_{wf}) &= \frac{2\chi}{\pi} \int\_0^{\xi\_c(t)} \lambda f(\overline{p}) \cosh \xi d\xi + \frac{p\_{\text{sc}}T}{2\zeta\pi^2 K\_{\text{sc}}(1+b/\overline{p})h} \left[ \frac{64\xi\_c(t)}{4\overline{\kappa}\eta^2 K\_{\text{m0}}(1+b/\overline{p})h} + \frac{\chi\overline{p}-h/2}{2\overline{\kappa}q\_0 w\_{\text{fl}}} + \frac{\ln(h/2w\_{\text{cr}})}{2\overline{\kappa}K\_{\text{l}}ww} \right] \cdot q\_{\text{sc}} \\ &+ \beta\_{\overline{\kappa}} \frac{MT2p\_{\text{sc}}^2}{4\overline{\kappa}pw\_{\text{f}}^2 \kappa\_{\text{sc}}^2 T\_{\text{sc}}^2} \left[ \frac{(x\eta - h/2)f(\overline{p})}{h^2} + \frac{(1/r\_w - 2/h)}{\pi^2} \right] \cdot q\_{\text{sc}}^2 \end{split} \tag{32}$$

#### *3.6. Equivalent Wellbore Radius Model*

The productivity equation of a complex well pattern under complex conditions, e.g., non-Darcy, can be obtained analogous to ordinary vertical wells under Darcy seepage conditions. Therefore, the equivalent diameter of a single fracture in a horizontal well with transverse fracturing can be obtained by combining Equation (32) with the ordinary vertical well productivity equation after considering the threshold pressure gradient, stress sensitivity effect, and slippage effect under generalized Darcy percolation conditions [40].

The pressure distribution equation of normal vertical wells is obtained from the generalized Darcy equation, stress sensitivity, and the slippage effect equation:

$$
\psi(p\_i) - \psi(p) = \frac{p\_{\rm sc} T \ln(R\_\varepsilon/r) \cdot q\_{\rm sc}}{2\pi K\_{\rm F0} Z\_{\rm sc} T\_{\rm sc} (1 + b/\overline{p}) h} + \lambda f(\overline{p}) (R\_\varepsilon - r) \tag{33}
$$

The pressure distribution equation of a normal vertical well can be expressed as follows:

$$
\psi(p\_i) - \psi(p\_{wf}) = \frac{p\_{\text{sc}} T \ln(R\_\text{c} / r\_{\text{equ}}) \cdot q\_{\text{sc}}}{2 \pi K\_{\text{FO}} Z\_{\text{sc}} T\_{\text{sc}} (1 + b/\overline{p}) h} + \lambda f(\overline{p}) (R\_\text{c} - r\_{\text{equ}}) \tag{34}
$$

Combining Equations (33) and (34), the equivalent diameter, *requ*, of a single fracture in the horizontal well is obtained.

#### *3.7. Productivity Model of Multi-fractured Horizontal Wells*

Multiplied vertical fractures in a horizontal well will interfere with each other, and the degree of mutual influence depends on the location of the fracture as well.

According to the superposition theory of the pseudo-pressure drop, MFHWs are equivalent to multiple vertical wells using the equivalent radius model [41,42]. Thus, the seepage flow of MFHW is transformed into the superposition of multi-fractured vertical wells. According to the pseudo-pressure superposition principle, the pseudo-pressure of each fracture at fracture *j* is as follows:

⎧

⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩

$$\begin{aligned} \psi(p\_i) - \psi(p\_{wf1}) &= \Delta\psi(p)\_{11}(q\_{\kappa1}) + \Delta\psi(p)\_{21}(q\_{\kappa2}) + \Delta\psi(p)\_{31}(q\_{\kappa3}) + \dots + \Delta\psi(p)\_{n1}(q\_{\kappa n})\\ \psi(p\_i) - \psi(p\_{wf2}) &= \Delta\psi(p)\_{12}(q\_{\kappa1}) + \Delta\psi(p)\_{22}(q\_{\kappa2}) + \Delta\psi(p)\_{32}(q\_{\kappa3}) + \dots + \Delta\psi(p)\_{n2}(q\_{\kappa n})\\ \psi(p\_i) - \psi(p\_{wf3}) &= \Delta\psi(p)\_{13}(q\_{\kappa1}) + \Delta\psi(p)\_{23}(q\_{\kappa2}) + \Delta\psi(p)\_{33}(q\_{\kappa3}) + \dots + \Delta\psi(p)\_{n3}(q\_{\kappa n})\\ &\vdots\\ \psi(p\_i) - \psi(p\_{wfn}) &= \Delta\psi(p)\_{1n}(q\_{\kappa1}) + \Delta\psi(p)\_{2n}(q\_{\kappa2}) + \Delta\psi(p)\_{3n}(q\_{\kappa3}) + \dots + \Delta\psi(p)\_{nn}(q\_{\kappa n}) \end{aligned} \tag{35}$$

Assuming that gas seepage in the horizontal wellbore is infinite and that pressure in the horizontal wellbore is balanced,

$$
\psi(p\_{wf1}) = \psi(p\_{wf2}) = \dots = \psi(p\_{wfn})\tag{36}
$$

Thus, the productivity of a fractured horizontal well is as follows:

$$q\_{\rm sc} = q\_{\rm sc1} + q\_{\rm sc2} + \dots + q\_{\rm sci} = \sum\_{i=1}^{n} q\_{\rm sci} \tag{37}$$

Rewriting Equation (35) for each fracture allows us to form a closed system equation, which can be solved to obtain the gas productivity of a multi-stage fractured horizontal well in a tight gas reservoir.

#### **4. Results and Analysis**

#### *4.1. Model Validation*

To verify the accuracy of the transient productivity forecast model, a numerical model of MFHW was established using CMG. By assuming the same parameters for a certain tight gas reservoir (Table 1), the daily production rate, cumulative production of the model established in this study, and CMG were obtained.


**Table 1.** Fundamental parameters of a gas field at Jilin Oilfield in China.

Figure 6 shows the daily production rate and cumulative production simulated by the model proposed in this paper and CMG, as well as the field production data. The data simulated by the proposed model and CMG are in good agreemen<sup>t</sup> with the field data. Meanwhile, the curves of the daily production rate simulated by the proposed model and CMG have a high degree of conformity and a consistent declining trend. At the end of simulation (1500 days), cumulative production simulated by the proposed model and CMG is 11,222 × 10<sup>4</sup> m<sup>3</sup> 11,404 × 10<sup>4</sup> m3, respectively. The results obtained by the two methods have an error of less than 1.62%. In this way, we demonstrated that the proposed model is accurate enough to simulate production of a multi-stage fractured horizontal well in tight gas reservoirs. During the field production process, the gas well is sometimes shut down. Therefore, when production begins again, the actual production rate is slightly higher than the simulated results because of the recovery of formation pressure facilitated by the shutdown.

**Figure 6.** Comparison of the gas production rate of the model in this study and the CMG numerical model.

Figure 7 shows the reservoir pressure distribution of the transient productivity model after 1800 days of production. The gas seepage shows an elliptical flow pattern around the fractures. With continuous production, the pressure wave continues to spread outwards and the discharge radius continues to expand, gradually forming an elliptical stimulated zone that considers the horizontal wellbore as the horizontal major axis and the fractures as the vertical minor axis.

**Figure 7.** Reservoir pressure distribution of the multi-stage fracture well production.

In this case, we introduce a new concept called the *contribution degree* to evaluate the influence of different factors to productivity growth. Contribution degree (a dimensionless quantity) is defined as the increase in production rate caused by certain factors divided by the production rate when the factor is neglected. The contribution degree of different factors to productivity growth at different times was analyzed, as illustrated in Figure 8. Except for the seepage effect, all of these factors have a negative effect on the output growth. In the early stage of production, the turbulence effect and the stress sensitive effect have considerable influence over gas production, while the effect of a threshold pressure gradient is insignificant. With continuous gas production, the effect of the threshold pressure gradient becomes significant, while the turbulence and stress sensitive effects weaken. It is practical to consider the turbulence effect due to high-speed seepage in the artificial fractures. In general, the increasing influence of contributing factors of production is as follows: (a) threshold pressure gradient; (b) stress

sensitivity; (c) turbulence effect; and (d) slippage effect. The slippage effect of gas only works in the low-pressure stage, and its influence is small.

**Figure 8.** A schematic of the contribution of different factors to productivity at different periods.

#### *4.2. The Influence of Seepage Mechanisms on Gas Production*

Different threshold pressure gradients have different effects on gas well productivity. Figure 9 shows that the larger is the threshold pressure gradient, the more significant is the impact on gas well productivity. When the threshold pressure gradient ranges from 0 to 0.6 MPa/m, it has relatively little effect on gas production, compared with a range around 0.6–0.8 MPa/m. While the gas production is primarily derived from fractures, there is no threshold pressure gradient in the initial stage; therefore, the daily production rate at the beginning is almost the same as that at 46 × 10<sup>4</sup> m<sup>3</sup>/d. With a gradual decrease of production until the stable period, or the final stage, the influence of the threshold pressure gradient becomes increasingly significant. When the production time is 500 days, the daily gas production (without considering the threshold pressure gradient) is 34 × 10<sup>4</sup> m<sup>3</sup>/d; however, a daily gas production with a threshold pressure gradient of 0.8 MPa/mis 10 × 10<sup>4</sup> m<sup>3</sup>/d. At the end of the production, the daily gas production (without considering the threshold gradient) is stable at about 33 × 10<sup>4</sup> m<sup>3</sup>/d, and the cumulative gas production is 5.2 × 10<sup>8</sup> m3. When the threshold pressure gradient is 0.8 MPa/m, the daily gas production is about 8.2 × 10<sup>4</sup> m<sup>3</sup>/d, and the cumulative gas production is 1.59 × 10<sup>8</sup> m3. Thus, the threshold pressure gradient has a considerable influence on productivity in the middle and final stages and hence cannot be ignored for predicting the production of tight gas reservoirs.

**Figure 9.** The effect of different slippage factors on productivity.

Because the pore size of the tight reservoir is extremely small, the seepage capacity of the medium is obviously sensitive to the pressure, and the deformation of the medium has considerable influences on the reservoir's properties. Figure 10 shows the effect of matrix stress sensitivity on gas productivity.

The figure indicates that the matrix stress sensitivity coefficient has negative effects on daily gas production; the larger is the pressure sensitivity coefficient, the smaller is the impact on productivity. When the stress sensitivity of the matrix is not considered, the gas production rate at the start is 43.7 × 10<sup>4</sup> m<sup>3</sup>/d, and the high-yield period is considerably long, reaching up to 100 days before entering the transition stage. When production time reaches 1000 days, production gradually transforms to the steady flow stage. At this time, the daily gas production remains stable at 8.5 × 10<sup>4</sup> m<sup>3</sup>/d, and the cumulative gas production at the end of production can reach 1.03 × 10<sup>8</sup> m3. When the coefficient of matrix stress sensitivity is 0.8 MPa−1, the initial production rate is only 6.5 × 10<sup>4</sup> m<sup>3</sup>/d, indicating that the stress sensitivity coefficient has particularly considerable impacts in the initial stage. When the production time reaches 100 days, the well productivity stabilizes at 1.8 × 10<sup>4</sup> m<sup>3</sup>/d.

**Figure 10.** The effect of different matrix stress sensitivity coefficients on productivity.

#### *4.3. The Influence of Formation Properties on Gas Production*

The matrix permeability of different tight reservoirs varies. Therefore, it is important to evaluate the impact of matrix permeability on productivity, as illustrated in Figure 11. Below, we discuss the influence of matrix permeability of 0.05, 0.1, 0.15, 0.2, and 0.25 × 10−<sup>3</sup> μm<sup>2</sup> on daily and cumulative gas production. Because gas is supplied by large fractures in the initial stage of production (without the participation of the matrix in the seepage process), the daily gas production at the initial stage of production remains virtually constant, i.e., 46 × 10<sup>4</sup> m<sup>3</sup>/d. At the middle and end of the production stage, gas is primarily supplied by pores in the matrix; therefore, the influence of matrix permeability on productivity gradually appears. As matrix permeability increases, the influence on gas wells becomes considerable. When the matrix permeability is 0.25 × 10−<sup>3</sup> μm2, the daily production rate is 20.3 × 10<sup>4</sup> m<sup>3</sup>/<sup>d</sup> at the final stage of production. When the matrix permeability is 0.05 × 10<sup>4</sup> m<sup>3</sup>/d, the daily production rate stabilizes at 9.5 × 10−<sup>3</sup> μm2, which is reduced by 53.2% compared to permeability of 0.25 × 10−<sup>3</sup> μm2. Moreover, the final cumulative gas production is only 1.83 × 10<sup>4</sup> m<sup>3</sup>/d/d, which is reduced by 88.5% compared to permeability of 0.2 × 10<sup>4</sup> m<sup>3</sup>/d. These results show that matrix permeability has considerable effects on productivity, and the influence of the matrix permeability should be considered in the production allocation and productivity prediction of new wells.

Figure 12 demonstrates that formation thickness has a considerable influence on gas productivity in the initial and stable production stages. A larger formation thickness induces larger gas productions at the initial and stable production stages. With the gradual increase in thickness, the previously increased range of stable productivity remains the same. Figure 12 shows the changing curves for daily production rate and cumulative production as a function of time—corresponding to different formation thicknesses, e.g., 5, 14, 23, 32, 25, 41, and 50 m. When formation thickness is 5 m, the initial production rate is 40 × 10<sup>4</sup> m<sup>3</sup>/d, and then is quickly reduced to the stable gas production rate. This indicates that the gas supply capacity of artificial fractures is insufficient, and the matrix rapidly begins to supply gas. Note that the daily gas production rate is 4.5 × 10<sup>4</sup> m<sup>3</sup>/d, and the cumulative gas production is 0.87 × 10<sup>8</sup> m<sup>3</sup> at the end of production. When the formation thickness is 50 m, it is large enough to be a sufficient gas source. Specifically, the initial stages of production can reach 41.2 × 10<sup>4</sup> m<sup>3</sup>/d, and the stable gas production in the final stage of production can reach 13.9 × 10<sup>4</sup> m<sup>3</sup>/d; the total gas production reaches 2.66 × 10<sup>8</sup> m3.

**Figure 11.** The effect of different matrix permeability on productivity.

**Figure 12.** The effect of different formation thicknesses on productivity.

#### *4.4. Influence of Fracture Length on Gas Production*

The primary function of artificial fractures is to increase formation conductivity, and hence effectively improve the production capacity of reservoirs. Therefore, the length of the artificial fracture has a direct impact on both daily gas production and stable production capacity. Figure 13 shows its impact on the daily production rate and cumulative production when the half-lengths of artificial fractures are 50, 100, 150, 200, and 250 m. It can be seen that, at the initial stage of production, the production rate changes little, i.e. it is almost 46 × 10<sup>4</sup> m<sup>3</sup>/<sup>d</sup> for all fracture lengths. The high-productivity period varies proportionally to the length of the artificial fractures, which might be caused by the fact that the length of these fractures can promote high-speed turbulence within them. When the half-length of the artificial fracture is 50 and 250 m, the production rate at the end of the stage is 8.7 × 10<sup>4</sup> m<sup>3</sup>/<sup>d</sup> and 10.9 × 10<sup>4</sup> m<sup>3</sup>/d, respectively. As the artificial fracture length increases, the drainage radius and area increase; therefore, the production capacity of reservoirs and the controlled reserves of fractured horizontal wells improve. When the gas supply is sufficient, the daily gas production rate must be increased during the stage of stable production.

**Figure 13.** The effect of different artificial fracture lengths on productivity.

In contrast to the fracture length, the artificial fracture permeability primarily affects daily production in the early stage, but has little influence on the seepage of the matrix in the final stage. Figure 14 shows the impact on productivity due to fractured horizontal wells, which correspond to artificial fracture permeability of 200, 600, 1000, 1400, and 1800 × 10−<sup>3</sup> μm2. Here, we see that the artificial fracture permeability has a significant influence on initial gas production. When the artificial fracture permeability is 200 ×10−<sup>3</sup> μm2, the initial gas production rate is only 6.45 × 10<sup>4</sup> m<sup>3</sup>/d; however the daily production rate can reach 51.38 × 10<sup>4</sup> m<sup>3</sup>/<sup>d</sup> for an artificial fracture permeability of 1800 × <sup>10</sup>−3μm2, i.e., an increase of 696%. This demonstrates that the higher artificial fracture permeability induces a higher gas seepage rate and hence a higher gas production rate. Moreover, because of the high permeability, the gas source cannot be fully supplied, which leads to a sharp decrease in gas production. When the production time is less than 100 days, the daily production rate drops to 18 × 10<sup>4</sup> m<sup>3</sup>/d. During the stable production period, the artificial fracture permeability has little effect on productivity because gas is supplied by the matrix; and fractures act as channels but cannot supply gas themselves. When the artificial permeability is 200 and 1800 × 10−<sup>3</sup> μm2, the final gas production rate is 2.98 and 9.28 × 10<sup>4</sup> m<sup>3</sup>/d, respectively. The cumulative production is correspondingly 0.56 × 10<sup>8</sup> m<sup>3</sup> and 1.79 × 10<sup>8</sup> m3.

**Figure 14.** The effect of different artificial fracture permeability on productivity.

We attribute the influence of pressure on productivity as primarily due to the effect of a pressure differential. When the bottom flow pressure is fixed, formation pressure is greater, production pressure difference is greater, the driving force of gas seepage is greater, and the gas production rate is likewise higher. When the formation pressure is fixed, the bottom flow pressure is greater, the production pressure differential is smaller, and the driving force of gas is smaller, making the gas production rate

smaller Figure 15 shows the daily and cumulative production rates when the formation pressure is 50, 55, 60, 65, and 70 MPa. The figure shows that the higher is the formation pressure, the greater is the initial production rate, in addition to the stable production rate at the final stage. When formation pressure is 50 MPa, the initial production rate is 31.2 × 10<sup>4</sup> m<sup>3</sup>/d. The gas production rate and cumulative gas production volume are 6.18 × 10<sup>4</sup> m<sup>3</sup>/<sup>d</sup> and 1.19 × 10<sup>8</sup> m3, respectively. When the formation pressure is 70 MPa, the initial production considerably increases to 137.1 × 10<sup>4</sup> m<sup>3</sup>/d. The production rate and cumulative production at the end are 26.3 × 10<sup>4</sup> m<sup>3</sup>/<sup>d</sup> and 5.15 × 10<sup>8</sup> m3, respectively.

**Figure 15.** The effect of different formation pressures on productivity.

Figure 16 shows the daily production rate and the cumulative production volume corresponding to different bottom flow pressures, i.e., 0, 10, 30, 40, and 50 MPa. The larger is the bottom hole flow pressure, the smaller are the initial and stable production rates. When the bottom hole flow pressure is zero (an ideal condition that is actually impossible), the initial production rate is 203 × 10<sup>4</sup> m<sup>3</sup>/d, and the production rate in the middle and final stage is 37 × 10<sup>4</sup> m<sup>3</sup>/d, while cumulative production volume at the end of production is 7.2 × 10<sup>8</sup> m3. When the bottom hole flow pressure is 40 MPa, the initial production rate is relatively small (64 × 104m3/d) at the initial stage, but production decreases quickly and enters into the low production period almost instantly. Note that the production rate in the middle and later stages is 12.4 × 10<sup>4</sup> m<sup>3</sup>/d, and the cumulative production volume at the end of production is 2.4 × 10<sup>8</sup> m3.

**Figure 16.** The effect of different WBHPs on productivity.
