A Step Forward to Closing the Loop between Static and Dynamic Reservoir Modeling

A Step Forward to Closing the Loop between Static and Dynamic Reservoir Modeling — The current trend for history matching is to find multiple calibrated models instead of a single set of model parameters that match the historical data. The advantage of several current workflows involving assisted history matching techniques, particularly those based on heuristic optimizers or direct search, is that they lead to a number of calibrated models that partially address the problem of the non-uniqueness of the solutions. The importance of achieving multiple solutions is that calibrated models can be used for a true quantification of the uncertainty affecting the production forecasts, which represent the basis for technical and economic risk analysis. In this paper, the importance of incorporating the geological uncertainties in a reservoir study is demonstrated. A workflow, which includes the analysis of the uncertainty associated with the facies distribution for a fluvial depositional environment in the calibration of the numerical dynamic models and, consequently, in the production forecast, is presented. The first step in the workflow was to generate a set of facies realizations starting from different conceptual models. After facies modeling, the petrophysical properties were assigned to the simulation domains. Then, each facies realization was calibrated separately by varying permeability and porosity fields. Data assimilation techniques were used to calibrate the models in a reasonable span of time. Results showed that even the adoption of a conceptual model for facies distribution clearly representative of the reservoir internal geometry might not guarantee reliable results in terms of production forecast. Furthermore, results also showed that realizations which seem fully acceptable after calibration were not representative of the true reservoir internal configuration and provided wrong production forecasts; conversely, realizations which did not show a good fit of the production data could reliably predict the reservoir behavior. Thus a statistical approach was confirmed to be the only way to reduce the uncertainty inherent to reservoir modeling and should be adopted as a standard in reservoir studies. LIST OF SYMBOLS u Porosity k Effective permeability (m) kabs Absolute permeability (m ) kr Relative permeability yk True model state vector Cpp Error covariance matrix AGM Adaptive Gaussian Mixture EnKF Ensemble Kalman Filter FOPT Field cumulative oil production FWPT Field cumulative water production OF Objective Function WBHP Well Bottom Hole Pressure (Pa) WOPR Well Oil Production Rate (m/day) WGOR Well Gas Oil Ratio WWCT Well water cut h Bandwidth parameter INTRODUCTION The main goal of a reservoir study is the integration of a large number of data, expressing all the available geological, geophysical, petrophysical and engineering information, into a model that, ideally, should describe with proper accuracy the structural and petrophysical properties of the reservoir (Cosentino, 2001; Benetatos and Viberti, 2010). In general, the ability of reproducing the historical dynamic behavior of a field is taken as an indicator that the model is a reliable representation of the real system so that it can be used for prediction purposes. Therefore, considerable effort is devoted to the calibration process, also known as history matching, in the belief that a calibrated model will provide reliable estimates of hydrocarbon recovery. The traditional manual history matching workflow consists in a sequential approach that begins with the matching of global or field parameters, e.g. field pressure, and follows with the adjustment of individual flow units or layers, to finally match the well and near-wellbore data (Mattax and Dalton, 1990). This approach has proved to be very flexible because engineers can tune the values of the reservoir parameters based on their own experience as well as good judgment; however, its disadvantage is being an extremely slow trial and error procedure. Furthermore, a good fit for the production data does not necessarily give a good estimation of the parameters of the reservoir and this might lead to errors in the prediction of the reservoir performance (Tavassoli et al., 2004). Additionally, an assessment of the overall uncertainty related to the calculated reserves and the associated production profiles is often required. When a single reservoir model is available a number of simulations 1202 Oil & Gas Science and Technology – Rev. IFP Energies nouvelles, Vol. 69 (2014), No. 7

Abstract -A Step Forward to Closing the Loop between Static and Dynamic Reservoir Modeling -The current trend for history matching is to find multiple calibrated models instead of a single set of model parameters that match the historical data.The advantage of several current workflows involving assisted history matching techniques, particularly those based on heuristic optimizers or direct search, is that they lead to a number of calibrated models that partially address the problem of the non-uniqueness of the solutions.The importance of achieving multiple solutions is that calibrated models can be used for a true quantification of the uncertainty affecting the production forecasts, which represent the basis for technical and economic risk analysis.In this paper, the importance of incorporating the geological uncertainties in a reservoir study is demonstrated.A workflow, which includes the analysis of the uncertainty associated with the facies distribution for a fluvial depositional environment in the calibration of the numerical dynamic models and, consequently, in the production forecast, is presented.The first step in the workflow was to generate a set of facies realizations starting from different conceptual models.After facies modeling, the petrophysical properties were assigned to the simulation domains.Then, each facies realization was calibrated separately by varying permeability and porosity fields.Data assimilation techniques were used to calibrate the models in a reasonable span of time.Results showed that even the adoption of a conceptual model for facies distribution clearly representative of the reservoir internal geometry might not guarantee reliable results in terms of production forecast.Furthermore, results also showed that realizations which seem fully acceptable after calibration were not representative of the true reservoir internal configuration and provided wrong production forecasts; conversely, realizations which did not show a good fit of the production data could reliably predict the reservoir behavior.Thus a statistical approach was confirmed to be the only way to reduce the uncertainty inherent to reservoir modeling and should be adopted as a standard in reservoir studies.

