# A Semianalytical Production Prediction Model and Dynamics Performance Analysis for Shale Gas Wells

Authors: Yingying Xu; Xiangui Liu; Zhiming Hu; Xianggang Duan; Jin Chang
Journal: Geofluids
 (2021)

---

## Abstract

The shale gas productivity model based on shale gas nonlinear seepage mechanism is an effective way to reasonably predict productivity. The incomplete gas nonlinear effects considered in the current production prediction models can lead to inaccurate production prediction. Based on the conventional five-zone compound flow model, comprehensive gas nonlinearities were considered in the improved compound linear flow model proposed in the paper and a semianalytical solution for productivity was obtained. The reliability of the productivity model was verified by the field data, and then, the 20-year production performance analysis of the gas well was studied. Ultimately, the key influencing factors of the fracture control stage and matrix control stage have been analyzed. Research indicated the following: (1) the EUR predicted by the productivity model is higher than the EUR that the comprehensive nonlinear effects are not considered, which demonstrated that the various nonlinear effects cannot be neglected during the production prediction to ensure the greater calculation accuracy; (2) during the early production stage of shale reservoir, the adsorbed gas is basically not recovered, and the cumulative adsorption contribution rate does not exceed 10%. The final adsorption gas contribution rate is 23.28%, and the annual adsorption rate can exceed 50% in the 20th year, showing that free gas and adsorbed gas are, respectively, important sources of the early stage of production and long-term stable production; (3) the widely ranged three-dimensional fracturing reformation of shale reservoirs and reasonable bottom hole pressure in the later matrix development process should be implemented to increase the effective early production of the reservoir and ensure the earlier gas production process of the matrix development. The findings of this study can help for better ensuring the prediction accuracy of the estimated ultimate recovery and understanding the main influencing factors of the dynamic performance of gas wells so as to provide a theoretical reference for production optimization and development plan formulation of the shale gas reservoirs.


---

## Body

## 1. Introduction

Nowadays, shale gas has occupied an important strategic position in the energy structure of worldwide natural gas production for its rich reserves and energy efficiency [1]. Productivity is adopted to evaluate the development effect of shale gas reservoirs and has become a key research hotspot of shale gas fields utilization. The shale gas production evaluation methods mainly include the decline curve method and the productivity model method. The decline curve equation [2–4] and the modified decline equation [5, 6] are adopted to fit and regress the production performance data of shale gas wells. Therefore, the key uncertain parameters in the empirical equation can be inverted by adjusting predefined parameters to predict the production of shale gas wells. Nevertheless, the method is limited to the production decline characteristics on the constant pressure production of shale gas wells and is easily affected by the irregular stability of production data, which makes it difficult to achieve the validity of this model in the shale gas development practice. On the contrary, productivity model based on shale gas complex flow mechanism is an effective approach to reasonably analyze the production performance of shale gas wells, which can avoid the limitations of shale gas production empirical equations and is of theoretical significance and application value to formulate shale reservoir fracturing design plans and guide the exploration and development of shale gas reservoirs.Different from conventional gas reservoirs, the dynamic production of shale gas wells is contributed by complex nonlinear effects such as multiple flow mechanism (continuous flow, slip flow, transition flow, or molecular free flow), gas PVT properties, and high-pressure adsorption at the same time [7, 8]. Nowadays, most productivity models [9–15] are lack of considering the comprehensive multiscale nonlinear effects, such as the impact of the real-gas effect and the nonlinearity of the supercritical desorption law on shale gas productivity, resulting in low productivity evaluation. Therefore, a unified cognition of the establishment of the shale gas productivity model has not yet been formed.In this paper, considering the lack of the nonlinear factors of the existing production model, the workflow is developed as follows: first, the nonlinearity of real-gas physical properties, high-pressure desorption of adsorbed gas, the non-Darcy flow mechanism such as slippage effect and diffusion, and the stress sensitivity of the secondary fracture network are comprehensively considered in the production prediction model. Subsequently, the semi-analytical solution of the productivity model is obtained based on Laplace transformation, canonical perturbation transformation, Stehfest numerical inversion, and iterative methods. Then, the production data of a shale gas well verifies the reliability of the productivity solution to confirm the importance of shale gas comprehensive nonlinear effects during the shale gas production. Next, the production performance of the shale gas well for 20 years can be illustrated. Finally, the key production influence factors of gas wells have been clarified, and specific guidelines for production optimization have been put forward, which is of certain theoretical reference to increase gas reservoir production rate and final gas recovery.


## 2. Methodology

### 2.1. Physical Model Description

Due to the irregular distribution of secondary fractures in the stimulated reservoir zone, the specific fracture heterogeneity characteristic parameters cannot be accurately obtained, and the fracture network is simplified to facilitate the calculation process. It is assumed that all the hydraulic fractures can be characterized by biwing transverse fracture [16] in this study and the secondary fracture network and matrix are treated into equivalent medium. Considering that there are some unstimulated reservoir zones between the hydraulic fractures, the physical model can be divided into the following five flow areas: hydraulic fractures, fracture network area 1, matrix area 2 between hydraulic fractures, and unstimulated matrix area 3 and matrix area 4. The flow behaviors in different areas are obviously different. A quarter of a single hydraulic fracture is selected in this model as illustrated in Figure 1:Figure 1
Simplified five-zone composite physical model of a multifractured horizontal shale gas well.The specific assumptions are as follows:(1)
The migration of single-phase methane gas in different seepage fields is considered as the one-dimensional flow. Shale gas flows into the inner zone 1 and 2 from the outer zone 3 and zone 4, respectively. The gas in the zone 1 from the zone 2 converges to the hydraulic fracture and the wellbore successively. The shale gas reservoir is of uniform thickness and the outer boundary is closed. The outer boundary distance is defined asye(2)
The adsorption characteristics of the adsorbed gas in the matrix 2, 3, and 4 can be explained by the supercritical Langmuir adsorption equation. The diffusion and slippage phenomenon of free gas in the matrix are considered(3)
The free gas in the inner fracture network 1 flows to the hydraulic fracture by viscous flow considering the stress sensitivity of the fracture network(4)
The hydraulic fractures in the inner zone are uniformly equal-length and symmetrically distributed, which fully penetrate the reservoir vertically. The half-length, width, cluster spacing of a single hydraulic fracture and the fractured half-width areyF, w, LF, d/2. The gas flow behavior in the hydraulic fractures follows Darcy flow law(5)
The gas flow mass exchange between different seepage areas is an unsteady flow behavior, and the gas production process is isothermal seepage. The gravity and capillary force are ignored, and the gas well is produced at constant pressure


### 2.2. Mathematical Model

#### 2.2.1. Governing Equation for Gas Flow in Matrix of Zone 4

