#### Abstract

Shale gas reservoir has been aggressively exploited around the world, which has complex pore structure with multiple transport mechanisms according to the reservoir characteristics. In this paper, a new comprehensive mathematical model is established to analyze the production performance of multiple fractured horizontal well (MFHW) in box-shaped shale gas reservoir considering multiscaled flow mechanisms (ad/desorption and Fick diffusion). In the model, the adsorbed gas is assumed not directly diffused into the natural macrofractures but into the macropores of matrix first and then flows into the natural fractures. The ad/desorption phenomenon of shale gas on the matrix particles is described by a combination of the Langmuir’s isothermal adsorption equation, continuity equation, gas state equation, and the motion equation in matrix system. On the basis of the Green’s function theory, the point source solution is derived under the assumption that gas flow from macropores into natural fractures follows transient interporosity and absorbed gas diffused into macropores from nanopores follows unsteady-state diffusion. The production rate expression of a MFHW producing at constant bottomhole pressure is obtained by using Duhamel’s principle. Moreover, the curves of well production rate and cumulative production *vs.* time are plotted by Stehfest numerical inversion algorithm and also the effects of influential factors on well production performance are analyzed. The results derived in this paper have significance to the guidance of shale gas reservoir development.

#### 1. Introduction

Shale gas is an unconventional natural gas with the main component of methane existing as adsorbed or free gas in the rich organic shale or its interlayers, the reserve of which is tremendous around the world. With the innovative techniques of horizontal well drilling as well as hydraulic fracturing in North American or other places of the world, the annual production of shale gas in the USA has increased from 564 × 10^{8} m^{3} in 2007 to 4382 × 10^{8} m^{3} in 2015. Moreover, well testing and production performance analysis in shale gas reservoir, especially for multiple fractured horizontal well (MFHW), is a hot research topic in recent years [1]. Compared to conventional gas reservoirs, shale gas reservoirs have special features on the accumulation conditions, occurrence manner, percolation mechanism, etc. For instance, shale is a self-sourcing reservoir within which the gas is stored either by compression in the pores and natural fracture as free gas or by adsorption on the surface of the solid material (either organic material or minerals) as adsorbed gas. According to Hill and Nelson [2], the adsorbed gas accounts for up to 85% of total gas, which is dependent on a variety of geologic and geochemical properties [3]. The pore network of shale gas reservoir is complex, which can be further divided into 4 categories: nanoinorganic pore, nanoorganic intragranular pore, macronatural fracture, and large-scaled hydraulic fractures. The first two pore types belong to matrix pores and the latter two belong to fractures [4–7], so the mechanisms of gas flow in shale reservoirs is very complicated. Shale gas reservoirs are generally described by the dual porosity model like Warren-Root model [8], in which the matrix blocks serve as the main storage space of gas and are separated by fractures.

Kucuk and Sawyer [9] first analyzed the transient pressure response of wells in shale gas reservoirs by using the dual porosity model (assuming cylindrical and spherical matrix geometries), and the corresponding analytical and numerical results were obtained. However, the analytical model presented in the paper did not take the effects of gas desorption and diffusion into account and the numerical model ignored the effects of gas diffusion.

On the basis of isothermal adsorption equation, Bumb and McKee [10] analyzed the transient pressure of a vertical well in coal bed methane (CBM) reservoir by introducing an additional compressibility coefficient in their model. Their method has been popularly used for production performance analysis of unconventional gas. Many researchers followed Bumb and McKee’s [10] approach to consider the ad/desorption effects in their studies; the results of which were proved to be much more realistic than previous studies [11–14].

Ertekin et al. [15] utilized a dual-mechanism model for tight gas reservoirs, which considered Darcy flow (driven by pressure gradient) and Fick diffusion (driven by concentration difference). Subsequently, Carlson and Mercer [16] investigated the gas flow behavior in shale gas reservoir by coupling the gas diffusion and desorption with Fick’s law and Langmuir isothermal adsorption equation. And then, using the same method, Chawathe et al. [17] analyzed the pressure and saturation distribution in the reservoir by numerical solution.

