**1. Introduction**

Unconventional reservoirs, especially the tight sand and shale reservoirs, have become more and more significant with the development of the horizontal wells and hydraulic fracturing technologies. The complex fracture network surrounding the multistage fractured horizontal wells (MFHWs) called

stimulated reservoir volume (SRV), has obvious influences on the production performance [1–3]. Currently, a lot of effort has been done to model the pressure transient and rate transient of MFHWs aiming to forecast the production accurately. Through numerical methods, analytical methods, and empirical methods [4–7], production data analysis in unconventional reservoirs, such as Barnett, Bakken, and Eagle, shows that the linear flow in formation, especially in hydraulic fractures (HFs) and SRV, dominates the production of MFHWs. Meanwhile, the advantage of analytic models is their simplicity and less parameters compared with numerical and empirical models, although they cannot characterize the over-complex heterogeneity properties and fracture system. Thus, considering the main flow regimes in the life of MFHWs in unconventional gas reservoirs, various linear flow models based on linear-flow assumptions have been developed. Ozkan et al. [8] and Brown et al. [9] reported a trilinear model to study the performance of MFHWs in unconventional reservoirs, which divides the formation into three main regions. In order to describe the SRV size between HFs more reasonably, Stalgorova and Matter [10] improved the trilinear model to a five-linear model by simplifying SRV in a region with limited width. On this basis, many scholars studied the production rate or pressure transient performances of MFHWs in unconventional reservoirs [11,12]. Although the models mentioned above can simulate the main linear-flow regimes in both SRV and unstimulated reservoir volume (USRV), the heterogeneous distribution of complex induced fractures are not considered in SRV region, which is the major distinct properties between SRV and USRV as shown in Figure 1.

**Figure 1.** Simplified physical model and three typical presented linear-flow models for a multi-fractured horizontal wells (MFHW) in unconventional reservoirs.

The distribution of induced fractures system in SRV affected by HFs and the pre-existing closed natural fractures (NFs) with self-similar characterization can be described by fractal theory [13–15]. Chang and Yortsos [16] quantified the heterogeneous permeability and porosity of NFs using fractal theory. Then, Cossio et al. [17] introduced the fractal permeability/porosity relation [18] into the fractal diffusivity equation (FDE), which proved to be more fundamental than the classical diffusivity equation. Wang et al. [19] combined FDE with dual-porosity media model in the trilinear model to describe oil flow in SRV region of MFHWs in tight reservoirs, and their model provided a suitable way to study the flow behaviors in heterogeneous induced fractures systems. Sheng et al. [20] extended the fractal trilinear model to shale gas reservoirs coupling multi-porosity media (porous kerogen, inorganic matter, and induced-fracture) with multiple gas transport mechanisms (viscous flow with slippage, Knudsen diffusion, and surface diffusion). However, they did not consider the unstimulated reservoir volume (USRV) between HFs, which cannot accurately describe the microseismic results. In addition, the fracture porosity and permeability distribution affected by the distribution of induced-fracture spacing and aperture, but there are limited reports about the calculation methods for the non-uniform distribution of induced fractures spacing and aperture using fractal theory [21].

