**3. Results**

#### *3.1. Derivation of Analytical Solution.*

The diffusivity equations in the mathematical model are all PDEs. In order to obtain the analytical solution, the first step is to transform the sets of equations into ODEs. In this work, we adopted the integral method other than the common Laplace Transform. By integrating over the spatial domain, the spatial dependence is eliminated, which is even feasible to irregular regions as will be demonstrated later. The pressure in different regions is assumed as average value in our model. Therefore, the pseudo-time is supposed to be independent with space. Integrating the equations with respect to spatial coordinates,

$$\begin{split} & \int\_{x\_1}^{x\_2} \int\_{\sigma}^{\tau\_f} \int\_{z\_o}^{\frac{\partial}{\partial x}} (\frac{\partial m\_2(p)}{\partial x}) dx dy dz + \int\_{x\_1}^{x\_2} \int\_{\sigma}^{\tau\_f} \int\_{z\_o}^{\frac{\partial}{\partial y}} (\frac{\partial m\_2(p)}{\partial y}) dx dy dz + \int\_{x\_1}^{x\_2} \int\_{\sigma}^{\tau\_f} \int\_{z\_o}^{\frac{\partial}{\partial z}} (\frac{\partial m\_2(p)}{\partial z}) dx dy dz \\ &= \frac{(\phi \mu\_r \epsilon\_n^\*)\_{\frac{\alpha}{2}}}{k\_{\Omega}} \frac{\partial}{\partial t\_{\alpha 2}} \int\_{\mathbf{x}\_1}^{\mathbf{x}\_2 \times \mathbf{y}\_1} \int\_{z\_o}^{\mathbf{x}\_2 \times \mathbf{y}\_2} \int\_{\mathbf{x}\_3}^{\mathbf{x}\_4} \int\_{z\_o}^{\mathbf{x}\_2} m\_2(p) dx dy dz \end{split} \tag{43}$$

In Equation (43), we move the pseudo-time outside the spatial integral since the pseudo-time is independent of the spatial coordinates. In order to ge<sup>t</sup> a simplified equation, the average pseudo-pressure and the effective pore volume is defined as,

$$
\overline{m}(p) = \frac{\iiint m(p)}{\iiint dx dy dz} = \frac{\iiint m(p)}{V\_b} \tag{44}
$$

$$V\_p = \phi V\_b \tag{45}$$

where *Vb* is the volume of the region and *Vp* is the pore volume of the region.

Equation (43) can be rewritten as,

$$\begin{split} & \int\_{0}^{\mathbf{x}\_{f}} \int\_{\mathbf{x}}^{\mathbf{z}\_{f}} \left( \left. \frac{\partial w\_{2}(\mathbf{p})}{\partial \mathbf{x}} \right|\_{\mathbf{x}} - \left. \frac{\partial w\_{2}(\mathbf{p})}{\partial \mathbf{x}} \right|\_{\mathbf{x}} \right) dydz + \int\_{\mathbf{x}\_{1}}^{\mathbf{x}\_{2}} \int\_{\mathbf{x}}^{\mathbf{z}\_{f}} \left( \left. \frac{\partial w\_{2}(\mathbf{p})}{\partial \mathbf{y}} \right|\_{\mathbf{x}\_{f}} - \left. \frac{\partial w\_{2}(\mathbf{p})}{\partial \mathbf{y}} \right|\_{\mathbf{0}} \right) d\mathbf{x} d\mathbf{z} + \int\_{\mathbf{x}\_{1}}^{\mathbf{x}\_{2}} \int\_{0}^{\mathbf{y}} \left( \left. \frac{\partial w\_{2}(\mathbf{p})}{\partial \mathbf{z}} \right|\_{\mathbf{z}\_{\ell}} - \left. \frac{\partial w\_{2}(\mathbf{p})}{\partial \mathbf{y}} \right|\_{\mathbf{z}\_{0}} \right) d\mathbf{x} d\mathbf{y} \\ &= \frac{(\phi\_{1} \omega\_{\text{s}}^{\*})\_{2} V\_{b\_{2}}}{k\_{\text{s}}} \frac{d\mathbf{m}\_{2}(\mathbf{p})}{dt\_{\ell}} \end{split} \tag{46}$$

