Infill Well Location Optimization Method based on Remaining Oil Recoverable Potential Evaluation

1. Introduction

Suitable infill well location is the basis of realizing efficient development of oil fields. The design of infill well location is based on the determination of infill method, that is, based on the known number of infill wells, the injection-production well type of infill wells and the injection-production well type conversion measures of existing wells, the position of infill wells in a certain infill area, so that the well location can better adapt to the geological characteristics and development conditions of the reservoir [1,2].

The design method of infill well location mainly includes the selection of infill well location and the optimization of infill well location. The infill well location optimization method is to design a variety of infill well location schemes, predict and compare the reservoir development effects of different schemes, and then select the best infill well location scheme [3,4,5,6,7]. The process of determining the infill well location is simple and easy to operate. However, the method of infill well location optimization can only select the best scheme from the limited number of infill well location schemes designed by reservoir engineers, and can not optimize the infill well location, and there is the possibility of not finding the best scheme.

In recent years, the method of infill well location optimization has become a research hotspot in production optimization [8,9,10,11,12]. The method transforms the well location problem into a mathematical problem. With reservoir development effect indicators (cumulative oil production, recovery degree, economic net present value, etc.) as the optimization objective [13,14,15,16,17], and with the characterization parameters of infill well location as the variables, the optimization of infill well location is carried out under certain constraints. The optimal infill well location is affected by many factors such as reservoir heterogeneity, fluid factors, production factors, and economic parameters [18,19,20,21,22]. The relationship between development effect and infill well location is a complex nonlinear, resulting in the infill well location optimization is extremely complex, non-linear and multi-peak [18,22,23,24]. At present, this kind of problem is usually solved by coupling reservoir numerical simulators and optimization algorithms. This method uses the algorithm solution mechanism to update the infill well position automatically. However, the problem of well location optimization is often time-consuming [25,26,27]. Therefore, to find the optimal infilling well location efficiently and accurately, some researchers use the way of production potential map to determine. Production potential mapping [28,29] is a method developed by Cruz and Horne to simplify the optimization of infill well placement. A well is arranged in all or part of the grid of the reservoir, and a reservoir numerical simulator is run to predict the cumulative oil production of the reservoir at each well location, and the kriging interpolation method is used to obtain the cumulative oil production map, that is, the production potential map. When the production potential map is applied to optimize the well location, all the grids are sorted according to the production potential from the largest to the smallest, the grids without drilling are selected in turn, and whether there are other wells in the control area of the selected well location is judged. If there are other wells, the well location is abandoned to continue the selection. Based on the calculation formula of production potential proposed by Liu et al [30], the potential map not only contains static oil and gas reserve parameters but also includes parameters characterizing oil and gas flow capacity. Ding et al. considered the influence of gas cap and bottom water and expanded the application scope of the optimization method of well location [31]. Molina et al. proposed a more commonly used formula for calculating the development potential index [32]. Karim et al. further combined the static and dynamic characteristics and carried out the corresponding development potential evaluation [10,33,34,35]. Geng et al. considering the nonlinear relationship between remaining oil diversion ability and water saturation, introduced dominant potential abundance to quantitatively characterize the distribution of remaining oil potential in high water cut reservoirs [36]. Ding et al. proposed to quantitatively characterize the distribution potential of remaining oil in the high water cut period by using superior production potential, and compared it with the abundance of reserves and superior potential [37].

Therefore, in light of the existing challenges in evaluating remaining oil potential in ultra-high water cut reservoirs, this paper first establishes a quantitative evaluation method based on the quantitative characterization of recoverable reserves of remaining oil, water flooding conditions, and reservoir heterogeneity. It also considers the varying degrees of influence different indicators have on remaining oil potential to achieve a comprehensive evaluation of the distribution of recoverable potential of remaining oil. Based on the evaluation results, the potential area for infill wells is determined. The infill well location optimization method is then used to optimize the heel coordinate, well length, well inclination, and azimuth angle of the infill wells, enabling the precise design of infill well locations.

2. Establishment of Remaining Oil Recoverable Potential Evaluation Method

After long-term water injection and erosion in ultra-high water cut reservoirs, the distribution of underground oil and water becomes complex. Parameters such as reservoir physical properties, pore structure, wettability, and seepage characteristics change significantly. The static and dynamic heterogeneity of reservoirs is further enhanced, resulting in a sharp increase in the difference in oil-water diversion capacity. The variation in water flooding degree in different parts leads to a significant difference in oil-water diversion capacity. The conventional quantitative characterization method for remaining oil potential can no longer effectively reflect this disparity. In ultra-high water cut reservoirs, the distribution of remaining oil is highly dispersed. Its development potential is influenced by multiple factors, including reservoir heterogeneity, the scale of recoverable remaining oil reserves, water flooding status, and oil-water diversion capacity.

