*Article* **Productivity Models of Infill Complex Structural Wells in Mixed Well Patterns**

#### **Liang Sun \*, Baozhu Li and Yong Li**

Research Institute of Petroleum Exploration & Development, PetroChina, Beijing 100083, China; lbz@petrochina.com.cn (B.L.); liyongph@petrochina.com.cn (Y.L.)

**\*** Correspondence: sunliang328@petrochina.com.cn

Received: 14 May 2019; Accepted: 27 May 2019; Published: 31 May 2019

**Abstract:** The mathematical models of productivity calculation for complex structural wells mainly focus on the single well or the regular well pattern. Previous research on the seepage theory of complex structural wells and vertical wells in mixed well pattern is greatly insufficient. Accordingly, this article presents a methodology of evaluating the productivity of infill complex structural wells in mixed well patterns. On the basis of the mirror-image method and source–sink theory, two semi-analytical models are established. These models are applied to the productivity prediction of an infill horizontal well inhorizontal-vertical well pattern and an infill multilateral well inmultilateral-vertical well pattern, respectively, in which the interference of other wells, the randomicity of well patterns, and the pressure drawdown along the horizontal laterals are taken into account. The semi-analytical models' results are consistent with those calculated by the Eclipse reservoir simulator with the relative error of less than 15%. Results indicate that the bottom hole flowing pressure decreases logarithmically while the wellbore flow rate increases monotonically from the toe to the heel of the horizontal well. Due to the pseudo-hemispherical flow at each endpoint and the pseudo-linear flow at the center of the horizontal well, the drainage area at each endpoint is relatively larger than that at the center. The radial inflow at each endpoint of the horizontal segment is considerably greater than that at the center, which presents the U-shape distribution. The proposed methodology enhances and promotes the theory of productivity evaluation for complex structural wells in mixed well patterns.

**Keywords:** complex structural well; mixed well pattern; productivity evaluation; semi-analytical model; well location optimization

#### **1. Introduction**

Complex structural wells including horizontal wells and multilateral wells have become a popular alternative for the development of oil and gas fields around the world because of their high flow efficiency due to larger contact area made with the reservoir and lower pressure drawdown at the same liquid volume [1]. As the process of oilfield development enters into the intermediary and later phase, the mixed well patterns of complex structural wells and conventional vertical wells have been widely applied to the implementation of adjustment plans [2]. The infilling horizontal wells or multilateral wells in mixed well patterns are of great significance to optimization of development strategies. Due to the complexity of seepage mechanism near the horizontal wellbore, the coupling between reservoir flow and wellbore conduit flow, and the interference of other wells in mixed well patterns [3], the original mathematical models based on time invariant flow are no longer applicable as a result of the change of flow regimes. Therefore, it is necessary to establish new productivity models of infill complex structural wells in mixed well patterns.

The productivity evaluation of complex structural wells under different reservoir conditions is an important topic in the field of complex structural wells. At present, studies have been conducted on the methods of production calculation for horizontal wells [4–18] and multilateral wells [19–23], which are mainly divided into analytical methods and semi-analytical methods. The analytical model aims to directly build a calculation formula based on ideal assumptions. In the semi-analytical model, each branch of the complex structural well is divided into several infinitesimal sections, thus the productivity can be obtained by solving the system of linear equations combined with fluid flow rate and pressure drawdown for each infinitesimal in the wellbore. However, these methods mainly focus on building models of single well or regular well patterns, and little research has been conducted on the seepage theory of mixed well patterns of vertical wells and complex structural wells [24–30]. In view of the field problems involving the productivity prediction of infill wells in irregular well patterns, there are limitations in current mathematical models [31–35]. On the one hand, they are only suitable for the productivity calculation of the entire well pattern not for single infill wells; on the other hand, they are merely applied to the regular five-spot pattern, seven-spot pattern, and nine-spot pattern not for the irregularly mixed well pattern. More importantly, the pressure drop caused by wall friction and fluid acceleration along the horizontal lateral is not comprehensively taken into account during the coupling of reservoir seepage and wellbore conduit flow. In order to overcome the deficiencies of the existing models, Ye et al. [36] presented a productivity evaluation model for infill horizontal wells considering the interference of other wells and the wellbore friction, which is suitable for mixed horizontal injection and production patterns. Although this method did not take into consideration more complicated conditions, such as the mixed well pattern including multilateral-vertical wells, it provides us with an effective approach to solve these problems.

The objective of this article here is to present two semi-analytical models for the productivity evaluation of infill complex structural wells in mixed well patterns on the basis of the mirror-image method and source–sink theory. The first model is suitable for the infill horizontal well in horizontal-vertical well pattern, and the other for the infill multilateral well in multilateral-vertical well pattern. Then, the two models are applied to the study the seepage mechanism in terms of the bottom hole flowing pressure and the distribution of wellbore flow and radial flow along the horizontal segment. The main feature of this methodology is that the models take into account the interference of other wells, the randomicity of well patterns, and the pressure drawdown along the horizontal laterals. The application in the optimization of infill well location also indicates the significant practical value of the proposed models.

#### **2. Productivity Model**

#### *2.1. Productivity Model of Mixed Horizontal-Vertical Well Pattern*

#### 2.1.1. Reservoir Flow Model

As is shown in Figure 1, in the deployment of mixed horizontal-vertical well pattern, we assume that the number of vertical producers including *Pz* <sup>1</sup>, *Pz* <sup>2</sup>, ··· , *<sup>P</sup><sup>z</sup> <sup>m</sup>* is *mz*, the number of horizontal producer including *P<sup>s</sup>* <sup>1</sup>, *Ps* <sup>2</sup>, ··· , *Ps <sup>m</sup>* is *ms* , the number of vertical injectors including *I z* 1, *I z* <sup>2</sup>, ··· , *I z <sup>n</sup>* is *nz*, the number of horizontal injectors including *I s* 1,*I s* <sup>2</sup>, ··· ,*I s <sup>n</sup>* is *ns* , and *Pnew* is the infill horizontal producer. The top and bottom boundaries of the reservoir are both closed, and the surrounding area is infinite. In order to analyze the well performance, the vertical interval and horizontal interval are both divided into *N* segments (Figure 2) [2,36].

