Well Placement Optimization for Fractured Reservoirs: Coupling StoSAG and EDFM

Well placement optimization is a signi ﬁ cant task in oil ﬁ eld development which aims to ﬁ nd 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 signi ﬁ cant 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 ﬂ ow discrete fracture numerical simulation model, i.e., the embedded discrete fracture method (EDFM) is used to model the natural fractures as its computational e ﬃ ciency and ﬂ exibility 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 ﬁ nd 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 ﬂ ooding naturally fractured reservoir production optimization substantially increased by using StoSAG.


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][3][4][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][8][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 fracturematrix (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][14][15][16][17][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][20][21]. Then, it is extended to the flow performance analysis of both naturally fracture reservoirs and hydraulic fractured tight oil reservoirs [22][23][24]. Some applications of EDFM have been reported. Yao et al. coupled the EDFM and dualcontinuum 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][34][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 Sto-SAG 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].

Matrix-matrix connection
Fracture-matrix connection Fracture-fracture connection: type I Fracture-fracture connection: type II

Geofluids
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 twophase 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    3 Geofluids 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.

Numerical Simulation Model.
In order to perform the well placement optimization for naturally fractured reser-voir, 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     5 Geofluids 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.

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 N e is the geological model realization number to describe the reservoir uncertainty. ucontains the placement    6 Geofluids information for all wells. ∇ u Jðu, m i Þ is the stochastic approximation of the simplex gradient. ∇ u Jðu, m i Þ is obtained by Then, the gradient becomes The superscript sign "+" on a matrix denotes the Moore-   4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Ln (k) Figure 12: Log-permeability distribution example 3.

Geofluids
N s represents the number of control perturbations. Each control perturbationû l,j ,j = 1, 2, ⋯N s at iteration l is generated from the distribution Nðu l , C U Þ; C U is a predefined covariance matrix.
where u is a N u -dimensional column vector which contains all well placement information; ndenotes the n th time step for the reservoir simulation; N t is the total number of time steps; the time at the end of the n th time step is denoted by t n ; t n is the n th time step size; b is the annual discount rate;    10 Geofluids N P and N I denote the number of producers and injectors, respectively; r o is the oil revenue, in $/STB; c w and c wi denote the disposal cost of produced water and the cost of water injection in units of $/STB, respectively; q o,j and q w,j , respectively, denote the average oil production rate and the average water production rate for the n th time step, in units of STB/ day; and q wi,k denote the average water injection rate at the k th injector for the n th 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 N e geological realizations: We consider only bound constraints and let u low and u up 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 k = 0, 1, 2 ⋯ until convergence, where x 0 is the initial guess and x k is the estimate of the optimal control parameter at the k th iteration; α k is the step size.

Workflow.
The workflow of well placement optimization using StoSAG and EDFM is shown in Figure 2. Firstly, N e geological realizations should be generated. In robust optimization, N e is commonly set to be bigger than 1. The initial well placement u 0 is used to calculate the initial objective function value J E ðu 0 Þ. 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.

Case Study
In  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 $2:66 × 10 7 for case 5, and the lowest NPV is $2:6 × 10 7 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 $2:85 × 10 7 for case 1, and the lowest NPV is $2:55 × 10 7 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.  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 $ 4:03 × 10 7 for case 6 and the lowest NPV is $3:89 × 10 7 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 Permeability (mD)  13 Geofluids be seen that the highest NPV is $3:54 × 10 7 for case 4 and the lowest NPV is $3:44 × 10 7 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 14 Geofluids 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 $1:02 × 10 8 for by using StoSAG. The final NPV for EnOpt and PSO is $9:56 × 10 7 and $9:65 × 10 7 , respectively. Also, when using StoSAG, the highest cumulative oil production can be achieved in 7200 days.

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.