**6. Conclusions**

We derived a new steady-state productivity model for a multi-fractured horizontal well in a reservoir with pseudo-elliptical boundary. We summarize some important findings as follows:


**Author Contributions:** M.C. established the model and wrote the manuscript, and H.X. supervised this research and revised the manuscript. C.W. validated the model and conducted the sensitivity analysis. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Science and Technology Major Project (No. 2017ZX05019005-006) and Heilongjiang Natural Science Foundation (LH2019F004). Special thanks to the anonymous reviewers and editors for their valuable comments.

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

#### **Appendix A. Analytical Solution for the Inner SRV Region**

The appendix is an optional section that can contain details and data supplemental to the main text. For example, explanations of experimental details that would disrupt the flow of the main text, but nonetheless remain crucial to understanding and reproducing the research shown; figures of replicates for experiments of which representative data is shown in the main text can be added here if brief, or as Supplementary data. Mathematical proofs of results not central to the paper can be added as an appendix.

According to the governing equation and boundary conditions, Fourier sine transformation are used to deal with Equation (4), which is defined as

$$F\_{\mathbf{x}}[f(\mathbf{x}\_{\mathcal{D}})] = \int\_{0}^{\mathbf{x}\_{\mathcal{d}\mathcal{D}}} f(\mathbf{x}\_{\mathcal{D}}) \sin(\boldsymbol{\beta}\_{\mathcal{V}} \mathbf{x}\_{\mathcal{D}}) d\mathbf{x}\_{\mathcal{D}} = \overline{f}(\boldsymbol{\beta}\_{\mathcal{V}}) \tag{A1}$$

$$F\_y[f(y\_D)] = \int\_0^{y\_D} f(y\_D) \sin(\gamma\_v y\_D) dy\_D = \overline{f}(\gamma\_v) \tag{A2}$$

$$\beta\_{\nu} = \frac{\nu \pi}{\chi\_{\text{t}\text{D}}}, \nu = 0, 1, 2, 3 \cdots, \text{ } \gamma\_{\nu} = \frac{\nu \pi}{\mathcal{Y}^{\text{d}\text{D}}}, \text{ } \nu = 0, 1, 2, 3 \cdots \text{ }, \tag{A3}$$

Taking the Fourier transformation of Equations (4) and (5), the first term on the left hand side (LHS) of the original equation is given as

$${}\_{3}F\_{y}\left[F\_{x}\left(\frac{\partial^{2}p\_{mD}}{\partial x\_{D}^{2}}\right)\right] = \frac{\beta\_{v}p\_{cDx}[1-(-1)^{v}]}{\gamma\_{v}}[1-\cos(\gamma\_{v}y\_{cD})] - \beta\_{v}^{2}\hat{\overline{p}}\_{mD} \tag{A4}$$

Likewise, the second term on the LHS is given as

$${}\_{x}F\_{y}\left[F\_{y}\left(\frac{\partial^{2}p\_{mD}}{\partial y\_{D}^{2}}\right)\right] = \frac{\gamma\_{\nu}\eta\_{\varepsilon Dy}\left[1 - (-1)^{\nu}\right]}{\beta\_{\nu}}[1 - \cos\left(\beta\_{\nu}\mathbf{x}\_{\varepsilon D}\right)] - \gamma\_{\nu}^{2}\hat{\tilde{p}}\_{mD} \tag{A5}$$

The third term on the LHS is given as

$$\begin{cases} \int\_{0}^{L\_{fDm}} q\_{fDm}(\mathbf{u}\_{D}) \delta(\mathbf{x}\_{D} - \mathbf{x}\_{ofDm} - \mathbf{u}\_{D}) \delta(y\_{D} - y\_{ofDm}) d\mu\_{D} \\ = \int\_{0}^{L\_{fD}} q\_{fD}(\mathbf{u}\_{D}) [\sin[\beta\_{\mathrm{fl}}(\mathbf{x}\_{ofD} + \mathbf{u}\_{D})] \sin(\gamma\_{\mathrm{mf}} y\_{ofD})] d\mu\_{D} \end{cases} \tag{A6}$$

Finally, the Fourier solution for Equation (4) is given as

