# A Mathematical Pressure Transient Analysis Model for Multiple Fractured Horizontal Wells in Shale Gas Reservoirs

Authors: Yan Zeng; Qing Wang; Zhengfu Ning; Hongliang Sun
Journal: Geofluids
 (2018)

---

## Abstract

Multistage fractured horizontal wells (MFHWs) have become the main technology for shale gas exploration. However, the existing models have neglected the percolation mechanism in nanopores of organic matter and failed to consider the differences among the reservoir properties in different areas. On that account, in this study, a modified apparent permeability model was proposed describing gas flow in shale gas reservoirs by integrating bulk gas flow in nanopores and gas desorption from nanopores. The apparent permeability was introduced into the macroseepage model to establish a dynamic pressure analysis model for MFHWs dual-porosity formations. The Laplace transformation and the regular perturbation method were used to obtain an analytical solution. The influences of fracture half-length, fracture permeability, Langmuir volume, matrix radius, matrix permeability, and induced fracture permeability on pressure and production were discussed. Results show that fracture half-length, fracture permeability, and induced fracture permeability exert a significant influence on production. A larger Langmuir volume results in a smaller pressure and pressure derivative. An increase in matrix permeability increases the production rate. Besides, this model fits the actual field data relatively well. It has a reliable theoretical foundation and can preferably describe the dynamic changes of pressure in the exploration process.


---

## Body

## 1. Introduction

Shale gas is known as a key resource to meet the increasing world energy demand because of its rich reserves and extensive distribution [1]. To be more specific, many pieces of research have shown that shale gas reservoirs are characterized by extremely low permeability values. Consequently, hydraulic fracturing simulations and horizontal drilling are developed for commercial exploitation. In shale gas reservoirs, the complex pore structure, macro/microfracture network distribution, and microscale reservoir properties make the gas flow characteristics become more complicated than those in conventional reservoirs. Gas flow in organic matter, natural fractures, secondary fractures, and hydraulic fractures is controlled by different mechanisms [2]. To be more specific, 20%–80% of shale gas exists in the adsorbed state on the surface of organic particulates. Desorption of adsorbed gas in organic matter affects gas transport mechanisms significantly [3, 4]. Understanding the pressure- and rate-transient behaviors of multistage fractured horizontal wells in shale gas reservoirs is of great importance for production forecasting, well placement, and configuration optimization [5].In recent studies [6–8], gas transport mechanisms in organic matter can be described by continuum flow, slip flow viscous flow, the Knudsen diffusion, and transition flow. Wang and Li [9] made a comparison of collision frequency in gas-gas and gas-solid interactions. They concluded that the collisions between gas molecules and the solid wall account for a large percentage of the total collisions. Therefore, the Knudsen diffusion should be considered when studying gas transport in nanopores. Darabi et al. [10] also reported that, under typical shale gas reservoir conditions, the Knudsen diffusion dominates gas transport and the contribution to cumulative production can reach 20%. Adsorbed gas desorption from the surface of organic matter also plays a significant role in gas transport [11, 12]. Swami et al. [13] developed a modified model to investigate the contribution of adsorbed gas to the final recovery. They found that desorption of adsorbed gas can increase the pore diameter, reduce tortuosity, and cause extra slippage at the solid boundary. Yu and Sepehrnoori [14] analyzed the shale gas development history for a period of 30 years in North America and concluded that different shale samples demonstrate quite different adsorption capacities.Mayerhofer et al. [15] studied the difference between conventional double wings symmetric fractures and hydraulic fractures of Barnet shale and were the first to propose the simulated reservoir volume (SRV) conception. In SRV, multistage hydraulic fracturing connects existing natural fractures, which generates a large fracture network [16]. Kucuk and Sawyer [17] firstly developed an analytical model to investigate transient pressure in shale gas reservoirs. However, the model ignored diffusion flow and the desorption effect. Lee and Brockenbrough [18] analyzed production data coming from vertical wells and developed an analytical solution to describe transient pressure. Ozkan et al. [19] proposed a dual-mechanism, that is to say, a dual-porosity horizontal well model for shale gas reservoirs, which includes hydraulic fractures, as well as the inner region and the outer region. The model can adopt the diffusive flow mechanism but fail to take adsorbed gas into consideration. Stalgorova and Mattar [20–22] extended the trilinear flow model to a five-region model for MFHWs in homogeneous shale gas reservoirs. In order to investigate the effect of natural fractures on matrix permeability, Apaydin et al. [23] combined composite blocks into an analytical trilinear flow model. Zhao et al. [24, 25] analyzed the pressure response and production performance of MFHWs in shale gas reservoirs and they introduced the source function theory to characterize SRV. Liu et al. [26] discussed the desorption and diffusion effects in the analytical model and studied the productivity-decline law. Wang [27] considered the stress sensitive effect of natural fractures and hydraulic fracture angles in the semianalytical model. His research implied that the stress sensitive effect has a significant impact on bottom hole flowing pressure.In recent years, many general semianalytical models, considering microscopic seepage mechanisms, have been proposed with the rapid development of computer technology [24, 25, 28–30]. These models can describe predominant flow regimes for shale well but fail to consider all of the typical characteristics and possible situations in shale gas reservoirs. In the present study, a novel apparent permeability model was put forward by taking into consideration viscous flow, slip flow, transition flow, gas adsorption/desorption [31], and the poromechanical response [32, 33]. The novelty of this study concentrates on the improvement of the semianalytical model to conduct pressure- and rate-transient analysis in SRV, by considering microscale percolation mechanisms and heterogeneity. The rest of this study can be described in the following way. In Section 2, the apparent permeability model for gas transport in shale nanopores is introduced. The semianalytical model for MFHWs in a heterogeneous shale gas reservoir, by considering microscale percolation mechanisms in a dual-porosity formation, is introduced in Section 3. Log-log dimension pressure, dimension pressure derivative, and production type curves are plotted firstly, and then a sensitivity analysis of formation parameters is conducted in Section 4. An actual field case is studied in Section 5. Finally, some conclusions are drawn in Section 6.


## 2. Apparent Permeability Model for Gas Transport in Shale Nanopores

### 2.1. Viscous Flow

The total shale gas flux is composed of a viscous flow flux and a Knudsen diffusion flux component. Viscous flow is strongly prominent when the short molecule free path is far less than the pore diameter. In other words, the transport mechanism is governed by viscous flow in macropores (where the pore diameter > 50 nm). The Hagen-Poiseuille equation can describe the molar flux [34] as stated below:(1)JV=-ϕτr2P8μRTdPdl,where JV is the viscous flow flux, mol/(m2∗s); ϕ is porosity, dimensionless; τ is tortuosity, dimensionless; r is the pore radius, m; μ is viscosity, Pa∗s; R is the universal gas constant, J/(mol∗K); T is temperature, K; P is pressure, Pa; and l is the shale gas transport distance, m.Gas viscosity gradually deviates from the traditional viscosity definition for high Knudsen numbers. Karniadakis et al. [35] modified the viscosity definition in the Knudsen layer by considering the rarefaction effect.(2)μeff=11+αKnμ,where μeff is the viscosity considering rarefaction effect, Pa∗s; α is the rarefaction coefficient, dimensionless; and Kn is the Knudsen number, dimensionless.According to the Darcy equation, the apparent permeability of viscous flow can be expressed as [36]:(3)Kv=Jv-effμVstddP/dl=1+αKnϕτr2PVstd8RT=1+αKnK∞,where Kv is the apparent permeability, taking into account the effect of rarefaction, m2; Vstd is the gas mole volume at standard conditions, 22.414 ∗ 10−3 m3/mol; and K∞ is the permeability of viscous flow, m2.


### 2.2. The Knudsen Diffusion

When there is shale gas transport through micropores under low pressure, the molecule free path is long and equal to the pore diameter. The Knudsen diffusion is prominent because collisions between molecules and the micropores wall are more frequent than the intermolecular collisions. The flux of the Knudsen diffusion can be expressed as [34](4)JK=-ϕτMDKRT∂C∂l,where JK is the Knudsen diffusion flux, mol/(m2∗s); DK is the Knudsen diffusion coefficient, m2/s; and C is the gas concentration, m3/mol. The Knudsen diffusion coefficient is obtained as(5)DK=23r8RTπM,where M is the gas molar mass, kg/mol. The rough wall has a significant effect on the Knudsen diffusion. Darabi et al. [10] proposed a function to describe the roughness effect as follows:(6)Deff=δDf-2DK=λ¯2rDf-2DK,where Deff is the Knudsen diffusion coefficient which takes into consideration the roughness effect, m2/s; λ¯ is the mean free path of gas molecules, m; and Df is the fractal dimension of the rough wall, dimensionless.By combining (4) and (6), the apparent permeability of the Knudsen diffusion can be obtained and written as follows:(7)KK=JK-effμVstdMdC/dl=23ϕτrμVstdMλ¯2rDf-28MπRT.


### 2.3. Weight Factor

In the actual shale gas reservoirs, different transport mechanisms exist at the same time. Therefore, the total shale gas flux should be a weighted summation of the viscous flow flux and the Knudsen diffusion flux, based on their different contributions. The weight factors of the viscous flow and the Knudsen diffusion flow are defined as the ratios of intermolecular collisions and molecular/pore-wall collisions to total collisions, respectively. Based on the definitions for a gas molecule free path [11] and collision numbers [37], the weight factors of viscous flow and Knudsen diffusion can be approximated as follows:(8)wV=w=11+Kn/2wK=1-w=11+2/Kn.Therefore, the total apparent permeability in nanopores is(9)K=wKV+1-wKK=w1+αKnK∞+1-w23ϕτrμVstdMλ¯2rDf-28MπRT.With a decrease in pressure, the gas molecule free path increases immediately. The weight factor of the Knudsen diffusion also increases. The piece of research carried out by Wu et al. [38] showed that, under the 10 nm radius condition, viscous flow dominates when the pressure is larger than 4.02 ∗ 105 Pa; otherwise, the Knudsen diffusion dominates.