Using the initial condition and boundary conditions, Equation (46) is simplified to

$$-\int\_{0}^{\mathbf{x}\_{f}}\int\_{z\_{0}}^{z\_{e}}\left(k\_{i2}\frac{\partial m\_{2}(p)}{\partial y}\bigg|\_{x\_{1}}\right)dydz = (V\_{p}\mu\_{i}\mathbf{c}\_{fi}^{\*})\_{2}\frac{d\overline{m}\_{2}(p)}{dt\_{d2}}\tag{47}$$

According to the equivalent Darcy's law,

$$q\_2 = \frac{T\_{sc} Z\_{sc}}{2T P\_{sc}} \int\_0^{x\_f} \int\_{z\_0}^{z\_e} k\_{i2} \frac{\partial m\_2(p)}{\partial y} \bigg|\_{x\_1} dy dz \tag{48}$$

Define

$$C\_2 = \frac{T\_{\rm sc} Z\_{\rm sc}}{2T P\_{\rm sc}} (\mu\_i c\_{ti}^\*)\_2 \tag{49}$$

Therefore, Equation (48) can be rewritten as,

$$C\_2 V\_{p2} \frac{d\overline{m\_2}(p)}{dt\_{d2}} = -q\_2 \tag{50}$$

where *Vp*2 is the pore volume of region 2, *q*2 is the flow rate in region 2.

For region 1, we also use the integration method to deal with the equation,

$$\begin{split} \frac{1}{2} \int\_{0}^{T} \left\lVert \frac{\partial w\_{1}(p)}{\partial x} \right\rVert\_{\mathbf{x}\_{1}} - \frac{\partial m\_{1}(p)}{\partial x} \Big\rVert\_{\mathbf{y}} dydz &+ \int\_{w/2}^{T} \int\_{0}^{\delta} \left\lVert \frac{\partial m\_{1}(p)}{\partial y} \right\rVert\_{\mathbf{x}\_{f}} - \frac{\partial m\_{1}(p)}{\partial y} \Big\lVert\_{0} \right] dxdz + \int\_{w/2}^{T} \int\_{0}^{T} \left\lVert \frac{\partial m\_{1}(p)}{\partial z} \right\rVert\_{\mathbf{z}\_{\ell}} - \frac{\partial m\_{1}(p)}{\partial y} \Big\lVert\_{\mathbf{z}\_{0}} \Big\lVert \mathbf{x}\_{0} \rVert\_{\mathbf{x}\_{0}} \\ &= \frac{(\phi\_{1} \mu\_{c})\_{1} V\_{b\_{1}}}{k\_{\mathrm{l}}} \frac{V\_{\mathrm{l}s\_{1}}(p)}{dt\_{\mathrm{l}1}} \end{split} \tag{51}$$

After applying the initial condition and boundary conditions, Equation (51) can be rewritten as,

$$-\int\_{0}^{\chi\_{f}}\int\_{z\_{0}}^{z\_{c}}\left(k\_{i1}\frac{\partial m\_{1}(p)}{\partial \mathbf{x}}\bigg|\_{\mathbf{x}\_{1}} - k\_{i1}\frac{\partial m\_{1}(p)}{\partial \mathbf{x}}\bigg|\_{\mathbf{x}\_{0}}\right)dydz = (V\_{p}\mu\_{i}c\_{ti})\_{1}\frac{d\overline{m}\_{i}(p)}{dt\_{d1}}\tag{52}$$

Note that,

$$q\_1 = \frac{T\_{\rm sc} Z\_{\rm sc}}{2T P\_{\rm sc}} \int\_0^{\overline{\chi}\_f} \int\_{z\_0}^{z\_\epsilon} \left( k\_{i1} \left. \frac{\partial m\_1(p)}{\partial \mathbf{x}} \right|\_{\mathbf{x}\_0} \right) dydz \tag{53}$$

Defining

$$\mathbf{C}\_{1} = \frac{T\_{\rm sc} Z\_{\rm sc}}{2T P\_{\rm sc}} (\mu\_{i} c\_{ti})\_{1} \tag{54}$$

Substituting Equations (30), (45), (49), and (50) into Equation (48) results in,

$$C\_1 V\_{p1} \frac{d\overline{m}\_1(p)}{dt\_{a1}} = -q\_1 + q\_2 \tag{55}$$

