**2. Method**

Production of shale gas requires horizontal wells with multistage hydraulic fracturing. The propagation of fractures is uncertain and branch fractures will be created around the main hydraulic fractures, resulting in a stimulated region around each main hydraulic fracture. In order to model such a complex system, the reservoir is simplified as a three-region system, where the first region is the hydraulic fracture region which is regarded as the sole connection to the well. The second one represents the stimulated region with aggregated volume of all the microfractures and the third one is the adjacent ultra-low-permeability matrix directly connected to the stimulated region. Figure 1a shows the schematic 3D model. Three regions are contained in the model: region 1 is the hydraulic fracture, region 2 (darker color) with higher permeability around each hydraulic fracture, and region 3 (lighter color) with lower permeability connected to the region 2. The arrows represent the flow directions. For this model, the flow directions are parallel, and the system is symmetrical with respect to the hydraulic fracture and horizontal-well. Thus, it is feasible to use one quarter of the reservoir shown in Figure 1b to replace the whole reservoir in order to simplify the derivation process. According to Anderson et al. [33], they verified that when the permeability of region 2 is less than or equal to 500 nD, the contribution of the region beyond fractures can be neglected for the 20-years production. For the case that the distance from the fracture face to permeability boundary (*x*1) is less than the half distance between fractures (*x*2), the contribution would be smaller. Meanwhile, Stalgorova et al. [21] also set numerical models to illustrate the negligible contribution of region beyond fractures, and they obtained that the di fference of 20-years production is negligible after comparing the results of numerical simulation with and without region beyond fractures. In the work of Heidari et al. [22], they also did not take the region beyond fractures into account. Therefore, the contribution of the region beyond the tips of the hydraulic fractures is assumed to be negligible in this work.

**Figure 1.** Schematic of a multistage fractured horizontal well with two regions around the hydraulic fracture. (**a**) The 3D model. (**b**) The plan view.

In this work, our analytical model is derived under the following assumptions:


#### *2.1. Gas Adsorption*/*Desorption E*ff*ect*

In contrast to conventional gas reservoirs, gas adsorption is an important feature of shale gas reservoir. The Langmuir isotherm adsorption equation [34] is widely used to calculate the shale gas adsorption and its expression is as follows:

$$V = V\_L \frac{P}{P\_L + P} \tag{1}$$

where *V* is the volume of the adsorbed gas, and *P* is pressure. *VL* and *PL* stand for Langmuir volume and pressure, respectively. When considering the gas adsorption, the effect of adsorption on compressibility of the reservoir is essential. According to Bumb et al. [35], the new compressibility equation can be expressed as,

$$\mathcal{C}\_{I}"=\mathcal{C}\_{f}+\mathcal{C}\_{w}\mathcal{S}\_{w}+\mathcal{C}\_{\mathcal{S}}(1-\mathcal{S}\_{w})+\mathcal{C}\_{\mathcal{S}^{d}}\tag{2}$$

$$C\_{\mathcal{S}^d} = \frac{0.031214 \rho\_m V\_L \overline{B\_{\mathcal{S}}} P\_L}{\phi \left(\overline{p} + P\_L\right)^2} \tag{3}$$

where *Cf* is rock compressibility, *Cw* is water compressibility, *Cg* is free gas compressibility and *Cgd* is adsorbed gas compressibility, *Sw* is water saturation, ρ*m* is matrix density, *Bg* is gas reservoir volume factor, and φ is the porosity.

Another important change is that the compressibility factor is modified by King [36],

$$z^\* = \frac{z}{(1 - S\_{\rm av}) + 0.031214 \rho\_{\rm m} \frac{zT}{\phi} \Big(\frac{p\_{\rm sc}}{T\_{\rm sc}}\Big) V\_L \frac{1}{p + p\_L}}\tag{4}$$

where *Psc* is the standard condition pressure, *Tsc* is the standard condition temperature, *z* is compressibility factor, *T* is the reservoir temperature.

#### *2.2. High-Velocity Non-Darcy Flow E*ff*ect*