**Figure 1.** The deployment of mixed horizontal-vertical well pattern.

**Figure 2.** Schematic of the segmented horizontal well.

The coordinates of an arbitrary point *M* in the vertical segment *k* (1 ≤ *k* ≤ *N*) of the *u*th vertical producer can be expressed as follows:

$$\begin{cases} \mathbf{x}\_z^p(\boldsymbol{\mu}, k, t) = \mathbf{x}\_z^p \\ y\_z^p(\boldsymbol{\mu}, k, t) = y\_{z\prime}^p \\ \vdots \\ \mathbf{z}\_z^p(\boldsymbol{\mu}, k, t) = z\_z^p + (k + t - 1)L/N, \qquad 0 \le t \le 1 \end{cases} \tag{1}$$

The coordinates of an arbitrary point *M* in the horizontal segment *k* (1 ≤ *k* ≤ *N*) of the *u*th horizontal producer can be expressed as follows:

$$\begin{cases} \mathbf{x}\_s^p(u,k,t) = \mathbf{x}\_s^p + (k+t-1)L/N, \\ y\_s^p(u,k,t) = y\_s^p, \\ \vdots \\ \mathbf{z}\_s^p(u,k,t) = \mathbf{z}\_s^p, \mathbf{x}\_s & 0 \le t \le 1 \end{cases} \tag{2}$$

*Processes* **2019**, *7*, 324

The coordinates of an arbitrary point *M* in the vertical segment *k* (1 ≤ *k* ≤ *N*) of the *v*th vertical injector can be expressed as follows:

$$\begin{cases} x\_z^i(v, k, t) = x\_{z'}^i \\ y\_z^i(v, k, t) = y\_{z'}^i \\ \vdots \\ z\_z^i(v, k, t) = z\_z^i + (k + t - 1)L/N, & 0 \le t \le 1 \end{cases} \tag{3}$$

The coordinates of an arbitrary point *M* in the horizontal segment *k* (1 ≤ *k* ≤ *N*) of the *v*th horizontal injector can be expressed as follows:

$$\begin{cases} \mathbf{x}\_s^i(\upsilon, k, t) = \mathbf{x}\_s^i + (k + t - 1)L/N, \\ y\_s^i(\upsilon, k, t) = y\_{s'}^i \\ \vdots \\ z\_s^i(\upsilon, k, t) = z\_{s'}^i & 0 \le t \le 1 \end{cases} \tag{4}$$

The coordinates of an arbitrary point *M* in the horizontal segment *k* (1 ≤ *k* ≤ *N*) of the infill horizontal producer can be expressed as follows [37]:

$$\begin{cases} \begin{array}{l} \chi(new,k,t) = \mathfrak{x}\_{new}^{p} + (k+t-1)L/N, \\ y(new,k,t) = y\_{new}^{p} \\ \vdots \\ z(new,k,t) = \mathfrak{z}\_{new}^{p} \end{array} \end{cases} \tag{5}$$

where *x p <sup>z</sup>* (*u*, *k*, *t*), *y p <sup>z</sup>* (*u*, *k*, *t*), *z p <sup>z</sup>* (*u*, *k*, *t*) are the coordinates of an arbitrary point in the segment *k* of the *u*th vertical producer; *x p <sup>z</sup>* , *y p <sup>z</sup>* , *z p <sup>z</sup>* are the coordinates of the left end (well heel) of the *u*th vertical producer; *x p <sup>s</sup>*(*u*, *k*,*t*), *y p <sup>s</sup>*(*u*, *k*,*t*), *z p <sup>s</sup>*(*u*, *k*,*t*) are the coordinates of an arbitrary point in the segment *k* of the *u*th horizontal producer; *x p <sup>s</sup>*, *y p <sup>s</sup>*, *z p <sup>s</sup>* are the coordinates of the left end (well heel) of the *u*th horizontal producer; *x<sup>i</sup> <sup>z</sup>*(*v*, *k*,*t*), *y<sup>i</sup> <sup>z</sup>*(*v*, *k*,*t*), *zi <sup>z</sup>*(*v*, *k*,*t*) are the coordinates of an arbitrary point in the segment *k* of the *v*th vertical injector; *xi <sup>z</sup>*, *y<sup>i</sup> <sup>z</sup>*, *zi <sup>z</sup>* are the coordinates of the left end (well heel) of the *v*th vertical injector; *xi <sup>s</sup>*(*v*, *k*,*t*), *yi <sup>s</sup>*(*v*, *k*,*t*), *zi <sup>s</sup>*(*v*, *k*,*t*) are the coordinates of an arbitrary point in the segment *k* of the *v*th horizontal injector; *xi <sup>s</sup>*, *yi <sup>s</sup>*, *zi <sup>s</sup>* are the coordinates of the left end (well heel) of the *v*th horizontal injector; *x*(*new*, *k*, *t*), *y*(*new*, *k*, *t*), *z*(*new*, *k*, *t*) are the coordinates of an arbitrary point in the segment *k* of the infill horizontal producer; *x p new*, *y p new*, *z p new* are the coordinates of the left end (well heel) of the infill horizontal producer. *L* is the length of the entire horizontal section, m; *N* is the total number of the segments, dimensionless.

The mirror-image method and source–sink theory have been widely applied to the boundary effect reduction. As is shown in Figure 3a, assuming the constant pressure boundary is a mirror, we can project the producer (we call sink in the method) into an image injector (we call source in the method) at the symmetric coordinate to counteract the constant pressure boundary effect. As is shown in Figure 3b, in the same way, we can project the producer into an image producer at the symmetric coordinate to counteract the closed boundary effect [2].

**Figure 3.** The mirror-image method and source–sink theory applied in different boundaries: (**a**) constant pressure boundary; (**b**) closed boundary.