Porosity k
Effective permeability (m 2 ) k abs Absolute permeability (m 2 ) k r Relative permeability y t

INTRODUCTION
The main goal of a reservoir study is the integration of a large number of data, expressing all the available geological, geophysical, petrophysical and engineering information, into a model that, ideally, should describe with proper accuracy the structural and petrophysical properties of the reservoir (Cosentino, 2001;Benetatos and Viberti, 2010).In general, the ability of reproducing the historical dynamic behavior of a field is taken as an indicator that the model is a reliable representation of the real system so that it can be used for prediction purposes.Therefore, considerable effort is devoted to the calibration process, also known as history matching, in the belief that a calibrated model will provide reliable estimates of hydrocarbon recovery.The traditional manual history matching workflow consists in a sequential approach that begins with the matching of global or field parameters, e.g.field pressure, and follows with the adjustment of individual flow units or layers, to finally match the well and near-wellbore data (Mattax and Dalton, 1990).This approach has proved to be very flexible because engineers can tune the values of the reservoir parameters based on their own experience as well as good judgment; however, its disadvantage is being an extremely slow trial and error procedure.Furthermore, a good fit for the production data does not necessarily give a good estimation of the parameters of the reservoir and this might lead to errors in the prediction of the reservoir performance (Tavassoli et al., 2004).Additionally, an assessment of the overall uncertainty related to the calculated reserves and the associated production profiles is often required.When a single reservoir model is available a number of simulations are performed varying the input parameters and evaluating discrepancies from a reference scenario.At this point, a great number of production profiles are simulated and the overall uncertainty can be evaluated.This task may diversify depending on the development stage of the field.For undeveloped fields, without a production history to match, the aim is to define a significant set of production profiles.Sensitivity studies, known as "risk analysis", are often carried out based on different reservoir realizations, i.e. distribution of static properties, to assess the uncertainty affecting the estimate of the hydrocarbon originally in place.However, when it comes to production forecasts, a reference case is selected among the reservoir realizations.Then, only one or few parameters of the dynamic model are varied at a time.The best approach is to concentrate on the parameters that are deemed to have a significant impact on the reservoir behavior; however, dependencies among the considered parameters exist and should be taken into account.Even if a production history exists, the evaluation of the uncertainty affecting the production forecast is a very difficult task.Having fixed the static and dynamic data of the reference scenario during the history matching phase, the only parameters that can be varied in the forecast phase are the well count, location and completion and the facilities constraints.Nevertheless, because the history match is an ill-posed problem and the solution is not unique, the representativeness of uncertainty in the forecast phase is definitely not guaranteed.Several optimal solutions, or reservoir history matches, are needed to truly evaluate the uncertainty associated to a calculated production forecast and thus to the recoverable reserves (Selberg et al., 2006).The possibility to obtain multiple calibrated models partially addresses the problem of the nonuniqueness of the solutions and at the same time the set of suboptimal models can provide a more reliable prediction of the reservoir performance, since the overall uncertainty is reduced by integrating geological and production data on a wide range of scenarios.To this end assisted history matching techniques can prove of great help (Fig. 1).