### 2.4. Sorption-Induced Swelling Response

Organic matter is described by weak strength and strong sensitivity to stress change. With an increase in reservoir pressure, adsorbed gas begins to desorb. Gas desorption results in shrinkage of the organic matrix and to an increased effective hydraulic diameter [39]. According to the solid deformation theory [40] and the Langmuir isotherm equation [41], the relationship between the degree of solid deformation and the reservoir pressure can be written as follows:(10)Δε=VLρsRTEVstdln⁡PL+Pini-ln⁡PL+P+cPPini-P,where Δε is the solid deformation degree, dimensionless; VL is the Langmuir volume, m3/kg; ρs is the density of the shale matrix, kg/m3; E is the shale matrix Young modulus, Pa; PL is the Langmuir pressure, Pa; Pini is the initial pressure, Pa; and cP is the rock compressibility, Pa−1.Meanwhile, based on Seidle’s model [42], the approximation relation between the porosity and the solid deformation degree is given by(11)ϕϕini=1+1+2ϕiniΔε.Assume that the pore volumes are proportional to the gas flow channels. In accordance with the capillary model, the effective hydraulic diameter can be obtained as(12)reff=r1+1+2ϕiniVLρsRTEVstdlnPL+PiniPL+P+cPPini-P.


### 2.5. Gas Adsorption/Desorption

Because of the large surface area and oil-wet characteristic, nanopores have a strong adsorbed gas capacity. The mass balance equation, which considers the adsorbed gas, can be expressed as(13)ρcgϕ1+1-ϕϕ∂V∂ρ∂P∂t=1r∂∂rrKμρ∂P∂r,where ρ is the shale gas density, kg/m3; cg is the gas compressibility, Pa−1; and V is the mass of adsorbed gas per unit volume of shale, kg/m3. The apparent permeability, which considers desorption of absorbed gas, can be expressed as(14)Kapp=1+1-ϕϕ∂V∂ρK.Cui et al. [12] defined the effective adsorption porosity as(15)ϕad=1-ϕ∂V∂ρ.According to the Langmuir isotherm equation, the effective adsorption porosity can be expressed as(16)ϕad=ρsMVstd1-ϕρcgVmPLPL+P2.Finally, the apparent permeability model for gas transport in shale nanopores can be expressed as(17)Kapp=1+ϕadϕKϕ=1+1+2ϕiniΔεϕiniϕad=ρsMVstd1-ϕρcgVmPLPL+P2Δε=VLρsRTEVstdln⁡PL+PiniPL+P+cPPini-PK=w1+αKnϕτreff2PVstd8RT+1-w23ϕτreffμVstdλ¯2rDf-28MπRTreff=r1+1+2ϕiniΔε.


## 2.1. Viscous Flow

The total shale gas flux is composed of a viscous flow flux and a Knudsen diffusion flux component. Viscous flow is strongly prominent when the short molecule free path is far less than the pore diameter. In other words, the transport mechanism is governed by viscous flow in macropores (where the pore diameter > 50 nm). The Hagen-Poiseuille equation can describe the molar flux [34] as stated below:(1)JV=-ϕτr2P8μRTdPdl,where JV is the viscous flow flux, mol/(m2∗s); ϕ is porosity, dimensionless; τ is tortuosity, dimensionless; r is the pore radius, m; μ is viscosity, Pa∗s; R is the universal gas constant, J/(mol∗K); T is temperature, K; P is pressure, Pa; and l is the shale gas transport distance, m.Gas viscosity gradually deviates from the traditional viscosity definition for high Knudsen numbers. Karniadakis et al. [35] modified the viscosity definition in the Knudsen layer by considering the rarefaction effect.(2)μeff=11+αKnμ,where μeff is the viscosity considering rarefaction effect, Pa∗s; α is the rarefaction coefficient, dimensionless; and Kn is the Knudsen number, dimensionless.According to the Darcy equation, the apparent permeability of viscous flow can be expressed as [36]:(3)Kv=Jv-effμVstddP/dl=1+αKnϕτr2PVstd8RT=1+αKnK∞,where Kv is the apparent permeability, taking into account the effect of rarefaction, m2; Vstd is the gas mole volume at standard conditions, 22.414 ∗ 10−3 m3/mol; and K∞ is the permeability of viscous flow, m2.


## 2.2. The Knudsen Diffusion

When there is shale gas transport through micropores under low pressure, the molecule free path is long and equal to the pore diameter. The Knudsen diffusion is prominent because collisions between molecules and the micropores wall are more frequent than the intermolecular collisions. The flux of the Knudsen diffusion can be expressed as [34](4)JK=-ϕτMDKRT∂C∂l,where JK is the Knudsen diffusion flux, mol/(m2∗s); DK is the Knudsen diffusion coefficient, m2/s; and C is the gas concentration, m3/mol. The Knudsen diffusion coefficient is obtained as(5)DK=23r8RTπM,where M is the gas molar mass, kg/mol. The rough wall has a significant effect on the Knudsen diffusion. Darabi et al. [10] proposed a function to describe the roughness effect as follows:(6)Deff=δDf-2DK=λ¯2rDf-2DK,where Deff is the Knudsen diffusion coefficient which takes into consideration the roughness effect, m2/s; λ¯ is the mean free path of gas molecules, m; and Df is the fractal dimension of the rough wall, dimensionless.By combining (4) and (6), the apparent permeability of the Knudsen diffusion can be obtained and written as follows:(7)KK=JK-effμVstdMdC/dl=23ϕτrμVstdMλ¯2rDf-28MπRT.


## 2.3. Weight Factor

In the actual shale gas reservoirs, different transport mechanisms exist at the same time. Therefore, the total shale gas flux should be a weighted summation of the viscous flow flux and the Knudsen diffusion flux, based on their different contributions. The weight factors of the viscous flow and the Knudsen diffusion flow are defined as the ratios of intermolecular collisions and molecular/pore-wall collisions to total collisions, respectively. Based on the definitions for a gas molecule free path [11] and collision numbers [37], the weight factors of viscous flow and Knudsen diffusion can be approximated as follows:(8)wV=w=11+Kn/2wK=1-w=11+2/Kn.Therefore, the total apparent permeability in nanopores is(9)K=wKV+1-wKK=w1+αKnK∞+1-w23ϕτrμVstdMλ¯2rDf-28MπRT.With a decrease in pressure, the gas molecule free path increases immediately. The weight factor of the Knudsen diffusion also increases. The piece of research carried out by Wu et al. [38] showed that, under the 10 nm radius condition, viscous flow dominates when the pressure is larger than 4.02 ∗ 105 Pa; otherwise, the Knudsen diffusion dominates.


## 2.4. Sorption-Induced Swelling Response

Organic matter is described by weak strength and strong sensitivity to stress change. With an increase in reservoir pressure, adsorbed gas begins to desorb. Gas desorption results in shrinkage of the organic matrix and to an increased effective hydraulic diameter [39]. According to the solid deformation theory [40] and the Langmuir isotherm equation [41], the relationship between the degree of solid deformation and the reservoir pressure can be written as follows:(10)Δε=VLρsRTEVstdln⁡PL+Pini-ln⁡PL+P+cPPini-P,where Δε is the solid deformation degree, dimensionless; VL is the Langmuir volume, m3/kg; ρs is the density of the shale matrix, kg/m3; E is the shale matrix Young modulus, Pa; PL is the Langmuir pressure, Pa; Pini is the initial pressure, Pa; and cP is the rock compressibility, Pa−1.Meanwhile, based on Seidle’s model [42], the approximation relation between the porosity and the solid deformation degree is given by(11)ϕϕini=1+1+2ϕiniΔε.Assume that the pore volumes are proportional to the gas flow channels. In accordance with the capillary model, the effective hydraulic diameter can be obtained as(12)reff=r1+1+2ϕiniVLρsRTEVstdlnPL+PiniPL+P+cPPini-P.


## 2.5. Gas Adsorption/Desorption

Because of the large surface area and oil-wet characteristic, nanopores have a strong adsorbed gas capacity. The mass balance equation, which considers the adsorbed gas, can be expressed as(13)ρcgϕ1+1-ϕϕ∂V∂ρ∂P∂t=1r∂∂rrKμρ∂P∂r,where ρ is the shale gas density, kg/m3; cg is the gas compressibility, Pa−1; and V is the mass of adsorbed gas per unit volume of shale, kg/m3. The apparent permeability, which considers desorption of absorbed gas, can be expressed as(14)Kapp=1+1-ϕϕ∂V∂ρK.Cui et al. [12] defined the effective adsorption porosity as(15)ϕad=1-ϕ∂V∂ρ.According to the Langmuir isotherm equation, the effective adsorption porosity can be expressed as(16)ϕad=ρsMVstd1-ϕρcgVmPLPL+P2.Finally, the apparent permeability model for gas transport in shale nanopores can be expressed as(17)Kapp=1+ϕadϕKϕ=1+1+2ϕiniΔεϕiniϕad=ρsMVstd1-ϕρcgVmPLPL+P2Δε=VLρsRTEVstdln⁡PL+PiniPL+P+cPPini-PK=w1+αKnϕτreff2PVstd8RT+1-w23ϕτreffμVstdλ¯2rDf-28MπRTreff=r1+1+2ϕiniΔε.


## 3. Productivity Prediction Model for MFHWs

### 3.1. Physical Model

