#### Abstract

Shale reservoirs have the characterizations of low porosity, low permeability, and hydrocarbon organic matter self-generation and self-storage, resulting in its complex flow mechanisms. Compared with fractured vertical wells, multiple-fractured horizontal wells are widely used due to their advantages of effectively increasing the well control range and further expanding the drainage area. To further study the multiscale flow mechanisms of shale gas, a flow model was established that considered viscous flow in microfractures and inorganic pores, the diffusion of Knudsen in nanoscale porosity, the coexistence of slippage, adsorption-desorption effects under infinity, and closed outer boundary conditions; based on the continuous point source solution, a multiple-fractured horizontal well flow model was established and solved by MATLAB programming. Then, the effects of various factors were investigated. The results show that the Knudsen diffusion and slippage coefficients mainly affect the apparent permeability of the matrix pores. The more the Knudsen diffusion and slippage coefficients are, the earlier the turbulent flow occurs and the higher the gas production is.

#### 1. Introduction

The revolution of shale gas in the United States of America (USA) is changing the energy structure of the world [1, 2]. Based on its various forms in the reservoir, the shale gas can be divided into free shale gas, adsorbed shale gas, and dissolved shale gas. Most of the shale gas is deposited in pores and microfractures in a free state, and the proportion of the gas adsorbed on the surface of organic matter and clay mineral particles [3] could generally reach more than 50% [4]. The dissolved gas exists in liquid hydrocarbons or adsorbed on the surface of other substances in the kerogen [5–7] as shown in Figure 1. Through tank degassing experiments, Cheng et al. [8] observed the diffusion of gas from the kerogen or clay particles to the surface of the pores.

Scholars have committed to the research on shale gas flow models since a decade ago. In 2009, a mathematical model was developed for horizontal wells with triple media branches [9], and the bottom well pressure solution of the horizontal well with the triple media branches in Laplace space was solved by using the point source function theory. Soon, a pressure calculation formula for any point in a stratum with fractures at different inclinations in a fractured horizontal well was developed based on the same theory [10]. Then, Ozkan et al. [11] applied the trilinear flow model to fractured horizontal wells, which is suitable for shale reservoirs with a matrix permeability of only micron-scale Darcy and nanoscale Darcy. This model assumes that the artificial vertical fractures are formed after fracturing, and the surrounding of the fractures is considered to be a dual medium.

In the same year, Stalgorova and Mattar [12] proposed a trilinear flow model to analyze the production dynamics of conventional reservoirs and fractured horizontal wells in unconventional gas reservoirs. In this model, the entire reservoir was divided into three regions (the external reservoir, the internal reservoir, and the hydraulic fracture). The flow of gas between the three regions was assumed to be linear in this model. The results indicated that the most effective way to increase the productivity of unconventional gas reservoirs is to increase the fracture density. The effects of fracture permeability on productivity are not obvious. Then, a new trilinear flow model [13] was proposed based on Ozkan’s model. Similarly, the model divided the entire reservoir into three parts: the fractured zone, the unreformed zone, and the artificial fracture. The results studied by this model points out that due to the multiple factors such as rock properties and ground stress, the hydraulic fractures formed during the fracturing process are not necessarily perpendicular to the horizontal wellbore but form a certain fracture network structure. In 2014, Dejam et al. [14] established a fractured horizontal well flow model for shale gas reservoirs based on the triple-porosity four-linear flow model, which considers the effects of gas adsorption and desorption effect of gas in the matrix. The model mainly studies the quasisteady state turbulence of the matrix to the microfracture medium and the large-fracture medium. In 2017, the unsteady flow regimes of a slightly compressible fluid under the linear and radial pre-Darcy flow conditions are modeled and the corresponding highly nonlinear diffusivity equations are solved analytically by aid of a generalized Boltzmann transformation technique [15]. The influence of pre-Darcy flow on the pressure diffusion for homogeneous porous media is studied in terms of the nonlinear exponent and the threshold pressure gradient. In addition, the pressure gradient, flux, and cumulative production per unit area are compared with the classical solution of the diffusivity equation based on the Darcy flow. In 2018, a two-phase flowback model is developed with multiscale diffusion mechanisms, and it is found that the two-phase flowback and the flow consistency between the matrix and fracture network have significant influences on cumulative gas production. The multiscale diffusion mechanisms in different zones should be carefully considered in the flowback model [16]. In 2019, a model considering convective flow, gas diffusion, and surface diffusion was established to investigate gas transport mechanism, and a common practice in modeling shale gas permeability is to use Knudsen diffusion coefficient when calculating diffusive flux, but the use of Knudsen diffusion coefficient would be incorrect if the shale gas flow regime is lying in either the transition diffusion or Fick’s diffusion, in which case the diffusion coefficient must correspond to that regime [17]. Then, a series of shale gas adsorption and desorption experiments are conducted [18]. Desorption and adsorption curves are not coincident, with the former located above the latter, which suggests that adsorption hysteresis also occurs in shale gas.

