Integration of structural uncertainty and experimental design for uncertainty quantification

During the history matching, attentions are shifted from the structural map and focus are on the dynamic properties such as relative permeability data, saturation end points and Pressure-VolumeTemperature (PVT) information. However, structural maps can inherit huge uncertainty which can affect the dynamic behaviour of the reservoir. The objective of this study is the presentation of an approach for quantifying uncertainties in production forecast by incorporating uncertainty inherent in structural maps. Three static models (a “Low”, “Mid” and “High” cases) were developed, history matched using traditional method and used for production forecast. Two (2) possible alternatives (replicated and conventional) were considered for uncertainty analysis. This yielded four (4) reservoir proxies using experimental design techniques. Uncertainty was quantified using validated proxies with Latin hypercube sampling in a Monte Carlo simulator. Conventional and replicated approach yielded non unique answer. Replicated method which assume equal probability for the three model cases and incorporated them in the design offeres a more realistic answer. The analysis revealed with 90% certainty that at least a normalized reserve figure of 0.91 will be realized and 10% confidence that at most reserve figure of 1.039 will be realized. The method is recommendable for developing new noncomplex reservoirs/fields where subsurface data acquisition is a challenge.


INTRODUCTION
The use of uncertainty and risk analysis in the petroleum industry has increased since the first oil crisis.Uncertainty can emanate from numerous sources and directly impact the physical reservoir description (Vanegas et al., 2006).In the geophysical domain uncertainties are linked to acquisition, processing, and interpretation of seismic data (Vincent et al., 1999;Corre et al., 2000).Some of these uncertainties are due to errors in picking horizons, difference between several interpretations, uncertainties and errors in depth conversion, uncertainties in seismic-to-well-tie, uncertaintes in pre-processing and migration, and uncertainties in the amplitude map of the top reservoir.For efficient reservoir management, it is important to understand these sources for uncertainties to be properlly accounted for during decision making since the actual reservoir description cannot be achieved even at the end of the field life (Mohaghegh et al., 2006).
In spite of this, static models, perfect or imperfect, have been helpful in supporting future development activities of hydrocarbon reservoirs.The decision to sanction an oil or gas field development are taking with consideration of the possible upside and downside outcome (Smith et al., 1993).These quantities require understanding of the uncertainties of the reservoir's static architecture and of its dynamic behaviour during production (Friedmann et al., 2003).Recent research and applied reservoir characterization and modelling projects have focused on the application of design of experiment (DoE) and Response Surface techniques as a basis for history matching and uncertainty analysis (Almeida et al., 2003;Amudo et al., 2008, Dai et al., 2014a, b, c).Detailed review of DoE can be found in Morris (2000).
Assumption of single deterministic static model for history matching is a common practice in the industry.After the history matching process, experimental design and Monte Carlo simulation are usually combined to generate important probabilities that correspond to topdown figures with their corresponding values of parameters.In most cases the extracted parameter values give different figures when validated using the numerical simulation.This observation are due to assumptions such as the distribution functions for different uncertainty factors and the integrity of the response surface models, an input to Monte Carlo simulation which combined these factors.
Generation of multiple realizations of the geological models allows quantification and minimization of uncertainties associated with geological information.However, this geostatistical method is inefficient due to difficulty in conditioning the entire models with dynamic data (Alhuthali et al., 2006).These multiple geostatistical realizations can be ranked to reduce the number of models.The idea of ranking realizations has been around for a while and was first published in the context of stochastic reservoir modeling as pointed out by Ballin et al. (1992).The application of most of the available ranking measures when dealing with multiple reservoir models can be time consuming.In these situations, the best approach is to minimize the number of geological models without compromising the evaluation process.This study presents an approach for quantifying uncertainties in production forecastsusing three (3) models obtained from structural map uncertainty and provides new approach for generating representative reserves probability figures for investment decision making.