Recently, many studies related to well pressure and rate transient analysis have been presented [5, 18–33]. El-Banbi [34] utilized a linear flow model to analyze the production performance of a MFHW by ignoring the gas desorption and diffusion. In this model, the reservoir was simplified by continuous medium with slab dual porosity and Warren-Root model [19, 20, 35]. However, most of these studies did not consider desorption or diffusion effects, such as “trilinear” model [24] and other models [36, 37].

In fact, a large number of core electron macroscope scanning images show that there are many organic pores existing in the shale matrix, the scale of which ranges from nanometer to macrometer. Previous studies considered the matrix pores as a whole without separating nanometer pores and macropores. Song [38, 39] proposed a “triple porosity” conceptual model to describe the gas storage mechanism, which was equivalent to a dual porosity model combined with gas adsorption. Thereafter, Zhao et al. [40] adopted the model to analyze the pressure transient and rate decline behavior of MFHW in shale gas reservoirs. Recently, a comprehensive mathematical model of a MFHW in a box-shaped reservoir is presented by Zhao et al. [32]. In this model, the reservoir was divided into a dual porosity system, which included natural fractures and nanopores. The gas flow in natural fracture system obeyed the Darcy flow theory, while the Knudsen diffusion and slippage effect were considered for the matrix low permeability system.

According to Dehghanpour and Shirdel [41], the dual porosity models have been widely used in naturally fractured shale gas reservoirs with multiple hydraulic fractures. In these models, only nanopores existed in the matrix blocks which were surrounded by a fracture network. Gas flow in the nanopores obeys the Fickian diffusion law and the storage capability of the fracture system is very low, the theoretical productivity of the well from such conceptual models is quite low, and the decline rate of well productivity is very large, which is not consistent with some practical data ([5, 24, 27]).

Most physical models proposed by previous investigators did not take into account the nanopores and the macropores in the matrix system [32, 40, 42]. According to the NMR scanning results of the core from YY32 well from the YanChang formation in P.R. China (Figure 1), two peaks are the most important characteristics for many shale gas reservoirs, which means that there are mainly two types of pores in the matrix system: the nanopores and macropores.

In this paper, a new conceptual model is developed to describe the gas flow in shale gas reservoir. The flow path is presented by into three pore types: natural fractures, macropores, and nanopores. It is assumed that the adsorbed gas first desorbs from the organic particle surface into the nanopore system. It then diffuses into the macropore system. Thereafter, the gas in the macropores is transported into the natural fractures with transient interporosity flow. Based on this physical model, the production performance of a MFHW in rectangular outer-boundary shale gas reservoirs is analyzed.

#### 2. Multiple Transport Mechanisms of Shale Gas Reservoirs

##### 2.1. Gas Ad/Desorption Mechanism

A large amount of shale gas is adsorbed on the surface of rock particles. As the reservoir pressure depletes, adsorbed gas desorbs from the matrix surface and becomes free gas. This feature differentiates shale gas reservoirs different from conventional reservoirs. The isothermal adsorption curves have the same shape as Langmuir isothermal curve [10], which is frequently used to describe ad/desorption process of shale gas.

The occurrence time of gas desorption behavior also depends on the amount of adsorbed gas at initial conditions. As shown in Figure 2, if the adsorbed gas is statured under initial conditions, the amount of adsorbed gas is at point B and the gas quickly desorbs from organic surfaces as pressure decreases. However, if the adsorbed gas is unsaturated, the amount of adsorbed gas would be at point C (if reservoir conditions are at point C and there is no free gas in matrix pores) and the gas would not desorb from organic surfaces until reservoir pressure decreases to the pressure at point A, which is called critical desorption pressure. The difference between the critical desorption pressure and the initial reservoir pressure determines the starting time of gas desorption [1].