Figure1 shows the diagram of MFHWs in a shale gas reservoir. The reservoir is divided into seven contiguous regions: two upper-reservoir regions (regions 6 and 7), two outer-reservoir regions (regions 4 and 5), two inner-reservoir regions (regions 2 and 3), and a hydraulic region (region 1). Not all reservoir zones between primary hydraulic fractures are stimulated by induced fractures. 1D linear flow is assumed within unstimulated region and the flow direction depends on the location. As shown in Figure 1, a region of higher permeability around each fracture, the so-called SRV area, is introduced to represent fracture branching (region 2). The SRV region occupies part of the space between fractures and the flow in the other part (region 3) towards region 2, parallel to the wellbore. Regions 4 and 5 are in the outer reservoir beyond the tips of the hydraulic fractures and they connect with regions 2 and 3 in x direction. Because of the higher permeability in region 2, the pressure decreases in this area first. Therefore, the flow direction in regions 4 and 5 is towards regions 2 and 3, perpendicular to the wellbore. Although the main extending direction of the hydraulic fractures is horizontal (x direction in Figure 1), we also need to consider fractures height in vertical direction (z direction in Figure 1). Based on this consideration, regions 6 and 7 are defined in the upper reservoir beyond the tips of the hydraulic (connecting with regions 2 and 3 in z direction) and the flow direction is vertical to the wellbore, towards the lower reservoir.Figure 1
Diagram of MFHWs and of the seven-flow-region model.The main assumptions considered for this model are listed below:(1)
The shale reservoir has been fractured and has a constant reservoir thickness. All hydraulic fractures are symmetric and arranged uniformly along the horizontal well. Length and conductivity of the hydraulic fractures are the same.(2)
Natural fracture permeability is stress-dependent and shale gas flow in the matrix is driven by concentration difference.(3)
The initial pressure throughout the reservoir is uniform.(4)
Hydraulic fractures penetrate the entire pay formation. Fracture interference is ignored. After SRV, large fracture networks can be generated in region 2 and it has higher natural fracture permeability compared to the other regions.(5)
The shale gas flow in the reservoir is considered to be an isothermal flow and the effect of gravity and capillary pressure are negligible.(6)
Nonsimulated areas are described by a dual-porosity system consisting of a shale matrix and natural fractures.(7)
The only route for shale gas to reach the horizontal wellbore is through hydraulic fractures. The flow rate of MFHW is the sum of the flow rate from every hydraulic fracture.


### 3.2. Mathematical Model

According to the mass conservation equation, the motion equation, and the state equation, the partial differential formulas for shale gas flow in a matrix system and a fracture system can be obtained as follows [19, 27]:(18)1rmD2∂∂rmDrmD2∂ψmD∂rmD=1ηmD∂ψmD∂tD(19)∂2ψfD∂xD2-βD∂ψfD∂xD2=exp⁡βDψfD2Kmapp,avehref2Kfihfrm∂ψmD∂rmDrmD=1+1ηfD∂ψfD∂tD,where βD is a characteristic dimensionless parameter of stress sensitivity and relates to the permeability modulus γ. β is defined as β=(μZ/2P)γ, where Kf=Kfiexp⁡(-γΔPf).Pedrosa’s substitution [43] is used to weaken the nonlinearity of (19). By setting ψfD=(1/γD)ln⁡1-γDξD, (19) can be simplified into (20)∂2ξfD∂xD2=2Kmapp,avehref2Kfihfrm∂ψmD∂rmDrmD=1+1ηfD∂ξfD∂tD.To solve (18) and (20), these two equations were converted to the Laplace domain. And then the definite boundary continuity is introduced to develop an integral seepage flow differential equation. The final equation can be written as(21)∂2ψ¯fD∂xD2=fssψ¯fD,where fs=1/ηfD-(λ/5s)1-s/ηmDcoth⁡s/ηmD.A series of variables and dimensionless variables are defined in Table1.Table 1
Definitions of variables.
Variable
Expression
Dimensionless expression
Pseudo-pressure, psi2/cp
ψ
=
∫
P
i
n
i
P
2
P
μ
Z
d
P
ψ
D
=
K
r
e
f
h
r
e
f
2704
q
s
c
T
ψ
i
n
i
-
ψ
Permeability modulus factor, psi−1∗hr
β
=
μ
Z
38
P
γ
β
D
=
2704
q
s
c
T
K
r
e
f
h
r
e
f
β
Diffusivity, ft2/hr
η
r
e
f
=
2.637
×
1
0
-
4
K
r
e
f
ϕ
r
e
f
c
t
r
e
f
μ
r
e
f
η
m
D
=
η
m
η
r
e
f
η
m
=
2.637
×
1
0
-
4
K
m
,
a
v
e
ϕ
m
a
a
v
e
c
m
a
v
e
μ
m
a
v
e
η
f
D
=
η
f
η
r
e
f
η
f
=
2.637
×
1
0
-
4
K
f
ϕ
f
c
f
,
a
v
e
μ
f
,
a
v
e
η
F
D
=
η
F
η
r
e
f
η
F
=
2.637
×
1
0
-
4
K
F
ϕ
F
c
F
,
a
v
e
μ
F
,
a
v
e
Adsorption coefficient, dimensionless
a
a
v
e
=
c
m
,
a
v
e
+
(
1
-
ϕ
m
)
/
ϕ
m
T
P
s
c
Z
a
v
e
/
T
s
c
P
m
,
a
v
e
Z
s
c
V
L
P
L
/
(
P
L
+
P
m
,
a
v
e
)
2
c
m
,
a
v
e
Time, hr
t
D
=
η
r
e
f
h
r
e
f
2
t
Length, ft
x
1
D
=
x
f
h
r
e
f
x
e
D
=
x
e
h
r
e
f
y
1
D
=
y
1
h
r
e
f
y
2
D
=
y
2
h
r
e
f
z
1
D
=
h
F
2
h
r
e
f
z
2
D
=
h
h
r
e
f
w
D
=
w
F
h
r
e
f
r
m
D
=
r
r
m


#### 3.2.1. Region 7

In this upper-reservoir region, the governing diffusivity equation for shale gas flow in thez-direction becomes(22)∂2ψ¯fD7∂zD2-f7ssψ¯fD7=0,where(23)f7s=1ηfD7-λ75s1-sηmD7coth⁡sηmD7.At the top of the reservoir(z=z2), the boundary condition is no-flow:(24)∂ψ¯fD7∂zDzD=z2D=0.On the interfaces(z=z1) of regions 7 and 3 and regions 7 and 5, pressure continuity condition can be obtained as(25)ψ¯fD7zD=z1D=ψ¯fD3zD=z1D=ψ¯fD5zD=z1D.Equation (22) is solved simultaneously with (24) and (25). The final solution for (22) can be given as(26)∂ψ¯fD7∂zDzD=z1D=f7sstanh⁡f7ssz1D-z2Dψ¯fD3zD=z1D=f7sstanh⁡f7ssz1D-z2Dψ¯fD5zD=z1D.


#### 3.2.2. Region 6

Consistent with region 7, the governing diffusivity equation for shale gas flow in thez-direction becomes(27)∂2ψ¯fD6∂zD2-f6ssψ¯fD6=0,where (28)f6s=1ηfD6-λ65s1-sηmD6coth⁡sηmD6.The boundary condition at the top of the reservoir(z=z2) is(29)∂ψ¯fD6∂zDzD=z2D=0.The boundary conditions on the interfaces of regions 6 and 2 and regions 6 and 4 are(30)ψ¯fD6zD=z1D=ψ¯fD2zD=z1D=ψ¯fD5zD=z1D.The solution for (27) can be given as(31)∂ψ¯fD6∂zDzD=z1D=f6sstanh⁡f6ssz1D-z2Dψ¯fD2zD=z1D=f6sstanh⁡f6ssz1D-z2Dψ¯fD4zD=z1D.


#### 3.2.3. Region 5

In this region, shale gas flow takes place in thex-direction and the z-direction. So the governing diffusivity equation for region 5 is(32)∂2ψ¯fD5∂xD2+∂2ψ¯fD5∂zD2-f5ssψ¯fD5=0,where (33)f5s=1ηfD5-λ55s1-sηmD5coth⁡sηmD5.Along the horizontal well(z=0), there is a no-flow boundary:(34)∂ψ¯fD5∂zDzD=0=0.Flux across the interface of regions 5 and 7(z=z1) is continuous; therefore,(35)Kf5∂ψ¯fD5∂zDzD=z1D=Kf7∂ψ¯fD7∂zDzD=z1D.Substituting (34) and (35) into (32), each of its terms from 0 to z1D is integrated with respect to zD. The governing diffusivity can be simplified as(36)∂2ψ¯fD5∂xD2+Kf7Kf5z1D∂ψ¯fD7∂zDzD=z1D-f5ssψ¯fD5=0.In this region, shale gas flow in thex-direction and with a closed boundary condition (x=xe) is expressed as (37)∂ψ¯fD5∂xDxD=xeD=0.Pressure continuity condition between regions 5 and 3 (atx=x1) is expressed as(38)ψ¯fD5xD=x1D=ψ¯fD2xD=x1D.By solving (36), the solution can be obtained as(39)∂ψ¯fD5∂xDxD=x1D=α5tanh⁡α5x1D-xeDψ¯fD3xD=x1D,where (40)α5=f5ss-Kf7Kf5z1Df7sstanh⁡f7ssz1D-z2D.


#### 3.2.4. Region 4

In a similar way, shale gas flow in this region takes place in thex-direction and the z-direction. The governing diffusivity equation for region 4 is(41)∂2ψ¯fD4∂xD2+∂2ψ¯fD4∂zD2-f4ssψ¯fD4=0,where (42)f4s=1ηfD4-λ45s1-sηmD4coth⁡sηmD4.The same method is also used in region 5, where the governing diffusivity equation can be simplified as(43)∂2ψ¯fD4∂xD2+Kf4Kf6z1D∂ψ¯fD6∂zDzD=z1D-f4ssψ¯fD4=0.No-flow condition at the outer-reservoir boundary(x=xe) is given as(44)∂ψ¯fD4∂xDxD=xeD=0.On the interface of regions 4 and 2, the boundary condition is(45)ψ¯fD4xD=x1D=ψ¯fD2xD=x1D.The solution for (41) is given as(46)∂ψ¯fD4∂xDxD=x1D=α4tanh⁡α4x1D-xeDψ¯fD2xD=x1D,where (47)α4=f4ss-Kf6Kf4f6sstanh⁡f6ssz1D-z2D.


#### 3.2.5. Region 3