CLOSED LOOP SIMULATION
Historically, during a reservoir study the work performed by geoscientists in each discipline was not integrated.The results were handed over from discipline to Uncertainty assessment using multiple calibrated models.
discipline without any active interconnection and with limited information exchange.Each discipline involved in the creation of the reservoir model had to provide data and highly accurate results in order to minimize the uncertainties during the construction of the reservoir model.This process was based on the claim that the independent accuracy of each discipline leads to an overall accuracy of the final model.However, this approach showed strong limitations in the case inconsistencies arose during the data processing, requiring a thorough and consistent re-evaluation of all the model parameters.
In particular, the history matching workflow was limited to the calibration of the dynamic model by changing few parameters of a given static model according to a deterministic approach, so the workflow was unidirectional.
As an example, local variations of permeability might be introduced in the model by the engineers regardless of the facies and the petrophysical properties distribution initially defined in the geological model, or variations of porosity could be adopted in order to pseudo-adjust a structural problem.However, in recent years, this trend has changed and the construction of the static model is more often conditioned by the dynamic data, closing the loop that had remained open for many years (Fig. 2).This approach guaranties a closer and faster collaboration between the different disciplines involved in a reservoir study because the feedback coming from the history match can be potentially incorporated into the structural, static and dynamic model.In order to account for geological uncertainties and integrate dynamic data to geological model several algorithms of different nature have been proposed, such as: spectral decomposition (Reynolds et al., 1996); pilot point method (RamaRao et al., 1995) gradual deformation method (Hu, 2000), probability perturbations methods (Caers, 2003) and spectral co-simulation perturbation (Le Ravalec-Dupin and Da Veiga, 2001).History matching is often performed on the basis of gradient-based methods or evolutionary strategies (Schulze-Riegert, et al., 2001;Schulze-Riegert and Haase, 2003;Zitzler, 1999).

APPLICATION OF THE ENKF/AGM TO ADDRESS FACIES DISTRIBUTION UNCERTAINTY
A workflow, which includes the analysis of the uncertainty associated with the facies distribution of a fluvial depositional environment, is presented below.This type of depositional environment was selected because it is particularly representative of complex internal geometries which are found in several reservoirs worldwide.The Ensemble Kalman Filter (EnKF) modified according to Adaptive Gaussian Mixture (AGM) filter in order to loosen up the EnKF requirement of a Gaussian prior distribution of the petrophysical properties was used for data assimilation.

Model Description
A synthetic reservoir model was set up for the purpose of this research.completes the reservoir seal.No active aquifer is present.Undersaturated oil was considered, so no gas cap exists at the initial state.A total of 11 producing wells were drilled and most of them were located at the top of the structure.Water was injected through 5 wells for pressure support.The reservoir domain and well locations are shown in Figure 3. Two different facies, channel sands and floodplain, typical of a fluvial depositional environment were modeled.The facies were distributed using an object modeling technique.The channel parameters were set as follows: -a net-to-gross of 40% was set in all the domain; -an orientation between 20 and 40 compass degrees; -a triangular probability distribution was used to define the amplitude, wavelength, width and thickness of the channels.The values used for each parameter distribution are summarized in Table 1.The facies modeling process was constrained to the well facies data.The facies distribution at each well is presented in Figure 4.The facies distribution in the 3D model varies in both horizontal and vertical directions, however, in the following, for simplicity's sake, only the first layer is shown, the remaining layers showing similar behavior.
The petrophysical properties were distributed in the model using a sequential Gaussian simulation algorithm for each facies type.A spherical variogram and a normal probability distribution were used for porosity and log-permeability for both facies and for all layers simultaneously.The first step was to generate a porosity field for the channel and floodplain facies.The correlation length for both the major and minor directions of the variogram was set to 2 000 m for the flood plain facies and to 5 000 m for the channel sands.The mean values of porosity were set to 3% and 16% with 2% and 6% standard deviation for floodplain and channel sands, respectively (Fig. 5).Furthermore, in order to squeeze the tails to acceptable physical values, an exponential transformation was applied to both the fields.Water saturation and pressure field at the initial state.Facies data at the wells.
The log-permeability (XY direction) field was generated using the sequential Gaussian co-simulation based upon cokriging algorithm (Deutsch, 2002) with a correlation factor to the porosity field of 0.8.The mean values of permeability were set to 7 mD and to 250 mD with a standard deviation of 0.05 and 0.4 in the logarithmic "space" for the floodplain and channel sands, respectively.The vertical permeability distribution was obtained by multiplying the horizontal permeability by the anisotropy factor, set equal to 0.1.
The facies modeling workflow and the procedure to assign the petrophysical properties are illustrated in Figure 6.
The petrophysical properties of the top layer as resulting from the facies distribution and characterization are shown in Figure 7.
The modeling workflow described above is intrinsically characterized by a high uncertainty level since the number of control points, i.e. well data, is too low in respect to the degree of freedom of the whole model.Thus the proposed History Matching (HM) workflow Porosity and permeability fields for layer 1 (true case).
will involve not only the uncertainty related to the model petrophysical characterization, formalized as porosity and permeability distributions, but it will also involve different conceptual facies models and realizations in order to better integrate the facies/petrophysical modeling with HM through consolidated tools.