Langmuir isothermal adsorption equation is expressed as
where is the gas volumetric concentration, m^{3}/m^{3}; is the Langmuir adsorption volume of shale gas, representing the limit adsorbence of unit reservoir, m^{3}/m^{3}; is the Langmuir pressure, which is the pressure when gas adsorbence reaches 50% of limit adsorbence, Pa; and is the pressure in matrix system, Pa.

The mass flow rate of desorbed gas from shale reservoir of volume in unit time is
where is the mass of desorption gas in unit time, kg/s; is the time, s; is the shale matrix volume, m^{3}; and is the shale gas density at standard conditions, kg/m^{3}.

If gas desorption happens instantaneously and all desorbed gas transforms into free gas, then (2) can be introduced directly into the continuity equation of the matrix or fracture system, resulting in the steady-state adsorption-desorption model ([10–13]; Civan, 2010; [37]; Wang et al., 2013).

##### 2.2. Gas Flow in Natural Fractures and Macropore System

Gas flow in natural and hydraulic fractures and macropores is continuous flow. According to previous analysis (Civan, 2010; [37]), Darcy’s flow equation used for conventional gas reservoir can be adopted to analyze gas flow behavior in fracture and macropore systems in shale reservoirs, that is,
where presents the natural fracture system and the presents the macropore system; is the gas deviation factor (dimensionless), *Z* = 1 for ideal gas; is the gas constant (8.314 J/(mol ⋅ K)); is the reservoir temperature, K; is the permeability in system, m^{2}; is the gas molar mass, kg/mol; is the pressure in system, Pa.

##### 2.3. Gas Diffusion in Nanopores

Gas flow behavior in shale is different from flow in fractures because the pores are of nanoscale. In these nanopores, it is erroneous to assume molecular continuous Darcy flow. Shale gas flow in nanopores includes not only viscous flow but also diffusion. In some ultratight shale reservoirs, the flow of gas molecules occur by only diffusion because the gas in the nanopores exists only as absorbed gas and not free gas.

Fick’s law of diffusion is generally used to describe shale gas migration in nanopores and at the surface of kerogen [15, 40, 43–46]. Under the influence of concentration gradient, desorbed gas migrates from high concentration area to low concentration area, and the diffusion stops when gas concentration becomes uniform (as illustrated in Figure 3). Assume that gas molecular diffusion satisfies Fick’s law, and then gas diffusion flux could be expressed as
where is the Fick diffusion coefficient, m^{2}/s.

If gas concentration in matrix changed with coordinates, meaning nonequivalent concentration at any time in the matrix, unsteady state diffusion occurs, which can be described by Fick’s second law. Otherwise, the Fick’s first law is used when the concentration is assumed to be the same at all locations in a given time (pseudosteady-state diffusion). It is assumed in this paper that the desorption of the adsorbed gas desorbed from the nanopores into the macropores follows unsteady-state diffusion [15].

To simplify the hypothetical shale gas reservoir model, it is assumed that the matrix element is spherical with radius ; each matrix element is composed of organics and clay grains, and diffusion is the only gas flow mechanism in the matrix (Figure 4). Therefore, Fick’s second law can be applied to describe gas diffusion from matrix elements into macrofractures or macropores.

According to the unsteady-state diffusion theory, gas concentration in the matrix is a function of time and space. Therefore, the following mathematical expression can be used to describe gas concentration change in the matrix,

The initial condition of matrix unsteady-state diffusion should satisfy the equation below. where is the initial reservoir pressure, Pa.

The inner boundary condition can be expressed as

Since the outer boundary of matrix is connected with the fracture system or macropore system, pressure at the matrix surface equals that of fracture system, then, the outer boundary can be represented as where is the radius of a spherical matrix element, m.

Using the above model and relevant boundary conditions, the distribution of gas concentration in the shale matrix can be solved. Then, the following equation can be used to determine gas diffusion flow rate within matrix of volume ,