Therefore, by thoroughly considering the factors affecting the potential of remaining oil recovery and taking into account the significant reservoir heterogeneity, a quantitative evaluation system is devised. This system comprehensively captures reservoir heterogeneity, the scale of recoverable remaining oil reserves, water flooding conditions, and oil-water diversion capacity. It delves deeply into the analysis of remaining oil recovery potential in ultra-high water cut reservoirs, offering insights for the comprehensive exploitation of remaining oil with varying potential.

2.1. Evaluation Index Construction

The unit conductance coefficient considers the influence of parameters related to remaining oil and reservoir properties on the flow capacity of remaining oil. It is represented by the absolute permeability and crude oil viscosity of grid units within a certain range near the reservoir grid unit, rather than solely representing the attributes of individual grid units. This characterization reflects the ability of grid units to exchange remaining oil with neighboring grid units. A higher unit conductance coefficient indicates a stronger oil exchange capability between grid units and neighboring ones, implying a greater potential for remaining oil recovery. Its calculation formula is as follows:

$$CO{I}_{i,j}=\sqrt{0.5\left[{\left(K/{\mu }_o\right)}_{i-\frac{1}{2},j}+{\left(K/{\mu }_o\right)}_{i+\frac{1}{2},j}\right]\cdot 0.5\left[{\left(K/{\mu }_o\right)}_{i,j-\frac{1}{2}}+{\left(K/{\mu }_o\right)}_{i,j+\frac{1}{2}}\right]}$$

(1)

$${\left(K/{\mu }_o\right)}_{i\pm \frac{1}{2},j}=\frac{{\left(K/{\mu }_o\right)}_{i\pm 1,j}{\left(K/{\mu }_o\right)}_{i,j}}{{\left(K/{\mu }_o\right)}_{i\pm 1,j}+{\left(K/{\mu }_o\right)}_{i,j}}$$

(2)

$${\left(K/{\mu }_o\right)}_{i,j\pm \frac{1}{2}}=\frac{{\left(K/{\mu }_o\right)}_{i,j\pm 1}{\left(K/{\mu }_o\right)}_{i,j}}{{\left(K/{\mu }_o\right)}_{i,j\pm 1}+{\left(K/{\mu }_o\right)}_{i,j}}$$

(3)

where COI represents the unit conductance coefficient; K denotes the reservoir permeability in 10⁻³μm²; μ_o is the oil viscosity, mPa⋅s; subscripts i and j represent the X and Y direction indices of the grid unit, respectively.

The mobility of residual oil is influenced by the current oil saturation and the residual oil saturation of the reservoir, which are affected by the permeability characteristics of the reservoir. Traditional methods for characterizing the potential scale of residual oil primarily use residual oil abundance or recoverable residual oil abundance. Residual oil abundance only reflects the scale of residual oil reserves per unit area, while recoverable residual oil abundance simply considers residual oil saturation, ignoring the diversion ability of mobile oil at different water saturation levels. Residual oil advantage potential abundance, building upon recoverable reserves abundance, further considers the nonlinear relationship between oil-water diversion capability and water saturation during high water saturation periods. This method characterizes both the recoverable reserves scale of underground residual oil and the differences in diversion capabilities of residual oil within the reservoir during high water saturation periods. Higher residual oil advantage potential abundance indicates greater recoverable potential of residual oil.

$$VO{I}_{i,j}=\frac{{\mu }_w}{{\mu }_o}{e}^a{(c-S_w)}^2{e}^{b{S}_{\mathrm{w}}}\left[S_{oi,j}\left(t\right)-S_{or}\right]\cdot {\varphi }_{i,j}\cdot h_{i,j}\cdot NT{G}_{i,j}$$

(4)

where VOI represents the unit remaining oil dominant potential abundance; μ_w is water viscosity, mPa·s; S_w denotes the water saturation; S_o signifies the oil saturation; S_or stands for the residual oil saturation; ø represents the porosity; h indicates the thickness of the grid unit, m; NTG represents the net-to-gross ratio; a, b, and c are the fitting coefficients of the nonlinear phase-permeability relationship during the high-water-cut period.