History matching to estimate the critical parameters of formation and hydraulic fracturing and rate prediction of MFHWs plays important roles in unconventional reservoirs based on the numerical and analytic models. Samandarli et al. [22] presented a semianalytic method to estimate reservoir parameters by history matching the production data of hydraulically fractured shale gas wells using a least absolute value regression. Correction for gas properties changing with time and desorbed gas at low pressure was considered in the model. According to the algebraic equations, bilinear flow regime and linear flow regime were used to match both synthetic data and field data by doing regression on matrix permeability, fracture permeability, and length. Hategan and Hawkes [23] used a simple analytical model based on the solution of the pseudo-steady state equation to analyze the production of MFHWs in shale gas reservoirs, and thought most shale gas reservoirs fit di fferent type curves by normalizing to reservoir pressure, permeability, and number of fracture stages. Considering the e ffect of pressure-dependent natural-fracture permeability on production from shale gas wells, Cho et al. [24] modified the trilinear flow model to match the production data of Haynesville and Barnett shale gas wells. The reservoir properties and the correlation coe fficients used for the pressure-dependent permeability were matched ignoring the non-uniqueness issues caused by the large number of parameters. Ogunyomi et al. [25] developed an approximate analytical solution to the double-porosity model for MFHWs in unconventional oil reservoirs. Due to the solution obtained in real-time space, the model was used to match the field data with the least-square method mainly by changing the parameters related to time, and forecast the production rates. Clarkson et al. [26] combined analytical, semi-analytical, and empirical methods for history matching and production forecasting for MFHWs in tight and shale gas reservoirs. Linear-to-boundary (LTB) model [27], composite model [9,10], and semi-analytical model (SAM) [26] were used to match the field data by fitting the parameters of SRV and unstimulated region. Chen et al. [28] proposed a comprehensive semi-analytical model for MFHWs considering the shale gas flow mechanisms. The unknowns including the well length, fracture number, conductivity, and length were determined by history matching. Yao et al. [7] developed an analytical multi-linear model based on linear flow and radial flow assumption. Permeability and length of di fferent regions, and fracture length were estimated by history matching with several field data. Therefore, the semianalytic models coupling with history matching methods can be applied for not only obtaining some critical parameters for MFHWs but also predicting the production rate.

The objective of this paper is development and calibration of a semianalytic model for shale wells with nonuniform distribution of induced fractures based on ensemble smoother with multiple data assimilation (ES-MDA). The 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 multi-media (porous kerogen, inorganic matter, and fracture system) are taken into account and the model is verified. The e ffects of the fractal dimensions and tortuosity index of induced fractures on MFHWs performances are analyzed. Then, we employ the ES-MDA method with the presented model to reduce uncertainty in the forecasting of gas production rate for MFHWs in unconventional gas reservoirs using a synthetic case for the tight gas reservoir and a real filed case for the shale gas reservoir. The matching results indicate that the proposed approach could enrich the application of the semianalytic model in the practical field. Finally, some conclusions are drawn and presented based on the results and analysis.

#### **2. Fractal Semianalytic Model for Mfhws in Unconventional Gas Reservoirs**

As shown in Figure 2, the basic workflow of a new semianalytic model construction in this study consisted of two main parts: firstly, for SRV, which dominates the productivity of MFHWs in unconventional gas reservoirs, the fractal fracture spacing distribution (FFSD) was taken into account to describe its heterogeneous properties instead of classical idealized dual-media model; secondly, based on the linear flow assumption, combining with the linear models for unstimulated reservoir volume (USRV) (original formation), a coupled fractal multi-linear flow model (FMFM) for MFHWs

was developed. The detailed derivations of the FMFM for MFHWs in unconventional gas reservoirs are presented in following sections.

**Figure 2.** Development of physical and mathematical models. (**a**) Map of micro-seismic monitoring results; (**b**) simplified physical model presented in this study; (**c**) multi-linear mathematical model presented in this study (modified from Sheng et al. [20]); (**d**) idealized homogeneous dual-media model (IDMM) for stimulated reservoir volume (SRV) (modified from Warren and Root [29]); and (**e**) fractal heterogeneous dual-media model (FDMM) for SRV.

#### *2.1. Fractal Fracture Network Permeability and Porosity in SRV*

Hydraulic fracturing creates the fracture network surround the horizontal well-bore and hydraulic fracture. Due to differences in situ stress and fracture perforation position, the distribution of induced-fracture properties including induced-fracture spacing and aperture are heterogeneous [21], which affects the distribution of permeability and porosity of fracture network directly. Therefore, the fractal dimensions for both induced-fracture spacing and aperture (*dfs* and *dfa*) are first discussed in this section. Then the fracture network permeability and porosity are proposed (details are shown in Appendix A). As shown in Figure 2c, using the fractal theory; thus, the fracture permeability can be obtained by

