Transient-Flow Modeling of Vertical Fractured Wells with Multiple Hydraulic Fractures in Stress-Sensitive Gas Reservoirs

Massive hydraulic fracturing of vertical wells has been extensively employed in the development of low-permeability gas reservoirs. The existence of multiple hydraulic fractures along a vertical well makes the pressure profile around the vertical well complex. This paper studies the pressure dependence of permeability to develop a seepage model of vertical fractured wells with multiple hydraulic fractures. Both transformed pseudo-pressure and perturbation techniques have been employed to linearize the proposed model. The superposition principle and a hybrid analytical-numerical method were used to obtain the bottom-hole pseudo-pressure solution. Type curves for pseudo-pressure are presented and identified. The effects of the relevant parameters (such as dimensionless permeability modulus, fracture conductivity coefficient, hydraulic-fracture length, angle between the two adjacent hydraulic fractures, the difference of the hydraulic-fracture lengths, and hydraulic-fracture number) on the type curve and the error caused by neglecting the stress sensitivity are discussed in detail. The proposed work can enrich the understanding of the influence of the stress sensitivity on the performance of a vertical fractured well with multiple hydraulic fractures and can be used to more accurately interpret and forecast the transient pressure.


Introduction
Gas flow in porous media has recently attracted much attention and stimulated great interest in the development of oil and gas reservoirs [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15]. In the past decades, many gas reservoirs with stress-sensitive permeability have been discovered and developed around the world. During the development process of stress-sensitive reservoirs, the continuously decreasing formation pressure leads to the decrease of the formation permeability. Therefore, it is critical to investigate the effect of the stress sensitivity of the permeability on the production. Recently, lots of models have been established to study the performance of various wells in stress-sensitive reservoirs [16][17][18].
On the other hand, in order to obtain economic benefit, hydraulic fracturing has been widely applied in the development of low-permeability gas reservoirs. A great variety of seepage models of vertical fractured wells have been proposed and the characteristic of pressure response has been studied in detail. Gringarten et al. [19] established a seepage model of a vertical fractured well with two symmetrical infinite-conductivity hydraulic fractures, which can be used to identify the linear flow. Later, by considering the effect of the fluid flow within the hydraulic fractures, a seepage model of a vertical fractured well with two symmetrical finite-conductivity hydraulic fractures was established [20,21], which can identify the bilinear flow. Some researchers pointed out the existence of two asymmetrical hydraulic fractures in practice and established some seepage models of vertical fractured wells with two asymmetrical hydraulic fractures [22][23][24][25][26][27]. In fact, massive hydraulic fracturing usually creates fractured network or several main hydraulic fractures near the wellbore instead of two symmetrical/asymmetrical hydraulic fractures [28][29][30][31]. Recently, seepage models of vertical fractured wells with multiple hydraulic fractures have received more attention. Shaoul et al. [32] proposed a numerical model of a vertical fractured well with multiple hydraulic fractures and investigated the pressure behavior of a vertical fractured well. Restrepo and Tiab [33] and Wang et al. [34] established semi-analytical models of a vertical fractured well with multiple infinite-conductivity hydraulic fractures. Ren and Guo [35] and Luo and Tang [36] proposed semi-analytical models of a vertical fractured well with multiple finite-conductivity hydraulic fractures. Compared with the pressure behavior of vertical fractured wells with two hydraulic fractures, the pressure response of vertical fractured wells with multiple hydraulic fractures is very complex. Although the effect of relevant parameters on the pressure response of vertical fractured wells with multiple hydraulic fractures has been studied, and the characteristics of the type of curves is not clearly understood. More important, little work has focused on vertical fractured wells with multiple hydraulic fractures in stress-sensitive gas reservoirs. In particular, the effect of relevant parameters on the error caused by neglecting the stress sensitivity has not been recognized clearly.
In this paper, we propose a seepage model of vertical fractured wells with multiple finiteconductivity hydraulic fractures. Type curves are identified and analyzed. The effects of relevant parameters on type curves and the error caused by neglecting the stress sensitivity are discussed in detail. This work can enrich the understanding of the influence of stress sensitivity on the performance of vertical fractured wells with multiple hydraulic fractures.

