Regular Article
Generalized MultiScale Stochastic Reservoir Opportunity Index for enhanced well placement optimization under uncertainty in green and brownfields
^{1}
Department of Petroleum Engineering, Amirkabir University of Technology, Hafez Ave. No. 424, Tehran, Iran
^{2}
Faculty of Economics and Business Administration, Ghent University, Tweekerkenstraat 2, 9000 Gent, Belgium
^{3}
Technology and Operations Management, Vlerick Business School, Reep 1, 9000 Gent, Belgium
^{4}
UCL School of Management, University College London, 1 Canada Square, London E14 5AA, United Kingdom
^{*} Corresponding author: m.ahmady@aut.ac.ir
Received:
19
September
2020
Accepted:
16
March
2021
Well placement planning is one of the challenging issues in any field development plan. Reservoir engineers always confront the problem that which point of the field should be drilled to achieve the highest recovery factor and/or maximum sweep efficiency. In this paper, we use Reservoir Opportunity Index (ROI) as a spatial measure of productivity potential for greenfields, which hybridizes the reservoir static properties, and for brownfields, ROI is replaced by Dynamic Measure (DM), which takes into account the current dynamic properties in addition to static properties. The purpose of using these criteria is to diminish the search region of optimization algorithms and as a consequence, reduce the computational time and cost of optimization, which are the main challenges in well placement optimization problems. However, considering the significant subsurface uncertainty, a probabilistic definition of ROI (SROI) or DM (SDM) is needed, since there exists an infinite number of possible distribution maps of static and/or dynamic properties. To build SROI or SDM maps, the kmeans clustering technique is used to extract a limited number of characteristic realizations that can reasonably span the uncertainties. In addition, to determine the optimum number of clustered realizations, HigherOrder Singular Value Decomposition (HOSVD) method is applied which can also compress the data for large models in a lowerdimensional space. Additionally, we introduce the multiscale spatial density of ROI or DM (D_{2}ROI and D_{2}DM), which can distinguish between regions of high SROI (or SDM) in arbitrary neighborhood windows from the local SROI (or SDM) maxima with low values in the vicinity. Generally, we develop and implement a new systematic approach for well placement optimization for both green and brownfields on a synthetic reservoir model. This approach relies on the utilization of multiscale maps of SROI and SDM to improve the initial guess for optimization algorithm. Narrowing down the search region for optimization algorithm can substantially speed up the convergence and hence the computational cost would be reduced by a factor of 4.
Coauthors emails: Forough.Vaseghi@UGent.be; m_sharifi@aut.ac.ir; Mario.Vanhoucke@UGent.be
© F. Vaseghi et al., published by IFP Energies nouvelles, 2021
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Nomenclature
D_{x}, D_{y}, D_{z}: Dimensions of the Reservoir
DM_{producer}: Dynamic Measure for producer
DM_{injector}: Dynamic Measure for injector
K_{ro} : Oil Relative Permeability, fraction
K_{rw} : Water Relative Permeability, fraction
ROI: Reservoir Opportunity index
S_{o} : Oil Saturation, fraction
ξ : Reservoir Opportunity Index Tensor
: Loworder Approximation of Tensor
D_{2} ROI: Multiscale spatial density of ROI
D_{2} DM: Multiscale spatial density of DM
SROI: Statistical Reservoir Opportunity Index
SDM: Statistical Dynamic Measure
1 Introduction
Identifying the optimum locations for production and/or injection wells in greenfields as well as infill wells in mature fields is a very challenging and important problem. Several factors including petrophysical property distribution, geological and structural heterogeneity, economic factors, operational constraints as well as the geological uncertainty influence the solution of this problem.
In some of studies conducted on well placement optimization, different types of derivativefree optimization algorithms, such as genetic algorithm, particle swarm optimization algorithm, and simultaneous perturbation stochastic approximation (Bangerth et al., 2006; Beckner and Song, 1995; Emerick et al., 2009; Onwunalu and Durlofsky, 2010) as well as gradientbased optimization algorithms, such as adjoint methods (Forouzanfar et al., 2010; Sarma et al., 2008; Zandvliet et al., 2008) have been used. Nonetheless, in the field applications, the use of these methods for a model having a large number of cells would be computationally very timeconsuming and costly. Therefore, the practical application of optimization algorithms is very limited.
Due to the significant computational cost needed for lots of simulation runs, some researchers have used the reduced order models or the statistical surrogates. For example, some have used proxy modeling, in which a proxy model substitutes the primary objective. Some of proxy models used in production optimization are least squares, kriging, radial basis functions and multiple regression techniques (Artus et al., 2006; Pan and Horne, 1998; Zubarev, 2009). Also, Bittencourt and Horne (1997), Guyaguler et al. (2002) have used hybrid algorithms by combining derivativefree optimization algorithms with surrogate models. In some studies, the combination of optimization algorithms and artificial neural network as a proxy for narrowing down the search region of optimization (Akin et al., 2010) has been used.
Additionally, recently there is a tendency towards considering the geological uncertainty in well placement optimization. The optimization under uncertainty has much higher computational cost due to need for running the simulations at each step on several realizations of reservoir model. Hence, some techniques have been developed to avoid or at least decrease the need for running many simulations. Among them, Reservoir Opportunity Index (or Simulation Opportunity Index) allows integrating the spatial maps of different static properties rather than visual inspection for finding suitable well locations. Another technique is the use of cumulative oil production maps, namely the quality maps, for example, see Liu and Jalali (2006), Da Cruz et al. (1999).
Among the studies considering the geological uncertainty, one can point to Lyons and Nasrabadi (2013) who have used an ensemble Kalman filter and optimization under timedependent uncertainty. Moreover, Taware et al. (2012) have used streamlinebased quality maps and calculated a dynamic measure based on the total time of flight to find the optimum well locations. In fact, they have created multiple geologic realizations to consider geological uncertainty and proposed a probabilistic definition of dynamic measure. Wang et al. (2012) have used a retrospective optimization framework under uncertainty.
The use of reservoir opportunity index for designing the optimum development plan improved the traditional well placement techniques with little computations needed (Molina and Rincon, 2009; Ghazali and Razib, 2011; VarelaPineda et al., 2014; Saputra and Ariaji, 2015). In these studies, the reservoir opportunity index is defined by combining the key elements influencing the well productivity such as the flow capacity, the mobile hydrocarbon pore volume, and the pressure.
Insuasty et al. (2017) have used HigherOrder Singular Value Decomposition (HOSVD) associated with kmeans clustering to handle the uncertainties within the production optimization problem. HOSVD was used to reduce the realizations dimensionality and hence it helped to compress the stored data. Thereafter, they performed the kmeans clustering based on a flowbased dissimilarity measure for each couple of reservoir realizations. Finally, the robust optimization under uncertainty was conducted on a few numbers of realizations clustered out of several realizations. To combat geological uncertainty and reduce the high number of scenarios, Meira et al. (2015) have used a mathematical function to model the representativeness of a subset of models with respect to the full set that characterizes the problem, and then an optimization tool to identify the representative models of any problem. Also, Shirangi and Durlofsky (2016) have introduced a general framework, based on clustering, for selecting a representative subset of realizations and quantified the performance of various realizationselection methods by computing the difference between the flow response for the subset and the full set.
In this paper, we have developed a method for well placement optimization based on the reservoir opportunity index (in undeveloped or greenfields) and the dynamic measure (in developed or brownfields) to achieve the maximum hydrocarbon recovery. In fact, our method offers faster convergence features due to the reduction of the algorithm search region. The regions with higher values of reservoir opportunity index or dynamic measure are considered as more favorable areas for drilling, i.e., sweet spots. Hence, the search regions of optimization algorithms are defined in the neighborhood of obtained sweets spots.
Theoretically, there is an infinite number of possible geological realizations, and the optimal point for drilling can be different in different realizations. From the perspective of probabilistic, the optimal point, in average, in all realizations, is close to the optimum location. To take the geological uncertainty into account, it is necessary to simulate a single scenario on several different realizations in each step of optimization, which will make this process much more timeconsuming. Therefore, to avoid running excessive simulation runs on all geological realizations, a clustering approach is needed to pick an affordable number of representative realizations for the optimization process. If the number of realizations used in optimization is low, that probabilistic interval may not adequately cover the uncertainty range. Nonetheless, the kmeans clustering technique allows selecting a limited number of realizations, which provide key information about the uncertainty and variability of properties in various realizations.
Moreover, we have defined the multiscale spatial density of ROI or DM (D_{2}ROI and D_{2}DM) for each of the representative realizations, because there is a risk of unsuccessful development scenario due to the geological uncertainties and likely well deviation from the target in drilling. We reason that those cells must be nominated for well placement for which ROI average over their neighborhood has higher values.
For large reservoir models, due to the large number of realizations, information should be saved within a limited volume. Hence, we use HigherOrder Singular Value Decomposition (HOSVD) to compress static and dynamic data. Moreover, HOSVD allows specifying the suitable number of clusters, such that most of the variability statistics are preserved after performing the kmeans clustering technique.
On the whole, we have generally addressed these main issues in our well placement optimization problem:
Defining Statistical Reservoir Opportunity Index (SROI) and/or Statistical Dynamic Measure (SDM) to account for the geological uncertainty and improving the convergence features of optimization algorithms by narrowing down the search region identified by SROI and/or SDM.
Introducing the multiscale spatial density of ROI or DM (D_{2}ROI and D_{2}DM) for each of the representative realizations to identify the locations with higher ROI average over their neighborhood.
Clustering a reduced representative set of realization to ensure spanning the uncertainty range within the optimization process.
Storing the data within more limited memory volume using HOSVD.
2 Methodology
2.1 Well placement optimization in greenfields
The following procedure has been developed for well placement optimization in a typical greenfield. 1000 realizations of porosity and permeability distribution within a synthetic reservoir model were propagated to resemble the full range of geological uncertainty.
Step 1. ROI map calculation for a large ensemble of realizations
1000 maps of Reservoir Opportunity Index (ROI) are generated based on 1000 realizations of porosity and permeability distribution. Equation (1) is used to calculate the ROI:
All these twodimensional ROI maps are gathered within a threedimensional tensor ξ in which the third dimension of this tensor shows the number of realizations.
Step 2. Reducedorder reconstruction of 3D ROI tensors using HOSVD
Thereafter, we apply HigherOrder Singular Value Decomposition (HOSVD) on 3D ROI tensor (HOSVD – Appendix A). It is required to reconstruct ROI tensor from a few numbers of the basis vectors of original tensor. These basis vectors are selected based on the threshold value used to cut off the tensor singular values. It is required to restrict the approximation error of ROI tensor at a reasonable level. By applying five different threshold values after HOSVD, we obtain five different sets of basis vectors for approximating the original ROI tensor.
HOSVD unfolds threedimensional ROI tensor into three matrices, and then decomposes each of these matrices through the Singular Value Decomposition (SVD) method. For data compression, original tensor is approximated by matrix reconstruction with a few basis vectors corresponding to the largest singular values. We can control the number of selected singular values and hence the number of basis vectors by the threshold value. By increasing the tolerance, the number of selected basis vectors for matrix reconstruction is reduced. This results in a tensor , which is a loworder approximation of the initial tensor. If the tensor approximation error is small, then the number of selected basis vectors suggests a clue of optimum number of clusters needed.
Step 3. Categorizing the realizations into an optimum number of clusters using kmeans clustering
After that, we use the kmeans clustering to select the appropriate number of clusters within 1000 realizations (Clustering – Appendix B). We categorize the realizations based on a distance measure in kmeans clustering such that the realizations with the similar characteristics are assigned to the same cluster.
In the tensor , consisting of 1000 matrices, an arbitrary matrix (here the first matrix) is considered as the reference matrix; and the Euclidean distance between reference matrix and other matrices is calculated using the definition of Frobenius norm. Therefore, the distance vector x (1 × 1000) provides a distance measure for clustering. The clustering output would be a vector L (1 × 1000), which includes the cluster index to which each realization has been assigned:
We conduct kmeans clustering for five cases (100 clusters, 60 clusters, 20 clusters, 10 clusters, and 5 clusters) based on 5 different threshold values applied in tensor reconstruction. Subsequently, the reservoir simulation under the same production scenario for 1000 realizations is performed. For each of the five cases defined in steps 2 and 3, a representative realization is determined for each cluster. At each cluster, a realization with the closest cumulative oil production (total field oil production, FOPT) value to the average of cumulative production values within that cluster is selected as the representative realization. Thereafter, for each of 5 cases, the statistical measures (including mean, standard deviation and the specific quantiles e.g. P10, P50, and P90) of cluster representatives FOPT values are compared with full reference case (1000 realizations). We calculate a Euclidean norm to quantify the difference of statistical measures between each clustered case and the reference case. The minimum number of clusters with difference norm less than a predefined threshold would be the optimum. Therefore, from now, we work with the representative realizations of clusters.
Step 4. Calculation of ROI multiscale density for representative realizations
Then, we calculate the multiscale density of ROI for each of the representative realizations. The common practice for the use of ROI maps or reservoir quality maps is to select well sweet spots from the cells with the highest ROI. However, there is a risk of unsuccessful development scenario due to the geological uncertainties and likely well deviation from the target in drilling. We reason that those cells must be nominated for well placement for which ROI average over their neighborhood has higher values. For each cell, the nthorder neighborhood is defined as a set containing the cells in the set of (n − 1)thorder neighborhood in addition to all other cells sharing an interface with cells in (n − 1)thorder neighborhood. Based on this, for each cell, we define nthorder ROI density as the arithmetic average of ROI values for cells in the corresponding nthorder neighborhood. We utilize the secondorder ROI density for well placement, which is defined as:
Step 5. Calculation of the statistical measures maps of second order ROI density (SROI)
Thereafter, the mean (expected value) and standard deviation maps of second order ROI density maps (SROI) over the representative realizations are calculated. Areas with the maximum probability of productivity (sweet spots) are identified by their high mean value and low standard deviation. The neighborhood of selected sweet spots is considered as the likely search regions for the optimization algorithm.
Step 6. Performing robust optimization with an evolutionary algorithm
With the search regions determined in the previous step, the Particle Swarm Optimization (PSO) algorithm is conducted to optimize the production well location under uncertainty. It means that well location is optimized over the ensemble of clustered representative realizations. In this robust optimization, the objective function is the average of cumulative oil production. The maximum number of iterations was set to 20 and the population size (swarm size) was 5. Results of PSO with localized regions around the sweet spots are compared against those obtained with a wide search region. Workflow showing overall steps in our proposed method is illustrated in Figure 1.
Fig. 1 Workflow showing overall steps. 
2.2 Well placement optimization in brownfields
The following procedure is developed for the optimization of infill well locations in a typical brownfield. This procedure has been implemented for the same synthetic reservoir model used in greenfield development optimization. It is assumed that the field has been produced for a while and now it is desired to find the optimum locations of production and injection wells.
Step 1. Selecting the appropriate number of clusters
This step is similar to step 3 from the greenfield procedure, with the difference that the simulation conditions have been changed, i.e., the time of production in simulation has been increased. Then, as in step 3 in the greenfield, the best number of clusters and the representative realizations of the clusters have been selected.
Step 2. Creating the separate Injector DM (IDM) and Producer DM (PDM) maps for clusters representatives
In order to find the best place for production and injection well in a brownfield, it is necessary to define a dynamic measure; because the distribution of pressure, saturation and other reservoir parameters have been altered after production. After simulation with initial development scenario in the greenfield, the outputs are required for the dynamic measure definition. These outputs for infill producers include permeability, oil relative permeability, oil saturation, pressure, and the bulk oil velocity, and for injectors, include the permeability, pressure, the bulk oil velocity, and the mobility ratio.
The suitable location of a production well can be different from the suitable location of an injection well. Hence, in DM definition proposed for the infill producers, the reservoir pressure, the oil saturation, the permeability and oil relative permeability are directly related. While there should be an inverse relationship between the producer DM with the bulk oil velocity, because low oil velocity and high oil saturation are indicators of undrained regions:
On the other hand, DM for the infill injectors is defined with a direct relationship with permeability an inverse relationship with the pressure, the bulk oil velocity, and the mobility ratio. This is because water sweep efficiency is in an inverse relationship with the mobility ratio. Moreover, drained and lowpressure areas with high permeability (resulting in high injectivity) are ideal locations for injection wells:
where
Step 3. Calculation of injector and producer DM multiscale density for representative realizations
After constructing the dynamic measure maps for each of the representative realizations, the maps of the secondorder injector/producer DM density are computed.
Step 4. Calculation of the statistical measures maps of second order DM density (SDM)
Afterward, two maps of statistical measures of second order DM density are created, which include the mean and Standard Deviation Maps (SDM) over the whole representative realizations. This allows recognizing the areas with the maximum potential of hydrocarbon production from the producer SDM maps. In addition, the areas with high potential of favorable sweep efficiency are identified by the injector SDM maps.
Step 5. Performing the optimal screening on candidate search regions
In order to increase the chance of finding the best well locations, optimal screening is performed to select the potential search regions for the optimization. Those cells are picked as the promising search regions for the optimization algorithm, which correspond to the injector/producer SDM mean value higher than a minimum threshold and the injector/producer SDM standard deviation value less than a maximum limit.
Step 6. Conducting the optimization algorithm to find the best well locations
Any optimization algorithm (here Particle Swarm Optimization (PSO) algorithm) is utilized to find the best locations for production and injection infill wells within the search areas provided in the previous step. The locations are optimized under uncertainty by taking the representative realizations into account. In this robust optimization, the objective function is the cumulative production of oil. The maximum number of iterations was set to 20 and the population size (swarm size) was 5.
3 Application
We have implemented our approach for a synthetic reservoir model.
3.1 Synthetic reservoir model description
A simple synthetic model with 40 × 40 × 1 cells is created with the cell dimensions of 40 ft × 40 ft × 6.5 ft. Distributions of permeability and porosity in the model are based on the Unconditional Sequential Gaussian Simulation (USGS). A view of the reservoir and the distribution of permeability and porosity are shown in Figure 2. It should be noted that in all the maps presented in this study, darkcolored areas represent lower amounts of the property (the lowest amount for violet areas), and bright areas indicate high levels of property (the highest amount for red areas). Also, the general characteristics of the reservoir are given in Table 1. For the reservoir rock, three types of rock (sand, shale and shaly sand) based on the reservoir quality index have been considered. In Table 2, the rock type definition based on the Reservoir Quality Index (RQI) is shown. The distribution of the rock type in the model for the representative clusters is shown in Figure C1 (Appendix C):
General characteristics of the reservoir.
The definition of rock type based on the reservoir quality index.
Fig. 2 Porosity (left) and permeability field (log10 – right) for synthetic 2D model. For all the maps in the paper, the lighter areas correspond to the higher values and the darker areas correspond to the lower values of the parameter. 
3.2 Procedural steps for greenfield development optimization
Here we follow the procedure explained in the methodology section.
Step 1: After constructing 1000 realizations of permeability, porosity, and then calculating ROI maps, a ROI tensor of size 40 × 40 × 1000 is formed which 1000 is the number of realizations and 40 is the number of grids in the x and y dimension. A number of realizations of permeability, porosity and ROI are shown in Figures C2–C4 (Appendix C).
Step 2: Then the HigherOrder Singular Value Decomposition (HOSVD) with five different tolerances (1, 2, 4, 8, 16) has been applied on this tensor. Higher tolerance means picking lower singular values resulting in more memory saving against less accurate approximation of original tensor. The number of retained singular values in the realization direction of tensor specifies the number of realizations clusters. The decomposed core tensor and the number of clusters corresponding to each of these five tolerances are shown in Table C1 (Appendix C). The error of the lowrank approximation defined by equation (8) is provided in Table C2 (Appendix C) for each of the five tolerances:
Step 3: In this step, the clustering operation for each of these five cases for reduced tensor has been performed. Afterward, by simulation on the 1000 realizations, a realization with the closest cumulative oil production value to the average value in each cluster, has been selected as the cluster representative. The simulation conditions are shown in Table 3. The probability density and the cumulative density functions of cumulative oil production for cluster representatives are plotted in Figure 3 for each of the five tolerances and the case of whole 1000 realizations (reference case). Also, the mean, standard deviation, P10, P50, and P90 values of cumulative oil for the reference case and five tolerance cases are presented in Table 4. All statistical measures including mean, standard deviation, P10, P50, and P90 values, as well as the probability density and cumulative density functions, are compared between each tolerance case and the reference case. It can be observed that in cases of 5 clusters, 10 clusters, and 20 clusters, the cost of lost data is more than other cases and these three cases are much far from the reference case. As a result, the case of 60 clusters, which has a good fit in terms of statistical measures with the reference case, has been selected as optimum with acceptable lowrank approximation error. Some representative realizations of permeability, porosity, and reservoir opportunity index are exhibited in Figures C5–C7 (Appendix C).
Fig. 3 PDF and CDF graphs for the five cases and reference case (greenfield). 
Simulation conditions to select the appropriate number of clusters (greenfield).
Comparison of mean, standard deviation. P10, P50, P90 values for 5 cases with the reference case (greenfield).
Step 4: In this step, the secondorder density of ROI has been calculated for representative realizations of the reservoir opportunity index. Some of these maps are shown in Figure C8 (Appendix C).
Step 5: In order to find the most promising areas, the mean and standard deviation maps of secondorder ROI density over 60 representatives have been created. The region corresponding to I ∈ [21, 28] J ∈ [33, 40] has been selected as the most probably productive area and the feed for the optimization algorithm. The mean and standard deviation SROI maps are shown in Figure C9 (Appendix C).
Step 6: Finally, optimization has been performed with 60 representative realizations of the clusters to select the best location for drilling production well in the greenfield. In Table C3 (Appendix C), the search regions are given for both cases of overall search and selective search. As we mentioned before, the objective function is defined as the average of cumulative oil production obtained for all representative realizations.
Due to the small size of the model and the lack of aquifer, the reservoir depletes very fast. To maintain the reservoir pressure at a higher level, such that production takes longer, a water injection well has been placed at the position (I, J) = (1, 1). The production time was set to 5 years. Also, the production and the injection wells have been set to constant pressure mode of control with bottomhole pressure values of 2000 psia and 6000 psia respectively.
3.3 Procedural steps for brownfield development optimization
Step 1: To select the appropriate number of clusters in the brownfield, simulation has been done on 1000 realizations with one production well at the optimum location identified from the previous section. The simulation conditions have been given in Table C4 (Appendix C). Similar to Step 3 of the greenfield optimization procedure, realizations are clustered in the five different cases (100 clusters, 60 clusters, 20 clusters, 10 clusters, and 5 clusters). In each cluster, a realization having closer cumulative oil production value to the average of that cluster was selected as the representative of that cluster. PDF and CDF plots of the cumulative oil production values for representatives of clusters are shown in Figure 4 for each of the five cases of clustering compared with 1000 realizations as the reference case. Also, the mean, the standard deviation, P10, P50, and P90 values are compared between the reference case and the five clustering cases as provided in Table 5. Based on this comparison, in cases of 5 clusters, 10 clusters, and 20 clusters, most of the statistical information is lost. Therefore, these cases cannot be suitable for a reasonable and representative classification of 1000 realizations. However, the case of 60 clusters, exhibits a good fit with reference case in terms of statistical measures and shape of PDF and CDF plots. Hence, the number of 60 clusters is identified as the minimum number of clusters which allows preservation of the original statistical information with an acceptable level of approximation. Some of the representative realizations of permeability are shown in Figure C10 (Appendix C).
Fig. 4 PDF and CDF graphs for the five cases and reference case (brownfield). 
Comparison of mean, standard deviation. P10, P50, P90 values for five cases with reference case (brownfield).
Step 2: The simulation time with optimum production well at greenfield condition is 5 years. The positions of production and injection wells are shown in Figure C11 (Appendix C). After the simulation, using the dynamic properties including pressure and saturation values at the start of the 6th year, the injection and production dynamic measures (PDM and IDM) have been calculated for each of the 60 representative realizations (the representative of the clusters). PDM and IDM maps for some of the realizations are demonstrated in Figures C12 and C13 (Appendix C).
Step 4: In terms of geological uncertainty, the mean value map and the standard deviation value map can define the confidence degree of the dynamic measures. These SDM maps have been made for 60 realizations (representative of clusters). Areas with higher mean value and lower standard deviation value are regarded as areas with higher productivity potential. These regions have been considered as the sweet spots. The mean and standard deviation maps of PDM and IDM and the sweet spots locations are shown in Figures C16 and C17 (Appendix C).
Step 5: In this step, the optimal screening technique has been used to select the most promising areas of the reservoir having a higher potential of productivity (or higher potential of sweep efficiency). Those cells with the mean value of more than 0.5 and the standard deviation value of less than 0.8, have been considered as the most likely areas of search for the brownfield optimization. The area of I ∈ [26, 40] and J ∈ [1, 40] has been selected as the optimization search region for the production well. Also, the region of I ∈ [24, 40] and J ∈ [27, 40], accordingly, has been selected as the optimization search region for injection well. These search regions are represented in Figures C18 and C19 (Appendix C).
Step 6: Similar to the greenfield optimization, we define and compare the results of two cases; in the first case, the whole model has been considered as the optimization search region, and in the second case, the search region has been limited to the more probable area for drilling wells. In Tables 6 and 7, the optimization algorithm search regions are given for both cases of overall search and selective search to find best location of producer and injector respectively. The simulation time for brownfield optimization is 5 years, which is restarted from a 5years simulation with the greenfield optimum strategy. The pressure at the production and injection wells has been set to 2000 psia and 6000 psia respectively. The aim is to find the optimum location for both production and injection infill wells in two different optimization problems. The maximum number of iterations in this optimization has been set to 20 and the population size (swarm size) is 5.
Optimization algorithm search region (brownfieldproduction well).
Optimization algorithm search region (brown fieldinjection well).
4 Results and validation
4.1 Results and validation for greenfield development optimization
As can be seen in Figure C20 (Appendix C), for the first case when the search region includes all grid cells, the objective function has reached its maximum value of 88916.8 at iteration number 15. The optimum well location found in this case is at (I, J) = (21, 34). One the other hand, in the second case that more promising area (sweet spots) has been considered as the PSO search region, in only two iterations, the objective function has reached the maximum amount of 89 235.6 that is higher than the value found in the first case after 15 iterations. This proves the impact of a good initial guess on both convergence rate and the quality of optimization results. The optimum well position obtained in this case is at (I, J) = (26, 38). Therefore, the grid cell (26, 38) was selected as the best location for production well in the greenfield case, as demonstrated in Figure C21 (Appendix C). The number of iterations and simulations performed for each case are given in Table C5 (Appendix C).
To validate the approach for the greenfield, the mean value map of reservoir opportunity index over 60 realizations has been compared to the mean value map of cumulative oil (FOPT) over 60 realizations after 5 years of production to create cumulative oil maps, 121 grid cells have been selected from 1600 grid cells and considered as a location for production well in the reservoir. The selection of these grid cells was such that the distance between two consecutive cells in the horizontal and vertical directions is 4. Afterward, production of reservoir is simulated for each case (121 cases) for each representative cluster (60 realizations). After simulation and calculation of the mean value map of cumulative oil maps over these 60 realizations, the value of cumulative oil for other grid cells is calculated through the interpolation method. As can be seen in Figure 5, the harmony between field cumulative oil production map and reservoir opportunity index map and also the highest value of cumulative oil in the area of I ∈ [21, 28] J ∈ [33, 40] confirm the selection of search region from reservoir opportunity index maps.
Fig. 5 Comparison of mean value map of reservoir opportunity index over 60 realizations with mean value map of cumulative oil over 60 realizations after 5 years of production. 
4.2 Results and validation for brownfield development optimization
As represented in Figure C22 (Appendix C), when the optimization algorithm searches the whole grid for the optimum producer location, the objective function has converged to its maximum value of 126 012 after 16 iterations. The optimum location of the production well was (I, J) = (6, 27). While with a selective search on higherpotential areas, the optimization algorithm has reached the maximum value of 129 530 in iteration number 3, when the optimized location was (I, J) = (39, 23). However, for the injector location optimization, the results have been shown in Figure C23 (Appendix C). With whole grid search, the objective function has reached 134 150 in iteration number 14. The optimized location of the injection well was (I, J) = (30, 19). However, with the selective search of sweetspot areas, objective function converges to the value of 136 770 after only four iterations. The optimized location of the injection well was (I, J) = (32, 28). The number of iterations and simulations performed for both producer and injector are given in Tables C6 and C7 (Appendix C) respectively.
It can be concluded that the optimization with a highpotential search region can converge in much fewer optimization to a higher value of objective function. This guarantees the superiority of developed methodology over the common optimization approaches due to saving in simulation time plus more reliable results of optimization.
To validate the method for the brownfield, the mean and standard deviation value maps of dynamic measure over 60 realizations have been compared to the mean value map of cumulative oil over 60 realizations after 10 years of production in this field. For this purpose, as we can see in Figure 6, the highest value of cumulative oil in the area of I ∈ [30, 40] and J ∈ [1, 40] confirm the region obtained from optimal screening the mean and standard deviation value maps of dynamic measure.
Fig. 6 Comparison of mean value map and standard deviation map of dynamic measure over 60 realizations with mean value map of cumulative oil over 60 realizations after 10 years (producer). 
4.2.1 Comparison between deterministic and optimistic approaches
In the deterministic approach, the best location is identified just based on the expected value and standard deviation maps of dynamic measure without optimization, but in the optimistic approach, after identifying the best search region and then optimization, the more potential location will find. In Table C8 (Appendix C) the results of the identified locations and corresponding cumulative oil for two approaches are given.
It can be seen that there is no significant difference between deterministic and optimistic approaches, and this is an indication that by using the expected value and standard deviation maps of dynamic measure, we can reach an acceptable solution for the well placement problem. These results make the methodology suitable for more complex and large fields, which are more time consuming to do optimization.
5 Conclusion
In this study, a new methodology has been developed for well placement optimization under uncertainty for both green and brownfields. The idea is to improve the initial guess of optimization algorithm. Another word, our approach limits the optimization search region to the higherpotential areas (sweetspots). We believe our methodology speeds up the optimization process, while results in higher (or at least comparable) values of objective function. In this regard, we defined the Statistical Reservoir Opportunity Index (SROI) for greenfields and also Statistical Dynamic Measure (SDM) for the mature fields. Moreover, a multiscale density function of SROI and SDM has been defined which is used to identify the sweet spots. Multiscale SROI/SDM density ensures the flexibility on target placement and reduces the risk created by geological uncertainty. With the search region specified after optimal screening on multiscale SROI/SDM density, the convergence time was roughly reduced to 0.25 of corresponding time for whole grid search. However, in both green and brownfield problems, the objective function has reached to higher levels compared to the whole grid search. Therefore, the main advantage of this method is reducing the cost and the time of optimization under geological uncertainty. This makes our approach to be suitable for largescale fields specially, in presence of geological uncertainty. To handle the curse of geological uncertainty, we utilize the kmeans clustering technique to cluster the realizations with the similar characteristics in one group. Thereafter, the optimization is conducted with a fewer number of realizations which allows to decrease the time and computational cost. To preserve the spatialgeological information during the clustering process, the HigherOrder Singular Value Decomposition (HOSVD) has been successfully used.
For future research efforts, we suggest consideration other parameters like cost and other types of risks and uncertainty when creating scenarios, defining the dynamic measure and also the objective function. Another suggestion for future studies could be investigating the impact of physical properties of fluids on the result of well placement optimization for more heterogenous and complex models, because phenomena like fingering, early breakthrough and high water cut are more likely to happen in these types of reservoirs and the effect of these extreme conditions on the well placement is inevitable and should be taken into account. For instance, parameters such as permeability variation coefficient can be added in the formula of dynamic measure.
References
 Akin S., Kok M.V., Uraz I. (2010) Optimization of well placement geothermal reservoirs using artificial intelligence, Comput. Geosci. 36, 776–785. [Google Scholar]
 Artus V., Durlofsky L.J., Onwunalu J., Aziz K. (2006) Optimization of nonconventional wells under uncertainty using statistical proxies, Comput. Geosci. 10, 389–404. [Google Scholar]
 Bangerth W., Klie H., Wheeler M., Stoffa P., Sen M. (2006) On optimization algorithms for the reservoir oil well placement problem, Comput. Geosci. 10, 303–319. [Google Scholar]
 Beckner B., Song X. (1995) Field development planning using simulated annealingoptimal economic well scheduling and placement, in: SPE Annual Technical Conference and Exhibition, Dallas, Texas, October 1995. [Google Scholar]
 Bittencourt A.C., Horne R.N. (1997) Reservoir development and design optimization, in: SPE Annual Technical Conference and Exhibition, San Antonio, Texas, October 1997. [Google Scholar]
 Da Cruz P.S., Horne R.N., Deutsch C.V. (1999) The quality map: a tool for reservoir uncertainty quantification and decision making, SPE Reserv. Eval. Eng. 7, 6–14. [Google Scholar]
 Emerick A.A., Silva E., Messer B., Almeida L.F., Szwarcman D., Pacheco M.A.C. (2009) Well placement optimization using a genetic algorithm with nonlinear constraints, in: SPE Reservoir Simulation Symposium, The Woodlands, Texas, February 2009. [Google Scholar]
 Forouzanfar F., Li G., Reynolds A.C. (2010) A twostage well placement optimization method based on adjoint gradient, in: SPE Annual Technical Conference and Exhibition, Florence, Italy, September 2010. [Google Scholar]
 Ghazali A.K., Razib A.R. (2011) Optimizing development strategy and maximizing field economic recovery through simulation opportunity index, in: SPE Reservoir Characterisation and Simulation Conference and Exhibition, Abu Dhabi, UAE, October 2011. [Google Scholar]
 Guyaguler B., Horne R.N., Rogers L., Rosenzweig J.J. (2002) Optimization of well placement in a Gulf of Mexico waterflooding project, SPE Reser. Eval. Eng. 5, 229–236. [Google Scholar]
 Insuasty E., Van den Hof P.M., Weiland S., Jansen J.D. (2017) Flowbased dissimilarity measures for reservoir models: a spatialtemporal tensor approach, Comput. Geosci. 21, 645–663. [Google Scholar]
 Liu N., Jalali Y. (2006) Closing the loop between reservoir modeling and well placement and positioning, in: Intelligent Energy Conference and Exhibition, Amsterdam, The Netherlands, April 2006. [Google Scholar]
 Lyons J., Nasrabadi H. (2013) Well placement optimization under timedependent uncertainty using an ensemble Kalman filter and a genetic algorithm, J. Pet. Sci. Eng. 109, 70–79. [Google Scholar]
 Meira L., Coelho G., Santos A.A.S., Schiozer D. (2015) Selection of representative models for decision analysis under uncertainty, Comput. Geosci. 88, 67–82. [Google Scholar]
 Molina A.R., Rincon A.A. (2009) Exploitation plan design based on opportunity index analysis in numerical simulation models, in: Latin American and Caribbean Petroleum Engineering Conference, Cartagena de Indias, Colombia, May 2009. [Google Scholar]
 Onwunalu J.E., Durlofsky L.J. (2010) Application of a particle swarm optimization algorithm for determining optimum well location and type, Comput. Geosci. 14, 183–198. [Google Scholar]
 Saputra W., Ariaji T. (2015) A new simulation opportunity index based software to optimize vertical well placements, public hearing of final project presentation, Department of Petroleum Engineering ITB. [Google Scholar]
 Sarma P., Durlofsky L.J., Aziz K. (2008) Computational techniques for closed–loop reservoir modeling with application to a realistic reservoir, Petrol. Sci. Technol. 26, 1120–1140. 10th European Conference on the Mathematics of Oil Recovery. [Google Scholar]
 Pan Y., Horne R.N. (1998) Improved methods for multivariate optimization of field development scheduling and well placement design, in: SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, September 1998. [Google Scholar]
 Shirangi M.G., Durlofsky L.J. (2016) A general method to select representative models for decision making and optimization under uncertainty, Comput. Geosci. 96, 109–123. [Google Scholar]
 Taware S.V., Park H.Y., DattaGupta A., Bhattacharya S., Tomar A., Kumar M. (2012) Well placement optimization in a mature carbonate waterflood using streamlinebased quality maps, in: SPE Oil and Gas India Conference and Exhibition, Mumbai, India, March 2012. [Google Scholar]
 VarelaPineda A., Hutheli A.H., Mutairi S.M. (2014) Development of mature fields using reservoir opportunity index: a case study from a Saudi field, in: SPE Saudi Arabia Section Technical Symposium and Exhibition, AlKhobar, Saudi Arabia, April 2014. [Google Scholar]
 Wang H., EcheverríaCiaurri D., Durlofsky L., Cominelli A. (2012) Optimal well placement under uncertainty using a retrospective optimization framework, SPE J. 17, 112–121. [Google Scholar]
 Zandvliet M., Handels M., Essen G., Brouwer R., Jansen J. (2008) Adjointbased wellplacement optimization under production constraints, SPE J. 13, 392–399. [Google Scholar]
 Zubarev D.I. (2009) Pros and cons of applying proxymodels as a substitute for full reservoir simulations, in: SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, October 2009. Journal of Petroleum Technology. 62, 4142. [Google Scholar]
