Closed-loop field development optimization with multipoint geostatistics and statistical assessment

12/01/2017 ∙ by Mehrdad G Shirangi, et al. ∙ Stanford University 0

Closed-loop field development (CLFD) optimization is a comprehensive framework for optimal development of subsurface resources. CLFD involves three major steps: 1) optimization of full development plan based on current set of models, 2) drilling new wells and collecting new spatial and temporal (production) data, 3) model calibration based on all data. This process is repeated until the optimal number of wells is drilled. This work introduces an efficient CLFD implementation for complex systems described by multipoint geostatistics (MPS). Model calibration is accomplished in two steps: conditioning to spatial data by a geostatistical simulation method, and conditioning to production data by optimization-based PCA. A statistical procedure is presented to assess the performance of CLFD. Methodology is applied to an oil reservoir example for 25 different true-model cases. Application of a single-step of CLFD, improved the true NPV in 64 cases. The full CLFD procedure (with three steps) improved the true NPV in 96 of cases, with an average improvement of 37



There are no comments yet.


page 9

page 11

page 12

page 13

page 16

page 17

page 18

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.


Closed-loop field development (CLFD) optimization is a comprehensive framework for optimal development of subsurface resources. CLFD involves three major steps: 1) optimization of full development plan based on current set of models, 2) drilling new wells and collecting new spatial and temporal (production) data, 3) model calibration based on all data. This process is repeated until the optimal number of wells is drilled. This work introduces an efficient CLFD implementation for complex systems described by multipoint geostatistics (MPS). Model calibration is accomplished in two steps: conditioning to spatial data by a geostatistical simulation method, and conditioning to production data by optimization-based PCA. A statistical procedure is presented to assess the performance of CLFD. Methodology is applied to an oil reservoir example for 25 different true-model cases. Application of a single-step of CLFD, improved the true NPV in 64%–80% of cases. The full CLFD procedure (with three steps) improved the true NPV in 96% of cases, with an average improvement of 37%.

1 Introduction

For optimal operation of oil and gas resources, engineers typically build physics-based flow simulation models. These models need to be calibrated with observed spatial and temporal data so that they can be reliable for purposes such as optimization and decision making. The common approach for management of these systems is closed-loop optimization which involves three major steps, repeated throughout the project life: 1) optimization of decision parameters, 2) implementation and operation for a time period, 3) and model calibration for consistency with (new) observed data. For different types of reservoir problems, various approaches have been proposed for Steps 1 and 3, i.e., the optimization and the model calibration step. In addition, implementations of closed-loop optimization may differ on type of decision parameters included in these two steps.

Closed-loop field development (CLFD) optimization, introduced by Shirangi and Durlofsky (2015), is a comprehensive reservoir management framework. In the optimization step of CLFD, the number, type, location and controls of new wells (together with controls of the existing wells) are optimized. The model calibration step of CLFD involves integration/assimilation of new spatial data (hard data from new wells) and temporal data (production data from existing wells). The impact of CLFD research is significant as drilling new wells is one of the most expensive parts of reservoir operations. In Shirangi and Durlofsky (2015), the CLFD framework is presented with two-point geostatistical models. In recent years, however, there has been significant research on multipoint geostatistics (MPS) which is able to model complex geological features such as channels and deltaic fans.

In this work, a closed-loop field development methodology is developed for reservoir models described by multipoint geostatistics. A two-step model calibration procedure is presented where we employ an efficient geostatistical algorithm to integrate spatial data, and a recent model calibration algorithm for integration of temporal data (production history). A statistical procedure is then presented and applied to assess the performance of this closed-loop optimization method. In this assessment, the project outcome (true NPV) is treated as a random variable and the improvement in true NPV is investigated after each CLFD step.

Before the introduction of CLFD framework, most of the reservoir optimization research has focused on developing improved approaches for optimal operation of existing wells, referred to as closed-loop reservoir management (CLRM). Compared to CLFD, CLRM does not involve drilling new wells. In optimization step of CLRM, referred to as production optimization or well control optimization, continuous operational settings of existing wells are optimized. The model calibration of CLRM only includes new temporal data, and integration of new spatial data is not typically considered. CLRM has been investigated extensively (see, e.g., Jansen et al (2009)). Most papers on CLRM investigated the application of particular history matching and optimization approaches for water-flooding operations (Aitokhuehi and Durlofsky, 2005; Sarma et al, 2006; Chen et al, 2009). CLRM has also been applied to SAGD operations (Regtien, 2010), and for the management of geological carbon storage operations (Cameron and Durlofsky, 2014).

In the field development planning (FDP) problem, encountered in CLFD, decision parameters may include the number of new wells, well type (producer or injector), well locations, drilling sequence, and well settings. Most papers on FDP considered the optimization of only a subset of these parameters, mainly the location of new wells (referred to as the well placement problem). Many of the methods developed for this problem entail stochastic (global) search methods such as GAs (Ozdogan and Horn, 2006), PSO (Onwunalu and Durlofsky, 2010; Arnold et al., 2016), evolution strategy with covariance matrix adaptation (CMA-ES) (Bayer and Finkel, 2007; Bayer et al, 2008; Bouzarkouna et al, 2012), and differential evolution (DE) (Bayer et al, 2010; Nwankwor et al, 2013). Local optimization methods such as SPSA (Jahandideh and Jafarpour, 2016; Shirangi and Mukerji, 2012), and pattern search techniques (Wilson and Durlofsky, 2013; Cameron and Durlofsky, 2012) have also been applied for well placement optimization.

Recent research investigated the more complex problem of jointly optimizing various parameters in FDP. Bellout et al (2012) introduced a nested approach for the joint optimization of well location and controls in which the outer well placement optimization is solved using a pattern search optimization method, while the inner well control problem is solved using gradient-based SQP. Li and Jafarpour (2012) presented an alternating iterative solution of the decoupled well placement and control subproblems, where each subproblem is solved in turn. Zhang et al (2017) presented an ensemble-based method for solving the joint optimization of well location and controls. Isebor et al (2014a, b) developed a formulation based on a hybrid of PSO, a global stochastic search algorithm, and mesh adaptive direct search (MADS), a local pattern search method. This PSO-MADS procedure can simultaneously optimize the number and type (e.g., injector or producer) of new wells and the drilling sequence, in addition to well locations and controls. This algorithm will be used in this work for the CLFD optimization step.