$$k\_f(y) = k\_{i-f}(y)n(y)\frac{A\_p}{A\_I}\frac{\text{dx}}{\text{dL}} = k\_w \left(\frac{y}{y\_w}\right)^{3\left(d\_{fs}-2\right)+\left(d\_{fs}-2\right)-0},\tag{1}$$

where *kf* is the fracture permeability, m2; μ is the fluid viscosity, Pa·s; *Ar* is the sectional-area of a REV, m2; *Ap* is the sectional-area of induced-fracture, m2; d*L* is the flow length in the induced-fracture, m; *kw* is the fracture permeability at the reference point, m2; *yw* is the position of reference point, m; *n* is the number of fracture sites per unit thickness, m<sup>−</sup>1; *dfs* is the fractal dimension of induced-fracture spacing; *dfa* is the fractal dimension of induced-fracture aperture; and θ is the tortuosity index of induced fracture [30].

Moreover, the fracture porosity is given by [21]

$$\phi\_f(y) = \phi\_{i-f} \frac{b\_f(y)}{b\_f(y) + s\_f(y)} = \phi\_w(\frac{y}{y\_w})^{d\_{fs} + d\_{fs} - 4},\tag{2}$$

where φ*f* is the fracture porosity; φ*<sup>i</sup>*−*<sup>f</sup>* is the single induced-fracture porosity; φ*w* is the fracture porosity at the reference point; *sf* is induced-fracture spacing, m; and *bf* is induced-fracture aperture, m.

Then, Equations (1) and (2) are the basic forms of fracture permeability and porosity considering the fractal distribution in SRV of MFHWs in unconventional reservoirs. When *yw* = 0.1 m and *kw* = 10 × 10−<sup>15</sup> m2, the fracture permeability distribution with different fractal dimensions of induced-fracture spacing and aperture (*dfs* and *dfa*) and tortuosity index of induced-fracture (θ) are shown in Figure 3. The fracture permeability *kf* decreases if *dfs* < 2 and increases if *dfs* > 2 with the

increase of the distance from the reference point *yw* as shown in Figure 3a. Meanwhile, the fracture permeability *kf* decreases if *dfa* < 2 and increases if *dfa* > 2 with the increase of the distance from the reference point *yw* as shown in Figure 3b. Figure 3c shows that the tortuosity index of induced-fracture θ decreases when θ < 0 and increases when θ > 0 with the increase of the distance from the reference point *yw* [31]. The reasons can be illustrated as follows: If *dfs* < 2 or *dfa* < 2, as the distance from the reference point increases, the induced-fracture number or aperture decreases and the induced-fracture spacing increases; if *dfs* = 2 or *dfa* = 2, as the distance from the reference point increases, the induced-fracture number or aperture and the induced-fracture spacing are constants; if *dfs* > 2 or *dfa* > 2, as the distance from the reference point increases, the induced-fracture number or aperture increases and the induced-fracture spacing decreases; if θ > 0, the tortuosity of pores is more complex than that when θ < 0. Therefore, for SRV, the fractal dimensions of nature fracture are usually reported from 1.3 to 1.7 [32].

**Figure 3.** Fracture permeability distribution with different fractal dimensions of induced-fracture spacing and aperture (*dfs* and *dfa*) and tortuosity index of induced-fracture (θ). (**a**) Fracture permeability distribution in SRV with different values of *dfs* when *dfa* = 2 and θ = 0; (**b**) fracture permeability distribution in SRV with different values of *dfa* when *dfs* = 2 and θ = 0; and (**c**) fracture permeability distribution in SRV with different values of θ when *dfs* = 2 and *dfa* = 2.

#### *2.2. Transient Di*ff*usivity Equations in SRV And USRV*