#### 3. Solutions of a MFHW in a Rectangular Shale Gas Reservoir

##### 3.1. Physical Model

In this section, continuous point source solution for any point in a rectangular shale gas reservoir with closed upper and lower boundaries is derived. The physical model is shown in Figure 5. The assumptions for this model are (1) the gas reservoir is homogeneous and anisotropic with horizontal and vertical permeability and . The width, length, and height of the reservoir are , , and ; (2) a MFHW with length is located in the middle of the reservoir and parallel to the upper and lower boundaries. All hydraulic fractures are symmetrical and perpendicular to the wellbore; the half length of the th fracture is ; the fractures are fully penetrating; (3) the well is perforated at the intersections between the fractures and the wellbore, which means that fluid flows from the reservoir to the fractures and then from the fractures to the wellbore. Infinite conductivity fractures are assumed with no pressure drop across the fracture face and along the horizontal wellbore length; (4) each fracture is evenly divided into 2*N* units with unified flow rate [32, 47].

##### 3.2. Continuous Point Source Solution

During the development of shale gas reservoirs, natural fractures are the main paths for gas flowing from the matrix to hydraulic fractures and wellbores (Figure 6). Therefore, flow models for these natural fractures are the key to deriving the point source function. The PDE of pseudopressure in cylindrical coordinates of fracture systems for different models are derived in Laplace domain (detailed derivation in the Appendix), which is where is the parameter group of the corresponding model under different percolation mechanism, dimensionless.

For unsteady-state diffusion and transient interporosity flow, the expression of is where can be expressed as

According to the mirror image and superposition principles (the schematic is shown in Figure 7) and Newman [48] as well as Ozkan and Raghavan [49], the instantaneous point source solution at the observation point generated by unit source sink at point is [1]

In the above equation, represents the pressure response at any point of the gas reservoir caused by each image well due to closed boundary. where

Introducing the following Poisson summation formula, which is [32, 49] Then, the continuous point source solution for rectangular gas reservoirs can be simplified and expressed as [32, 49] where and

Then, the pressure drop in rectangular reservoir generated by the intensity of the continuous point source is

Introducing (17) into (20) yields

The above expression is the pressure drop at any point of the reservoir generated by constant intensity in rectangular reservoir. Noted that it can clearly see from the above equation that some infinite summations of eigenfunctions are existed, which generally have slow convergence characteristics. Since the solutions presented in this work are in the Laplace transform domain, we utilized the same approach proposed by Ozkan and Raghavan [49] to improve the computation efficient of source function given in (21). The main routine to overcome this problem is to separate in the infinite-acting and boundary-domained flow contributions in the solution and then recast the slow converging summations and computationally difficult integrals into computationally more efficient forms. In the previous paper [32], the detailed expressions during different time are given, and many authors has also reported them [49, 50].

##### 3.3. Solution of the MFHW Producing with Constant Bottomhole Pressure

As is shown in Figure 5, a horizontal well with *M* fully penetrated fractures is drilled in a rectangular outer-boundary reservoir. According to the source function theory, the pressure drop at arbitrary point caused by each discrete element can be obtained by integrating the continuous point source function along the fracture element, which is [40, 51]

Due to the fractures are parallel to the *x-*axis and all of them are fully penetrated, the above equation can be converted into the following formulation:

Taking (21) into (23) and arranging it yield where is the continuous line source strength, which has the following relationship with the continuous point source :

Defining the following dimensionless production rate and pseudopressure, we have and

Introducing (26) and (27) into (24), the dimensionless pseudopressure drop generated by the ^{th} unit at any reservoir location can be expressed as

Due to the horizontal well is composed of many discrete fracture units, the pressure at any point in the reservoir should be the superposition of the pressure drop caused by all these units at the same point. For the observation point , according to the superposition principle, there is [32, 40]