The gas in the matrix of zone 4 flows into the matrix of zone 2 along they direction. Considering the supercritical desorption of adsorbed gas in the shale matrix, the diffusion, and slippage of free gas, the outer boundary is closed and the pressure on the inner boundary is continuous.According to the law of mass conservation, the governing equation for gas flow in matrix [1] can be obtained:
(1)ρgφ4Ct4∂P4∂t=∂∂yρgK4aμ4∂P4∂y,where the matrix comprehensive compression factor is
(2)Ct4=Cg4+Cd4+Cf4.Gas compression factor is(3)Cg4=1ρg∂ρg∂P4.As the conventional Langmuir equation cannot fit the high temperature and high-pressure isotherm adsorption law, the high-pressure isotherm adsorption model [17] is adopted, and desorption gas compression factor [18]is transformed as equation (4):
(4)Cd4=VL1−φ4TZPscφ4ZscTscPLP4P4+PL21−ρgρa.Javadpour apparent permeability model [19] is used to describe the gas diffusion and slip in the matrix:
(5)Kia=cDgiμi+FiKi,whereslip speed correction factor is(6)Fi=1+8πRTMμi×10−3Pavg×106ri2α−1.The hydraulic flow radius is(7)ri=8Kiφi.The pseudopressure and pseudotime [18]are used to linearize the high-pressure physical property parameters in the governing equation so as to solve the equations easily.The pseudopressure is(8)ψ=∫0P2PμZdP.The pseudotime is(9)ta=∫0tμiCtiμPCtPdt.According to the definition of pseudopressure and pseudotime and the definition of the dimensionless parameters in Table1, Eq. (1) can be rewritten as follows:
(10)∂ψ4D∂ta=1η4D∂2ψ4D∂y2.Table 1
Dimensionless parameters definition.
Dimensionless parameterDefinitionDimensionless parameterDefinitionDimensionless pseudopressureψD=ψe−ψ/ψi−ψwfDimensionless production1qD=TscKFAcwψe−ψwfPscqscTDimensionless pseudotimetaD=KFtaμφ1Ct1+φ2Ct2+φFCtFAcwStorage capacity ratiowi=ψiCtiφ1Ct1+φ2Ct2+φFCtFDimensionless lengthxD=2xLfyD=yAcwZone3-1mass transfer coefficientλ13=12K1LF2K3AcwZone2-4Mass transfer coefficientλ24=12K4LF2K2AcwZone1-Fmass transfer coefficientλ1F=12K1LF2KFAcwDimensionless conductivity in zone 3η3D=KFφ1Ct1+φ2Ct2+φFCtFφ3Ct3K3aDimensionless conductivity in zone 4η3D=KFφ1Ct1+φ2Ct2+φFCtFφ4Ct4K4aDimensionless pressureγfD=γfψe−ψwfDimensionless formation conductivityRCD=K1dK2LFThe initial matrix pressure of zone 4 is equal to the original formation pressure:(11)P4t=0=Pe,(12)Outerboundary:∂P4∂yy=ye=0,(13)Innerboundary:P4y=yF=P2y=yF.The dimensionless forms of Eq. ((11)–(13)) are as follows:
(14)ψ4DyD,0=0,ψ4DyFD,taD=ψ2D,∂ψ4DyD,taD∂yDyD=yeD=0.Then, the pseudopressure in the dimensionless flow equations can be transformed as equation (15) in Laplace form [20]:
(15)ζLD=LζD=∫0∞ζDe−stDdtD.Next, both sides of the dimensionless flow control equation (10) are multiplied by e−stDsimultaneously and integrated by the dimensionless time from zero to infinity, which turns out as the equation (16):
(16)∫0∞e−staD∂ζ4D∂taDdtaD=1η4D∫0∞e−staD∂2ζ4D∂yD2dtaD.The integral term at the left side of equation (16) is simplified as equation (17):
(17)∫0∞e−staD∂ζ4D∂taDdtaD=∫0∞e−staDdζ4D=e−staDζ4D0∞+s∫0∞e−staDζ4DdtaD=sζL4D.The integral term at the right side of equation (16) is simplified as equation (18):
(18)∫0∞e−staD∂2ζ4D∂yD2dtaD=∂2∫0∞e−staDζ4DdtaD∂yD2=∂2ζL4DyD,s∂yD2.Eq. (17) and Eq. (18) are substituted into Eq. (16) to obtain the dimensionless flow equation (19) of the outer matrix 4 in the Laplace space:
(19)∂2ζL4D∂yD2=sη4DζL4D.Boundary conditions can be expressed in Laplace form as equations ((20), (21)):
(20)ζL4DyFD,s=ζL2D,(21)∂ζL4DyD,s∂yDyD=yeD=0.The dimensionless pseudopressure solution of zone 4 in Laplace space is expressed in equation (22) based on equation (19) and boundary conditions ((20), (21)):
(22)ζL4DyD,s=ζL2Dcoshsη4DyeD−yDcoshsη4DyeD−yFD.


#### 2.2.2. Governing Equation for Gas Flow in the Matrix of Zone 3

Considering the gas desorption, diffusion, and slippage from matrix zone 3 into zone 1 along they direction, the dimensionless gas percolation equation in the matrix of zone 3 can be derived in equation (23):
(23)∂ψ3D∂taD=1η3D∂2ψ3D∂yD2.As the outer boundary is closed and the pressure on the inner boundary is continuous, the dimensionless initial and boundary conditions are reflected in equations ((24)–(26)) as follows:
(24)ψ3DyD,0=0,(25)ψ3DyFD,taD=ψ1D,(26)∂ψ3DyD,taD∂yDyD=yeD=0.Similar to the solution process of dimensionless flow equation in the matrix of zone4, the dimensionless pseudo pressure solution of matrix zone 3 in Laplace space can be derived in equation (27) based on Laplace transformation:
(27)ζL3DyD,s=ζL1Dcoshsη3DyeD−yDcoshsη3DyeD−yFD.


#### 2.2.3. Governing Equation for Gas Flow in Matrix of Zone 2

Considering the unsteady gas flow exchange between matrix zone 4 and matrix zone 2, the dimensionless gas seepage equation in the matrix of zone 2 is established in equation (28):
(28)∂2ψ2D∂xD2=3w2λ12∂ψ2D∂taD−6λ24yFD∂ψ4D∂yDy=yFD.The outer boundary is closed, the inner boundary pressure is continuous, and the initial and boundary conditions are expressed in equations ((29)–(31)):The initial condition is(29)ψ2DxD,0=0.Outer boundary is(30)∂ψ2DxD,taD∂xDxD=1=0.Inner boundary is(31)ψ2DdD2,taD=ψ1D.The dimensionless pseudopressure solution of the matrix in zone 2 in Laplace space can be deduced in equation (32) based on Laplace transformation:
(32)ζL2DxD,s=ζL1Dcoshα21−xDcoshα21−dD/2,where
(33)α2=3w2sλ12+6sη4Dλ24yFDtanhsη4DyeD−yFD.


#### 2.2.4. Governing Equation for Gas Flow in Fracture Network of Zone 1

Considering the impact of fracture network stress sensitivity on gas viscous flow in zone 1, the flow rate on the outer boundary between zone 1 and zone 2 is continuous, and the pressure on the inner boundary is continuous.An exponential stress sensitivity empirical model [21] is adopted in this paper to characterize the impact of the stress sensitivity of secondary fracture on gas flow, which is expressed in equation (34):
(34)Kf=Kfie−γψi−ψf.The dimensionless gas flow equation in fracture network zone 1 can be depicted in expression (35):
(35)∂2ψ1D∂xD2−γ1D∂ψ1D∂xD2=3w1λ1Feγ1Dψ1D∂ψ1D∂taD−6eγ1Dψ1Dλ13yFD∂ψ3D∂yDy=yFD.The boundary conditions are as followed in equations ((36)–(38)):Initial condition is(36)ψ1DxD,0=0.Outer boundary is(37)ψ1D0,taD=ψFD.Inner boundary is(38)∂ψ1DxD,taD∂xDxD=dD/2=dDeγ1Dψ1D2RCD∂ψ2DxD,taD∂xDxD=dD/2.Due to the strong nonlinearity of the stress-sensitivity term in the dimensionless gas flow equation, the pedrosa variable [22] is adopted to replace the pseudopressure term. The variable term expression is
(39)ψD=−1γDLn1−γD⋅ζD,where ζD is the perturbation transformation function.The equation (39) is substituted into the dimensionless equation (35), and equation (35) can be derived into equation (40) in zone 1:
(40)∂2ζ1∂xD2=3w1λ1F⋅∂ζ1∂taD−6λ13yFD⋅11−γ1Dζ1D∂ζ3∂yD.The term1/1−γ1Dζ1D in equation (40) can be expanded into Taylor form in equation (41):
(41)11−γ1Dζ1D=1+γ1Dζ1D+γ1Dζ1D2+γ1Dζ1D3+.⋯⋯Because the dimensionless pseudopermeability modulusγ1D ≪1, the term γ1Dζ1D can be considered as zero approximately; then, the equation (40) can be simplified as equation (42):
(42)∂2ζ1D∂xD2=3w1λ1F⋅∂ζ1D∂taD−6λ13yFD∂ζ3D∂yD.The initial and boundary conditions are simplified as follows:Initial conditions are(43)ζ1DxD,0=0.Outer boundary is(44)ζ1D0,taD=ζFD1−e−γ1Dγ1D.Inner boundary is(45)∂ζ1DxD,taD∂xDxD=dD/2=dD2RCD∂ζ2DxD,taD∂xDxD=dD/2.Similarly, combining the Laplace transform, the pseudopressure solution of equation (42) can be expressed in equation (46):
(46)ζL1DxD,s=ζLF1−e−γ1Dγ1Dc3sinhα1xD−dD/2+coshα1xD−dD/2−c3sinhα1dD/2+coshα1dD/2,where
(47)c3=λ12λ1Fα2α1tanhα2dD2−1,α1=3w1sλ1F+6sη3Dλ13yFDtanhsη3DyeD−yFD.


#### 2.2.5. Governing Equation for Gas Flow in Inner Zone Hydraulic Fracture