Based on the trilinear model, a newly coupled FMFM was developed for MFHWs in unconventional gas reservoirs as shown in Figure 2c, which followed the basic linear flow assumptions and some simplifying assumptions: (1) the reservoir formation consisted of five regions including the hydraulic fracture, the effective SRV and three USRV regions; (2) the fractal triple-media model (FTMM) was applied to describe the heterogeneous characteristics of induced fractures in SRV, which consisted of three media (fracture network, kerogen, and inorganic matter); (3) idealized dual-media model (IDMM) was used to characterize the homogeneous properties of formation in USRV, which included porous kerogen and inorganic matter; (4) porous kerogen and inorganic matter are uniform and continuous media in both SRV and USRV regions; (5) for shale gas reservoirs, viscous flow, Knudsen di ffusion, surface di ffusion, and slip factor were taken into account; (6) single-phase fluid flowed into the horizontal wellbore from reservoir only by the hydraulic fractures, and frictional loss within the horizontal well and e ffects of gravity and capillary forces were ignored in the reservoir; (7) the hydraulic fractures had identical properties, and were evenly distributed along the horizontal well; and (8) there were the no-flow boundaries parallel to the hydraulic fractures at the middle of two fractures in closed reservoir.

#### 2.2.1. Di ffusivity Equations of FTMM in SRV (Region II)

Since the induced-fracture system, porous kerogen and inorganic matter coexist in SRV, there are two types of cross-flow flux happening: one is from porous kerogen to inorganic matter, and another one is from inorganic matter to induced-fracture system, both of which can be described by Warren–Root pseudo-steady model. Thus, combining Equations (1) and (2), the di ffusivity equation of fluid flow in the induced-fracture system can be expressed as

$$\begin{array}{c} \frac{\partial \left[\phi\_{w}(y/y\_{w})^{d\_{f}+d\_{f}-4}\left(p\_{f}\Lambda l/ZRT\right)\right]}{\partial t} \ = \frac{\partial}{\partial y} \left[\frac{p\_{f}M}{ZRT\mu\_{g}}k\_{w}(\frac{y}{y\_{w}})^{3(d\_{f}-2)+(d\_{f}-2)-0}\frac{\partial p\_{f}}{\partial y}\right] \\ \quad + \frac{1}{x\_{f}}\frac{p\_{f}M}{ZRT\mu\_{g}}k\_{m}\frac{\partial p\_{m}}{\partial x}\bigg|\_{x=x\_{f}} + q\_{m-f}' \end{array} \tag{3}$$

where *pf* is the induced-fracture pressure in SRV, Pa; *Z* is the gas compressibility factor; *R* is the universal gas constant, 8.314 J/(K·mol); μ*g* is the gas viscosity, Pa·s; *km* is the apparent permeability of the inorganic matter, m2; *pm* is the inorganic matter pressure, Pa; *xf* is the hydraulic fracture half-length, m; subscript, 3 represents region III; and *q <sup>m</sup>*−*f* is the mass transfer from inorganic matter to the induced-fracture system, kg/m3·s.

For porous kerogen, if the Knudsen di ffusion, surface di ffusion, viscous flow and slippage e ffect are included, the apparent permeability can be given by (Sheng et al., 2019b)

$$k\_{k} = \frac{\phi\_{k}}{\tau\_{k}} \frac{2r\_{k}}{3} \sqrt{\frac{8ZRT}{\pi M}} \mathcal{C}\_{\mathcal{S}} \mu\_{\mathcal{S}} + \varepsilon\_{k\mathcal{S}} (1 - \phi\_{k}) D\_{\mathcal{S}} \frac{c\_{\mu\mathcal{Z}} ZRT}{\left(p\_{L} + p\_{k}\right)^{2}} \mu\_{\mathcal{S}} + \frac{\phi\_{k}}{\tau\_{k}} \frac{r\_{k}^{2}}{8} \left[1 + \sqrt{\frac{8ZRT}{M}} \frac{\mu\_{\mathcal{S}}}{p\_{k} r\_{k}} \left(\frac{2}{f} - 1\right)\right] \tag{4}$$

