Regular Article
A new approach with multiple realizations for image perturbation using cosimulation and probability perturbation method
University of Campinas, PO Box 6122, Campinas, São Paulo
13083970, Brazil
^{*} Corresponding author: goncalo@dep.fem.unicamp.br
Received:
13
December
2017
Accepted:
18
September
2018
History matching is an inverse problem with multiple possible answers. The petrophysical properties of a reservoir are highly uncertain because data points are scarce and widely scattered. Some methods reduce uncertainty in petrophysical characterization; however, they commonly use a single matched model as a reference, which may excessively reduce uncertainty. Choosing a single image may cause the model to converge to a local minimum, yielding less reliable history matching. This work improves on the history matching presented by Oliveira et al. ((2017a) J. Petrol. Sci. Eng. 153, 111–122) using a benchmark model (UNISIMIH based on the Namorado field in Brazil). We use a new approach for a Probability Perturbation Method and image perturbation using CoSimulation. Instead of using a single image as the reference, a set of best images is used to increase variability in the properties of the reservoir model while matching production data with history data. This approach mitigates the risk of the potentially excessive reduction of uncertainties that can happen when using a single model. Our methodology also introduces a new objective function for water breakthrough, improving model quality because of the importance of matching the water breakthrough in the process. Our proposed methodology for image perturbation uses the UNISIMIH, which comprises 25 wells and has 11 years of history data. Our methodology made the process of calibration more effective than the history matching by Oliveira et al. ((2017a) J. Petrol. Sci. Eng. 153, 111–122). Crossinfluence was minimized, making the history matching more objective and efficient, and consequently, the production forecasts more reliable.
© G. Soares Oliveira et al., published by IFP Energies nouvelles, 2018
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
An essential step in reservoir modeling is the characterization of petrophysical properties in the reservoir. Due to scarce data and reservoir complexity, a large uncertainty is present in the process, resulting in wide production variability. In order to minimize the misfit between production and history data a calibration is necessary with a process known as history matching or data assimilation. The resulting matched models can then be used production forecast and decision analysis.
History matching is an iterative process that changes the reservoir model with the objective of decreasing the misfit between production and history data. This misfit can be measured using the concept of objective function, that represents the normalized misfit between simulated and history data, allowing the comparison between different dynamic data types (e.g. pressure and rates).
While deterministic history matching creates a single model, probabilistic history matching creates a set of models that, although fitted to the history data, may provide different production forecasts. The central idea of probabilistic approach is to use observed data to condition the process of reducing reservoir parameter uncertainties, such as porosity, permeability, for example. Due to the multiple possible answers to match the data, a single model (deterministic approach) is no longer a goal and instead, a probabilistic history matching is preferable in order to allow uncertainty quantification which will help reservoir management and decision making.
A single objective function, normally the sum or average of all objective functions is not enough to capture specific behavior, and instead a multiobjective function is preferable. This approach allows the evaluation of each objective function and well individually (Almeida et al., 2014; Caeiro et al., 2014; Christie et al., 2013; Hajizadeh et al., 2011; Hutahaean et al., 2015; Mesquita et al., 2015; SchulzeRiegert and Kroschem, 2007). This improves understanding of reservoir behavior and enables the identification of necessary changes to history match a particular objective function (Oliveira et al., 2017a).
Several works have presented reservoir image perturbation methods that improve the geological consistency, the best known of which are the Probability Perturbation Method (Hoffman, 2005; Hoffman and Caers, 2003, 2005) for categorical properties such as facies, Gradual Deformation Method (Ding and Roggero, 2010; GervaisCouplet et al., 2007; Hu, 2000; Roggero and Hu, 1998) and CoSimulation (MataLima, 2008; MataLima and Soares, 2007; Oliveira, 2014) for continuous properties such as porosity and permeability.
These methods have the particularity of allowing a local perturbation within defined nonoverlapping regions, avoiding the perturbation of the entire reservoir equally. Using the local perturbation has three main advantages: (1) the reservoir images have an increase in variability, (2) the method is more flexible to achieve the match and (3) the convergence is twice as fast than using global perturbation (Hoffman and Caers, 2005; Oliveira, 2014; Oliveira et al., 2017b).
Oliveira et al. (2017a) proposed a methodology that was able to deal with multiple petrophysical properties by locally reducing the uncertainty in facies, porosity, and permeability simultaneously. It was noticed that using a single image could generate multiple solutions but with similar characteristics (local minimum) making the history matching process less effective because of the crossinfluence between regions and properties. Occasionally, it was necessary to undo the perturbation to correct the misfit by selecting a new reservoir image to perturb the model again.
For continuous properties perturbation, there are several works, using the gradual deformation method to combine multiple reservoir images to generate a new image (Ding and Roggero, 2010; GervaisCouplet and Le Ravalec, 2018; GervaisCouplet et al., 2007). However in cosimulation this approach has not been used, and for Probability Perturbation (PPM), the method is actually limited to use a single image.
Oliveira et al. (2017a) proposed using a multivariate sensitivity analysis to identify the regions and properties that impact each objective function as well as how to perturb an image to history match dynamic data.
The purpose of probabilistic history matching is to efficiently explore the search space and find different models that fit the history data by combining the uncertain parameters. However, in image perturbation, whether for categorical or continuous properties, a single reservoir image is usually chosen as the reference and conditioned in the following iterations. This contradicts the concept of probabilistic history matching because the models created will tend to a set of solutions with a low variability that is highly dependent on the chosen image.
The objective of this work is to improve the efficiency of history matching by avoiding the convergence to a local minimum. We achieve this using a set of models as input for image perturbation rather than a single image. We adapt the probability perturbation method and cosimulation for use in our proposed methodology.
This work also uses an indicator to quantify the misfit for water time breakthrough to further improve history matching.
Oliveira et al. (2017a) efficiently explored the advantage of using multiobjective functions with a multivariate sensitivity analysis and a local perturbation method. The multiobjective functions allow problems diagnosis by identifying the most critical dynamic data for the history matching. On the other hand, the multivariate sensitivity analysis supports the user to understand how the reservoir image (which regions and properties) should be perturbed and the local perturbation (either using PPM or CoSGS) allows the perturbation of specific regions and properties.
This work starts by presenting the necessary modifications to the Probability Perturbation Method (Sect. 2.1) and CoSimulation (Sect. 2.2) in order to use multiple reservoir images. Following in Section 3 is a brief explanation of the general workflow developed in Oliveira et al. (2017a) to update the reservoir image and the most relevant topics for the method to quantify match quality, using the concept of NQDS (Sect. 3.1) and the multivariate sensitivity analysis for petrophysical properties (Sect. 3.2). In Section 4, the proposed method and modifications are applied to a complex case study, with a complete history matching of the 25 wells. A comparison between history matching and production forecast is done for the cases using a single and multiple images.
2 Image perturbation methods
During the process of history matching, a problem occasionally occurs: previously historymatched wells get worse throughout the iterative process. The two most likely causes are: (1) the limits definition of the influential region that affect each well and (2) the crossinfluence that may exist (a single region may affect wells in another region that is further away). As history matching is an inverse and nonlinear problem with multiple answers, to match a well, it is not enough to take into account only the regions and properties local to wells. As time progresses, the area of influence increases and adjacent wells start to impact well productions so it is possible that history matching one well may affect another.
Another influential factor is the defined region. Regions that take into account the influential area and reservoir behavior are difficult to define because streamlines are timedependent. Many cases end up using geometric regions like Voronoi, however, there is no guarantee that the region covers the entire influential area, or that the perturbed area is larger than necessary.
To deal with this problem, we use a different approach to perturb the reservoir image for categorical and continuous properties. A set of the best historymatched images is used to identify the pattern and characteristics that fit production data with the history data for each well. These similarities are searched within the defined region to avoid the perturbation of cell blocks that are far from the wells that we want to match. Although a region is defined, only some of the cell blocks within are perturbed and not the entire defined region.
The modifications made to the Probability Perturbation Method and CoSimulation are explained in Sections 2.1 and 2.2, followed by the application to UNISIMIH model.
2.1 Probability perturbation method
The Sequential Indicator Simulation (Journel and Alabert, 1990) used in this work for facies modeling uses the location of welllog data, the histogram of the property, and the modeled variogram to calculate a conditioned probability based on geological data P(AB). This probability represents the chance of facies A occurring, taking into account geological and petrophysical information B. It is then possible to create a set of models, which, because of the heterogeneity and different facies spatial distribution, provides different reservoir production answers that may, or may not, be a better fit to history data.
The Probability Perturbation Method integrates dynamic data (D) with probability P(AB) creating a posterior probability P(AB, D) to reproduce the models that best fit the history data. In the method developed by Hoffman (2005), a single image is chosen, based on the information of dynamic data (model match quality) to calculate the probability P(AD). The weight given to the selected image ranges between 0 and 1 and is chosen arbitrarily.
In our proposed methodology, we adapt the work of Oliveira et al. (2017a) to use a set of n models to calculate the probability P(AD). This identifies similarities within the chosen set. Considering this, P(AD) can be calculated using equations (1) and (2):(1) (2)where determines the occurrence of facies at location of model among . If the occurrence of facies A is high for the chosen reservoir images in a certain block, the probability will tend to have values near 1 to represent its high chance of occurrence.
To calculate the posterior probability , a conditional independence is assumed between B and D, and equation (3) is used as in Hoffman (2005):(3)
2.2 Image perturbation using cosimulation
CoSimulation was first used to simulate properties with few observed data but correlated with other properties that had more information (Deutsch, 2002; Yamamoto and Landim, 2013). Simulating porosity using seismic data, or permeability using porosity, are wellknown procedures in the industry. Lately, cosimulation has been used to reduce global or regional uncertainty in petrophysical properties, by using, as secondary information, a reservoir image with the spatial distribution that is to be reproduced (MataLima, 2008; MataLima and Soares, 2007).
Two inputs are necessary for cosimulation: a reservoir image of the property with the characteristics we want to preserve, and a correlation coefficient ranging from 0 to 1 that defines the influence of the chosen image in the perturbation.
The usual method to reduce uncertainty using CoSimulation chooses a single reservoir image as the reference. The correlation coefficient is usually given to a predefined geometric region, which is normally chosen despite reservoir dynamic behavior.
The approach proposed here consists in choosing a set of models with a good match and in using the similarities between them to generate possible new images that better fit the history data by perturbing only the influential cell blocks within a region instead of the entire region. This can be achieved using statistics (such as mean, variance, and coefficient of variation) to identify where, from all images, the variability is the smallest. Cells with high variability may then be defined as noncrucial to history matching for particular production data because there are no specific characteristics required to match the production data (e.g., low local permeability) and for this reason, the initial variability should be kept. By conditioning only cells with low variability, we minimize the number of conditioned cells, increasing variability in petrophysical characterization.
For each region, a set of n best models is chosen and equation (4) defines the secondary image (SI):(4)where is the property value at location and is the number of models used to build the secondary information. To evaluate variability, standard deviation () and Coefficient of Variation (CV) are used:(5) (6)
Cells with CV (Coefficient of Variation) below a given cutoff (Cut Value) are given a high correlation coefficient, while others receive a low value. The high value for the correlation coefficient can be a constant value defined by the user or function of the match quality as used in Oliveira (2014); while for the low correlation coefficient, a low value can be arbitrarily chosen (in this case, 0.05). We assume as a high value and as a low value for the correlation coefficient between the secondary information and the new reservoir image. The Sequential Gaussian CoSimulation requires each block to have a value for the secondary information and a correlation coefficient, which can be defined using equation (7):(7)
3 Methodology
The general methodology used in this work follows that fully detailed in Oliveira et al. (2017a) for reducing multiproperties and is briefly explained in this paper. Oliveira et al. (2017a) treat geostatistical images to achieve a set of models with good history matching; this process is integrated into a general workflow for history matching and uncertainty reduction (Fig. 1).