For gas flow in the hydraulic fracture, high-velocity non-Darcy effect is considered in this study. Forchheimer [37] proposed that Darcy's law is inadequate to describe high-velocity gas flow without adding an inertial effect, which is proportional to the square of the flow velocity. To account for the non-Darcy flow effect, an inertial term must be included. The Forchheimer's flow equation is given as,

$$
\overrightarrow{v} - \nabla p = \frac{\mu}{k}\overrightarrow{v} + \beta\rho \left|\overrightarrow{v}\right|\overrightarrow{v} \tag{5}
$$

In order to reduce the nonlinearity, the equivalent permeability is introduced to obtain the extended Darcy's law [37],

$$-\nabla p = \frac{\mu}{k\_{c\eta}}\vec{v}\tag{6}$$

Substituting Equation (6) into Equation (5), the equivalent permeability of hydraulic fracture yields,

$$k\_{\rm eq}^F = \frac{\mu k\_{\rm iF}}{\mu + \beta \rho \left| \overrightarrow{\boldsymbol{\upsilon}} \right| k\_{\rm iF}} \tag{7}$$

where

$$\beta = \frac{4.1 \times 10^{11}}{(k\_{\text{iF}})^{1.5}} \tag{8}$$

among which, β is the non-Darcy flow coefficient, *kiF* is the hydraulic fracture permeability, *kFeq* is the equivalent permeability of hydraulic fracture, *v* is the flow velocity, and ρ is the gas density.

#### *2.3. Non-Darcy Flow in the Matrix*

For the nano-porous media in shale reservoir, Darcy's law has difficulty describing the actual gas behavior. The gas flow can be classified as four flow conditions, such as continuum flow, slip flow, transition flow, and free-molecular flow. According to the previous publications [38], the four flow conditions can be qualified by the Knudsen number. However, the Knudsen number varies from 10−<sup>3</sup> to 1 in most shale reservoirs. In order to represent the non-Darcy gas flow in matrix, the apparent permeability is presented by the following general form:

$$k\_{\rm n}^{\rm m} = k\_{\rm co}^{\rm m} f(Kn) \tag{9}$$

where *km*∞ is the intrinsic permeability of the porous medium, *f*(*Kn*) is the correlation term expressed as a function of the Knudsen number, which is modeled as [38],

$$f(Kn) = (1 + aKn) \left( 1 + \frac{4Kn}{1 + Kn} \right) \tag{10}$$

in which,

$$a = \frac{128}{15\pi^2} \tan^{-1} \left( 4Kn^{0.4} \right) \tag{11}$$

Meanwhile, for a capillary tube of radius, *r*, the intrinsic permeability can be derived [39],

$$k^m\_{\phi\phi} = \frac{r^2}{8} \tag{12}$$

The Knudsen number *Kn* is defined as the ratio of the molecular mean free path length and the pore radius in shale matrix,

$$Kn = \frac{\lambda}{r} \tag{13}$$

The mean free path can also be calculated,

$$
\lambda = \frac{\mu\_{\mathcal{S}}}{P} \sqrt{\frac{\pi RT}{2M}} \tag{14}
$$

Substituting Equation (14) into Equation (13) on basis of the real gas condition, we can obtain:

$$\mathbf{K}n = \frac{\mu\_{\mathcal{S}} z}{rP} \sqrt{\frac{\pi RT}{2M}}\tag{15}$$

where μ*g* is the gas viscosity, *z* is compressibility factor, *r* is the pore radius, *P* is the reservoir pressure, *R* is the universal gas constant, *T* is the reservoir temperature, and *M* is the gas molecular weight.

#### *2.4. Derivation of Linearized Gas Di*ff*usivity Equation.*

For the flow of shale gas, the gas diffusivity equation will be nonlinear which makes deriving the analytical solution difficult. On one hand, with the reduction in average reservoir pressure, the gas properties like the gas viscosity (μ), gas compressibility (*Ct*), and gas compressibility factor (*z*) will change with the pressure. On the other hand, when incorporating the significant mechanisms in shale gas reservoir like gas adsorption and non-Darcy flow effect, the permeability is a variate with pressure rather than a constant.