$$\begin{split} & \frac{\beta\_{v} p\_{vDn} [1 - (1-)^{v}][1 - \cos(\gamma\_{v} y\_{vD})]}{\gamma\_{v}} - \beta\_{v}^{2} \hat{p}\_{mD} + \frac{\gamma\_{v} p\_{vDy} [1 - (1-)^{v}][1 - \cos(\beta\_{v} x\_{D})]}{\beta\_{v}} - \gamma\_{v}^{2} \hat{p}\_{mD} \\ & + 2\pi \sum\_{m=1}^{N\_{f}} \int\_{0}^{L\_{fDm}} q\_{fDm}(u\_{D}) \Big{(}\sin[\beta\_{v} (x\_{ofDm} + u\_{D})] \sin(\gamma\_{v} y\_{ofDm}) \Big{)} du\_{D} = 0 \end{split} \tag{A7}$$

The inversion Fourier sine transformation is written as

$$\begin{split} \hat{\vec{p}}\_{mD} &= \quad p\_{\epsilon D \mathbf{x}} \frac{\beta\_{\nu} [1 - (-1)^{\nu}][1 - \cos(\gamma\_{\nu} y\_{\rm{cD}})]}{\gamma\_{\nu} (\beta\_{\nu}^{2} + \gamma\_{\nu}^{2})} + p\_{\epsilon D y} \frac{\gamma\_{\nu} [1 - (-1)^{\nu}][1 - \cos(\beta\_{\nu} x\_{\rm{cD}})]}{\beta\_{\nu} (\beta\_{\nu}^{2} + \gamma\_{\nu}^{2})} \\ &+ 2\pi \sum\_{m=1}^{N\_{f}} \int\_{0}^{L\_{fDm}} q\_{fDm}(\boldsymbol{\mu} \boldsymbol{\nu}) \Big[ \frac{\sin[\beta\_{\nu} (\mathbf{x}\_{\rm{cf}Dm} + \boldsymbol{u}\_{D})] \sin(\gamma\_{\nu} y\_{\rm{cD}} \boldsymbol{\mu})}{\beta\_{\nu}^{2} + \gamma\_{\nu}^{2}} \Big] d\boldsymbol{u}\_{D} \end{split} \tag{A8}$$

where

$$N(\boldsymbol{\beta}\_{\boldsymbol{\nu}}) = \int\_{0}^{\mathbf{x}\_{\ell D}} \sin^{2}(\boldsymbol{\beta}\_{\boldsymbol{\nu}} \mathbf{x}\_{D}) d\mathbf{x}\_{D} = \frac{\mathbf{x}\_{\ell D}}{2} \Big( 1 - \frac{\sin(\beta\_{\boldsymbol{\nu}} \mathbf{x}\_{\ell D}) \cos(\beta\_{\boldsymbol{\nu}} \mathbf{x}\_{\ell D})}{\beta\_{\boldsymbol{\nu}} \mathbf{x}\_{\ell D}} \Big) \tag{A9}$$

$$N(\boldsymbol{\gamma}\_{\boldsymbol{\gamma}}) = \int\_{0}^{\mathcal{Y}\_{\mathrm{d}D}} \sin^{2}(\boldsymbol{\gamma}\_{\boldsymbol{\gamma}} \boldsymbol{y}\_{D}) d\boldsymbol{y}\_{D} = \frac{y\_{\mathrm{\ell D}}}{2} \Big( 1 - \frac{\sin(\boldsymbol{\gamma}\_{\boldsymbol{\gamma}} \boldsymbol{y}\_{\mathrm{\ell D}}) \cos(\boldsymbol{\gamma}\_{\boldsymbol{\gamma}} \boldsymbol{y}\_{\mathrm{\ell D}})}{\boldsymbol{\gamma}\_{\boldsymbol{\gamma}} \boldsymbol{y}\_{\mathrm{\ell D}}} \Big) \tag{A10}$$

Accordingly, the original solution is expressed in forms of Fourier solution,

$$p\_{mD} = \sum\_{\nu=0}^{\infty} \frac{\sin(\beta\_{\nu} \mathbf{x}\_{D})}{N(\beta\_{\nu})} \left[ \sum\_{\nu=0}^{\infty} \frac{\sin(\gamma\_{\nu} y\_{D})}{N(\gamma\_{\nu})} \hat{\vec{p}}\_{mD}(\beta\_{\nu\nu} \gamma\_{\nu}) \right] \tag{A11}$$

Equation (A11) could be further expressed as