During the high-water-cut period, the reservoir has been extensively affected by water injection, leading to significant water-flooded areas where the late-stage development potential is limited. Based on Leverett's concept of fractional flow, neglecting the influence of capillary forces and gravity, the formula for calculating the fractional flow of residual oil (i.e., oil saturation) is as follows:

$$f_o=1-f_w=1-\frac{1}{1+\frac{{\mu }_w}{{\mu }_o}\frac{K_{ro}}{K_{rw}}}=1-\frac{1}{1+\frac{{\mu }_w}{{\mu }_o}{e}^a{(c-S_w)}^2{e}^{b{S}_{\mathrm{w}}}}$$

(5)

where f_o represents the oil cut of the unit; f_w represents the water cut of the unit; K_ro represents the relative permeability of oil; K_rw represents the relative permeability of water.

For locations with different physical properties in the reservoir, the relative permeabilities of oil and water vary. Even if the oil saturation remains the same, the oil phase's ability to flow differs. As the oil saturation decreases in a reservoir grid unit, the residual oil's flow rate continuously decreases. Generally, stronger residual oil flow capacity indicates greater remaining oil recovery potential.

2.2. Establishment of Evaluation Model

To characterize the remaining oil recovery potential during the high-water-cut period, the combination of the three evaluation criteria is weighted and multiplied together. This combined metric is defined as the residual oil recovery potential index, expressed mathematically as follows. The advantage of combining these three indices by multiplication is that when any one of them indicates non-recoverable conditions, the residual oil recovery potential index can capture that characteristic. A higher value of the residual oil recovery potential index represents a higher potential for residual oil recovery in the area unit.

$$ROI=VO{I}^{\alpha }\cdot CO{I}^{\beta }{f}_o^{\gamma }$$

(6)

where ROI is the residual oil recovery potential index; α, β, and γ are weighting coefficients used to emphasize the importance of each criterion in the residual oil recovery potential index. A higher weighting coefficient implies greater importance of the corresponding criterion in the residual oil recovery potential index.

To objectively determine the importance of each evaluation index and characterize the residual oil recovery potential during the high-water-cut period, this study introduces a projection pursuit model to project multidimensional evaluation index values into one-dimensional projection data [38]. An accelerated genetic algorithm is employed to optimize the projection index function to seek the optimal projection direction of the evaluation index system [39]. Based on the magnitudes of each component in the optimal projection direction, the objective weights of each evaluation index are determined, effectively avoiding subjective biases from manual assignment and enhancing the objectivity of evaluating the residual oil recovery potential during the high-water-cut period.

Taking the actual reservoir geological model grid as the unit, the spatial distribution of remaining oil is determined and the characteristic parameters are extracted, and the three evaluation index values of the remaining oil recoverable potential evaluation in the ultra-high water cut stage are calculated. Take the logarithm at both ends of Eq.(6):

$$RO{I}^{*}=\lg ROI=\alpha \lg VOI+\beta \lg COI+\gamma \lg f_o$$

(7)

The sample data set x_k={lgVOI_k, lgCOI_k, lgf_ok} composed of m grid units and three evaluation indexes is normalized to eliminate the influence of dimensional difference or numerical range difference of each index on the evaluation results. Since the above three indicators are positive, the normalization formula is as follows:

$$\begin{array}{l}VO{I}_k*=\frac{\lg VO{I}_k-\lg VO{I}_{\min }}{\lg VO{I}_{\max }-\lg VO{I}_{\min }}\\ CO{I}_k*=\frac{\lg CO{I}_k-\lg CO{I}_{\min }}{\lg CO{I}_{\max }-\lg CO{I}_{\min }}\\ f_{ok}*=\frac{\lg f_{ok}-\lg f_{\mathrm{omin}}}{\lg f_{\mathrm{omax}}-\lg f_{\mathrm{omin}}}\end{array}$$

(8)

Assuming that d={d₁, d₂, d₃} is the projection direction vector, the normalized one-dimensional projection value Z_k of the evaluation sample data x_k* in this direction is:

$$\begin{array}{cc}{Z}_k=d_{1}VO{I}_k*+d_{2}CO{I}_k*+d_3{f}_{ok}*& k=1,2,\mathrm{...}m\end{array}$$

(9)

Usually, the optimal projection direction is obtained by maximizing the standard deviation matrix S_Z and the local density matrix D_Z of the projection value Z_k, which preserves as much useful information in the original data as possible. Therefore, the projection index function can be expressed as:

$$Q\left(d\right)=S_Z{D}_Z$$

(10)

where