On the basis of the method of images and the principle of superposition, we take the closed boundaries of the top and the bottom of the reservoir as mirrors. Thus, it is transformed into the problem of infinite well rows for horizontal producers, which is easy to solve. The potential of the producing segments of the infill horizontal well at any point *M* (*x*, *y*, *z*) in the infinite formation is:

$$\Phi(\mathbf{x}, y, z) = \sum\_{k=1}^{N} \Phi\_k(\mathbf{x}, y, z) = \frac{1}{4\pi\Lambda L\_k} \sum\_{k=1}^{N} \left[ q\_r(k)\varphi\_k(\mathbf{x}, y, z) \right] + \text{C.} \tag{6}$$

With:

<sup>ϕ</sup>*i*(*x*, *<sup>y</sup>*, *<sup>z</sup>*) = <sup>+</sup><sup>∞</sup> *n*=−∞ ξ*k*[*x*(*new*, *k*, *t*), *y*(*new*, *k*, *t*), 2*nh* + *z*(*new*, *k*, *t*)]+ ξ*k*[*x*(*new*, *k*, *t*), *y*(*new*, *k*, *t*), 2*nh* − *z*(*new*, *k*, *t*)]+ *mz u*=1 *q p <sup>u</sup>*,*z*(*k*) *qr*(*k*) <sup>ξ</sup>*k*[*<sup>x</sup> p <sup>z</sup>* (*u*, *k*, *t*), *y p <sup>z</sup>* (*u*, *k*, *t*), 2*nh* + *z p <sup>z</sup>* (*u*, *k*, *t*)]+ *mz u*=1 *q p <sup>u</sup>*,*z*(*k*) *qr*(*k*) <sup>ξ</sup>*k*[*<sup>x</sup> p <sup>z</sup>* (*u*, *k*, *t*), *y p <sup>z</sup>* (*u*, *k*, *t*), 2*nh* − *z p <sup>z</sup>* (*u*, *k*, *t*)]+ *ms u*=1 *q p <sup>u</sup>*,*s*(*k*) *qr*(*k*) <sup>ξ</sup>*k*[*<sup>x</sup> p <sup>s</sup>*(*u*, *k*, *t*), *y p <sup>s</sup>*(*u*, *k*, *t*), 2*nh* + *z p <sup>s</sup>*(*u*, *k*, *t*)]+ *ms u*=1 *q p <sup>u</sup>*,*s*(*k*) *qr*(*k*) <sup>ξ</sup>*k*[*<sup>x</sup> p <sup>s</sup>*(*u*, *k*, *t*), *y p <sup>s</sup>*(*u*, *k*, *t*), 2*nh* − *z p <sup>s</sup>*(*u*, *k*, *t*)]− *nz v*=1 *qi <sup>v</sup>*,*z*(*k*) *qr*(*k*) <sup>ξ</sup>*k*[*xi <sup>z</sup>*(*v*, *k*, *t*), *y<sup>i</sup> <sup>z</sup>*(*v*, *k*, *t*), 2*nh* + *zi <sup>z</sup>*(*v*, *k*, *t*)]− *nz v*=1 *qi <sup>v</sup>*,*z*(*k*) *qr*(*k*) <sup>ξ</sup>*k*[*xi <sup>z</sup>*(*v*, *k*, *t*), *y<sup>i</sup> <sup>z</sup>*(*v*, *k*, *t*), 2*nh* − *zi <sup>z</sup>*(*v*, *k*, *t*)]− *ns v*=1 *qi <sup>v</sup>*,*s*(*k*) *qr*(*k*) <sup>ξ</sup>*k*[*xi <sup>s</sup>*(*v*, *k*, *t*), *yi <sup>s</sup>*(*v*, *k*, *t*), 2*nh* + *z<sup>i</sup> <sup>s</sup>*(*v*, *k*, *t*)]− *ns v*=1 *qi <sup>v</sup>*,*s*(*k*) *qr*(*k*) <sup>ξ</sup>*k*[*x<sup>i</sup> <sup>s</sup>*(*v*, *k*, *t*), *yi <sup>s</sup>*(*v*, *k*, *t*), 2*nh* − *zi <sup>s</sup>*(*v*, *<sup>k</sup>*, *<sup>t</sup>*)] , *n* = 0,±1,±2, ··· (7)

where *q p <sup>u</sup>*,*z*(*k*) <sup>≈</sup> *Qp u*,*z <sup>N</sup>* , *qi <sup>v</sup>*,*z*(*k*) <sup>≈</sup> *<sup>Q</sup><sup>i</sup> v*,*z <sup>N</sup>* , *q p <sup>u</sup>*,*s*(*k*) <sup>≈</sup> *Qp u*,*s <sup>N</sup>* , *qi <sup>v</sup>*,*s*(*k*) <sup>≈</sup> *Qi v*,*s <sup>N</sup>* ; *<sup>Q</sup><sup>p</sup> <sup>u</sup>*,*<sup>z</sup>* is the production of the *u*th vertical producer, m3/d; *Q<sup>i</sup> <sup>v</sup>*,*<sup>z</sup>* is the injection rate of the *<sup>v</sup>*th vertical injector; *Qp <sup>u</sup>*,*<sup>s</sup>* is the production of the *u*th horizontal producer, m3/d; *Q<sup>i</sup> <sup>v</sup>*,*<sup>s</sup>* is the injection rate of the *v*th horizontal injector, m3/d.

For:

$$\mathbb{E}\_{\mathbb{R}}[\mathbf{x}(new, k, t), y(new, k, t), z(new, k, t)] = \ln \frac{r\_1(new, k) + r\_2(new, k) + \frac{L}{N}}{r\_1(new, k) + r\_2(new, k) - \frac{L}{N}},$$

$$\begin{aligned} r\_1(n\kappa w, k) &= \sqrt{\left[\mathbf{x}(n\kappa w, k, t=0) - \mathbf{x}\right]^2 + \left[y(n\kappa w, k, t=0) - y\right]^2 + \left[z(n\kappa w, k, t=0) - z\right]^2}, \\ r\_2(n\kappa w, k) &= \sqrt{\left[\mathbf{x}(n\kappa w, k, t=1) - \mathbf{x}\right]^2 + \left[y(n\kappa w, k, t=1) - y\right]^2 + \left[z(n\kappa w, k, t=1) - z\right]^2}. \end{aligned}$$