Review of static modeling
The summary of the methodology adopted for the static modeling and uncertainty analysis is shown in Figure 1.The only structural map available is that generated for a most likely case.This was Olatunde et al. 75 used for building the base case model.Two additional models, the "High" and "Low" caseswere built by using ±10 ft uncertainty on the original structural map.Geocellular frameworks were built using the three structural maps.The cell size (ΔX*ΔY*ΔZ) used is approximately 50×50×1ft which was assumed to have captured all the intra-reservoir heterogeneities.A sedimentology concept was developed for the distribution of volume of shale (Vshale) at different layers.Based on the established concept, spherical variogram model (major range of 2500, minor range of 1000, vertical distant of 100 ft, Azimuth of 60° and sill and nugget of approximately 1 and 0) was adequate for Vsh distribution at different para-sequences/zones. Since the facies model is absent for this study, the distribution of porosity was achieved using Sequential Gaussian Simulation.The simulation was constrained by porosity variogram and corresponding Vshale distribution map for various zones.
Figure 2 shows average property maps to illustrate the variability of the properties within the reservoir for a realization of the reservoir.The volume of shale distribution shows an increase from North to South.At the crestal part, the Vsh ranges between 0 and 10% which justified the location of the drainage points.The water saturation map shows increasing water saturation from the north towards the south with minimum of 19.2% at the crest and value as high as 67% towards the south.Both porosity and net-to-gross properties show high values in the upper structure and decrease from north to south.Permeability distribution is similar to the porosity distribution with minimum and maximum values of 13 and 2154 mD, respectively.

Dynamic modelling
The model dimensions in the X-, Y-and Z-direction are 197, 75 and 9 m, respectively.This amounts to 132,975 active cells on a fine scale equivalent to 98*37*9 grid dimension (32,634 3D grid blocks) on coarse scale.The upscaling from fine to coarse scale was done on the basis of available computer for simulation and to reduce the computational time from 5 h to about 1 h using a fine and coarse grids respectively.However, this was able to preserve certain geologic features in certain layers in the static model.The models were initialized under static condition of 4951 psia, 170°F, API gravity of 35°API and average viscosity of 1.3 cp at reference depth of 11142 ft.
Figure 3a is a cross section showing the direction of water encroachment.This is known by the relative positions of the oil water contact (OWC) at the start (November, 1993) and end (November, 2008) of history.So water encroaches from the flanks.Material balance calculation was carried out using MBAL software and with Fetkovich analytical model, the reservoir energy was captured.Analytical plot in Figure 3b shows clearly that the observed production cannot be matched without the aquifer in place.The aquifer influence function in Figure 3c revealed the ratio of the aquifer radius to the reservoir radius (aquifer size) of 8 which is an indication of appreciable aquifer which is likely to react later to voidage caused by oil, water and gas production through the wells.The energy plot in Figure 3d revealed the drive mechanism to be majorly a combination of water influx and pore volume expansion.Other aquifer characteristics including permeability, encroachment angle and tank radius are shown in Table 1.History matching of the three models was achieved using experimental design assisted traditional method.

Relative permeability to water
There was no relative permeability measurement available.The measurement from analogous reservoir was modified and used for history matching (Arinkoola et al., 2015).Keeping end points constant, a low and a high case water relative permeability were derived using some multipliers on the original curve.

Critical water saturation
Critical water saturation was globally increased by 25% to match the water breakthrough time.The connate water saturations obtained from the static model were relatively low.

Oil viscosity
The key oil properties available after carefully reviewing the PVT and sampling depth data for this reservoir include; gas-oil-ratio (GOR), gas gravity, oil API gravity and oil formation volume factor (FVF).PVT model was developed using empirical correlations.Oil viscosity represents huge uncertainty.Error of ±10% was assumed for the resulted oil viscosity property.This gives rise to three viscosity models, minimum viscosity model (degrading original viscosity by 10%), base case and maximum viscosity model (increase original viscosity by 10%).

Vertical/horizontal permeability
KV/KH is one of the major uncertainties since there is no core measurement for kv/kh ratio supervision.At the beginning, the vertical permeability model (PERMZ) was calculated using the vertical/horizontal permeability ratio of 0.1.This was guided from the information about a neighboring field with KV/KH ratio ranges from 0 to 0.25.However, the minimum KV/KH that preserves history matching ranges between 0.05 and 0.12.