A multiple-fractured horizontal well (MFHW) is extensively applied to develop the shale gas reservoirs, which makes the production of unconventional reservoirs economically and practically feasible [19]. After massive fracturing, the complex fracture network will exist around the horizontal wellbore. To model the flow behavior more rigorously, the fractures within the network should be represented explicitly rather than idealized as dual-porosity media around the wellbore. A mathematical model considering the stress sensitivity of the reservoir permeability, Darcy flow, diffusion, and adsorption and desorption in shale gas reservoirs is developed, and the numerical nonlinear production decline equations are derived and obtained. Subsequently, the model is verified by a simplified model, and production decline curves for a MFHW with discrete fracture networks in shale gas reservoirs are plotted [20]. An analytical solution is developed for the shale gas productivity of a multiple-fractured horizontal well based on a diffusion model and a trilinear flow pattern. The shale gas reservoir is divided into three flow regions: hydraulic-fracture region, microfracture network or dual-porosity region, and pure-matrix region. For the pure-matrix region, a transient diffusion equation is solved based on our previous diffusivity model developed for the shale matrix. For the microfracture network region, a modified dual-porosity model is proposed wherein both the free and adsorbed gases in the shale matrix flow into the microfracture network through a pseudosteady diffusion process. A dimensionless solution is obtained for the bottom-hole pressure in the Laplace domain considering the skin effect. An analytical solution is obtained for the gas production rate in a real-time domain through a partial Taylor series simplification and Laplace inverse transform. This analytical solution is compared with the field data of the shale gas produced from a fractured horizontal well located in southwestern China, and a good agreement is observed [21]. A semianalytical model of finite-conductivity multiple-fractured horizontal wells in LC gas reservoirs is established based on the Laplace-space superposition principle and fracture discrete method. The proposed model is validated against a commercial numerical simulator. Type curves are obtained to study pressure characteristics and identify flow regimes. The effects of some parameters on type curves are discussed [22].

Different models are featured in different aspects, so there are many problems at present. (1) The migration mechanisms of shale gas in the pores are not clear. In addition, the equation describing the flow of shale gas in the pores is mostly derived under the conditions of single capillary and ideal gas. (2) The effects of flow law are not considered comprehensively in the existing models, and the models do not fit well with the actual conditions; it is difficult for the model to fully describe the Darcy flow, Knudsen diffusion, adsorption-desorption, Fick diffusion, and other mechanisms. (3) The research results and methods of shale gas in microcosmic are difficult to be used in macroscopic research. The research results between microcosmic flow and macroscopic flow are not closely related.

First, the shale gas multiscale flow mechanisms are presented. Then, a flow model was established that considered viscous flow in microfractures and inorganic pores, the diffusion of Knudsen in nanoscale porosity, the coexistence of slippage, and adsorption-desorption effects under infinity and closed outer boundary conditions; based on the continuous point source solution, a multiple-fractured horizontal well flow model was established and solved by MATLAB programming, and the effects of various factors were investigated.

#### 2. Shale Gas Multiscale Flow Mechanisms

Unlike conventional reservoirs, shale gas reservoirs have the characteristics of low porosity, low permeability, self-generation, and self-storage of hydrocarbon organics. The spatial structure of the reservoir is a typical multiscale feature. The migration of shale gas in the reservoir is a typical multiscale flow usually including the following three stages: (1) free gas flows from pores to the wellbore; (2) after the pore pressure drops, the adsorbed gas on the pore wall surface begins to desorb into the pores and then become free gas; and (3) dissolved gas in organic matter diffuses to the pore surface [23–25].

##### 2.1. Flow in Fractures and Macropores

The flow of shale gas in natural fractures and large matrix pores satisfies Darcy’s law, which is is the gas flow velocity in m/s, is the medium permeability in , is the gas viscosity in mPa·s, and is the gradient operator, .