Model calibration is the third major step in CLFD. In model calibration (also referred to as history matching, or data assimilation), the goal is to learn from data through generating one or more computational models that are consistent with prior geological information and new spatial data, and provide flow simulation results that closely match observed (temporal) production data. The model calibration methods typically involve a challenging minimization problem as this problem is usually ill-posed and the number of unknown model parameters can be very large. The ill-posedness of model calibration can be mitigated by reducing the number of parameters through an appropriate parameterization such as TSVD (Shirangi, 2011, 2014; Shirangi and Emerick, 2016; Bjarkason et al, 2017; Dickstein et al, 2017), and PCA (Sarma et al, 2008; Vo and Durlofsky, 2016). Vo and Durlofsky (2015, 2014) presented a differentiable PCA-based parameterization (O-PCA) that enables application of efficient gradient-based approaches for model calibration of complex channelized systems. This optimization-based principal component analysis (O-PCA) approach is used in the CLFD model calibration step in this work.

The CLFD framework has also been applied in recent work by other researchers. Morosov and Schiozer (2016) applied the CLFD framework to a realistic example, while Hanea et al (2017) applied the CLFD framework in simpler settings where only drilling sequence is optimized. Compared with the original work (Shirangi and Durlofsky, 2015; Shirangi, 2013), they used different reservoir modeling and different model calibration and optimization approaches in their CLFD application. Both papers, reported that CLFD significantly improved true-model NPV in most cases; however, in at least one case in each work, the true-model NPV decreased slightly after the application of CLFD. Most recently, Hidalgo et al (2017) applied the CLFD algorithm to a realistic reservoir case and reported 40.8% improvement in (true) NPV.

In Shirangi and Durlofsky (2015), reservoir properties were represented in terms of two-point spatial statistics. In this work, the CLFD framework is extended to complex channelized models, described by a training image. A training image contains a range of valuable information such as orientation, proportions, affinity, trend that need to be reproduced (Tahmasebi et al, 2014; Mariethoz and Caers, 2014; Mustapha and Dimitrakopoulos, 2011). This extension involves particular treatments for integration of spatial and temporal data in the model calibration step. This work also presents the results of a significant computational experiment, involving millions of reservoir simulation runs, to statistically assess the performance of this new CLFD approach and other treatments in terms of improvement in true NPV.

This paper proceeds as follow. In the next section, the CLFD framework and treatments for the model calibration and the optimization steps are described. In the section after that, computational results are presented for a two-dimensional binary channelized reservoir case. Results of statistical assessment of CLFD performance is also presented in this section. Finally, conclusions and suggestions for future work are provided.

2 Methodology: closed-loop field development optimization

Closed-loop field development (CLFD) (Fig. 1) has three major steps as follows. 1) Optimization: optimize the “full development plan” (the number, type, location and controls of all new wells together with controls of existing wells), for a set of models describing the current state of knowledge, 2) operation: drill and complete new wells, produce from existing wells, and collect new spatial and temporal data, 3) model calibration (history matching): update models for consistency with new data. These steps are repeated until drilling additional wells does not increase the expected NPV, at which step, the CLRM workflow is applied. We now describe the CLFD workflow in more detail.

Figure 1: Schematic and notations for CLFD framework with =2.

We consider the oil field development optimization problem where parameters of new wells are determined. In practical settings, the number of drilling rigs is limited and new wells are drilled sequentially, one or a few at a time. We let denote the number of rigs (in this work, we set =2), and we let denote the duration (in days) that it takes to drill and complete a new well and start production/injection. The initial field development plan is determined by optimizing the (full) development plan over a prior set of geological models, generated by use of a geostatistical simulation method. We let , , denote the time that the th group of wells is drilled. The first group of wells is drilled at .