EnKF and AGM
A sequential data assimilation technique is a process which aims at estimating and predicting (analysis step) an unknown true state y t k of the system by integrating the forward model (system dynamics) F in time (k) using measurements, whenever available, to reinitialize the model before the integration continues.The model equation can be written as: where q is the unknown model error over one time step.
The standard Kalman Filter is an optimal sequential data assimilation method for linear dynamics (i.e.F is a linear operator) under the hypothesis that the prior distributions of the state vector components and of the measurement errors are all Gaussian and unbiased, i.e. white noise is assumed (Evensen, 2009).Unfortunately, the majority of real world problems are highly nonlinear (and petroleum engineering problems belong to such category).In order to overcome the standard KF limitation, subsequent variants of the original scheme were proposed, such as the Ensemble Kalman Filter (EnKF) introduced by Evensen in 1994 (Evensen, 1994), which uses an ensemble representation of the solution states.Sets of ensemble realizations are generated using Monte-Carlo sampling for the initial state, model noise and measurement noise.Ensemble members are then forwarded in time by solving the nonlinear state equations and are analyzed by an approximate Kalman filter scheme.The ensemble covariance is used as an approximation of the true covariance, thus avoiding explicit evolution and storage of the covariance as well.The proposed method can therefore be used with realistic highly nonlinear models on large domains and it is also well suited for distributed computing environments.When applied to the history matching problem the ensemble Kalman filter is able to estimate a large number of model variables and to assimilate a wide range of data types and can be easily integrated with existing reservoir simulators.
In order to loosen up the EnKF requirement of a Gaussian prior distribution the Adaptive Gaussian Mixture filter (AGM) was introduced by Stordal et al. in 2011(Stordal et al., 2011, 2012;Valenstrand et al., 2012).Gaussian mixture filters are based on the assumption that at each time-step k the p density f y k d 1:kÀ1 À Á can be approximated by a Gaussian kernel density estimator: where are scalar weights which satisfy P N e i¼1 x i kÀ1 ¼ 1.In the construction of the N e Gaussian kernels a bandwidth parameter h 2 ð0; 1 is introduced.The choice of the bandwidth parameter h determines the magnitude of the Kalman filter update step.This is treated as a design parameter which can be beneficial to reduce the risk of filter divergence.
In the proposed scheme, Stordal assigns the same covariance matrix to each kernel.In this way, each ensemble member is updated linearly, proportionally to the EnKF update.This filter suffers from a degeneracy phenomenon and the most efficient way to get around it is resampling, i.e. drawing new particles according to the distribution of the ensemble and then reassigning them the uniform weights.In order to reduce weight collapse and consequent resampling, which is time consuming, an adaptive parameter a can be introduced, extending the scheme introduced in by Hoteit et al. (2008).
In this framework, the EnKf can be viewed as a mixture filter, obtained by using Gaussian kernels and forcing the weights to be uniform.

Methodology
The main goal of the case study discussed below was to include the uncertainty associated to the facies distribution, both in terms of conceptual models parameters and petrophysical characterization, in the calibration of the numerical model and, consequently, in the production forecast.The proposed workflow is presented in Figure 8.
The first step in the workflow was to generate a set of facies realizations starting from different conceptual models.Here a conceptual model is intended as a given set of geometrical parameters that are used to stochastically generate different facies realizations (Fig. 9).For example, straight channels oriented in the north-south direction are conceptually different from meandering channels oriented in the east-west direction.
After facies modeling the petrophysical properties were assigned to the simulation domain.The number of petrophysical realizations needed for each facies realization depends on the dimensions of the model and on the number of historical parameters to be matched.In the analyzed cases, porosity and logpermeability values have to be assessed for each cell of the domain.In general, no more than 200 ensemble members are required for a medium complexity problem in which the number of parameters -to be assimilated at each assimilation step -does not exceed 100.In this case, 60 ensemble members were needed.Then, each facies realization was calibrated separately with the aid of the AGM algorithm.The optimal value of the bandwidth parameter h is known to be problem dependent and the optimization requires a trial-and-error procedure.In this work, the bandwidth parameter h was set to 0.5, from a preprocessing analysis.After completion of all the assimilation steps, simulations were run again the final porosity and permeability fields.The quality   of the match for each facies realization was computed using the following misfit function:

Conceptually similar Conceptually different
where the component OF w obs represents the contribution to the total objective function of the observation obs obtained at the well w: the term d k À Á w obs is the ensemble mean of the observation obs obtained at the well w at the time k, i.e.: and the coefficient r 2 obs;w is the variance of each observation at each well during the whole history matching: In this way, the different scale sizes of the observations were nearly normalized.
Eventually, the posterior distribution of the field cumulative oil and water production were compared against the prior distributions.