This region is adjacent to regions 5 and 7, so shale gas flow direction takes place along thex-direction, y-direction, and z-direction. The governing diffusivity equation becomes(48)∂2ψ¯fD3∂yD2+∂2ψ¯fD3∂xD2+∂2ψ¯fD3∂zD2-f3ssψ¯fD3=0,where (49)f3s=1ηfD3-λ35s1-sηmD3coth⁡sηmD3.For thex-direction, the boundary condition used is a no-flow condition for a horizontal well (z=0) and a flux continuity condition is used on the interface of regions 5 and 3 (x=x1):(50)∂ψ¯fD3∂xDxD=0=0Kf3∂ψ¯fD3∂xDxD=x1D=Kf5∂ψ¯fD5∂xDxD=x1D.For thez-direction, (51) can be obtained:(51)∂ψ¯fD3∂zDzD=0=0Kf3∂ψ¯fD3∂zDzD=z1D=Kf7∂ψ¯fD7∂zDzD=z1D.By substituting (50) and (51) into (48), the final governing diffusivity equation becomes(52)∂2ψ¯fD3∂yD2+Kf7Kf3z1D∂ψ¯fD7∂zDzD=z1D+Kf5Kf3x1D∂ψ¯fD5∂xDxD=x1D-f3ssψ¯fD3=0.Between two fractures, there is a no-flow boundary (aty=y2):(53)∂ψ¯fD3∂yDyD=y2D=0.On the interface of regions 3 and 2, the pressure continuity condition is written as(54)ψ¯fD3yD=y1D=ψ¯fD2yD=y1D.The solution for (48) is given as(55)∂ψ¯fD3∂yDyD=y1D=α3tanh⁡α3y1D-y2Dψ¯fD2yD=y1D,where (56)α3=f3ss-Kf7Kf3z1Df7sstanh⁡f7ssz1D-z2D-Kf5Kf3x1Dα5tanh⁡α5x1D-xeD.


#### 3.2.6. Region 2

Similarly, as region 2 is adjacent to hydraulic fracture, the governing diffusivity equation is(57)∂2ψ¯fD2∂yD2+∂2ψ¯fD2∂xD2+∂2ψ¯fD2∂zD2-f2ssψ¯fD2=0.Consistent with region 3, the simplified equation can be written as(58)∂2ψ¯fD2∂yD2+Kf6Kf2z1D∂ψ¯fD6∂zDzD=z1D+Kf4Kf2x1D∂ψ¯fD4∂xDxD=x1D-f2ssψ¯fD2=0,where (59)f2s=1ηfD2-λ25s1-sηmD2coth⁡sηmD2.At the interface of regions 2 and 3, the boundary condition is characterized by flux continuity:(60)Kf2∂ψ¯fD2∂yDyD=y1D=Kf3∂ψ¯fD3∂yDyD=y1D.And between the hydraulic fracture level and region 2, the boundary condition is(61)ψ¯fD2yD=wD/2=ψ¯FDyD=wD/2.Solving (58) to (61),(62)∂ψ¯fD2∂yDyD=wD/2=-β2ψ¯FDyD=wD/2,where (63)β2=α2exp⁡-α2wD/2-β1exp⁡α2wD/2exp⁡-α2wD/2+β1exp⁡α2wD/2β1=exp⁡-2α2y1DKf2α2+Kf3α3tanh⁡α3y1D-y2DKf2α2-Kf3α3tanh⁡α3y1D-y2Dα2=f2ss-Kf6Kf2z1Df6sstanh⁡f6ssz1D-z2D-Kf4Kf2x1Dα4tanh⁡α4x1D-xeD.


#### 3.2.7. Region 1

In this region, the diffusivity equation is given by(64)∂2ψ¯pFD∂xD2+2FCD∂ψ¯fD2∂yDyD=wD/2-sηFDψ¯pFD=0,where FCD=KFwD/Kf2. The boundary condition can be stated by the expressions:(65)∂ψ¯pFD∂xDxD=x1D=0∂ψ¯pFD∂xDxD=x1D=-πFCDs.The solution can be obtained:(66)ψ¯pFD=πcosh⁡αFxD-x1DαFFCDssinh⁡αFx1D,where (67)αF=2β2FCD+sηFD.Assuming that pressure distribution along the horizontal well is uniform, the dimensionless down-hole pressure can be solved by settingxD=0 [44]:(68)ψ¯wD=ψ¯pFDxD=0=παFFCDstanh⁡αFx1D.Shale gas flow in hydraulic fractures is not completely linear. The nonlinear streamlines result in an additional pressure drop. Mukherjee and Economides [45] presented a formulation to take this additional pressure drop into consideration.(69)Sc=Kf1hKFwfln⁡h2rw-π2.Based on the Duhamel principle and the superposition principle, the dimensionless flowing bottom hole pressure (FHBP), which takes into consideration the wellbore storage effect and the skin effect, can be expressed as [30](70)ψ¯wDSc,CD=Sc+sψ¯wDs1+CDssψ¯wD+Sc.The FHBP for MFHWs, when considering the stress sensitivity effect, can be obtained as(71)ψ¯wDSc,CD,tD=-1γDln⁡1-γDψ¯wDSc,CD.The relationship between the dimensionless rate under the condition of constant pressure and dimensionless pressure caused by the constant rate is described as [46](72)q¯D=1s2ψ¯wD.


## 3.1. Physical Model

Figure1 shows the diagram of MFHWs in a shale gas reservoir. The reservoir is divided into seven contiguous regions: two upper-reservoir regions (regions 6 and 7), two outer-reservoir regions (regions 4 and 5), two inner-reservoir regions (regions 2 and 3), and a hydraulic region (region 1). Not all reservoir zones between primary hydraulic fractures are stimulated by induced fractures. 1D linear flow is assumed within unstimulated region and the flow direction depends on the location. As shown in Figure 1, a region of higher permeability around each fracture, the so-called SRV area, is introduced to represent fracture branching (region 2). The SRV region occupies part of the space between fractures and the flow in the other part (region 3) towards region 2, parallel to the wellbore. Regions 4 and 5 are in the outer reservoir beyond the tips of the hydraulic fractures and they connect with regions 2 and 3 in x direction. Because of the higher permeability in region 2, the pressure decreases in this area first. Therefore, the flow direction in regions 4 and 5 is towards regions 2 and 3, perpendicular to the wellbore. Although the main extending direction of the hydraulic fractures is horizontal (x direction in Figure 1), we also need to consider fractures height in vertical direction (z direction in Figure 1). Based on this consideration, regions 6 and 7 are defined in the upper reservoir beyond the tips of the hydraulic (connecting with regions 2 and 3 in z direction) and the flow direction is vertical to the wellbore, towards the lower reservoir.Figure 1
Diagram of MFHWs and of the seven-flow-region model.The main assumptions considered for this model are listed below:(1)
The shale reservoir has been fractured and has a constant reservoir thickness. All hydraulic fractures are symmetric and arranged uniformly along the horizontal well. Length and conductivity of the hydraulic fractures are the same.(2)
Natural fracture permeability is stress-dependent and shale gas flow in the matrix is driven by concentration difference.(3)
The initial pressure throughout the reservoir is uniform.(4)
Hydraulic fractures penetrate the entire pay formation. Fracture interference is ignored. After SRV, large fracture networks can be generated in region 2 and it has higher natural fracture permeability compared to the other regions.(5)
The shale gas flow in the reservoir is considered to be an isothermal flow and the effect of gravity and capillary pressure are negligible.(6)
Nonsimulated areas are described by a dual-porosity system consisting of a shale matrix and natural fractures.(7)
The only route for shale gas to reach the horizontal wellbore is through hydraulic fractures. The flow rate of MFHW is the sum of the flow rate from every hydraulic fracture.


## 3.2. Mathematical Model

According to the mass conservation equation, the motion equation, and the state equation, the partial differential formulas for shale gas flow in a matrix system and a fracture system can be obtained as follows [19, 27]:(18)1rmD2∂∂rmDrmD2∂ψmD∂rmD=1ηmD∂ψmD∂tD(19)∂2ψfD∂xD2-βD∂ψfD∂xD2=exp⁡βDψfD2Kmapp,avehref2Kfihfrm∂ψmD∂rmDrmD=1+1ηfD∂ψfD∂tD,where βD is a characteristic dimensionless parameter of stress sensitivity and relates to the permeability modulus γ. β is defined as β=(μZ/2P)γ, where Kf=Kfiexp⁡(-γΔPf).Pedrosa’s substitution [43] is used to weaken the nonlinearity of (19). By setting ψfD=(1/γD)ln⁡1-γDξD, (19) can be simplified into (20)∂2ξfD∂xD2=2Kmapp,avehref2Kfihfrm∂ψmD∂rmDrmD=1+1ηfD∂ξfD∂tD.To solve (18) and (20), these two equations were converted to the Laplace domain. And then the definite boundary continuity is introduced to develop an integral seepage flow differential equation. The final equation can be written as(21)∂2ψ¯fD∂xD2=fssψ¯fD,where fs=1/ηfD-(λ/5s)1-s/ηmDcoth⁡s/ηmD.A series of variables and dimensionless variables are defined in Table1.Table 1
Definitions of variables.
Variable
Expression
Dimensionless expression
Pseudo-pressure, psi2/cp
ψ
=
∫
P
i
n
i
P
2
P
μ
Z
d
P
ψ
D
=
K
r
e
f
h
r
e
f
2704
q
s
c
T
ψ
i
n
i
-
ψ
Permeability modulus factor, psi−1∗hr
β
=
μ
Z
38
P
γ
β
D
=
2704
q
s
c
T
K
r
e
f
h
r
e
f
β
Diffusivity, ft2/hr
η
r
e
f
=
2.637
×
1
0
-
4
K
r
e
f
ϕ
r
e
f
c
t
r
e
f
μ
r
e
f
η
m
D
=
η
m
η
r
e
f
η
m
=
2.637
×
1
0
-
4
K
m
,
a
v
e
ϕ
m
a
a
v
e
c
m
a
v
e
μ
m
a
v
e
η
f
D
=
η
f
η
r
e
f
η
f
=
2.637
×
1
0
-
4
K
f
ϕ
f
c
f
,
a
v
e
μ
f
,
a
v
e
η
F
D
=
η
F
η
r
e
f
η
F
=
2.637
×
1
0
-
4
K
F
ϕ
F
c
F
,
a
v
e
μ
F
,
a
v
e
Adsorption coefficient, dimensionless
a
a
v
e
=
c
m
,
a
v
e
+
(
1
-
ϕ
m
)
/
ϕ
m
T
P
s
c
Z
a
v
e
/
T
s
c
P
m
,
a
v
e
Z
s
c
V
L
P
L
/
(
P
L
+
P
m
,
a
v
e
)
2
c
m
,
a
v
e
Time, hr
t
D
=
η
r
e
f
h
r
e
f
2
t
Length, ft
x
1
D
=
x
f
h
r
e
f
x
e
D
=
x
e
h
r
e
f
y
1
D
=
y
1
h
r
e
f
y
2
D
=
y
2
h
r
e
f
z
1
D
=
h
F
2
h
r
e
f
z
2
D
=
h
h
r
e
f
w
D
=
w
F
h
r
e
f
r
m
D
=
r
r
m