$$\begin{split} p\_{\text{mD}} &= \underbrace{\frac{4w\_{\text{rD}}}{\text{r}\_{\text{rD}}g\_{\text{rD}}}}\_{\text{r-1}} \underbrace{\frac{\sin(\gamma\_{\text{r}}y\_{\text{r}})[1-\cos(\gamma\_{\text{r}}y\_{\text{r}})]}{\mathcal{V}^{\text{r}}} \underbrace{\sum\_{\nu=1}^{\infty} \sin(\beta\_{\text{r}}x\_{\text{rD}}) \frac{\beta\_{\text{r}}[1-(-1)^{\nu}]}{(\beta\_{\text{r}}^{2}+\nu\_{\text{r}}^{2})}}}\_{\text{V}} \\ &+ \underbrace{\frac{4w\_{\text{r}O}}{\text{s.t.}} \sum\_{\nu=1}^{\infty} \frac{\sin(\beta\_{\text{r}}x\_{\text{rD}})[1-\cos(\beta\_{\text{r}}x\_{\text{rD}})]}{\beta\_{\text{s}}} \underbrace{\sum\_{\nu=1}^{\infty} \sin(\gamma\_{\text{r}}y\_{\text{rD}}) \frac{\gamma\_{\text{r}}[1-(-1)^{\nu}]}{(\gamma\_{\text{r}}^{2}+\beta\_{\text{r}}^{2})}}\_{\text{V}\_{2}} \\ &\left(\sum\_{\nu=1}^{\infty} \sin(\beta\_{\text{r}}x\_{\text{rD}}) \sin(\beta\_{\text{r}}(x\_{\text{rD}}+u\_{\text{rD}})) \underbrace{\sum\_{\nu=1}^{\infty} \sin(\beta\_{\text{r}}x\_{\text{rD}}) \sin(\beta\_{\text{r}}x\_{\text{rD}}+u\_{\text{rD}})}\_{\text{V}\_{2}}\right)}\_{\text{V}\_{3}} \end{split} \tag{A12}$$

where

+ 8π *xeD yeD*

$$V\_1 = \sum\_{\nu=1}^{\text{co}} \frac{\beta\_{\nu} \sin(\beta\_{\nu} \mathbf{x}\_{\text{D}})}{\beta\_{\nu}^2 + \gamma\_{\nu}^2} + \sum\_{\nu=1}^{\text{co}} (-1)^{\nu=1} \frac{\beta\_{\nu} \sin(\beta\_{\nu} \mathbf{x}\_{\text{D}})}{\beta\_{\nu}^2 + \gamma\_{\nu}^2} = \frac{\mathbf{x}\_{\text{cD}}}{2} \frac{\sinh[\nu \pi (\mathbf{x}\_{\text{cD}} - \mathbf{x}\_{\text{D}})/y\_{\text{cD}}] + \sinh(\nu \pi \mathbf{x}\_{\text{D}}/y\_{\text{cD}})}{\sinh(\nu \pi \mathbf{x}\_{\text{cD}}/y\_{\text{cD}})} \tag{A13}$$

$$V\_2 = \frac{y\_{eD}}{2} \frac{\sinh\left[\nu\pi(y\_{eD} - y\_D)/\mathbf{x}\_{eD}\right] + \sinh(\nu\pi y\_D/\mathbf{x}\_{eD})}{\sinh(\nu\pi y\_{eD}/\mathbf{x}\_{eD})}\tag{A14}$$

<sup>ν</sup>=1

6-------------------------------78-------------------------------9

*V*3

$$V\_{\eth} = \frac{\mathbf{x}\_{\text{cD}} y\_{\text{cD}}}{4\upsilon\pi} \frac{\cosh\left[\upsilon\pi \left(y\_{\text{cD}} - \left|y\_{\text{D}} - y\_{\text{of}D}\right|\right)/\mathbf{x}\_{\text{cD}}\right] - \cosh\left[\upsilon\pi \left[y\_{\text{cD}} - \left(y\_{\text{D}} + y\_{\text{of}D}\right)\right]/\mathbf{x}\_{\text{cD}}\right]}{\sinh\left(\upsilon\pi y\_{\text{cD}}/\mathbf{x}\_{\text{cD}}\right)}\tag{A15}$$

Correspondingly, Equation (A12) is finalized by

$$p\_{\rm nD}(\mathbf{x}\_{\rm D}, y\_{\rm D}) = p\_{\rm tDx}lB(\mathbf{x}\_{\rm D}, y\_{\rm D}) + p\_{\rm tDy}lD(\mathbf{x}\_{\rm D}, y\_{\rm D}) + 2\sum\_{m=1}^{N\_f} \int\_{0}^{\mathcal{L}\_{f\rm Dm}} q\_{f\rm Dm}(\mathbf{u}\_{\rm D}) \delta p\_{\rm tD}(\mathbf{x}\_{\rm D}, y\_{\rm D}, \boldsymbol{\mu}\_{\rm D})d\mathbf{u}\_{\rm D} \tag{A16}$$