If setting the observation point on each fracture element, then we can obtain linear equations and another auxiliary equation of mass conservation for flux can be also obtained; the detailed expressions can be referred to the paper of Zhao et al. [32].

When well producing at a constant bottomhole pressure, the dimensionless rate can be defined as follows:

According to previous studies [52], there is the following relation formula between dimensionless pseudopressure and rate, where the definition of dimensionless pseudopressure in (31) is

The above equation is the production rate solution when well producing at a constant pressure.

#### 4. Results Analysis

##### 4.1. Model Validation

In this section, we use the solution of conventional dual-porosity reservoirs to prove the correctness of our model. Comparing to the conventional model, the difference between them is that our mode considers the gas desorption and diffusion from the nanopores into macropores. The model presented in this paper can simplify into the conventional dual-porosity model by neglecting the shale gas desorption and diffusion processes. For conventional reservoirs, approximation expressions for in (11) can be found in the literatures [8, 9, 53–55]. The expression of introduced by Warren and Root for blocked matrix with pseudosteady-state interporosity flow is

Using de Swaan’s approach [53], the solutions for Warren and Root’s model can be extended to the solutions of transient interporosity flow, which are for horizontal slab blocks and for spherical blocks.

The above three types of matrix block are the most popular idealization model for naturally fractured reservoirs; they have been widely used in the petroleum industry since they are proposed. So, the correctness of our model and method can be validated by comparing the result to the conventional model with spherical blocks.

In our model, if the gas diffusion and desorption are neglected, the following relationship must be established, which is

Substituting the above equation into the expression of in (11), which becomes

It can be clearly seen that (11) can be simplified into the same solution for spherical blocks described by (35) for naturally fractured reservoirs by taking the above equation into it, which prove the correctness of our solution.

##### 4.2. Well Production Performance Analysis

On the basis of the above solution, the semianalytical solution in Laplace domain can be transformed into real space using Stehfest numerical inversion algorithm [56]. Therefore, when the well is produced at a constant bottomhole pressure, the rate decline semilog type curves of production rate and cumulative production *vs* time can be obtained. The input parameters for simulation sensitivities for the following models are shown in Table 1.

For a well producing at constant bottomhole pressure in a rectangular gas reservoir, the production rate and cumulative production vs. time curves in semilog coordinates with different fracture numbers are shown in Figure 8. From the plot, it can be seen that the early time production rate is increasing with the increasing of fracture number, which is mainly because the larger the *M* is, the bigger the drainage area of the well will be during the early flow period. And when the pressure wave transports to the outer region formation, this effect will be gradually weaken. However, the cumulative production with different *M* under the same production time in later production period will be the same due to the reserve of them are equivalent.

Figure 9 shows the influence of Langmuir volume on the well production performance. Due to the mechanism, physical model in this paper has divided the pores in matrix system into the following two subspace: macropores and nanopores. The adsorbed gas will be desorbed from the particle surface in the matrix and then flow into the macropores by Fick diffusion. When the pressure in the macropore system is lower than the Langmuir pressure, desorption phenomenon will happen. So, when the well producing with the constant bottomhole pressure of 4 MPa, the pressure in the vicinity of the well will be lower than the Langmuir pressure, which will cause the adsorbed gas to desorb from matrix surface. So, we can find from Figure 9 that the production rate is larger with the increasing of Langmuir volume. But this effect is very small. When well is producing for a long time, the cumulative production of well with bigger will be larger than the others.

Figure 10 shows the effects of permeability in natural fracture system on the well production performance. Because the natural fracture system is the main flow path of gas flow from the formation into hydraulic fractures, it has a significant influence on well production performance. According to the productivity equation of wells, the well production rate has a proportional relationship with the natural fracture permeability. And due to the constant control volume of our model, the bigger production rate at early flow period will lead to a bigger rate decline in the middle flow period. According to the cumulative production curves on Figure 10, we can clearly see that the larger the is, the larger the cumulative production will be in our given time period.