In order to deal with this problem, the pseudo-pressure and pseudo-time instead of the pressure and time are adopted to linearize the equations. We set a general real gas diffusivity equation in three-dimensional Cartesian coordinate system as example. When the non-Darcy effect is coupled, the

*k*(*p*) can be calculated by Equation (7) or Equation (9). Certainly, the *ct*(*p*) can be replaced by *ct*<sup>∗</sup>(*p*) to represent the gas adsorption effect.

$$\frac{\partial}{\partial \mathbf{x}} \Big( k(p) \frac{p}{\mu(p) z(p)} \frac{\partial p}{\partial \mathbf{x}} \Big) + \frac{\partial}{\partial y} \Big( k(p) \frac{p}{\mu(p) z(p)} \frac{\partial p}{\partial y} \Big) + \frac{\partial}{\partial z} \Big( k(p) \frac{p}{\mu(p) z(p)} \frac{\partial p}{\partial z} \Big) = \phi \mathbf{c}\_t(p) \frac{p}{z(p)} \frac{\partial p}{\partial t} \tag{16}$$

where *k*(*p*) is the pressure-dependent permeability, μ(*p*) is the pressure-dependent viscosity, *ct*(*p*) is the pressure-dependent compressibility.

Considering the pressure dependence of the permeability, we define a general modified pseudo-pressure transformation. Surely, *k*(*p*) can be calculated by Equation (7) or Equation (9) to couple with the non-Darcy effect. *z*(*p*) can be replaced with *z*<sup>∗</sup>(*p*) calculated by Equation (4) to consider the gas adsorption.

$$m(p) = \frac{1}{k\_i} \int\_{p\_b}^{p} k(p) \frac{p}{\mu(p)z(p)} dp \tag{17}$$

where *m*(*p*) is the pseudo-pressure, *ki* is the intrinsic permeability, *z*(*p*) is the pressure-dependent compressibility factor.

Substituting Equation (17) into Equation (16), the right side of the partial differential equation is still nonlinear.

$$\frac{\partial^2 m(p)}{\partial x^2} + \frac{\partial^2 m(p)}{\partial y^2} + \frac{\partial^2 m(p)}{\partial z^2} = \frac{\phi \mu(p) c\_t(p)}{k(p)} \frac{\partial m(p)}{\partial t} \tag{18}$$

To linearize Equation (18), a general modified pseudo-time transformation is defined [24].

$$t\_d = \frac{1}{k\_i} \int\_0^t k(p) \frac{(\mu c\_t)\_i}{\mu(p)c\_t(p)} dt \tag{19}$$

where *ta* is the pseudo-time, μ*i* is the initial viscosity, *cti* is the initial compressibility.

After substituting the Equation (19) into Equation (18), the general linear partial differential equation is derived as,

$$\frac{\partial^2 m(p)}{\partial x^2} + \frac{\partial^2 m(p)}{\partial y^2} + \frac{\partial^2 m(p)}{\partial z^2} = \frac{\phi \mu\_i c\_{ti}}{k\_i} \frac{\partial m(p)}{\partial t\_a} \tag{20}$$

Considering the permeability and compressibility in three regions are different, therefore, the pseudo-time in three regions can be shown as follows.

$$t\_{aF} = \frac{1}{k\_{iF}} \int\_0^t k\_{c\eta}^F(p) \frac{(\mu c\_t)\_{iF}}{\mu\_F(p)c\_{tF}(p)} dt \tag{21}$$

$$t\_{d1} = \frac{1}{k\_{i1}} \int\_0^t k\_1(p) \frac{(\mu c\_t)\_{i1}}{\mu\_1(p)c\_{t1}(p)} dt \tag{22}$$

$$t\_{d2} = \frac{1}{k\_{i2}} \int\_0^t k\_{i\}^m(p) \frac{(\mu c\_t^\*)\_{i2}}{\mu\_2(p) c\_{t2}^\*(p)} dt \tag{23}$$