The gas flow in the hydraulic fracture is linear flow, and the dimensionless gas flow governing equation (48) in the hydraulic fracture is established:
(48)∂2ψFD∂yD2=wF∂ψFD∂taD−λ1Feγ1Dψ1D3∂ψ1D∂xDxD=dD/2.Considering the constant bottom hole pressure of the gas well and the outer boundary is closed, the initial and boundary conditions are as follows:The initial condition is(49)ψFDyD,0=0.Outer boundary is(50)∂ψFDyD,taD∂taDyD=yFD=0.Inner boundary is(51)ψFD0,taD=1.Then, the dimensionless governing equations (48) of the hydraulic fracture can be simplified and solved by using pedrosa variable substitution and Laplace transform in equation (52):
(52)ζLFDyD,s=coshα4yD−yFDcoshα4yFD,where
(53)α4=wFs−λ1Fα33γfD1−e−γfD,α3=α1c3coshα1xD−dD/2+α1sinhα1xD−dD/2−c3sinhα1dD/2+coshα1dD/2.Due to the number of hydraulic fractures for the gas well, the total dimensionless production rate at the bottom hole of the shale gas well in Laplace space can be derived in equation (54):
(54)qLD=−N2π∂ζLFD∂yDyD=0=Nα42πstanhα4·yFD.Based on the Stehfest numerical inversion [23], the semianalytical solution of the dimensionless production in the real space is obtained. Then, combining with the Newton iteration method, the production solution in the real space at constant pressure can be derived.


### 2.3. Model Verification and Analysis

The production data of a multistage fractured horizontal well in the Weiyuan shale gas block can be adopted to testify the validity of the proposed production model. The relevant geological parameters and horizontal well parameters are shown in the following Table2. By adjusting the relevant parameters of the model, the semianalytical model proposed in this paper can historically match the dynamic production data of the gas well. (∗ indicates model fitting parameters).Table 2
Reservoir geological parameters and completion parameters.
Parameters nameValueParameters nameValueInitial formation pressure/MPa45The stimulated fracture stages25Formation temperature/K392.4∗The hydraulic fracture half-length/m100Formation thickness/m30∗ The hydraulic fracture width/m0.003Matrix porosity0.038∗eternal distance/m300Matrix permeability/mD0.000447Horizontal well length/m1500Fracture network porosity0.032Bottom hole flow pressure /MPa5Fracture network permeability/mD0.851Rock density/(kg/m3)2600Hydraulic fracture permeability /(mD)750Stress sensitivity coefficient /MPa-10.13Langmuir volume/(m3/t)3.7Langmuir pressure /MPa8.5∗The hydraulic fracture zone half width/m5Fracture cluster/stage3It can be seen from Figure2 that the production process of this gas well lasted 1140 days, and the daily gas production data was quite volatile, but the overall monotonous decline law was obvious. At 800 days, the daily gas production rate dropped from 230,000 m3/day to 10,000 m3/day, and the cumulative gas production was 71 million cubic meters, basically entering the stable production stage of the gas well. Obviously, the gas production curve calculated by the productivity model is basically consistent with the typical production curve. Furthermore, the cumulative gas production calculated by the model is 77.9 million cubic meters, and the actual gas production for 1140 days of the gas well reaches 76.4 million cubic meters, with a relative error of 1.986%, so the reliability of the model can be considered to be testified. As is shown in Figure 3, the contribution of nonlinear effects to production is not obvious in the first year of the gas well production but cannot be ignored with the increase of the production time. Then, the EUR of the gas well for 20 years is 137.9 million cubic meters, which is 23.46% higher than the 111.7 million calculated without considering the comprehensive nonlinear effect, which shows that the production can be underestimated seriously without considering any of the abovementioned nonlinear effects. This phenomenon is attributed to the thinning effect of gas molecules in the micronano pores of the matrix in the low-pressure reservoir at the later stage of shale gas development. The collision frequency between the gas molecules and the pore wall becomes the dominant kinematic mechanism of gas motivation, so microscopic mass transfer methods such as Knudsen diffusion, desorption, and slip flow are beginning to play a significant part in enhancing gas permeability greatly, which can delay the deceleration rate of daily gas production and strengthen the gas supply degree in the later stage. The semianalytical model closely relies on on-site production information and considers more multinonlinearities, which greatly improves the calculation accuracy so as to reasonably evaluate the medium and long-term productivity of gas wells and predict the EUR value.Figure 2
Model verification results.
(a)
Fitting curve of daily gas production(b)
Fitting curve of cumulative gas productionFigure 3
Prediction results of gas production.
(a)
The model calculation result of the daily gas production(b)
The model calculation result of the cumulative gas production


## 2.1. Physical Model Description

Due to the irregular distribution of secondary fractures in the stimulated reservoir zone, the specific fracture heterogeneity characteristic parameters cannot be accurately obtained, and the fracture network is simplified to facilitate the calculation process. It is assumed that all the hydraulic fractures can be characterized by biwing transverse fracture [16] in this study and the secondary fracture network and matrix are treated into equivalent medium. Considering that there are some unstimulated reservoir zones between the hydraulic fractures, the physical model can be divided into the following five flow areas: hydraulic fractures, fracture network area 1, matrix area 2 between hydraulic fractures, and unstimulated matrix area 3 and matrix area 4. The flow behaviors in different areas are obviously different. A quarter of a single hydraulic fracture is selected in this model as illustrated in Figure 1:Figure 1
Simplified five-zone composite physical model of a multifractured horizontal shale gas well.The specific assumptions are as follows:(1)
The migration of single-phase methane gas in different seepage fields is considered as the one-dimensional flow. Shale gas flows into the inner zone 1 and 2 from the outer zone 3 and zone 4, respectively. The gas in the zone 1 from the zone 2 converges to the hydraulic fracture and the wellbore successively. The shale gas reservoir is of uniform thickness and the outer boundary is closed. The outer boundary distance is defined asye(2)
The adsorption characteristics of the adsorbed gas in the matrix 2, 3, and 4 can be explained by the supercritical Langmuir adsorption equation. The diffusion and slippage phenomenon of free gas in the matrix are considered(3)
The free gas in the inner fracture network 1 flows to the hydraulic fracture by viscous flow considering the stress sensitivity of the fracture network(4)
The hydraulic fractures in the inner zone are uniformly equal-length and symmetrically distributed, which fully penetrate the reservoir vertically. The half-length, width, cluster spacing of a single hydraulic fracture and the fractured half-width areyF, w, LF, d/2. The gas flow behavior in the hydraulic fractures follows Darcy flow law(5)
The gas flow mass exchange between different seepage areas is an unsteady flow behavior, and the gas production process is isothermal seepage. The gravity and capillary force are ignored, and the gas well is produced at constant pressure


## 2.2. Mathematical Model

### 2.2.1. Governing Equation for Gas Flow in Matrix of Zone 4

The gas in the matrix of zone 4 flows into the matrix of zone 2 along they direction. Considering the supercritical desorption of adsorbed gas in the shale matrix, the diffusion, and slippage of free gas, the outer boundary is closed and the pressure on the inner boundary is continuous.According to the law of mass conservation, the governing equation for gas flow in matrix [1] can be obtained:
(1)ρgφ4Ct4∂P4∂t=∂∂yρgK4aμ4∂P4∂y,where the matrix comprehensive compression factor is
(2)Ct4=Cg4+Cd4+Cf4.Gas compression factor is(3)Cg4=1ρg∂ρg∂P4.As the conventional Langmuir equation cannot fit the high temperature and high-pressure isotherm adsorption law, the high-pressure isotherm adsorption model [17] is adopted, and desorption gas compression factor [18]is transformed as equation (4):
(4)Cd4=VL1−φ4TZPscφ4ZscTscPLP4P4+PL21−ρgρa.Javadpour apparent permeability model [19] is used to describe the gas diffusion and slip in the matrix:
(5)Kia=cDgiμi+FiKi,whereslip speed correction factor is(6)Fi=1+8πRTMμi×10−3Pavg×106ri2α−1.The hydraulic flow radius is(7)ri=8Kiφi.The pseudopressure and pseudotime [18]are used to linearize the high-pressure physical property parameters in the governing equation so as to solve the equations easily.The pseudopressure is(8)ψ=∫0P2PμZdP.The pseudotime is(9)ta=∫0tμiCtiμPCtPdt.According to the definition of pseudopressure and pseudotime and the definition of the dimensionless parameters in Table1, Eq. (1) can be rewritten as follows:
(10)∂ψ4D∂ta=1η4D∂2ψ4D∂y2.Table 1
Dimensionless parameters definition.
Dimensionless parameterDefinitionDimensionless parameterDefinitionDimensionless pseudopressureψD=ψe−ψ/ψi−ψwfDimensionless production1qD=TscKFAcwψe−ψwfPscqscTDimensionless pseudotimetaD=KFtaμφ1Ct1+φ2Ct2+φFCtFAcwStorage capacity ratiowi=ψiCtiφ1Ct1+φ2Ct2+φFCtFDimensionless lengthxD=2xLfyD=yAcwZone3-1mass transfer coefficientλ13=12K1LF2K3AcwZone2-4Mass transfer coefficientλ24=12K4LF2K2AcwZone1-Fmass transfer coefficientλ1F=12K1LF2KFAcwDimensionless conductivity in zone 3η3D=KFφ1Ct1+φ2Ct2+φFCtFφ3Ct3K3aDimensionless conductivity in zone 4η3D=KFφ1Ct1+φ2Ct2+φFCtFφ4Ct4K4aDimensionless pressureγfD=γfψe−ψwfDimensionless formation conductivityRCD=K1dK2LFThe initial matrix pressure of zone 4 is equal to the original formation pressure:(11)P4t=0=Pe,(12)Outerboundary:∂P4∂yy=ye=0,(13)Innerboundary:P4y=yF=P2y=yF.The dimensionless forms of Eq. ((11)–(13)) are as follows:
(14)ψ4DyD,0=0,ψ4DyFD,taD=ψ2D,∂ψ4DyD,taD∂yDyD=yeD=0.Then, the pseudopressure in the dimensionless flow equations can be transformed as equation (15) in Laplace form [20]:
(15)ζLD=LζD=∫0∞ζDe−stDdtD.Next, both sides of the dimensionless flow control equation (10) are multiplied by e−stDsimultaneously and integrated by the dimensionless time from zero to infinity, which turns out as the equation (16):
(16)∫0∞e−staD∂ζ4D∂taDdtaD=1η4D∫0∞e−staD∂2ζ4D∂yD2dtaD.The integral term at the left side of equation (16) is simplified as equation (17):
(17)∫0∞e−staD∂ζ4D∂taDdtaD=∫0∞e−staDdζ4D=e−staDζ4D0∞+s∫0∞e−staDζ4DdtaD=sζL4D.The integral term at the right side of equation (16) is simplified as equation (18):
(18)∫0∞e−staD∂2ζ4D∂yD2dtaD=∂2∫0∞e−staDζ4DdtaD∂yD2=∂2ζL4DyD,s∂yD2.Eq. (17) and Eq. (18) are substituted into Eq. (16) to obtain the dimensionless flow equation (19) of the outer matrix 4 in the Laplace space:
(19)∂2ζL4D∂yD2=sη4DζL4D.Boundary conditions can be expressed in Laplace form as equations ((20), (21)):
(20)ζL4DyFD,s=ζL2D,(21)∂ζL4DyD,s∂yDyD=yeD=0.The dimensionless pseudopressure solution of zone 4 in Laplace space is expressed in equation (22) based on equation (19) and boundary conditions ((20), (21)):
(22)ζL4DyD,s=ζL2Dcoshsη4DyeD−yDcoshsη4DyeD−yFD.