$$\begin{array}{l}{S}_Z=\sqrt{{\displaystyle \sum _{i=1}^m\frac{{\left(Z_i-\overline{Z}\right)}^2}{m-1}}}\\ D_Z={\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^m\left(R-r_{ij}\right)\times }}u\left(R-r_{ij}\right)\\ r_{ij}=\left|Z_i-Z_j\right|\end{array}$$

(11)

When the sample data set of each evaluation index is given, the projection direction vector d is the only variable that affects the projection index function Q(d). Therefore, when the projection index function takes the maximum value, the corresponding projection direction vector reflects the best projection direction of the data. The optimization problem of the best projection direction is transformed into the problem of solving the maximization of the projection index function, that is:

$$\begin{array}{l}\begin{array}{cc}\max \text{:}& Q\left(a\right)\end{array}=S_Z{D}_Z\\ \begin{array}{cc}\mathrm{s}.\mathrm{t}.& {\displaystyle \sum _{j=1}^3{a}_j^2=1}\end{array}\end{array}$$

(12)

The optimization problem is a typical complex nonlinear problem, which is solved by an accelerated genetic algorithm based on real number coding [40]. The standard deviation matrix S_Z and the local density matrix D_Z in the projection pursuit model are used as the initial population of the accelerated genetic algorithm. The accelerated genetic algorithm is used to perform global optimization in the full parameter space until the optimal solution that maximizes the projection index function is found. That is the optimal projection direction of the projection pursuit model. According to the projection pursuit model, each component of the optimal projection direction reflects the degree of influence of each evaluation index on the recoverable potential of the remaining oil in ultra-high water cut reservoirs.

3. Construction of Infill Well Location Optimization Method

The precise optimization of infill well placement involves identifying the optimal parameters, such as toe coordinates, wellbore length, inclination angle, azimuth, etc., within regions of identified infill potential. This optimization is conducted based on the assessment results of residual oil potential, considering factors such as the number of infill wells, the assignment of injector-producer pairs, and the conversion of existing well assignments. The objective is to enhance reservoir development effectiveness by strategically locating infill wells in areas with high infill potential.

3.1. Construction of Optimization Mathematical Model

3.1.1. Optimization Variable

The optimization problem of infill well placement aims to find the optimal locations for these wells, with characteristic parameters representing the variables to be optimized. Assuming the wellbore is a straight line segment, six variables (x_h, y_h, z_h, L, θ, ϕ) are used to represent the position of each well, where x_h, y_h, z_h are the coordinates of the wellhead, L is the length of the wellbore, θ is the azimuth angle, and ϕ is the inclination angle. The specific parameters are illustrated in Figure 1. For inclined wells, each well can be represented by six variables (x_h, y_h, z_h, L, θ, ϕ). For vertical wells, each well can be represented by four variables (x_h, y_h, z_h, L), and for horizontal wells, each well can be represented by five variables (x_h, y_h, z_h, L, θ).

3.1.2. Objective Function

Performance indicators represent the quantities to be optimized in the optimization problem. Common performance indicators for infill well optimization include cumulative oil production, oil recovery factor, and net present value (NPV). The objective function is the mathematical expression of the performance indicators. In this study, we use the maximization of NPV as the objective function for infill well optimization. NPV comprehensively evaluates the economic factors of reservoir production during the development process, including production revenue, injection cost, water treatment cost, and the impact of horizontal section drilling costs on development effectiveness. Additionally, NPV, as an optimization objective function, has been widely applied in optimization problems throughout the oil and gas field development process.

$$\max NPV={\displaystyle \sum _{t=1}^T\left[\frac{1}{{\left(1+b\right)}^t}\left({\displaystyle \sum _{p=1}^{N_p}{C}_o{Q}_o^{\left(t\right)}-{\displaystyle \sum _{p=1}^{N_p}{C}_w{Q}_w^{\left(t\right)}-{\displaystyle \sum _{I=1}^{N_I}{C}_I{Q}_I^{\left(t\right)}}}}\right)\right]}-N_{\inf }{C}_{d}L$$

(13)

where T represents the total production time in years; t represents the specific year; C_o is the oil price, USD/m³; C_w is water production cost, USD /m³; C_I is water injection cost, USD /m³; C_d is the drilling cost, USD /m; Q_o is the oil production rate, m³/d; Q_w is the water production rate, m³/d; Q_I is the injection rate, m³/d; N_P is the number of production wells; N_I is the number of injection wells; N_inf is the number of infill wells; L is the wellbore length, m; b is the discount rate.

3.1.3. Constraint Condition