where *kk* is the apparent permeability of the porous kerogen, m2; φ*k* is the porous kerogen porosity; τ*k* is the porous kerogen tortuosity; *rk* is the pore size in porous kerogen, m; *Cg* is the gas compressibility, Pa−1; ε*ks* is proportion of solid kerogen volume in total interconnected matrix; *Ds* is the surface di ffusion coe fficient, m<sup>2</sup>/s; *<sup>c</sup>*μ*s* is the Langmuir volume on the porous kerogen, mol/m3; *pL* is the Langmuir pressure, Pa; *pk* is the porous kerogen pressure, Pa; and *f* is fraction of molecules striking pore wall which are di ffusely reflected.

Here, the Warren–Root pseudo-steady-state model is used to describe the cross-flow between porous kerogen and inorganic matter. Therefore, the di ffusivity equation for the porous kerogen can be obtained by

$$\frac{\partial \left[\phi\_k(p\_{2k}M/ZRT) + \varepsilon\_{ks}(1-\phi\_k)c\_s\right]}{\partial t} = \sigma\_k \frac{p\_{2k}Mk\_k}{ZRT\mu\_\mathcal{S}}(p\_{2m} - p\_{2k}),\tag{5}$$

where *<sup>c</sup>*μ*s* is the adsorbed gas volume on the porous kerogen, mol/m3, which follows Langmuir isotherm equation; σ*k* is the pseudo-steady-state shape factor for the porous kerogen, <sup>1</sup>/m2; and subscript, 2 represents region II.

Similarly, if the Warren–Root pseudo-steady-state model is also used to describe the cross-flow between inorganic matter and induced-fracture system, the di ffusivity equation for the porous kerogen can be expressed as

$$\frac{\partial(\phi\_{m}p\_{2m}M/ZRT)}{\partial t} = \sigma\_{k}\frac{p\_{2k}Mk\_{k}}{ZRT\mu\_{\mathcal{S}}}(p\_{2k} - p\_{2m}) - \sigma\_{m}\frac{p\_{2m}Mk\_{m}}{ZRT\mu\_{\mathcal{S}}}(p\_{f} - p\_{2m}),\tag{6}$$

where φ*m* is the inorganic matter porosity; σ*m* is the pseudo-steady-state shape factor for the inorganic matter, <sup>1</sup>/m2; *km* is the apparent permeability of the inorganic matter, m2, which includes the Knudsen di ffusion, viscous flow, and the slippage e ffect are included, and can be given by

$$k\_{\rm mf} = \frac{\phi\_m}{\tau\_m} \frac{2r\_m}{3} \sqrt{\frac{8ZRT}{\pi M}} C\_{\rm g} \mu\_{\rm g} + \frac{\phi\_m}{\tau\_m} \frac{r\_m^2}{8} \left[1 + \sqrt{\frac{8ZRT}{M}} \frac{\mu\_{\rm g}}{p\_m r\_m} \left(\frac{2}{f} - 1\right)\right] \tag{7}$$

where τ*m* is the inorganic matter tortuosity; *rm* is the pore size in the inorganic matter, m.

Outer boundary conditions: flux continuity between SRV (region II) and USRV (region IV)

$$\frac{k\_w}{\mu\_\mathcal{S}} \left(\frac{y}{y\_w}\right)^{3(d\_{fs}-2)+(d\_{fs}-2)-\theta} \left.\frac{p\_f M}{ZRT} \frac{\partial p\_f}{\partial y}\right|\_{y=y\_f} = \frac{k\_m}{\mu\_\mathcal{S}} \frac{p\_{4m}M}{ZRT} \left.\frac{\partial p\_{4m}}{\partial y}\right|\_{y=y\_f}.\tag{8a}$$

Inner boundary conditions: pressure continuity between SRV (region II) and HF (region I)

$$\left.p\_f\right|\_{y=y\_{\text{iv}}} = \left.p\_{HF}\right|\_{y=y\_{\text{iv}}}.\tag{8b}$$