### 3.2.1. Region 7

In this upper-reservoir region, the governing diffusivity equation for shale gas flow in thez-direction becomes(22)∂2ψ¯fD7∂zD2-f7ssψ¯fD7=0,where(23)f7s=1ηfD7-λ75s1-sηmD7coth⁡sηmD7.At the top of the reservoir(z=z2), the boundary condition is no-flow:(24)∂ψ¯fD7∂zDzD=z2D=0.On the interfaces(z=z1) of regions 7 and 3 and regions 7 and 5, pressure continuity condition can be obtained as(25)ψ¯fD7zD=z1D=ψ¯fD3zD=z1D=ψ¯fD5zD=z1D.Equation (22) is solved simultaneously with (24) and (25). The final solution for (22) can be given as(26)∂ψ¯fD7∂zDzD=z1D=f7sstanh⁡f7ssz1D-z2Dψ¯fD3zD=z1D=f7sstanh⁡f7ssz1D-z2Dψ¯fD5zD=z1D.


### 3.2.2. Region 6

Consistent with region 7, the governing diffusivity equation for shale gas flow in thez-direction becomes(27)∂2ψ¯fD6∂zD2-f6ssψ¯fD6=0,where (28)f6s=1ηfD6-λ65s1-sηmD6coth⁡sηmD6.The boundary condition at the top of the reservoir(z=z2) is(29)∂ψ¯fD6∂zDzD=z2D=0.The boundary conditions on the interfaces of regions 6 and 2 and regions 6 and 4 are(30)ψ¯fD6zD=z1D=ψ¯fD2zD=z1D=ψ¯fD5zD=z1D.The solution for (27) can be given as(31)∂ψ¯fD6∂zDzD=z1D=f6sstanh⁡f6ssz1D-z2Dψ¯fD2zD=z1D=f6sstanh⁡f6ssz1D-z2Dψ¯fD4zD=z1D.


### 3.2.3. Region 5

In this region, shale gas flow takes place in thex-direction and the z-direction. So the governing diffusivity equation for region 5 is(32)∂2ψ¯fD5∂xD2+∂2ψ¯fD5∂zD2-f5ssψ¯fD5=0,where (33)f5s=1ηfD5-λ55s1-sηmD5coth⁡sηmD5.Along the horizontal well(z=0), there is a no-flow boundary:(34)∂ψ¯fD5∂zDzD=0=0.Flux across the interface of regions 5 and 7(z=z1) is continuous; therefore,(35)Kf5∂ψ¯fD5∂zDzD=z1D=Kf7∂ψ¯fD7∂zDzD=z1D.Substituting (34) and (35) into (32), each of its terms from 0 to z1D is integrated with respect to zD. The governing diffusivity can be simplified as(36)∂2ψ¯fD5∂xD2+Kf7Kf5z1D∂ψ¯fD7∂zDzD=z1D-f5ssψ¯fD5=0.In this region, shale gas flow in thex-direction and with a closed boundary condition (x=xe) is expressed as (37)∂ψ¯fD5∂xDxD=xeD=0.Pressure continuity condition between regions 5 and 3 (atx=x1) is expressed as(38)ψ¯fD5xD=x1D=ψ¯fD2xD=x1D.By solving (36), the solution can be obtained as(39)∂ψ¯fD5∂xDxD=x1D=α5tanh⁡α5x1D-xeDψ¯fD3xD=x1D,where (40)α5=f5ss-Kf7Kf5z1Df7sstanh⁡f7ssz1D-z2D.


### 3.2.4. Region 4

In a similar way, shale gas flow in this region takes place in thex-direction and the z-direction. The governing diffusivity equation for region 4 is(41)∂2ψ¯fD4∂xD2+∂2ψ¯fD4∂zD2-f4ssψ¯fD4=0,where (42)f4s=1ηfD4-λ45s1-sηmD4coth⁡sηmD4.The same method is also used in region 5, where the governing diffusivity equation can be simplified as(43)∂2ψ¯fD4∂xD2+Kf4Kf6z1D∂ψ¯fD6∂zDzD=z1D-f4ssψ¯fD4=0.No-flow condition at the outer-reservoir boundary(x=xe) is given as(44)∂ψ¯fD4∂xDxD=xeD=0.On the interface of regions 4 and 2, the boundary condition is(45)ψ¯fD4xD=x1D=ψ¯fD2xD=x1D.The solution for (41) is given as(46)∂ψ¯fD4∂xDxD=x1D=α4tanh⁡α4x1D-xeDψ¯fD2xD=x1D,where (47)α4=f4ss-Kf6Kf4f6sstanh⁡f6ssz1D-z2D.


### 3.2.5. Region 3

This region is adjacent to regions 5 and 7, so shale gas flow direction takes place along thex-direction, y-direction, and z-direction. The governing diffusivity equation becomes(48)∂2ψ¯fD3∂yD2+∂2ψ¯fD3∂xD2+∂2ψ¯fD3∂zD2-f3ssψ¯fD3=0,where (49)f3s=1ηfD3-λ35s1-sηmD3coth⁡sηmD3.For thex-direction, the boundary condition used is a no-flow condition for a horizontal well (z=0) and a flux continuity condition is used on the interface of regions 5 and 3 (x=x1):(50)∂ψ¯fD3∂xDxD=0=0Kf3∂ψ¯fD3∂xDxD=x1D=Kf5∂ψ¯fD5∂xDxD=x1D.For thez-direction, (51) can be obtained:(51)∂ψ¯fD3∂zDzD=0=0Kf3∂ψ¯fD3∂zDzD=z1D=Kf7∂ψ¯fD7∂zDzD=z1D.By substituting (50) and (51) into (48), the final governing diffusivity equation becomes(52)∂2ψ¯fD3∂yD2+Kf7Kf3z1D∂ψ¯fD7∂zDzD=z1D+Kf5Kf3x1D∂ψ¯fD5∂xDxD=x1D-f3ssψ¯fD3=0.Between two fractures, there is a no-flow boundary (aty=y2):(53)∂ψ¯fD3∂yDyD=y2D=0.On the interface of regions 3 and 2, the pressure continuity condition is written as(54)ψ¯fD3yD=y1D=ψ¯fD2yD=y1D.The solution for (48) is given as(55)∂ψ¯fD3∂yDyD=y1D=α3tanh⁡α3y1D-y2Dψ¯fD2yD=y1D,where (56)α3=f3ss-Kf7Kf3z1Df7sstanh⁡f7ssz1D-z2D-Kf5Kf3x1Dα5tanh⁡α5x1D-xeD.


### 3.2.6. Region 2

Similarly, as region 2 is adjacent to hydraulic fracture, the governing diffusivity equation is(57)∂2ψ¯fD2∂yD2+∂2ψ¯fD2∂xD2+∂2ψ¯fD2∂zD2-f2ssψ¯fD2=0.Consistent with region 3, the simplified equation can be written as(58)∂2ψ¯fD2∂yD2+Kf6Kf2z1D∂ψ¯fD6∂zDzD=z1D+Kf4Kf2x1D∂ψ¯fD4∂xDxD=x1D-f2ssψ¯fD2=0,where (59)f2s=1ηfD2-λ25s1-sηmD2coth⁡sηmD2.At the interface of regions 2 and 3, the boundary condition is characterized by flux continuity:(60)Kf2∂ψ¯fD2∂yDyD=y1D=Kf3∂ψ¯fD3∂yDyD=y1D.And between the hydraulic fracture level and region 2, the boundary condition is(61)ψ¯fD2yD=wD/2=ψ¯FDyD=wD/2.Solving (58) to (61),(62)∂ψ¯fD2∂yDyD=wD/2=-β2ψ¯FDyD=wD/2,where (63)β2=α2exp⁡-α2wD/2-β1exp⁡α2wD/2exp⁡-α2wD/2+β1exp⁡α2wD/2β1=exp⁡-2α2y1DKf2α2+Kf3α3tanh⁡α3y1D-y2DKf2α2-Kf3α3tanh⁡α3y1D-y2Dα2=f2ss-Kf6Kf2z1Df6sstanh⁡f6ssz1D-z2D-Kf4Kf2x1Dα4tanh⁡α4x1D-xeD.


### 3.2.7. Region 1