Original oil-water contact
The first well drilled in 1973 penetrated the sand without encountering the water-up-to (WUP) depth.However, the water-upto from two simultaneously drilled wells indicated slightly varying depths.The average Oil-Down-To (ODT) from these wells is approximately 11239 TVD ftss.The base case match was obtained at 11240.5TVD ftss, for the low case model, a good match was obtained with ODT at 11245.5TVD ftss while for the high case, the ODT at 11238TVD ftss was adequate to obtain a good pressure and saturation matches.
A good match on water cut and pressure was obtained at reservoir level.Figures 4 to 6 shows the outcome of the history match for the three models.Although the water breakthrough time for high case was earlier than actually observed, the overall water cut matches, are relatively good.

Uncertainty analyses
The uncertainty analysis was performed by treating the two extreme models (low and high cases) as experimental replicates.Thus the experiment was performed three times using a full factorial 2-level experimental design.This approach can be likened to situation when a random error is introduced in the experiment.Usually, random error in any numerical simulation is always assumed zero since repeating numerical simulation time without number does not change the fundamental equation with which calculations are made except if the input data is changed.With this method the three reservoir models are assumed to have equal likelihood of representing the actual reservoir.A total of 96 experiments were conducted and analyzed considering uncertain parameters shown in Table 2.
Table 2 shows the range of values (low and high as multipliers on base case values) of parameters used in the design.The design matrix consisted of the parameter values in coded form as well as the actual response values in stock-tank barrels after experimentation is shown in Table 3.In all the experiments, water cut and GOR were constrained.To select the best model that fit the experimental data, linear, quadratic, polynomial, cubic and factorial models were evaluated for their significance (α = 5%) using Fstatistics, t-statistics and correlation coefficient criteria.This was done using the analysis of variance (ANOVA) package available in  Figure 7 shows half normal plot that shows interactively, how the factors were selected for analysis during the analysis using the design expert.
From this plot, the main factors, horizontal permeability (A), vertical permeability (B), water relative permeability

Conventional method
Assumption of single deterministic static model for history matching is a common practice in the industry.After the history matching such model uncertainty are quantified on it and cdf is generated and extraction of important probabilities that correspond to topdown figures are made with the parameter values that correspond to them.Most often, extracted parameter values give different figures when simulated numerically.The discrepancy could be attributed to the assumption of the distribution function but largely to quality of the response surface developed.This work established comparison between the four possible reservoir proxies in terms of quantified uncertainty figures.In addition to the proxy Equation 1, the three other reservoir proxies were developed using the CCD method.Central composite design belongs to the family of response surface methods.Also in this family are 3-level algorithms consist of Box-Behnken, Full factorial, D-optimal e.tc.Central composite design is a 5-level design consists of a factorial design with corners at ±1 of the cube, augmented by additional "star" and "centre" points which allows for efficient estimation of the secondorder polynomial equations.Because it requires minimum number of experiments and its ability to mirror beyond the specified range of parameter using the "star" point denoted by ±α CCD is selected for response surface modeling ahead of others.
The parameter used excluded the water relative permeability (KRW) being insignificant from previous analysis; and the shows the new data summary for CCD experiment.CCD is an improved algorithm due to its ability to mirror and integrate factor values outside their defined ranges.Table 5 shows the design matrix for CCD comprises in total 75 experiments using "high", "mid" and "low" case models.

Analysis of variance
As shown in Table 6, the Model F-value of 116.05 implies the model is significant.Thus the chosen linear and interaction factorsare more than adequate for a good model.There is only a 0.01% probability that a "Model F-Value" this large could occur due to noise.Values of "Prob > F" less than 0.0500 indicate models terms are significant.
In this case A, B, D, E, AE and BE are significant only significant model terms for the response.Thus factor C (KRW) is not signficant on the model.The regression equation relating ultimate recovery (dependent variable) and the history match variables (independent variables) is Equation (1).The response was transformed into logarithmic form to address nonlinearity issues.
(1)  Similarly, after performing ANOVA, the significant parameters were used to formulate the response surfaces which are the regression equations relating ultimate recovery and the independent variables.
The model equations are: (2) (3) (4) The values of coefficients in the response surface equations are presented in Table 7. Figure 8 is the crossplots of the actual experimental values against the predicted values.If the prediction by these models were a perfect fit of the experimental data then all of the points would lie on the .On all these plots, it is clear that the dependent and independent variables are not completely linearly related.This is responsible for some of the data points observed not to aligned completely on the straight line.There are more data points on the replicate cross plot because it is a conglomerate of the three models.The variation of the data points on the remaining other cross plots however indicated some degree of uncertainty in the three models.With some degree of confidence nevertheless, the developed response surface modelsare excellent representation of reservoir and agreed amenable for sampling within the experimental range of parameters for further analysis.