To ensure that the optimized infill well locations fall within feasible ranges, constraint conditions need to be incorporated into the optimization process. The constraints for infill well placement optimization include feasible infill range constraints, minimum well spacing constraints, wellbore length constraints, azimuth angle range constraints, and inclination angle range constraints.

• Feasible infill range constraints

The feasible infill range constraints are incorporated into the optimization of infill wells, ensuring that the root and toe coordinates of the optimized wells do not exceed the feasible infill area. The upper and lower limits of the infill range are represented by x_max, y_max, z_max, and x_min, y_min, z_min. The constraints that the root coordinates (x_h, y_h, z_h) must satisfy are as follows:

$$\begin{array}{l}{x}_{\min }\le x_{\mathrm{h}}\le x_{\max }\\ y_{\min }\le y_{\mathrm{h}}\le y_{\max }\\ z_{\min }\le z_{\mathrm{h}}\le z_{\max }\end{array}$$

(14)

The toe coordinates (x_t, y_t, z_t) can be calculated based on the representative parameters of the well position:

$$\begin{array}{l}{x}_{\mathrm{t}}=x_{\mathrm{h}}+Lcos\theta \sin \phi \\ y_{\mathrm{t}}=y_{\mathrm{h}}+L\sin \theta \sin \phi \\ z_{\mathrm{t}}=z_{\mathrm{h}}+Lcos\phi \end{array}$$

(15)

The toe of the infill well needs to be within the feasible infill area. The constraints for the toe coordinates are as follows:

$$\begin{array}{l}{x}_{\min }\le x_{\mathrm{t}}\le x_{\max }\\ y_{\min }\le y_{\mathrm{t}}\le y_{\max }\\ z_{\min }\le z_{\mathrm{t}}\le z_{\max }\end{array}$$

(16)

• Minimum well spacing constraint

The well spacing refers to the distance between two wellbores. When the well spacing is too small, it can lead to interference between wells. To prevent this from happening, a minimum distance between wells should be set during the optimization of infill well locations in the well pattern. If i and j represent two wells, M is a point on the straight line of well i's borehole, and N is a point on the straight line of well j's borehole, then the distance between points M and N can be denoted by D_i,j.

$${\mathrm{D}}_{i,j}\ge {\mathrm{D}}_{\min },i,j\in N_{\mathrm{T}}i\ne j$$

(17)

where N_T is the total number of wells in the reservoir.

• Well length constraint

The wellbore length needs to be within a certain range. When the minimum wellbore length is L_min and the maximum length is L_max, the wellbore length constraint satisfies a linear constraint condition, which can be expressed as:

$$L_{\min }\le L\le L_{max}$$

(18)

• Orientation angle range constraint

When optimizing for deviated or horizontal wells, the orientation angle of the infill wells must fall within a feasible range, influenced by actual production conditions and technical constraints. When the minimum orientation angle is θ_min and the maximum angle is θ_max, the constraint on the orientation angle falls within linear constraints.

$${\theta }_{\min }\le \theta \le {\theta }_{max}$$

(19)

The wellbore inclination angle, representing the angle between the wellbore and the vertical line, is influenced by drilling techniques and geological factors. When optimizing for deviated wells, the inclination angle must fall within a feasible range. With a minimum angle φ_min and maximum angle φ_max, the inclination angle range constraint conforms to linear constraints.

$${\phi }_{\min }\le \phi \le {\phi }_{max}$$

(20)

Combining the analysis above, based on the determination of infill potential areas, maximizing the net present value (NPV) is set as the objective function. The optimization variables consist of a matrix of well network characteristic parameters (well location, wellbore length, well deviation angle, azimuth angle). Together with constraints including the infill range, wellbore length, azimuth angle, and minimum well spacing, the mathematical model for optimizing infill well locations is established to achieve the precise well design. The mathematical model is represented by Eq.(21).