In this region, the diffusivity equation is given by(64)∂2ψ¯pFD∂xD2+2FCD∂ψ¯fD2∂yDyD=wD/2-sηFDψ¯pFD=0,where FCD=KFwD/Kf2. The boundary condition can be stated by the expressions:(65)∂ψ¯pFD∂xDxD=x1D=0∂ψ¯pFD∂xDxD=x1D=-πFCDs.The solution can be obtained:(66)ψ¯pFD=πcosh⁡αFxD-x1DαFFCDssinh⁡αFx1D,where (67)αF=2β2FCD+sηFD.Assuming that pressure distribution along the horizontal well is uniform, the dimensionless down-hole pressure can be solved by settingxD=0 [44]:(68)ψ¯wD=ψ¯pFDxD=0=παFFCDstanh⁡αFx1D.Shale gas flow in hydraulic fractures is not completely linear. The nonlinear streamlines result in an additional pressure drop. Mukherjee and Economides [45] presented a formulation to take this additional pressure drop into consideration.(69)Sc=Kf1hKFwfln⁡h2rw-π2.Based on the Duhamel principle and the superposition principle, the dimensionless flowing bottom hole pressure (FHBP), which takes into consideration the wellbore storage effect and the skin effect, can be expressed as [30](70)ψ¯wDSc,CD=Sc+sψ¯wDs1+CDssψ¯wD+Sc.The FHBP for MFHWs, when considering the stress sensitivity effect, can be obtained as(71)ψ¯wDSc,CD,tD=-1γDln⁡1-γDψ¯wDSc,CD.The relationship between the dimensionless rate under the condition of constant pressure and dimensionless pressure caused by the constant rate is described as [46](72)q¯D=1s2ψ¯wD.


## 3.2.1. Region 7

In this upper-reservoir region, the governing diffusivity equation for shale gas flow in thez-direction becomes(22)∂2ψ¯fD7∂zD2-f7ssψ¯fD7=0,where(23)f7s=1ηfD7-λ75s1-sηmD7coth⁡sηmD7.At the top of the reservoir(z=z2), the boundary condition is no-flow:(24)∂ψ¯fD7∂zDzD=z2D=0.On the interfaces(z=z1) of regions 7 and 3 and regions 7 and 5, pressure continuity condition can be obtained as(25)ψ¯fD7zD=z1D=ψ¯fD3zD=z1D=ψ¯fD5zD=z1D.Equation (22) is solved simultaneously with (24) and (25). The final solution for (22) can be given as(26)∂ψ¯fD7∂zDzD=z1D=f7sstanh⁡f7ssz1D-z2Dψ¯fD3zD=z1D=f7sstanh⁡f7ssz1D-z2Dψ¯fD5zD=z1D.


## 3.2.2. Region 6

Consistent with region 7, the governing diffusivity equation for shale gas flow in thez-direction becomes(27)∂2ψ¯fD6∂zD2-f6ssψ¯fD6=0,where (28)f6s=1ηfD6-λ65s1-sηmD6coth⁡sηmD6.The boundary condition at the top of the reservoir(z=z2) is(29)∂ψ¯fD6∂zDzD=z2D=0.The boundary conditions on the interfaces of regions 6 and 2 and regions 6 and 4 are(30)ψ¯fD6zD=z1D=ψ¯fD2zD=z1D=ψ¯fD5zD=z1D.The solution for (27) can be given as(31)∂ψ¯fD6∂zDzD=z1D=f6sstanh⁡f6ssz1D-z2Dψ¯fD2zD=z1D=f6sstanh⁡f6ssz1D-z2Dψ¯fD4zD=z1D.


## 3.2.3. Region 5

In this region, shale gas flow takes place in thex-direction and the z-direction. So the governing diffusivity equation for region 5 is(32)∂2ψ¯fD5∂xD2+∂2ψ¯fD5∂zD2-f5ssψ¯fD5=0,where (33)f5s=1ηfD5-λ55s1-sηmD5coth⁡sηmD5.Along the horizontal well(z=0), there is a no-flow boundary:(34)∂ψ¯fD5∂zDzD=0=0.Flux across the interface of regions 5 and 7(z=z1) is continuous; therefore,(35)Kf5∂ψ¯fD5∂zDzD=z1D=Kf7∂ψ¯fD7∂zDzD=z1D.Substituting (34) and (35) into (32), each of its terms from 0 to z1D is integrated with respect to zD. The governing diffusivity can be simplified as(36)∂2ψ¯fD5∂xD2+Kf7Kf5z1D∂ψ¯fD7∂zDzD=z1D-f5ssψ¯fD5=0.In this region, shale gas flow in thex-direction and with a closed boundary condition (x=xe) is expressed as (37)∂ψ¯fD5∂xDxD=xeD=0.Pressure continuity condition between regions 5 and 3 (atx=x1) is expressed as(38)ψ¯fD5xD=x1D=ψ¯fD2xD=x1D.By solving (36), the solution can be obtained as(39)∂ψ¯fD5∂xDxD=x1D=α5tanh⁡α5x1D-xeDψ¯fD3xD=x1D,where (40)α5=f5ss-Kf7Kf5z1Df7sstanh⁡f7ssz1D-z2D.


## 3.2.4. Region 4

In a similar way, shale gas flow in this region takes place in thex-direction and the z-direction. The governing diffusivity equation for region 4 is(41)∂2ψ¯fD4∂xD2+∂2ψ¯fD4∂zD2-f4ssψ¯fD4=0,where (42)f4s=1ηfD4-λ45s1-sηmD4coth⁡sηmD4.The same method is also used in region 5, where the governing diffusivity equation can be simplified as(43)∂2ψ¯fD4∂xD2+Kf4Kf6z1D∂ψ¯fD6∂zDzD=z1D-f4ssψ¯fD4=0.No-flow condition at the outer-reservoir boundary(x=xe) is given as(44)∂ψ¯fD4∂xDxD=xeD=0.On the interface of regions 4 and 2, the boundary condition is(45)ψ¯fD4xD=x1D=ψ¯fD2xD=x1D.The solution for (41) is given as(46)∂ψ¯fD4∂xDxD=x1D=α4tanh⁡α4x1D-xeDψ¯fD2xD=x1D,where (47)α4=f4ss-Kf6Kf4f6sstanh⁡f6ssz1D-z2D.


## 3.2.5. Region 3

This region is adjacent to regions 5 and 7, so shale gas flow direction takes place along thex-direction, y-direction, and z-direction. The governing diffusivity equation becomes(48)∂2ψ¯fD3∂yD2+∂2ψ¯fD3∂xD2+∂2ψ¯fD3∂zD2-f3ssψ¯fD3=0,where (49)f3s=1ηfD3-λ35s1-sηmD3coth⁡sηmD3.For thex-direction, the boundary condition used is a no-flow condition for a horizontal well (z=0) and a flux continuity condition is used on the interface of regions 5 and 3 (x=x1):(50)∂ψ¯fD3∂xDxD=0=0Kf3∂ψ¯fD3∂xDxD=x1D=Kf5∂ψ¯fD5∂xDxD=x1D.For thez-direction, (51) can be obtained:(51)∂ψ¯fD3∂zDzD=0=0Kf3∂ψ¯fD3∂zDzD=z1D=Kf7∂ψ¯fD7∂zDzD=z1D.By substituting (50) and (51) into (48), the final governing diffusivity equation becomes(52)∂2ψ¯fD3∂yD2+Kf7Kf3z1D∂ψ¯fD7∂zDzD=z1D+Kf5Kf3x1D∂ψ¯fD5∂xDxD=x1D-f3ssψ¯fD3=0.Between two fractures, there is a no-flow boundary (aty=y2):(53)∂ψ¯fD3∂yDyD=y2D=0.On the interface of regions 3 and 2, the pressure continuity condition is written as(54)ψ¯fD3yD=y1D=ψ¯fD2yD=y1D.The solution for (48) is given as(55)∂ψ¯fD3∂yDyD=y1D=α3tanh⁡α3y1D-y2Dψ¯fD2yD=y1D,where (56)α3=f3ss-Kf7Kf3z1Df7sstanh⁡f7ssz1D-z2D-Kf5Kf3x1Dα5tanh⁡α5x1D-xeD.


## 3.2.6. Region 2

Similarly, as region 2 is adjacent to hydraulic fracture, the governing diffusivity equation is(57)∂2ψ¯fD2∂yD2+∂2ψ¯fD2∂xD2+∂2ψ¯fD2∂zD2-f2ssψ¯fD2=0.Consistent with region 3, the simplified equation can be written as(58)∂2ψ¯fD2∂yD2+Kf6Kf2z1D∂ψ¯fD6∂zDzD=z1D+Kf4Kf2x1D∂ψ¯fD4∂xDxD=x1D-f2ssψ¯fD2=0,where (59)f2s=1ηfD2-λ25s1-sηmD2coth⁡sηmD2.At the interface of regions 2 and 3, the boundary condition is characterized by flux continuity:(60)Kf2∂ψ¯fD2∂yDyD=y1D=Kf3∂ψ¯fD3∂yDyD=y1D.And between the hydraulic fracture level and region 2, the boundary condition is(61)ψ¯fD2yD=wD/2=ψ¯FDyD=wD/2.Solving (58) to (61),(62)∂ψ¯fD2∂yDyD=wD/2=-β2ψ¯FDyD=wD/2,where (63)β2=α2exp⁡-α2wD/2-β1exp⁡α2wD/2exp⁡-α2wD/2+β1exp⁡α2wD/2β1=exp⁡-2α2y1DKf2α2+Kf3α3tanh⁡α3y1D-y2DKf2α2-Kf3α3tanh⁡α3y1D-y2Dα2=f2ss-Kf6Kf2z1Df6sstanh⁡f6ssz1D-z2D-Kf4Kf2x1Dα4tanh⁡α4x1D-xeD.


## 3.2.7. Region 1