Reservoir State Vector
The state vector contains the variables of interest.It describes the state of the dynamic system and represents its degrees of freedom.The state vector y k j mainly consists of three data types: static reservoir properties, dynamic reservoir properties and measured values: where j refers to the jth ensemble member and k is the time step index for time t k .
In this case study, the static properties were represented by porosity and logarithmic transformation of the horizontal and vertical permeabilities for each gridblock, that is: The dynamic variables v k j were the pressure and water saturation values for each gridblock: The number of elements of the state vector corresponding to the observed data d k j can vary depending on the measurement availability at a given assimilation time step.Bottom hole pressures were observed at both producing and injector wells, while block pressure, water-cut and oil production rates were recorded at the production wells only: The maximum length of the state vector is: A total of 100 petrophysical realizations were created for each facies realization.The variogram parameters and probability density functions were equal to those used for the generation of the true case.The variogram parameters and the values of the probability density functions are summarized in Table 5 and Table 6, respectively.
Examples of petrophysical realizations from different facies models are shown in Figure 11.

Assimilation of Production Data
Two different scenarios were considered; in each scenario different geostatistical realizations (although generated with the same input parameters) and different production histories were used.The scope was that of evaluating the efficiency of the workflow under two different conditions.Scenario 1 represents a conservative situation in which the true value of the field cumulative oil production at the end of the forecast phase is very close to the medium value of the prior distribution.Conversely, in scenario 2 the prior field cumulative oil production was biased relative to the true case.
A total of twelve years of historical production data were considered available.Four of them were used for history matching; the remaining eight were used to check the accuracy of the forecast phase.
Four different types of observations were considered during the assimilation step: Well Bottom-Hole Pressure (WBHP), Well Block Pressure (WBP), well water cut (WWCT) and Well Oil Production Rate (WOPR).The well block pressure was used as an approximation of the well pressure measurements under static conditions, when the wells are shut-in.Although the wells were supposed to be produced at a constant rate, the actual WOPR values varied because in some simulation models wells were unable to respect the target value.
The data used for assimilation were 318 WBHP, 66 WBP (pseudo static pressure), 264 WCT and 264 WOPR values.All the data was subject to noise and the noise was assumed to be white.The imposed standard deviations are: 3 bar for pressure; 1% for water-cut; and 0.0001 sm 3 /d for oil production rate.
During the history matching phase producers are set to honor the historical oil production rates, with a minimum operational BHP.For the forecast phase a constant field liquid rate constraint of 1150 m 3 /day was imposed to the producers.The injector is set to replace the volume of fluid produced.The well production rate was reduced by a factor of 0.9 whenever the water cut exceed 0.9.

Facies Probability
The porosity mean was converted to a probability map that described the channel occurrence probability in each gridblock of the model.The transformation from mean porosity values to channel probability was expressed as a piecewise function (Fig. 12).
Defining / ch and / fp , the mean porosity values of channels and floodplain, respectively, and r ch and r fp their corresponding standard deviations, the transformation function was defined as: Such map is needed to construct the channel probability map which helps assessing the presence of channels through its key parameters such as the mean value and standard deviation.An example will be shown in Figure 27 in Section 2.10.

Results: Scenario 1. History Matching Objective Function Evaluation
The mean values, the standard deviation and the ranking of the objective function of each facies realization The heat map relative to the assimilation phase for all the facies realizations is shown in Figure 14.A logarithmic color scale was used to qualitative distinguish different OF values.Results show that the imposed Well Oil Production Rates (WOPR) were respected fairly well in all the cases, which means that the production wells did not reach the minimum bottom hole pressure.The quality of the water cut match varied depending on the well and facies realization.Production wells PROD1, PROD2, PROD7 and PROD8 showed a good match for almost all the cases.The C3 and B3 facies realizations showed a better overall  match of the water cut than the other facies realizations.In general, the bottom hole pressure and well block pressure showed the worst OF values.In the column orientation, facies realization A3 and C1 showed a much more accurate response in terms of pressure than the other realizations.