During the development of shale gas reservoirs, natural fractures and matrix systems exhibit stress-sensitive effects as formation pressures decrease. The permeability is a function of formation pressure at this time, which is is the medium system pressure in Pa.

##### 2.2. Flow in Nanoscale Pores

In shale reservoirs, nanoscale pores are mainly organic pores, of which the velocity of gas molecules is not zero at the wall surface, and both the slip flow and the Knudsen effect exist in pores, as shown in Figure 2.

There would be slippage when the fluid flowed in porous media [26]. Therefore, the following formula was proposed to calculate the apparent gas permeability of reservoir rocks: where is the apparent permeability in , is the equivalent liquid permeability in , is the slip factor which depends on the gas type and pore throat characteristics in Pa, and is the average pressure of the flow path in Pa.

There are many models of visual permeability in shale nanoscale pores, which are summarized in Table 1, and the slip coefficient expressions of various scholars are obtained. It is noteworthy that these models only consider the slippage mechanism. The Knudsen diffusion terms must be added if these models are applied to nanoscale pores.

In this paper, the expression of the coefficient of permeability correction [27] considering slippage effect and Knudsen is as follows:

is the sparse function related to ; it can be expressed as follows:

When the flow pattern of shale gas in the nanoscale pores is slippage, , , , and .

##### 2.3. Gas Desorption on Organic Surface

A large amount of shale gas is absorbed on the surface of organic matter under actual formation conditions, and desorption occurs mainly during mining [28–32]. There are many factors affecting the desorption of shale gas, such as organic matter content, organic matter maturity, temperature, and pressure [33–35]. The relationship between various factors is complicated.

In the establishment of the flow model, the instantaneous desorption model is used to characterize the shale gas. The expression of the Langmuir isotherm equation is as follows:

The mass flow of desorption per unit volume of shale gas reservoir per unit of time during extraction is
where is the shale gas mass flow desorbed from a unit volume of reservoir rock in kg/s and is the shale gas density under standard conditions in kg/m^{3}.

Substituting Equation (6) into Equation (7),

If the Langmuir equation is expressed in pseudopressure form, the above equation can be changed to

At this point, Equations (8) and (9) can be used to express the instantaneous amount of desorption gas, and then the instability flow model of shale gas wells can be studied.

##### 2.4. Diffusion in Porous Kerogen (Organic Matter)

There are two forms of diffusion: steady state diffusion and nonsteady state diffusion. When the gas consistence in the porous kerogen is independent of the spatial position and only related to time, it is called pseudosteady state diffusion and can be described by Fick’s first diffusion law; when the gas concentration in organic matter is related to both time and space coordinates, it is called nonsteady state diffusion and can be described by Fick’s second diffusion law.

###### 2.4.1. Pseudosteady State Diffusion

The derivative of shale gas consistence in the matrix with respect to time is proportional to the consistence difference of gas inside and outside the organic matter, which can be expressed as

Therefore, the diffusion gas flow from the unit volume of shale organic matter to the nanoscale pores is
where is the molar concentration of shale gas in organic matter under pseudosteady state conditions in mol/m^{3}; is the molar concentration of gas at organic wall in mol/m^{3}; is the nanoscale pore pressure; is the shape factor in 1/m^{2}, determined by the shape of the organic matter; is the Fick diffusion coefficient in m^{2}/s; is the time in s; is the Fick mass flow in kg/(m^{3}·s); and is the geometric factor; it is determined by the organic shapes that are listed in Table 2.

###### 2.4.2. Nonsteady State Diffusion

The diffusion of dissolved gas in organic matter is usually an unsteady state diffusion. If there is a stable and continuous diffusion source, the pseudosteady diffusion rule will be met. Based on theory of uniformity, this paper assumes that the organic matter is spheres of equal size with radius of . The physical model is shown in Figure 3.

According to the unsteady diffusion theory, the dissolved gas consistence is a function of time and spatial position. For the spherical matrix model, it is assumed that the gas consistence on the outer surface of the organic matter is in dynamic equilibrium with the free gas in the nanoscale pores. The gas consistence change in the organic matter can be described by the following mathematical model:
where is the volumetric consistence of shale gas in organic matter under unsteady conditions in mol/m^{3} and is the radial coordinates of spherical organic rock mass in m.