STEP 1: Define each uncertainty distribution and the parameterization limits ensuring that the combination of the uncertain attributes encompasses the history data.

STEP 2: Define an acceptable misfit limit for each objective function allowing the normalization of different production data using the Normalized Quadratic Difference with Sign (NQDS) function (detailed in Sect. 3.1).

STEP 3: Decide whether history matching or uncertainty reduction is preferable. Uncertainty reduction is designated for the objective of recharacterizing the distribution function while history matching has the objective of finding the best matched models, normally with an optimization algorithm.

STEP 4: Generate a set of models by combining different levels of uncertainty.

STEP 5: Simulate models.

STEP 6: Evaluate models variability using the NQDS function for each production data.

STEP 7: The iterative process concludes if the models exhibit acceptable behavior.

STEP 8: Review simulation model to evaluate if the problem is in the numeric model, for example, problems in productivity index, evaluate boundary conditions, aquifer, etc.

STEP 9: Perform the necessary changes highlighted in Step 8.

STEP 10: When, along the process, few or no models that match the data, a new parameterization is considered by identifying new influential attributes.

STEP 11: A reparameterization is done on the influential attributes.

STEP 12: Perform the uncertainty reduction on the attributes or the history matching using an optimization algorithm.

STEP 13: When a sufficient percentage of models has an acceptable match, generate new models to be filtered in the next step.