Forecast Uncertainty
Results in terms of observed data (WOPR, WBHP, and WWCT) of the initial and calibrated models for well PROD3 are shown in Figure 15.The vertical line separates the assimilation and forecast periods.The true curve for each observed quantity is plotted as a reference.The gray and blue lines are the data obtained from the initial (prior) and final (posterior) ensembles, respectively.Case B3 provides a better match for WOPR and WWCT at well PROD3.Conversely, case B3 shows a less accurate reproduction of the reservoir pressure if compared to cases A3 and B3.The calibrated models honor fairly well the WWCT during the assimilation period leading to a sensible reduction in the uncertainty compared to the initial models in the forecast phase.Similar behavior is observed for the WBHP.
The heat map associated to the forecast phase for all the facies realizations is shown in Figure 16.During the forecast phase, the total liquid rate was imposed, which means that the oil rate varies as a function of the water cut.Consequently, the values of the OF for the WOPR are well correlated to the OF values of WWCT and are much higher than in the assimilation phase.The realizations A3, B3, C1 and C3 show better pressure matches than the rest of the realizations.
The values of the OF, based on the cumulative oil production at the end of the forecast phase, for a progressively larger number of facies realizations are shown in 1 000 2 000 3 000 4 000 Time (days) 0 1 000 2 000 3 000 4 000 Time (days) 0 1 000 2 000 3 000 4 000 Time (days) 0 1 000 2 000 3 000 4 000 Time (days) 0 1 000 2 000 3 000 4 000 Time (days) 0 1 000 2 000 3 000 4 000 Time (days) 0 1 000 2 000 3 000 4 000 Time (days) 0 1 000 2 000 3 000 4 000 Time  The correlation between the OF rank of the facies realizations in the forecast phase and in the assimilation phase is shown in Figure 18.The facies realizations that lie on the diagonal line x = y represent the realizations that did not change their ranking from one phase to the following.The realizations with a better rank in the forecast phase than in the assimilation phase are located in the upper part of the plot (y > x).Conversely, the realizations that undergo a ranking decrease in the forecast phase relative to the assimilation phase are located in the lower portion of the plot (y < x).This plot clearly shows the non-uniqueness of the solution because more than one reservoir configuration can provide a satisfactory reservoir response if compared with the observed behavior.Furthermore and even more importantly, realizations which seem fully acceptable after calibration are not necessarily representative of the true reservoir internal configuration and provide wrong production forecasts; on the other hand, realizations which have a relatively low ranking after history matching can reliably predict the reservoir behavior.
The field cumulative oil production (FOPT) distribution obtained from the prior and posterior ensembles for each conceptual facies model is shown in Figure 19.The colored regions identify the [P10, P90] intervals and the dotted colored line indicates the mode (P50).Cumulative oil production at the end of the forecast phase grouped OF ascending.Forecast vs assimilation OF ranking.
The continuous line represents the true value.The uncertainty of the cumulative oil production is efficiently reduced in all the cases, but it is worth noting that the facies models of type B show the best performance in reducing such uncertainty.In all models, the FOPT mean is very close to the true value for both the prior and posterior distributions.
The distribution of the FOPT is vital information for the evaluation of a development plan for a given reservoir.However, another important aspect is the field cumulative water production (FWPT), mainly because it has a great impact on the required surface equipment and facilities and, as a consequence, on the economic feasibility of a project.The FWPT distribution obtained from the prior and posterior ensembles for each conceptual facies model is shown in Figure 20.Similarly to the FOPT, the FWPT uncertainty shows an evident reduction.The conceptual model B shows the larger reduction relative to the prior distribution.Conversely, model A is less effective to reduce the uncertainty.

Facies and Petrophysical Properties
The evolutions of the mean and standard deviation of the ensemble permeability for cases A3, B3 and C1 are shown in Figures 21, 22 and 23, respectively.The petrophysical properties generation was not constrained at the well locations, thus the corresponding standard deviation is greater than zero.The petrophysical distribution of model B3 seems to reproduce in a more realistic way the petrophysical properties of the true case.
Similarly, the evolutions of the mean and standard deviation of the ensemble porosity for cases A3, B3 and C1 are shown in Figure 24, 25 and 26, respectively.Prior (gray) vs posterior (blue) cumulative oil production distribution at the end of the forecast phase for each conceptual models group, A, B and C, respectively.Black solid line is the true value and dotted lines the mean values of prior and posterior distributions.Outlined regions represent, instead, the amplitude interval between the P10 to P90 of the corresponding distributions.
Case A3 shows unrealistic high values of porosity in the southern channel, which is a clear indication that a wider channel exists, i.e. more volume is needed in that zone, as confirmed by comparing it to the true case.Case B3 shows a fairly good match of the porosity field and lower values of standard deviation.Localized spots of higher values of standard deviation exist in zones with large porosities, indicating a higher uncertainty in these zones.Case B3, indeed, is characterized by the same facies model of the true case.As expected, the assimilation process results in the most satisfying petrophysical characterization even though the reference porosity and logpermeability distribution were not included in the ensemble (Fig. 22,25).
The mean and standard deviations of the channel probability field are shown in Figure 27.Red zones in mean probability map represent high probabilities of channel occurrence and blue zones low probability of channel occurrence.This was obtained by applying Equation ( 12).The map seems to capture fairly well   the presence of channels in the central region of the reservoir, which is related to the presence of several wells providing the dynamic observations used to calibrate the models.The obtained facies probability map could be used as soft data for the generation of a new set of facies models in an iterative procedure so as to better represent the reservoir internal geometry.Eventually, the red areas of the standard deviation map show the zones of the reservoir where the information on the probability of channel occurrence is less reliable.True, mean and standard deviation of the channel probability field considering all the cases.