Appendix A
A.1 Singular Value Decomposition (SVD)
The singular value decomposition of a matrix A (m × n) is the factorization of A into the product of three matrices:
where the columns of U (m × m) and V (n × n) are orthonormal and the matrix S (m × n) is diagonal with positive real entries and elements on its diagonal range from large values to low values. The columns of V in the singular value decomposition, called the right singular vectors of A, always form an orthogonal set with no assumptions on A. The columns of U are called the left singular vectors and they also form an orthogonal set. Elements on the original diagonal of matrix S are called eigenvalues of matrix A.
A can be decomposed into a sum of min{m, n} rank one matrices as
In many applications, the data matrix A is close to a matrix of low rank and it is useful to find a low rank matrix which is a good approximation to the data matrix.
A.2 Best rankk approximations
If we have a representation of matrix A as the sum of several elements that are arranged in order of importance, to have a best rankk approximation degree of matrix A, we should keep the first k elements as the most important of them:
A_{k} is the best rankk approximation to A. It is clear that A_{k} has rank k and is expressed as the sum of k rankone matrices.
The five steps to get the matrix A_{k} are summarized as follows:
Step 1: Calculate the singular value decomposition of matrix A (A = USV^{T}) where U is an orthogonal matrix m × m, S is a nonnegative diagonal matrix m × n and V^{T} is an orthogonal matrix n × n.
Step 2: Conserve only the first k singular vector of the right singular vectors. is equivalent to the first k row of the matrix V^{T}.
Step 3: Preserve only the first k singular vector of left singular vectors. U_{k} (m × k) is equivalent to the first k column of the matrix U.
Step 4: Preserving only first k eigenvalue of the matrix S. S_{k} (k × k) is equivalent to k the first column and k the first row of the matrix S.
Step 5: Then the rankk approximation is defined as follows:
Usually, if the number of large eigenvalues is small, k is equal to the number of these eigenvalues. Also, according to a rule of thumb, k must be chosen such that the sum of k upper eigenvalue at least c times is larger than the sum of the other eigenvalues (c is a domaindependent constant, for example 10).
A.3 HigherOrder Singular Value Decomposition (HOSVD)
Higherorder singular value decomposition or multilinear singular value decomposition is a generalization of singular value decomposition for higher order tensors (arrays of more than two dimensions) and is one of the tensor decomposition methods that with helping optimization the least squares of objective function by algorithms, decomposes the input matrix or tensor into the matrix elements.
Higherorder singular value decomposition of threedimensional tensor A involves the computation of the singular value decomposition of the three matrices, which are referred to as matrix modes. As a result, singular values with more importance and larger values are preserved. The number of singular values can be controlled by the tolerance. Finally, core tensor S and approximation tensor are built. The six steps of the HOSVD are summarized as follows:
Step 1: Build tensor A.
Step 2: The matrix unfoldings of tensor A (m × n × p) where we matricize the tensor in all three modes, creating three new matrices (1mode (U^{(1)}), 2mode, (U^{(2)}), 3mode (U^{(3)})):