Before a gas reservoir is put into production, the pressure at each point in the reservoir is the original formation pressure , and the gas concentration at each point in the organic matter is related to its internal pressure. Furthermore, it is assumed that the rate of change in gas concentration at the center of the organic matter is zero, and the gas concentration at the outer boundary of the organic matter is the same as the gas concentration in the nanoscale pores. Therefore, the initial boundary conditions for the nonsteady state diffusion in the organic matter satisfy the following process: where is the reservoir original formation pressure in MPa, is the nanoscale pore pressure connected to the outer surface of organic matter in MPa, and is the spherical organic radius in m.

According to the above model and the corresponding definite conditions, the dissolved gas concentration distribution in the organic rock mass can be obtained, and then the mass diffusion of gas per unit volume of organic matter can be determined according to

is the Fick mass flow in kg/(m^{3}·s).

#### 3. Establishment and Point Source Solutions for Flow Model

##### 3.1. Characterization of the Microscopic Flow Model

Strictly speaking, the shale reservoirs should be divided into four parts of organic matter that existed which are microfracture systems, inorganic pores, nanoscale pores, and gas dissolved organic matter [36, 37]. A comprehensive flow model is proposed in this section that considers multiple porous media, including the viscous flow in microfractures and inorganic pores, the diffusion of Knudsen in nanoscale porosity, the coexistence of slippage, adsorption-desorption effects, and the diffusion of dissolved gases in organic matter [38, 39].

###### 3.1.1. Dissolved Gas Diffusion in Porous Kerogen

The diffusion of dissolved gas in porous kerogen can be described by the nonsteady state Fick diffusion law:

###### 3.1.2. Flow in Nanoscale Pores in Kerogen

There are three sources of gas in organic pores, namely, free gas in pores, desorption gas on the surface of organic matter, and dissolved gas diffused from kerogen. The flow equation in the pores of spherical kerogen is as follows:
where is the radial radius of kerogen in nanoscale pores in m, is the gas density in kerogen nanoscale pores in kg/m^{3}, is the gas permeation velocity in kerogen nanoscale pores in m/s, is the ratio of kerogen volume to rock volume, is the inherent permeability of kerogen nanoscale pores in mD, is the gas mass rate desorbed from organic pore in kg/s, is the mass rate of gas diffused from organic matter in kg/s, is the Langmuir volume in kerogen in kg/s, and is the Langmuir pressure in kerogen in MPa.

The correction coefficients of permeability can be taken as follows: , , , and .

###### 3.1.3. Flow in Matrix Macropore

The flow equation of matrix macropore in spherical coordinates can be written as

###### 3.1.4. Flow in Microfracture System

The mass conservation equation of the microfracture system is

##### 3.2. The Solution of Microscopic Flow Model

The general equation of flow for the corresponding microfracture system under various flow models is as follows: where is the corresponding parameter and is the Laplace variable.

The developed microscopic flow model divides shale reservoirs into four parts: the microfracture system, inorganic pores, nanoscale porosity pores, and organic matter with dissolved gas. A flow model that is a comprehensive consideration of porous media was established, including viscous flow in microfractures and inorganic pores, Knudsen diffusion in nanoscale porosity, coexistence of slippage, adsorption-desorption effects, and dissolved gas diffusion in organic matter. The solution processes of the various microscopic flow mechanism model are in the appendix.

The dimensionless variables are defined in Table 3.

The resulting expression is as follows:

So far, we have gotten the expressions of dimensionless pseudopressure of microfractures under various flow mechanisms of shale gas reservoirs and given their unified form in spherical coordinates.

##### 3.3. Continuous Point Source Solution of Circular Boundary Shale Gas Reservoir

###### 3.3.1. Flow Model

In this section, the continuous point source solution at any point in an anisotropic circular gas reservoir with a closed upper and lower boundary is derived. The physical model is shown in Figure 4. The main assumptions in the model are as follows: (1) the gas reservoir is homogeneous, the thickness is , the initial pressure is ; (2) the permeability of the gas reservoir in the horizontal and vertical directions are and , respectively; (3) the middle distance of the gas reservoir is the lower boundary , where exists a cylindrical source sink. The radius of the microelement cylinder is and the height is . (4) The source and sink strength of the cylinder is , and (5) the gas flow in the gas reservoir meets the Darcy flow, and the effects of capillary pressure and gravity on the flow are ignored.