### 2.2.2. Governing Equation for Gas Flow in the Matrix of Zone 3

Considering the gas desorption, diffusion, and slippage from matrix zone 3 into zone 1 along they direction, the dimensionless gas percolation equation in the matrix of zone 3 can be derived in equation (23):
(23)∂ψ3D∂taD=1η3D∂2ψ3D∂yD2.As the outer boundary is closed and the pressure on the inner boundary is continuous, the dimensionless initial and boundary conditions are reflected in equations ((24)–(26)) as follows:
(24)ψ3DyD,0=0,(25)ψ3DyFD,taD=ψ1D,(26)∂ψ3DyD,taD∂yDyD=yeD=0.Similar to the solution process of dimensionless flow equation in the matrix of zone4, the dimensionless pseudo pressure solution of matrix zone 3 in Laplace space can be derived in equation (27) based on Laplace transformation:
(27)ζL3DyD,s=ζL1Dcoshsη3DyeD−yDcoshsη3DyeD−yFD.


### 2.2.3. Governing Equation for Gas Flow in Matrix of Zone 2

Considering the unsteady gas flow exchange between matrix zone 4 and matrix zone 2, the dimensionless gas seepage equation in the matrix of zone 2 is established in equation (28):
(28)∂2ψ2D∂xD2=3w2λ12∂ψ2D∂taD−6λ24yFD∂ψ4D∂yDy=yFD.The outer boundary is closed, the inner boundary pressure is continuous, and the initial and boundary conditions are expressed in equations ((29)–(31)):The initial condition is(29)ψ2DxD,0=0.Outer boundary is(30)∂ψ2DxD,taD∂xDxD=1=0.Inner boundary is(31)ψ2DdD2,taD=ψ1D.The dimensionless pseudopressure solution of the matrix in zone 2 in Laplace space can be deduced in equation (32) based on Laplace transformation:
(32)ζL2DxD,s=ζL1Dcoshα21−xDcoshα21−dD/2,where
(33)α2=3w2sλ12+6sη4Dλ24yFDtanhsη4DyeD−yFD.


### 2.2.4. Governing Equation for Gas Flow in Fracture Network of Zone 1

Considering the impact of fracture network stress sensitivity on gas viscous flow in zone 1, the flow rate on the outer boundary between zone 1 and zone 2 is continuous, and the pressure on the inner boundary is continuous.An exponential stress sensitivity empirical model [21] is adopted in this paper to characterize the impact of the stress sensitivity of secondary fracture on gas flow, which is expressed in equation (34):
(34)Kf=Kfie−γψi−ψf.The dimensionless gas flow equation in fracture network zone 1 can be depicted in expression (35):
(35)∂2ψ1D∂xD2−γ1D∂ψ1D∂xD2=3w1λ1Feγ1Dψ1D∂ψ1D∂taD−6eγ1Dψ1Dλ13yFD∂ψ3D∂yDy=yFD.The boundary conditions are as followed in equations ((36)–(38)):Initial condition is(36)ψ1DxD,0=0.Outer boundary is(37)ψ1D0,taD=ψFD.Inner boundary is(38)∂ψ1DxD,taD∂xDxD=dD/2=dDeγ1Dψ1D2RCD∂ψ2DxD,taD∂xDxD=dD/2.Due to the strong nonlinearity of the stress-sensitivity term in the dimensionless gas flow equation, the pedrosa variable [22] is adopted to replace the pseudopressure term. The variable term expression is
(39)ψD=−1γDLn1−γD⋅ζD,where ζD is the perturbation transformation function.The equation (39) is substituted into the dimensionless equation (35), and equation (35) can be derived into equation (40) in zone 1:
(40)∂2ζ1∂xD2=3w1λ1F⋅∂ζ1∂taD−6λ13yFD⋅11−γ1Dζ1D∂ζ3∂yD.The term1/1−γ1Dζ1D in equation (40) can be expanded into Taylor form in equation (41):
(41)11−γ1Dζ1D=1+γ1Dζ1D+γ1Dζ1D2+γ1Dζ1D3+.⋯⋯Because the dimensionless pseudopermeability modulusγ1D ≪1, the term γ1Dζ1D can be considered as zero approximately; then, the equation (40) can be simplified as equation (42):
(42)∂2ζ1D∂xD2=3w1λ1F⋅∂ζ1D∂taD−6λ13yFD∂ζ3D∂yD.The initial and boundary conditions are simplified as follows:Initial conditions are(43)ζ1DxD,0=0.Outer boundary is(44)ζ1D0,taD=ζFD1−e−γ1Dγ1D.Inner boundary is(45)∂ζ1DxD,taD∂xDxD=dD/2=dD2RCD∂ζ2DxD,taD∂xDxD=dD/2.Similarly, combining the Laplace transform, the pseudopressure solution of equation (42) can be expressed in equation (46):
(46)ζL1DxD,s=ζLF1−e−γ1Dγ1Dc3sinhα1xD−dD/2+coshα1xD−dD/2−c3sinhα1dD/2+coshα1dD/2,where
(47)c3=λ12λ1Fα2α1tanhα2dD2−1,α1=3w1sλ1F+6sη3Dλ13yFDtanhsη3DyeD−yFD.


### 2.2.5. Governing Equation for Gas Flow in Inner Zone Hydraulic Fracture

The gas flow in the hydraulic fracture is linear flow, and the dimensionless gas flow governing equation (48) in the hydraulic fracture is established:
(48)∂2ψFD∂yD2=wF∂ψFD∂taD−λ1Feγ1Dψ1D3∂ψ1D∂xDxD=dD/2.Considering the constant bottom hole pressure of the gas well and the outer boundary is closed, the initial and boundary conditions are as follows:The initial condition is(49)ψFDyD,0=0.Outer boundary is(50)∂ψFDyD,taD∂taDyD=yFD=0.Inner boundary is(51)ψFD0,taD=1.Then, the dimensionless governing equations (48) of the hydraulic fracture can be simplified and solved by using pedrosa variable substitution and Laplace transform in equation (52):
(52)ζLFDyD,s=coshα4yD−yFDcoshα4yFD,where
(53)α4=wFs−λ1Fα33γfD1−e−γfD,α3=α1c3coshα1xD−dD/2+α1sinhα1xD−dD/2−c3sinhα1dD/2+coshα1dD/2.Due to the number of hydraulic fractures for the gas well, the total dimensionless production rate at the bottom hole of the shale gas well in Laplace space can be derived in equation (54):
(54)qLD=−N2π∂ζLFD∂yDyD=0=Nα42πstanhα4·yFD.Based on the Stehfest numerical inversion [23], the semianalytical solution of the dimensionless production in the real space is obtained. Then, combining with the Newton iteration method, the production solution in the real space at constant pressure can be derived.