STEP 14: The filter will discard models that have an objective function outside the limits defined by a cut value.

STEP 15: Use the filtered models for forecasting or risk analysis.
Fig. 1. General workflow for history matching and uncertainty reduction (adapted from Oliveira et al., 2017a). 
The image perturbation method, used to improve reservoir models behavior is done in Step 12. In order to do that, the workflow presented in Figure 2, adapted from Oliveira et al. (2017a) is used:

STEP 12.1: Each loop should have an objective. At this point the objective function related to the mismatched production data should be chosen.

STEP 12.2: In order to match a specific objective function, it is necessary to understand which petrophysical properties most affect the dynamic data and in which region of the reservoir. If the influence between petrophysical properties and production is not well known, the sensitivity analysis described in more details in Oliveira et al. (2017a) and briefly explained in Section 3.2 should be done.

STEP 12.3: Perform the multivariate sensitivity analysis for petrophysical properties (detailed in Sect. 3.2).

STEP 12.4: Identify which attribute is more influential in the history matching process and follow the flowchart properly.

STEP 12.5: Perform the Probability Perturbation Method as in Section 2.1.

STEP 12.6 and 12.7: Perform Collocated CoSimulation as in Section 2.2.

STEP 12.8: Generate a new set of reservoir images.
Fig. 2. Flowchart for uncertainty reduction in geostatistical realizations (adapted from Oliveira et al., 2017a). 
The main contribution and changes applied from Oliveira et al. (2017a) in this work consists in considering a set of reservoir images to define the secondary information for probability perturbation method and cosimulation as explained in Sections 2.1 and 2.2.
3.1 Normalized quadratic difference and water breakthrough objective function
In order to quantify the quality of each scenario, the concept of Normalized Quadratic Difference with Sign (NQDS) is used as in Almeida et al. (2014) and Mesquita et al. (2015). We quantify the misfit between production and history data for bottomhole pressure, oil, water, and liquid rates for production wells, and bottomhole pressures and injection rates for injector wells. The NQDS can be expressed by equations (8) and (9):(8) (9)where, is the simulated data and is the history data for an observed data point for data type.
The main purpose is not to match the history data perfectly, but to establish an acceptable misfit between models to be considered matched. This is calculated using the definition of tolerance (Tol) and, for some data types (e.g. water rate), a constant C. Constant C is included to avoid the exclusion of models which may produce some water when there is no water on history data, resulting in high values of NQDS (even if water rate is small and the model has satisfactory match). This tolerance and constant are defined in STEP 2 of the general workflow for history matching and uncertainty reduction. Taking this into account, values of NQDS between −1 and 1 have an acceptable misfit and the objective function is considered matched.
To evaluate the variability and quality of the models the NQDS values are plotted (Fig. 3), showing whether the mean behavior of the objective function is above or below the history data.
Fig. 3. Examples of production curves and NQDS values (adapted from Oliveira et al, 2017b). 
The particular case of P3 from Figure 3 does not have any matched models despite having models above and below the acceptance limit; this usually happens for oil and water rates and is caused by incorrect water breakthrough times. For this reason, this particular case was treated using an objective function for Water Breakthrough Time (WBT):(10)where WBT_{Sim} is the time (days) that water takes to erupt in the production well in the simulation model, WBT_{Hist} is the time at which the production well starts to produce water in the history data, Tolerance is used to normalize the misfit and is the acceptable difference, in days, between the start of water production in simulation and history data.
3.2 Sensitivity analysis for petrophysical properties
Petrophysical characterization is crucial for reservoir production behavior; however, the order of magnitude of the uncertainty can be from thousand to millions due to the number of grid blocks and possible property values (e.g., porosity, permeability, etc) to occur in each cell. Understanding the reservoir and how perturbation should be done is essential in data assimilation.
In order to do that Oliveira et al. (2017a) proposed a method by chancing properties’ spatial distribution by using porosity and permeability multipliers on a chosen base case (for example, the model with best global match). The reservoir is split into different regions and the perturbation is done randomly in each region, creating a new set of models that will be used to correlate the property multiplier with the dynamic data response. These models are only used for the sensitivity analysis.
The creation of these models and the corresponding correlation between the petrophysical characterization and the dynamic data can be seen in Figures 4 and 5, respectively.
Fig. 4. Example of porosity models created for sensitivity analysis (from Oliveira et al., 2017a). 
Fig. 5. Plot between attribute and NQDS (from Oliveira et al., 2017a). 
4 Application and results
The proposed methodology was applied to a complex synthetic reservoir used for history matching and uncertainty reduction, which uses real well data from the Namorado Field, Campos Basin. The benchmark case study is known as UNISIMIH (Avansi and Schiozer, 2015). The model has a cornerpoint grid with 36.739 active cells, each measuring 100 × 100 × 8 m.
The production strategy has 14 producer and 11 injector wells. There is a total of 11 years of history data with information for oil, water, liquid and water injection rates, and bottomhole pressures. During the history matching, liquid rates were informed for the producers while injection rates were informed for the injector wells.
In this specific work, we focus only on the petrophysical properties and, for this reason, all other uncertain attributes have been reduced to a single level. Avansi and Schiozer (2015) extensively described all the uncertainties present in the benchmark case and a brief description of the geostatistical modeling is given as follows.
All the petrophysical properties are correlated somehow and they are facies, porosity and permeability. For facies a Sequential Indicator Simulation (SIS) was used, and nettogross is directly related to the facies type. Within each facies, the porosity was simulated independently using the Sequential Gaussian Simulation. To simulate permeability, instead of using a deterministic approach as in Avansi and Schiozer (2015), a cross plot is used for a probabilistic approach while guaranteeing that the relationship between properties is respected as in (Oliveira, et al., 2017a).
An initial set of 150 reservoir images were created to study the effects of petrophysical properties on production data. Variability in production data for each model was evaluated using the production curves and plotting the NQDS values for each objective function.
Table 1 lists the tolerance and constant values used to calculate the NQDS and WBT. The same tolerance and constant values were used for all wells. However, these can vary depending on how reliable the history data is or how rigorous we want the history match to be.
Tolerances and constants for each objective function.
As seen in Figures 6–10, the high variability in the initial set of models is due to the uncertainty present in petrophysical properties. Because the injection rates and liquid production are not honored, field average pressure has a large mismatch. For this reason, the first iteration seeks to improve the injection rates and liquid production for each well. We performed a preliminary test in the first iteration to compare the proposed methodology, which uses a set of reservoir images to perform the image perturbation, with the more common approach of using a single image.
Fig. 6. Average field pressure for the initial set of models. 
Fig. 7. Cumulative field liquid production for the initial set of models. 
Fig. 8. Cumulative field water injection for the initial set of models. 
Fig. 9. Cumulative field oil production for the initial set of models. 
Fig. 10. Cumulative field water production for the initial set of models. 
4.1 Preliminary test
The common method for image perturbation is compared with the proposed methodology. The common method referred to as “Approach A” uses a single image for the image perturbation, while “Approach B” refers to the proposed methodology and uses a set of images for each region to locally perturb the reservoir image. The defined region used to look for a pattern between images was the Voronoi, defining a region for each well, due to simplicity in using and to use the same regions as in Oliveira et al. (2017a). In order to choose the best models, to be used in the perturbation, all objective functions were taken into account for each well inside the Voronoi region. In order to be with a good match, every objective function from the well should be within an acceptable range of NQDS value. To improve the injection rate for the wells, we reduced the uncertainty for facies and porosity, leaving only permeability to perturb together with the production wells, to avoid creating preferential flux paths without excluding the producer match. Starting with the same initial set of images, each case is described in Table 2.
Test between the common and proposed methodology for image perturbation.
The initial set of reservoir images is the same for both approaches. The sequential uncertainty reduction in reservoir images was performed, and in each loop, a new set of 150 reservoir images was created.
A major impact of the proposed methodology is that it seeks to identify a similar spatial pattern and perturb only those cells that may influence the history matching; for this reason, fewer cells are conditioned within a certain region in approach B than in approach A (Fig. 11).
Fig. 11. Correlation coefficient for the first iteration for approaches A and B. 
The number of images used in the process was previously stabilised taking into account the number of models with reasonable match that could be used to build the reference image. The number of images chosen was 5, less than that, the approach would start to be no different than using a single image.
Figure 12 shows that when calculating the percentage of occurrence of facies type 1 (from the set of 150 models) in the cell after the perturbation, the spatial pattern of facies distribution is highly predictable when using a single image for perturbation; a geological structure is easily identified (left side). However, the same is not true when using a set of reservoir images to perform the Probability Perturbation Method, the probability changes more smoothly. The prior probability changes but the probability goes smoothly through the reservoir, and it is difficult to predict the limit of occurrence for each facies (right side).
Fig. 12. Percentage of occurrence of a facies in the set of 150 reservoir images. 
The coefficient of variation was calculated for porosity. Figure 13 shows that, for image perturbation, variability is higher when using a set of reservoir images than for only a single image. Using a single image in cosimulation and a correlation coefficient for the entire region excessively reduces uncertainty, affecting not only the conditioned region but also the surrounding area, because of continuity.
Fig. 13. Coefficient of variation for porosity for approach A and B for the generated set. 
Most importantly, not only the variability in the distribution of spatial properties, but also the effects on production data must be considered. To assess these effects, the NQDS for all objective functions (Fig. 14) are compared. Figure 14 shows that both methods improved the injection rate.
Fig. 14. NQDS functions (water rate and BHP) for injector wells (initial set, perturbation using 1, or 5 images). 
For INJ022, for example, the water rate significantly improved using five images compared to a single image. However, the values of the correlation coefficient at the top of reservoir suggest little perturbation in the region (highlighted region in Fig. 15). On the other hand, analysis beneath the top of the reservoir, for the proposed method, identified geological structures that are crucial to the history matching. In Figure 15B, a channel is highlighted beneath the top of the reservoir with a high value of correlation coefficient, while the remaining cells have a low value.
Fig. 15. Correlation coefficient with highlighted channel beneath the top of the reservoir. 
4.2 Complete history matching
The history matching of the reservoir followed the proposed methodology. The following step is to history match the producers.
First, using the sensitivity analysis described in Oliveira et al. (2017a), the regions that most affected each producer were identified and the producerinjector groups were defined.
Second, the best images from each loop were chosen, (lowest values of NQDS) for each producerinjector group. These models were used as secondary information for the Probability Perturbation Method and CoSimulation inputs.
Each iteration had an objective requiring several loops to reach a desirable match. Each loop generated a set of 150 reservoir images. Table 3 summarizes how many loops were required for each objective.
Iterations and required loops for each objective to history match the reservoir model.
As stated before, the first iteration matched the injectivity and productivity of the wells. This not only improved average reservoir pressure but also the bottomhole pressure of producers and injectors.
The following four iterations matched all the producers of the model. Because of the complexity of the case study and the crossinfluence between the dynamic data and each property within the region, adjustments were necessary to improve the history matching for some objective functions. These adjustments are local recharacterizations, for example, the history matching for PROD09 was very complex. In this specific case, a downscaling further divided the region and then a sensitivity analysis was performed. The downscaling was necessary to better understand how the heterogeneity within the region affected well production. In fact, we were able to identify that a barrier was necessary to decrease water production while the remaining area should be kept highly permeable to ensure well productivity.
A set of 26 models was finally selected using a cut value of 15 for the Normalized Quadratic Function. The history matches for the field and wells were successfully achieved as shown in Figures 16–18.
Fig. 16. Field cumulative liquid production, for the first iteration and the filtered models. 
Fig. 17. Field cumulative oil production, for the first iteration and the filtered models. 
Fig. 18. Field cumulative water injection, for the first iteration and the filtered models. 
Most importantly, not only were all models created through the process of history matching geologically consistent and realistic, but they also reliably forecasted production, the ultimate purpose of model calibration (Figs. 19–21).
Fig. 19. Field cumulative liquid production forecast. 
Fig. 20. Field cumulative oil production forecast. 
Fig. 21. Field cumulative water injection forecast. 
To highlight the advantages and contribution of our proposed methodology, we compared the history matching resulting from this work (Approach B) and that using the same model but with a single image to define the secondary information (Oliveira et al., 2017a), approach A. The comparison showed that our methodology presented improvements for both history matching and production forecasting, also the resulting reservoir images have more variability using the proposed methodology than using a single image for image perturbation (Fig. 22).
Fig. 22. Coefficient of variation for permeability using Approach A and Approach B. 
Fig. 23. NQDS for pressure and liquid rates in the first and final iterations for approaches A and B. 
Fig. 24. NQDS for oil and water rates in the first and final iterations for approaches A and B. 
Using an indicator to quantify the misfit for water breakthrough time efficiently improved the production data match. The water breakthrough in wells was corrected (Fig. 25) and symmetry in NQDS values for oil and water rates prevailed (Fig. 24). Comparing the NQDS values for all objective functions, we see that approach A excessively reduced the production variability (Figs. 23–25). On the other hand, because approach B tries to keep variability in petrophysical characterization, it affects the variability in production data. The best scenario would be to have symmetrical NQDS value above and below a perfect match.
Fig. 25. NQDS water breakthrough time in the first and final iterations for approaches A and B. 
Figure 25 shows that the water breakthrough time should be considered when history matching. Because of the large amount of history data, the water production may appear to match while the water breakthrough does not, significantly affecting production forecast. Increasing the number of objective functions will result in minimizing the risk of making an incorrect production forecast.
The mismatch in water breakthrough significantly affects the match of field production and consequently the production forecast as represented in Figure 26. Approach A converged to a local minimum and reduced uncertainty excessively, whereas Approach B generated a reliable production forecast.
Fig. 26. Cumulative oil production forecasts for the initial set of models, the matched models using the classic and the proposed approach. 
5 Conclusion
The proposed methodology to reduce uncertainty in reservoir characterization efficiently achieved good history matching and a reliable production forecast for a challenging benchmark case (UNISIMIH). Using a set of models, instead of a single matched model, to identify the characteristics and geological events to match the dynamic data, we improved the quality in history matching while increasing variability in the spatial distribution of petrophysical properties. This contribution to probabilistic history matching avoids convergence to a local minimum, while searching for different answers.
Using a single image for image perturbation has demonstrated that matching the petrophysical characteristics (facies, porosity, and permeability) to the history data may easily converge to a set of similar models with little variability. This can potentially excessively reduce uncertainty in the first steps of history matching, resulting in future difficulty when exploring other possible solutions.
Another disadvantage of using a single image for the probability perturbation method is that from the posterior probability function it is possible to easily highlight the limits of the geological events such as channels, proving that the created models are very similar (unique solution). However, this does not happen when using a set of images to calculate the posterior probability.
The usual method to perform the cosimulation assigns a correlation coefficient to the entire predefined region. Parameterizing a predefined region, usually geometric, is a simple and efficient approach. However, it excessively reduces uncertainty either through locally perturbing a region that does not affect the history match of the dynamic data, or affecting nearby wells. On the other hand, regions that factor in reservoir behavior, such as streamlines, may be difficult to parameterize because the region of influence is timedependent. Our proposed approach has the advantage of identifying the local characteristics that match the data, without conditioning an entire region.
The inclusion of the misfit for water breakthrough time has been extremely useful along with the other objective functions and we highly recommend its use.
For future works, other uncertain attributes could be included in the historymatching process to better meet requirements in the industry.
Acknowledgments
The authors are grateful to PETROBRAS, the research network SIGER 3 (Grant Agreement No. 0050.0100204.16.9), the National Agency of Petroleum, Natural Gas and Biofuels (ANP), the Center of Petroleum Studies (CEPETROUNICAMP) particularly the UNISIM research group, the Department of Energy of the School of Mechanical Engineering of the University of Campinas (DEFEMUNICAMP), the Foundation CMG and National Council for Scientific and Technological Development (CNPq), for supporting this research, and also, to the Computer Modelling Group (CMG) and Schlumberger for software licenses.
References
 Almeida F., Davolio A., Schiozer D. (2014) A New Approach to Perform a Probabilistic and MultiObjective History Matching, SPE170623MS, SPE Annual Technical Conference and Exhibition, 27–29 October. [Google Scholar]
 Avansi D., Schiozer D. (2015) UNISIMI: synthetic model for reservoir development and management applications, Int. J. Model. Simul. Pet. Ind. 9, 21–30. [Google Scholar]
 Caeiro M., Demyanov V., Soares A. (2014) Multiobjective History Matching of a Deltaic Reservoir with NonStationary Geostatistical Modelling, European Conference on the Mathematics of Oil Reservoir. [Google Scholar]
 Christie M., Eydinov D., Demyanov V., Talbot J., Arnold D., Shelkov V. (2013) MultiObjective Algorithms in History Matching of a Real Field, Matching of a Real Field, SPE 163580MS. [Google Scholar]
 Deutsch C. (2002) Geostatistical Reservoir Modeling, Applied Geostatistics Series, Oxford University Press, p. 376. [Google Scholar]
 Ding D.Y., Roggero F. (2010) History matching geostatistical model realizations using a geometrical domain based parameterization technique, Math. Geosci. 42, 4, 413–432. [Google Scholar]
 GervaisCouplet V., Le Ravalec M. (2018) Identifying influence areas with connectivity analysis – application to the local perturbation of heterogeneity distribution for history matching, Comput Geosci 22, 3–28, doi: https://doi.org/10.1007/s105960179663y. [Google Scholar]
 GervaisCouplet V., Gautier Y., RavalecDupin L., Roggero F. (2007) History Matching Using Local Gradual Deformation, SPE107173, EUROPEC/EAGE Conference and Exhibition. [Google Scholar]
 Hajizadeh Y., Christie M., Demyanov V. (2011) Towards multiobjective history matching: Faster convergence and uncertainty quantification. SPE 141111. [Google Scholar]
 Hoffman B. (2005) Geologically consistent history matching while perturbing facies, PhD Dissertation, Stanford University. [Google Scholar]
 Hoffman B., Caers J. (2003) Geostatistical history matching using a regional probability perturbation method, SPE 84409MS, SPE Annual Technical Conference and Exhibition, 5–8 October. [Google Scholar]
 Hoffman B., Caers J. (2005) Regional Probability perturbation for history matching, J. Petrol. Sci. Eng. 46 (1–2), 53–71. [CrossRef] [Google Scholar]
 Hu L.Y. (2000) Gradual deformation and iterative calibration of Gaussianrelated stochastic models, Math. Geol. 32, 1, 87–108. [Google Scholar]
 Hutahaean J., Demyanov V., Christie M. (2015) Impact of model parameterization and objective choices on assisted history matching and reservoir forecasting, SPE 176389MS, SPE/IATMI Asia Pacific Oil & Gas Conference and Exhibition, 20–22 October. [Google Scholar]
 Journel A.G., Alabert F.G. (1990) New method for reservoir mapping, J. Petrol. Technol. 42, 212–218. [CrossRef] [Google Scholar]
 MataLima H., Soares A. (2007) Geostatistical history matching with direct sequential transformation of images. Petrol. Geostat.. http://earthdoc.eage.org/publication/publicationdetails/?publication=7865 [Google Scholar]
 MataLima H. (2008) Reservoir characterization with iterative direct sequential cosimulation: integrating fluid dynamic data into stochastic model, J. Petrol. Sci. Eng., Elsevier Science 63, 3–4, 59–72. [CrossRef] [Google Scholar]
 Mesquita F., Davolio A., Schiozer D. (2015) A systematic approach to uncertainty reduction with a probabilistic and multiobjective history matching, SPE174359MS, EUROPEC 2015, 1–4 June. [Google Scholar]
 Oliveira G.S. (2014) Integração da Modelagem Geostatística com o Processo de Ajuste de Histórico. Master Degree Thesis, Instituto Superior Técnico, Lisboa. [Google Scholar]
 Oliveira G., Schiozer D.J., Maschio C. (2017a) History matching by integrating regional multiproperty image perturbation methods with a multivariate sensitivity analysis, J. Petrol. Sci. Eng. 153, 111–122, doi: http://dx.doi.org/10.1016/j.petrol.2017.03.031. [CrossRef] [Google Scholar]
 Oliveira G., Soares A., Schiozer D., Maschio C. (2017b) Reducing uncertainty in reservoir parameters combining history matching and conditioned geostatistical realization, J. Petrol. Sci. Eng. 156, 75–90, doi: https://dx.doi.org/10.1016/j.petrol.2017.05.003. [CrossRef] [Google Scholar]
 Roggero F., Hu L.Y. (1998) Gradual deformation of continuous geostatistical models for history matching, Paper SPE 49004, SPE Annual Technical Conference and Exhibition. [Google Scholar]
 SchulzeRiegert R., Kroschem M. (2007) Multiobjective optimization with application to model validation and uncertainty quantification, SPE 105313, SPE Middle East Oil and Gas Show and Conference, 11–14 March. [Google Scholar]
 Yamamoto J., Landim P. (2013) Geostatística: conceitos e aplicações, Oficina de Textos, p. 215. ISBN 9788579750779. [Google Scholar]