According to the flow theory of oil and gas reservoirs, the continuity equation of gas flow is
where is the density of gas in crack systems in kg/m^{3}, is the radial flow velocity of gas in m/s, is the flow velocity of gas in the vertical direction in m/s, is the porosity of the fracture system, is the radial and vertical coordinates, and is the production time in s.

The flow equations of gas flow in the transverse and longitudinal directions are as follows:
where is the permeability of the reservoir fracture system in the horizontal direction in m^{2}, is the permeability of the reservoir fracture system in the vertical direction in m^{2}, is the crack system pressure in Pa, and is the gas viscosity in Pa·s.

According to the state equation of real gas, the density of gas can be expressed as follows: where is the gas molar mass in kg/mol, is the absolute temperature of the gas in K, and is the gas compression factor.

According to the definition of rock compressibility, the following relationship exists between the porosity of the rock fracture system and the pressure of the fracture system: where is the fracture pore compression factor.

For the inner boundary conditions, it is assumed that gas enters mainly from the side of the cylinder, and as the radius of the microelement tends to 0, it can be considered as a continuous point sink, and its inner boundary expression is [38]

Among them, the volume coefficient of gas is as follows:

The upper and lower boundaries of the gas reservoir are all closed. Therefore, the upper and lower boundary conditions can be described as

The outer boundary conditions can be described as

Introducing the following proposed pressure,

Substituting the gas density, volumetric coefficient, and pseudopressure definition into the flow equation of the fracture system,

For the above continuity equation, the viscosity and compressibility of the gas are all functions of pressure. In order to obtain the analytical solution of the equation, the viscosity and compressibility of the part are usually regarded as the values in the original state.

Several dimensionless variables are defined in Table 4.

Substituting the above dimensionless quantities into Equations (30) and (31), the flow equation of the fracture system becomes as follows:

The inner boundary condition becomes as follows:

The upper and lower boundary conditions become as follows:

The outer boundary condition becomes as follows:

###### 3.3.2. The Solution of the Flow Model

Since there are many derivation formulas for solving similar flow models (appendix), this section does not carry out long derivation but directly gives the general solution of the model; as shown in Equation (38), the solution process requires Laplace transform and Fourier cosine transform. where , , and .

###### 3.3.3. The Continuous Point Source Solution under Different Outer Boundary Conditions

Above, we have solved the hybrid outer boundary model. However, the outer boundary of a gas reservoir has the following two situations: infinite and closed boundaries. Therefore, the continuous point source solutions were analyzed in these two cases.

*(1) Infinity Outer Boundary*. When the outer boundary of a gas reservoir is infinite, the outer boundary equation is as follows:

Therefore,

Substituting Equation (40) into Equation (39), we obtain a continuous point source solution with closed top and bottom boundaries and infinite outer boundary:

The source intensity for a continuous point is a constant. Thus, the above equation can be written as follows:

*(2) Closed Outer Boundary*. When the outer boundary is closed, it is known from the closed boundary conditions:

When the outer boundary is closed, the continuous point source solution is as follows:

#### 4. Study on Unsteady Flow Rules of Circular Boundary Multiple-Fractured Horizontal Wells

##### 4.1. The Establishment of a Physical Model

A model diagram of a multiple-fractured horizontal well in a circular shale gas reservoir is shown in Figure 5. The following assumptions are given: (1)The gas reservoir is homogeneous, horizontal, and equal in thickness. The upper and lower boundaries are closed. The lateral direction is infinite or closed and the reservoir thickness is (2)The horizontal well is located in the center of the gas reservoir and its length is (3)The pressure drop caused by fluid flowing in the wellbore is ignored, and the hydraulic fractures supply fluid to the wellbore, regardless of the flow of the reservoir fluid directly to the wellbore(4)The number of hydraulic fractures is , and it is a fully cracked fracture with infinite flow conductivity. The fracture is regarded as a 2D surface source rather than a 3D volume source(5)Fracture cracks are symmetrical bifurcations distributed at equal intervals along the horizontal wellbore(6)The flow process is an isothermal flow, regardless of the effects of fluid gravity and capillary forces

##### 4.2. Mathematical Model Establishment and Solution