In Equation (8a,b), *yw* is assumed as hydraulic fracture half-width, m; *yf* is the half-width of SRV in y-direction, m; and subscript, 4 represents region IV and HF represents hydraulic fracture.

#### 2.2.2. Di ffusivity Equations of FTMM in USRV (Region III, Region IV and Region V)

The fluid flow in both region III and region V can be described as linear flow in x-direction. Thus, according to Equations (5) and (6), the di ffusivity equations of porous kerogen and inorganic matter in these two regions can be expressed as, respectively,

$$\begin{split} \frac{\frac{\partial \left[\phi\_{k}(p\_{\rm{M}}M/ZRT) + M\varepsilon\_{\rm{k}}(1-\phi\_{k})\varepsilon\_{\rm{s}}\right]}{\partial t} &= \sigma\_{k}\frac{p\_{\rm{M}}Mk\_{\rm{k}}}{ZRT\mu\_{\rm{g}}}(p\_{\rm{3m}} - p\_{\rm{3k}})\\ \frac{\partial \left(\phi\_{\rm{m}}p\_{\rm{3m}}M/ZRT\right)}{\partial t} &= \frac{\partial}{\partial x}\left(\frac{p\_{\rm{3m}}Mk\_{\rm{m}}}{ZRT\mu\_{\rm{g}}}\frac{p\_{\rm{3m}}}{\partial x}\right) + \sigma\_{k}\frac{p\_{\rm{3k}}Mk\_{\rm{k}}}{ZRT\mu\_{\rm{g}}}(p\_{\rm{3k}} - p\_{\rm{3m}}) \end{split} \tag{9}$$

$$\begin{array}{c} \frac{\partial \left[\phi\_{k}\left(p\_{55}M/ZRT\right) + M\varepsilon\_{k}\left(1-\phi\_{k}\right)c\_{s}\right]}{\partial t} = \sigma\_{k}\frac{p\_{55}Mk\_{k}}{ZRT\mu\_{\mathcal{S}}}\left(p\_{55} - p\_{5k}\right) \\ \frac{\partial \left(\phi\_{m}p\_{5m}M/ZRT\right)}{\partial t} = \frac{\partial}{\partial x}\left(\frac{p\_{5m}Mk\_{m}}{ZRT\mu\_{\mathcal{S}}}\frac{\partial p\_{5m}}{\partial x}\right) + \sigma\_{k}\frac{p\_{5k}Mk\_{k}}{ZRT\mu\_{\mathcal{S}}}\left(p\_{5k} - p\_{5m}\right) \end{array}, \text{for region V,} \tag{10}$$

where subscript 5 represents region V.

⎧

⎪⎪⎪⎨

⎪⎪⎪⎩

⎧

⎪⎪⎪⎨

⎪⎪⎪⎩

Outer boundary conditions: no-flow conditions for region III and region V

$$\left.\frac{\partial p\_{3\text{in}}}{\partial \mathbf{x}}\right|\_{\mathbf{x}=\mathbf{x}\_{\mathbf{f}}} = \mathbf{0},\tag{11a}$$

$$\left.\frac{\partial p\_{5\text{nr}}}{\partial \mathbf{x}}\right|\_{\mathbf{x}=\mathbf{x}\_{\text{c}}} = 0.\tag{11b}$$

Inner boundary conditions: pressure continuity for region III—SRV (region II), and region V—region IV

$$\left.p\_{\mathfrak{M}}\right|\_{\mathfrak{x}=\mathfrak{x}\_{f}} = \left.p\_{f}\right|\_{\mathfrak{x}=\mathfrak{x}\_{f}'} \tag{11c}$$

$$\left.p\_{5m}\right|\_{\mathbf{x}=\mathbf{x}\_f} = \left.p\_{4m}\right|\_{\mathbf{x}=\mathbf{x}\_f}.\tag{11d}$$

In Equation (11a–d), *xe* is the reservoirs half-size in x-direction, m; subscript 4 represents region IV. Whereas, the fluid flow in y-direction in region IV, and the di ffusivity equation can be given by

$$\begin{cases} \frac{\partial \left[\phi\_{k}(p\_{4k}M/ZRT) + M\varepsilon\_{k}(1-\phi\_{k})\varepsilon\_{s}\right]}{\partial t} = \sigma\_{k}\frac{p\_{4k}Mk\_{\mathrm{k}}}{ZRT\mu\_{\mathrm{g}}}(p\_{4m} - p\_{4k})\\ \frac{\partial \left(\phi\_{m}p\_{4m}M/ZRT\right)}{\partial t} = \frac{\partial}{\partial y}\left(\frac{p\_{4m}Mk\_{\mathrm{w}}}{ZRT\mu\_{\mathrm{g}}}\frac{\partial p\_{4m}}{\partial y}\right) + \sigma\_{k}\frac{p\_{4k}Mk\_{\mathrm{k}}}{ZRT\mu\_{\mathrm{g}}}(p\_{4k} - p\_{4m}) + \frac{1}{\chi\_{f}}\frac{p\_{4m}M}{ZRT\mu\_{\mathrm{g}}}k\_{\mathrm{m}}\frac{\partial p\_{4m}}{\partial x}\Big|\_{x=x\_{f}} \end{cases} \tag{12}$$

Outer boundary conditions: no-flow conditions for region IV

$$\left.\frac{\partial p\_{4m}}{\partial y}\right|\_{y=y\_{\varepsilon}} = 0.\tag{13a}$$

Inner boundary conditions: pressure continuity for region IV—SRV (region II)

$$\left.p\_{4m}\right|\_{y=y\_f} = \left.p\_f\right|\_{y=y\_f}.\tag{13b}$$

In Equation (13a,b), *ye* is the half-size of HF spacing in y-direction, m.

#### 2.2.3. Di ffusivity Equations of FTMM in HF (Region I)

Linear flow occurs in the hydraulic fractures with closed tip and constant production in x-direction, and the di ffusivity equation can be expressed as

$$\frac{\partial(\phi\_{\rm F}\mathrm{p\_{F}M/ZRT})}{\partial t} = \frac{\partial}{\partial \mathbf{x}} \bigg| \frac{\mathbf{p\_{F}M/ZRT}}{\mu\_{\mathcal{S}}} \mathbf{k\_{F}} \frac{\partial p\_{\rm F}}{\partial \mathbf{x}} \bigg) + \frac{1}{y\_{\rm w}} \frac{\mathbf{p\_{F}M}}{ZRT\mu\_{\mathcal{S}}} \mathbf{k\_{w}} \bigg(\frac{y}{y\_{\rm w}}\bigg)^{3(d\_{fz}-2)+(d\_{fz}-2)-0} \frac{\partial p\_{f}}{\partial y} \bigg|\_{y=y\_{\rm w}},\tag{14}$$

where φ*F* is the HF porosity; *pF* is the HF pressure, Pa; and *kF* is the HF permeability, m2.

Outer boundary conditions: no-flow conditions for region I

$$\left.\frac{\partial p\_F}{\partial x}\right|\_{x=x\_f} = 0.\tag{15a}$$

Inner boundary conditions: constant production without wellbore storage and skin e ffect

$$\left.\frac{\partial p\_F}{\partial x}\right|\_{x=0} = -\frac{q\_f ZRT\mu\_g}{2Mk\_F h y\_w p\_F \Big|\_{x=0}},\tag{15b}$$

where *qf* is constant well production from HF, kg/s.

For solving the FMFM as shown in Appendix B, we can obtain the pseudo-pressure in well bottom-hole in Laplace domain as

$$\overline{\boldsymbol{\varrho}\_{fD}} = \overline{\boldsymbol{\varrho}\_{FD}}\Big|\_{\mathbf{x}\_{D}=0} = \frac{\boldsymbol{\pi}k\_{w}\mathbf{x}\_{f}}{Zk\_{F}\mathbf{y}\_{w}\mathbf{s}\cdot\tan\left(\sqrt{\frac{k\_{f}}{k\_{F}}\left(\boldsymbol{\alpha}\_{F}\mathbf{s}-\frac{\mathbf{x}\_{f}}{\boldsymbol{\mathcal{Y}}\_{w}}\mathbf{a}\_{1}\right)}\right)\sqrt{\frac{k\_{f}}{k\_{F}}\left(\boldsymbol{\alpha}\_{F\mathbf{s}}-\frac{\mathbf{x}\_{f}}{\boldsymbol{\mathcal{Y}}\_{w}}\mathbf{a}\_{1}\right)}},\tag{16}$$

$$\begin{array}{l} \text{where}, \ a\_{1} = \sqrt{a\_{2}} \underline{y}\_{\text{w}\boldsymbol{\mathcal{D}}}^{\frac{1+m\_{1}-2\boldsymbol{m}\_{1}}{2}} \left( \frac{\boldsymbol{I}}{\frac{2+m\_{1}-\boldsymbol{m}\_{1}}{2}} \boldsymbol{Y}\_{\boldsymbol{w}\boldsymbol{\mathcal{D}}}^{\frac{2+m\_{1}-2\boldsymbol{m}\_{1}}{2}} \boldsymbol{Y}\_{\boldsymbol{w}\boldsymbol{\mathcal{D}}}^{-\frac{1}{2}}}{\frac{1-\boldsymbol{m}\_{1}}{2} \left( \boldsymbol{I} \right)} \frac{\boldsymbol{I}}{\frac{2+m\_{1}-\boldsymbol{m}\_{1}}{2}} \boldsymbol{Y}\_{\boldsymbol{w}\boldsymbol{\mathcal{D}}}^{-\frac{1}{2}} \boldsymbol{Y}\_{\boldsymbol{w}\boldsymbol{\mathcal{D}}}^{\frac{2+m\_{1}-2\boldsymbol{m}\_{1}}{2}} \boldsymbol{Y}\_{\boldsymbol{w}\boldsymbol{\mathcal{D}}}^{-\frac{1}{2}}} \right) \\ \text{variables are shown in Appendix B.} \end{array}$$