Step 3: The application of SVD in all three new matrices.

Step 4: Keep the kmost important singular values for each matrix (krank of each mode).

Step 5: Build the core tensor S that reduces dimensionality:

Step 6: Build the approximation tensor:
Appendix B
B.1 Clustering
Clustering (unsupervised classification) means the classification of data, observations and vectors into different groups or clusters based on a similarity criterion and has been used in many topics and disciplines and it is very useful in data analysis. In fact, clustering analysis means organizing a set of patterns (usually as a vector of dimensions or a point in a multidimensional space) based on similar properties into clusters. That is, patterns in one cluster are more similar to patterns in other clusters.
In the clustering process, an appropriate similarity criterion must be selected after displaying the patterns. It is common to use distance criteria to calculate the dissimilarity between the two patterns. The distance helps us move through the data space and form clusters. By calculating the distance between two data, we can find out how close the two data are to each other and put them in a cluster. There are various mathematical functions to calculate the distance. The most popular method for continuous properties is the Euclidean distance method obtained from the following equation:
After selecting the similarity criterion, the clustering operation is performed.
B.2 kmeans clustering
kmeans clustering, which is considered a partitioning clustering technique, is a basic method for many other clustering methods. Various algorithms have been proposed for it but all of them have a repetitive procedure that is used for a fixed number of clusters. In a simple form of this method, points are selected at random equal to the number of clusters required. Then the data are assigned to one of these clusters according to their similarity and thus new clusters are obtained. By repeating the same procedure, new centers are calculated for each iteration by averaging the data and reassign the data to new clusters. This process continues until there is no change in the data. So the algorithm is defined as follows:
Step 1: Obtain points as the centers of the clusters where these points are actually the mean points belonging to each cluster.
Step 2: Assign each data to a cluster with the least distance to the center of that cluster.
Step 3: After belonging all the data to clusters, a new point is calculated for each cluster as the center (mean points belonging to each cluster).
Steps 2 and 3 are repeated until no changes are made to the cluster centers.
Appendix C
Fig. C1 Distribution of rock type based on reservoir quality index over realizations. 
Fig. C2 1000 permeability realizations used for generating ROI map under geologic uncertainty. 
Fig. C3 1000 porosity realizations used for generating ROI map under geologic uncertainty. 
Fig. C4 ROI maps for 1000 realizations [0.07 (darker areas) – 0.99 (lighter areas)]. 
The decomposed core tensor in different tolerances.
Comparison of approximation error in tensor decomposition in different tolerances.
Fig. C5 Permeability representatives for 60 clusters (greenfield). 
Fig. C6 Porosity representatives for 60 clusters (greenfield). 
Fig. C7 ROI representatives for 60 clusters [0.07 (darker areas) – 0.99 (lighter areas)]. 
Fig. C8 Secondorder density of ROI for representative of 60 clusters. 
Fig. C9 Mean value map and standard deviation value map of secondorder density over 60 realizations (greenfield). 
Optimization algorithm search region (green field).
Simulation conditions to select the appropriate number of clusters (brownfield).
Fig. C10 Permeability representatives for 60 clusters (brownfield). 
Fig. C11 Position of injection well and production well obtained from optimization in the greenfield. 
Fig. C12 Dynamic measure representatives (production well) for 60 clusters [0.19 (darker areas) – 9.87 (lighter areas)]. 
Fig. C13 Dynamic measure representatives (injection well) for 60 clusters [0.09 (darker areas) – 9.88 (lighter areas)]. 
Step 3: After calculating dynamic measures for 60 realizations (representatives of clusters), the secondorder density maps of the IDM and PDM have been computed. Some of these IDM and PDM density maps are shown in Figures C14 and C15 (Appendix C).
Fig. C14 Secondorder density of dynamic measure representatives (production well) for 60 clusters. 
Fig. C15 Secondorder density of dynamic measure representatives (injection well) for 60 clusters. 
Fig. C16 Mean and standard deviation value map of secondorder density of dynamic measure (production well) over 60 realizations. 
Fig. C17 Mean and standard deviation value map of secondorder density of dynamic measure (injection well) over 60 realizations. 
Fig. C18 Sweet spots for production well using optimal screening technique. 
Fig. C19 Sweet spots for injection well using optimal screening technique. 
Fig. C20 Optimization results for case 1 (whole grid search) and case 2 (selective sweetspot search) – greenfield. 
Fig. C21 The position obtained for production well – greenfield. 
Number of iterations and simulations performed for each case – greenfield.
Fig. C22 Optimization results for case 1 and case 2 (brownfield – production well). 
Fig. C23 Optimization results for case 1 and case 2 (brownfield – injection well). 
Number of iterations and simulations performed for each case – brownfield (producer).
Number of iterations and simulations performed for each case – brownfield (injector).
Comparison between deterministic and optimistic approaches.
All Tables
Comparison of mean, standard deviation. P10, P50, P90 values for 5 cases with the reference case (greenfield).
Comparison of mean, standard deviation. P10, P50, P90 values for five cases with reference case (brownfield).
Comparison of approximation error in tensor decomposition in different tolerances.
Simulation conditions to select the appropriate number of clusters (brownfield).
Number of iterations and simulations performed for each case – brownfield (producer).
Number of iterations and simulations performed for each case – brownfield (injector).
All Figures
Fig. 1 Workflow showing overall steps. 