where *Vp*1 is the pore volume of region 1, *q*1 is the flow rate in region 1.

For hydraulic fracture region, we also use the integration method to deal with the di ffusion equation,

$$\begin{split} \int\_{0}^{T} \int\_{0}^{z} \left( \left. \frac{\partial m\_{f}(p)}{\partial x} \right|\_{\frac{\pi}{2}} - \left. \frac{\partial m\_{f}(p)}{\partial x} \right|\_{0} \right) dy dz + \int\_{0}^{w/2} \int\_{0}^{z} \left( \left. \frac{\partial m\_{f}(p)}{\partial y} \right|\_{x\_{f}} - \left. \frac{\partial m\_{f}(p)}{\partial y} \right|\_{0} \right) dx dz + \int\_{0}^{w/2} \int\_{0}^{z} \left( \left. \frac{\partial m\_{f}(p)}{\partial z} \right|\_{z\_{d}} - \left. \frac{\partial m\_{f}(p)}{\partial y} \right|\_{z\_{0}} \right) dx dy \\ = \frac{(\phi\_{f} \mu\_{i} c\_{i})\_{f} V\_{b\_{f}^{\*}}}{k\_{\text{J}}} \frac{d \overline{m}\_{f}(p)}{dt\_{d}} \end{split} \tag{56}$$

After applying the initial condition and boundary conditions, Equation (56) can be rewritten as,

$$\int\_{0}^{\frac{\pi}{2}} \int\_{z\_{0}}^{z\_{\varepsilon}} \int\_{0}^{z} \left(k\_{\text{if}} \frac{\partial m\_{\text{F}}(p)}{\partial x}\Big|\_{\frac{\pi}{2}}\right) dydz + \int\_{0}^{w/2} \int\_{z\_{0}}^{z} \left(-k\_{\text{if}} \frac{\partial m\_{\text{F}}(p)}{\partial y}\Big|\_{0}\right) \text{dx}dz = \left(\phi\mu\_{\text{if}}c\_{\text{if}}\right)\_{\text{F}} V\_{\text{bf}} \frac{d\overline{m}\_{\text{F}}(p)}{dt\_{\text{aff}}}\tag{57}$$

The flow rate in hydraulic fracture can be expressed as,

$$q\_F = \frac{T\_{\rm sc} Z\_{\rm sc}}{2TP\_{\rm sc}} \int\_0^{\chi\_f} \int\_0^{z\_{\rm c}} \left( k\_{\rm iF} \frac{\partial m\_F(p)}{\partial \mathbf{x}} \Big|\_0 \right) dy dz \tag{58}$$

Defining

$$\mathbf{C}\_{\rm F} = \frac{T\_{\rm sc} Z\_{\rm sc}}{2T P\_{\rm sc}} (\mu\_i c\_{ti})\_{\rm F} \tag{59}$$

Substituting Equations (42), (48), (58), and (59) into Equation (57) results in,

$$\mathcal{L}\_F V\_{pF} \frac{d\overline{m}\_F(p)}{dt\_{aF}} = -q\_F + q\_1 \tag{60}$$

The next step is to substitute the average pseudo-pressure with the relationship between pressure and the flow rate. Since it is assumed that gas flows sequentially from region 2 to region 1 then to hydraulic fracture, a general analytical solution for one-dimensional linear gas flow is derived to solve the problem (details are shown in Appendix A), which is given by:

$$\overline{m}(p) = m(p\_{wf}) + \frac{4}{\pi^2} [m(p\_i) - m(p\_{wf})] \sum\_{n=1}^{\infty} \frac{q\_{D\_n}}{\left(2n - 1\right)^2} \tag{61}$$

where *m*(*p*) is the average pseudo-pressure, *qDn* is the dimensionless production from the nth mode.

Combining the assumptions, the average pseudo-pressure in each region can be expressed as,

$$\overline{m}\_r(p) = m(p\_{wf}) + \frac{4}{\pi^2} [m(p\_i) - m(p\_{wf})] \sum\_{n=1}^{\infty} \frac{q\_{DF\_n}}{(2n-1)^2} \tag{62}$$

