#### Abstract

Well placement optimization is a significant task in oil field development which aims to find the optimal well locations by maximizing the net present value (NPV) or other objective function. It is a highly nonlinear problem involving large number of variables. Despite lots of work has been done on conventional reservoirs, the optimization tool for naturally fractured reservoir (NFR) is still rare. Naturally fractured reservoirs represent significant amounts of oil reserves. The well placement optimization for NFR is challenging due to the high heterogeneity of matrix and fracture system. In this work, a computer-aided well placement optimization method is established for NFR based on the recent advances. The two phases flow discrete fracture numerical simulation model, i.e., the embedded discrete fracture method (EDFM) is used to model the natural fractures as its computational efficiency and flexibility to handle fracture. Then, stochastic simplex approximate gradient (StoSAG) is employed to obtain the approximate gradient by combing the EDFM. The steepest ascent algorithm is used to find the optimal well placement. A series of numerical case studies are set up to examine the performance of the proposed approach. The NPV for water flooding naturally fractured reservoir production optimization substantially increased by using StoSAG.

#### 1. Introduction

Determining the location of wells is crucial during field development process because it can affect the final NPV significantly [1]. Lots of automated well placement optimization methods has been investigated in previous studies [2–5]. However, little has suggested an effective way to assess the well location optimization for naturally fractured reservoir. In recent years, naturally fractured reservoirs receive great attention as its significant amounts of oil reserves [6]. In order to improve the efficiency of reservoir development and enhance oil recovery, well location for NFR should be carefully arranged and optimized.

Well placement optimization is a highly nonlinear problem involving the reservoir response. Reservoir simulation simulators are the common tools to achieve the rates of oil and water. In the optimization process, the reservoir simulation model may require thousands of runs. Thus, the simulator should be computational efficiency. For NFR, simulators are developed based on two classical types of fracture models: dual-continuum method and explicit discrete-fracture and matrix model. Dual-continuum method is widely used in most commercial reservoir simulators such as Eclipse and CMG. Dual-continuum method divides the reservoir domain into fracture and matrix [7–9]. Commonly, it can keep applicability when the fractures are denser and the representative elementary volume (REV) exists [10]. It also loses accuracy in flow calculation when a number of large fractures locates in the reservoir [11]. The fracture and matrix system are coupled by transfer functions. Fracture-matrix flow is controlled mainly by matrix propertied, and the shape factor needs to be determined. Unfortunately, the shape factor is difficult to calculate when considering capillarity, gravity [12]. Explicit discrete-fracture and matrix model or discrete fracture-matrix (DFM) has grown in popularity during recent years. The model deals with every fracture explicitly. Thus, it can capture the fracture geometries and accurately characterize the flow exchange between fracture and matrix. Many DFMs have been developed [13–18]. The unstructured grid is always chosen to discrete the calculation domain. The grid size near the fractures should keep small to exactly simulate the flow between matrix and fracture. However, the local grid refinement leads to large computational load. A recently popular DFM—embedded discrete fracture model (EDFM)—attracts much attention and shows some advantages. The EDFM is firstly proposed as a hierarchical modeling method to deal with multiple length scales in naturally fractured formations [19–21]. Then, it is extended to the flow performance analysis of both naturally fracture reservoirs and hydraulic fractured tight oil reservoirs [22–24]. Some applications of EDFM have been reported. Yao et al. coupled the EDFM and dual-continuum method to inverse multiscale fractures hierarchically using dynamic production data [25, 26]. Yao et al. [27] optimized the fracturing parameters in shale gas reservoir. Yu et al. [28] simulated the pressure response of well interference in tight oil reservoirs with complex–fracture geometries. Alessio et al. [29] employed the EDFM for computing fracture-fracture and matrix-fracture transmissibilities as an upscaling tool. In this study, the EDFM will be used as the simulating tool.