In the same way, we can obtain the expressions of ξ*k*[*x*(*u*, *k*,*t*), *y*(*u*, *k*,*t*), *z*(*u*, *k*,*t*)] and ξ*k*[*x*(*v*, *k*, *t*), *y*(*v*, *k*, *t*), *z*(*v*, *k*, *t*)].

Once the distributions of the flow rates of *qr*(*k*), *q p <sup>u</sup>*,*z*(*k*), *qi <sup>v</sup>*,*z*(*k*), *q p <sup>u</sup>*,*s*(*k*), *q<sup>i</sup> <sup>v</sup>*,*s*(*k*) are known, all terms of flow pressure *pw f*(*k*) of the infill horizontal producer can be calculated. Hence Equation (6) is simplified to a linear equation system, and can be arranged in the form as follows:

$$
\begin{bmatrix}
\begin{matrix}
\varphi\_{11} & \varphi\_{12} & \varphi\_{13} & \cdots & \varphi\_{1N} \\
\varphi\_{21} & \varphi\_{22} & \varphi\_{23} & \cdots & \varphi\_{2N} \\
\varphi\_{31} & \varphi\_{32} & \varphi\_{33} & \cdots & \varphi\_{3N} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
\varphi\_{N1} & \varphi\_{N2} & \varphi\_{N3} & \cdots & \varphi\_{NN}
\end{matrix}
\end{bmatrix}
\begin{bmatrix}
q\_{r}(1) \\
q\_{r}(2) \\
q\_{r}(3) \\
\vdots \\
q\_{r}(N)
\end{bmatrix} = \frac{4\pi K\Delta L}{\mu\_{\text{o}}} \begin{bmatrix}
p\_{c} - p\_{wf}(1) \\
p\_{c} - p\_{wf}(2) \\
p\_{c} - p\_{wf}(3) \\
\vdots \\
p\_{c} - p\_{wf}(N)
\end{bmatrix} \tag{8}
$$

where ϕ*ij* is the value of ϕ*<sup>j</sup>* at the midpoint of the segment *i* of the infill horizontal producer, dimensionless; *qr*(*i*) is the radial flow rate of the segment *i* of the infill horizontal producer, m3/d; *pe* is the reservoir pressure, MPa; *pwf*(*i*) is the flowing pressure of the segment *i* of the infill horizontal producer, MPa; *<sup>K</sup>* is the reservoir permeability, 10−<sup>3</sup> <sup>μ</sup>m2; <sup>μ</sup>*<sup>o</sup>* is the oil viscosity, mPa·s; *<sup>h</sup>* is the net pay thickness, m.

The wellbore flow rate along the horizontal well can be expressed as follows:

$$q\_l(k) = \sum\_{j=k}^{N} q\_{l'}(j),\tag{9}$$

where *ql*(*k*) is the wellbore flow rate of the segment *k* of the infill horizontal producer, m3/d.

#### 2.1.2. Wellbore Flow Model

The wellbore pressures of each segment along the horizontal segment are not independent of each other, instead, they are related to each other via wellbore hydraulics. More specifically, the pressure difference between two adjacent segment midpoints is dependent on the radial flow into the two segments, the local pressure, and the fluid property [2,36].

The pressure drop along the horizontal well is caused by the wall friction and the fluid acceleration [3,11,12,38], and can be calculated by the correlation of pressure drop in the horizontal interval under the condition of open-hole completion proposed by Liu et al. [12]:

$$\begin{array}{l} \Delta p\_{wf}(i) = \frac{l}{N} \Big\{ \frac{2f\rho}{\pi^2 D^5} \Big[ 2q\_l(i) - q\_r(i) \Big]\_N^L \Big] + \frac{16\alpha q\_r(i)}{\pi^2 D^4} \Big[ 2q\_l(i) - q\_r(i) \Big]\_N^L \Big\} \\\ p\_{wf}(i) = p\_{wf}(i-1) + 0.5 \Big( \Delta p\_{wf}(i-1) + \Delta p\_{wf}(i) \Big) \Big. \quad (1 \le i \le N+1) \end{array} \tag{10}$$

where *pwf* is the flowing pressure of the well heel, *pwf*(0) = *pwf*, Δ*pwf*(0) = Δ*pwf*(*N* + 1) = 0, MPa; ρ is the fluid density, g/cm3; *f* is the wall friction factor, and can be calculated as follows:

$$\begin{aligned} f &= \frac{64}{N\_{\text{Re}}}, & N\_{\text{Re}} &\leq 2100\\ \frac{1}{\sqrt{f}} &= 1.14 - 2\lg\left(\frac{\tau}{D} + \frac{21.25}{N\_{\text{Re}}^{0.9}}\right), & N\_{\text{Re}} &\geq 2100 \end{aligned}$$

By combing Equations (9) and (10), the following system of equations can be obtained:

$$f[q\_r(i), p\_{wf}(i)] = 0.\tag{11}$$

As expected, the total number of Equations (8) and (11) is equal to 2*N*, which is the same as that of unknowns. Thus, the model has a unique solution.

#### 2.1.3. Solution Procedure

The iterative algorithm is an efficient way to solve this problem, and the solution procedure is shown in Figure 4.

**Figure 4.** The flow chart of solution procedure.

#### *2.2. Productivity Model of Mixed Multilateral-Vertical Well Pattern*

#### 2.2.1. Reservoir Flow Model

As is shown in Figure 5, in the deployment of mixed multilateral-vertical well pattern, we assume that the number of vertical producers including *Pz* <sup>1</sup>, *Pz* <sup>2</sup>, ··· , *Pz <sup>m</sup>* is *mz*, the number of vertical injectors including *I z* 1,*I z* <sup>2</sup>, ··· ,*I z <sup>n</sup>* is *nz*, and *Pnew* is the infill multilateral well. The top and bottom boundaries of the reservoir are both closed, and the surrounding area is infinite. In order to analyze the well performance, the vertical interval is divided into *N* segments.

