**5. Conclusions**

In this paper, a semianalytic fractal multi-linear flow model (FMFM) for MFHWs in unconventional gas reservoirs with consideration of the heterogeneous distribution of the properties of induced-fracture system is proposed. Fractal dimensions of induced-fracture spacing and aperture and tortuosity index of induced-fracture system are included based on fractal theory to describe the properties of SRV region. For shale gas reservoirs, gas transport mechanisms (viscous flow with slippage, Knudsen di ffusion, and surface di ffusion) among various medium including porous kerogen, inorganic matter, and fracture system are take into account in FMFM. Then, combining with ES-MDA, the FMFM can be applied not only for prediction of gas production rate for MFHWs in unconventional gas reservoirs but also for main uncertainty parameters matching. Some findings can be drawn as follows:


**Author Contributions:** Q.Z. contributed to the conception and drafting the article; S.J. contributed to the study design and manuscript editing. X.W. contributed to programming and analysis; Y.W. contributed to the literature research; Q.M. contributed to the manuscript revision. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (No.51904279), Fundamental Research Funds for the Central Universities (China University of Geosciences, Wuhan) (No. CUGGC04), and Foundation of Key Laboratory of Tectonics and Petroleum Resources (China University of Geosciences) (No. TPR-2019-03).

**Acknowledgments:** The authors thank Guanglong Sheng, School of Petroleum Engineering, Yangtze University.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Derivation of Fractal Permeability and Porosity of Induced-Fracture**

As shown in Figure 2c, using the fractal theory, the total number of induced-fracture, which is relevant to the fractal dimension of induced-fracture spacing *dfs*, can be expressed as [37]

$$\int\_{0}^{L\_{x,y}} n(L\_{x,y}) dL\_{x,y} \propto L\_{x,y}^{d\_{fs}} \tag{A1}$$

where *Lx,y* is the position of system, m. Then, Chang and Yortsos [16] presented the expression of the distribution of fracture site as

$$m\left(L\_{x,y}\right) = aL\_{x,y}^{d\_{fs}-D} \tag{A2}$$

where *D* is the Euclidean dimension, which is equal to 2 in Cartesian coordinate system; α is a constant related to fracture porosity.

For the fracture site, we can obtain

$$\text{res}\_f(L\_{\mathbf{x},\mathbf{y}})\mathfrak{n}(L\_{\mathbf{x},\mathbf{y}}) + b\_f(L\_{\mathbf{x},\mathbf{y}})\mathfrak{n}(L\_{\mathbf{x},\mathbf{y}}) = A\_{SRV\prime} \tag{A3}$$

where *ASRV* is the hydraulic fracture half-length, m.

Substituting Equation (A2) into Equation (A3), Equation (A3) can be rewritten as

$$1 + \frac{b\_f(L\_{x,y})}{s\_f(L\_{x,y})} = \frac{x\_f}{\alpha s\_f(L\_{x,y})} L\_{x,y}^{2-d\_{fs}}.\tag{A4}$$

In the actual unconventional reservoirs, the value of induced-fracture aperture is commonly far less than that of induced-fracture spacing (*bf* << *sf*), and thus Equation (A4) becomes

$$\left. \mathrm{s}\_{f}(L\_{\mathrm{x},y}) \right|\_{\frac{b\_{f}(L\_{\mathrm{x},y})}{\nu\_{f}(L\_{\mathrm{x},y})} \to 0} = \frac{\mathrm{x}\_{f}}{\alpha} L\_{\mathrm{x},y}^{2-d\_{f\*}}.\tag{A5}$$

If the induced-fracture spacing at the plane of HF is known, which can be chosen as a reference point, the FFSD in Cartesian coordinate system can be obtained by Equation (A5)

$$s\_f(y) = s\_{fw} \left(\frac{y}{y\_w}\right)^{2-d\_{fs}},\tag{A6}$$

where *sfw* is the induced-fracture spacing at reference point, m.

For the distribution of induced-fracture aperture (FFAD), the basic form can also be obtained according to Equation (A6) [21]

$$b\_f(y) = b\_{fw} \left(\frac{y}{y\_w}\right)^{d\_{fa}-2},\tag{A7}$$

where *bfw* is the induced-fracture aperture at reference point, m.

Bai and Elsworth [38] reported the expression of a fracture with a slab shape as

$$k\_{i-f} = \frac{b\_f(y)^2}{12}.\tag{A8}$$

In a representative elementary volume (REV) of fracture network, the fluid flux in the fracture system is given by