$$\begin{array}{l}\max NPV={\displaystyle \sum _{t=1}^T\left[\frac{1}{{\left(1+b\right)}^t}\left({\displaystyle \sum _{p=1}^{N_p}{C}_o{Q}_o^{\left(t\right)}-{\displaystyle \sum _{p=1}^{N_p}{C}_w{Q}_w^{\left(t\right)}-{\displaystyle \sum _{I=1}^{N_I}{C}_I{Q}_I^{\left(t\right)}}}}\right)\right]}-N_{\inf }{C}_{d}L\\ \mathrm{s}.\mathrm{t}.x_{\min }\le x_{\mathrm{h}}\le x_{\max }\\ y_{\min }\le y_{\mathrm{h}}\le y_{\max }\\ z_{\min }\le z_{\mathrm{h}}\le z_{\max }\\ x_{\min }\le x_{\mathrm{h}}+Lcos\theta \sin \phi \le x_{\max }\\ y_{\min }\le y_{\mathrm{h}}+L\sin \theta \sin \phi \le y_{\max }\\ z_{\min }\le z_{\mathrm{h}}+Lcos\phi \le z_{\max }\\ {\mathrm{D}}_{i,j}\ge {\mathrm{D}}_{\min },i,j\in N_{\mathrm{T}}(i\ne j)\\ L_{\min i}\le L_i\le L_{maxi}\\ {\theta }_{\min i}\le {\theta }_i\le {\theta }_{maxi}\\ {\phi }_{\min i}\le {\phi }_i\le {\phi }_{maxi}\end{array}$$

(21)

3.2. Solution of optimization Mathematical Model

After water injection, oil reservoirs often exhibit increased heterogeneity in parameters such as porosity, permeability, and oil saturation. Within feasible densification ranges, different densification well locations result in varying reservoir production parameters and consequently different reservoir development outcomes. Therefore, the optimization of densification well locations involves multiple factors such as reservoir heterogeneity, fluid heterogeneity, and production systems, making it a complex optimization problem.

The optimization problem of densification well locations exhibits the following characteristics: (1) It involves a large number of optimization variables, resulting in a high-dimensional optimization problem. For oil reservoirs of varying scales, the number of optimization wells can range from a few to hundreds, corresponding to several hundred optimization variables. The higher the number of optimization variables, the more challenging the optimization problem becomes to solve. (2) The solution space is rugged, with multiple peaks in the optimization problem. As the complexity of the optimization problem increases, the ruggedness of the solution space also increases. When solving optimization problems with multiple peaks, the rougher the solution space, the more prone it is to getting trapped in local optima, making the solution process more complex. (3) The optimization problem is nonlinear, making it difficult to obtain derivative information. In the optimization of densification well locations, parameters such as permeability, porosity, oil saturation, reservoir fluids, well spacing, and displacement pressure gradients within injection and production units vary at different well locations. This leads to different reservoir development outcomes at different well locations, and the cumulative oil production of the reservoir exhibits a complex nonlinear relationship with reservoir and production parameters, making it difficult to obtain derivatives. Based on the analysis above, obtaining derivative information during the solution process of the optimization problem of densification well locations is challenging. Therefore, gradient-free optimization algorithms are used to solve this problem.

PSO-MADS is a typical hybrid optimization algorithm that couples particle swarm optimization (PSO) [41,42] algorithm with mesh adaptive direct search (MADS) [43] algorithm, belonging to the class of global-local stochastic search algorithms. It combines the global optimization capability of PSO with the local optimization capability of MADS. The PSO-MADS algorithm benefits from the global search provided by PSO and the local search provided by MADS. The algorithm switches between PSO and MADS based on specific conditions and settings.

The PSO-MADS algorithm begins with an initial particle swarm, which includes one or more custom initial values, and then applies PSO for iterations. If the global best position improves, indicating a successful PSO search step where a new global best position with an objective function better than the previous one is generated, PSO continues to run. When the search step is not successful (i.e., no improvement is provided), the algorithm proceeds to the MADS polling step. In this step, the best position of the cluster serves as the center, and polling is conducted around the best solution found by PSO. As long as the polling step is successful, yielding better feasible points regarding the objective function, the polling continues. The MADS polling step runs until an unsuccessful iteration occurs, at which point the grid size parameter decreases, and the hybrid algorithm returns to PSO search. PSO search continues around the best solution generated by the current MADS, repeating the above steps. This alternation between PSO and MADS continues until the optimal solution is found and the search terminates.

4. Field Application

4.1. Reservoir Model Description

An oilfield in Bohai is a large low-amplitude drape structure developed on the background of a Paleogene paleo-uplift. Due to the multiple constraints of structure, fault, and lithology, the oil-water system is complex, the reservoir types are diverse, and lithologic structural reservoirs and structural reservoirs are developed. According to the internal fault segmentation and sand body distribution range, the oilfield is divided into three blocks in the north, south, and west. The southern area is located in the southeast of the oilfield, and the structure is generally gentle. The NmIL sand body is an important main production layer in the southern area. The drilling data and 3-D seismic data reveal that the reservoir plane changes rapidly and the vertical thickness changes greatly, showing the sedimentary characteristics of typical multi-stage channel sand body superposition. The reservoir type is a lithologic structural edge water reservoir. From the perspective of reservoir physical properties, the average porosity of the reservoir in the southern area is 32 %, and the average permeability is 3689 mD. It is a high-porosity and high-permeability reservoir with a large reservoir thickness, generally around 9 m. The current injection-production well pattern in the southern area of the oilfield is relatively perfect. There are 62 production wells and 29 injection wells in the NmIL sand body. The oil recovery is 31.5 % and the comprehensive water cut is 95.2 %. It has entered the development stage of the ultra-high water cut period (Figure 3).