**Figure 5.** The mixed well pattern of multilateral-vertical wells.

*Processes* **2019**, *7*, 324

The coordinates of an arbitrary point *M* in the vertical segment *k* (1 ≤ *k* ≤ *N*) of the *u*th vertical producer can be expressed as follows:

$$\begin{cases} \mathbf{x}\_z^p(\boldsymbol{\mu}, k, t) = \mathbf{x}\_{z\prime}^p\\ \mathbf{y}\_z^p(\boldsymbol{\mu}, k, t) = \mathbf{y}\_{z\prime}^p\\ \vdots\\ \mathbf{z}\_z^p(\boldsymbol{\mu}, k, t) = \mathbf{z}\_z^p + (k + t - 1)L/\mathcal{N}, \qquad 0 \le t \le 1 \end{cases} \tag{12}$$

The coordinates of an arbitrary point *M* in the vertical segment *k* (1 ≤ *k* ≤ *N*) of the *v*th vertical injector can be expressed as follows:

$$\begin{cases} \mathbf{x}\_z^i(v, k, t) = \mathbf{x}\_{z'}^i \\ y\_z^i(v, k, t) = y\_{z\prime}^i \\ \vdots \\ z\_z^i(v, k, t) = z\_z^i + (k + t - 1)L/N, \quad 0 \le t \le 1 \end{cases} \tag{13}$$

where *x p <sup>z</sup>* (*u*, *k*, *t*), *y p <sup>z</sup>* (*u*, *k*, *t*), *z p <sup>z</sup>* (*u*, *k*, *t*) are the coordinates of an arbitrary point in the segment *k* of the *u*th vertical producer; *x p <sup>z</sup>* , *y p <sup>z</sup>* , *z p <sup>z</sup>* are the coordinates of the left end (well heel) of the *u*th vertical producer; *xi <sup>z</sup>*(*v*, *k*, *t*), *y<sup>i</sup> <sup>z</sup>*(*v*, *k*, *t*), *z<sup>i</sup> <sup>z</sup>*(*v*, *k*, *t*) are the coordinates of an arbitrary point in the segment *k* of the *v*th vertical injector; *xi <sup>z</sup>*, *y<sup>i</sup> <sup>z</sup>*, *zi <sup>z</sup>* are the coordinates of the left end (well heel) of the *v*th vertical injector.

Figure 6 is the three-dimensional (3D) schematic of the infill multilateral well, we assume that each branch is symmetrically distributed on the same horizontal plane, the length of which is identical. The *xoy* coordinate system is established by taking the subpoint of the main hole in the *xoy* plane as the origin of coordinates, and the direction of the main hole as the *z* axis. Each branch is distributed counterclockwise around the *z* axis with the *x* axis as the starting direction.

**Figure 6.** 3D schematic of infill multilateral well.

*M*0(*x*, *y*, *z*0) is the bottom-hole coordinate of the main wellbore, and *Mi*(*xi*,0, *yi*.0, *zi*,0) is the starting coordinate of the *i*th branch with the length of *Li* (1 ≤ *i* ≤ *M*). *M* is the number of branches. The *i*th branch of the infill multilateral well is divided into *Ni* (1 ≤ *i* ≤ *M*) segments, thus the infinitesimal length of the *i*th branch is Δ*Li* = *Li*/*Ni*.

*Processes* **2019**, *7*, 324

Accordingly, the coordinates of an arbitrary point *M* in the horizontal segment *i* (1 ≤ *i* ≤ *M*) of the infill multilateral well can be expressed as follows:

$$\begin{cases} x(i,j) = x\_{i,0} + \Delta L\_i[\sum\_{k=1}^{j-1} (\sin \theta\_{ik} \cdot \cos \alpha\_{ik}) + t \sin \theta\_{i,j} \cdot \cos \alpha\_{ik}]\_\prime \\\ y(i,j) = y\_{i,0} + \Delta L\_i[\sum\_{k=1}^{\prime} (\sin \theta\_{ik} \cdot \sin \alpha\_{ik}) + t \sin \theta\_{i,j} \cdot \sin \alpha\_{ik}]\_\prime \\\ z(i,j) = z\_{i,0} + \Delta L\_i[\sum\_{k=1}^{\prime} \cos \theta\_{ik} + t \cos \theta\_{i,j}]\_\prime \quad (1 \le i \le M, 0 \le j \le N\_i, 0 \le t \le 1) \end{cases} \tag{14}$$

where *x*(*i*, *j*), *y*(*i*, *j*), *z*(*i*, *j*) are the coordinates of an arbitrary point in the *i*th branch of the infill multilateral producer; *xi*,0, *yi*,0, *zi*,0 are the coordinates of the left end (well heel) in the *i*th branch of the infill multilateral producer; θ is the deviation angle of the *i*th segment of the infill multilateral producer; α is the azimuth angle of the *i*th segment of the infill multilateral producer.

As mentioned above, we can also transform this problem into the one of infinite well rows by using the method of images and the principle of superposition. Thus the potential of the producing segments of the infill multilateral well at any point *M*(*x*, *y*, *z*) in the infinite formation is:

$$f(\mathbf{x}, y, z) = \sum\_{i=1}^{M} f\_i(\mathbf{x}, y, z) = \frac{1}{4\pi} \sum\_{i=1}^{M} \sum\_{j=1}^{N\_j} \left[ q\_{l'}(i, j) q\_{i, j}(\mathbf{x}, y, z) \right] + \text{C.} \tag{15}$$

With:

<sup>ϕ</sup>*i*,*j*(*x*, *<sup>y</sup>*, *<sup>z</sup>*) = <sup>+</sup><sup>∞</sup> *n*=−∞ [ 1 Δ*Li* ξ*i*,*<sup>j</sup> xi*,*j*, *yi*,*j*, 2*nh* + *zi*,*<sup>j</sup>* + <sup>1</sup> Δ*Li* ξ*i*,*<sup>j</sup> xi*,*j*, *yi*,*j*, 2*nh* − *zi*,*<sup>j</sup>* + 1 Δ*L mz u*=1 *q p <sup>u</sup>*,*z*(*k*) *qr*(*i*,*j*) <sup>ξ</sup>*k*[*<sup>x</sup> p <sup>z</sup>* (*u*, *k*, *t*), *y p <sup>z</sup>* (*u*, *k*, *t*), 2*nh* + *z p <sup>z</sup>* (*u*, *k*, *t*)]+ 1 Δ*L mz u*=1 *q p <sup>u</sup>*,*z*(*k*) *qr*(*i*,*j*) <sup>ξ</sup>*k*[*<sup>x</sup> p <sup>z</sup>* (*u*, *k*, *t*), *y p <sup>z</sup>* (*u*, *k*, *t*), 2*nh* − *z p <sup>z</sup>* (*u*, *k*, *t*)]− 1 Δ*L nz v*=1 *qi <sup>v</sup>*,*z*(*k*) *qr*(*i*,*j*) <sup>ξ</sup>*k*[*x<sup>i</sup> <sup>z</sup>*(*v*, *k*, *t*), *yi <sup>z</sup>*(*v*, *k*, *t*), 2*nh* + *zi <sup>z</sup>*(*v*, *k*, *t*)]− 1 Δ*L nz v*=1 *qi <sup>v</sup>*,*z*(*k*) *qr*(*i*,*j*) <sup>ξ</sup>*k*[*x<sup>i</sup> <sup>z</sup>*(*v*, *k*, *t*), *yi <sup>z</sup>*(*v*, *k*, *t*), 2*nh* − *z<sup>i</sup> <sup>z</sup>*(*v*, *k*, *t*)] (16)

where *q p <sup>u</sup>*,*z*(*k*) <sup>≈</sup> *Qp u*,*z <sup>N</sup>* , *qi <sup>v</sup>*,*z*(*k*) <sup>≈</sup> *Qi v*,*z <sup>N</sup>* ; *<sup>Q</sup><sup>p</sup> <sup>u</sup>*,*<sup>z</sup>* is the production of the *u*th vertical producer, m3/d; *Qi <sup>v</sup>*,*<sup>z</sup>* is the injection rate of the *v*th vertical injector, m3/d.

Once the distributions of the flow rates of *qr*(*i*, *j*), *q p <sup>u</sup>*,*z*(*k*), *qi <sup>v</sup>*,*z*(*k*) are known, all terms of flow pressure *pw f*(*i*, *j*) of the infill multilateral producer can be calculated. Hence Equation (15) is simplified to a linear equation system, and can be arranged in the form as follows:

$$A \begin{bmatrix} q\_r(1,1) \\ q\_r(1,2) \\ \vdots \\ q\_r(1,N\_1) \\ \vdots \\ q\_r(M\_rN\_M) \end{bmatrix} = \frac{4\pi\sqrt{k\_lh\_\upsilon}}{\mu\_o} \begin{bmatrix} P\_c - P\_{wf}(1,1) \\ P\_c - P\_{wf}(1,2) \\ \vdots \\ P\_c - P\_{wf}(1,N\_1) \\ \vdots \\ P\_c - P\_{wf}(1,N\_M) \end{bmatrix} \tag{17}$$

With:

$$A = \begin{bmatrix} \phi\_{(1,1)(1,1)} - \phi\_{\varepsilon(1,1)} & \cdots & \phi\_{(1,1)(1,N\_{1})} - \phi\_{\varepsilon(1,N)} & \phi\_{(1,1)(2,1)} - \phi\_{\varepsilon(2,1)} & \cdots & \phi\_{(1,1)(M\_{\mathrm{Mul}})} - \phi\_{\varepsilon(M\mathrm{Mul})} \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ \phi\_{(1,N\_{1})(1,1)} - \phi\_{\varepsilon(1,1)} & \cdots & \phi\_{(1,N\_{1})(1,N\_{1})} - \phi\_{\varepsilon(1,N)} & \phi\_{(1,N\_{1})(2,1)} - \phi\_{\varepsilon(2,1)} & \cdots & \phi\_{(1,N\_{1})(M\_{\mathrm{Mul}})} - \phi\_{\varepsilon(M\mathrm{Mul})} \\ \phi\_{(2,1)(1,1)} - \phi\_{\varepsilon(1,1)} & \cdots & \phi\_{(2,1)(1,N\_{1})} - \phi\_{\varepsilon(1,N)} & \phi\_{(2,1)(2,1)} - \phi\_{\varepsilon(2,1)} & \cdots & \phi\_{(2,1)(M\_{\mathrm{Mul}})} - \phi\_{\varepsilon(M\mathrm{Mul})} \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ \phi\_{(M\_{M})(1,1)} - \phi\_{\varepsilon(1,1)} & \cdots & \phi\_{(M\_{M})(M\_{\mathrm{M}})} - \phi\_{\varepsilon(1,N)} & \phi\_{(M\_{M})(2,1)} - \phi\_{\varepsilon(2,1)} & \cdots & \phi\_{(M\_{M})(M\_{\mathrm{Mul}})} - \phi\_{\varepsilon(M\mathrm{Mul})} \end{bmatrix}$$

where *qr*(*i*, *j*) is the radial flow rate of the *i*th branch of the infill multilateral producer, 1 ≤ *i* ≤ *M*, <sup>1</sup> <sup>≤</sup> *<sup>j</sup>* <sup>≤</sup> *Ni*, m3/d; *Pe* is the reservoir pressure, MPa; *Pwf*(*i*, *j*) is the flow pressure of the *i*th branch of the infill multilateral producer, 1 ≤ *i* ≤ *M*, 1 ≤ *j* ≤ *Ni*, MPa; *kh* is the reservoir horizontal permeability, <sup>10</sup>−<sup>3</sup> <sup>μ</sup>m2; *kv* is the reservoir vertical permeability, 10−<sup>3</sup> <sup>μ</sup>m2; <sup>μ</sup>*<sup>o</sup>* is the oil viscosity, mPa·s.

#### 2.2.2. Wellbore Flow Model

#### 1. Pressure drop model of the horizontal lateral

The pressure drop along the horizontal interval is caused by the wall friction and the fluid acceleration. According to Equation (11), the pressure drop of the *i*th branch of the multilateral well along the horizontal interval can be calculated by:

$$F\_{i,j}'[q\_r(i,j), p\_{wf}(i,j)] = 0.\tag{18}$$

#### 2. Pressure drop model of the deviated segment

The branch holes enter the reservoir through the deviated segment. Although each deviated segment does not contribute to the production, the pressure drop in the deviated segments cannot be ignored. According to fluid seepage theory in the elbow, the pressure drop can be calculated by the following formula:

$$
\Delta p\_s = \mathcal{C}\_w \frac{4 \times 10^{-6} \rho \lambda\_c R\_c Q\_W^2}{\pi d\_w^5} + \frac{\rho \mathcal{g} R\_c}{10^6} \,\mathrm{}\tag{19}
$$

where Δ*Ps* the pressure drop along the deviated segment, MPa; *Rc* is curvature radius of the deviated segment, m; ρ is the fluid density, kg/m3; *dw* is the wellbore diameter of the deviated segment, m; *Cw* is the similarity factor of the wellbore and deviated segment, dimensionless; λ*<sup>c</sup>* is the friction factor of the deviated segment, dimensionless.

#### 3. Pressure drop model of the vertical segment

According to the law of conservation of energy, the equation of pressure gradient in the pipe with inclination can be expressed as follows:

$$\frac{dP}{dz} = \rho \lg \sin \theta + \rho v \frac{dv}{dz} + f \frac{\rho v^2}{2d},\tag{20}$$

where *dP*/*dz* is the pressure loss per unit length, MPa/m; ρ*g*sinθ is the pressure drop caused by the fluid gravity, MPa/m; <sup>ρ</sup>*v* <sup>×</sup> *dv*/*dz* is the pressure drop caused by the fluid acceleration, MPa/m; *f* <sup>×</sup> <sup>ρ</sup>*v*2/*2d* is the pressure drop caused by the friction, MPa/m; ρ is the fluid density, kg/m3; *d* is the pipe diameter, m. The friction factor *f* can be calculated by:

 $f = \frac{64}{N\_{\text{Re}}}, N\_{\text{Re}} \le 2100$ 
 $\frac{1}{\sqrt{f}} = 1.14 - 2\lg\left(\frac{\tau}{d} + \frac{21.25}{N\_{\text{Re}}^{0.9}}\right), N\_{\text{Re}} \ge 2100$ 

where *NRe* is the Reynolds number, dimensionless.

#### 4. Internal pressure drop model of the branch hole

The flowing pressure at the heel of the *i*th branch can be calculated by simultaneous Equations (19) and (20):

$$P\_{wf}(i) = P\_{mb} - \Delta P\_{vi} - \Delta P\_{wi}(1 \le i \le M), \tag{21}$$

where *Pmb* is the flowing pressure of the main wellbore, MPa; Δ*Pvi* is the pressure drop of the vertical segment of the *i*th branch, MPa; Δ*Pwi* is the pressure drop of the deviated segment of the *i*th branch, MPa. By combining Equations (20) and (24), the following system of equations can be obtained:

> *Fi*,*<sup>j</sup>* 1 *qr*(*i*, *j*), *pw f*(*i*, *j*) 2 = 0,(1 ≤ *i* ≤ *M*, 1 ≤ *j* ≤ *Ni* . (22)

As expected, the total number of equations of Equations (17) and (22) is equal to 2(*M* × *Ni*), which is the same as that of unknowns. Therefore, the solution of the coupling model is unique.

#### 2.2.3. Solution Procedure

The coupling model can also be solved by using the iterative algorithm, and the solution procedure is shown in Figure 7.

**Figure 7.** The flow chart of solution procedure.

#### **3. Results and Discussion**

#### *3.1. Model Validation*

The semi-analytical models can be applied to the productivity prediction of the infill horizontal well and the multilateral well in mixed well patterns. In order to verify the reliability of the models, two cases are considered, and the results calculated by the semi-analytical models are compared with

those calculated by the Eclipse software (Schlumberger). In the first case, we selected a typical mixed well pattern of horizontal-vertical wells from a thin carbonate reservoir in the Middle East. There are three vertical producers (*P*11, *P*12, and *P*13), one horizontal producer (*P*14), one vertical injector (*I*11), one horizontal injector (*I*12), and *P*1*new* is the infill horizontal producer. In the other case, the reservoir is a thick, multi-layer reservoir in the Middle East, multilateral wells are adopted to improve the degree of the reservoir development. The mixed well pattern of multilateral-vertical well includes two vertical producers (*P*<sup>21</sup> and *P*22), one vertical injector (*I*21), and *P*2*new* is the infill multilateral producer with two branch holes.

The productivity models can be solved by VB modular programming of the object-oriented technology. The reservoir properties and wellbore geometry are listed in Tables 1 and 2. By inputting the same model parameters, the productivity evaluation of the two wells were conducted respectively with the semi-analytical models and the Eclipse simulator. The calculation results are shown in Table 3.


**Table 1.** Productivity calculation parameters of *P*1*new*.

**Table 2.** Productivity calculation parameters of *P*2*new*.


Table 3 indicates that the semi-analytical models have a high accuracy with the relative error of less than 15%. The results calculated by the models are consistent with those by the Eclipse numerical simulator, which verify the reliability of the productivity models. By taking into account the actual conditions, namely the interference of other wells in the well pattern, the coupling between the reservoir flow and wellbore flow, and the pressure drop along the horizontal lateral, the accuracy of the productivity prediction by the semi-analytical models is greatly improved.


**Table 3.** Results comparison of two methods.

#### *3.2. Model Application*

#### 3.2.1. Study on Seepage Mechanism of Horizontal Well

In the first case of model validation, we can further analyze the bottom hole flowing pressure, the distribution of wellbore flow rate and radial flow rate along the horizontal segment. Figures 8 and 9 show that the bottom hole flowing pressure decreases logarithmically while the wellbore flow rate increases monotonically from the toe to the heel of the horizontal well. The reservoir flow coupled with variable mass wellbore flow results in this near wellbore dynamics. On one hand, the mass flow rate of fluid from the toe to the heel increases gradually. In this case, the fluid velocity along the main flow direction also increases with an accelerating pressure drop. The radial inflow of reservoir fluid along the horizontal wellbore disturbs the boundary layer of the main stream tube and affects its velocity profile, thus altering the wall friction resistance determined by the velocity distribution. These factors lead to the bottom hole flowing pressure decreasing in the form of logarithmic function. On the other hand, the radial inflow performance affects the distribution of pressure and pressure drop, and vice versa. Due to the pseudo-hemispherical flow at each endpoint (well heel and well toe) and the pseudo-linear flow at the center of the horizontal segment, the drainage area at each endpoint is relatively larger than that at the center. The radial flow rate at each endpoint of the horizontal segment is considerably greater than that at the center, and decreases quickly toward the intermediate section, which generally presents a U-shape distribution.

**Figure 8.** Distribution of bottom hole flowing pressure along the horizontal lateral.

**Figure 9.** Distribution of wellbore flow and radial flow along the horizontal lateral.

#### 3.2.2. Optimization of Infilling Well Location

In the intermediary and later phase of oilfield development, infill wells are usually needed to improve the well pattern and the effect of tapping potential of remaining oil. Therefore, the key is the optimization of infill well location. The productivity models provide a shortcut for this target. Based on the understanding of geological conditions and development status, the productivity of infill well can be accurately predicted, and in combination with the findings of remaining oil distribution, the infill well location could be optimized.

Taking the mixed five-spot well pattern of vertical injectors and horizontal producer as an example, we conducted a study on the productivity variation of horizontal well at different locations in the well pattern (Figure 10). The basic reservoir parameters refer to the first case of model validation. Figure 11 shows that the lateral shifting of horizontal well has some effects on the productivity, and the productivity at the center of the well pattern is relatively higher than that at other locations. However, this does not mean that the central location is the optimum selection for the horizontal well. In the adjustment of well pattern, geological condition and development factors, such as the distribution of flow field and remaining oil, should also be considered.

**Figure 10.** Mixed five-spot well pattern of vertical-injectors and horizontal producer (-1 –-7 show the code of different horizontal locations).

**Figure 11.** Productivity variation of horizontal well in mixed well pattern.

We made further study on the flow field distribution in the mixed well pattern so as to provide references for the optimization of infilling well location. Figures 12 and 13 indicate that the equipotential lines and stream lines are symmetrically distributed along the *x*-axis as the horizontal well moves laterally. When the horizontal well is located at the center of the well pattern, both the equipotential lines and stream lines present a uniform distribution; when the horizontal well moves toward the right side, both of the two types of the lines on the right side become dense, while those on the left side are sparsely distributed, and vice versa. The more the horizontal well location changes, the more obvious this phenomenon is. The flow pattern almost presents a linear flow as the horizontal well locates at the line of the two vertical wells on the same side, which indicates a relatively higher swept efficiency. Therefore, in combination with the findings of productivity variation and flow fluid distribution, the infilling well location could be optimized. This improves the effect of well pattern adjustment considerably.

**Figure 12.** Equipotential line pattern for lateral shifting of the horizontal well. ((**a**–**c**) show the equipotential line pattern of different lateral location.)

**Figure 13.** Stream line pattern for lateral shifting of the horizontal well. ((**a**–**c**) show the stream line pattern of different lateral location.)

#### *3.3. Discussion*

In this study, the semi-analytical models for productivity evaluation of infill complex structural wells are built with the method of images and source–sink theory. Compared with the existing models, the advantages of the semi-analytical model lie in the following three aspects. Firstly, the models are suitable for mixed well patterns of vertical wells and complex structure wells; Secondly, both the interference of other wells in the well pattern and the pressure drawdown along the horizontal segment are taken into account, which makes the models more reliable. In addition, the models provide some insight into the seepage mechanism of complex structural wells such as the distribution of wellbore flow and radial flow along the horizontal segment. However, the semi-analytical models also have some limitations, the basic assumptions of the models are steady-state flow and closed top and bottom boundaries, and the productivity evaluation of different reservoir boundaries under the condition of unsteady-state flow still needs to be further studied.

#### **4. Conclusions**

In this paper, on the basis of the mirror-image method and the source–sink theory, two semi-analytical models are established to predict the productivity of infill horizontal wells and multilateral wells in mixed well patterns respectively, in which the interference of other wells, the randomicity of well pattern, and the pressure drawdown along the horizontal lateral are taken into account.

The application of the productivity models verifies the reliability and practicability. The results calculated by the models are consistent with those by the Eclipse numerical simulator, and the semi-analytical models have a high accuracy with the relative error of less than 15%.

The results indicate that the bottom hole flowing pressure decreases logarithmically while the wellbore flow rate increases monotonically from the toe to the heel of the horizontal segment. Due to the pseudo-hemispherical flow at each endpoint and the pseudo-linear flow at the center of the horizontal segment, the drainage area at each endpoint is relatively larger than that at the center. The radial flow rate at each endpoint of the horizontal segment is considerably greater than that at the center, which generally presents a U-shape distribution.

This study also proposes a practical and efficient approach for the study on the horizontal seepage mechanism, and the optimization of infilling well locations.

**Author Contributions:** Conceptualization, L.S.; Data curation, L.S.; Formal analysis, L.S.; Funding acquisition, Y.L.; Methodology, L.S.; Project administration, B.L.; Supervision, Y.L.; Validation, L.S.; Writing – review & editing, L.S., Y.L.

**Funding:** This research was founded by the National Science and Technology Major Project of China, grant number 2016ZX05015-002, and the National Natural Science Foundation of China, grant number 2017ZX05030-001. **Acknowledgments:** The authors are grateful for the technical and financial support from the Research Institute of Exploration and Development, PetroChina, Beijing 100083, China.

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

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Processes* Editorial Office E-mail: processes@mdpi.com www.mdpi.com/journal/processes

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18