**Nomenclature**



#### **Appendix A Derivation of the Average Pseudo-Pressure**

The solution to the 1-D linear oil flow with constant-pressure boundary condition is given by [27]:

$$p\_D = 1 + \frac{4}{\pi} \sum\_{n=1}^{\infty} \frac{(-1)^{l}}{2n - 1} e^{-\left(\frac{(2n-1)\pi}{2}\right)^2 l\_D} \cos\left(\frac{(2n-1)\pi}{2} x\_D\right) \tag{A1}$$

For the shale gas case, the pressure and time should be replaced by pseudo-pressure and pseudo-time shown in Equations (17) and (19). Following the same step, we derive the general solution applied for the 1-D linear gas flow with constant-pressure boundary condition.

$$m\_D(p) = 1 + \frac{4}{\pi} \sum\_{n=1}^{60} \frac{(-1)^n}{2n - 1} e^{-\left(\frac{(2n-1)\pi}{2}\right)^2 t\_{dD}} \cos\left(\frac{(2n-1)\pi}{2} \mathbf{x}\_D\right) \tag{A2}$$

where

*mD* = *<sup>m</sup>*(*pi*)−*<sup>m</sup>*(*p*) *<sup>m</sup>*(*pi*)−*<sup>m</sup>*(*pw f*) is the dimensionless reservoir pseudo-pressure. *taD* = *kta* is the dimensionless pseudo-time.

φμ*ctL xD* = *xe*−*x xe*−*xw f* is the dimensionless distance in the x-direction.

Similarly, the production rate can be obtained as,

$$
\eta\_D = 2e^{-\left(\frac{(2n-1)\pi}{2}\right)^2 t\_{dD}} \tag{A3}
$$

According to the definition of volume-weighted average pseudo-pressure,

$$
\overline{m}(p) = \frac{\int\_V m(p)dV}{\int\_V dV} \tag{A4}
$$

Therefore,

$$\overline{m}\_{\rm D}(p) = \frac{\int\_{1}^{\chi\_{\rm D}} 1 + \frac{4}{\pi} \sum\_{n=1}^{\infty} \frac{\left(\frac{-1}{2n}\right)^{n} e^{-\left(\frac{(2n-1)\pi}{2}\right)^{2} t\_{\rm dD}}}{\int\_{1}^{\chi\_{\rm D}} d\chi\_{\rm D}}} \times \frac{\left(\frac{(2n-1)\pi}{2} \chi\_{\rm D}\right) d\chi\_{\rm D}}{} \tag{A5}$$

Performing the integration results in,

$$\overline{m}\_{\rm D}(p) = 1 + \frac{8}{\pi^2} \sum\_{n=1}^{\infty} \frac{(-1)^n}{\left(2n-1\right)^2} e^{-\frac{(2n-1)^2 \pi^2}{4} t\_{\rm dD}} \sin\frac{(2n-1)\pi}{2} \tag{A6}$$

For odd values of n, sin (<sup>2</sup>*n*−<sup>1</sup>)<sup>π</sup> 2 is positive unity and (−<sup>1</sup>)*<sup>n</sup>* is negative unity. As a result, the product of these two coefficients is always negative. Equation (A6) can be rewritten as,

$$\overline{m}\_{\mathbb{D}}(p) = 1 - \frac{8}{\pi^2} \sum\_{n=1}^{\infty} \frac{1}{\left(2n - 1\right)^2} e^{-\frac{(2n-1)^2 \pi^2}{4} t\_{d\mathbb{D}}} \tag{A7}$$

Substituting the Equation (A3) into the Equation (A7), we can obtain that,

$$\overline{m}\_{\rm o}(p) = 1 - \frac{4}{\pi^2} \sum\_{n=1}^{\infty} \frac{qn}{(2n-1)^2} \tag{A8}$$

Transforming the Equation (A8) into dimensional form, Equation (A8) can be rewritten as,

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

#### **Appendix B Solution of the System of ODEs**

To solve the systems of ODEs in Equation (71), we extract the coefficient matrix as,

$$A\_{ll} = \begin{pmatrix} -\left(\frac{1}{\tau\_2} + \frac{T\_{\tau 1}}{\tau\_1 T\_{\tau 1 F}}\right) & \frac{T\_{\tau 2 1}}{\tau\_1 T\_{\tau 1 F}} & 0\\ \frac{1}{\tau\_1} & -\left(\frac{1}{\tau\_1} + \frac{T\_{\tau 1 F}}{\tau\_F f}\right) & \frac{T\_{\tau 1 F}}{\tau\_F f} \\ 0 & \frac{1}{\tau\_F} & -\frac{1}{\tau\_F} \end{pmatrix}\_{3 \times 3} \tag{A10}$$

For the 3 × 3 matrix, we can obtain three eigenvalues λ1, λ2, λ3 and three eigenvectors ξ1, ξ2, ξ3 where

$$\begin{aligned} \, \, \, \xi\_1 = \begin{pmatrix} a\_1 \\ a\_2 \\ a\_3 \end{pmatrix} \xi\_2 = \begin{pmatrix} a\_4 \\ a\_5 \\ a\_6 \end{pmatrix} \xi\_3 = \begin{pmatrix} a\_7 \\ a\_8 \\ a\_9 \end{pmatrix} \end{aligned} \tag{A11}$$

Combining the eigenvalues and eigenvectors, we obtain the solution of the ODEs,