$$
\overline{m}\_1(p) = \overline{m}\_\mathbb{F}(p) + \frac{4}{\pi^2} [m(p\_i) - \overline{m}\_\mathbb{F}(p)] \sum\_{n=1}^\infty \frac{q\_{D1\_n}}{(2n-1)^2} \tag{63}
$$

$$\overline{m}\_2(p) = \overline{m}\_1(p) + \frac{4}{\pi^2} [m(p\_i) - \overline{m}\_1(p)] \sum\_{n=1}^{\infty} \frac{q\_{D2\_n}}{(2n-1)^2} \tag{64}$$

We define the productivity index (*J*) and transmissibility (*Tr*1*F* and *Tr*21) as,

$$J = \frac{\pi^2}{4} \frac{q\_{\text{iF}}}{m(p\_i) - m(p\_{wf})} \tag{65}$$

$$T\_{r1F} = \frac{\pi^2}{4} \frac{q\_{i1}}{m(p\_i) - \overline{m}\_1(p)}\tag{66}$$

$$T\_{r21} = \frac{\pi^2}{4} \frac{q\_{l2}}{m(p\_i) - \overline{m}\_2(p)}\tag{67}$$

where *<sup>m</sup>*1 (*p*), *<sup>m</sup>*2 (*p*) is the average pseudo-pressure in hydraulic fracture region, region 1, and region 2, respectively. *qiF* = *TscZsc* 2*TPsc kiFAF m*(*pi*) *LF* is the initial production rate from the hydraulic fracture region. *qi*1 = *TscZsc* 2*TPsc ki*1*A*1 *<sup>m</sup>*1 (*p*) *L*1 is the initial production rate from the region 1. *qi*2 = *TscZsc* 2*TPsc ki*2*A*2 *<sup>m</sup>*2 (*p*) *L*2 is the initial production rate from the region 2.

Substituting Equations (62)–(64) into Equations (50), (55), and (60), there are six ordinary differential terms on the left side of the system of equations, such as *dqFn dtaF* , *dqFn dta*1 , *dqFn dta*2 , *dq*<sup>1</sup>*n dta*1 , *dq*<sup>1</sup>*n dta*2 , *dq*<sup>2</sup>*n dta*2 . In order to solve the system of ordinary differential equations, an approximate pseudo-time . *ta* is introduced. Therefore, the Equations (50), (55), and (60) can be rewritten as follows.

$$\sum\_{n=1}^{\infty} \frac{dq\_{2\_n}}{d\tilde{t}\_n} = -\left(\frac{T\_{r21}}{\mathbb{C}\_2 V\_{p2}} + \frac{T\_{r21}}{\mathbb{C}\_1 V\_{p1}}\right) \sum\_{n=1}^{\infty} \left(2n-1\right)^2 q\_{2\_n} + \frac{T\_{r21}}{\mathbb{C}\_1 V\_{p1}} \sum\_{n=1}^{\infty} \left(2n-1\right)^2 q\_{1\_n} \tag{68}$$

$$\sum\_{n=1}^{\infty} \frac{dq\_{1\_n}}{\tilde{d}\tilde{t}\_a} = \frac{T\_{r1\mathcal{F}}}{\mathbb{C}\_1 V\_{p1}} \sum\_{n=1}^{\infty} \left( 2n - 1 \right)^2 q\_{2\_n} - \left( \frac{T\_{r1\mathcal{F}}}{\mathbb{C}\_1 V\_{p1}} + \frac{T\_{r1\mathcal{F}}}{\mathbb{C}\_\mathcal{F} V\_{p\mathcal{F}}} \right) \sum\_{n=1}^{\infty} \left( 2n - 1 \right)^2 q\_{1\_n} + \frac{T\_{r1\mathcal{F}}}{\mathbb{C}\_\mathcal{F} V\_{p\mathcal{F}}} \sum\_{n=1}^{\infty} \left( 2n - 1 \right)^2 q\_{\mathcal{F}\_n} \tag{69}$$

$$\sum\_{n=1}^{\infty} \frac{dq\_{\mathbb{F}\_n}}{d\overline{t}\_a} = -\frac{I}{\mathbb{C}\_{\mathbb{F}} V\_{p\mathbb{F}}} \sum\_{n=1}^{\infty} (2n-1)^2 q\_{\mathbb{F}\_n} + \frac{I}{\mathbb{C}\_{\mathbb{F}} V\_{p\mathbb{F}}} \sum\_{n=1}^{\infty} (2n-1)^2 q\_{1\_n} \tag{70}$$

