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

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.


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][3][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][10][11][12][13][14][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.

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: The specific assumptions are as follows: (1) The migration of single-phase methane gas in different seepage fields is considered as the onedimensional 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 as y e The gas in the matrix of zone 4 flows into the matrix of zone 2 along the y 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: where the matrix comprehensive compression factor is Gas compression factor is 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): Javadpour apparent permeability model [19] is used to describe the gas diffusion and slip in the matrix: where slip speed correction factor is The hydraulic flow radius is 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 The pseudotime is According to the definition of pseudopressure and pseudotime and the definition of the dimensionless parameters in Table 1, Eq. (1) can be rewritten as follows: The initial matrix pressure of zone 4 is equal to the original formation pressure: Outer boundary : Inner boundary : The dimensionless forms of Eq. ((11)-(13)) are as follows: Then, the pseudopressure in the dimensionless flow equations can be transformed as equation (15) in Laplace form [20]: Next, both sides of the dimensionless flow control equation (10) are multiplied by e −stD simultaneously and integrated by the dimensionless time from zero to infinity, which turns out as the equation (16):

Geofluids
The integral term at the left side of equation (16) is simplified as equation (17): The integral term at the right side of equation (16) is simplified as equation (18): 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: Boundary conditions can be expressed in Laplace form as equations ((20), (21)): The dimensionless pseudopressure solution of zone 4 in Laplace space is expressed in equation (22) based on equation (19) and boundary conditions ((20), (21)):

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 the y direction, the dimensionless gas percolation equation in the matrix of zone 3 can be derived in equation (23): 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:  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): Dimensionless pressure 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 Outer boundary is Inner boundary is The dimensionless pseudopressure solution of the matrix in zone 2 in Laplace space can be deduced in equation (32) based on Laplace transformation: where

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): The dimensionless gas flow equation in fracture network zone 1 can be depicted in expression (35): The boundary conditions are as followed in equations ((36)-(38)): Initial condition is Outer boundary is Inner boundary is 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 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: The term 1/1 − γ 1D ζ 1D in equation (40) can be expanded into Taylor form in equation (41): 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): The initial and boundary conditions are simplified as follows: Initial conditions are Outer boundary is Inner boundary is Similarly, combining the Laplace transform, the pseudopressure solution of equation (42) can be expressed in equation (46): where

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: 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 Outer boundary is Inner boundary is 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): where 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): 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.

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 Table 2. 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).
It can be seen from Figure 2 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 m 3 /day to 10,000 m 3 /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 6 Geofluids accuracy so as to reasonably evaluate the medium and longterm productivity of gas wells and predict the EUR value.

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 Figures 4 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. It can be seen from Figure 7 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   7 Geofluids 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.

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 ( Figure 8). 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 highconductivity 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 m 3 /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     Free gas proportion Adsorbed gas proportion Figure 5: Adsorbed gas and free gas contribution rate versus production time. 8 Geofluids 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.

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 Figure 9 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   Figure 6: The adsorbed gas annual contribution rate versus the production time curve.  Declining period Stable production period Figure 8: The daily gas production and production time of the gas well curve. 9 Geofluids 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.

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 Figure 10, 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,  10 Geofluids 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.

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. Figure 11 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.

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 pro-duction 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 Figure 12, 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.
To 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 Adsorbed gas-without stress sensitivity Adsorbed gas-stress sensitivity coefficient: 0.14 MPa-1 Free gas-without stress sensitivity Free gas-stress sensitivity coefficient: 0.14 MPa-1 (b) The curve of adsorbed gas, free gas cumulative production verse time under different secondary fracture stress sensitivity Figure 11: Effect of stress sensitivity of secondary fractures on shale gas production.
11 Geofluids 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.

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 20year 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 20 th 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 12 Geofluids the matrix to promote the production of free gas and adsorbed gas in the matrix Nomenclature A cw : Wellbore crossflow area, m 2 C di : Modified supercritical desorption gas compression coefficient, Pa -1 C gi : The gas compression coefficient, Pa -1 C ti : The matrix comprehensive compressibility coefficient, Pa -1 C fi : The reservoir compressibility coefficient, Pa -1 P sc : The standard atmospheric pressure, Pa T sc : The standard temperature, K T: The reservoir temperature, K Z sc : The ideal gas Compression factor, dimensionless Z: The gas compression factor, dimensionless Bgi: Original volume coefficient, dimensionless Bg: The gas volume coefficient at reservoir condition, dimensionless VL: The Langmuir volume, m 3 /m 3 PL: The Langmuir pressure, Pa P e : The seepage field pressure, Pa F: The slip velocity correction factor, dimensionless μðiÞ: The gas viscosity of the different seepage field, Pa.s (Subscript i = 2, 3, 4 represents the three matrix regions, respectively) P avg : The average gas pressure, Pa D: The diffusion coefficient, m 2 /s Kia: The apparent matrix area permeability, m 2 Ki: Darcy permeability, m 2 ri: The circular capillary radius, m C tF : The hydraulic fracture comprehensive compressibility coefficient, Pa -1 x e : Effective horizontal length, m ye: Well spacing, m w: The width of a hydraulic fracture, m h: The reservoir thickness, m P: Formation pressure, Pa Pi: Original formation pressure/Pa ta: The pseudo time, s y F : The fracture half length, m K: The reservoir permeability, mD N: The number of hydraulic fractures, m K 1 : The fracture network permeability after considering the stress sensitivity effect, m 2 K 1i : The initial fracture network permeability, m 2 K F : The hydraulic fracture permeability, m 2 .

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.