*2.2. Variables Definitions*

To introduce the mathematical model conveniently and concisely, the dimensionless variables are defined, respectively, by

$$p\_D = \frac{2\pi k\_m h (p\_i - p)}{q \mu B}, \; \xi\_D = \frac{\xi}{L\_{ref}}, \; q\_D = \frac{q}{q\_{ref}}, \; \mathcal{C}\_{fD} = \frac{k\_f w\_f}{k\_m L\_{ref}} \tag{1}$$

When the flow in the fracture satisfies non-Darcy effect owing to high velocity, the corresponding dimensionless variables are defined as

$$(q\_{\rm DND})\_f = \frac{k\_f \beta}{\mu w\_f h} q\_{ref'} \text{ and } q\_{\rm cD} = \frac{q\_c}{q\_{ref}} \tag{2}$$

When the fracture conductivity behaves pressure sensitivity, the corresponding dimensionless variables are defined as

$$
\gamma\_{fD} = \frac{q\mu B}{2\pi k\_m h} \gamma\_f \tag{3}
$$

where *pm* is the pressure in the reservoir, Pa; *pf* is the pressure in the fracture, Pa; *pi* is the pressure on the elliptical boundary, Pa; *pex* is the pressure on the width of inner SRV region, Pa; *pey* is the pressure

on the length of inner SRV region, Pa; (*xofm*, *yofm*) is the location of *m*-th fracture tip, m; *xe* is the width of inner SRV region, m; *y*e is the length of inner SRV region, m; *xR* is the width of the whole drainage region, m3; *h* is the reservoir thickness, m; *wf* is the fracture width, m; *Lf* is the fracture length, m; *rw* is the radius of wellbore, m; *Lref* is the reference length, m; *qfm* is the function of inflow distribution along the *m*-th fracture, m<sup>2</sup>/s; *qwm* is the production rate, m<sup>3</sup>/s; *qc* is the influx rate, m<sup>3</sup>/s; μ is the viscosity, Pa s; *B* is the volume factor; γ*f* is the permeability modulus, Pa−1; and β is the Beta factor, m<sup>−</sup>1. Here, the subscript D represents dimensionless variables. + − − − = ∂ ∂ + ∂ ∂ ¦³ =

## **3. Mathematical Model**

#### *3.1. Fluid Flow in the Inner SRV Region*

The SRV region is regarded as a rectangular reservoir with constant-pressure boundary containing a set of hydraulic fractures, as shown in Figure 3. The diffusion equation governing fluid flow is given as the following dimensionless form:

$$\frac{\partial^2 p\_{\rm mD}}{\partial \mathbf{x}\_D^2} + \frac{\partial^2 p\_{\rm mD}}{\partial y\_D^2} + 2\pi \sum\_{m=1}^{N\_f} \int\_0^{\mathcal{L}\_{f\rm Dm}} q\_{f\rm Dm}(\mathbf{u}\_D) \delta(\mathbf{x}\_D - \mathbf{x}\_{\rm of\rm Dm} - \mathbf{u}\_D) \delta(y\_D - y\_{\rm of\rm Dm}) d\mathbf{u}\_D = 0 \tag{4}$$

**Figure 3.** Schematic of MFHW in the inner region with constant-pressure outer boundary.

The outer boundary conditions are satisfied by

$$\begin{cases} \begin{aligned} p\_{\textit{mD}}(\mathbf{x}\_{\textit{D}} = \mathbf{x}\_{\textit{eD}}, y\_{\textit{D}}) &= p\_{\textit{mD}}(\mathbf{x}\_{\textit{D}} = \mathbf{0}, y\_{\textit{D}}) = p\_{\textit{eD}\mathbf{x}} \\ p\_{\textit{mD}}(\mathbf{x}\_{\textit{D}}, y\_{\textit{D}} = \mathbf{0}) &= p\_{\textit{mD}}(\mathbf{x}\_{\textit{D}}, y\_{\textit{D}} = y\_{\textit{eD}}) = p\_{\textit{eD}y} \end{aligned} \tag{5}$$

where there are *Nf* fractures in the SRV region, the dimensionless length of *m*-th fracture is *LfDm*, and the dimensionless coordinate of *m*-th fracture tip is (*xofDm*, *yofDm*). The pressure on the outer boundary parallel to the fracture is *peDx*, and the pressure on the outer boundary perpendicular to the fracture is *peDy*.

Using Fourier transformation and inverse transformation (Appendix A), the dimensionless pressure distribution caused by *Nf* fracture is given by

$$p\_{\rm nD}(\mathbf{x}\_{\rm D}, y\_{\rm D}) = p\_{\rm eD\pi} lB(\mathbf{x}\_{\rm D}, y\_{\rm D}) + p\_{\rm eDy} lD(\mathbf{x}\_{\rm D}, y\_{\rm D}) + 2\sum\_{m=1}^{N\_f} \int\_{0}^{L\_{f\rm Dm}} q\_{f\rm Dm}(\mathbf{u}\_{\rm D}) \delta p\_{\rm nD}(\mathbf{x}\_{\rm D}, y\_{\rm D}, u\_{\rm D}) du\_{\rm D} \tag{6}$$

#### *3.2. Fluid Flow in the Outer Region*

The fluids stored in the outer region are depleted and then flow into the inner SRV region through the interface (denoted by the yellow block in Figure 4). The streamlines satisfy the hyperbolic shape (orange solid line in Figure 4), which are the counterpart of elliptical-shaped isopotential line (black dashed line). Once the fluids enter the inner region, the streamlines are in a linear shape around the interface. Therefore, the coordinate systems are divided into two subsets: 1D elliptical coordinate in the outer region (blue region), and 1D linear coordinate (yellow region).

**Figure 4.** Elliptical coordinate used in the outer elliptical region.

The outer boundary satisfies the constant pressure of *peD*.

In the outer elliptic region, as shown in Figure 4, there are two sets of coordinates: the Cartesian coordinates (*x*',*y*') and the elliptic coordinates (ξ,η). The relationship is given by

$$\mathbf{x}' = (\mathbf{x}\_{\mathfrak{e}}/2)\cosh\xi\cos\eta,\ y' = (\mathbf{x}\_{\mathfrak{e}}/2)\sinh\xi\sin\eta \tag{7}$$

We used the elliptic coordinates to simulate the flowing process. The outer and inner boundaries of the elliptic region, represented by ξ in the elliptic coordinates, could be expressed using spatial variables in the Cartesian coordinates:

$$\begin{cases} \, \, \mathbb{X}\_{\varepsilon} = \mathbb{X}(\mathbf{x}' = \mathbf{x} \mathbf{\hat{x}} / 2, \mathbf{y}' = \mathbf{0}) = \ln[(\mathbf{x} \mathbf{\hat{x}} + \sqrt{\mathbf{x}\_{R}^{2} - \mathbf{x}\_{\varepsilon}^{2}}) / \mathbf{x}\_{\varepsilon}] \\\ \, \, \mathbb{X}\_{w} = \xi(\mathbf{x}' = \mathbf{0}, \mathbf{y}' = \mathbf{0}) = \ln 1 = 0 \end{cases} \tag{8}$$

The corresponding pressure gradient is written as the following dimensionless expression,

$$\left(\frac{\partial p\_{mD}}{\partial \xi\_D}\right)\_{\xi\_D = 0} = \frac{p\_i - p\_{cy}}{\xi\_\varepsilon - \xi\_w} = -\frac{p\_{cyD}}{\ln[(\chi\_{RD} + \sqrt{\chi\_{RD}^2 - \chi\_{\varepsilon D}^2})/\chi\_{\varepsilon D}]} \tag{9}$$

Note that the flux inflow in Equation (9) is continuous on the interface between the outer elliptical region and inner SRV region. Meanwhile, the flowing is modeled using linear coordinates. As a result, the continuous condition is written in the form of the average integral, which is given by

$$\underbrace{\frac{1}{\mathbf{x}\_{\mathrm{cD}}} \int\_{0}^{\pi} \left. \frac{\partial p\_{mD}}{\partial \dot{\xi}\_{D}} \right|\_{\xi\_{D} = 0} d\eta\_{D}}\_{\text{elliptic flow}} = \underbrace{\frac{1}{\mathbf{x}\_{\mathrm{cD}}} \int\_{0}^{\mathbf{x}\_{\mathrm{cD}}} \left. \frac{\partial p\_{mD}}{\partial y\_{D}} \right|\_{y\_{D} = y\_{\mathrm{cD}}} d\mathbf{x}\_{\mathrm{D}}}\_{\text{linear flow}} \tag{10}$$

It could be expressed as the function of *peDx* and *peDy* (Appendix B), which is written as follows:

$$A\_2 p\_{\varepsilon Dx} + B\_2 p\_{\varepsilon yD} = \sum\_{m=1}^{N\_f} q\_{\varepsilon Dm} \text{C}\_{2,m} \tag{11}$$

Likewise, the continuity condition of flux inflow between outer linear region and inner SRV region contributes to

$$A\_1 p\_{eDx} + B\_1 p\_{cyD} = \sum\_{m=1}^{N\_f} q\_{wDm} \mathbb{C}\_{1,m} \tag{12}$$

We can obtain the dimensionless pressure on the interface by combing Equations (11) and (12),

⎧⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎩ *peDx* = *N*&*f <sup>m</sup>*=1 *qwDm <sup>B</sup>*2*C*1,*m*<sup>−</sup>*B*1*C*2,*m A*1*B*2−*A*2*B*<sup>1</sup> *peDy* = *N*&*f <sup>m</sup>*=1 *qwDm <sup>A</sup>*1*C*2,*m*<sup>−</sup>*A*2*C*1,*m A*1*B*2−*A*2*B*<sup>1</sup> (13)

where the relationship between influx distribution and production rate of *m*-th fracture is given as.

Substituting Equation (13) into Equation (6) could yield the pressure distribution caused by *Nf* fractures, which is the function with regard to influx distribution along the fracture.

#### *3.3. Fluid Flow in the Fracture*

For hydraulic fracture, the flow pattern is widely assumed to be 1D linear and incompressible [36]. As shown in Figure 5, the governing equation of fluid flow on the *m*-th fracture is described as the following dimensionless form:

$$\frac{\partial^2 p\_{fDm}}{\partial \mathbf{x}\_{Dm}^2} - \frac{2\pi}{\mathbb{C}\_{fDm}} q\_{fDm}(\mathbf{x}\_{Dm}) + \frac{2\pi}{\mathbb{C}\_{fDm}} [q\_{wfDm} \delta(\mathbf{x}\_{Dm} - \mathbf{x}\_{wDm})] = 0 \tag{14}$$

**Figure 5.** Schematic of the flow pattern in the transverse fracture.

The corresponding boundary conditions are described as

$$\left.\frac{\partial p\_{fDm}}{\partial \mathbf{x}\_{Dm}}\right|\_{\mathbf{x}\_{Dm}=0} = \left.\frac{\partial p\_{fD}}{\partial \mathbf{x}\_{Dm}}\right|\_{\mathbf{x}\_{Dm}=L\_{fDm}} = 0\tag{15}$$

Integrating Equation (14) from 0 to *xDm* with respect to *xDm* based on Equation (15) twice, the result is given by

$$p\_{f\text{Dm}}(\mathbf{x}\_{\text{Dm}}) = p\_{f\text{Dm}}(0) + \frac{2\pi}{\mathbb{C}\_{f\text{Dm}}} \int\_{0}^{\mathbf{x}\_{\text{Dm}}} d\zeta \int\_{0}^{\zeta} q\_{f\text{Dm}}(\xi) d\xi - \frac{2\pi}{\mathbb{C}\_{f\text{Dm}}} q\_{\text{tr}\text{Dm}} H(\mathbf{x}\_{\text{Dm}}, \mathbf{x}\_{\text{wDm}}) \tag{16}$$

where <sup>δ</sup>(*<sup>x</sup>*, *<sup>x</sup>*0) = 2 <sup>∞</sup>, *x* = *x*0 0, *x* - *x*0 is the Dirac function, *<sup>H</sup>*(*<sup>x</sup>*, *<sup>x</sup>*0) = 2 1, *x* ≥ *x*0 0, *x* < *x*0 is the Heaviside function, and *<sup>G</sup>*(*<sup>x</sup>*, *<sup>x</sup>*0) = 2 *x* − *x*0, *x* ≥ *x*0 0, *x* < *x*0 is the integral of Heaviside function.

As shown in Figure 5, the fluid flow is regarded as a 1D linear pattern in the region far from wellbore, and becomes a 1D radial pattern around the wellbore. It is di fferent from the vertical fracture, where the vertical fracture is in lateral contact with the wellbore and only linear flow is considered. The additional pressure loss owing to the convergence of fluids into the horizontal wellbore should be considered. The e ffect of flow convergence is regarded as a skin factor, namely convergence flow skin:

$$S\_{\rm c,m} = \frac{h\_D}{L\_{fDm}} \frac{1}{\mathcal{C}\_{fDm}} \left[ \ln \left( \frac{h\_D}{2r\_{wD}} \right) - \frac{\pi}{2} \right] \tag{17}$$

Equation (16) is rewritten as the final expression by

$$p\_{f\text{Dm}}(\mathbf{x}\_{\text{Dm}}) = p\_{f\text{Dm}}(0) + q\_{\text{pDm}}\mathbf{S}\_{\text{c,W}} + \frac{2\pi}{\mathbf{C}\_{f\text{Dm}}} \int\_{0}^{\mathbf{x}\_{\text{Dm}}} d\zeta \int\_{0}^{\zeta} q\_{f\text{Dm}}(\xi) d\xi - \frac{2\pi}{\mathbf{C}\_{f\text{Dm}}} q\_{\text{pDm}} \mathbf{G}(\mathbf{x}\_{\text{Dm}}, \mathbf{x}\_{\text{uDm}}) \tag{18}$$

#### 3.3.1. Model of a Conductivity Fracture with Non-Darcy Flow

Non-Darcy flow may occur in the fracture when the inertial forces may no longer be neglected compared with viscous forces. It is very common near production wells where local velocities can be very high. Extra pressure drop is required to overcome the inertial forces.

To account for the e ffect of extra pressure loss caused by high-velocity flow, which is proportional to the square of the velocity, the Forchheimer flow equation is presented [27]:

$$
\nabla p = \frac{\mu}{k} v + \beta v |v| \tag{19}
$$

The apparent fracture permeability in Equation (19) is defined as

$$k\_{f, \mu p} = k\_f \frac{1}{1 + k\_f \beta(\left|v\right| / \mu)}\tag{20}$$

The corresponding apparent permeability is rewritten by

$$k\_{f,app} = k\_f \frac{1}{1 + (q\_{DND})\_f \left| q\_{cD}(\mathbf{x}\_D) \right|} \tag{21}$$

Substituting Equation (21) into Equation (18) yields the following form:

$$\mathcal{L}\_{f\text{D,app}}(q\_{c\text{Dm}})\frac{\partial p\_{f\text{Dm}}}{\partial \mathbf{x}\_{\text{Dm}}} + 2\pi q\_{c\text{Dm}} + 2\pi q\_{\text{wDm}}H(\mathbf{x}\_{\text{Dm}}, \mathbf{x}\_{\text{wDm}}) = 0 \tag{22}$$

where the apparent conductivity is written as *Cf <sup>D</sup>*,*app*(*qcDm*) = *Cf Dm*(*xD*) <sup>1</sup>+(*qDND*)*f* · + + + +*qcDm* + + + + , and the flux on the cross

section is given as *qcDm* = −*Cf Dm* 2π <sup>∂</sup>*pf Dm* ∂*xDm* . The corresponding integral form is given by

$$q\_{\rm cDm} = q\_{\rm wDm} H(\mathbf{x}\_{\rm Dm}, \mathbf{x}\_{\rm wDm}) - \int\_0^{\rm x\_{\rm Dm}} q\_{fDm}(\mathbf{x}') d\mathbf{x}' \tag{23}$$

Note that Equation (22) presents a general model for a uniform-conductivity under the non-Darcy-flow condition.

#### 3.3.2. Model of a Conductivity Fracture with Pressure Sensitivity E ffect

The variation of permeability caused by a stress change can be expressed as a function of pore pressure. For hydraulic fracture, a general model presented by Zhang et al. [31] is introduced to describe the relationship between fracture permeability and pore pressure next:

$$k\_f(p\_f) = k\_{\min} + [k\_f(p\_i) - k\_{f\min}] \cdot \exp[-\gamma\_f \cdot (p\_i - p\_f)]\tag{24}$$

According to the definitions of the dimensionless parameters, the pressure-sensitive conductivity of a fracture is

$$\frac{\mathbb{C}\_{fD}(p\_{fD})}{\mathbb{C}\_{fDi}} = \left(1 - \frac{\mathbb{C}\_{fD\text{min}}}{\mathbb{C}\_{fDi}}\right) \cdot \exp\left[-\gamma\_{fD} \cdot p\_{fD}\right] + \frac{\mathbb{C}\_{fD\text{min}}}{\mathbb{C}\_{fDi}}\tag{25}$$

Substituting Equation (25) into Equation (14) yields the following equation governing flow in the m-th fracture flowing:

$$\frac{\partial}{\partial \mathbf{x}\_{Dm}} \Big( \mathbf{C}\_{fDm}(p\_{fDm}) \frac{\partial p\_{fDm}}{\partial \mathbf{x}\_{Dm}} \Big) - 2\pi q\_{fDm}(\mathbf{x}\_{Dm}) + 2\pi q\_{wDm} \delta(\mathbf{x}\_{Dm}, \mathbf{x}\_{wDm}) = 0 \tag{26}$$

Note that Equation (26) presents a general model for a uniform-conductivity with the pressure sensitivity e ffect.

#### **4. Semi-Analytical Solution for Coupled Model**

## *4.1. Dimension Transformation*

Compared with the equation of uniform conductivity fracture described by Equations (14), (22) and (26) are subject to nonlinear equations. In the case of considering non-Darcy flow, the fracture conductivity is a function of flowing rate; in the case of considering pressure-sensitivity e ffect, the fracture conductivity is a function of pressure. In essence, both flow rate and pressure chang with the spatial variable. Here, we attempt to solve the problem in the same way as the linear treatment in Equation (14) by a dimension transformation, which is defined as

$$\begin{cases} \ \mathcal{L}\_{Dm}(\mathbf{x}\_{Dm}) = \mathcal{C}\_{fDm} \cdot \int\_{0}^{\mathbf{x}\_{Dm}} \frac{1}{\mathcal{C}\_{fDm}(\mathbf{x}\_{D}')} d\mathbf{x}\_{D}' \\\ \mathcal{C}\_{fDm} = L\_{fDm} / \int\_{0}^{L\_{fDm}} \frac{1}{\mathcal{C}\_{fDm}(\mathbf{x}\_{D}')} d\mathbf{x}\_{D}' \end{cases} \tag{27}$$

Then, the fracture-flow equations (Equations (22) and (26)) can be written as

$$p\_{f\rhd\mathfrak{m}}(\xi\_{\square\mathfrak{m}}) = p\_{f\rhd\mathfrak{m}}(\xi\_{\square\mathfrak{m}} = 0) + q\_{\mathfrak{m}\triangleright\mathfrak{m}}\xi\_{\square,\mathfrak{m}} + \frac{2\pi}{\mathfrak{C}\_{f\rhd\mathfrak{m}}} \left( \int\_{0}^{\xi\_{\square\mathfrak{m}}} d\zeta \int\_{0}^{\zeta} \widetilde{q}\_{f\rhd\mathfrak{m}}(\zeta) d\zeta - q\_{\mathfrak{u}\triangleright\mathfrak{m}}G(\xi\_{\square\mathfrak{m}}, \xi\_{\mathsf{u}\triangleright\mathfrak{m}}) \right) \tag{28}$$

Equation (28) has the same form as the fracture-flow equations for a uniform-conductivity fracture with constant conductivity given by Equation (18). The di fference is that the variable ξ*D* is a function of the variable *xD* and depends on the distribution of conductivity *CfD*(*xD*).