## 2.2.1. Governing Equation for Gas Flow in Matrix of Zone 4

The gas in the matrix of zone 4 flows into the matrix of zone 2 along they direction. Considering the supercritical desorption of adsorbed gas in the shale matrix, the diffusion, and slippage of free gas, the outer boundary is closed and the pressure on the inner boundary is continuous.According to the law of mass conservation, the governing equation for gas flow in matrix [1] can be obtained:
(1)ρgφ4Ct4∂P4∂t=∂∂yρgK4aμ4∂P4∂y,where the matrix comprehensive compression factor is
(2)Ct4=Cg4+Cd4+Cf4.Gas compression factor is(3)Cg4=1ρg∂ρg∂P4.As the conventional Langmuir equation cannot fit the high temperature and high-pressure isotherm adsorption law, the high-pressure isotherm adsorption model [17] is adopted, and desorption gas compression factor [18]is transformed as equation (4):
(4)Cd4=VL1−φ4TZPscφ4ZscTscPLP4P4+PL21−ρgρa.Javadpour apparent permeability model [19] is used to describe the gas diffusion and slip in the matrix:
(5)Kia=cDgiμi+FiKi,whereslip speed correction factor is(6)Fi=1+8πRTMμi×10−3Pavg×106ri2α−1.The hydraulic flow radius is(7)ri=8Kiφi.The pseudopressure and pseudotime [18]are used to linearize the high-pressure physical property parameters in the governing equation so as to solve the equations easily.The pseudopressure is(8)ψ=∫0P2PμZdP.The pseudotime is(9)ta=∫0tμiCtiμPCtPdt.According to the definition of pseudopressure and pseudotime and the definition of the dimensionless parameters in Table1, Eq. (1) can be rewritten as follows:
(10)∂ψ4D∂ta=1η4D∂2ψ4D∂y2.Table 1
Dimensionless parameters definition.
Dimensionless parameterDefinitionDimensionless parameterDefinitionDimensionless pseudopressureψD=ψe−ψ/ψi−ψwfDimensionless production1qD=TscKFAcwψe−ψwfPscqscTDimensionless pseudotimetaD=KFtaμφ1Ct1+φ2Ct2+φFCtFAcwStorage capacity ratiowi=ψiCtiφ1Ct1+φ2Ct2+φFCtFDimensionless lengthxD=2xLfyD=yAcwZone3-1mass transfer coefficientλ13=12K1LF2K3AcwZone2-4Mass transfer coefficientλ24=12K4LF2K2AcwZone1-Fmass transfer coefficientλ1F=12K1LF2KFAcwDimensionless conductivity in zone 3η3D=KFφ1Ct1+φ2Ct2+φFCtFφ3Ct3K3aDimensionless conductivity in zone 4η3D=KFφ1Ct1+φ2Ct2+φFCtFφ4Ct4K4aDimensionless pressureγfD=γfψe−ψwfDimensionless formation conductivityRCD=K1dK2LFThe initial matrix pressure of zone 4 is equal to the original formation pressure:(11)P4t=0=Pe,(12)Outerboundary:∂P4∂yy=ye=0,(13)Innerboundary:P4y=yF=P2y=yF.The dimensionless forms of Eq. ((11)–(13)) are as follows:
(14)ψ4DyD,0=0,ψ4DyFD,taD=ψ2D,∂ψ4DyD,taD∂yDyD=yeD=0.Then, the pseudopressure in the dimensionless flow equations can be transformed as equation (15) in Laplace form [20]:
(15)ζLD=LζD=∫0∞ζDe−stDdtD.Next, both sides of the dimensionless flow control equation (10) are multiplied by e−stDsimultaneously and integrated by the dimensionless time from zero to infinity, which turns out as the equation (16):
(16)∫0∞e−staD∂ζ4D∂taDdtaD=1η4D∫0∞e−staD∂2ζ4D∂yD2dtaD.The integral term at the left side of equation (16) is simplified as equation (17):
(17)∫0∞e−staD∂ζ4D∂taDdtaD=∫0∞e−staDdζ4D=e−staDζ4D0∞+s∫0∞e−staDζ4DdtaD=sζL4D.The integral term at the right side of equation (16) is simplified as equation (18):
(18)∫0∞e−staD∂2ζ4D∂yD2dtaD=∂2∫0∞e−staDζ4DdtaD∂yD2=∂2ζL4DyD,s∂yD2.Eq. (17) and Eq. (18) are substituted into Eq. (16) to obtain the dimensionless flow equation (19) of the outer matrix 4 in the Laplace space:
(19)∂2ζL4D∂yD2=sη4DζL4D.Boundary conditions can be expressed in Laplace form as equations ((20), (21)):
(20)ζL4DyFD,s=ζL2D,(21)∂ζL4DyD,s∂yDyD=yeD=0.The dimensionless pseudopressure solution of zone 4 in Laplace space is expressed in equation (22) based on equation (19) and boundary conditions ((20), (21)):
(22)ζL4DyD,s=ζL2Dcoshsη4DyeD−yDcoshsη4DyeD−yFD.


## 2.2.2. Governing Equation for Gas Flow in the Matrix of Zone 3

Considering the gas desorption, diffusion, and slippage from matrix zone 3 into zone 1 along they direction, the dimensionless gas percolation equation in the matrix of zone 3 can be derived in equation (23):
(23)∂ψ3D∂taD=1η3D∂2ψ3D∂yD2.As the outer boundary is closed and the pressure on the inner boundary is continuous, the dimensionless initial and boundary conditions are reflected in equations ((24)–(26)) as follows:
(24)ψ3DyD,0=0,(25)ψ3DyFD,taD=ψ1D,(26)∂ψ3DyD,taD∂yDyD=yeD=0.Similar to the solution process of dimensionless flow equation in the matrix of zone4, the dimensionless pseudo pressure solution of matrix zone 3 in Laplace space can be derived in equation (27) based on Laplace transformation:
(27)ζL3DyD,s=ζL1Dcoshsη3DyeD−yDcoshsη3DyeD−yFD.


## 2.2.3. Governing Equation for Gas Flow in Matrix of Zone 2

Considering the unsteady gas flow exchange between matrix zone 4 and matrix zone 2, the dimensionless gas seepage equation in the matrix of zone 2 is established in equation (28):
(28)∂2ψ2D∂xD2=3w2λ12∂ψ2D∂taD−6λ24yFD∂ψ4D∂yDy=yFD.The outer boundary is closed, the inner boundary pressure is continuous, and the initial and boundary conditions are expressed in equations ((29)–(31)):The initial condition is(29)ψ2DxD,0=0.Outer boundary is(30)∂ψ2DxD,taD∂xDxD=1=0.Inner boundary is(31)ψ2DdD2,taD=ψ1D.The dimensionless pseudopressure solution of the matrix in zone 2 in Laplace space can be deduced in equation (32) based on Laplace transformation:
(32)ζL2DxD,s=ζL1Dcoshα21−xDcoshα21−dD/2,where
(33)α2=3w2sλ12+6sη4Dλ24yFDtanhsη4DyeD−yFD.


## 2.2.4. Governing Equation for Gas Flow in Fracture Network of Zone 1