In this region, the diffusivity equation is given by(64)∂2ψ¯pFD∂xD2+2FCD∂ψ¯fD2∂yDyD=wD/2-sηFDψ¯pFD=0,where FCD=KFwD/Kf2. The boundary condition can be stated by the expressions:(65)∂ψ¯pFD∂xDxD=x1D=0∂ψ¯pFD∂xDxD=x1D=-πFCDs.The solution can be obtained:(66)ψ¯pFD=πcosh⁡αFxD-x1DαFFCDssinh⁡αFx1D,where (67)αF=2β2FCD+sηFD.Assuming that pressure distribution along the horizontal well is uniform, the dimensionless down-hole pressure can be solved by settingxD=0 [44]:(68)ψ¯wD=ψ¯pFDxD=0=παFFCDstanh⁡αFx1D.Shale gas flow in hydraulic fractures is not completely linear. The nonlinear streamlines result in an additional pressure drop. Mukherjee and Economides [45] presented a formulation to take this additional pressure drop into consideration.(69)Sc=Kf1hKFwfln⁡h2rw-π2.Based on the Duhamel principle and the superposition principle, the dimensionless flowing bottom hole pressure (FHBP), which takes into consideration the wellbore storage effect and the skin effect, can be expressed as [30](70)ψ¯wDSc,CD=Sc+sψ¯wDs1+CDssψ¯wD+Sc.The FHBP for MFHWs, when considering the stress sensitivity effect, can be obtained as(71)ψ¯wDSc,CD,tD=-1γDln⁡1-γDψ¯wDSc,CD.The relationship between the dimensionless rate under the condition of constant pressure and dimensionless pressure caused by the constant rate is described as [46](72)q¯D=1s2ψ¯wD.


## 4. Discussion and Analysis

### 4.1. Apparent Permeability in the Shale Matrix

Figure2 shows the variation of the apparent permeabilities against the average pressure. The input parameters are listed in Table 2. Based on Figure 2, one can observe that, in the early period, with a decrease of the average pressure, the bulk permeability Kapp,bulk decreases dramatically, but the apparent permeability caused by gas desorption Kapp,ad increases greatly. In the later period, Kapp,bulk increases slightly with an increase of the pressure. Kapp,ad also increases in a similar way with an increase in pressure. The explanation is that, at high pressure, the Knudsen diffusion increases significantly. As can be seen in Figure 2, viscous flow dominates at high pressure and the Knudsen diffusion dominates at low pressure. Figure 2 also indicates that the pore diameter has a great impact on the apparent permeability. With a decrease of the pressure, nanopores become larger and larger under the influence of matrix shrinkage, so the apparent permeability increases gradually. However, there is an inflection point between the apparent permeability and the average pressure. Furthermore, pore diameter is bigger and the time of the inflection point for apparent permeability occurring appears later. The apparent permeability decreases faster after the inflection point. The results presented in Figure 2 agree well with those presented by Wu et al. [38].Table 2
Basic parameters.
Variable
Value
ϕ, dimensionless
0.05
r, m
10−8
P
i
n
i, MPa
25
T, K
333
k
B, J/K
1.38∗ 10−23
ρ
s, t/m3
2.56
c
p, MPa−1
4.35∗ 10−4
V
L, m3/t
18.66
E, MPa
4000
τ, dimensionless
4.3
μ, Pa∗s
1.75∗ 10−5
P
f
i
n, MPa
1
M, kg/mol
0.016
V
s
t
d, m3/mol
22.4∗ 10−3
ρ
g, t/m3
0.085
c
g, MPa−1
6.89∗ 10−2
P
L, MPa
2.72
α, dimensionless
2.0Figure 2
Variation ofKapp,ad, Kv, Kapp,Kn, and Kapp,t with average pressure.
(a)
(b)
(c)


### 4.2. Effects of Hydraulic Fracture Half-Length

Figures3 and 4 show the effect of hydraulic fracture half-length on pseudo-pressure, pseudo-pressure derivative, and dimensionless production. The input parameters are listed in Table 3. According to the results, fracture half-length affects the whole production process. The longer the hydraulic fracture half-length, the more significant the pressure drop and the higher the production rate. The contact area between the hydraulic fracture and the shale gas reservoir is larger when the fracture half-length is longer. Therefore, when gas flow resistance is reduced, more shale gas can flow to the fracture. Furthermore, larger flow rate results in significant pressure drop and a larger production.Table 3
Model parameters.
Variable
Value
Reservoir size inx-direction, xe, ft
500
Distance to no-flow boundary,y2, ft
250
Reservoir thickness,h, ft
250
Fracture width,wf, ft
0.01
Reservoir temperature,T, K
333
Gas volumetric factor,Bg
0.004
Initial matrix compressibilitycm, psi−1
10−3
Fracture porosity,ϕf
0.45
Fracture permeability,Kf, mD
103
Hydraulic fracture compressibility,cF, psi−1
10−3
Hydraulic fracture permeability,KF, mD
105
Well length,Lw, ft
1500
Fracture half-length,xf, ft
400
SRV width,y1, ft
250
Fracture height,hf, ft
0.001
Matrix radius,rm, ft
5
Initial viscosity,μini, cp
0.018
Initial matrix porosity,ϕm
0.05
Matrix permeability,Km, mD
10−6
Fracture compressibilitycf, psi−1
10−3
Hydraulic fracture porosity,ϕF
0.38Figure 3
Effect of hydraulic fracture half-length on pressure.Figure 4
Effect of hydraulic fracture half-length on production.


### 4.3. Effect of Hydraulic Fracture Permeability

Figures5 and 6 present the curve of the relationships of dimensionless time and pseudo-pressure and dimensionless production rate under different hydraulic fracture permeability. The figures show that hydraulic fracture permeability dominates all of the production period except for when we have early pure wellbore storage regime and boundary flow regime. The results are quite different from those observed in a single porosity tight reservoir where hydraulic fracture permeability only exerts an influence on early time pressure behavior [47]. The figures also indicate that, with an increase in fracture permeability, the pseudo-pressure drop and the production rate growth tend to slow down. Therefore, there is an optimal fracture permeability for a certain reservoir.Figure 5
Effect of hydraulic fracture permeability on pressure.Figure 6
Effect of hydraulic fracture permeability on production.


### 4.4. Effect of the Langmuir Volume.

The effects of the Langmuir volume on pressure and production are plotted in Figures7 and 8. The Langmuir volume is defined as the maximum volume of shale gas that can be absorbed to a unit weight organic matrix at infinite pressure. The Langmuir volume is a significant parameter for gas adsorption and desorption. As illustrated in Figures 7 and 8, the larger the Langmuir volume is, the lower the pressure drop and the production will be. This is because a larger Langmuir volume results in a greater amount of absorbed gas which provides an extra contribution for well reserves. However, the Langmuir volume does not affect the linear flow stage as only shale gas in hydraulic fractures contributes to production at that period. It can also be seen that, with an increase of the Langmuir volume, the dimensionless rate decreases slowly. This results in longer production time.Figure 7
Effect of Langmuir volume (VL) on pressure.Figure 8
Effect of Langmuir volume (VL) on production.


### 4.5. Effect of Matrix Radius

Figures9 and 10 show the effect of matrix block radius on pseudo-pressure, pseudo-pressure derivative, and dimensionless production rate. The shape and size of the matrix blocks exert a very significant effect on the interporosity flow stage. In this model, the value of the matrix block radius determines the magnitudes of the interporosity flow coefficient and the storativity ratio. According to the results, the matrix block radius mainly affects the transitional stage and the stage of bilinear flow. That is to say, the larger the matrix block radius, the more significant the pressure drop of the transitional stage. Consequently, the stage of bilinear flow occurs later. When the matrix block radius is large enough, bilinear flow disappears. Meanwhile, the magnitude of the matrix block radius also affects the production in the early and middle stages. With an increasing matrix block radius, the production rate gradually decreases, and the production decline rate decreases as well. In the late stage of production, the production rates under different matrix block radii are approximately the same [30].Figure 9
Effect of matrix radius on pressure.Figure 10
Effect of matrix radius on production.


### 4.6. Effect of Matrix Permeability

Figures11 and 12 demonstrate how matrix permeability affects pseudo-pressure and dimensionless production rate. As can be seen in these two figures, as the matrix permeability decreases, pressure drop in the transient period becomes more significant and the fracture linear flow period becomes longer. Dimensionless production rate increases slowly with an increase of matrix permeability. The results indicate that the rock has more capacity to transfer gas molecules from the matrix to the fracture system. Under ultra-low matrix permeability conditions, the pressure variation is more sensitive.Figure 11
Effect of matrix permeability on pressure.Figure 12
Effect of matrix permeability on production.


### 4.7. Effect of Induced Fracture Permeability

Figures13 and 14 represent the effect of induced fracture permeability on pressure and production. The induced fracture permeability of shale reservoirs is an extremely important reservoir parameter, as shale gas transport to hydraulic fracture should take place through induced fracture systems. Therefore, the induced fracture permeability determines the scale of gas well productivity. The results show that the higher the induced fracture permeability, the more significant the pressure drop during the transitional period. With a continuous increase in induced fracture permeability, the duration of the stage of bilinear flow becomes shorter and shorter until it completely disappears. Figure 13 indicates that the induced fracture permeability has a more significant influence on pressure than matrix radius and matrix permeability. Meanwhile, the induced fracture permeability also affects the scale of production rate; that is, the higher the permeability, the larger the scale. However, as far as production rate is concerned, the secondary fracture permeability only affects its scale, but not its decline rate.Figure 13
Effect of induced fracture permeability on pressure.Figure 14
Effect of induced fracture permeability on production.


## 4.1. Apparent Permeability in the Shale Matrix