Figure 11 shows the effects of permeability in macropores on well production rate and cumulative production rate curves. The bigger represents a stronger supplement of fluid from matrix system into natural fractures. Because the fluid flow from matrix system into natural fractures follows the transient interporosity flow, this supplement will happen immediately after the well opens. The bigger the is, the larger the well production rate will be, and this trend will be lasting for a long time. So, for shale gas reservoir development, the permeability of the matrix system is a very important parameter to evaluate the well production performance.

Figure 12 shows the effects of macropores porosity () on well production performance. The value of presents the amount of free gas in macropore system of shale matrix system. With the increase of the , the quantity of free gas flows from macropores into natural fracture system by interporosity will increase. So, when well produces at a constant bottomhole pressure, the higher the porosity is, the larger the well production rate will be. Although the increment of production rate is not quite obvious, for the situation listed in this paper, the increment is still nearly 5000 m^{3}/d. If this increment for conventional gas reservoir could be negligible, it is very considerable for shale gas reservoir due to the long production period.

Figure 13 shows the effects of gas desorption time () on the well production performance. As we know, the gas desorption time presents the gas desorption rate of adsorbed gas desorbed into free gas. The smaller the is, the faster the desorption process happens. It can be clearly seen from Figure 12 that the smaller the gas desorption time is, the larger the well production time will be, and the corresponding cumulative production will be larger too.

#### 5. Discussion

This paper established a new mechanism model to describe the gas flow in multiscaled shale gas reservoir, which divided the pores in matrix systems into two subtypes: macropores and nanopores. The gas flow in nanopores obeys the Fick’s second diffusion and in macropores follows the interporosity flow as conventional dual-porosity model. On the basis of the above flow mechanism, the production performance of a MFHW in a rectangular outer-boundary shale gas reservoir is analyzed by semianalytical solution method. Comparing to the solutions of conventional naturally fractured reservoir, the correctness of our model is certified by neglecting the gas desorption and diffusion properties. After that, the effects of corresponding sensitive parameters on well production rate and cumulative production rate curves are analyzed. According to the results, we can conclude that the fracture number only affects the well performance in early flow period; the permeability of natural fractures mainly affects the early performance due to the finite drainage volume in our physical model. Comparing to the permeability of natural fracture system, the effects of permeability in macropore system on early and middle flow periods are obvious. So, for some high rate well in field, the macropore may be an important storage space, and which cannot be neglected or combined with other system to analyze.

#### Appendix

#### Continuity Equation in Natural Fracture System

The continuity equation of shale gas flow in macrofracture system is given in the cylindrical coordinate system by [32]

For transient interporosity flow from matrix to fractures, the continuity equation for flow in matrix is

The interporosity flow rate is expressed as

Gas diffusion from matrix to fractures in unit reservoir volume is

For unsteady-state gas diffusion from matrix to fractures, there is

For transient interporosity flow model, define the following dimensionless variables:

The other variable groups used in the above equations are defined as where the expression of pseudopressure is defined as where is the dimensionless interporosity flow coefficient from nanopore into macropore system; is the dimensionless interporosity flow coefficient from macropore into natural fracture system; is the sorption time constant which can be measured in the lab; is the storability ratio, dimensionless; and is the characteristic length, m; in this paper, we can choose the horizontal well length as the reference length.

And then the continuity equation in the natural fracture and macropore system in the matrix system can be prescribed as the following in Laplace domain, which are where .

According to the mathematical model of gas diffusion in the nanopores described in (5)–(8) and combining the defined dimensionless variables with them, we can obtained the dimensionless diffusion equation in the Laplace domain is [1]

Define the following group as

Substituting the above equation into (A.11) yields

The general solution is

Substituting (A.16) into (A.14) yields

According to the inner boundary condition of matrix system, we have