$$w = A\_r \frac{k\_f(y)}{\mu} \frac{\mathrm{d}p}{\mathrm{d}x} = n(y)A\_p \frac{k\_{i-f}(y)}{\mu} \frac{\mathrm{d}p}{\mathrm{d}L} \tag{A9}$$

where *v* is flux of induced-fracture system in a REV, m<sup>3</sup>/s; μ is the fluid viscosity, Pa·s; and d*p* is the pressure-difference of a REV in y direction, Pa.

Thus, the fracture permeability and porosity can be obtained as Equations (1) and (2).

#### **Appendix B. Solution for the FMFM Model**

For solving the FMFM, the basic dimensionless parameters are defined in Table A1.


**Table A1.** Basic parameters used in FMFM.

Subscript *i* represents the initial conditions; *f*, *m*, *k* represents fracture, inorganic matter, and porous kerogen respectively.

Therefore, transforming the FMFM into Laplace space, the dimensionless forms and solutions for HF, SRV, and USRV can be obtained as follows:

#### (1) Dimensionless forms and solutions for USRV.

According to Equations (9) and (10), the dimensionless diffusivity equations in Laplace domain for region III and V can be given by

$$\begin{cases} \overline{\rho\_{3kD}} = \frac{\sigma\_k k\_x \overline{x\_f^2}}{\sigma\_k k\_x \overline{x\_f^2} + k\_w s \omega\_k} \overline{\rho\_{3mD}}\\ \omega\_m \overline{\rho\_{3mD}} = \frac{k\_m}{k\_w} \frac{\partial^2 \overline{\rho\_{3m}}}{\partial \overline{x\_D^2}} + \sigma\_k \frac{k\_L}{k\_w} \mathbf{x}\_f^2 (\overline{\rho\_{3kD}} - \overline{\rho\_{3mD}})\\ \vdots \end{cases}, \text{ for region III,} \tag{A10a}$$

$$\begin{cases} \overline{\varphi\_{5\text{k}D}} = \frac{\sigma\_{k}k\_{k}\overline{x\_{f}^{2}}}{\alpha\_{k}k\_{k}\overline{x\_{f}^{2}} + k\_{a}\kappa\omega\_{k}}\overline{\varphi\_{5\text{m}D}}\\ \omega\_{m}\overline{\varphi\_{5\text{m}D}} = \frac{k\_{\text{m}}}{\overline{k}\_{\text{w}}}\frac{\partial^{2}\overline{\varphi\_{5\text{m}}}}{\partial x\_{D}^{2}} + \sigma\_{k}\frac{k\_{\text{k}}}{\overline{k}\_{\text{w}}}x\_{f}^{2}(\overline{\varphi\_{5\text{k}D}} - \overline{\varphi\_{5\text{m}D}}) \end{cases}, \text{ for region V.} \tag{A10b}$$

Solving Equation (A10a,b) with boundary conditions, we can obtain