The distribution of remaining oil in the plane of the NmIL sand body is quite different, showing the characteristics of overall dispersion and local enrichment. The remaining oil in the central well area is more dispersed, the remaining oil in the edge is locally enriched, and the remaining oil is difficult to tap. Therefore, according to the evaluation method of the recoverable potential of remaining oil, this section first clarifies the area with encryption potential, and on this basis, the precise optimization of the infill well location is carried out.

4.2. Determination of Encryption Potential Area

Taking the NmIL sand body of the main production layer in the southern area of the oilfield as the research object, the evaluation of the remaining oil exploitation potential of the ultra-high water cut reservoir is carried out. On the basis of reservoir numerical simulation, the grid units divided by the reservoir numerical simulation model are taken as the research object, and the oil saturation distribution data and related attribute parameters are extracted. According to Eqs.(1), (4), and (5), the conductivity coefficient, the abundance of remaining oil dominant potential and the oil cut of each grid unit are calculated respectively, as shown in Figure 4. It can be seen from Figure 4 that for the locations with different physical conditions in the reservoir, under the condition that the regional remaining oil saturation is the same, the remaining oil recoverable potential evaluated by the three evaluation indexes is quite different. For example, for the low well pattern control area at the edge of the sand body, the remaining oil saturation is generally high, but the local unit conduction coefficient and the remaining oil dominant potential abundance distribution are different. Even due to the regional difference in oil phase diversion ability, there is a big difference in the oil content distribution of local units. Therefore, there will be some contradictions or problems when only a single index is used to describe the remaining oil recoverable potential. It is difficult to effectively guide the potential tapping of remaining oil in an ultra-high water cut period.

Based on the calculation of three evaluation index values of each grid unit of the NmIL sand body, three evaluation index values of each grid unit are taken as a set of sample data, and 1000 grid units are randomly selected from all grid units to form an evaluation index sample data set. The sample dimension is 1000, and the number of evaluation indexes is 3. After the index data is normalized, the projection index function is optimized by an accelerated genetic algorithm based on real number coding. The specific algorithm parameters are set as follows: the population size is 400, the random number required for the mutation direction is 20, the crossover probability is 0.8, the mutation probability is 0.2, the number of accelerations is 7, and the number of accelerations after two generations of evolution is 2. Then the optimal projection direction vector d = (0.5106, 0.0538, 0.8581) of the remaining oil recoverable potential evaluation index system is obtained.

According to the projection pursuit model, the size of each component of the optimal projection direction essentially reflects the influence degree of each index on the evaluation of the remaining oil recoverable potential of the ultra-high water cut reservoir, that is, the weight. The order from large to small is the unit oil content, the unit remaining oil advantage potential abundance, and the unit conduction coefficient, of which the first two indicators are the main influencing factors. After normalization, the corresponding weights of each evaluation index are obtained: α=0.35895, β=0.03782, γ=0.60323. Substituted into Eq.(6), the expression of the remaining oil recoverable potential index is obtained:

$$ROI=VO{I}^{0.35895}\cdot CO{I}^{0.03782}\cdot f_o^{0.60323}$$

(22)

The remaining oil recoverable potential index of each grid unit is calculated by Eq. (22), which can comprehensively characterize the influence of reservoir physical properties, recoverable reserves abundance, and oil-water diversion capacity in different regional positions on the remaining oil production potential in the ultra-high water cut stage. Then the remaining oil recoverable potential map of the NmIL sand body in the ultra-high water cut stage is drawn, as shown in Figure 5.