where *taF*, *taF*, *ta*2 are the pseudo-time in the hydraulic fracture region, region 1, and region 2 respectively.

Therefore, we can obtain the linearized gas diffusivity equation in three regions respectively. Firstly, we consider the gas adsorption and gas slippage in matrix and the high-velocity non-Darcy flow in hydraulic fracture region. Thus, we use a new total compressibility coupling the gas adsorption effect in diffusivity equation of region 2. Finally, through the linearization by the modified pseudo-pressure

and pseudo-time, the non-Darcy flow effect for matrix and fracture are included in the governing linear equations in different regions.

#### *2.5. Model Description in Di*ff*erent Regions.*

#### 2.5.1. Model Description in Matrix Region (Region 2)

The system of equations based on the conceptual model is presented as follows. For the lowpermeability matrix region (region 2), the diffusivity equation for gas flow is derived as,

$$\frac{\partial^2 m\_2(p)}{\partial \mathbf{x}^2} + \frac{\partial^2 m\_2(p)}{\partial y^2} + \frac{\partial^2 m\_2(p)}{\partial z^2} = \frac{(\phi \mu\_i \mathbf{c}\_{i1}^\*)\_2}{k\_{i2}} \frac{\partial m\_2(p)}{\partial t\_{i2}}\tag{24}$$

where *<sup>m</sup>*2(*p*) is the pseudo-pressure in region 2, *x*, *y*, *z* are the Cartesian coordinates, *ta*2 is the pseudo-time for gas flow, *<sup>c</sup>*<sup>∗</sup>*ti*2 and μ*i*2 are the initial modified total compressibility and viscosity in region 2, *ki*2 and φ*i*2 are the permeability and porosity in region 2, respectively.

The initial condition for the region is that the pseudo-pressure of the region is equal to initial pseudo-pressure at *t* = 0.

$$m\_{\mathbf{2}}(p)(\mathbf{x}, y, z, 0) = m(p\_{\bar{i}}) \tag{25}$$

where *m*(*pi*) is the initial pseudo-pressure.

The boundary conditions are defined as no-flow at top and bottom of the reservoir. Additionally, the both ends of y-direction can also be regarded as no-flow boundaries.

$$\left. \frac{\partial m\_2(p)}{\partial y} \right|\_{y=0, x\_f} = 0 \tag{26}$$

$$\left.\frac{\partial m\_2(p)}{\partial z}\right|\_{z=z\_0, z\_\varepsilon} = 0\tag{27}$$

Due to the plane of symmetry between adjacent fractures, the location *x* = *x*2 is also a no-flow boundary.

$$\left.\frac{\partial m\_2(p)}{\partial \mathbf{x}}\right|\_{\mathbf{x}=\mathbf{x}\_2} = 0\tag{28}$$

The continuity of flux and pressure across the boundaries between the regions is assumed. In Ogunyomi's [27] and Qiu's [29] work, only one is chosen as a boundary condition. Considering the average pseudo-pressure will be adopted, the continuity of flux is chosen as the last boundary condition which can be given by,

$$k\_{l2} \frac{\partial m\_2(p)}{\partial \mathbf{x}} \bigg|\_{\mathbf{x} = \mathbf{x}\_1} = k\_{l1} \frac{\partial m\_1(p)}{\partial \mathbf{x}} \bigg|\_{\mathbf{x} = \mathbf{x}\_1} \tag{29}$$

#### 2.5.2. Model Description in Stimulated Region Volume (Region 1)

For the stimulated region volume (region 1), the diffusivity equation for gas flow is also expressed as,

$$\frac{\partial^2 m\_1(p)}{\partial \mathbf{x}^2} + \frac{\partial^2 m\_1(p)}{\partial y^2} + \frac{\partial^2 m\_1(p)}{\partial z^2} = \frac{(\phi \mu\_i c\_{ii})\_1}{k\_{i1}} \frac{\partial m\_1(p)}{\partial t\_{a1}} \tag{30}$$