$$\left. \frac{\partial \overline{\Psi^{3\text{c}}} \mathbf{u} \mathbf{D}}{\partial \mathbf{x}\_{\text{D}}} \right|\_{\mathbf{x}\_{\text{D}}=1} = \sqrt{\mathbf{u}} \mathbf{c} \left| \frac{\left\{ \mathbf{K}\_{-\left(1/2\right)} \left( \sqrt{\mathbf{u}} \mathbf{x}\_{\text{c}} \mathbf{x}\_{\text{D}} \right) \mathbf{I}\_{-\left(1/2\right)} \left( \sqrt{\mathbf{u}} \mathbf{y} \right) - \mathbf{I}\_{-\left(1/2\right)} \left( \sqrt{\mathbf{u}} \mathbf{x}\_{\text{c}} \mathbf{y} \right) \mathbf{K}\_{-\left(1/2\right)} \left( \sqrt{\mathbf{u}} \mathbf{y} \right)}}{\left\{ \mathbf{I}\_{\left(1/2\right)} \left( \sqrt{\mathbf{u}} \mathbf{y} \right) \mathbf{K}\_{-\left(1/2\right)} \left( \sqrt{\mathbf{u}} \mathbf{x}\_{\text{c}} \mathbf{y} \right) + \mathbf{I}\_{-\left(1/2\right)} \left( \sqrt{\mathbf{u}} \mathbf{x}\_{\text{c}} \mathbf{y} \right) \mathbf{K}\_{\left(1/2\right)} \left( \sqrt{\mathbf{u}} \mathbf{y} \right) \right\}}} \left| \overleftarrow{\mathbf{q}} \mathbf{y} \right|\_{\mathbf{x}\_{\text{D}}=1'} \tag{A11a} \right| \mathbf{x}\_{\text{c}} = 1' \mathbf{y}$$

$$\left. \frac{\partial \overline{\Psi\_{\rm SmD}}}{\partial \mathbf{x}\_{\rm D}} \right|\_{\mathbf{x}\_{\rm D} = 1} = \sqrt{a\_{5}} \left( \frac{\left| \mathcal{K}\_{-\left(1/2\right)} \left( \sqrt{a\_{5}} \mathbf{x}\_{\rm d} \right) \mathcal{I}\_{-\left(1/2\right)} \left( \sqrt{a\_{5}} \right) - \mathcal{I}\_{-\left(1/2\right)} \left( \sqrt{a\_{5}} \mathbf{x}\_{\rm d} \right) \mathcal{K}\_{-\left(1/2\right)} \left( \sqrt{a\_{5}} \right) \right|}{\left| \mathcal{I}\_{\left(1/2\right)} \left( \sqrt{a\_{5}} \right) \mathcal{K}\_{-\left(1/2\right)} \left( \sqrt{a\_{5}} \mathbf{x}\_{\rm d} \right) + \mathcal{I}\_{-\left(1/2\right)} \left( \sqrt{a\_{5}} \mathbf{x}\_{\rm d} \right) \mathcal{K}\_{\left(1/2\right)} \left( \sqrt{a\_{5}} \right) \right|} \right) \overline{\Psi\_{\rm AmD}} \Big|\_{\mathbf{x}\_{\rm D} = 1'} \tag{A11b}$$

where *a*3 = *a*5 = *kwkm* ω*ms* + λ*k* − λ2*k* <sup>ω</sup>*ks*+λ*k* . According to Equation (3), the dimensionless diffusivity equations in Laplace domain for region IV can be given by

$$\begin{cases} \begin{aligned} \overline{\varphi\_{4kD}} &= \frac{\alpha k k x\_f^2}{\sigma\_k k\_k x\_f^2 + k\_w \kappa \omega\_k} \overline{\varphi\_{4mD}} \\\ \omega\_m s \overline{\varphi\_{4mD}} &= \frac{k\_m}{k\_w} \frac{1}{y\_{wD}^2} \frac{\partial^2 \overline{\varphi\_{4mD}}}{\partial y\_D^2} + \sigma\_k \frac{k\_l}{k\_w} x\_f^2 (\overline{\varphi\_{4kD}} - \overline{\varphi\_{4mD}}) + \frac{k\_l}{k\_w} \frac{\overline{\varphi\_{4mD}}}{\partial x\_D} \end{aligned} &\text{A12} \end{cases} \tag{A12}$$

Solving Equation (A12) with boundary conditions, we can obtain

$$
\begin{bmatrix}
\frac{\partial \overline{\Psi\_{\mathsf{U}\mathsf{U}\mathsf{U}\mathsf{D}}}}{\partial \mathsf{y}\_{\mathsf{D}}}\Big|\_{\substack{\mathfrak{Y}\mathsf{D}-\mathfrak{Y}\mathsf{U}\mathsf{D}}} = \sqrt{\mathsf{q}\_{4}} \begin{vmatrix}
\left|\mathcal{K}\_{-\left(1/2\right)}\left(\sqrt{\mathsf{q}\_{4}}y\_{\mathsf{U}\mathsf{D}}\right)I\_{-\left(1/2\right)}\left(\sqrt{\mathsf{q}\_{4}}y\_{\mathsf{D}}\right)I\_{-\left(1/2\right)}\left(\sqrt{\mathsf{q}\_{4}}y\_{\mathsf{D}}\right)\right|\_{-\left(1/2\right)}\left(\sqrt{\mathsf{q}\_{4}}y\_{\mathsf{D}}\right)
\end{vmatrix}\right|\_{\substack{\mathfrak{Y}\mathsf{D}-\mathfrak{Y}\mathsf{D}}
}\Big|\_{\substack{\mathfrak{Y}\mathsf{D}-\mathfrak{Y}\mathsf{D}\mathsf{U}
\end{bmatrix}}}\tag{A13}
$$

 equations in Laplace domain

where *a*4 = *a*5 − *kwkm <sup>y</sup>*<sup>2</sup>*wD* √*<sup>a</sup>*5tan <sup>h</sup>(*ylD* − *yeD*) <sup>√</sup>*<sup>a</sup>*5. DimensionlessformsandsolutionsforSRV.

(2) Accordingto(3),(5), and(6),thedimensionlessdiffusivity

Equationsfor region SRV can be given by

⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩<sup>ω</sup>*fsy*θ−<sup>2</sup>(*df a*<sup>−</sup><sup>2</sup>) *D* ϕ*f D* = 1*<sup>y</sup>*<sup>2</sup>*wD* <sup>∂</sup><sup>2</sup>ϕ*f D* <sup>∂</sup>*y*<sup>2</sup>*D* + 1*<sup>y</sup>*<sup>2</sup>*wD* 3(*df a*<sup>−</sup><sup>2</sup>)+(*df s*<sup>−</sup><sup>2</sup>)−<sup>θ</sup> *yD* ∂ϕ*f D* ∂*yD* + *y*θ−<sup>2</sup>(*df a*<sup>−</sup><sup>2</sup>) *D kmkw* ∂ϕ<sup>3</sup>*m* ∂*xD* ++++*xD*=<sup>1</sup> <sup>+</sup>*y*θ−<sup>2</sup>(*df a*<sup>−</sup><sup>2</sup>) *D* <sup>λ</sup>*m*<sup>ϕ</sup>2*mD* − ϕ*f D* ϕ2*kD* = λ*k* λ*k*+*s*ω*k* ϕ2*mD* <sup>ω</sup>*m<sup>s</sup>*ϕ2*mD* = <sup>λ</sup>*k*(ϕ2*kD* − ϕ2*mD*) − <sup>λ</sup>*m*<sup>ϕ</sup>2*mD* − ϕ*f D* . (A14)

Solving Equation (A14) with boundary conditions, we can obtain

∂ϕ*f D* ∂*yD* ++++*yD*=*ywD* = √*<sup>a</sup>*2 *y* <sup>1</sup>+*m*1−2*n*<sup>1</sup> 2 *wD* ⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ *I* −1−*m*<sup>1</sup> <sup>2</sup>+*m*1−*n*<sup>1</sup> ⎛⎜⎜⎜⎜⎝ 2 √*<sup>a</sup>*2 <sup>2</sup>+*m*1−*n*<sup>1</sup> *y* <sup>2</sup>+*m*1−2*n*<sup>1</sup> 2 *wD* ⎞⎟⎟⎟⎟⎠− *h*1*h*2 *K* −1−*m*<sup>1</sup> <sup>2</sup>+*m*1−*n*<sup>1</sup> ⎛⎜⎜⎜⎜⎝ 2 √*<sup>a</sup>*2 <sup>2</sup>+*m*1−*n*<sup>1</sup> *y* <sup>2</sup>+*m*1−2*n*<sup>1</sup> 2 *wD* ⎞⎟⎟⎟⎟⎠ *y* <sup>1</sup>−*n*<sup>1</sup> 2 *wD* ⎛⎜⎜⎜⎜⎝*I* <sup>1</sup>−*n*<sup>1</sup> <sup>2</sup>+*m*1−*n*<sup>1</sup> ⎛⎜⎜⎜⎜⎝ 2 √*<sup>a</sup>*2 <sup>2</sup>+*m*1−*n*<sup>1</sup> *y* <sup>2</sup>+*m*1−2*n*<sup>1</sup> 2 *wD* ⎞⎟⎟⎟⎟⎠+ *h*1*h*2 *K* <sup>1</sup>−*n*<sup>1</sup> <sup>2</sup>+*m*1−*n*<sup>1</sup> ⎛⎜⎜⎜⎜⎝ 2 √*<sup>a</sup>*2 <sup>2</sup>+*m*1−*n*<sup>1</sup> *y* <sup>2</sup>+*m*1−2*n*<sup>1</sup> 2 *wD* ⎞⎟⎟⎟⎟⎠⎞⎟⎟⎟⎟⎠ ⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠<sup>ϕ</sup>*FD*+++*yD*=*ywD* , (A15) where *a*2 = *<sup>y</sup>*<sup>2</sup>*wD* ⎧⎪⎪⎪⎨⎪⎪⎪⎩ω*fs* − *kmkw* √*a*3⎛⎜⎜⎜⎜⎜⎜⎜⎜⎝ 2*K*−( 12 )( <sup>√</sup>*<sup>a</sup>*3*xeD*)*<sup>I</sup>*−( 12 )( <sup>√</sup>*<sup>a</sup>*3)−*I*−( 12 )( <sup>√</sup>*<sup>a</sup>*3*xeD*)*<sup>K</sup>*−( 12 )( <sup>√</sup>*<sup>a</sup>*3)3 2*I*( 12 )( <sup>√</sup>*<sup>a</sup>*3)*<sup>K</sup>*−( 12 )( <sup>√</sup>*<sup>a</sup>*3*xeD*)+*<sup>I</sup>*−( 12 )( <sup>√</sup>*<sup>a</sup>*3*xeD*)*<sup>K</sup>*( 12 )( <sup>√</sup>*<sup>a</sup>*3)3 ⎞⎟⎟⎟⎟⎟⎟⎟⎟⎠+λ*m* <sup>ω</sup>*ms*+ <sup>λ</sup>*k*<sup>ω</sup>*k<sup>s</sup>* <sup>λ</sup>*k*+<sup>ω</sup>*k<sup>s</sup>* <sup>ω</sup>*ms*+ <sup>λ</sup>*k*<sup>ω</sup>*k<sup>s</sup>* <sup>λ</sup>*k*+<sup>ω</sup>*k<sup>s</sup>* +λ*m* ⎫⎪⎪⎪⎬⎪⎪⎪⎭, *m*1 = *df a* + *df s* − 4, *n*1 = 3*df a* − 2 + *df s* − 2 − θ, *h*1 = √*<sup>a</sup>*2 *b*2 *y <sup>m</sup>*1−*n*1 2 *f D I* −1−*m*<sup>1</sup> <sup>2</sup>+*m*1−*n*<sup>1</sup> ! 2 √*<sup>a</sup>*2 2+*m*1−*n*<sup>1</sup> *y*<sup>2</sup>+*m*1−2*n*<sup>1</sup> 2 *f D* " − *I* <sup>1</sup>−*n*<sup>1</sup> <sup>2</sup>+*m*1−*n*<sup>1</sup> ! 2 √*<sup>a</sup>*2 2+*m*1−*n*<sup>1</sup> *y*<sup>2</sup>+*m*1−2*n*<sup>1</sup> 2 *f D* ", *h*2 = *K*<sup>1</sup>−*n*<sup>1</sup> <sup>2</sup>+*m*1−*n*<sup>1</sup> ! 2 √*<sup>a</sup>*2 2+*m*1−*n*<sup>1</sup> *y*<sup>2</sup>+*m*1−2*n*<sup>1</sup> 2 *f D* " + √*<sup>a</sup>*2 *b*2 *y <sup>m</sup>*1−*n*1 2 *f D K* −1−*m*<sup>1</sup> <sup>2</sup>+*m*1−*n*<sup>1</sup> ! 2 √*<sup>a</sup>*2 2+*m*1−*n*<sup>1</sup> *y*<sup>2</sup>+*m*1−2*n*<sup>1</sup> 2 *f D* ", *b*2 = √*a*5! -*K*−(1/2)( √*<sup>a</sup>*5*xeD*)*<sup>I</sup>*−(1/2)( √*<sup>a</sup>*5)−*I*−(1/2)( √*<sup>a</sup>*5*xeD*)*<sup>K</sup>*−(1/2)( <sup>√</sup>*<sup>a</sup>*5) -*<sup>I</sup>*(1/2)( √*<sup>a</sup>*5)*<sup>K</sup>*−(1/2)( √*<sup>a</sup>*5*xeD*)+*<sup>I</sup>*−(1/2)( √*<sup>a</sup>*5*xeD*)*<sup>K</sup>*(1/2)( <sup>√</sup>*<sup>a</sup>*5) ".

(3) Dimensionless forms and solutions for HFs.

According to Equation (14), the dimensionless diffusivity equations in Laplace domain for HFs can be given by

$$
\omega\_F \mathbf{s} \overline{\mathbf{q}\_{FD}} = \frac{k\_F}{k\_W} \frac{\partial^2 \overline{\mathbf{q}\_{FD}}}{\partial \mathbf{x}\_D^2} + \frac{\mathbf{x}\_f^2}{y\_w^2} \left. \frac{\partial \overline{\mathbf{q}\_{FD}}}{\partial y\_D} \right|\_{\mathbf{y} = y\_{wD}} \text{ \tag{A16}}
$$

Solving Equation (A16) with boundary conditions, we can obtain the pseudo-pressure in well bottom-hole in Laplace domain as Equation (16).
