Abstract
Closedloop 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 optimizationbased PCA. A statistical procedure is presented to assess the performance of CLFD. Methodology is applied to an oil reservoir example for 25 different truemodel cases. Application of a singlestep 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 physicsbased 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 closedloop 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 closedloop optimization may differ on type of decision parameters included in these two steps.
Closedloop 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 twopoint 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 closedloop field development methodology is developed for reservoir models described by multipoint geostatistics. A twostep 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 closedloop 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 closedloop 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 waterflooding 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 (CMAES) (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 gradientbased 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 ensemblebased 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 PSOMADS 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 illposed and the number of unknown model parameters can be very large. The illposedness 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 PCAbased parameterization (OPCA) that enables application of efficient gradientbased approaches for model calibration of complex channelized systems. This optimizationbased principal component analysis (OPCA) 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 truemodel NPV in most cases; however, in at least one case in each work, the truemodel 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 twopoint 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 twodimensional 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: closedloop field development optimization
Closedloop 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.
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 truemodel and adding random noise to the true data. The truemodel 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
(1) 
where
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(2) 
Computing NPV value for each model involves a flow simulation run. In this work, the waterflooding process with twophase flow is considered where NPV is computed as
(3) 
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 PSOMADS optimization algorithm (Isebor et al, 2014b). Nonlinear constraint here is the minimum well distance constraint which is handled through the filter method in PSOMADS. The PSOMADS 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.,
(4) 
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.

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 .

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

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

Compute in Eq. 4 by replacing with and with .

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 flowbased and permeabilitybased 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 (historymatched) 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 optimizationbased principal component analysis (OPCA) 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 multiscale crosscorrelation simulation (MSCCSIM) geostatistical algorithm (Tahmasebi et al, 2014) for generating these conditional realizations. In the second step, the OPCAbased RML is applied for conditioning realizations to production data. The schematic of model calibration procedure is presented in Fig. 2.
We now briefly describe these procedures. The OPCA 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 MSCCSIM algorithm is very efficient and can generate these models in a few minutes. Then the centered matrix of realizations, , is computed,
(5) 
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:
(6) 
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 OPCA RML, a calibrated model is generated by minimizing the following objective function.
(7) 
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 MSCCSIM). Note that the hard data mismatch term does not appear since hard data are already honored in the realizations and thus in the OPCA 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 PSOMADS 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 PSOMADS.
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 AutomaticDifferentiationbased General Purpose Research Simulator (ADGPRS) (Younis, 2011). The existing OpenMPbased parallelized version of ADGPRS (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.
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 truemodel NPV. In other words, assimilating new data and reoptimizing decision variables over the remaining span of project life, in cases has resulted in decreasing the ultimate outcome, i.e., the NPV of the truemodel.
It is essential to note that the outcome of any closedloop 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 closedloop optimization framework for multiple cases of truemodel. The procedure can be outlined as follows.

Select candidate models, randomly, from the prior set.

For each of the truemodel cases, perform the closedloop optimization (CLFD, here) after an optimization over the prior models.

Compute all performance measures for each of cases.

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 truemodel 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 twodimensional 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 MSCCSIM 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.
Parameter  Value 

$10 MM  
$90 STB  
$10 STB  
$10 STB  
Prod BHP range  psi 
Inj BHP range  psi 
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 PSOMADS 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 truemodel (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.
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 30day 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 MSCCSIM algorithm, and then applying OPCAbased 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 truthcase NPV from CLFD is $646.2 MM, which is 36.6% higher than the NPV from optimization (using OSV) over prior realizations.
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.
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.
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.
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 truemodel (Fig. 4). This indicates the importance of integrating production data to reduce the uncertainty in the geological description.
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 truthcase models. Specifically, we consider 24 additional models, and we repeat the CLFD application for each of these truemodels. 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 truemodel 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
(8) 
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 (125)  339  392  0.640  12.5 
Median (125)  –  –  0.168  9 
3.2.2 Assessment of singlestep CLFD performance
We now assess the performance of CLFD framework in terms of improvement of truemodel 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 truemodel NPV improvement for truemodel as
(9) 
where the numerator is the NPV improvement for truemodel 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 truemodel 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 singlestep 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  OSVprior (=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 (125)  392  431  463  510 
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 
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 multisteps 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
(10) 
where relpr 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 multistep 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.
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 
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 truemodel . 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 truemodel, and therefore they demonstrate the overall performance of CLFD and present an assessment of project value under optimal development.
OSVprior (=0)  =180 (=1)  =360 (=2)  =540 (=3)  

521  558  561  588  
397  449  464  504  
246  326  358  415  
392  431  463  510 
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 PSOMADS 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 truemodel improvement.
4 Conclusions and future work
In this work, a closedloop 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 twostep 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 gradientbased RML procedure using OPCA parameterization.
The new methodology is applied 25 times for 25 different truemodel 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 multistep 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 truemodel cases to assess the performance of any closedloop 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 COEOR 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).
Acknowledgments
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.
References
 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 highreliabilitybased 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 LevenbergMarquardt 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 metamodels. Computational Geosciences 16(1):75–92
 Bukshtynov et al (2015) Bukshtynov V, Volkov O, Durlofsky LJ, Aziz K (2015) Comprehensive framework for gradientbased optimization in closedloop reservoir management. Computational Geosciences 19(4):877–897
 Caers (2003) J. Caers. History matching under trainingimagebased 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. AlKhoury, J. Bundschuh (Eds.), Computational Models for CO Geosequestration & 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 ensemblebased closedloop 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 https://doi.org/10.1007/s1059601796944

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
https://doi.org/10.1016/j.advwatres.2017.10.011  Gill et al (2005) Gill PE, Murray W, Saunders MA (2005) SNOPT: An SQP algorithm for largescale 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 SPE182719MS presented at SPE Reservoir Simulation Conference, Houston, Texas, DOI https://doi.org/10.2118/182719MS
 Hidalgo et al (2017) Hidalgo DM, Emerick AA, Couto P, et al (2017) Closedloop field development under geological uncertainties: application in a Brazilian benchmark case. Paper OTC28089MS presented at Offshore Technology Conference, Rio de Janeiro, Brasil.
 Isebor et al (2014a) Isebor OJ, Durlofsky LJ, Ciaurri DE (2014a) A derivativefree 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 fielddevelopment optimization with derivativefree 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) Closedloop reservoir management. Paper SPE 119098 presented at SPE Reservoir Simulation Symposium, The Woodlands, Texas.
 Kitanidis (1995) Kitanidis PK (1995) Quasilinear geostatistical theory for inversing. Water Resources Research 31(10):2411–2419
 Li and Jafarpour (2012) Li L, Jafarpour B (2012) A variablecontrol well placement optimization for improved reservoir development. Computational Geosciences 16(4):871–889
 Mariethoz and Caers (2014) Mariethoz G, Caers J (2014) MultiplePoint 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 SPE180094MS presented at the SPE Europec featured at 78th EAGE Conference, Vienna, Austria. DOI https://doi.org/10.2118/180094MS
 Mustapha and Dimitrakopoulos (2011) Mustapha H, and Dimitrakopoulos R (2011) HOSIM: A highorder stochastic simulation algorithm for generating threedimensional complex geological patterns. Computers & geosciences, 37(9), pp.12421253. DOI https://doi.org/10.1016/j.cageo.2010.09.007

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 timedependent 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 realtime reservoir management using adjointbased 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 ClosedLoop Reservoir Optimization under Uncertainty. PhD Thesis, Stanford University.
 Shirangi (2013) Shirangi MG (2013) Closedloop 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 https://doi.org/10.1016/j.petrol.2013.11.025
 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) Closedloop field development under uncertainty by use of optimization with sample validation. SPE Journal 20(05):908 – 922, DOI https://doi.org/10.2118/173219PA
 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 https://doi.org/10.1016/j.cageo.2016.08.002
 Shirangi and Emerick (2016) Shirangi MG, Emerick AA (2016) An improved TSVDbased LevenbergMarquardt algorithm for history matching and comparison with GaussNewton. Journal of Petroleum Science and Engineering 143:258–271, DOI https://doi.org/10.1016/j.petrol.2016.02.026
 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 https://doi.org/10.2118/182642PA
 Tahmasebi et al (2014) Tahmasebi P, Sahimi M, Caers J (2014) MSCCSIM: accelerating patternbased geostatistical simulation of categorical variables using a multiscale 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 lowdimensional 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 PCAbased 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 https://doi.org/10.1016/j.jcp.2016.07.011
 Wilson and Durlofsky (2013) Wilson KC, Durlofsky LJ (2013) Optimization of shale gas field development using direct search techniques and reducedphysics 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 ensemblebased method. Computers & Chemical Engineering 101: 193–209
 Zhou and Tchelepi (2012) Zhou Y, Tchelepi H (2012) Multicore and GPU parallelization of a general purpose reservoir simulator. Proceedings of the 13th European Conference on the Mathematics of Oil Recovery, Biarritz, France.
Comments
There are no comments yet.