where *<sup>m</sup>*1(*p*) is the pseudo-pressure in region 1, *x*, *y*, *z* are the Cartesian coordinates, *ta*1 is the pseudo-time for gas flow, *cti*1 and μ*i*1 are the initial total compressibility and viscosity in region 1, *ki*1 and φ1 are the permeability and porosity in region 1, respectively.

For the whole model, the initial condition remains the same.

$$m\_1(p)(\mathbf{x}, y, z, 0) = m(p\_i) \tag{31}$$

The flow in region 1 is also the 1D linear in x-direction, therefore, the outer boundary conditions are identical with those of region 2.

$$\left. \frac{\partial m\_1(p)}{\partial y} \right|\_{y=0, x\_f} = 0 \tag{32}$$

$$\left.\frac{\partial m\_1(p)}{\partial z}\right|\_{z=z\_0, z\_d} = 0\tag{33}$$

According to the continuity hypothesis, the boundary condition is expressed as,

$$k\_{l1} \frac{\partial m\_1(p)}{\partial \mathbf{x}} \bigg|\_{\mathbf{x} = x\_1} = k\_{l2} \frac{\partial m\_2(p)}{\partial \mathbf{x}} \bigg|\_{\mathbf{x} = x\_1} \tag{34}$$

$$k\_{i1} \frac{\partial m\_1(p)}{\partial \mathbf{x}} \Big|\_{\mathbf{x} = \frac{\mathbf{y}}{2}} = k\_{i\mathbf{F}} \frac{\partial m\_{\mathbf{F}}(p)}{\partial \mathbf{x}} \Big|\_{\mathbf{x} = \frac{\mathbf{y}}{2}} \tag{35}$$

2.5.3. Model Description in Hydraulic Fracture Region

For the hydraulic fracture region, the diffusivity equation for gas flow is also expressed as,

$$\frac{\partial^2 m\_F(p)}{\partial \mathbf{x}^2} + \frac{\partial^2 m\_F(p)}{\partial y^2} + \frac{\partial^2 m\_F(p)}{\partial z^2} = \frac{(\phi \mu\_i c\_{i\bar{t}})\_F}{k\_{i\bar{F}}} \frac{\partial m\_F(p)}{\partial t\_{a\bar{F}}} \tag{36}$$

where *mF*(*p*) is the pseudo-pressure in hydraulic fracture region, *x*, *y*, *z* are the Cartesian coordinates, *taF* is the pseudo-time for gas flow, *ctiF* and μ*iF* are the initial total compressibility and viscosity in hydraulic fracture region, *kiF* and φ*F* are the permeability and porosity in hydraulic fracture region, respectively.

For the whole model, the initial condition remains the same.

$$m\_F(p)(x, y, z, 0) = m(p\_i) \tag{37}$$

With the assumption of the constant bottom-hole pressure, at the location of *y* = 0, the pressure is the same as the bottom-hole pressure at any time.

$$m\_F(p)(\mathbf{x}, y = 0, z\_\prime t\_\mathbf{a}) = m(p\_{wf}) \tag{38}$$

The flow in hydraulic fracture region is also the 1D linear in the y-direction, therefore, the outer boundary conditions is expressed as,

$$\left.\frac{\partial m\_F(p)}{\partial x}\right|\_{x=0} = 0\tag{39}$$

$$\left.\frac{\partial m\_F(p)}{\partial y}\right|\_{y=x\_f} = 0\tag{40}$$

$$\left.\frac{\partial m\_F(p)}{\partial z}\right|\_{z=z\_0, z\_\epsilon} = 0\tag{41}$$

According to the continuity hypothesis, the boundary condition is expressed as,

$$k\_{\rm iF} \frac{\partial m\_F(p)}{\partial \mathbf{x}} \Big|\_{\mathbf{x} = \frac{\mathbf{w}}{2}} = k\_{i1} \frac{\partial m\_1(p)}{\partial \mathbf{x}} \Big|\_{\mathbf{x} = \frac{\mathbf{w}}{2}} \tag{42}$$