Another important issue for well placement optimization is the chosen of optimization algorithm. Two kinds of optimization algorithms are implemented to find the optimal well placement: (1) gradient-based methods and (2) gradient-free methods. Gradient-based methods mainly include SPSA, finite difference, adjoint method, and descent method. These gradient-based techniques are easy to trap in local optima, and the gradient is difficult to calculate. On the other hand, derivative-free techniques do not require the calculation of derivatives, and they can achieve global search. Many gradient-free methods are used like particle swarm optimization (PSO), genetic algorithm (GA), differential evolution (DE), and covariance matrix adaptation evolution strategy (CMA-ES). Ensemble-based optimization increases popularity recently due to its ability to capture uncertainty of multiple reservoir realizations [30]. Ensemble-based optimization (EnOpt) was first introduced by Lorentzen et al. [31] and Nwaozo [32]. Then, this method is greatly improved and used in reservoir development field [33–35]. In 2017, Fonseca et al. [36] found that not all cases can get the optimal value in the process of robust optimization using this method. Based on this observation, a stochastic simplex approximate gradient (StoSAG) method was proposed. StoSAG improves the EnOpt gradient formula in two aspects, using the initial variable and initial function value to replace the average value of random perturbation position and corresponding function value, respectively. After that, the method is widely used in the field of well location optimization and injection production optimization [37, 38]. In 2019, the researchers of Alamos National Laboratory of the United States proved the advantages of StoSAG in robust optimization theoretically, and the optimal value obtained by StoSAG was significantly higher than that obtained by EnOpt, which provided a strict theoretical basis for the popularization and application of StoSAG. Four kinds of stochastic gradient calculation criteria, which are StoSAG, f- StoSAG, sf-StoSAG, and ss-StoSAG, are proposed. The results show that the optimal injection production control variables can be obtained by using four kinds of gradients. The NPV calculated is greatly improved than EnOpt method [39].

Here, we developed the well placement optimization tool for naturally fractured reservoirs by coupling the EDFM and StoSAG. As far as we know, this is the first well placement optimization tool by introducing the EDFM method. It has broad application prospects. In the tool, the classical EDFM of two-phase flow problem is adopted to simulate the naturally fractures. The detailed geological characteristics of each fracture is kept. The pressure and saturation are solved by Newton-Raphson iteration by carefully setting the time steps to guarantee the convergence. The standard StoSAG is chosen as the optimization algorithm to search for the optimization well placements. The well can connect with the fracture and matrix domain freely. A series of example are test from simple to complex to show the validation of the workflow. Specially, the robust optimization results are presented. This is the first tool to optimize the well placement by coupling the StoSAG and EDFM. The paper is organized as follows. In Section 2, the theoretical model for well placement optimization is presented. Then, we show the optimization model and results. Finally, discussion is given.

#### 2. Methodology

##### 2.1. Numerical Simulation Model

In order to perform the well placement optimization for naturally fractured reservoir, the numerical simulation model should be carefully prepared. In this work, the embedded discrete fracture method (EDFM) is adopted as the numerical simulation tool for well placement optimization. Figure 1 shows four kinds of connections when using EDFM. The three kinds of NNCs are the fracture-matrix connection in the same matrix grid, the fracture-fracture connection in the same matrix grid, and fracture-fracture connection in different matrix grids. By defining three kinds of NNC in preprocessing program, the in-house numerical simulation code can be called to perform the calculation. Previous research results show that pressure distribution, saturation distribution, and the well flow response agree with each other for EDFM, DFM, and LGR. Easy implementation, general applicability, and high computational efficiency are also demonstrated compared to DFM and LGR.

##### 2.2. Optimization Algorithm

Ensemble-based methods show advantages in gradient-based optimization. The motivation is that the real gradient is not always available. The general formulation for StoSAG search direction is given by the following equation [39, 40]: where is the geological model realization number to describe the reservoir uncertainty. contains the placement information for all wells. is the stochastic approximation of the simplex gradient. is obtained by

Then, the gradient becomes

The superscript sign “+” on a matrix denotes the Moore-Penrose pseudo-inverse, where

represents the number of control perturbations. Each control perturbation , at iteration is generated from the distribution ; is a predefined covariance matrix.

##### 2.3. Well Placement Optimization Tool

The objective function commonly used in well placement optimization problem is the NPV, which is defined as where is a -dimensional column vector which contains all well placement information; denotes the time step for the reservoir simulation; is the total number of time steps; the time at the end of the time step is denoted by ; is the time step size; is the annual discount rate; and denote the number of producers and injectors, respectively; is the oil revenue, in $/STB; and denote the disposal cost of produced water and the cost of water injection in units of $/STB, respectively; and , respectively, denote the average oil production rate and the average water production rate for the time step, in units of STB/day; and denote the average water injection rate at the injector for the time step, in units of STB/day.

To account for geological uncertainty, robust optimization is performed. The problem is to maximize the expectation for life cycle NPV which is approximated by the average NPV over geological realizations:

We consider only bound constraints and let and denote the lower bound and upper bound for the well placement variable, respectively. Then, the problem can be expressed as

The logarithm transformation to each element of the control vector is used to search the solution of the well placement optimization problem. The steepest descent optimization algorithm is used in which the new search position is updated as for until convergence, where is the initial guess and is the estimate of the optimal control parameter at the iteration; is the step size.

##### 2.4. Workflow