In the text 
Fig. 2 Porosity (left) and permeability field (log10 – right) for synthetic 2D model. For all the maps in the paper, the lighter areas correspond to the higher values and the darker areas correspond to the lower values of the parameter. 

In the text 
Fig. 3 PDF and CDF graphs for the five cases and reference case (greenfield). 

In the text 
Fig. 4 PDF and CDF graphs for the five cases and reference case (brownfield). 

In the text 
Fig. 5 Comparison of mean value map of reservoir opportunity index over 60 realizations with mean value map of cumulative oil over 60 realizations after 5 years of production. 

In the text 
Fig. 6 Comparison of mean value map and standard deviation map of dynamic measure over 60 realizations with mean value map of cumulative oil over 60 realizations after 10 years (producer). 

In the text 
Fig. C1 Distribution of rock type based on reservoir quality index over realizations. 

In the text 
Fig. C2 1000 permeability realizations used for generating ROI map under geologic uncertainty. 

In the text 
Fig. C3 1000 porosity realizations used for generating ROI map under geologic uncertainty. 

In the text 
Fig. C4 ROI maps for 1000 realizations [0.07 (darker areas) – 0.99 (lighter areas)]. 

In the text 
Fig. C5 Permeability representatives for 60 clusters (greenfield). 

In the text 
Fig. C6 Porosity representatives for 60 clusters (greenfield). 