variables are shown in Appendix B.

Based on Equation (16), if the wellbore storage and skin factor are not considered, the dimensionless flow rate can be derived as [33,34]

$$\overline{q\_D} = \frac{Zk\_F y\_{w^\*} \tan\left(\sqrt{\frac{k\_f}{k\_F} \left(\omega\_{FS} - \frac{x\_f}{y\_w} a\_1\right)}\right) \sqrt{\frac{k\_f}{k\_F} \left(\omega\_{FS} - \frac{x\_f}{y\_w} a\_1\right)}}{\pi s k\_w x\_f} \tag{17}$$

If we hope to obtain the solution in the real-time domain, here, the Stehfest algorithm is applied, and Equation (17) becomes

$$q\_D = \frac{\ln 2}{t\_D} \sum\_{i=1}^{N} (-1)^{\frac{N}{2} + 1} \times \left\{ \begin{array}{c} \left| \min(\frac{i \frac{N}{2}}{2}) \right| \frac{k \frac{\Psi}{2} (2k)!}{\left( \frac{N}{2} - k \right)! \operatorname{Id} \left( (k-1)! \left( \frac{i}{k} \right)! (2k-i)!} \right) \\\\ \times \frac{\operatorname{Zky}\_{\mathcal{Y}} \cdot \tan \left( \sqrt{\frac{k\_f}{k\_F} \left( \operatorname{wyr} - \frac{x\_f}{y\_w} a \right)} \right) \sqrt{\frac{k\_f}{k\_F} \left( \operatorname{wyr} - \frac{x\_f}{y\_w} a \right)}}{\pi k\_0 x\_f} \times \left( \frac{\ln 2}{t\_D} i \right) \end{array} \right\}. \tag{18}$$

#### **3. Model Verification and Discussion**