The workflow of well placement optimization using StoSAG and EDFM is shown in Figure 2. Firstly, geological realizations should be generated. In robust optimization, is commonly set to be bigger than 1. The initial well placement is used to calculate the initial objective function value . In order to compute the objective function, the EDFM simulator is called. Then, the iteration step is performed. For a certain iteration *k*, the stochastic simplex approximate gradient is calculated using the search direction. In our work, the ensemble size is set to be 10 for all cases. Figure 3 shows the reservoir model with three wells. The number of grids is 100 for the reservoir without fracture, while the number of grids is 119 for the reservoir with fracture. Note that if there is no fracture in the reservoir, the well placement is located at the center of the matrix grid. However, when fractures exist, the well placements can locate at the centers for both matrix and fracture grids.

**(a) Reservoir without fracture**

**(b) Reservoir with fractures**

#### 3. Case Study

In this section, some synthetic cases are presented to test the new workflow. The examples are set from simply to complex. The 2D models are firstly discussed. Then, the 3D model is presented. A number of wells and fractures also increase for different examples which can test the performance of the new workflow.

##### 3.1. Example 1: 2D Model with Inclined Single Fracture

A water flooding example is considered which is a 2D heterogeneous model. The model size is 400 m ×400 m ×10 m with 20 ×20 × 1 uniform grid. The size for each grid is 20 m ×20 m ×10 m. The horizontal log-permeability distribution is presented in Figure 4. the fracture width is set to be 0.001 m. The matrix porosity is homogenous and equal to 0.1. The initial pressure is set to be 30 MPa. The reservoir lifetime is 2000 days. One production wells and one injection wells are placed in the reservoir. The production well is operated at fixed bottom hole pressure of 10 MPa, and the injection wells are operated with BHP of 40 MPa. To optimize the NPV, the oil price is set equal to USD 60/stb; the water injection cost is USD 5/stb; the cost of disposing produced water is USD 5/stb; the annual discount rate is 0.1. The maximum number of step size cuts is set to be 5. The total maximum allowable iteration is 50.

The performance of the optimization algorithm is always dependent on the setting parameters. In the optimization process, different value of perturbation size and initial step size are taken to examine their effect on the objective function value and optimal placement. Six test cases are conducted, and the optimization results are also presented in Table 1(case 1 and case 4 have the same setting parameters). Figure 5 shows the NPVs after 50 iterations. It can be seen that all cases converge to a steady NPV value after a series of exploring. Compared the initial well placement and the final solution, the NPV increases substantially. At the initial iterations, the NPV increases rapidly, and then the NPV increases slowly. Lastly, the NPV trends to be a constant for different cases. Figure 6 shows the final converged NPV. The highest NPV is for case 5, and the lowest NPV is for case 2. Despite that the final converged NPV value is different, the difference is very small. The NPV evaluation curves follow similar paths which start from small to large value monotonously. Taking a closer look at the six curves calculated from different cases in Figure 5, we can observe that the red curve shows a relatively high converged performance than other for curves. In the initial iteration stage, case 5 has a rapid search efficiency and converges to the highest NPV finally. Though case 6 has a slow search efficiency during the first 25 iterations, it finds a relatively high local optimum. Observing the curves of case 1, case 2, and case 3, the initial stage search efficiency becomes higher if we set a smaller perturbation size. On the other hand, by comparing case 4, case 5, and case 6, we can see that the initial step size also has great effect for the search path. Overall, the results demonstrate that StoSAG generates optimal well placement, so the stochastic simplex approximate gradient can be used in well placement optimization problem. Figure 7 shows the optimal well placement overlapped on water saturation field in 2000 days. The injection well moves along the +*x* and -*y* direction, and the production well moves along -*x* and + *y* from the initial position. The line connecting two wells trends to be perpendicular to the fracture. The distance between two wells is very close for different cases.

##### 3.2. Example 2: 2D Model with Vertical Single Fracture

A water flooding example is considered shown in Figure 8. The basic model parameter is the same as example 1 except the fracture. A fracture is located at the center of the reservoir. The length of the fracture is 200 m. The fracture width is set to be 0.001 m. Six test cases are conducted, and the optimization results are presented in Table 2. Figure 9 shows the NPVs after 50 iterations. It can be seen that the first five cases are converged to a steady NPV value after 50 iteration. Case 6 is trapped to a local optimum. In this example, case 1 demonstrated an extraordinary ability in research efficiency. For most of the iteration, its NPV is higher than others. Figure 10 is the maximum NPV after 50 iterations. The highest NPV is for case 1, and the lowest NPV is for case 6. Figure 11 shows the optimal well placement overlapped on water saturation field in 2000 days. The injection well moves along the *x* direction, and the production well moves along -*x* from the initial position. The line connecting two wells trends to be perpendicular to the fracture for the first five cases. The distance between two wells is very close for the first five cases.

