© F. Vaseghi et al., published by IFP Energies nouvelles, 2021

Nomenclature

D_x, D_y, D_z: Dimensions of the Reservoir

D₂ : Second-order Density

DM_producer: Dynamic Measure for producer

DM_injector: Dynamic Measure for injector

e: Calculation error

K: Permeability, md

K_ro : Oil Relative Permeability, fraction

K_rw : Water Relative Permeability, fraction

L : Clustering Output

MR: Mobility Ratio

P : Reservoir Pressure, psia

ROI: Reservoir Opportunity index

RQI: Reservoir Quality Index

S_o : Oil Saturation, fraction

V_o : Oil Velocity

V_w : Water Velocity

X : Distance criteria

ϕ : Porosity, fraction

ξ : Reservoir Opportunity Index Tensor

: Low-order Approximation of Tensor

μ_o : Oil Viscosity, cp

μ_w : Water Viscosity, cp

D₂ ROI: Multiscale spatial density of ROI

D₂ 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 derivative-free 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 gradient-based 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 time-consuming 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 derivative-free 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 time-dependent uncertainty. Moreover, Taware et al. (2012) have used streamline-based 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; Varela-Pineda 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 Higher-Order Singular Value Decomposition (HOSVD) associated with k-means 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 k-means clustering based on a flow-based 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 realization-selection 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 time-consuming. 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 k-means 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₂ROI and D₂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 Higher-Order 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 k-means 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₂ROI and D₂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:

(1)

All these two-dimensional ROI maps are gathered within a three-dimensional tensor ξ in which the third dimension of this tensor shows the number of realizations.

Step 2. Reduced-order reconstruction of 3D ROI tensors using HOSVD

Thereafter, we apply Higher-Order 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 three-dimensional 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 low-order 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 k-means clustering

After that, we use the k-means clustering to select the appropriate number of clusters within 1000 realizations (Clustering – Appendix B). We categorize the realizations based on a distance measure in k-means 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:

(2)

We conduct k-means 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 multi-scale 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 nth-order neighborhood is defined as a set containing the cells in the set of (n − 1)th-order neighborhood in addition to all other cells sharing an interface with cells in (n − 1)th-order neighborhood. Based on this, for each cell, we define nth-order ROI density as the arithmetic average of ROI values for cells in the corresponding nth-order neighborhood. We utilize the second-order ROI density for well placement, which is defined as:

(3)

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.

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:

(4)

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 low-pressure areas with high permeability (resulting in high injectivity) are ideal locations for injection wells:

(5)

where

(6)

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 second-order 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, dark-colored 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):

(7)

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 Higher-Order 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 low-rank approximation defined by equation (8) is provided in Table C2 (Appendix C) for each of the five tolerances:

(8)

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 low-rank approximation error. Some representative realizations of permeability, porosity, and reservoir opportunity index are exhibited in Figures C5–C7 (Appendix C).

Step 4: In this step, the second-order 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 second-order 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 bottom-hole 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).

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 5-years 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.

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.

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 higher-potential 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 sweet-spot 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 high-potential 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.

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 higher-potential areas (sweet-spots). 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 multi-scale density function of SROI and SDM has been defined which is used to identify the sweet spots. Multi-scale 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 multi-scale 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 large-scale fields specially, in presence of geological uncertainty. To handle the curse of geological uncertainty, we utilize the k-means 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 spatial-geological information during the clustering process, the Higher-Order 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

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:

(A1)

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

(A2)

(A3)

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 rank-k 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 rank-k approximation degree of matrix A, we should keep the first k elements as the most important of them:

(A4)

A_k is the best rank-k approximation to A. It is clear that A_k has rank k and is expressed as the sum of k rank-one 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 non-negative 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 rank-k approximation is defined as follows:

(A5)

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 domain-dependent constant, for example 10).

A.3 Higher-Order Singular Value Decomposition (HOSVD)

Higher-order singular value decomposition or multi-linear 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.

Higher-order singular value decomposition of three-dimensional 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 (1-mode (U⁽¹⁾), 2-mode, (U⁽²⁾), 3-mode (U⁽³⁾)):

(A6)

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

• Step 4: Keep the k-most important singular values for each matrix (k-rank of each mode).

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

(A7)

• Step 6: Build the approximation tensor:

(A8)

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:

(B1)

After selecting the similarity criterion, the clustering operation is performed.

B.2 k-means clustering

k-means 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

Step 3: After calculating dynamic measures for 60 realizations (representatives of clusters), the second-order 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).

All Tables

All Figures