where

$$\begin{cases} \begin{aligned} {}^{B}I &= \frac{2}{\pi} \sum\_{\nu=1}^{\infty} \frac{1}{\nu} \sin\left(\frac{\nu \pi y\_{D}}{\mathcal{Y}\_{dD}}\right) \frac{[1+(-1)^{\nu-1}]}{\sinh(\nu \pi x\_{D}/\mathcal{Y}\_{dD})} \Big(\sinh\frac{\nu \pi (x\_{dD} - \mathbf{x}\_{D})}{\mathcal{Y}\_{dD}} + \sinh\frac{\nu \pi x\_{D}}{\mathcal{Y}\_{dD}}\Big) \\\ {}^{B}I &= \frac{2}{\pi} \sum\_{\nu=1}^{\infty} \frac{1}{\nu} \sin\left(\frac{\nu \pi x\_{D}}{\mathcal{X}\_{dD}}\right) \frac{[1+(-1)^{\nu-1}]}{\sinh(\nu \pi y\_{dD}/\mathcal{X}\_{dD})} \Big(\sinh\frac{\nu \pi (y\_{D} - y\_{D})}{\mathcal{X}\_{dD}} + \sinh\frac{\nu \pi y\_{D}}{\mathcal{X}\_{dD}}\Big) \\\ \delta p\_{\nu D} &= \sum\_{\nu=1}^{\infty} \frac{4}{\nu} \frac{\sinh(\wp \pi\_{\text{D}}) \sin[\wp \pi (x\_{\nu (Dm} + u\_{D)})]}{\sinh(\nu \pi y\_{\text{D}}/\mathcal{X}\_{\text{dD}})} \cosh n\pi \frac{y\_{\text{D}} - [y\_{\text{D}} \pm y\_{\text{y}Dm}]}{\mathbf{x}\_{\text{eD}}} \end{aligned} \tag{A17}$$

#### **Appendix B. Analytical Solution for the Outer Region**

On the interface between the inner SRV region and outer elliptical region, the continuity condition is given as

$$\underbrace{\int\_{0}^{\pi} \left. \frac{\partial p\_{mD}}{\partial \xi\_{D}} \right|\_{\xi\_{D} = 0} d\eta\_{D}}\_{\text{elliptir flow}} = \underbrace{\int\_{0}^{\chi\_{\text{cD}}} \left. \frac{\partial p\_{mD}}{\partial y\_{D}} \right|\_{y\_{D} = y\_{\text{cD}}} d\mathbf{x}\_{D}}\_{\text{linear flow}} \tag{A18}$$

In detail, the right hand side (RHS) of Equation (A18) is expressed as

$$RHS = -p\_{cDn}A\_2 - p\_{cDy} \left(\frac{1}{\mathbf{x}\_{cD}\ln[\mathbf{x}\_{RD} + \sqrt{\mathbf{x}\_{RD}^2 - \mathbf{x}\_{cD}^2}/\mathbf{x}\_{cD}]} + B\_2\right) + \sum\_{m=1}^{N\_f} q\_{wDm} \mathbf{C}\_{2,m} \tag{A19}$$

where

$$\begin{cases} A\_{2} = \frac{4}{\pi x\_{\mathrm{r}D}} \sum\_{\nu=1}^{\infty} \frac{1}{\nu} \frac{[1+(-1)^{\nu-1}]}{\sinh(\nu \pi x\_{\mathrm{r}D}/y\_{\mathrm{r}D})} \Big(\cosh\frac{\nu \pi x\_{\mathrm{r}D}}{y\_{\mathrm{r}D}} - 1\Big) \\\ B\_{2} = \frac{-1}{x\_{\mathrm{r}D} \ln((\pi x\_{\mathrm{r}D} + \sqrt{\pi\_{\mathrm{RD}}^{2} - x\_{\mathrm{r}D}^{2}})/x\_{\mathrm{r}D})} - \frac{2}{\pi x\_{\mathrm{r}D}} \sum\_{\nu=1}^{\infty} \frac{1}{\nu} \frac{[1+(-1)^{\nu-1}]}{\sinh(\nu \pi y\_{\mathrm{r}D}/x\_{\mathrm{r}D})} \Big(\cosh\frac{\nu \pi y\_{\mathrm{r}D}}{x\_{\mathrm{r}D}} - 1\Big) \\\ C\_{2,m} = -\frac{8}{\pi \pi\_{\mathrm{f}} f \nu u} \sum\_{\nu=1}^{\infty} \frac{1}{\nu^{2}} \frac{[1-(-1)^{\nu}]}{\sinh(\nu \pi y\_{\mathrm{r}D}/x\_{\mathrm{r}D})} \sin\left(\frac{\nu \pi L\_{f\mathrm{Dm}}}{2x\_{\mathrm{r}D}}\right) \sin\left(\frac{\nu \pi (x\_{\mathrm{r}fD\mathrm{m}} + 0.5L\_{f\mathrm{Dm}})}{x\_{\mathrm{r}D}}\right) \sinh \nu \pi t \frac{y\_{\mathrm{r}D\mathrm{m}}}{x\_{\mathrm{r}D}} \end{cases} \tag{A20}$$