Results: Scenario 2. History Matching Objective Function Evaluation
The mean values, the standard deviation and the ranking of the objective function of each facies realization are shown in the boxplot and table of Figure 28.
The heat map relative to the assimilation phase for all the facies realization is shown in Figure 29.A logarithmic color scale was used to qualitative distinguish different OF values.Results show that not always the imposed WOPR were respected, which means that the production wells in some cases reached the minimum bottom hole pressure imposed at the wells.The water cut match quality varies depending on the well and facies realization.However, wells PROD8 and PROD6 present the best and the worst WWTC match, respectively.The C1, B3 and B1 facies realizations show a better overall match with respect to the other facies realizations.

Forecast Uncertainty
Results in terms of observed data (WOPR, WBHP, and WWCT) of the initial and calibrated models for well PROD3 are shown in Figure 30, for cases A2, B3 and C1 respectively.Results are shown for the models of each conceptual facies model with the best OF value, i.e.A2, B3 and C1.The vertical line separates the assimilation and forecast periods.The true curve for each observed data is plotted as a reference.The gray and blue lines are the data obtained from the initial (prior) and final (posterior) ensemble, respectively.Cases A2 and B3 match fairly well the historical water cut values both in assimilation and forecast phase.The WOPR have almost perfect matches in the assimilation phase and acceptable matches for the forecast phase in all the models except for well PROD2 in case C1.Cases B3 and C1 match accurately the pressures in the assimilation and forecast phases.Conversely, case A2 fails to reproduce the bottom hole pressure during the assimilation phase in both wells.The heat map associated to the forecast phase for all the facies realization is shown in Figure 31.Low values of WBHP for the injection wells are obtained; in these cases the injection pressure reached the maximum operational value so the wells change to a pressure controlled condition.The oil rates show higher mismatch than the rest of observations.In general, the realization B3 shows the lowest mismatch values.
The values of the OF, based on the cumulative oil production at the end of the forecast phase, for a progressively larger number of facies realizations are shown in Figure 32.The corresponding evolutions of the P90, P50 (or median) and P10 are indicated by a red, green and blue line, respectively.The FOPT uncertainty increases when cases of higher OF values are considered.Cases with lower OF values exhibit a FOPT distribution in which the true value lies outside the P10-P90 range.On the other hand, when all the available models are accounted for, a wider FOPT (P10-P90) range is obtained, but the true case more realistically lies within the P10-P90 range.
The correlation between the OF rank of the facies realizations in the forecast phase and in the assimilation phase is shown in Figure 33.The facies realizations that lie in the diagonal line x = y represent the realizations that did not change their assimilation ranking from one phase to the following.The realization with a better rank in the forecast phase than in the assimilation phase are located in the upper part of the plot (y > x).Conversely, the realizations that experiment a ranking decreasing in the forecast phase with respect to the assimilation phase are located in the lower of the plot (y < x).In this scenario as well, realizations which seem fully acceptable after calibration are not necessarily representative of the true reservoir internal configuration and provide wrong production forecasts; on the other hand, realizations which had a relatively low ranking after history matching can reliably predict the reservoir behavior.
The field cumulative oil production (FOPT) and the field cumulative water production (FWPT) distribution obtained from the prior and posterior ensembles are presented in Figure 34 to Figure 37.The colored regions identify the [P10, P90] intervals and the dotted colored line represents the median (P50).The continuous line represents the true value.A remarkable increase is observed both in the precision and accuracy of the estimation, i.e. a reduction of the distribution spread from the initial to the final models and a shift of the mean of the distribution towards the true value.Thus, the uncertainty associated to the cumulative oil and water production is efficiently reduced.It is worth noting that although case B and case C generally show better OF values they fail to match the true FOPT value within the uncertainty range P10-P90 (Fig. 34).