where *qDFn* , *qD*<sup>1</sup>*n* , and *qD*<sup>2</sup>*n* are the initial production rate from the nth mode in hydraulic fracture region, region 1, and region 2, respectively. . *ta* is an approximate pseudo-time, which is defined as the average of the *ta*1, *ta*2, and *taF*.

Then we define three time-constant parameters as,

$$
\pi\_F = \frac{\mathbb{C}\_F V\_{pF}}{J} \tag{71}
$$

$$\pi\_1 = \frac{\mathbb{C}\_1 V\_{p1}}{T\_{r1F}} \tag{72}$$

$$
\tau\_2 = \frac{C\_2 V\_{p2}}{T\_{f21}} \tag{73}
$$

We can rewrite this set of ODEs in matrix form,

$$
\begin{pmatrix}
\frac{dq\_{2u}}{d\overline{t}\_{u}}\\\
\frac{dq\_{1u}}{d\overline{t}\_{d}}\\\
\frac{dq\_{F\_{0}}}{d\overline{t}\_{d}}
\end{pmatrix} = (2u - 1)^{2} \begin{pmatrix}
\frac{1}{\tau\_{1}} & -\left(\frac{1}{\tau\_{1}} + \frac{T\_{r1F}}{\tau\_{F}T}\right) & \frac{T\_{r1F}}{\tau\_{F}}\\\
0 & \frac{1}{\tau\_{F}} & -\frac{1}{\tau\_{F}}
\end{pmatrix} \begin{pmatrix}
q\_{2u} \\\ q\_{1u} \\\ q\_{F\_{0}}
\end{pmatrix} \tag{74}
$$

where the initial conditions for the system of equations are,

$$q\_F(t\_a = 0) = q\_{iF} \tag{75}$$

$$q\_1(t\_a = 0) = 0\tag{76}$$

$$q\_2(t\_a = 0) = 0\tag{77}$$

After solving Equation (70), we can ge<sup>t</sup> the *n*th flow rate in combination with the initial conditions. The flow rate is the summation of them. By converting the summation to an indefinite integral, the analytical solution in real-time space can be derived as follows (details are shown in Appendix B),

$$\begin{split} q\_{F} &= \beta\_{3} q\_{\text{iF}} a\_{3} e^{\lambda\_{1} \overleftarrow{t\_{a}}} - \beta\_{2} q\_{\text{iF}} a\_{6} e^{\lambda\_{2} \overleftarrow{t\_{a}}} + \beta\_{1} q\_{\text{iF}} a\_{9} e^{\lambda\_{3} \overleftarrow{t\_{a}}} + \frac{\beta \text{gas} \mu \sqrt{\pi}}{4 \sqrt{|\lambda\_{1} \overleftarrow{t\_{a}}|}} \text{erfc} \Big( 3 \sqrt{|\lambda\_{1} \overleftarrow{t\_{a}}|} \Big) \\ &- \frac{\beta\_{2} \text{as} q\_{\text{iF}} \sqrt{\pi}}{4 \sqrt{|\lambda\_{2} \overleftarrow{t\_{a}}|}} \text{erfc} \Big( 3 \sqrt{|\lambda\_{2} \overleftarrow{t\_{a}}|} \Big) + \frac{\beta\_{1} \text{ao} q\_{\text{iF}} \sqrt{\pi}}{4 \sqrt{|\lambda\_{3} \overleftarrow{t\_{a}}|}} \text{erfc} \Big( 3 \sqrt{|\lambda\_{3} \overleftarrow{t\_{a}}|} \Big) \end{split} \tag{78}$$

In Equation (78), the defined parameters are,

$$\beta\_1 = \frac{a\_1(a\_1a\_5 - a\_2a\_4)}{(a\_1a\_9 - a\_3a\_7)(a\_1a\_5 - a\_2a\_4) - (a\_1a\_8 - a\_2a\_7)(a\_1a\_6 - a\_3a\_4)}\tag{79}$$