Physical Model
Complex geometry patterns of hydraulic fractures for vertical fractured wells can be described by the multiple-fracture model which is shown in Figure 1. Some assumptions of the present model are given as follows: (1) The gas reservoir is a homogeneous and laterally infinite formation with a constant thickness. The top and bottom boundaries of the gas reservoir are assumed to be impermeable.
(2) The finite-conductivity hydraulic fractures emanate from the wellbore in the horizontal plane (as shown in Figure 1b) and completely penetrate the formation in the vertical plane. The length of each hydraulic fracture may be different.
(3) The vertical fractured well produces only from the hydraulic fractures and the production rate of the vertical fractured well keeps constant.
(4) The permeabilities of both the reservoir and hydraulic fractures change with pressure. (5) The gas reservoir has a constant temperature and uniform initial pressure.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 2 of 24 two symmetrical infinite-conductivity hydraulic fractures, which can be used to identify the linear flow. Later, by considering the effect of the fluid flow within the hydraulic fractures, a seepage model of a vertical fractured well with two symmetrical finite-conductivity hydraulic fractures was established [20,21], which can identify the bilinear flow. Some researchers pointed out the existence of two asymmetrical hydraulic fractures in practice and established some seepage models of vertical fractured wells with two asymmetrical hydraulic fractures [22][23][24][25][26][27]. In fact, massive hydraulic fracturing usually creates fractured network or several main hydraulic fractures near the wellbore instead of two symmetrical/asymmetrical hydraulic fractures [28][29][30][31]. Recently, seepage models of vertical fractured wells with multiple hydraulic fractures have received more attention. Shaoul et al. [32] proposed a numerical model of a vertical fractured well with multiple hydraulic fractures and investigated the pressure behavior of a vertical fractured well. Restrepo and Tiab [33] and Wang et al. [34] established semi-analytical models of a vertical fractured well with multiple infinite-conductivity hydraulic fractures. Ren and Guo [35] and Luo and Tang [36] proposed semi-analytical models of a vertical fractured well with multiple finite-conductivity hydraulic fractures. Compared with the pressure behavior of vertical fractured wells with two hydraulic fractures, the pressure response of vertical fractured wells with multiple hydraulic fractures is very complex. Although the effect of relevant parameters on the pressure response of vertical fractured wells with multiple hydraulic fractures has been studied, and the characteristics of the type of curves is not clearly understood. More important, little work has focused on vertical fractured wells with multiple hydraulic fractures in stress-sensitive gas reservoirs. In particular, the effect of relevant parameters on the error caused by neglecting the stress sensitivity has not been recognized clearly.
In this paper, we propose a seepage model of vertical fractured wells with multiple finite-conductivity hydraulic fractures. Type curves are identified and analyzed. The effects of relevant parameters on type curves and the error caused by neglecting the stress sensitivity are discussed in detail. This work can enrich the understanding of the influence of stress sensitivity on the performance of vertical fractured wells with multiple hydraulic fractures.

Physical Model
Complex geometry patterns of hydraulic fractures for vertical fractured wells can be described by the multiple-fracture model which is shown in Figure 1. Some assumptions of the present model are given as follows: (1) The gas reservoir is a homogeneous and laterally infinite formation with a constant thickness. The top and bottom boundaries of the gas reservoir are assumed to be impermeable.
(2) The finite-conductivity hydraulic fractures emanate from the wellbore in the horizontal plane (as shown in Figure 1b) and completely penetrate the formation in the vertical plane. The length of each hydraulic fracture may be different.
(3) The vertical fractured well produces only from the hydraulic fractures and the production rate of the vertical fractured well keeps constant.
(4) The permeabilities of both the reservoir and hydraulic fractures change with pressure. (5) The gas reservoir has a constant temperature and uniform initial pressure.  c c p  in this paper.
With the assumption of the uniform initial pressure, the initial condition is expressed as

Flow Model for Stress-Sensitive Gas Reservoirs
The reservoir permeability changes with reservoir pressure and can be described as follows [37]: where ψ is the pseudo-pressure, which is defined as The governing equation of gas flow in stress-sensitive porous media is given as follows [37]: Both gas viscosity µ and total compressibility c t in Equation (3) are functions of pressure. In order to simplify the proposed model and obtain an analytical solution, many research papers [26,37,38] treat the values of µ and c t in Equation (3) as constants, which could make the simulation results be accurate enough for engineering requirements. Following the suggestion by Wang [37], the values of µ and c t are evaluated at the initial condition, i.e., µ = µ(p i ) and c t = c t (p i ) in this paper.
With the assumption of the uniform initial pressure, the initial condition is expressed as The gas reservoir is assumed to be laterally infinite, and thus, the outer boundary condition is The inner boundary condition for a point source can be written as Introducing the dimensionless variables listed in Table 1, Equations (3)-(6) are rewritten as

Nomenclature Definition
Dimensionless pseudo − pressure The point source model including Equations (7)-(10) shows strong non-linearity which can be alleviated by introducing the following expression [16]: The perturbation technique, Laplace transform, and superposition principle [39] were employed to obtain the pressure distribution in the gas reservoir with a vertical fractured well with multiple hydraulic fractures producing at a constant rate (see Appendix B): Appl. Sci. 2019, 9, 1359 5 of 23