Considering the impact of fracture network stress sensitivity on gas viscous flow in zone 1, the flow rate on the outer boundary between zone 1 and zone 2 is continuous, and the pressure on the inner boundary is continuous.An exponential stress sensitivity empirical model [21] is adopted in this paper to characterize the impact of the stress sensitivity of secondary fracture on gas flow, which is expressed in equation (34):
(34)Kf=Kfie−γψi−ψf.The dimensionless gas flow equation in fracture network zone 1 can be depicted in expression (35):
(35)∂2ψ1D∂xD2−γ1D∂ψ1D∂xD2=3w1λ1Feγ1Dψ1D∂ψ1D∂taD−6eγ1Dψ1Dλ13yFD∂ψ3D∂yDy=yFD.The boundary conditions are as followed in equations ((36)–(38)):Initial condition is(36)ψ1DxD,0=0.Outer boundary is(37)ψ1D0,taD=ψFD.Inner boundary is(38)∂ψ1DxD,taD∂xDxD=dD/2=dDeγ1Dψ1D2RCD∂ψ2DxD,taD∂xDxD=dD/2.Due to the strong nonlinearity of the stress-sensitivity term in the dimensionless gas flow equation, the pedrosa variable [22] is adopted to replace the pseudopressure term. The variable term expression is
(39)ψD=−1γDLn1−γD⋅ζD,where ζD is the perturbation transformation function.The equation (39) is substituted into the dimensionless equation (35), and equation (35) can be derived into equation (40) in zone 1:
(40)∂2ζ1∂xD2=3w1λ1F⋅∂ζ1∂taD−6λ13yFD⋅11−γ1Dζ1D∂ζ3∂yD.The term1/1−γ1Dζ1D in equation (40) can be expanded into Taylor form in equation (41):
(41)11−γ1Dζ1D=1+γ1Dζ1D+γ1Dζ1D2+γ1Dζ1D3+.⋯⋯Because the dimensionless pseudopermeability modulusγ1D ≪1, the term γ1Dζ1D can be considered as zero approximately; then, the equation (40) can be simplified as equation (42):
(42)∂2ζ1D∂xD2=3w1λ1F⋅∂ζ1D∂taD−6λ13yFD∂ζ3D∂yD.The initial and boundary conditions are simplified as follows:Initial conditions are(43)ζ1DxD,0=0.Outer boundary is(44)ζ1D0,taD=ζFD1−e−γ1Dγ1D.Inner boundary is(45)∂ζ1DxD,taD∂xDxD=dD/2=dD2RCD∂ζ2DxD,taD∂xDxD=dD/2.Similarly, combining the Laplace transform, the pseudopressure solution of equation (42) can be expressed in equation (46):
(46)ζL1DxD,s=ζLF1−e−γ1Dγ1Dc3sinhα1xD−dD/2+coshα1xD−dD/2−c3sinhα1dD/2+coshα1dD/2,where
(47)c3=λ12λ1Fα2α1tanhα2dD2−1,α1=3w1sλ1F+6sη3Dλ13yFDtanhsη3DyeD−yFD.


## 2.2.5. Governing Equation for Gas Flow in Inner Zone Hydraulic Fracture

The gas flow in the hydraulic fracture is linear flow, and the dimensionless gas flow governing equation (48) in the hydraulic fracture is established:
(48)∂2ψFD∂yD2=wF∂ψFD∂taD−λ1Feγ1Dψ1D3∂ψ1D∂xDxD=dD/2.Considering the constant bottom hole pressure of the gas well and the outer boundary is closed, the initial and boundary conditions are as follows:The initial condition is(49)ψFDyD,0=0.Outer boundary is(50)∂ψFDyD,taD∂taDyD=yFD=0.Inner boundary is(51)ψFD0,taD=1.Then, the dimensionless governing equations (48) of the hydraulic fracture can be simplified and solved by using pedrosa variable substitution and Laplace transform in equation (52):
(52)ζLFDyD,s=coshα4yD−yFDcoshα4yFD,where
(53)α4=wFs−λ1Fα33γfD1−e−γfD,α3=α1c3coshα1xD−dD/2+α1sinhα1xD−dD/2−c3sinhα1dD/2+coshα1dD/2.Due to the number of hydraulic fractures for the gas well, the total dimensionless production rate at the bottom hole of the shale gas well in Laplace space can be derived in equation (54):
(54)qLD=−N2π∂ζLFD∂yDyD=0=Nα42πstanhα4·yFD.Based on the Stehfest numerical inversion [23], the semianalytical solution of the dimensionless production in the real space is obtained. Then, combining with the Newton iteration method, the production solution in the real space at constant pressure can be derived.


## 2.3. Model Verification and Analysis

The production data of a multistage fractured horizontal well in the Weiyuan shale gas block can be adopted to testify the validity of the proposed production model. The relevant geological parameters and horizontal well parameters are shown in the following Table2. By adjusting the relevant parameters of the model, the semianalytical model proposed in this paper can historically match the dynamic production data of the gas well. (∗ indicates model fitting parameters).Table 2
Reservoir geological parameters and completion parameters.
Parameters nameValueParameters nameValueInitial formation pressure/MPa45The stimulated fracture stages25Formation temperature/K392.4∗The hydraulic fracture half-length/m100Formation thickness/m30∗ The hydraulic fracture width/m0.003Matrix porosity0.038∗eternal distance/m300Matrix permeability/mD0.000447Horizontal well length/m1500Fracture network porosity0.032Bottom hole flow pressure /MPa5Fracture network permeability/mD0.851Rock density/(kg/m3)2600Hydraulic fracture permeability /(mD)750Stress sensitivity coefficient /MPa-10.13Langmuir volume/(m3/t)3.7Langmuir pressure /MPa8.5∗The hydraulic fracture zone half width/m5Fracture cluster/stage3It can be seen from Figure2 that the production process of this gas well lasted 1140 days, and the daily gas production data was quite volatile, but the overall monotonous decline law was obvious. At 800 days, the daily gas production rate dropped from 230,000 m3/day to 10,000 m3/day, and the cumulative gas production was 71 million cubic meters, basically entering the stable production stage of the gas well. Obviously, the gas production curve calculated by the productivity model is basically consistent with the typical production curve. Furthermore, the cumulative gas production calculated by the model is 77.9 million cubic meters, and the actual gas production for 1140 days of the gas well reaches 76.4 million cubic meters, with a relative error of 1.986%, so the reliability of the model can be considered to be testified. As is shown in Figure 3, the contribution of nonlinear effects to production is not obvious in the first year of the gas well production but cannot be ignored with the increase of the production time. Then, the EUR of the gas well for 20 years is 137.9 million cubic meters, which is 23.46% higher than the 111.7 million calculated without considering the comprehensive nonlinear effect, which shows that the production can be underestimated seriously without considering any of the abovementioned nonlinear effects. This phenomenon is attributed to the thinning effect of gas molecules in the micronano pores of the matrix in the low-pressure reservoir at the later stage of shale gas development. The collision frequency between the gas molecules and the pore wall becomes the dominant kinematic mechanism of gas motivation, so microscopic mass transfer methods such as Knudsen diffusion, desorption, and slip flow are beginning to play a significant part in enhancing gas permeability greatly, which can delay the deceleration rate of daily gas production and strengthen the gas supply degree in the later stage. The semianalytical model closely relies on on-site production information and considers more multinonlinearities, which greatly improves the calculation accuracy so as to reasonably evaluate the medium and long-term productivity of gas wells and predict the EUR value.Figure 2
Model verification results.
(a)
Fitting curve of daily gas production(b)
Fitting curve of cumulative gas productionFigure 3
Prediction results of gas production.
(a)
The model calculation result of the daily gas production(b)
The model calculation result of the cumulative gas production


## 3. Dynamic Production Performance Analysis

The 20-year development performance of the gas well based on the productivity model proposed in this paper has been studied. It can be seen from Figures4 and 5 that when the average formation pressure drops from 45 MPa to 33 MPa at the production time of 1 year, the initial production is mainly contributed by free gas from fracture systems, while the adsorbed gas is basically not used, and the proportion of EUR does not exceed 40%. As the production time increases, due to the decline of reservoir pressure, the production of adsorbed gas gradually increases to result in the high cumulative adsorbed gas proportion. When the average reservoir pressure has dropped to about 50% of the initial formation pressure at the eighth year of production, the cumulative adsorption gas contribution rate reaches 15%, and the EUR recovery rate is 75%. In the 10th year, the reservoir pressure drops to 21.8 MPa, and the cumulative adsorption gas contribution rate is close to 18%. The utilization of the reserves can reach 80% of EUR. The final average formation pressure drop is about 67%, the calculated cumulative gas production can reach nearly 137.9 million cubic meters, and then, the contribution rate of adsorbed gas to EUR is 23.28%. Meanwhile, the annual adsorbed gas contribution rate curve in Figure 6 shows that the annual contribution rate of adsorbed gas production gradually increases with the production time. Specifically, the annual adsorbed gas production accounts for less than 10% of the gas output in the first year, and the annual adsorbed gas rate exceeds 40% in the 12th year. The adsorbed gas plays a more important role in the later gas production composition than free gas in the 20th year, and the final annual contribution rate of adsorbed gas can be 53.71%. It can be seen that the desorbed gas plays a key part in production composition mainly in the later stage of the production period, which is beneficial to long-term stable production of gas wells.Figure 4
Shale gas EUR recovery rate and average formation pressure versus production time.Figure 5
Adsorbed gas and free gas contribution rate versus production time.Figure 6
The adsorbed gas annual contribution rate versus the production time curve.It can be seen from Figure7 that when the average formation pressure is greater than 37 MPa, the cumulative gas production has a linear relationship with the average formation pressure, and the linear slope is about 0.050.9 billion cubic meters/MPa. During this period, the cumulative gas production volume curve and the free gas production curve are basically overlapped, which means the free gas is mainly produced from the shale gas reservoir. While the formation pressure drops to 37 MPa, the cumulative gas production curve starts to deviate from the straight section and gradually curves upward; then, the adsorbed gas begins to be utilized. The low average formation pressure makes the adsorption curve and the cumulative gas production curve easily turn up, whose slopes are, respectively, 0.062 million cubic meters/MPa and 0.0261 million cubic meters/MPa. It can explicitly explain that the adsorbed gas is a key part of later production composition in the shale reservoir.Figure 7
The relationship curve between cumulative gas production and average formation pressure.