$$
\beta\_2 = \frac{a\_1 a\_8 - a\_2 a\_7}{a\_1 a\_5 - a\_2 a\_4} \beta\_1 \tag{80}
$$

$$
\beta\_{\mathfrak{P}} = \frac{\beta\_2 a\_4 - a\_7}{a\_1} \beta\_1 \tag{81}
$$

where λ1, λ2, and λ3 are the eigenvalues of matrix in the Equation (74), *a*1 ∼ *a*9 are all the values in the eigenvector of matrix in the Equation (75), β1, β2, and β3 are all the coefficients.

In our model, shale gas in region 2 must firstly flow into region 1 and then into the hydraulic fracture region and finally into the wellbore. Therefore, the flow rate in hydraulic fracture is equal to that in the wellbore. From the final solution, we can find that the flow rate is related to six variables. Through fitting the production data, these variables can be obtained and then the solution can be used for production analysis and forecasting.

#### *3.2. Model Validation with Numerical Models*

In order to verify the derived analytical solution, a numerical model is built with the commercial Eclipse reservoir simulator for comparing with the previous physical model, which is one-quarter of the volume around one hydraulic fracture. The numerical model has 27 grid cells in the x-direction, 50 grid cells in the y-direction, and only one grid cell in the z-direction. In order to capture the transient flow towards the hydraulic fracture, local grid refinement in x-direction is constructed. The top view of the model is shown in Figure 2, where the first column of grids represents half-length of hydraulic fracture, and the horizontal well located in the first row of the grids along the x-direction. The blue region represents the region 2 in which the gas adsorption plays an important role, while the red one is the region 1. Table 1 summarizes the input parameters used in the numerical models, which include the reservoir conditions, hydraulic conductivity, gas adsorption, and non-Darcy flow parameters.

**Figure 2.** Reservoir grid for the numerical model. The red grids are the high-permeability region and the blue ones are the low-permeability region.


**Table 1.** Summary of input parameters for the synthetic numerical model.

In deriving the new analytical solution in this study, the pressure dependence of gas properties is considered using pseudo-pressure and pseudo-time. However, the results from the numerical models are gas rate versus real time. Therefore, the necessary step is to transform the numerical simulation results of gas rate with time into pseudo-time and then fit with the new model. Figure 3 presents the result analysis. According to the previous definition, the relationship between the pseudo-time and real time is shown in Figure 3a. The plot of 1/q vs. square root pseudo-time for regime diagnosis is shown in Figure 3b. The comparison of production rates obtained from numerical simulation and our analytical model is presented in Figure 3c. The blue dotted line represents the relationship between gas rate and pseudo-time, whereas the red one indicates that from the derived analytical solution. It can be seen that the results from the simulator and the analytical solution agree well with each other. Due to the high velocity gas flow in hydraulic facture region, the time constant in hydraulic fracture (0.001 days) is too short to be observed. Therefore, four flow regimes are identified in Figure 4. Regime 1 exhibits a half-slope straight line on the log-log plot and represents the transient linear flow in region 1. The permeability in this region is higher and hence the time constant in this region is shorter and nearly 22 days. Then an exponential curve of regime 2 indicates that the boundary of region 1 is reached, which is called the boundary-dominated flow in region 1 or inner-boundary-dominated flow. After that, the pressure propagates into the region 2. As for regime 3, it still presents an expected straight line with a half-slope signature. In our model, the permeability of region 2 is low-permeability matrix and thus the time constant is relatively longer (209 days). Regime 4 is the outer-boundary-dominated flow, which is a ffected by the boundary of region 2. Table 2 summarizes the four model parameters after fitting.

**Figure 3.** The result analysis for Case 1. (**a**) Comparison of the relationship between pseudo-time and time. (**b**) The plot 1/q vs. square root pseudo-time for regime diagnosis. (**c**) Comparison of the relationship between gas rate and pseudo-time. (Regime 1: transient linear flow in region 1 with slope -1/2. Regime 2: boundary-dominated flow for region 1. Regime 3: transient linear flow in region 2 with slope -1/2. Regime 4: boundary-dominated flow for region 2).

**Figure 4.** Reservoir grid for the numerical cases showing different shapes of stimulated region. The red grids are the stimulated region and the blue ones are the low-permeability region.