Flow Model for Stress-Sensitive Hydraulic Fractures
The dimensionless model for gas flow within stress-sensitive hydraulic fractures is given as (see Appendix C) According to the relationship between Cartesian coordinate system (x i , y i ) and polar coordinate system (r, θ) (see Figure 3), Equations (13)- (16) can be rewritten in the polar coordinate system as follows:

Flow Model for Stress-Sensitive Hydraulic Fractures
The dimensionless model for gas flow within stress-sensitive hydraulic fractures is given as (16) According to the relationship between Cartesian coordinate system   , i i x y and polar coordinate system   , r  (see Figure 3), Equations (13)- (16) can be rewritten in the polar coordinate system as follows:  Integrating Equation (17) from 0 to x with respect to r D yields Substituting Equations (18) and (19) into Equation (21) yields Integrating Equation (22) from 0 to r D with respect to r D yields where

Coupled Discrete Model
The flow model for gas reservoirs can associate with the flow model for hydraulic fractures by the following relationship: Substituting Equations (11) and (A37) into Equation (25), one can derive The zero-order perturbation solution of the ξ D is accurate enough in practice, and thus we can obtain Combining Equations (26) and (27) yields that Substituting Equation (28) into Equation (23) leads to Taking the Laplace transform of Equation (29), one can obtain Substituting Equation (12) into Equation (30), one can derive Employing the Laplace transform of Equations (19) and (20), one can obtain The coupled model consists of Equations (31) and (32), which can be solved by the numerical discrete method [40]. Each hydraulic fracture is discretized into some segments, and the flow rate in Appl. Sci. 2019, 9, 1359 7 of 23 each segment is considered to stay the same at a certain time. The discrete schematic of a hydraulic fracture of a vertical fractured well is shown in Figure 4.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 7 of 24 each segment is considered to stay the same at a certain time. The discrete schematic of a hydraulic fracture of a vertical fractured well is shown in Figure 4. And then, the discrete forms of Equations (31) and (32) can be obtained as follows Equations (33) and (34) compose The unknowns are which can be obtained by solving the linear equations, and then the wellbore storage and the skin near the wellbore can be taken into account by the following formula [41] With the aid of the numerical inversion method [42], And then, the discrete forms of Equations (31) and (32) can be obtained as follows Equations (33) and (34) compose The unknowns are ξ wHD (s) and q fk,vD (s) (1 ≤ k ≤ n, 1 ≤ v ≤ N k ), which can be obtained by solving the linear equations, and then the wellbore storage and the skin near the wellbore can be taken into account by the following formula [41] With the aid of the numerical inversion method [42], where Appl. Sci. 2019, 9, 1359 8 of 23 Finally, the dimensionless bottom-hole pseudo-pressure of a vertical fractured well with constant production rate is obtained by

Model Validation and Application
Because there is no analytical inversion solution of ξ wD (s), it is difficult to discuss the accuracy of the Laplace transform by comparing the analytical inversion solution with the numerical inversion solution [43]. In order to validate the proposed model and show how this model was used in practice, the drawdown test data of a vertical fractured well with four fracture wings were collected from the published literature [44]. The drawdown test data in the literature [44] were generated by the reservoir simulator which is based on the implicit finite difference method. The vertical fractured well produces at a constant production rate of 0.1639 m 3 /s. Basic parameters of the vertical fractured well are listed in Table 2. The proposed model is employed to simulate the bottom-hole pressure under the constant-rate-production condition. The simulated bottom-hole pseudo-pressure and bottom-hole pressure are compared with the results published in the literature [44]. It is shown from Figure 5 that the simulated bottom-hole pseudo-pressure and bottom-hole pressure by the proposed model are in good agreement with the bottom-hole pressure data published in the literature [44], which validates the proposed model.
Furthermore, the proposed model can be used to predict the bottom-hole pressure under the constant-rate-production condition. Based on the proposed model, it was easy to obtain the bottom-hole pressure of the vertical fractured well at any given time. For example, when the vertical fractured well produces at a constant production rate of 0.1639 m 3 /s, the bottom-hole pressures are 41.132 MPa, 37.676 MPa, 31.484 MPa at t = 100, 1000, 10000 hours, respectively.