In the text 
Fig. C7 ROI representatives for 60 clusters [0.07 (darker areas) – 0.99 (lighter areas)]. 

In the text 
Fig. C8 Secondorder density of ROI for representative of 60 clusters. 

In the text 
Fig. C9 Mean value map and standard deviation value map of secondorder density over 60 realizations (greenfield). 

In the text 
Fig. C10 Permeability representatives for 60 clusters (brownfield). 

In the text 
Fig. C11 Position of injection well and production well obtained from optimization in the greenfield. 

In the text 
Fig. C12 Dynamic measure representatives (production well) for 60 clusters [0.19 (darker areas) – 9.87 (lighter areas)]. 

In the text 
Fig. C13 Dynamic measure representatives (injection well) for 60 clusters [0.09 (darker areas) – 9.88 (lighter areas)]. 

In the text 
Fig. C14 Secondorder density of dynamic measure representatives (production well) for 60 clusters. 

In the text 
Fig. C15 Secondorder density of dynamic measure representatives (injection well) for 60 clusters. 

In the text 
Fig. C16 Mean and standard deviation value map of secondorder density of dynamic measure (production well) over 60 realizations. 

In the text 
Fig. C17 Mean and standard deviation value map of secondorder density of dynamic measure (injection well) over 60 realizations. 

In the text 
Fig. C18 Sweet spots for production well using optimal screening technique. 

In the text 
Fig. C19 Sweet spots for injection well using optimal screening technique. 

In the text 
Fig. C20 Optimization results for case 1 (whole grid search) and case 2 (selective sweetspot search) – greenfield. 

In the text 
Fig. C21 The position obtained for production well – greenfield. 

In the text 
Fig. C22 Optimization results for case 1 and case 2 (brownfield – production well). 

In the text 
Fig. C23 Optimization results for case 1 and case 2 (brownfield – injection well). 

In the text 