Factors interaction
The ANOVA table shows clearly that all the selected main factors have significant effects on the ultimate recovery except the water relative permeability curve (KRW).Nevertheless, the interaction of some of these factors exhibits significant effect on the model development.diagrams showing the main effects of OVISC, PERMX, PERMZ, Qwinj and their interaction on the ultimate recovery.A close observation of the base contour plot in all these diagrams shows that main factors OVISC, PERMZ, PERMX and Qwinj independently have influence on the response.For instance, as can be observed in Figure 9, as the viscosity increases from minimum to its maximum value, the amount of oil recoverable was observed to be decrease.This observation is expected because viscosity controls the mobility of oil and as such as oil viscosity increases, oil mobility reduces so also recovery.More recovery was obtained while keeping Qwinj constant and increasing PERMX from minimum value to maximum.This observation is also true indicating the efficiency of the water flooding as Qwinj was set at its maximum value.Similarly, a low recovery of oil was obtained keeping PERMX constant at minimum while Q winj is at minimum.With PERMX increases to its maximum value, a proportionate increase in recovery was observed as water injection rate increases.This was due to preferential increase in flow of oil aerially relative to water and gas.Thus more residual oil was contacted and swept more efficiently.Maximum recovery was obtainable at low OVISC and increasing Q winj .At high Q winj , the oil recovery decreases as the OVISC increases from minimum to maximum.This observation indicates that oil viscosity is a huge uncertainty and the result further revealed the need to acquire more realistic PVT data.

Monte Carlo simulation
Figure 12 shows the representative distributions from the Latin Hypercube sampling comprises of 100,000 samples using different proxy models.Figure 13 is a bar chart showing P10, P50 and P90 realizations using proxy Equations 1 to 4 compared with when the 3 history matched models was used directly for prediction without subjected them to experimental design (direct simulation).The bar chart shown in Figure 13 is normalized with regards to the P50 value (34.88 Mstb).This normalization was done for easy interpretation of the result.The P50 valueis the total forecast reserves using the deterministic figure obtained from direct forecast assuming base case/mid case model as shown in Figure 14.
The Monte carlo simulation was constrained on the mean value of reserves for all the models.In other words, the assumption of the distribution functions for different parameters was highly dependent on the P50 realization after the simulation using different models.The mean reserves was ensured approximated the deterministic figure of 34.88 Million barrels.Therefore the base case   normalized value equals one (unity).
The analysis was done by considered extreme percentiles: 10-percentile (denoting a 90% probability that the normalized reserves figures will at least greater than the 10-percentile normalized figure) and 90 percentiles (denoting a 10% probability that the normalized reserves figure will at most greater than the 90-percentile figures).These figures are useful for investment decisions.To assess the degree of uncertainty and quantify them, the disparity between the normalized P10 and P90 figures can provide managers with information on weither to acquire more data for detailed study or adoption of additional enhanced recovery techniques.
Considering the Mid/base case model, the P10 and P90 normalized value are 0.93 and 1.02 respectively.This figure indicates approximately disparity of three million stock tank barrel (3.2 Mstb) assuming no uncertainty (deterministic) in the structural map.For low case model, the P10 and P90 normalized value are 0.957 and 1.02 respectively.This figure indicates approximately disparity of 2.2 million stock tank barrel (2.2 Mstb) assuming the base case map was interpreted 10 ft below the base map (-10% uncertainty) .For high case model, the P10 and P90 normalized value are 0.961 and 1.042 respectively.This figure indicates approximately disparity of 2.8 million stock tank barrel (2.8 Mstb) had it been that the base case structural map was interpreted 10ft above the base map (+10% uncertainty).These figures represent uncertainty in the production forecast when compare with the base case.The top and base of the sand if not properly picked can introduce huge uncertainty.In some places, the uncertainty due to structural map can be as high as ±50 ft this can have a huge impact on economic decisions.
The convention in the industry is to use the base case model then quantify uncertainty using experimental design and Monte Carlo simulation.By this approach, the P10 and P90 normalized value are 0.89 and 1.077, respectively.This figure indicates approximately disparity of 6.5 million stock tank barrel (6.52 Mstb).However, the proposed approach (the replicates method) which we demonstrated by the integration of the three distict uncertainty models (low, mid and high)resulted to another composite uncertainty model upon which uncertainty analysis was based.In this method, the low and the high case models are regarded as replicates of the base case which resulted in model Equation 1.For FOPT Repl , the P10 and P90 normalized value are 0.91 and 1.04, respectively.This figure indicates approximately disparity of 4.5 million stock tank barrel (4.53 Mstb) assuming ±10 ft uncertaity in the structural map.