Facies and Petrophysical Properties
The evolutions of the mean and standard deviation of the ensemble permeability and porosity cases B3 and C1 are presented in Figure 38 to Figure 41.Case B3 and case C1 were selected because they were the models that showed the lowest objective function.In this case as well, the petrophysical distribution of model B3 seems to reproduce in a more realistic way the petrophysical properties of the true cases.Furthermore, the standard deviation for the B3 in both porosity and permeability is substantially lower than in the C1 case.
Prior vs posterior cumulative water production at the end of the forecast phase for each facies group.
Results also showed that even the adoption of a conceptual model for facies distribution clearly representative of the reservoir internal geometry does not guarantee reliable results in terms of production forecast, thus the statistical approach seems the only way to reduce the uncertainty inherent to reservoir modeling.In addition, the models that showed the best OF values in the assimilation phase were not necessarily able to describe the system behavior in the forecast phase.When the true value of FOPT lays close to the mean FOPT value of the prior distribution the uncertainty of the FOPT remains almost constant, independently of the number of facies realizations considered for the forecast.Conversely, when the prior distribution is biased in respect to the true value, the uncertainty of the model strongly depends on the number of realizations considered for the forecast.
Recently, progresses were made in the automatic update of the facies maps; however, these researches are still under development and not available for application by reservoir engineers in everyday life.Conversely, the proposed approach/workflow is based on wellconsolidated techniques and is virtually ready for use by means of available commercial tools, thus it constitutes a step forward for future technologies capable to automatically integrate different facies modeling in the history matching process.To this end the workflow could be extended to use the facies probability maps as soft data for the generation of new facies realizations to further improve the quality of the match in an iterative process.However, it should be mentioned that the computational cost for large models could be relevant. Figure1 Figure 2 Traditional workflow a) vs closed loop workflow b) for reservoir modeling.
Figure 3 Figure 5Porosity and permeability histograms for the channel and floodplain facies.

Figure 8
Figure 8Assisted history matching workflow.

Figure 9
Figure 9Conceptually different vs conceptually similar facies models.

2. 5
Figure 10Set of facies realizations based on the different conceptual models.
Figure 11Examples of petrophysical realizations from different facies models.
Figure 13Ranking of the objective functions for each facies realization.
Figure 15Prior and posterior observations vs time for well PROD3 for cases A3, B3 and C1.

Figure 17 .
Figure 17.The corresponding evolutions of the P90, P50 (or median) and P10 are indicated by a red, green and blue line, respectively.It is worth noting that the mean and P90 values remain almost constant independently of the considered realizations.In contrast, an overestimation of the WWCT values in the cases with high OF values leads to lower field cumulative oil values and, consequently, to lower P10 values.The correlation between the OF rank of the facies realizations in the forecast phase and in the assimilation phase is shown in Figure18.The facies realizations that lie on the diagonal line x = y represent the realizations that did not change their ranking from one phase to the following.The realizations with a better rank in the forecast phase than in the assimilation phase are located in the upper part of the plot (y > x).Conversely, the realizations that undergo a ranking decrease in the forecast phase relative to the assimilation phase are located in the lower portion of the plot (y < x).This plot clearly shows the non-uniqueness of the solution because more than one reservoir configuration can provide a satisfactory reservoir response if compared with the observed behavior.Furthermore and even more importantly, Figure16OF heat map of the forecast phase.
Figure 18 Figure 20Prior (gray) vs posterior (blue) cumulative water production distribution at the end of the forecast phase for each conceptual models group A, B and C, respectively.Black solid line is the true value and dotted lines the mean values of prior and posterior distributions.Outlined regions represent, instead, the amplitude interval between the P10 to P90 of the corresponding distributions.
Figure 21Comparison between the true ln(PERMx) field for layer 1 and ensemble mean in case A3 at different assimilation time (initial, at day 390, at day 870, at day 1 350).
Figure 22Comparison between the true ln(PERMx) field for layer 1 the ensemble mean in case B3 at different assimilation time (initial, at day 390, at day 870, at day 1 350).
Figure 23Comparison between the true ln(PERMx) field for layer 1 and ensemble mean in case C1 at different assimilation time (initial, at day 390, at day 870, at day 1 350).
Figure 25Comparison between the porosity field for layer 1 and ensemble mean in case B3 at different assimilation time (initial, at day 390, at day 870, at day 1 350).
Figure 28Objective function and ranking for each facies realization.
Figure31OF heat map of the forecast phase.