As shown in Figure 5, the direction of the wellbore along the horizontal well is the -axis, the vertical upward direction is the -axis direction, and the direction along the fracture surface is the -axis direction. For multistage fracturing horizontal wells, the flow process of fluids in the reservoir is quite complicated. Not only is the mutual interference between the seams and the fractures, but the flow on the fracture surface is also uneven. Therefore, the solution method for continuing to use the fracturing vertical wells does not work. We can first divide the crack into several units, and it is assumed that the flow on each section is evenly distributed, but the flow in the cracks on different units is different. By solving the pressure response expressions at each fracture cell and combining the bottom hole pressure with the normalized production conditions, the expressions of the bottom hole pressures in multistage fracturing horizontal wells in shale gas reservoirs are solved.

In this study, the fractures are infinitely inducing flow cracks that are completely fracturing. According to the continuous point source solution of the circular reservoir obtained, it is integrated along the longitudinal direction that we can obtain a continuous linear source solution. For a continuous line source with an intensity of , there is an infinite atmospheric boundary: where and .

The continuous line source intensity and continuous point source intensity have the following relationship:

The definition of the dimensionless distance is as follows:

Substituting Equations (46) and (47) into Equation (45), we can obtain a continuous linear source solution for an infinite boundary shale gas reservoir with a circular outer boundary:

In the same way, the solution of a continuous line source can be obtained when the gas reservoir boundary is closed:

The predecessors have basically formed a set of theoretical systems for the crack dispersion and superposition principle to solve the pressure response. The following is a brief introduction. In the hypothetical part of the model in this section, the number of fractures is , the length of horizontal wells is , and each fracture is equally divided into units as shown in Figure 6.

From the discrete fracture schematic, the node coordinates of each segment in the fracture can be obtained. The crack node coordinates in the negative direction of the -axis can be obtained by the following equation:

The coordinates of the crack node in the positive direction of the -axis can be obtained by the following equation: where is the abscissa of the -th crack cell in m, is the ordinate of the -th crack cell, and is the section -th crack length.

For the crack unit length, since each crack is divided into parts, the length of the -th crack cell is

As mentioned earlier, it has been assumed that the flow distribution of each fracture cell is uniform, so that the pseudopressure drop formed at any point () in the formation at the -th crack cell is

Since the cracks are all parallel to the -axis, the above curve integrals can be directly converted into coordinate integrals. When the outer boundary of the reservoir is infinite, Equation (49) is substituted into Equation (53), where the coordinates of the -th node are substituted into, we can get the following equation:

Since the flow is uniformly distributed on the discrete units of the fracture, assuming that the -th discrete unit is , the following relationship can be satisfied between the line flow and the unit flow:

Introducing the following dimensionless pressure and dimensionless production:

Substituting Equations (56) and (57) into Equation (54), where is the dimensionless length of the -th crack unit.

According to the superposition principle of the potential, the total pressure response generated at each position of the crack in the cracks at any position in the reservoir is

According to the previous assumption, where is the intersection point of the -th section of a fracture and the wellbore.

When a horizontal well is produced at a given production volume, the production of the well is numerically equal to the sum of the flow of all cracks. Therefore, we have the following equation:

Substituting Equation (61) into Equation. (57), there are the following normalization conditions:

In Figure 6, we separate the cracks into units and solve the pseudopressure at each unit, so we can have equations, which, plus the normalized equation for production, total equations, while the unknowns include the flow rate and bottom hole pressure of each fracture cell, which is also . The equations are closed and there is a unique solution. For ease of understanding, this set of equations is written in the following matrix form:

The expression of is as follows:

The above formula is the fracturing horizontal well production and pressure calculation formula under the condition of the infinite outer boundary. If the boundary of the gas reservoir is closed, the formula can be simplified to Equation (64):

##### 4.3. Gas Well Bottom Pressure Dynamic Curve Analysis

So far, we have obtained the unstable pressure of multiple-fractured horizontal wells in circular reservoirs considering a variety of complex flow mechanisms. According to the expression, the computer programming and Stehfest numerical inversion techniques [38] were used to obtain pressure and production dynamic curves of shale gas reservoir fractured.

Figure 7 shows the typical pressure response curves of multistage fracturing horizontal wells in shale reservoirs with different outer boundary conditions. Multiple concaves appear in the dynamic curve, reflecting the effects of various flow mechanisms on the pressure response. As the pressure wave propagates in the formation, different pore media rich in shale gas begin to flow, causing the pressure of the porous media to fall.

According to Figure 7, we can divide the flow of horizontal wells into the following stages:

*Stage I*: the wellbore storage and subsequent transitional flow stages. In the wellbore storage phase, the pressure and pressure derivative curves coincide with each other with a slope of 1