It can be seen from Figure 5 that when the remaining oil recoverable potential map is used to describe the remaining oil potential in the ultra-high water cut period, the remaining oil recoverable potential in different regions is different, and the dominant potential areas are prominent, which can more accurately and scientifically guide the adjustment of injection-production well pattern. (a) For the eastern part of the sand body adjacent to the edge water area, the reservoir permeability is generally high, and the local unit conductivity coefficient is high, resulting in serious water flooding, and the remaining oil advantage potential abundance is low. Therefore, the remaining oil recoverable potential in this area is low. (b) Similarly, for the high remaining oil saturation area in the northwest of the sand body, the unit conductivity coefficient and the abundance of remaining oil dominant potential are generally low, resulting in the low recoverable potential of the remaining oil in the region, and there are local remaining oil high recoverable potential areas, which can be used as the next infill well adjustment potential area. (c) The central well area has a high degree of well pattern control, and the injection production is relatively perfect. There is only a high recoverable potential area between local wells, and the potential for infill adjustment is low. (d) At present, the remaining oil high recoverable potential area of the main layer is generally located at the edge of the sand body, especially in the southwest and southeast of the sand body, the permeability is high, the oil phase flow ability is strong, the remaining oil recoverable potential is more prominent, and the well pattern control degree is low, which can be given priority as the next infill well adjustment potential area.

4.3. Infill Well Location Optimization

Based on the analysis of the remaining oil recoverable potential of the main production layer NmIL sand body in the southern area of the oilfield, the infill well location optimization method is applied to optimize the infill well location of the target block. In this optimization problem, the optimization variables are the coordinates of the heel end, the length of the well, the deviation angle, and the azimuth angle of each infill well. When conducting the optimization of the infill well location, the optimized infill well location should not change the well pattern infill mode.

The objective function is to seek the maximum NPV and related economic parameters are listed in Table 1. The well spacing between any two wells needs to be greater than 100 m, the well length does not exceed 250 m, the well deviation angle needs to be less than 50°, and the infill well range constraint is determined by the remaining oil potential analysis in the previous section. Therefore, the constraint conditions of the optimization problem include feasible infill range constraint, minimum well distance constraint, well length constraint, and well inclination constraint. The heel end coordinate constraint, well length constraint, and well inclination constraint are treated by the truncation method, the toe end coordinate constraint is treated by the direct rejection method, and the minimum well distance constraint is treated by the adaptive penalty function method. The PSO-MADS algorithm is used to solve the problem. The number of iterations is 50, the number of populations is 100, and the optimization is performed independently 10 times.

In order to analyze the difference between the remaining oil potential based infill well optimization method and the conventional infill well optimization method in development effect, this section adopts the conventional infill well optimization method and the new method respectively. The conventional infill well optimization method means that the infill well location is optimized throughout the reservoir. Based on the analysis of the remaining oil potential, the new method determines the region with encryption potential and then optimizes the region. As shown in Figure 6, the average iteration curve obtained by the two methods is that the initial infill well location of the conventional method is randomly given within the reservoir range, while the initial infill well location of the new method is given within the range of potential for infill. Therefore, the NPV obtained by the initial scheme of the two methods is different, and the NPV obtained by the new method is higher than that of the traditional method.

The relationship curve of water cut and oil recovery before and after optimization is shown in Figure 7. It can be seen that when the same water cut is reached, the oil recovery after optimization is higher than that before optimization, and the development effect of the target block is improved. Moreover, the reservoir development effect can be improved greatly by using the infill well location optimization method based on remaining oil potential

Table 2 shows the well location parameters of the optimized infill well of the new method. Figure 8 shows the remaining oil saturation distribution and streamline field distribution before and after optimization. It can be seen from the Figure that the streamline of the remaining oil enrichment area after optimization of the infilled well increases, and the efficiency direction of the single well increases, so that the remaining oil that cannot be displaced by the original streamline can be effectively used, and the displacement of crude oil is more balanced.

5. Conclusions

This study comprehensively addresses the impact of reservoir physical properties, recoverable reserve abundance, and oil-water diversion capacity on the production potential of remaining oil in ultra-high water cut reservoirs. It establishes evaluation indices including unit conductivity coefficient, unit remaining oil dominant potential abundance, and unit oil content index to assess the recoverable potential of remaining oil. By integrating the accelerated genetic algorithm and projection pursuit model, we determine the objective weight of each evaluation index and develop a quantitative evaluation method for remaining oil recoverable potential in ultra-high water cut reservoirs.

Building upon this foundation, we formulate a mathematical model for optimizing the well location of infill wells. This model optimizes the well location coordinates, well length, well inclination angle, and azimuth angle of infill wells. Field application demonstrates that relying solely on a single index, such as saturation, to describe the recoverable potential of remaining oil may lead to contradictions or issues. This approach often struggles to effectively guide the potential tapping of remaining oil during the ultra-high water cut period. In contrast, the remaining oil recoverable potential index enables the evaluation of remaining oil recoverable potential in different regions, offering clear guidance for infill well location optimization.