## 4. Shale Gas Production Influence Factors Analysis

The gas production decline curve shows that the actual production supply of the gas well mainly comes from the gas of the fracture in the stimulated zone and the mass transfer in the matrix near the fracture system (Figure8). After large-scale hydraulic fracturing of the shale reservoir, the induced secondary fracture network and natural fractures in the near-well area intersect each other, forming high-conductivity seepage channels. The shale gas in the fractures under the motivation of pressure difference will quickly flow into the wellbore from the fracture, resulting in a rapid increase in daily gas production in the early production stage. When the production time is about 38 days, the daily gas production reaches a peak of 325,900 m3/day. After the shale gas in the fracture flows into the wellbore, the gas well enters the gas supply in the matrix control stage. Caused by the low porosity and the low permeability of the shale reservoir, the production rate of shale gas in the matrix to the fractures is fairly slow, which makes the gas production from the matrix cannot supply the gas production from fracture system with high conductivity in time so that the gas production rate drops rapidly until the occurrence of shale gas desorption. Then, the obvious non-Darcy effects such as gas diffusion and slippage under low pressure will enhance the gas production to a certain extent. It is expected that the gas well will enter a long-lasting stable production period on the 800th day. Therefore, the initial production of the shale gas well is mainly affected by the gas flow in the fracture, while the shale gas production during the later production stage is obviously provided by the matrix long-term development. Based on the production prediction model in the paper, clarifying the influence factors for the high production of gas wells in the gas in the fracture control stage and gas supply in the matrix is helpful to propose theoretical optimization suggestions and guide measurements to enhance the reserve production rate and the final gas recovery level.Figure 8
The daily gas production and production time of the gas well curve.


### 4.1. Fracturing Complexity Index

Fracture Complexity Index (FCI) [24] reflects the complexity of hydraulic fractures and the degree of effective communication between fractures. It is defined as the ratio of the width (Xn) to the length (2Xf) of the hydraulic fracture network monitored by microseismic data. Then, the influence of different fracture network indexes on the gas well production is studied based on the model in the paper.The above Figure9 shows that when other factors remain unchanged, the fracture complexity index mainly affects the production performance of the gas well in the early and middle stages. The larger the fracture complexity index means the smoother daily production decline curve of the gas well, and the cumulative gas production of gas wells will increase significantly. When the FCI increased from 0.25 to 0.4, the EUR of gas wells in 20 years increased by 3.53%, and the output of free gas and adsorbed gas increased by 1.9% and 6.9% in 20 years. It can be seen that the higher FCI means more complex fractures formed by hydraulic fracturing and better reservoir fracturing effect, which greatly promote the production of free gas and the desorption of adsorbed gas in the fracture system.Figure 9
Effect of different FCI value on shale gas production.
(a)
The curve of gas production verse time under different Fracture Complexity Index(b)
The curve of adsorbed gas, free gas cumulative production verse time under different Fracture Complexity Index


### 4.2. Permeability of Secondary Fracture Network

In the actual gas production process, the secondary fractures play an important role in the gas mass transfer between the matrix and the hydraulic fractures. The effect of the different secondary fracture permeability on the gas production, free gas, and adsorption has been studied.From Figure10, it can be intuitively reflected that when the secondary fracture network permeability increases from 0.1 mD to 5 mD, the cumulative gas production curve of gas wells increases significantly, and the daily gas production decline rate decreases. The EUR, free gas volume, and adsorption gas volume of the 20-year gas well are increased by 11.71%, 6.56%, and 22.42%. This is because the high permeability of the secondary fracture network can denote the high conductivity of the secondary fractures to lower the gas flow resistance in the fracture network and fasten the spreading speed of the pressure drop funnel section. Finally, the gas flow area is widened so that the free gas in the fracture network is easily developed. Due to the continuous increase in the secondary fracture network permeability, the extremely low permeability of the matrix makes it far from sufficient to supply the gas in the fracture, resulting in a relatively slow growth rate of daily gas production and cumulative gas production.Figure 10
Effect of secondary fracture network permeability on shale gas production.
(a)
The curve of gas production verse time under different secondary fracture network permeability(b)
The curve of adsorbed gas, free gas cumulative production verse time under different secondary fracture network permeability


### 4.3. Stress Sensitivity of Secondary Fractures

When the effective formation pressure increases, secondary fractures tend to show obvious stress sensitivity. When other parameters remain unchanged, the influence of pressure sensitivity coefficients on the productivity of gas wells is studied.Figure11 shows that the daily gas production and cumulative gas production have a strong negative correlation with the stress sensitivity coefficient. The higher stress sensitivity coefficient can result in the smaller daily gas production and cumulative gas production of gas wells. The EUR, cumulative adsorbed gas volume, and cumulative free gas volume considering the max stress sensitivity are lower than that without the stress sensitivity by 31.33%, 66.58%, and 17.91%. It can be seen that the gas production reduction in the fracture without proppant leads to obvious stress shadow areas near the fracture, which cannot achieve the effective gas transmission between the fracture and the wellbore, resulting in a significant declination in gas production. Therefore, ceramsite proppants with higher strength and good pressure resistance should be selected during hydraulic fracturing to reduce the loss of fracture network permeability.Figure 11
Effect of stress sensitivity of secondary fractures on shale gas production.
(a)
The curve of gas production verse time under different secondary fracture stress sensitivity(b)
The curve of adsorbed gas, free gas cumulative production verse time under different secondary fracture stress sensitivity


### 4.4. Bottom Hole Pressure

The selection of reasonable bottom hole pressure has a significant impact on the production of gas wells, which can help delay the rapid attenuation of reservoir pore pressure, maintain the long-term opening of fracture systems, and alleviate the stress-sensitivity of the seepage field so as to increase the final cumulative gas production of a single well. Based on this model, the bottom hole pressure is set as 0.1 MPa, 5 MPa, 10 MPa, and 15 MPa, and the gas production characteristics of adsorbed gas and free gas are studied.As is shown in Figure12, the bottom hole pressure has a significant impact on the production in the middle and late stages of gas well. The daily gas production and cumulative gas production of gas wells increase significantly with the decrease of bottom hole pressure. Lowering bottom hole pressure increases the pressure gradient between the reservoir and bottom hole, effectively promoting gas production in fractures and early utilization of the gas in the matrix. The EUR, adsorbed gas volume, and free gas volume of the gas well at 5 MPa increased by 6.38%, 12.53%, and 3.54% compared to that at 10 MPa. However, excessively low bottom hole pressure promotes large effective stress, leading to the closure of secondary fractures, and the deformation of reservoir flow channels, so that the gas in the fracture systems cannot be effectively produced. Therefore, during the process of formulating a reasonable bottom hole pressure, it is necessary to consider the impact of the effective stress on the fracture conductivity to control the closure of fractures, alleviate the production rate reduction, and ensure the full utilization of the gas wells production, thereby increasing the accumulation of gas wells production and shale reservoir recovery rate.Figure 12
Effect of bottom hole pressure on shale gas production.
(a)
The curve of gas production verse time under different bottom hole pressure(b)
The curve of adsorbed gas, free gas cumulative production verse time under different bottom hole pressureTo sum up, the production decline curve of the shale gas well has the typical characteristics of high initial production, rapid decline rate, and long-term low production period in the later stage, forming a typical L-shaped production decline curve. Enhancing the shale gas reservoir production rate and obtaining the optimal productivity must rely on the gas supply in the fracture network control area and the matrix control area. First, three-dimensional hydraulic fracturing of the shale gas reservoirs should be implemented to improve the initial gas production. Reservoir rocks are transformed into granular form, and multilevel secondary fractures are induced under the effective stress difference between the hydraulic fractures and the reservoir to strengthen the communication between secondary fractures and natural fractures, maximize the fracturing control area in the reservoir, and increase the effective permeability of the reservoir, which can greatly increase the early productivity of shale gas wells. In addition, the gas supply in the matrix is the key guarantee for achieving a long-term stable production of gas wells in the later stage. The extremely low matrix permeability can result in several key issues that are not conducive to efficient production such as difficult pressure propagation and rapid decline of gas production rate in the matrix of the shale reservoir. Therefore, valid bottom hole pressure can be carried out to both relieve stress sensitivity and increase the effective pressure gradient between the fracture network and the matrix to promote the production of free gas and the desorption of adsorbed gas in the matrix to ensure the earlier occurrence of the gas well stable production period, which is beneficial to enhance the reservoir effective utilization rate and final recovery rate.