Figure2 shows the variation of the apparent permeabilities against the average pressure. The input parameters are listed in Table 2. Based on Figure 2, one can observe that, in the early period, with a decrease of the average pressure, the bulk permeability Kapp,bulk decreases dramatically, but the apparent permeability caused by gas desorption Kapp,ad increases greatly. In the later period, Kapp,bulk increases slightly with an increase of the pressure. Kapp,ad also increases in a similar way with an increase in pressure. The explanation is that, at high pressure, the Knudsen diffusion increases significantly. As can be seen in Figure 2, viscous flow dominates at high pressure and the Knudsen diffusion dominates at low pressure. Figure 2 also indicates that the pore diameter has a great impact on the apparent permeability. With a decrease of the pressure, nanopores become larger and larger under the influence of matrix shrinkage, so the apparent permeability increases gradually. However, there is an inflection point between the apparent permeability and the average pressure. Furthermore, pore diameter is bigger and the time of the inflection point for apparent permeability occurring appears later. The apparent permeability decreases faster after the inflection point. The results presented in Figure 2 agree well with those presented by Wu et al. [38].Table 2
Basic parameters.
Variable
Value
ϕ, dimensionless
0.05
r, m
10−8
P
i
n
i, MPa
25
T, K
333
k
B, J/K
1.38∗ 10−23
ρ
s, t/m3
2.56
c
p, MPa−1
4.35∗ 10−4
V
L, m3/t
18.66
E, MPa
4000
τ, dimensionless
4.3
μ, Pa∗s
1.75∗ 10−5
P
f
i
n, MPa
1
M, kg/mol
0.016
V
s
t
d, m3/mol
22.4∗ 10−3
ρ
g, t/m3
0.085
c
g, MPa−1
6.89∗ 10−2
P
L, MPa
2.72
α, dimensionless
2.0Figure 2
Variation ofKapp,ad, Kv, Kapp,Kn, and Kapp,t with average pressure.
(a)
(b)
(c)


## 4.2. Effects of Hydraulic Fracture Half-Length

Figures3 and 4 show the effect of hydraulic fracture half-length on pseudo-pressure, pseudo-pressure derivative, and dimensionless production. The input parameters are listed in Table 3. According to the results, fracture half-length affects the whole production process. The longer the hydraulic fracture half-length, the more significant the pressure drop and the higher the production rate. The contact area between the hydraulic fracture and the shale gas reservoir is larger when the fracture half-length is longer. Therefore, when gas flow resistance is reduced, more shale gas can flow to the fracture. Furthermore, larger flow rate results in significant pressure drop and a larger production.Table 3
Model parameters.
Variable
Value
Reservoir size inx-direction, xe, ft
500
Distance to no-flow boundary,y2, ft
250
Reservoir thickness,h, ft
250
Fracture width,wf, ft
0.01
Reservoir temperature,T, K
333
Gas volumetric factor,Bg
0.004
Initial matrix compressibilitycm, psi−1
10−3
Fracture porosity,ϕf
0.45
Fracture permeability,Kf, mD
103
Hydraulic fracture compressibility,cF, psi−1
10−3
Hydraulic fracture permeability,KF, mD
105
Well length,Lw, ft
1500
Fracture half-length,xf, ft
400
SRV width,y1, ft
250
Fracture height,hf, ft
0.001
Matrix radius,rm, ft
5
Initial viscosity,μini, cp
0.018
Initial matrix porosity,ϕm
0.05
Matrix permeability,Km, mD
10−6
Fracture compressibilitycf, psi−1
10−3
Hydraulic fracture porosity,ϕF
0.38Figure 3
Effect of hydraulic fracture half-length on pressure.Figure 4
Effect of hydraulic fracture half-length on production.


## 4.3. Effect of Hydraulic Fracture Permeability

Figures5 and 6 present the curve of the relationships of dimensionless time and pseudo-pressure and dimensionless production rate under different hydraulic fracture permeability. The figures show that hydraulic fracture permeability dominates all of the production period except for when we have early pure wellbore storage regime and boundary flow regime. The results are quite different from those observed in a single porosity tight reservoir where hydraulic fracture permeability only exerts an influence on early time pressure behavior [47]. The figures also indicate that, with an increase in fracture permeability, the pseudo-pressure drop and the production rate growth tend to slow down. Therefore, there is an optimal fracture permeability for a certain reservoir.Figure 5
Effect of hydraulic fracture permeability on pressure.Figure 6
Effect of hydraulic fracture permeability on production.


## 4.4. Effect of the Langmuir Volume.

The effects of the Langmuir volume on pressure and production are plotted in Figures7 and 8. The Langmuir volume is defined as the maximum volume of shale gas that can be absorbed to a unit weight organic matrix at infinite pressure. The Langmuir volume is a significant parameter for gas adsorption and desorption. As illustrated in Figures 7 and 8, the larger the Langmuir volume is, the lower the pressure drop and the production will be. This is because a larger Langmuir volume results in a greater amount of absorbed gas which provides an extra contribution for well reserves. However, the Langmuir volume does not affect the linear flow stage as only shale gas in hydraulic fractures contributes to production at that period. It can also be seen that, with an increase of the Langmuir volume, the dimensionless rate decreases slowly. This results in longer production time.Figure 7
Effect of Langmuir volume (VL) on pressure.Figure 8
Effect of Langmuir volume (VL) on production.


## 4.5. Effect of Matrix Radius

Figures9 and 10 show the effect of matrix block radius on pseudo-pressure, pseudo-pressure derivative, and dimensionless production rate. The shape and size of the matrix blocks exert a very significant effect on the interporosity flow stage. In this model, the value of the matrix block radius determines the magnitudes of the interporosity flow coefficient and the storativity ratio. According to the results, the matrix block radius mainly affects the transitional stage and the stage of bilinear flow. That is to say, the larger the matrix block radius, the more significant the pressure drop of the transitional stage. Consequently, the stage of bilinear flow occurs later. When the matrix block radius is large enough, bilinear flow disappears. Meanwhile, the magnitude of the matrix block radius also affects the production in the early and middle stages. With an increasing matrix block radius, the production rate gradually decreases, and the production decline rate decreases as well. In the late stage of production, the production rates under different matrix block radii are approximately the same [30].Figure 9
Effect of matrix radius on pressure.Figure 10
Effect of matrix radius on production.


## 4.6. Effect of Matrix Permeability

Figures11 and 12 demonstrate how matrix permeability affects pseudo-pressure and dimensionless production rate. As can be seen in these two figures, as the matrix permeability decreases, pressure drop in the transient period becomes more significant and the fracture linear flow period becomes longer. Dimensionless production rate increases slowly with an increase of matrix permeability. The results indicate that the rock has more capacity to transfer gas molecules from the matrix to the fracture system. Under ultra-low matrix permeability conditions, the pressure variation is more sensitive.Figure 11
Effect of matrix permeability on pressure.Figure 12
Effect of matrix permeability on production.


## 4.7. Effect of Induced Fracture Permeability

Figures13 and 14 represent the effect of induced fracture permeability on pressure and production. The induced fracture permeability of shale reservoirs is an extremely important reservoir parameter, as shale gas transport to hydraulic fracture should take place through induced fracture systems. Therefore, the induced fracture permeability determines the scale of gas well productivity. The results show that the higher the induced fracture permeability, the more significant the pressure drop during the transitional period. With a continuous increase in induced fracture permeability, the duration of the stage of bilinear flow becomes shorter and shorter until it completely disappears. Figure 13 indicates that the induced fracture permeability has a more significant influence on pressure than matrix radius and matrix permeability. Meanwhile, the induced fracture permeability also affects the scale of production rate; that is, the higher the permeability, the larger the scale. However, as far as production rate is concerned, the secondary fracture permeability only affects its scale, but not its decline rate.Figure 13
Effect of induced fracture permeability on pressure.Figure 14
Effect of induced fracture permeability on production.


## 5. Case Study

The MFHW is located in the transitional zone between the Sichuan Basin and the Yunnan-Guizhou Plateau (400–1,200 m in altitude) and mainly developed a marine and continental sedimentary facies. The specific stratum parameters are listed in Table4.Table 4
Parameters of the typical well.
Porosity
Viscosity
Compressibility
Temperature
0.032
1.47∗ 10−2 cp
4.8∗ 10−2 MPa−1
96.8°C
Wellhead pressure
Depth
Pressure coefficient
Horizontal length
39 MPa
2900 m
1.96
1438 m
ψ
i
n
i
ψ
w
D
Fracturing segments
Thickness
1.62∗ 105 MPa2/cp
1.62∗ 105 MPa4/cp
15
35 mAs seen in Figure15, the production data show that there is a relatively high degree of correlation between pressure and production; the production data present an obvious response to the pressure fluctuations caused by the test carried out in the early stage; therefore, the data have a certain degree of reliability. As shown by a comparative analysis on the actual production data and production curve, the measured data are deviated from the typical curve by a relatively significant extent. This is mainly because of the skin effect and wellbore storage effect on the early stage seepage. The fitting curve shows that the data points are concentrated in the stage of linear flow of matrix (see Figure 16), mainly because the extremely low permeability of matrix-fractures has extended the duration of the stage of nonsteady state linear flow.Figure 15
Production history.Figure 16
History matching of typical well.


## 6. Conclusion

In this study, an analytical model was presented to simulate shale gas flow through MFHW in a shale reservoir. Based on microscale percolation mechanisms in nanopores, a seven-linear-flow model for dual-porosity formations was established. This model takes viscous flow, the Knudsen diffusion, the sorption-induced swelling response, the adsorption/desorption, the matrix shrinkage and the stress sensitivity effect into consideration. The pressure transient analysis was carried out under conditions of constant rate. The following conclusions were deduced:(
1
) The shale gas transport mechanism in nanopores mainly includes bulk gas transport and gas adsorption/desorption. Bulk gas apparent permeability is a weighted summation of viscous flow permeability and the Knudsen diffusion permeability. Gas desorption can increase the total apparent permeability.(
2
) A semianalytical model was established for the pressure behavior of MFHWs in shale gas reservoirs. Based on the proposed model, flow regimes were identified into two top reservoir regions, two inner-reservoir regions, two outer-reservoir regions, and one artificial fracture region.(
3
) Larger fracture half-length results in a more significant pressure drop and a higher production rate. Fracture permeability dominates throughout the production period except during the early pure wellbore storage stage and the boundary flow stage. The effect of desorption cannot be neglected when simulating, as the Langmuir volume has a significant influence on the production rate. Compared with the matrix radius and matrix permeability, pressure and production are more sensitive to induced fracture permeability.(
4
) The semianalytical model fits relatively well with the actual production data, has a reliable theoretical foundation, and can preferably describe the dynamic changes of pressure in the exploration process of shale gas wells.



---

*Source: 8065949-2018-01-17.xml*