After drilling and completion of each group of new wells, new spatial data (hard data from well locations) and new temporal data (production data) are collected. The CLFD model calibration is performed at each ( to update models for consistency with all data. Here, synthetic observed data is generated by simulating the true-model and adding random noise to the true data. The true-model here, is a realization randomly selected from the prior set.

The optimization is then performed to determine decision parameters for all future wells (those to be delivered at and later steps), and controls of existing wells. Optimizing the full development plan at each step is to avoid a greedy approach where one would only optimize the next wells. CLFD workflow is shown in Fig. 1. The optimization and model calibration steps (together with notations) are discussed next, in turn.

2.1 CLFD optimization

In CLFD, we use , a set of models, to represent the current geological knowledge at step . At each CLFD optimization step, the expected net present value (NPV) of the project is optimized for a set of representative models selected from the full set. The problem is formulated as



is the vector of decision parameters,

is a set of representative models ( is a subset of ), and denote the lower and upper bound vectors, is the number of inequality constraints, and is the expected NPV defined as


Computing NPV value for each model involves a flow simulation run. In this work, the water-flooding process with two-phase flow is considered where NPV is computed as


where is the cost of drilling and completing a well, is the time of delivering well , is the vector of reservoir model which contains properties such as permeability values at grid blocks, is the number of simulation time steps, is the simulation time (in days), and is the annual discount rate. Variables and denote the number of producers and injectors, respectively, and , and indicate the oil price and the cost of handling produced and injected water (all in ). Variables and denote oil and water production rates for producer at simulation time step , is the water injection rate of injector (all in STBD), and is the size (in days) of simulation time step .

The vector of decision parameters, , specify the number of wells, their locations and controls, well type (producer/injector) and the sequence in which new wells are drilled. The optimization problem in (1) is solved through a parallelized PSO-MADS optimization algorithm (Isebor et al, 2014b). Nonlinear constraint here is the minimum well distance constraint which is handled through the filter method in PSO-MADS. The PSO-MADS algorithm converges when a maximum number of iterations or a minimum stencil length is reached. This algorithm is applied within the optimization with sample validation (OSV) procedure, discussed next.

In CLFD optimization at any , the NPV for each model in Eq. 3 is computed from time zero to the end of specified project life. This is to ensure that we can compare the value of the project at different steps. We let denote the optimal solution obtained from the optimization at step . Note that other economic measures such as rate of return could be incorporated in the objective function here (Shirangi et al, 2017).

2.2 Optimization with sample validation and selection of representative models

We now discuss the optimization with sample validation (OSV) procedure from Shirangi and Durlofsky (2015). Ideally, the optimization should be performed over the full set of models (here, we set = 50). Since the computational cost scales linearly with the number of realizations used, an optimization over the full set will be computationally expensive. Therefore, we perform the optimization over a set of representative models (e.g., ). We will discuss later how these models are selected.

The optimization over the representative subset, improves the initial objective value of to an optimal value of . We then apply this optimal solution to the full set of models. The NPV improvement for the full set of models is therefore - . The relative improvement, , is computed as the ratio of NPV improvement for the full set, divided by that computed for the representative subset, i.e.,


We require to accept as the optimal solution at step () and we set = 0.5. If this criterion is not satisfied, we increase and repeat the optimization, until 0.5 is satisfied, or the maximum number of iterations is reached. In this work, a sequence of models is considered in OSV, i.e., the optimization is first performed over models. If the criterion is not satisfied, the optimization is repeated for a (new) set of models selected based on the current solution. The OSV procedure is outlined as follows.

  1. At CLFD step , initialize the OSV iteration index to , and set where is the optimal solution from the previous CLFD optimization step. Also set the initial value of .

  2. Select representative models from the full set of realizations by use of flow simulation results based on .

  3. Solve the optimization problem in (1) with the optimal solution denoted by .

  4. Compute in Eq. 4 by replacing with and with .

  5. If , accept the optimal solution and set .
    Otherwise, set , increase , and go to 2.

We now discuss the selection of representative models in CLFD. There are various methods presented for realization selection. Here we use a CDF approach as explained in Shirangi (2017). In this approach, flow simulation is performed for the full set of models by use of current

, and the NPV values are computed. The cumulative distribution function (CDF) plot is then generated. The models are selected such that two of them correspond to P10 and P90 (for

) or P5 and P95 (for ), and the rest of models correspond to even increments in NPV percentile (e.g., for , the models correspond to P10, P30, P50, P70, and P90). Other selection methods such as equally weighted flow-based and permeability-based clustering (Shirangi and Durlofsky, 2016) could be applied and tested here.

2.3 CLFD model calibration for channelized models

In CLFD model calibration step, the realizations are updated with new hard data and production data. Conditioning to production data is an inverse problem where one or multiple models need to be generated such that they reproduce the observed flow data and honor prior geological information. In the context of inverse problem theory (Tarantola, 2005)

, generating multiple calibrated (history-matched) models correspond to a sampling of the posterior probability density function (pdf). A common approach for solving this problem is the randomized maximum likelihood (RML)

(Kitanidis, 1995; Oliver et al, 1996; Shirangi, 2014). With RML, a sample from the posterior pdf is generated by minimizing an objective function that quantifies the mismatch between observed and predicted data. This objective function also has a model mismatch term to preserve the prior geological information.

For model calibration of a channelized reservoir described by MPS, the optimization-based principal component analysis (O-PCA) parameterization (Vo and Durlofsky, 2014, 2015) is applied here. The CLFD model calibration for MPS models consists of two steps. In the first step, new conditional realizations are generated using a geostatistical simulation approach with hard data from all wells (including the most recent well). Here we use the multi-scale cross-correlation simulation (MS-CCSIM) geostatistical algorithm (Tahmasebi et al, 2014) for generating these conditional realizations. In the second step, the O-PCA-based RML is applied for conditioning realizations to production data. The schematic of model calibration procedure is presented in Fig. 2.

Figure 2: Schematic of CLFD model calibration step with multipoint geostatistics.

We now briefly describe these procedures. The O-PCA method requires generating realizations of the permeability field conditioned to hard data at well locations. We use in this work. New realizations must be generated at each CLFD step since new conditioning data become available as we proceed in time. The MS-CCSIM algorithm is very efficient and can generate these models in a few minutes. Then the centered matrix of realizations, , is computed,


where is the mean of the realizations. A truncated SVD of is then computed as , where . Given a random -dimensional vector , a new realization can be generated by solving the following optimization problem:


where is a regularization term that is specified such that the realization is generated consistent with the training image (see Vo and Durlofsky (2014) for details). For binary models here, =, where is the unity vector of same length as . The SNOPT algorithm (Gill et al, 2005) is applied to solve (6).

In the O-PCA RML, a calibrated model is generated by minimizing the following objective function.


where corresponds to a projected MPS realization, i.e., . Here is a realization that is unconditioned to production data, but conditioned to hard data ( is also generated by MS-CCSIM). Note that the hard data mismatch term does not appear since hard data are already honored in the realizations and thus in the O-PCA representation. Minimization of Eq. 7 is also performed by SNOPT (Gill et al, 2005). This minimization is to find the optimal -dimensional .

Note that for each trial , the vector is obtained from solving (6). Predicted data, , is then generated by performing a reservoir simulation run using this . The mismatch objective function in Eq. 7 is then computed. The gradient of with respect to is constructed through an adjoint solution, using the automatic differentiation framework (Bukshtynov et al, 2015)

. The gradient is then projected using the chain rule to obtain derivatives with respect to


2.4 Parallelized implementation

The CLFD algorithm can be parallelized in both the optimization and model calibration steps. In the PSO-MADS algorithm, the objective function for PSO particles and the stencil points in MADS algorithm, can be computed simultaneously, by accessing compute nodes. The number of stencil points is twice the number of decision parameters, and is typically selected to be a comparable number (e.g., or ). See Isebor et al (2014a) for more details on parallelization of PSO-MADS.

In the model calibration step, as the RML runs are independent of one another, each RML realization is generated on a separate compute node using distributed computing. All flow simulations are performed using Stanford’s Automatic-Differentiation-based General Purpose Research Simulator (AD-GPRS) (Younis, 2011). The existing OpenMP-based parallelized version of AD-GPRS (Zhou and Tchelepi, 2012) allows us to run each simulation on a computational node with 8 cores. This gives an average speedup of about a factor of 5 for each simulation. For CLFD results presented in this work, at each model calibration step, posterior RML models are generated simultaneously using 50 compute nodes (Fig. 3), providing a speedup factor of 250 compared with a sequential application on a single CPU.

Figure 3: Parallelized implementation of model calibration in CLFD. The models are generated simultaneously by accessing 50 nodes. During iterations of model calibration algorithm, each simulation run is parallelized using 8 cores through OpenMP implementation.

2.5 Statistical assessment of CLFD performance

As discussed in Introduction, recent work that applied modified versions of CLFD (Morosov and Schiozer, 2016; Hanea et al, 2017), reported that application of CLFD, in some cases, has resulted in slight reduction of true-model NPV. In other words, assimilating new data and re-optimizing decision variables over the remaining span of project life, in cases has resulted in decreasing the ultimate outcome, i.e., the NPV of the true-model.

It is essential to note that the outcome of any closed-loop optimization procedure can be treated as a random variable. Therefore, we are able to obtain a distribution of the outcome at each optimization step. This can be accomplished by repeating the closed-loop optimization framework for multiple cases of true-model. The procedure can be outlined as follows.

  1. Select candidate models, randomly, from the prior set.

  2. For each of the true-model cases, perform the closed-loop optimization (CLFD, here) after an optimization over the prior models.

  3. Compute all performance measures for each of cases.

  4. Present the distribution obtained from these runs for each performance measure.

We will define performance measures for CLFD and apply the above procedure in the computational results later in this paper. Note finally that it is essential to perform the optimization over the prior models for each true-model case as noted at Step 2 of the above procedure.

3 Computational results and discussion

In this section computational results are presented for a binary channelized system. In Part 1, detailed application of CLFD is presented. In Part 2, multiple CLFD runs are performed for performance assessment.

3.1 Part 1: CLFD for a channelized model

This example involves a reservoir model described on a two-dimensional uniform grid of dimensions with . A binary channelized training image from Vo and Durlofsky (2014) is taken as the prior geological description, from which a set of unconditional realizations are generated using the MS-CCSIM geostatistical algorithm (Tahmasebi et al, 2014). The realizations are not constrained to honor the sand/shale ratio observed in the training image. The sand permeability is 500 mD, while the shale permeability is 10 mD. The true permeability field, along with an initial guess for the well locations (), is shown in Fig. 4. Three prior realizations of the permeability field are shown in Fig. 5. Porosity is assumed to be constant and equal to 0.2. Initially, the reservoir is at irreducible water saturation which is equal to =0.2.

Reservoir life is 3000 days, which is divided into seven control steps, with the first five control steps of length 180 days, and the last two control steps of length 1050 days. In this example, we assume that two rigs are available (=2) and therefore (a maximum of) two wells can be drilled at each CLFD step. The optimal number of wells, however, is determined from the optimization. The drilling time is specified as =180 days. Therefore, the values are given by . The last optimization is performed at 540 days, which determines the decision parameters corresponding to well type and location of wells 9 & 10 and operational settings of existing wells (wells 1–8). Discount rate is specified to . Simulation and optimization parameters are presented in Table 1, and relative permeability curves are shown in Fig. 6.

Figure 4: True permeability field, with red indicating sand facies (permeability of 500 mD), and blue indicating shale facies (permeability of 10 mD). The initial well configuration is also shown with circles denoting producers and triangles denoting injectors, Example 1.
(a) Realization 1
(b) Realization 2
(c) Realization 3
Figure 5: Three prior realizations of permeability field, with red indicating sand facies (permeability of 500 mD), and blue indicating shale facies (permeability of 10 mD).
Figure 6: Oil and water relative permeability curves.
Parameter Value
$10 MM
$90 STB
$10 STB
$10 STB
Prod BHP range psi
Inj BHP range psi
Table 1: Optimization parameters for all examples

The number of decision parameters is 80. These correspond to 10 categorical variables for well types, 20 integer variables for well locations, and 50 continuous parameters for well settings. Well type can take values in

which corresponds to {producer, do not drill, injector} (for more details please refer to Isebor et al (2014b)). In each optimization, the PSO-MADS algorithm is applied with 60 PSO particles, and the minimum mesh size for MADS is specified to be 1% of the variable range. In the evaluation of a trial point during optimization, if the well distance constraint is violated, simulation run is not performed and a zero objective function value is assigned to that point. This treatment is to avoid frequent failed simulation jobs that happens for trial points with two wells located too close or at the same grid blocks. In addition, economic constraint for each producer is enforced by the simulator by specifying a maximum WOR (water to oil ratio).

First, deterministic optimization is applied using the true-model (shown in Fig 4). The optimal development plan and the final oil saturation map are shown in Fig. 7. The NPV from the initial guess is $441 MM while the optimal NPV is $783 MM. Next, CLFD using OSV is applied for this case.

Figure 7: Well configuration from deterministic optimization (using ), with red denoting producer, blue denoting injector, and the well numbers indicating the drilling sequence. Background shows final oil saturation. Note that two wells are drilled at a time.

CLFD is applied following the workflow shown in Fig. 1. Observed data for CLFD model calibration include production data (oil and water production rates at producers and water injection rate at injectors) measured at 30-day intervals and hard data from all existing wells. Observed production data is generated by adding random noise (measurement error) to the true data, where the true data corresponds to flow results of the true models when it is simulated with the optimal solution. Model calibration is performed every 180 days by first generating new geological realizations conditioned to hard data using the MS-CCSIM algorithm, and then applying O-PCA-based RML using all production data from time 0.

The progress of the true NPV with CLFD step is shown in Fig. 8. In the evaluation of the true NPV, a nominal strategy is applied where a producer is shut in when the cost of handling produced water exceeds the revenue from oil. The NPV from the deterministic optimization is also shown here for comparison. The final truth-case NPV from CLFD is $646.2 MM, which is 36.6% higher than the NPV from optimization (using OSV) over prior realizations.

Figure 8: Optimal expected NPV, and the corresponding NPV for the true-model, versus CLFD step. The number of realizations at each CLFD step is determined using OSV. The star shows the final true NPV from CLFD.

Fig. 9 presents the P10–P50–P90 results for NPV, determined by simulating all 50 realizations and then constructing the cdf, at each CLFD step. The expected NPV based on the current representative subset (which satisfies validation criterion of ) is also displayed. It is evident that the optimal expected NPV for the representative subset falls within the P10–P90 range. Note that in Fig. 9, the percentile values are determined for the updated models () based on current at each CLFD step . Therefore, the uncertainty range may not necessarily shrink as both the models and the solution are evolving.

Figure 9: P10, P50, P90 NPVs evaluated for the entire set of 50 realizations, along with the expected NPV for the representative subset, versus CLFD step.

Fig. 10 shows the evolution of the well configuration and the geological model for two realizations. It is evident that the well scenario involves three injectors at , but four injectors at later times. Each realization continues to show differences though the CLFD steps due to conditioning to new hard and production data.

(a) Realization 1 at with
(b) Realization 1 at with
(c) Realization 1 at with
(d) Realization 2 at with
(e) Realization 2 at with
(f) Realization 2 at with
Figure 10: Evolution of two RML realizations for different CLFD steps, with red indicating sand facies (permeability of 500 mD), and blue indicating shale facies (permeability of 10 mD). Current optimal well configuration and drilling sequence is also depicted. Solid white circles and triangles denote producers and injectors (drilled or in the process of being drilled), and yellow circles and triangles denote planned producers and injectors. Numbers indicate the drilling sequence (Example 1).

It may be worth comparing the final oil saturation from final CLFD solution and from the solution obtained by optimization over prior models through OSV (Fig. 11). It is evident that the solution from CLFD corresponds to improved sweep of oil in this case. Further, the CLFD solution also corresponds to higher oil production as shown in Fig. 12.

Figure 11: Final oil saturation from (a) optimization over prior models by OSV, and (b) final CLFD solution. Corresponding well configuration is also shown, with red denoting producer, blue denoting injector, and the well numbers indicating the drilling sequence.
Figure 12: Cumulative oil production from different approaches.

As discussed earlier, the CLFD model calibration step involves integrating production data from all previous wells (except the most recent wells) together with hard data from all wells including the most recent wells. Integrating both hard data and production data is required to achieve optimal CLFD performance. At each CLFD model calibration step, new realizations conditioned to hard data are generated. We let designate the mean of these realizations (conditioned to all available hard data) at each update step. Figs. 13(a)-13(d) show the evolution of the prior mean with CLFD step. We also compute the mean of posterior realizations (conditioned to both production and hard data). These are shown in Figs. 13(e)-13(h). The optimal development plan at each CLFD step is also shown. It is evident that the mean of the posterior realizations (conditioned to production data) more closely resembles the true-model (Fig. 4). This indicates the importance of integrating production data to reduce the uncertainty in the geological description.

(a) , 180 days
(b) , 360 days
(c) , 540 days
(d) , 720 days
(e) , 180 days
(f) , 360 days
(g) , 540 days
(h) , 720 days
Figure 13: Evolution of mean of () prior realizations (conditioned to hard data) and mean of () posterior realizations of facies distribution, for different CLFD steps. Current optimal well configuration and drilling sequence is also depicted. White circles and triangles denote producers and injectors, respectively. Wells with colored (red or blue) numbers are drilled, while outlined red circles and blue triangles denote planned producers and injectors. For the prior model (a-d) only the drilled wells are shown. Numbers indicate the drilling sequence.

In terms of computational cost, this CLFD application involved a total number of 356,500 simulation runs. As computations were performed on a cluster with access to maximum 400 compute cores, the parallelized computations is equivalent to 1,151 sequential simulation runs (calls to cluster) where each batch of simulations takes about 104 seconds on average. The total computational time is about 32 hours. Note that in practice, each CLFD step is performed at a different time and thus the computational cost will be a fraction of the total cost here.

3.2 Part 2: Statistical assessment of CLFD performance

In this section, our objective is to assess both the CLFD and the OSV performance through application with multiple truth-case models. Specifically, we consider 24 additional models, and we repeat the CLFD application for each of these true-models. We then assess the performance in terms of various measures as follows.

3.2.1 Assessment of optimization with sample validation performance

We first discuss the assessment of OSV performance, i.e., true NPV improvement obtained from optimization with sample validation (OSV) relative to optimization over representative models. We only consider the optimization over prior models () here. The results for five of the =25 cases are shown in Table 2 (the first row corresponds to the true-model discussed earlier). Only for 60% of cases, the true NPV obtained from OSV (third column of Table 2) is greater than that obtained from optimization with representative realizations (in three cases, was adequate to achieve ). These results may seem counter intuitive, as the OSV procedure, which involves optimization over a larger set to achieve , does not improve the true model NPV for about 40% of the cases. However, statistics as the mean and the median indicate significant improvement obtained through OSV, as discussed next.

We define as true NPV improvement obtained from OSV (for prior models) relative to optimization over prior models, as


where is the optimal solution obtained form the first optimization in OSV procedure (i.e., from optimization over representative prior models). Values of are shown in Table 2, with the CDF plot shown in Fig. 14. On average, OSV improves the true NPV by 64% over the NPV obtained from optimization with prior models, while the median improvement is 16.8%. In other words, although there is a chance that the true NPV may decrease from OSV (relative to optimization over models), the potential NPV improvement from OSV is significant in terms of average (64%) and median (16.8%) measures. This can also be observed in the CDF plot (Fig. 14) which demonstrates greater values for positive improvements. The OSV procedure is therefore statistically desirable.

The OSV procedure systematically determines adequate number of representative models in optimization over a large set of realizations. The average and median in this case are 12.5 and 9, respectively (Table 2). We conclude that for this problem, is typically not sufficient, while optimization over or realizations is more likely to satisfy the OSV validation criterion.

True model from opt with from OSV final
1 474 473 0.00 9
6 462 376 -0.19 16
11 486 461 -0.05 9
16 402 553 0.38 9
21 318 323 0.02 9
Avg (1-25) 339 392 0.640 12.5
Median (1-25) 0.168 9
Table 2: NPV values ($ MM) from optimization over prior realizations and by use of OSV (where is increased to satisfy ) for prior realizations, and the values of for five different true models. Average and median values are also provided.
Figure 14: Cumulative distribution function (CDF) for 25 values (relative improvement in true-model NPV from OSV) for (prior models).

3.2.2 Assessment of single-step CLFD performance

We now assess the performance of CLFD framework in terms of improvement of true-model NPV from a single CLFD step. The NPV values for 5 models at each CLFD step are compiled in Table 3. The average true NPV from 25 runs (also shown in Table 3) increases monotonically with CLFD step. The average true NPV at CLFD step 1 is 9.9% greater than optimization over prior (), and the increase from step 1 to 2 and from step 2 to 3 is 7.4% and 10.1%, respectively. At each CLFD step  (), we define the true-model NPV improvement for true-model as


where the numerator is the NPV improvement for true-model as a result of performing one step of CLFD (step ). The mean and percentile values of are compiled in Table 4. The mean improvement varies between 10–14%, and the P10 monotonically increase with CLFD step (Table 4). CLFD improved the true-model NPV for 64% of cases at , 76% of cases at , and 80% of cases at step . This performance improvement can also be observed visually in Fig. 15. Therefore, we observe that although application of a single-step CLFD may even reduce the true NPV of the project in some cases, the chance of increasing true NPV increases at later steps.

True model OSV-prior (=0) =180 (=1) =360 (=2) =540 (=3)
1 473 488 577 646
6 376 376 314 406
11 461 289 449 570
16 553 498 529 532
21 323 518 596 602
Avg (1-25) 392 431 463 510
Table 3: NPV values ( in $ MM) from OSV over prior models and at each CLFD step. Average values are also provided.
Figure 15: Cumulative distribution function (CDF) for 25 values (relative improvement in true-model NPV from one-step of CLFD run) for =1 (=180 days), =2 (=360 days), and =3 (=540 days).
CLFD step (time) mean P10 P50 P90 fraction positive
(180 days) 0.14 -0.18 0.11 0.55 0.64
(360 days) 0.10 -0.16 0.10 0.42 0.76
(540 days) 0.13 -0.01 0.06 0.46 0.80
Table 4: Statistics of (true NPV improvement after one step of CLFD) for 25 different true models

At ( days), the uncertainty reduction in model description through integration of spatial data (hard data from four wells) and temporal data (production data from two wells) is not as significant as later steps and by optimizing over the updated models, for almost 36% of cases, the new optimal solution resulted in a lower true NPV of the project. At later steps, however, more wells are drilled and more data results in improved model description. As a result there is a higher chance that the optimization over the updated models would result in a higher project outcome (true NPV).

3.2.3 Assessment of multi-steps CLFD performance

To further assess the overall net impact of performing multiple steps of CLFD, we define the relative improvement with respect to optimization over prior as


where rel-pr stands for “relative to prior”. For =1, in Eq. 10 provides the same quantity as that in Eq. 9. However, at later steps (=2 and =3), these two quantities are different.

Fig. 16 shows the CDF plot for values, and Table 5 provides the numeric percentiles. These results indicate that application of two CLFD steps, has improved the NPV of 68% of the models, with an average improvement of 23%. The full CLFD run (three CLFD steps in this case), improved the true NPV in 96% of cases with an average improvement of 37%.

These results together with those in the previous section (and Table 4), indicate that the probability of increasing true model NPV increases with CLFD step. In other words, there is a higher chance of improving true NPV through multi-step CLFD, compared with a single step of CLFD. In this case, we observe that application of one, two, and three consecutive steps of CLFD improved the NPV in 64%, 68%, and 96% of case, respectively. These results infer about an implicit correction mechanism in CLFD that decisions about planned wells (those to be drilled at later steps in the time horizon) are ultimately improved.

Figure 16: Cumulative distribution function (CDF) for 25 values (improvement in true-model NPV relative to optimization over prior models).
CLFD step (time) mean P10 P50 P90 fraction positive
(360 days) 0.23 -0.17 0.15 0.80 0.68
(540 days) 0.37 0.08 0.24 0.85 0.96
Table 5: Statistics of for 25 different true models

The average percentiles and the average expected NPV, for these 25 runs, versus the CLFD step is shown in Fig 17 and numerical values are provided in Table 6. Each point in the figure (and table) is an average of 25 values (results in Fig. 9 shows values for run =1), e.g., = where is the optimal solution at CLFD step for run , and is true-model . These results indicate a monotonic increase in all average values. Further, the difference between the average P10 and the average P90 shrinks with CLFD step, although this does not necessarily hold for each run. These results are not specific to a particular true-model, and therefore they demonstrate the overall performance of CLFD and present an assessment of project value under optimal development.

Figure 17: Average values of P10, P50, P90 NPVs, over 25 runs along with the average true NPV, versus CLFD step.
OSV-prior (=0) =180 (=1) =360 (=2) =540 (=3)
521 558 561 588
397 449 464 504
246 326 358 415
392 431 463 510
Table 6: Numerical values of average P10, P50, P90 NPVs, over 25 runs along with the average true NPV, versus CLFD step.

The computational experiment of performing 25 CLFD applications took about 789 hours (33 days). The total number of simulations is 9,333,900 which are performed through 30,136 calls to a computer cluster with access to maximum 400 cores.

3.3 Discussion

The results for performance assessment showed a high range of possible outcomes from CLFD runs. The highly nonlinear nature of the optimization step is perhaps the main cause of the wide range of variability in the ultimate outcome. Although, the application of PSO component in the PSO-MADS algorithm, provides some degree of global exploration, the optimization may ultimately converge to a locally optimal solution (for the representative subset) and there is no guarantee of finding the global optimum. Results presented in Onwunalu and Durlofsky (2010) show that for a problem where the location of a single well is optimized, the objective function surface is highly rough with many local optima. The field development optimization problem in CLFD is significantly more complex, and a significant number of local optima may exist. Another main cause to the wide range of final outcome, is the sparsity in spatial data and the inherent uncertainty in model description after model calibration. Despite all the complexity, application of the full CLFD improved the true NPV in 96% of cases, showing the effectiveness of the framework.

The statistical assessment can also be applied to test and compare the effect of different treatments. Different geostatistical algorithms or model calibration methods can be applied in CLFD. In order to compare two different treatments, multiple CLFD runs should be performed and results can be analyzed in terms of a performance measure, e.g., average true-model improvement.

4 Conclusions and future work

In this work, a closed-loop field development (CLFD) optimization framework is developed for models described by multipoint geostatistics (MPS). This framework involves three major steps: 1) optimization of the full development plan for the entire project life, 2) drilling and completing new wells and collecting new data, and 3) updating models based on all data. The optimization step is performed over multiple models through optimization with sample validation (OSV). A two-step model calibration procedure is presented for model calibration where spatial data (hard data at well locations) are integrated through a MPS simulation method, and production data are integrated through a gradient-based RML procedure using O-PCA parameterization.