## 4.1. Fracturing Complexity Index

Fracture Complexity Index (FCI) [24] reflects the complexity of hydraulic fractures and the degree of effective communication between fractures. It is defined as the ratio of the width (Xn) to the length (2Xf) of the hydraulic fracture network monitored by microseismic data. Then, the influence of different fracture network indexes on the gas well production is studied based on the model in the paper.The above Figure9 shows that when other factors remain unchanged, the fracture complexity index mainly affects the production performance of the gas well in the early and middle stages. The larger the fracture complexity index means the smoother daily production decline curve of the gas well, and the cumulative gas production of gas wells will increase significantly. When the FCI increased from 0.25 to 0.4, the EUR of gas wells in 20 years increased by 3.53%, and the output of free gas and adsorbed gas increased by 1.9% and 6.9% in 20 years. It can be seen that the higher FCI means more complex fractures formed by hydraulic fracturing and better reservoir fracturing effect, which greatly promote the production of free gas and the desorption of adsorbed gas in the fracture system.Figure 9
Effect of different FCI value on shale gas production.
(a)
The curve of gas production verse time under different Fracture Complexity Index(b)
The curve of adsorbed gas, free gas cumulative production verse time under different Fracture Complexity Index


## 4.2. Permeability of Secondary Fracture Network

In the actual gas production process, the secondary fractures play an important role in the gas mass transfer between the matrix and the hydraulic fractures. The effect of the different secondary fracture permeability on the gas production, free gas, and adsorption has been studied.From Figure10, it can be intuitively reflected that when the secondary fracture network permeability increases from 0.1 mD to 5 mD, the cumulative gas production curve of gas wells increases significantly, and the daily gas production decline rate decreases. The EUR, free gas volume, and adsorption gas volume of the 20-year gas well are increased by 11.71%, 6.56%, and 22.42%. This is because the high permeability of the secondary fracture network can denote the high conductivity of the secondary fractures to lower the gas flow resistance in the fracture network and fasten the spreading speed of the pressure drop funnel section. Finally, the gas flow area is widened so that the free gas in the fracture network is easily developed. Due to the continuous increase in the secondary fracture network permeability, the extremely low permeability of the matrix makes it far from sufficient to supply the gas in the fracture, resulting in a relatively slow growth rate of daily gas production and cumulative gas production.Figure 10
Effect of secondary fracture network permeability on shale gas production.
(a)
The curve of gas production verse time under different secondary fracture network permeability(b)
The curve of adsorbed gas, free gas cumulative production verse time under different secondary fracture network permeability


## 4.3. Stress Sensitivity of Secondary Fractures

When the effective formation pressure increases, secondary fractures tend to show obvious stress sensitivity. When other parameters remain unchanged, the influence of pressure sensitivity coefficients on the productivity of gas wells is studied.Figure11 shows that the daily gas production and cumulative gas production have a strong negative correlation with the stress sensitivity coefficient. The higher stress sensitivity coefficient can result in the smaller daily gas production and cumulative gas production of gas wells. The EUR, cumulative adsorbed gas volume, and cumulative free gas volume considering the max stress sensitivity are lower than that without the stress sensitivity by 31.33%, 66.58%, and 17.91%. It can be seen that the gas production reduction in the fracture without proppant leads to obvious stress shadow areas near the fracture, which cannot achieve the effective gas transmission between the fracture and the wellbore, resulting in a significant declination in gas production. Therefore, ceramsite proppants with higher strength and good pressure resistance should be selected during hydraulic fracturing to reduce the loss of fracture network permeability.Figure 11
Effect of stress sensitivity of secondary fractures on shale gas production.
(a)
The curve of gas production verse time under different secondary fracture stress sensitivity(b)
The curve of adsorbed gas, free gas cumulative production verse time under different secondary fracture stress sensitivity


## 4.4. Bottom Hole Pressure

The selection of reasonable bottom hole pressure has a significant impact on the production of gas wells, which can help delay the rapid attenuation of reservoir pore pressure, maintain the long-term opening of fracture systems, and alleviate the stress-sensitivity of the seepage field so as to increase the final cumulative gas production of a single well. Based on this model, the bottom hole pressure is set as 0.1 MPa, 5 MPa, 10 MPa, and 15 MPa, and the gas production characteristics of adsorbed gas and free gas are studied.As is shown in Figure12, the bottom hole pressure has a significant impact on the production in the middle and late stages of gas well. The daily gas production and cumulative gas production of gas wells increase significantly with the decrease of bottom hole pressure. Lowering bottom hole pressure increases the pressure gradient between the reservoir and bottom hole, effectively promoting gas production in fractures and early utilization of the gas in the matrix. The EUR, adsorbed gas volume, and free gas volume of the gas well at 5 MPa increased by 6.38%, 12.53%, and 3.54% compared to that at 10 MPa. However, excessively low bottom hole pressure promotes large effective stress, leading to the closure of secondary fractures, and the deformation of reservoir flow channels, so that the gas in the fracture systems cannot be effectively produced. Therefore, during the process of formulating a reasonable bottom hole pressure, it is necessary to consider the impact of the effective stress on the fracture conductivity to control the closure of fractures, alleviate the production rate reduction, and ensure the full utilization of the gas wells production, thereby increasing the accumulation of gas wells production and shale reservoir recovery rate.Figure 12
Effect of bottom hole pressure on shale gas production.
(a)
The curve of gas production verse time under different bottom hole pressure(b)
The curve of adsorbed gas, free gas cumulative production verse time under different bottom hole pressureTo sum up, the production decline curve of the shale gas well has the typical characteristics of high initial production, rapid decline rate, and long-term low production period in the later stage, forming a typical L-shaped production decline curve. Enhancing the shale gas reservoir production rate and obtaining the optimal productivity must rely on the gas supply in the fracture network control area and the matrix control area. First, three-dimensional hydraulic fracturing of the shale gas reservoirs should be implemented to improve the initial gas production. Reservoir rocks are transformed into granular form, and multilevel secondary fractures are induced under the effective stress difference between the hydraulic fractures and the reservoir to strengthen the communication between secondary fractures and natural fractures, maximize the fracturing control area in the reservoir, and increase the effective permeability of the reservoir, which can greatly increase the early productivity of shale gas wells. In addition, the gas supply in the matrix is the key guarantee for achieving a long-term stable production of gas wells in the later stage. The extremely low matrix permeability can result in several key issues that are not conducive to efficient production such as difficult pressure propagation and rapid decline of gas production rate in the matrix of the shale reservoir. Therefore, valid bottom hole pressure can be carried out to both relieve stress sensitivity and increase the effective pressure gradient between the fracture network and the matrix to promote the production of free gas and the desorption of adsorbed gas in the matrix to ensure the earlier occurrence of the gas well stable production period, which is beneficial to enhance the reservoir effective utilization rate and final recovery rate.


## 5. Summary and Conclusions

(1)
Based on the conventional five-zone composite flow model, the supercritical adsorption, diffusion, slippage stress sensitivity, and the gas high-pressure physical property are comprehensively considered in the improved flow model, which makes up for the lack of consideration of gas nonlinear effects in previous productivity models(2)
The proposed production model was verified by the field production data and illustrated that the EUR predicted by the productivity model is higher than the EUR that the comprehensive nonlinear effects are not considered, which emphasized that the gas nonlinear effect in the production prediction process cannot be ignored(3)
Based on the production prediction model, the 20-year production performance of the gas well is studied: in the first year of gas well production, the adsorbed gas is basically not recovered, and the cumulative adsorption contribution rate does not exceed 10%. The final adsorption gas contribution rate is 23.28%, and the annual adsorption rate can exceed 50% in the 20th year(4)
The cumulative gas production curve and the free gas production curve are basically overlapped when formation pressure is higher than the critical adsorption pressure; then, the gas cumulative gas production rises significantly when the formation pressure drops below the critical desorption pressure. Thus, free gas mainly affects the productivity in the early and midstage of production, and adsorbed gas is an important source of long-term stable production when the reservoir is in low pressure(5)
The production influence factors of the fracture network stage and the matrix stage are studied, including fracture network permeability, stress sensitivity, fracture complexity index, and bottom hole pressure. The optimization suggestions for improving the shale reservoir production are proposed: (1) in the process of hydraulic fracturing, the wide-range of three-dimensional fracturing reformation of shale reservoirs should be implemented to increase the effective early production of the reservoir; (2) reasonable bottom hole pressure can be used in the later matrix development process to increase the effective pressure gradient between the fracture network and the matrix to promote the production of free gas and adsorbed gas in the matrix



---

*Source: 9920122-2021-09-24.xml*