Thus, Equation (A18) is written as

$$A\_2 p\_{eDx} + B\_2 p\_{eyD} = \sum\_{m=1}^{N\_f} q\_{wDm} \mathcal{C}\_{2,m} \tag{A21}$$

Likewise, on the interface between the inner SRV region and outer linear region, the continuity condition is given as

$$\frac{1}{\mathcal{Y}^{\rm eD}} \int\_{0}^{y\_{\rm eD}} \left. \frac{\partial p\_{mD}}{\partial \mathbf{x}\_{\rm D}} \right|\_{\mathbf{x}\_{\rm D} = \mathbf{x}\_{\rm eD}} dy\_{\rm D} = \frac{1}{\mathcal{Y}^{\rm eD}} \int\_{0}^{y\_{\rm eD}} \left. \frac{\partial p\_{mD}}{\partial \mathbf{x}\_{\rm D}} \right|\_{\mathbf{x}\_{\rm D} = \mathbf{x}\_{\rm d}} dy\_{\rm D} \tag{A22}$$

In detail, the right hand side (RHS) of Equation (A22) is expressed as

$$RHS = -p\_{cDx} \left(\frac{1}{\mathbf{x}\_{RD} - \mathbf{x}\_{cD}} - A\_1\right) - p\_{cDy}B\_1 - \sum\_{m=1}^{N\_f} q\_{wDm}C\_{1,m} \tag{A23}$$

where

⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩

$$\begin{split} A\_{1} &= \frac{-1}{\left(\mathbf{x}\_{\mathrm{dO}} - \mathbf{x}\_{\mathrm{dO}}\right)} - \frac{1}{\pi y\_{\mathrm{dO}}} \sum\_{v=1}^{\infty} \frac{1}{v} \frac{\left[1 + (-1)^{v-1}\right]^{2}}{\sinh(v\pi x\_{\mathrm{dO}}/y\_{\mathrm{dO}})} \Big(\cosh\frac{v\pi x\_{\mathrm{dO}}}{y\_{\mathrm{dO}}} - 1\right) \\ B\_{1} &= \frac{2}{\pi y\_{\mathrm{dO}}} \sum\_{v=1}^{\infty} \frac{1}{v} \frac{\left[1 + \left(\frac{1 + (-1)^{v-1}}{\pi x\_{\mathrm{dO}}/x\_{\mathrm{dO}}}\right) \left(\cosh\frac{v\pi x\_{\mathrm{dO}}}{x\_{\mathrm{dO}}} - 1\right)}{\sinh(v\pi x\_{\mathrm{dO}}/x\_{\mathrm{dO}}) \sin\frac{v\pi L\_{\mathrm{dO}}}{v-1}} - 1 \end{split} \tag{A24}$$
 
$$\mathbf{C}\_{1,m} = -\frac{4x\_{\mathrm{dO}}}{\pi y\_{\mathrm{d}} L\_{f\Omega m}^{s}} \sum\_{v=1}^{\infty} \frac{1}{v^{2}} \frac{(-1)^{v}}{\sinh(v\pi y\_{\mathrm{d}O}/x\_{\mathrm{dO}})} \sin\left(\frac{v\pi L\_{f\Omega m}}{2x\_{\mathrm{dO}}}\right) \sin\left(\frac{v\pi (x\_{\mathrm{d}O\alpha} + 0.5L\_{f\Omega m})}{x\_{\mathrm{dO}}}\right) \sinh\frac{v\pi (y\_{\mathrm{d}O\alpha} - y\_{\mathrm{d}\Omega m})}{x\_{\mathrm{dO}}}$$

Thus, Equation (A22) is written as

$$A\_2 p\_{\rm rDx} + B\_2 p\_{\rm rydD} = \sum\_{m=1}^{N\_f} q\_{\rm wDm} \text{C}\_{2,m} \tag{A25}$$