SUMMARY AND CONCLUSION
The structural map uncertainty has been successfully used as a basis for generating stochastic models for production forecast in a case study of a Niger Delta marginal field.The three models developed ("Low", "Mid" and "High" cases) show marked variation indicating that the top structural map has uncertainty inherent in it.Experimental design and response surface methodology was sucessfully applied to quantify structural map and other petrophisical uncertainties.The uncertainty analysis was achieved considering three possible scenerio: conventional approach which utilized only the base case model, the second method which assumed equal likelihood for the three models and thus treated them as possible realizations, the replicate method that assumed extreme models as replicates of the base case model.
The result of assumption of uncertainty as small as ±10 ft on the base case map suggested that the map has inherent uncertainty.The top and base of the sand if not properly picked can introduce huge uncertainty and capable of impact decision on investment and overall management of the reservoir.From the analysis carried out, the degree of heterogeneity within the reservoir was captured by the disparity between the P10 and P90 reserves figures.In all the cases the discripancy varied between 2.2 and 6.5 million stock tank barrels with the upper extreme corresponds to the conventional method and the lower extreme the low case model.The difference suggested that the reserves figure can vary from one company to another depending on the analysis method.The replicate method that was introduced in this study however, indicates approximately the disparity of 4.5 million stock tank barrels.Although the conventional method overpredicted the reserves, it can be concluded that the reservoir is not too heterogeneous due to the proximity of these figures and injection of polymer instead of water injection can actually improve the recovery because the early water breakthrough recorded in some of the wells could actually be delayed.

Figure 1 .
Figure 1.Workflow for static and uncertainty analysis.

Figure 2 .
Figure 2. Average property maps for a reservoir realization.

Figure 3 .
Figure 3. Reservoir cross section and material balance plots.

Figure 4 .
Figure 4. History match plots of reservoir pressure, GOR, oil rate and water cut for base case model.

Figure 5 .
Figure 5. History match plots of reservoir pressure, GOR, oil rate and water cut for low case model.

Figure 6 .
Figure 6.History match plots of reservoir pressure, GOR, oil rate and water cut for high case model.

Figure 7 .
Figure 7. Half normal plot depicting the manual process that led to selection of model and significant model terms.

Figure 8 .
Figure 8. Cross plots of actual experimental data versus predicted value for all models.

FigureFigure 11 . 3 -
Figure 11.3-Dimensional diagram for the effect of water injection rate (Qwinj), PERMX and their interaction on ultimate recovery.

Figure 12 .
Figure 12.Oil reserves histograms from Monte Carlo simulation using different proxies.

Figure 13 .
Figure 13.Normalized bar chart for P10, P50 and P90 realizations from simulation and developed reservoir proxies.

Figure 14 .
Figure 14.Reprisentative models analogous to P10, P50 and P90 based on simulation approach from three history matched models.

Table 1 .
Aquifer properties for field pressure history match.
the Design Expert 6.0.Linear model augmented with the factorial model was selected having exhibited highest F-value, p-value (<0.0001) and correlation R-square of 99.6%.To select the impactful factors on the model, ANOVA was carried out.

Table 2 .
Uncertainty factors, their descriptions and range of values as multipliers. S/

Table 4 Table 3 .
Factorial design matrix and corresponding uncertainty recoveries for three uncertainty forecast models.

Table 4 .
Central Composite Design summary showing factors and range of values as multiplier on the base case value.

Table 5 .
Central Composite Design matrix and corresponding responses for three uncertainty forecast models.

Table 6 .
Analysis of variance table showing significance of model terms for replicate model.

Table 7 .
Value of coefficients in response surface correlations in equations 2, 3 and 4.