##### 3.3. Example 3: 2D Model with Multiple Fractures

In example 3, we set 12 fractures in the reservoir. The position is shown in Figure 12. The orientation of the fractures keeps consistent. The other parameters are the same as example 1. Three injection wells and production wells are arranged in the reservoir. Note that in this model, the number of whole optimization variables is 12 considering the well location coordinates in (*x*, *y*)-plane. Six test cases are conducted, and the optimization results are presented in Table 3. Figure 13 shows the NPVs after 50 iterations. It can be seen that the first five cases are converged to a steady NPV value after 50 iteration. In this example, case 6 demonstrated an extraordinary ability in search efficiency. For most of the iteration, its NPV is higher than others. Figure 14 is the maximum NPV after 50 iterations. The highest NPV is for case 6 and the lowest NPV is for case 1. Figure 15 shows the optimal well placement overlapped on water saturation field in 2000 days. The three production wells are located along the diagonal of the reservoir. Two production wells are located at the zone where the permeability is relatively high.

##### 3.4. Example 4: Robust Optimization

In example 4, we consider the robust production optimization. Here, 10 reservoir realizations are randomly chosen, which is used to characterize the reservoir geological uncertainty. Figure 16 shows the log-permeability distribution in the horizontal direction of 10 reservoir models. Like example 3, we set 12 fractures in the reservoir. Three injection wells and production wells are arranged in the reservoir. Table 4 shows the optimal well location. Figure 17 shows the NPVs after 50 iterations. Figure 18 is the maximum NPV after 50 iterations. It can be seen that the highest NPV is for case 4 and the lowest NPV is for case 3. Compared with the results of example 3, the uncertainty decreases the NPV greatly. On the other hand, four types of search methods (f-StoSAG, sf-StoSAG, StoSAG, and ss-StoSAG) for the steepest ascent optimization algorithm are used to optimize the well placement. Figure 19 presents the expected NPV of different methods. As expected, the average NPV generated from different search methods is higher than the initial average NPV. The f-StoSAG obtains a relatively high average NPV.

##### 3.5. Example 5: Modified Egg Model

The egg model has been widely used for well placement and control optimization. The geological parameters, fluids parameters, and production control parameters can be found in Jansen et al. (2014) [41]. The number of gridblocks is 25200 for which (60,60,7) is used to discretize the reservoir in *x*, *y*, and *z* directions, respectively. In this study, all grids are set to be active and will be considered in the simulation runs. The grid block size is set to 8 m ×8 m × 4 m. There are eight injection wells and four production wells. Because the model has no aquifer and no gas cap, primary production is almost negligible. The production wells are operated at fixed bottom hole pressure (BHP) with 39.5 MPa, and the injection wells are operated under a BHP constraint of 42 MPa. Total production time is 7200 days. We seek to optimize the well locations of 8 injection wells and 4 production wells. Figure 20 shows the fracture distribution for each layer. Figure 21 shows the permeability distribution of seven layers for egg model. After defining connections and calculating transmissibility in preprocessing code, the simulation is simple to performance using in-house simulators. Two other typical algorithms, particle-swarm optimization (PSO) and ensemble-based optimization (EnOpt), are both used to study their performance on well placement optimization. Figure 22 shows the NPVs with respect to number of iterations for different methods. Figure 23 shows the oil saturation (first layer) at the final simulation time for different optimization method. Figure 24 shows the cumulative oil production and water cut for different optimization method. As can be seen, the highest NPV is for by using StoSAG. The final NPV for EnOpt and PSO is and , respectively. Also, when using StoSAG, the highest cumulative oil production can be achieved in 7200 days.

#### 4. Conclusions

In this work, we use StoSAG for the well placement optimization. The computer-aided well placement optimization method is established for naturally fractured reservoirs based on the recent advances. The embedded discrete fracture method (EDFM) is used to model the natural fractures as its computational efficiency and flexibility. The stochastic simplex approximate gradient (StoSAG) is employed to obtain the approximate gradient by combing the EDFM. The steepest ascent algorithm is used to find the optimal well placement. A series of numerical case studies are set up to examine the performance of the proposed approach. We also demonstrate that f-StoSAG and StoSAG and sf-StoSAG and ss-StoSAG can achieve fairly close results. The workflow can be taken as an effective tool in well placement optimization for naturally fractured reservoirs.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This study was supported by the Natural Science Foundation of Shandong Province of China (ZR2019BEE030).