$$\mathbb{E}\begin{pmatrix} q\_{2n} \\ q\_{1n} \\ q\_{\mathbb{P}n} \end{pmatrix} = (2n-1)^2 \left[ \mathbb{C}\_1 \begin{pmatrix} a\_1 \\ a\_2 \\ a\_3 \end{pmatrix} \mathbb{s}^{(2n-1)^2 \lambda\_1 t\_s} + \mathbb{C}\_2 \begin{pmatrix} a\_4 \\ a\_5 \\ a\_6 \end{pmatrix} \mathbb{s}^{(2n-1)^2 \lambda\_2 t\_s} + \mathbb{C}\_3 \begin{pmatrix} a\_7 \\ a\_8 \\ a\_9 \end{pmatrix} \mathbb{s}^{(2n-1)^2 \lambda\_3 t\_s} \right] \tag{A12}$$

where *C*1, *C*2, and *C*3 are unknown constants.

> According to the initial assumption in Equations (72)–(74), we can calculate the value of *C*1, *C*2, and

$$\mathbf{C}\_3 \mathbf{C}\_1 = \beta\_3 q\_{\mathrm{iF}} \mathbf{C}\_2 = -\beta\_2 q\_{\mathrm{iF}} \mathbf{C}\_3 = \beta\_1 q\_{\mathrm{iF}} \tag{A13}$$

Where

$$\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)}$$

$$\beta\_2 = \frac{a\_1a\_8 - a\_2a\_7}{a\_1a\_5 - a\_2a\_4}\beta\_1$$

$$\beta\_3 = \frac{\beta\_2a\_4 - a\_7}{a\_1}\beta\_1$$

Therefore, we obtain the function of the final rate,

$$q\_{\rm Fn} = (2n - 1)^2 \Big( \beta\_3 a\_3 q\_{\rm \parallel \rm F} e^{\left(2n - 1\right)^2 \lambda\_1 t\_4} - \beta\_2 a\_6 q\_{\rm \parallel \rm F} e^{\left(2n - 1\right)^2 \lambda\_2 t\_4} + \beta\_1 a\_8 q\_{\rm \parallel \rm F} e^{\left(2n - 1\right)^2 \lambda\_3 t\_4} \Big) \tag{A14}$$

To obtain the total production rate, we summate the Equation (A14),

$$\begin{array}{ll} q\_F = & \sum\_{\substack{\alpha=0\\n=1\\n\equiv 1\end{array}}^{\infty} q\_{F\alpha} = \beta\_3 a\_3 q\_{\parallel F} e^{\lambda\_1 t\_d} - \beta\_2 a\_6 q\_{\parallel F} e^{\lambda\_2 t\_d} + \beta\_1 a\_9 q\_{\parallel F} e^{\lambda\_3 t\_d} \\ & + \sum\_{\substack{\alpha=0\\n\equiv 1\pmod{4}}}^{\infty} (2n-1)^2 \left(\beta\_3 a\_3 q\_{\parallel F} e^{\left(2n-1\right)^2 \lambda\_1 t\_d} - \beta\_2 a\_6 q\_{\parallel F} e^{\left(2n-1\right)^2 \lambda\_2 t\_d} + \beta\_1 a\_9 q\_{\parallel F} e^{\left(2n-1\right)^2 \lambda\_3 t\_d}\right) \end{array} \tag{A15}$$

Through converting the summation to an integral,

$$\begin{array}{lcl} q\_{\var|F} &=& \beta\_{3} a\_{3} q\_{\var|F} \epsilon^{\lambda\_{1} t\_{4}} - \beta\_{2} a\_{6} q\_{\var|F} \epsilon^{\lambda\_{2} t\_{4}} + \beta\_{1} a\_{9} q\_{\var|F} \epsilon^{\lambda\_{3} t\_{4}} \\ &+ \lim\_{z \to \ast \infty} \int\_{3}^{z} \left(\beta\_{3} a\_{3} q\_{\var|F} \epsilon^{-z^{2} t\_{4}} - \beta\_{2} a\_{6} q\_{\var|F} \epsilon^{\frac{\lambda\_{2}}{\lambda\_{1}} t\_{4}} + \beta\_{1} a\_{9} q\_{\var|F} \epsilon^{\frac{\lambda\_{3}}{\lambda\_{1}} t\_{4}}\right) \frac{dz}{2\sqrt{|\lambda\_{1}| t\_{4}}} \end{array} \tag{A16}$$

Adopting the complementary error function, the rate versus pseudo-time function can be

$$\begin{array}{ll} q\_{\overline{F}} = & \beta\_{\overline{3}} a\_{3} q\_{\overline{\mu}} e^{\lambda\_{1} t\_{a}} - \beta\_{\overline{2}} a\_{6} q\_{\overline{\mu}} e^{\lambda\_{2} t\_{a}} + \beta\_{1} a\_{9} q\_{\overline{\mu}} e^{\lambda\_{3} t\_{d}} + \frac{\beta\_{3} a\_{5} q\_{\overline{\mu}} \sqrt{\pi}}{4 \sqrt{|\lambda\_{1} t\_{d}|}} \text{erfc} \Big( 3 \sqrt{|\lambda\_{1} t\_{d}|} \Big) \\\ & - \frac{\beta\_{2} a\_{9} q\_{\overline{\mu}} \sqrt{\pi}}{4 \sqrt{|\lambda\_{2} t\_{d}|}} \text{erfc} \Big( 3 \sqrt{|\lambda\_{2} t\_{d}|} \Big) + \frac{\beta\_{1} a\_{9} q\_{\overline{\mu}} \sqrt{\pi}}{4 \sqrt{|\lambda\_{3} t\_{d}|}} \text{erfc} \Big( 3 \sqrt{|\lambda\_{3} t\_{d}|} \Big) \end{array} \tag{A17}$$