Type Curves and Flow Regimes
The bottom-hole pseudo-pressure of a vertical fractured well with multiple hydraulic fractures in stress-sensitive gas reservoirs was calculated by the proposed model. Type curves for the bottom-hole pseudo-pressure were plotted and the characteristics of the bottom-hole pseudo-pressure behavior was analyzed. The established model in this work was suitable for multiple hydraulic fractures with arbitrary fracture number, arbitrary fracture length, and arbitrary fracture orientation. In order to investigate the effect of the hydraulic-fracture distribution on the pseudo-pressure response, without loss of generality, we mainly focused on the vertical fractured well with four fracture wings which is shown in Figure 6. The first and second fracture wings were assumed to have the equal length   f X , and the third and fourth fracture wings were set to be the

Type Curves and Flow Regimes
The bottom-hole pseudo-pressure of a vertical fractured well with multiple hydraulic fractures in stress-sensitive gas reservoirs was calculated by the proposed model. Type curves for the bottom-hole pseudo-pressure were plotted and the characteristics of the bottom-hole pseudo-pressure behavior was analyzed. The established model in this work was suitable for multiple hydraulic fractures with arbitrary fracture number, arbitrary fracture length, and arbitrary fracture orientation. In order to investigate the effect of the hydraulic-fracture distribution on the pseudo-pressure response, without loss of generality, we mainly focused on the vertical fractured well with four fracture wings which is shown in Figure 6. The first and second fracture wings were assumed to have the equal length (X f ), and the third and fourth fracture wings were set to be the equal length (Y f ). Figure 7 shows the type curves for the bottom-hole pseudo-pressure of a vertical fractured well with multiple hydraulic fractures in stress-sensitive gas reservoirs. The type curves for the bottom-hole pseudo-pressure consist of the pseudo-pressure curve (x−axis: t D /C D , y−axis: ψ wD ) and the pseudo-pressure derivative curve (x−axis: t D /C D , y−axis: ψ wD · t D /C D ). It is observed that there are seven possible flow regions in the type curves as follows: Appl. Sci. 2019, 9, x FOR PEER REVIEW 10 of 24 Figure 6. Schematic of a vertical fractured well with four fracture wings. Figure 7 shows the type curves for the bottom-hole pseudo-pressure of a vertical fractured well with multiple hydraulic fractures in stress-sensitive gas reservoirs. The type curves for the bottom-hole pseudo-pressure consist of the pseudo-pressure curve ( x  axis: D D t C , y  axis: wD  ) and the pseudo-pressure derivative curve ( x  axis: . It is observed that there are seven possible flow regions in the type curves as follows: (1) Wellbore storage period (WSP): The produced gas comes from the storage gas in the wellbore, and the gas in the reservoir has not flowed into the wellbore. Both the pseudo-pressure curve and the pseudo-pressure derivative curve exhibit as straight lines with the slope being one.
(2) Transitional flow after wellbore storage period (TFAWSP): the gas in the reservoir begins to flow into the wellbore. This flow period is controlled by both the wellbore storage and skin near the wellbore. The pseudo-pressure derivative curve appears as a "hump".
(3) Bilinear flow period (BFP): Linear flows take place within the hydraulic fractures and perpendicular to the hydraulic-facture surfaces in the reservoir simultaneously. The pseudo-pressure derivative curve exhibits as a straight line with the slope being 1/4.
(4) Transitional flow after bilinear flow period (TFABFP): The duration of this flow period is dependent on the difference of the hydraulic-fracture lengths. The shorter the hydraulic-fracture length is, the earlier the BFP ends. Therefore, this flow period will occur when the BFP of the shorter hydraulic fracture ends and the longer hydraulic fracture still lies in the BFP. The    (1) Wellbore storage period (WSP): The produced gas comes from the storage gas in the wellbore, and the gas in the reservoir has not flowed into the wellbore. Both the pseudo-pressure curve and the pseudo-pressure derivative curve exhibit as straight lines with the slope being one.
(2) Transitional flow after wellbore storage period (TFAWSP): the gas in the reservoir begins to flow into the wellbore. This flow period is controlled by both the wellbore storage and skin near the wellbore. The pseudo-pressure derivative curve appears as a "hump".
(3) Bilinear flow period (BFP): Linear flows take place within the hydraulic fractures and perpendicular to the hydraulic-facture surfaces in the reservoir simultaneously. The pseudo-pressure derivative curve exhibits as a straight line with the slope being 1/4.
(4) Transitional flow after bilinear flow period (TFABFP): The duration of this flow period is dependent on the difference of the hydraulic-fracture lengths. The shorter the hydraulic-fracture length is, the earlier the BFP ends. Therefore, this flow period will occur when the BFP of the shorter hydraulic fracture ends and the longer hydraulic fracture still lies in the BFP. The (1) Wellbore storage period (WSP): The produced gas comes from the storage gas in the wellbore, and the gas in the reservoir has not flowed into the wellbore. Both the pseudo-pressure curve and the pseudo-pressure derivative curve exhibit as straight lines with the slope being one.
(2) Transitional flow after wellbore storage period (TFAWSP): the gas in the reservoir begins to flow into the wellbore. This flow period is controlled by both the wellbore storage and skin near the wellbore. The pseudo-pressure derivative curve appears as a "hump".
(3) Bilinear flow period (BFP): Linear flows take place within the hydraulic fractures and perpendicular to the hydraulic-facture surfaces in the reservoir simultaneously. The pseudo-pressure derivative curve exhibits as a straight line with the slope being 1/4.
(4) Transitional flow after bilinear flow period (TFABFP): The duration of this flow period is dependent on the difference of the hydraulic-fracture lengths. The shorter the hydraulic-fracture length is, the earlier the BFP ends. Therefore, this flow period will occur when the BFP of the shorter hydraulic fracture ends and the longer hydraulic fracture still lies in the BFP. The pseudo-pressure derivative curve appears as a straight line with the slope in the range of 1/4-1/2.
(5) Linear flow period (LFP): When the BFP of each hydraulic fracture ends, the LFP will appear. In this period, the linear flow perpendicular to the hydraulic-facture surfaces in the reservoir dominates in the area near the wellbore. The pseudo-pressure derivative curve appears as a straight line with the slope being 1/2. (6) Transitional flow after linear flow period (TFALFP): When the pressure wave continues to travel in the reservoir, the interference between the pressure waves from different hydraulic fractures will take place, which makes the slope of the pseudo-pressure derivative curve be greater than 1/2. It should be noted that the stress sensitivity begins to obviously affect the type curves in this period. As the magnitude of the γ D increases, the pseudo-pressure and its derivative curves are shifted up.
(7) Pseudo-radial flow period (PRFP): Compared with the previous flow period (i.e., the TFALFP), much slower growth of the pseudo-pressure drop was observed. The pseudo-pressure derivative without the effect of the stress sensitivity (i.e., γ D = 0) keeps a value of 0.5, while the magnitude of the pseudo-pressure derivative with the effect of the stress sensitivity (i.e., γ D > 0) increases with the time. Furthermore, the pseudo-pressure and its derivative increase with increasing the value of the γ D at a fixed time. Therefore, if the stress sensitivity of the gas reservoir is stronger, the larger pressure drop is needed to remain the constant rate of produced well.
In order to quantify the effect of the stress sensitivity on the pseudo-pressure, the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity is introduced as follows where ψ wDnl and ψ wDl are the bottom-hole pseudo-pressures with and without the effect of the stress sensitivity, respectively. Figure 8 shows the effect of the dimensionless permeability modulus (γ D ) on the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity. With the increase of the time, the relative difference δ first increases very slowly and finally increases rapidly. Considering the flow regimes, it was found that the impact of the stress sensitivity on the pseudo-pressure was negligibly small during and before the LFP, while after the LFP, this impact became more and more obvious with the increase of the time. Furthermore, the δ became larger with increasing the magnitude of the γ D at a fixed time.
appear. In this period, the linear flow perpendicular to the hydraulic-facture surfaces in the reservoir dominates in the area near the wellbore. The pseudo-pressure derivative curve appears as a straight line with the slope being 1/2.
(6) Transitional flow after linear flow period (TFALFP): When the pressure wave continues to travel in the reservoir, the interference between the pressure waves from different hydraulic fractures will take place, which makes the slope of the pseudo-pressure derivative curve be greater than 1/2. It should be noted that the stress sensitivity begins to obviously affect the type curves in this period. As the magnitude of the D  increases, the pseudo-pressure and its derivative curves are shifted up.
(7) Pseudo-radial flow period (PRFP): Compared with the previous flow period (i.e., the TFALFP), much slower growth of the pseudo-pressure drop was observed. The pseudo-pressure derivative without the effect of the stress sensitivity (i.e., D 0   ) keeps a value of 0.5, while the magnitude of the pseudo-pressure derivative with the effect of the stress sensitivity (i.e., D 0   ) increases with the time. Furthermore, the pseudo-pressure and its derivative increase with increasing the value of the D  at a fixed time. Therefore, if the stress sensitivity of the gas reservoir is stronger, the larger pressure drop is needed to remain the constant rate of produced well.
In order to quantify the effect of the stress sensitivity on the pseudo-pressure, the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity is introduced as follows

Sensitivity Analysis
Besides the dimensionless permeability modulus (γ D ), other parameters may have effects on the type curve and the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity. Therefore, it is critical to conduct the sensitivity analysis of relevant parameters. Figures 9 and 10 show the effects of the fracture conductivity coefficient (C fD ) on the type curve and the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity, respectively. It is obvious that the impact of the C fD on the type curve takes place in the early periods (i.e., from the TFAWSP to the LFP), where the values of the pseudo-pressure and its derivative increase with the decrease of the magnitude of the C fD . Furthermore, it was found from Figure 9 that as the value of the C fD increases, the duration of the BFP becomes shorter and the start time of the LFP becomes earlier. As shown in Figure 10, the δ was little affected by the C fD in all periods, and this was because the main impact of the C fD on the pseudo-pressure occurred in the early periods, but the effect of the stress sensitivity on the pseudo-pressure was negligibly small in the early periods.

Sensitivity Analysis
Besides the dimensionless permeability modulus   D  , other parameters may have effects on the type curve and the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity. Therefore, it is critical to conduct the sensitivity analysis of relevant parameters. Figures 9 and 10 show the effects of the fracture conductivity coefficient   fD C on the type curve and the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity, respectively. It is obvious that the impact of the fD C on the type curve takes place in the early periods (i.e., from the TFAWSP to the LFP), where the values of the pseudo-pressure and its derivative increase with the decrease of the magnitude of the fD C . Furthermore, it was found from Figure 9 that as the value of the fD C increases, the duration of the BFP becomes shorter and the start time of the LFP becomes earlier. As shown in Figure 10, the  was little affected by the fD C in all periods, and this was because the main impact of the fD C on the pseudo-pressure occurred in the early periods, but the effect of the stress sensitivity on the pseudo-pressure was negligibly small in the early periods.     Figures 11 and 12 show the effects of the length of the hydraulic fractures (L f ) on the type curve and the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity, respectively. It is seen from Figure 11 that the L f has an important effect on the magnitudes of the pseudo-pressure and its derivative in the intermediate and late flow periods (i.e., from the BFP to the PRFP). Decreasing the L f reduced the duration of the BFP and resulted in the increase of the pseudo-pressure and its derivative during and after the BFP. Furthermore, it was found from Figure 12 that the δ increased with decreasing the magnitude of the L f in the late flow period, indicating that the pseudo-pressure obtained by the conventional model without the effect of the stress sensitivity will result in a much bigger error when the length of the hydraulic fractures becomes shorter.  Figures 11 and 12 show the effects of the length of the hydraulic fractures   f L on the type curve and the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity, respectively. It is seen from Figure 11 that the f L has an important effect on the magnitudes of the pseudo-pressure and its derivative in the intermediate and late flow periods (i.e., from the BFP to the PRFP). Decreasing the f L reduced the duration of the BFP and resulted in the increase of the pseudo-pressure and its derivative during and after the BFP. Furthermore, it was found from Figure 12 that the  increased with decreasing the magnitude of the f L in the late flow period, indicating that the pseudo-pressure obtained by the conventional model without the effect of the stress sensitivity will result in a much bigger error when the length of the hydraulic fractures becomes shorter. Figure 11. Type curves for the bottom-hole pseudo-pressure of a vertical fractured well with multiple hydraulic fractures in stress-sensitive gas reservoirs with different values of the f L ( 4 n  , Appl. Sci. 2019, 9, x FOR PEER REVIEW 14 of 24 Figure 12.
Effect of the f L on the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity ( . Figures 13 and 14 show the effects of the angle between the two adjacent hydraulic fractures    on the type curve and the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity, respectively. As shown in Figure 13, with the decrease of the  , the start time of the TFALFP becomes earlier, and the magnitudes of the pseudo-pressure and its derivative in the TFALFP increase. This is because decreasing the  results in the earlier and stronger interference between the hydraulic fractures, which makes the pressure drop increase to remain the constant rate of the produced well. Whereas the impact of the  on the type curve is relatively unobvious in the PRFP. Furthermore, it is seen from Figure 14 that the  increases with decreasing the magnitude of the  during and after the TFALFP, indicating that the effect of the stress sensitivity on the pseudo-pressure becomes greater when the  decreases.  Figures 13 and 14 show the effects of the angle between the two adjacent hydraulic fractures (θ) on the type curve and the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity, respectively. As shown in Figure 13, with the decrease of the θ, the start time of the TFALFP becomes earlier, and the magnitudes of the pseudo-pressure and its derivative in the TFALFP increase. This is because decreasing the θ results in the earlier and stronger interference between the hydraulic fractures, which makes the pressure drop increase to remain the constant rate of the produced well. Whereas the impact of the θ on the type curve is relatively unobvious in the PRFP. Furthermore, it is seen from Figure 14 that the δ increases with decreasing the magnitude of the θ during and after the TFALFP, indicating that the effect of the stress sensitivity on the pseudo-pressure becomes greater when the θ decreases.    on the type curve and the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity, respectively. As shown in Figure 13, with the decrease of the  , the start time of the TFALFP becomes earlier, and the magnitudes of the pseudo-pressure and its derivative in the TFALFP increase. This is because decreasing the  results in the earlier and stronger interference between the hydraulic fractures, which makes the pressure drop increase to remain the constant rate of the produced well. Whereas the impact of the  on the type curve is relatively unobvious in the PRFP. Furthermore, it is seen from Figure 14 that the  increases with decreasing the magnitude of the  during and after the TFALFP, indicating that the effect of the stress sensitivity on the pseudo-pressure becomes greater when the  decreases.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 15 of 24 Figure 14. Effect of the  on the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity ( 4 n  , f 1 Introducing the coefficient of to quantify the difference of the hydraulic-fracture lengths, the effects of  on the type curve, and the relative differences between the pseudo-pressures with and without the effect of the stress sensitivity are shown in Figures 15  and 16, respectively. It is obvious that the difference of the hydraulic-fracture lengths becomes larger as the magnitude of the  increases. The range of the  is from zero to one. As shown in Figure 15, with the decrease of the  , the end of the BFP takes place later, the duration of the TFABFP becomes shorter, and the start time of the LFP occurs earlier. The reason for the existence of the TFABFP is that the difference of the hydraulic-fracture lengths leads to the inconsistency of the end times for the BFP for each hydraulic fracture. If  is equal to zero, the TFABFP may not even appear. Furthermore, the difference of the hydraulic-fracture lengths makes the pseudo-pressure and its derivative increase in the TFABFP. It is interesting to find in Figure 15 that during the TFALFP, the pseudo-pressure and its derivative increase as the  decreases. This may be because that for the distribution of the hydraulic fractures shown in Figure 6, the interference  Figure 14. Effect of the θ on the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity (n = 4, Introducing the coefficient of β = |X f − Y f |/(X f + Y f ) to quantify the difference of the hydraulic-fracture lengths, the effects of β on the type curve, and the relative differences between the pseudo-pressures with and without the effect of the stress sensitivity are shown in Figures 15 and 16, respectively. It is obvious that the difference of the hydraulic-fracture lengths becomes larger as the magnitude of the β increases. The range of the β is from zero to one. As shown in Figure 15, with the decrease of the β, the end of the BFP takes place later, the duration of the TFABFP becomes shorter, and the start time of the LFP occurs earlier. The reason for the existence of the TFABFP is that the difference of the hydraulic-fracture lengths leads to the inconsistency of the end times for the BFP for each hydraulic fracture. If β is equal to zero, the TFABFP may not even appear. Furthermore, the difference of the hydraulic-fracture lengths makes the pseudo-pressure and its derivative increase in the TFABFP. It is interesting to find in Figure 15 that during the TFALFP, the pseudo-pressure and its derivative increase as the β decreases. This may be because that for the distribution of the hydraulic fractures shown in Figure 6, the interference between hydraulic fractures enhances with the decrease of the β in the TFALFP. As shown in Figure 16, the β has an effect on the error caused by neglecting the stress sensitivity during and after the TFALFP. Decreasing the β results in the increase of the error. Figure 15, with the decrease of the  , the end of the BFP takes place later, the duration of the TFABFP becomes shorter, and the start time of the LFP occurs earlier. The reason for the existence of the TFABFP is that the difference of the hydraulic-fracture lengths leads to the inconsistency of the end times for the BFP for each hydraulic fracture. If  is equal to zero, the TFABFP may not even appear. Furthermore, the difference of the hydraulic-fracture lengths makes the pseudo-pressure and its derivative increase in the TFABFP. It is interesting to find in Figure 15 that during the TFALFP, the pseudo-pressure and its derivative increase as the  decreases. This may be because that for the distribution of the hydraulic fractures shown in Figure 6, the interference between hydraulic fractures enhances with the decrease of the  in the TFALFP. As shown in Figure 16, the  has an effect on the error caused by neglecting the stress sensitivity during and after the TFALFP. Decreasing the  results in the increase of the error.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 16 of 24 Figure 15. Type curves for the bottom-hole pseudo-pressure of a vertical fractured well with multiple hydraulic fractures in stress-sensitive gas reservoirs with different values of the  ( 4 n  , . Figures 17 and 18 show the effects of the hydraulic-fracture number   n on the type curve and the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity, respectively. The distribution of the hydraulic fractures was assumed to be uniform (i.e., the angle between the two arbitrary adjacent hydraulic fractures was equal). As shown in Figure 17, it is obvious that the n affects the pseudo-pressure behavior in all periods except for the WSP. Increasing the n leads to the decrease of the pseudo-pressure and its derivative. Due to the increase of the n , the angle between the two adjacent hydraulic fractures was reduced, and thus the interference between hydraulic fractures occurred earlier and the effect of this interference on the pseudo-pressure was more significant. Figure 18 indicates that a smaller value of the n leads to a bigger error caused by neglecting the stress sensitivity during and after the TFALFP.  Figure 16. Effect of the β on the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity (n = 4, Figures 17 and 18 show the effects of the hydraulic-fracture number (n) on the type curve and the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity, respectively. The distribution of the hydraulic fractures was assumed to be uniform (i.e., the angle between the two arbitrary adjacent hydraulic fractures was equal). As shown in Figure 17, it is obvious that the n affects the pseudo-pressure behavior in all periods except for the WSP. Increasing the n leads to the decrease of the pseudo-pressure and its derivative. Due to the increase of the n, the angle between the two adjacent hydraulic fractures was reduced, and thus the interference between hydraulic fractures occurred earlier and the effect of this interference on the pseudo-pressure was more significant. Figure 18 indicates that a smaller value of the n leads to a bigger error caused by neglecting the stress sensitivity during and after the TFALFP. and the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity, respectively. The distribution of the hydraulic fractures was assumed to be uniform (i.e., the angle between the two arbitrary adjacent hydraulic fractures was equal). As shown in Figure 17, it is obvious that the n affects the pseudo-pressure behavior in all periods except for the WSP. Increasing the n leads to the decrease of the pseudo-pressure and its derivative. Due to the increase of the n , the angle between the two adjacent hydraulic fractures was reduced, and thus the interference between hydraulic fractures occurred earlier and the effect of this interference on the pseudo-pressure was more significant. Figure 18 indicates that a smaller value of the n leads to a bigger error caused by neglecting the stress sensitivity during and after the TFALFP.

Conclusions
We have established a non-linear model of a vertical fractured well with multiple hydraulic fractures in gas reservoirs by incorporating the pressure-dependent permeability. Both transformed pseudo-pressure and perturbation techniques were employed to linearize the proposed model. Superposition principle and numerical discrete methods were used to obtain the semi-analytical solution. Type curves for pseudo-pressure were plotted and discussed in detail. The influence of relevant parameters on the type curve and the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity were analyzed. Some main conclusions are listed as follows: (1) The type curve of a vertical fractured well with multiple hydraulic fractures was identified by seven possible flow regions. Compared with the conventional model of a vertical fractured well with two symmetrical hydraulic-fracture wings, two transitional flow regimes (i.e., the TFABFP and TFALFP) were observed, which were caused by the difference of the hydraulic-fracture lengths and the interference between hydraulic fractures, respectively.
(2) As the time increased, the pressure drops increased, which made the permeability decrease. The impact of the stress sensitivity on the pseudo-pressure increased with the increase of  Figure 18. Effect of n on the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity (L f1 = L f2 = X f = 80 m, L f3 = L f4 = Y f = 80 m, θ = 360 • /n, C fD = 5 × 10 4 , S f = 10 −3 , C D = 0.1, γ D = 0.1).

Conclusions
We have established a non-linear model of a vertical fractured well with multiple hydraulic fractures in gas reservoirs by incorporating the pressure-dependent permeability. Both transformed pseudo-pressure and perturbation techniques were employed to linearize the proposed model. Superposition principle and numerical discrete methods were used to obtain the semi-analytical solution. Type curves for pseudo-pressure were plotted and discussed in detail. The influence of relevant parameters on the type curve and the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity were analyzed. Some main conclusions are listed as follows: (1) The type curve of a vertical fractured well with multiple hydraulic fractures was identified by seven possible flow regions. Compared with the conventional model of a vertical fractured well with two symmetrical hydraulic-fracture wings, two transitional flow regimes (i.e., the TFABFP and TFALFP) were observed, which were caused by the difference of the hydraulic-fracture lengths and the interference between hydraulic fractures, respectively.
(2) As the time increased, the pressure drops increased, which made the permeability decrease. The impact of the stress sensitivity on the pseudo-pressure increased with the increase of the time. Furthermore, the influence of the stress sensitivity on the pseudo-pressure was negligibly small in the early and intermediate flow periods but becomes very significant in the late flow period.
(3) Some relevant parameters, such as dimensionless permeability modulus, fracture conductivity coefficient, hydraulic-fracture length, angle between the two adjacent hydraulic fractures, the difference of the hydraulic-fracture lengths, and hydraulic-fracture number, not only affected the type curve, but also have an influence of the error caused by neglecting the stress sensitivity.

Conflicts of Interest:
The authors declare no conflict of interest.