According to the L’Hospital rule, the solution of parameter can be derived from (A.18)

Then, the following equation can be obtained from the inner boundary condition:

Similarly, the above equation can be simplified as

Combining with the outer boundary condition and taking into the corresponding equation, the parameter can be obtained as

Substituting (A.19) and (A.22) into (A.17), and then the solution of the gas concentration is

According to the Langmuir isotherm desorption equation, the following equation is established:

Define the following desorption coefficient as follows:

Then, (A.24) becomes

Substituting (A.26) into (A.23) yields

Taking the derivation with for the above equation, we have

Then, substituting (A.28) into (A.10), we have where .

Then, based on the previous research [8], the relationship of macropores in matrix system with pressure in natural fracture system can be solved. So, the uniform expression of the governing equation in natural system can be changed into where .

#### Nomenclature

*Latin*

c:_{mgi,} c_{fgi} | Gas compressibility in fracture and matrix system, respectively, Pa^{−1} |

: | Molar mass of gas in matrix, mol/m^{3} |

: | Volume of gas adsorbed per unit volume of the shale gas reservoir in equilibrium at pressure pm, sm^{3}/m^{3} |

: | Langmuir volume, sm^{3}/m^{3} |

: | Diffusion coefficient, m^{2}/s |

: | The parameters group of the corresponding model under different percolation mechanism, dimensionless |

f, mc: | Fractures and macropore media |

: | Formation thickness, m |

: | Fick mass diffusion flux (gas mass passing through unit acreage in unit time), kg/(m^{2} ⋅ s) |

: | Mass velocity of gas in media , kg/(m^{2} ⋅ s) |

: | Permeability of natural fracture system, m^{2} |

: | Permeability of macropore in matrix system, m^{2} |

: | Horizontal and vertical permeability of reservoir, respectively, m^{2} |

: | Well length, m |

: | The characteristic length, m |

: | Number of fractures |

: | Discrete number of each half fracture |

: | Molar mass of gas, kg/mol |

: | Pore pressure, Pa |

: | Langmuir pressure, Pa |

: | Initial reservoir pore pressure, Pa |

: | Pressure in fracture and matrix system, respectively, Pa |

: | Pressure at standard condition, |

: | Pressure in media , Pa |

: | Interporosity flow rate, kg/s |

: | The intensity of the cylinder source sink |

: | Mass flow rate of gas desorption from reservoir of volume , kg/s |

: | Well production rate, m^{3}/s |

: | Fick mass flow(gas mass passing through volume in the unit time), kg/s |

: | Inner diameter of sphere matrix element, m |

: | Radius of sphere matrix, m |

: | Laplace variables, dimensionless |

: | Production time, hour, and day |

: | Standard temperature, |

: | Reservoir temperature, K |

: | Shale matrix volume, m^{3} |

: | The half length of the th fracture, m |

: | Distance coordinates, m |

: | Effective reservoir width and length, m |

: | Dimensionless space coordinates, dimensionless |

: | Dimensionless space distance between the well with the bottom boundary, dimensionless |

: | Gas deviation factor, fraction. |

*Greek*

𝛌: | Interporosity flow coefficient, dimensionless |

: | Sorption time, constant, s |

: | Storability ratio, dimensionless |

Gas viscosity, Pa ⋅ s | |

: | Shale gas density at standard conditions, kg/m^{3} |

: | Gradient operator, |

: | Natural fracture and shale matrix porosity, respectively, fraction |

: | Desorption coefficient |

: | 3.141….. |

*Superscript*

−: | Laplace domain. |

*Subscripts*

D: | Dimensionless |

I: | Initial condition |

sc: | Standard state |

w: | Wellbore |

f: | Fracture system |

m: | Matrix system |

wD: | Dimensionless wellbore. |

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was supported by the National Basic Research Program of China (Grant No. 2014CB239205) and the National Natural Science Foundation of China (Grant Nos. 51704247, 51504202).