The new methodology is applied 25 times for 25 different true-model cases, and the OSV and the CLFD performance were assessed in terms of true NPV improvement. The median true NPV improvement from OSV was 17% and the mean improvement was 64%, indicating the importance of performing optimization on an adequately large set. The CLFD performance is assessed in terms of true NPV improvement. Results indicated that the true NPV in 64%, 68% and 96% of cases increased after one, two and three steps, respectively. These results showed that application of multi-step CLFD is much more likely to increase true NPV than a single step CLFD, particularly as richer data appears at later CLFD steps. The statistical assessment presented in this work, demonstrates the importance of performing multiple runs and considering several true-model cases to assess the performance of any closed-loop optimization framework.

There are several directions to investigate in future. In the model calibration step, a single training image was considered. It will be of interest to investigate the CLFD performance in cases where multiple training images with assigned probabilities exist. Application of CLFD optimization to geothermal reservoir, groundwater remediation (Bayer et al, 2008; Ghorbanidehno et al, 2017), and CO-EOR problems should also be investigated. In this work, a fixed project life was specified in the optimization step of CLFD and the expected NPV was optimized. It will be of interest to assess CLFD performance when the rate of return is incorporated in optimization objective and the economic project life is also optimized (Shirangi et al, 2017).


This research was conducted at Stanford University between 2015 and 2017. The author acknowledges helpful discussions with Prof. Louis Durlofsky. Access to the computational resources provided by the Stanford Center for Computational Earth & Environmental Science (CEES) is gratefully acknowledged. Financial support of this work was partially provided by the Stanford Smart Fields Consortium.


  • Aitokhuehi and Durlofsky (2005) Aitokhuehi I, Durlofsky LJ (2005) Optimizing the performance of smart wells in complex reservoirs using continuously updated geological models. Journal of Petroleum Science and Engineering 48(3):254–264
  • Arnold et al. (2016) D. Arnold, V. Demyanov, M. Christie, A. Bakay, and K. Gopa. Optimisation of decision making under uncertainty throughout field lifetime: A fractured reservoir example. Computers & Geosciences, 95:123–139, 2016.
  • Bayer and Finkel (2007) Bayer P, Finkel M (2007) Optimization of concentration control by evolution strategies: Formulation, application, and assessment of remedial solutions. Water resources research 43(2)
  • Bayer et al (2008) Bayer P, Bürger C, Finkel M (2008) Computationally efficient stochastic optimization using multiple realizations. Advances in Water Resources 31(2):399–417
  • Bayer et al (2010) Bayer P, de Paly M, Bürger CM (2010) Optimization of high-reliability-based hydrological design problems by robust automatic sampling of critical model realizations. Water Resources Research 46(5)
  • Bellout et al (2012) Bellout MC, Ciaurri DE, Durlofsky LJ, Foss B, Kleppe J (2012) Joint optimization of oil well placement and controls. Journal of Mathematical Modelling and Numerical Optimisation 16(4):1061–1079
  • Bjarkason et al (2017) Bjarkason EK, Maclaren OJ, O’Sullivan JP, O’Sullivan MJ, (2017) Randomized Truncated SVD Levenberg-Marquardt Approach to Geothermal Natural State and History Matching. preprint arXiv:1710.01031.
  • Bouzarkouna et al (2012) Bouzarkouna Z, Ding DY, Auger A (2012) Well placement optimization with the covariance matrix adaptation evolution strategy and meta-models. Computational Geosciences 16(1):75–92
  • Bukshtynov et al (2015) Bukshtynov V, Volkov O, Durlofsky LJ, Aziz K (2015) Comprehensive framework for gradient-based optimization in closed-loop reservoir management. Computational Geosciences 19(4):877–897
  • Caers (2003) J. Caers. History matching under training-image-based geological model constraints. SPE Journal, 8(3):218–226, 2003.
  • Cameron and Durlofsky (2014) D. A. Cameron and L. J. Durlofsky. Optimization and data assimilation for geological carbon storage. In: R. Al-Khoury, J. Bundschuh (Eds.), Computational Models for CO Geo-sequestration & Compressed Air Energy Storage. CRC Press, 355–-388, 2014.
  • Cameron and Durlofsky (2012) Cameron DA, Durlofsky LJ (2012) Optimization of well placement, CO injection rates, and brine cycling for geological carbon sequestration. International Journal of Greenhouse Gas Control 10:100–112
  • Chen et al (2009) Chen Y, Oliver DS, Zhang D, et al (2009) Efficient ensemble-based closed-loop production optimization. SPE Journal 14(04):634–645
  • Dickstein et al (2017) Dickstein F, Goldfeld P, Pfeiffer GT, Pinto RV (2017) Truncated conjugate gradient and improved LBFGS and TSVD for history matching. Computational Geosciences, DOI
  • Ghorbanidehno et al (2017)

    Ghorbanidehno H, Kokkinaki A, Kitanidis PK, Darve E (2017) Optimal estimation and scheduling in aquifer management using the Rapid Feedback Control Method. Advances in Water Resources, DOI
  • Gill et al (2005) Gill PE, Murray W, Saunders MA (2005) SNOPT: An SQP algorithm for large-scale constrained optimization. SIAM review 47(1):99–131
  • Hanea et al (2017) Hanea RG, Casanova P, Hustoft L, Bratvold RB, et al (2017) Drill and Learn: A Decision Making Workflow to Quantify Value of Learning. Paper SPE-182719-MS presented at SPE Reservoir Simulation Conference, Houston, Texas, DOI
  • Hidalgo et al (2017) Hidalgo DM, Emerick AA, Couto P, et al (2017) Closed-loop field development under geological uncertainties: application in a Brazilian benchmark case. Paper OTC-28089-MS presented at Offshore Technology Conference, Rio de Janeiro, Brasil.
  • Isebor et al (2014a) Isebor OJ, Durlofsky LJ, Ciaurri DE (2014a) A derivative-free methodology with local and global search for the constrained joint optimization of well locations and controls. Computational Geosciences (18):463–482
  • Isebor et al (2014b) Isebor OJ, Echeverría Ciaurri D, Durlofsky LJ (2014b) Generalized field-development optimization with derivative-free procedures. SPE Journal (19):891–908
  • Jahandideh and Jafarpour (2016) Jahandideh A, Jafarpour B (2016) Optimization of hydraulic fracturing design under spatially variable shale fracability. Journal of Petroleum Science and Engineering 138:174–188
  • Jansen et al (2009) Jansen JD, Douma SD, Brouwer DR, den Hof PMJV, Heemink AW (2009) Closed-loop reservoir management. Paper SPE 119098 presented at SPE Reservoir Simulation Symposium, The Woodlands, Texas.
  • Kitanidis (1995) Kitanidis PK (1995) Quasi-linear geostatistical theory for inversing. Water Resources Research 31(10):2411–2419
  • Li and Jafarpour (2012) Li L, Jafarpour B (2012) A variable-control well placement optimization for improved reservoir development. Computational Geosciences 16(4):871–889
  • Mariethoz and Caers (2014) Mariethoz G, Caers J (2014) Multiple-Point Geostatistics: Stochastic Modeling with Training Images. John Wiley & Sons.
  • Morosov and Schiozer (2016) Morosov AL, Schiozer DJ (2016) Field Development Process Revealing Uncertainty Assessment Pitfalls. Paper SPE-180094-MS presented at the SPE Europec featured at 78th EAGE Conference, Vienna, Austria. DOI
  • Mustapha and Dimitrakopoulos (2011) Mustapha H, and Dimitrakopoulos R (2011) HOSIM: A high-order stochastic simulation algorithm for generating three-dimensional complex geological patterns. Computers & geosciences, 37(9), pp.1242-1253. DOI
  • Nwankwor et al (2013)

    Nwankwor E, Nagar AK, Reid D (2013) Hybrid differential evolution and particle swarm optimization for optimal well placement. Computational Geosciences 17(2):249–268

  • Oliver et al (1996) Oliver DS, He N, Reynolds AC (1996) Conditioning permeability fields to pressure data. In: Proccedings of the European Conference for the Mathematics of Oil Recovery
  • Onwunalu and Durlofsky (2010) Onwunalu JE, Durlofsky L (2010) Application of a particle swarm optimization algorithm for determining optimum well location and type. Computational Geosciences 14(1): 183–198.
  • Ozdogan and Horn (2006) Ozdogan U, Horn RN (2006) Optimization of well placement under time-dependent uncertainty. SPE Journal 9(2):135–145
  • Regtien (2010) J. M. Regtien. Extending the smart fields concept to enhanced oil recovery. Paper SPE 136034 presented at SPE Russian Oil and Gas Conference, Moscow, Russia, 2010.
  • Sarma et al (2006) Sarma P, Durlofsky L, Aziz K, Chen W (2006) Efficient real-time reservoir management using adjoint-based optimal control and model updating. Computational Geosciences 10:3–36
  • Sarma et al (2008) Sarma P, Durlofsky LJ, Aziz K (2008) Kernel principal component analysis for efficient differentiable parameterization of multipoint geostatistics. Mathematical Geosciences 40:3–32
  • Shirangi (2017) Shirangi MG (2017) Advanced Techniques for Closed-Loop Reservoir Optimization under Uncertainty. PhD Thesis, Stanford University.
  • Shirangi (2013) Shirangi MG (2013) Closed-loop Field Development Optimization. Proceedings of the 26th Annual Meeting of Stanford Center for Reservoir Forecasting, Palo Alto, California.
  • Shirangi (2014) Shirangi MG (2014) History matching production data and uncertainty assessment with an efficinet TSVD parameterization algorithm. Journal of Petroleum Science and Engineering 113:54–71, DOI
  • Shirangi (2011) Shirangi MG (2011) History Matching Production Data with Truncated SVD Parameterization. Master’s thesis, University of Tulsa.
  • Shirangi and Durlofsky (2015) Shirangi MG, Durlofsky LJ (2015) Closed-loop field development under uncertainty by use of optimization with sample validation. SPE Journal 20(05):908 – 922, DOI
  • Shirangi and Durlofsky (2016) Shirangi MG, Durlofsky LJ (2016) A general method to select representative models for decision making and optimization under uncertainty. Computers & Geosciences 96:109–123. DOI
  • Shirangi and Emerick (2016) Shirangi MG, Emerick AA (2016) An improved TSVD-based Levenberg-Marquardt algorithm for history matching and comparison with Gauss-Newton. Journal of Petroleum Science and Engineering 143:258–271, DOI
  • Shirangi and Mukerji (2012) Shirangi MG, Mukerji T (2012) Retrospective optimization of well controls under uncertainty using kernel clustering. Proceedings of the 25th Annual Meeting of the Stanford Center for Reservoir Forecasting, Monterey, California.
  • Shirangi et al (2017) Shirangi MG, Volkov O, Durlofsky LJ (2017) Joint optimization of economic project life and well controls. SPE Journal, DOI
  • Tahmasebi et al (2014) Tahmasebi P, Sahimi M, Caers J (2014) MS-CCSIM: accelerating pattern-based geostatistical simulation of categorical variables using a multi-scale search in Fourier space. Computers & Geosciences 67:75–88
  • Tarantola (2005) Tarantola A (2005) Inverse Problem Theory and Methods for Model Parameter Estimation. SIAM, Philadelphia, USA
  • Vo and Durlofsky (2014) Vo HX, Durlofsky LJ (2014) A new differentiable parameterization based on principal component analysis for the low-dimensional representation of complex geological models. Mathematical Geosciences 46(7):775–813
  • Vo and Durlofsky (2015) Vo HX, Durlofsky LJ (2015) Data assimilation and uncertainty assessment for complex geological models using a new PCA-based parameterization. Computational Geosciences 19(4):747–767
  • Vo and Durlofsky (2016) Vo HX, Durlofsky LJ (2016) Regularized kernel PCA for the efficient parameterization of complex geological models. Journal of Computational Physics 322:859–881, DOI
  • Wilson and Durlofsky (2013) Wilson KC, Durlofsky LJ (2013) Optimization of shale gas field development using direct search techniques and reduced-physics models. Journal of Petroleum Science and Engineering 108:304–315
  • Younis (2011) Younis MR, 2011. Modern Advances in Software and Solution Algorithms for Reservoir Simulation. Ph.D. thesis, Stanford University.
  • Zhang et al (2017) Zhang Y, Lu R, Forouzanfar F, Reynolds AC (2017) Well placement and control optimization for WAG/SAG processes using ensemble-based method. Computers & Chemical Engineering 101: 193–209
  • Zhou and Tchelepi (2012) Zhou Y, Tchelepi H (2012) Multi-core and GPU parallelization of a general purpose reservoir simulator. Proceedings of the 13th European Conference on the Mathematics of Oil Recovery, Biarritz, France.