**Table 2.** Model parameters used in the analytical model for validation against numerical simulation.


Based on the output parameters, we can predict the values of physical parameters to further validate our model according to the following step-by-step procedure:

Step 1: Calculate the productivity index.

As defined above, we can calculate the productivity index *J* combined with the values of initial rate *qiF* initial pseudo-pressure *m*(*pi*) and pseudo-bottom-hole pressure *m*(*pw f*).

Step 2: Calculate the transmissibility between fractures and region 1 and region 1 and region 2. Among the output parameters, we can obtain the ratio of transmissibility (*Tr*21/*Tr*<sup>1</sup>*F*, *Tr*1*F*/*J*). Considering the calculated value in Step 1, the transmissibility can be calculated (*Tr*21, *Tr*1*F*).

Step 3: Calculate the pore volume of hydraulic fracture region, region 1, and region 2.

By transforming the formula that we defined in Equations (71)–(73), the pore volumes of different regions (*VpF*, *Vp*1, *Vp*2) can be obtained by:

$$V\_{pF} = \frac{f \tau\_F}{\mathcal{C}\_F} \ V\_{p1} = \frac{\tau\_1 T\_{r1F}}{\mathcal{C}\_1} \ V\_{p2} = \frac{\tau\_2 T\_{r21}}{\mathcal{C}\_2} \tag{82}$$

Following the above steps, we calculate the physical parameters in the numerical case as shown in Table 3 to compare with the given data. It shows that our model solution is correct within the accepted error bound. Among them, *Vg* and *Vc* express the given volume and calculated volume, respectively.


**Table 3.** Inferred parameter values according to the fitting results.

#### *3.3. Irregular Stimulated Region with One Hydraulic Fracture*

In this section, three numerical cases with only one hydraulic fracture are designed for investigating the applicability of the analytical model for irregular stimulated region in Figure 4. For the three

numerical cases, the input parameters are the same with the Case 1 except for the model dimensions. These three cases are identical with 31 grid cells in the x-direction, 51 grid cells in the y-direction, and only one grid cell in the z-direction, which represent the volume of 214 × 521 × 10 ft3. For Case 2, there is a rectangular stimulated region (region 1). As for Case 3, the stimulated region is irregular but symmetrical. In order to better describe the real condition, the stimulated region in Case 4 is designed to be neither regular nor symmetrical. The D-factor in hydraulic fracture is set as 0.0012 to represent the high-velocity non-Darcy flow.

In general, the results from the simulator and the analytical solution fit very well, as shown in Figure 5. There are also four flow regimes in the three cases. Comparing with Case 1 and Case 2, there are deviations for the cases with irregular stimulated region. The deviations are caused by the irregular inner boundary. Due to the irregular shapes of region 1, the time when the pressure propagates to inner boundaries changes with the distance away from the wellbore in di fferent parts of inner boundary. However, the analytical solution is derived based on the regular inner boundary conditions. For the same reason, the variable starting time in regime 2 results in the deviation in the early time of the third flow regime. The flow time in regime 3 is long enough, thus, the curves in this regime fit well with each other at late time. For the three cases, the outer boundaries are identical and therefore the curves in the fourth regime also agree well with each other. According to the results summarized in Table 4, we can observe that the estimated results shown in Figure 5 are acceptable within engineering accuracy. In conclusion, our analytical solution is feasible for both regular and irregular region conditions.

**Figure 5.** Comparison of the relationship of gas rate vs. pseudo-time from the numerical cases and new analytical solution. (**a**) Comparison of the relationship of gas rate vs. pseudo-time from the Case 2 and new model. (**b**) Comparison of the relationship of gas rate vs. pseudo-time from the Case 3 and new model. (**c**) Comparison of the relationship of gas rate vs. pseudo-time from the Case 4 and new model.


**Table 4.** Inferred parameter values according to fitting results.

#### *3.4. Irregular Stimulated Regions with Several Hydraulic Fractures*