*Stage II*: the early linear flow stage perpendicular to the fracturing fracture. After the completion of the well storage phase, the free gas in the reservoir flowed first in the direction perpendicular to the fracture surface, and there was no mutual interference between the fractures. In this flow stage, the slope of the pressure derivative curve is 1/2

*Stage III*: early radial flow stage. The pressure derivative curve shows a horizontal line with a slope of

*Stage IV*: the turbulence stage. This stage corresponds to the linear flow stage from the formation to the entire fracturing zone. The concave character appears on the pressure derivative curve

*Stage V*: the boundary control stage. The characteristics of the boundary flow are exhibited after the pressure wave spreads to the boundary. For the closed boundary, when the pressure wave spreads to the boundary, the pressure will increase faster; combined with the continuous line source solution Equation (49) under the closed boundary, the pressure and pressure derivative curves coincide and the slope is 1

The parameters used in the mechanistic model are listed in Table 5.

The pressure curves for different slippage coefficients under double logarithmic coordinates are shown in Figure 8. Making the slippage coefficient equals 1, 8, and 16, it can be seen that as the slippage coefficient increases and the turbulence trough appears earlier and earlier. According to the knowledge of flow mechanics, when the shale gas desorbed from the organic pore, slippage appears, so the slippage coefficient affects the turbulence stage the most. The greater the slippage coefficient is, the more obvious is the slippage effect, and the larger the apparent permeability of the matrix is, so that the turbulence of matrix pores to the fracture system is more serious, and eventually, the time of turbulence appears earlier. Naturally, a larger slippage coefficient will increase the production of shale gas wells.

Figure 9 shows the influence of different Knudsen diffusion coefficients on the pressure response of multiple-fractured horizontal wells. It is not difficult to see from the figure that the influence of Knudsen diffusion coefficient (m^{2}/s) on the pressure dynamics and the gas slippage effect on the well test curve of gas wells are similar. With the increase of (m^{2}/s), the turbulence troughs appear earlier and the pressure and pressure derivative curves are getting lower and lower, so the production of gas wells will also increase.

#### 5. Summary and Conclusions

(1)A flow model for shale gas reservoirs considering flow in microfractures and inorganic pores, the diffusion of Knudsen in nanoscale porosity, the coexistence of slippage, established adsorption-desorption effects, and the continuous point source solution of the flow mechanisms model were obtained under closed and infinity outer boundaries(2)Based on the continuous point source solution, a circular boundary multiple-fractured horizontal well flow model was established, and the gas well bottom pressure dynamic curve and the relevant sensitivity factors are analyzed with the help of MATLAB programming(3)The results of the sensitivity analysis show that the Knudsen diffusion and slippage coefficients mainly affect the apparent permeability of the matrix pores. The more the Knudsen diffusion and slippage coefficients are, the earlier the turbulent flow occurs and the higher the gas production is

#### Appendix

#### A. The Diffusion of Dissolved Gas in Porous Kerogen

The diffusion of dissolved gas in porous kerogen can be described by the nonsteady state Fick diffusion law:

The internal boundary conditions:

The outer boundary conditions:

The initial conditions:

The following dimensionless variables were defined as follows:

Firstly, the aforementioned equations become dimensionless, then the Laplace transformation was conducted:

Ordering , we can get

#### B. Flow in Nanoscale Pores in Kerogen

The flow equation in the pores of spherical kerogen is as follows:

According to the research of Javadpour [23] where

The internal boundary conditions of the porous kerogen are

The outer boundary conditions:

The initial conditions:

The following dimensionless variables were defined:

Firstly, the aforementioned equations become dimensionless, and then the Laplace transformation was conducted:

The general form of the equation:

#### C. Flow in Matrix Macropore

The macropore flow control equation of the matrix under the spherical coordinates can be written as follows:

The internal boundary conditions:

The outer boundary conditions:

The initial conditions:

The following dimensionless variables were defined:

Dimensionless form:

Laplace transformation:

The general form of the equation: where

#### D. Flow in Microfracture System

The mass conservation equation for natural fracture systems is as follows:

The internal boundary conditions:

The outer boundary conditions:

The initial conditions:

The following dimensionless variables were defined:

Dimensionless form:

Substituting internal and outer boundary conditions into the above general form: where

#### Data Availability

This article does not have any underlying data.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.