All Tables
Iterations and required loops for each objective to history match the reservoir model.
All Figures
Fig. 1. General workflow for history matching and uncertainty reduction (adapted from Oliveira et al., 2017a). 

In the text 
Fig. 2. Flowchart for uncertainty reduction in geostatistical realizations (adapted from Oliveira et al., 2017a). 

In the text 
Fig. 3. Examples of production curves and NQDS values (adapted from Oliveira et al, 2017b). 

In the text 
Fig. 4. Example of porosity models created for sensitivity analysis (from Oliveira et al., 2017a). 

In the text 
Fig. 5. Plot between attribute and NQDS (from Oliveira et al., 2017a). 

In the text 
Fig. 6. Average field pressure for the initial set of models. 

In the text 
Fig. 7. Cumulative field liquid production for the initial set of models. 

In the text 
Fig. 8. Cumulative field water injection for the initial set of models. 

In the text 
Fig. 9. Cumulative field oil production for the initial set of models. 

In the text 
Fig. 10. Cumulative field water production for the initial set of models. 

In the text 
Fig. 11. Correlation coefficient for the first iteration for approaches A and B. 

In the text 
Fig. 12. Percentage of occurrence of a facies in the set of 150 reservoir images. 

In the text 
Fig. 13. Coefficient of variation for porosity for approach A and B for the generated set. 

In the text 
Fig. 14. NQDS functions (water rate and BHP) for injector wells (initial set, perturbation using 1, or 5 images). 

In the text 
Fig. 15. Correlation coefficient with highlighted channel beneath the top of the reservoir. 

In the text 
Fig. 16. Field cumulative liquid production, for the first iteration and the filtered models. 

In the text 
Fig. 17. Field cumulative oil production, for the first iteration and the filtered models. 

In the text 
Fig. 18. Field cumulative water injection, for the first iteration and the filtered models. 

In the text 
Fig. 19. Field cumulative liquid production forecast. 

In the text 
Fig. 20. Field cumulative oil production forecast. 

In the text 
Fig. 21. Field cumulative water injection forecast. 

In the text 
Fig. 22. Coefficient of variation for permeability using Approach A and Approach B. 

In the text 
Fig. 23. NQDS for pressure and liquid rates in the first and final iterations for approaches A and B. 

In the text 
Fig. 24. NQDS for oil and water rates in the first and final iterations for approaches A and B. 

In the text 
Fig. 25. NQDS water breakthrough time in the first and final iterations for approaches A and B. 

In the text 
Fig. 26. Cumulative oil production forecasts for the initial set of models, the matched models using the classic and the proposed approach. 

In the text 