In order to further validate the applicability in irregular region of the new derived model, three more numerical cases of multifractured horizontal wells are designed, as shown in Figure 6. These three cases have the identical dimension with 163 grid cells in the x-direction, 63 grid cells in the y-direction, and only one grid cell of 10 ft in the z-direction. The length of the horizontal well is 762 ft with 10 hydraulic fractures equally spaced along the x-direction. For Case 5, the stimulated regions are all the same regular and symmetrical regions. As a comparison, the stimulated regions in Case 6 are all irregular but symmetrical regions. Meanwhile, there is no interference between fractures in Case 5 and Case 6 and the length of the hydraulic fractures is 352 ft in Case 5 and 464 ft in Case 6. As for Case 7, all the stimulated regions around each hydraulic fracture are di fferent, irregular, and asymmetric. The length of 10 hydraulic fractures ranges from 288 to 432 ft and there is interference between fractures. The other input parameters are the same with Case 1.

**Figure 6.** Reservoir grid for multifractured horizontal well (MFHWs).

Figure 7 shows the comparison of gas rate versus pseudo-time obtained from numerical simulations and the new analytical solutions. It can be seen that the results agree very well. Four regimes are also identified. The linear flow time in region 1 based on the fitting results is 9, 7, and 5 days for the three cases, respectively. Then it displays the boundary flow of region 1, which lasts until about the 100th day. The time constant in region 2 is 72, 79, and 91 days, respectively. Finally, the boundary flow of region 2 lasts for hundreds of days. We can obtain that the decline rate of gas production is

faster during regimes 1 and 2 and thus the regimes 3 and 4 can last for a longer time. Therefore, for the shale gas reservoir, it is crucial to enhance the volume of stimulated regions to keep a longer high production period. Meanwhile, the contribution of gas adsorption mainly reflects in regimes 3 and 4. It helps to extend the stable production period. According to the output parameters after fitting with the numerical cases, we can calculate the volume of the hydraulic fracture region, region 1, and region 2, as shown in Table 5. Considering that the relative errors meet the requirement of engineering accuracy, our new model is also suitable for the multifractured horizontal well.

**Figure 7.** Comparison of the relationship of gas rate vs. pseudo-time from the numerical cases and new analytical solution. (**a**) Comparison of the relationship of gas rate vs. pseudo-time from the Case 5 and new model. (**b**) Comparison of the relationship of gas rate vs. pseudo-time from the Case 6 and new model. (**c**) Comparison of the relationship of gas rate vs. pseudo-time from the Case 7 and new model.


#### *3.5. Application to Field Case*

The previous section demonstrated the accuracy of the derived analytical solution for production analysis. In this section, we apply the method to history matching and forecasting of field data in a shale gas reservoir. The gas well is selected for its relatively long production history and availability of pressure data. The main work flow is as follows:


We chose a horizontal well called Well B-15 with multiple fracture stages from Barnett shale [40]. The general data of the well for analysis is listed in Table 6. The well was produced under constant pressure and even for short early variable p/q the constant pressure was also applicable [41]. Firstly, we transformed the time into pseudo-time and the relationship between pseudo-time and time is shown in Figure 8a. The plot of 1/q vs. square root pseudo-time for regime diagnosis is shown in Figure 8b. Then the log-log diagnostic plot of gas rate versus pseudo-time was obtained, which exhibits a half-slope straight line in Figure 8c. According to the previous analysis, the flow time in matrix is longer and the decline rate in region 2 is slower. Therefore, the half-slope linear flow was diagnosed as the regime 3. Following the work flow, the next step was to match our model to the production data. The results of history matching and forecasting is shown in Figure 8c. The red marks in Figure 8c represent the production data and the green ones represent the analytical solution for history matching. It indicates a good matching. The six model parameters obtained from history matching are summarized in Table 7. Based on the parameters, we forecast the gas rate represented by the black markers in Figure 8c. According to the step-by-step procedure, we calculated the volume of the hydraulic fracture region as 89 ft3, the volume of stimulated region as 84.2 MMscf, and the volume of region 2 as 1784.2 MMscf. It is essential for increasing high production period by extending the volume of the stimulated region and decreasing the decline rate in the stimulated region.

**Figure 8.** Summary of production profiles for the field example. (**a**) The relationship between time and pseudo-time. (**b**) The plot 1/q vs. square root pseudo-time for regime diagnosis. (**c**) The results of history-matching and forecasting on log-log plot.


**Table 6.** Summary of the reservoir parameters for analysis.

**Table 7.** Model parameters used in the analytical model for application in the field example.

