 Home
 Conferences
 Conference Proceedings
 Conferences
ECMOR X  10th European Conference on the Mathematics of Oil Recovery
 Conference date: 04 Sep 2006  07 Sep 2006
 Location: Amsterdam, Netherlands
 ISBN: 9789073781474
 Published: 04 September 2006
1  50 of 78 results


Modeling of Nucleation Fronts during Depletion of GasSaturated Porous Media
Authors S. Zaleski, M. Chraibi and F. FrancoOil production through cold depletion leads to degassing of the light species and the formation of a bubbly phase, sometimes called the foamy oil effect. This bubbly phase is particularly observed with heavy oils, combining high viscosity and asphaltènes. We have modeled depletion experiments on laboratoryscale cores using a onedimensional model at the Darcy scale, describing the multiphase flow of the oil and gas. The oil and gas phases move with a classical relative permeability model.
Initially, the gas phase is in the form of very small bubbles, which are physically viewed as stabilized in crevices of the solid phase or as small germs surrounded by a coating skin.
When surface tension is taken into account in the phase equilibrium a Gibbs radius appears
so that bubbles grow rapidly above the Gibbs radius and collapse below it.
As a result, the solution of the partial differential equations describing mass conservation display nonlinear fronts that connect two regions of almost constant concentration and velocity, in a manner analogous to shock waves. We describe an asymptotic theory that allows understanding the formation of the fronts and connects them to the dynamics of bubble activation at nucleation sites.



Numerical Analysis of Foam Motion in Porous Media Using a New Stochastic Bubble Population Model
Authors P. L. J. Zitha, D. X. Du and F. J. VermolenWe present a numerical analysis of the stochastic population balance (SPB) theory for foam motion in porous media. The theory condenses into a set of nonlinear partial differential equations in the saturation, pressure and bubble density. We solve the equations using the IMPES method and perform sensitivity analysis. Finally, we compare the saturation profiles obtained numerically with those obtained from CT scan foam experiments.



Underground Storage of H2 and H2CO2CH4 Mixtures
Authors M. Panfilov, G. Gravier and S. FillacierIn conditions of a natural irreversible reduction in worldwide oil production, the hydrogen is examined as one of the basic sources of the renewable energy able to replace the fossil energies. The underground storage of H2 is considered as a most perspective way to store huge amounts of hydrogen. Various industrial methods usually produces the H2 mixtures with a dominating presence of CO2 and CH4, so the objective of the present research was to analyse various hydrodynamic aspects related to the problem of storing the hydrogen mixtures in porous reservoirs. The storage of H2+CO2+CH4 mixtures like manufactured or civil gas produced by the industrial coal gasification technique represents a double industrial interest, by ensuring, in addition, to capture CO2. When stored in an aquifer, this gas mixture is subjected to an intense bacterial activity. The bacteria consume H2 and CO2 with generating new gases like CH4. Such transformation of the stored gas composition was observed in European civil gas storages. The new model of multicomponent gas injection in reservoir accompanied with chemical reactions and a bacteria population growth is developed, which was next analysed numerically and by using the mathematical methods of the theory of population dynamics. The characteristic parameters like the constants of Monod and Michaelis are obtained by fitting the results of modelling with the insitu experimental data. Various scenarios of the system evolution are revealed. In particular, an irregular oscillating regime was detected, in which various individual gases can create their owns space accumulations which move in time all along the reservoir. These results are confirmed by insitu observations. The results of modelling yield optimal regimes of gas injection. The third part deals with using CO2 as a cushion gas for a H2 storage, which may be examined as a new technique of storing large amounts of CO2 in order to reduce the climate changes in the atmosphere. Within the zone of contact between the gases the same chemicalbacteria interactions are observed which leads to a natural reduction of CO2 and a possibility to perform a permanent CO2 injection into the reservoir. Using the mathematical model developed, we have estimated the characteristic rate and the amount of CO2 reduction in time, which determine the main parameters of CO2 injection in practice.



Analytical Modeling of CO2 Sequestration and Enhanced Coalbed Methane Recovery (ECBM)
Authors C. J. Seto, K. Jessen, T. La Force and F. M. Orr Jr.Injection of CO2 into deep unminable coal seams is an option for geological storage of CO2. In many industrial settings, pure CO2 streams are expensive to obtain and a mixture of CO2 and N2 would be less expensive. New analytical solutions are presented for twophase, fourcomponent flow with volume change on mixing in adsorbing systems. Analytical solutions have been reported previously for singlephase threecomponent gas flow (Zhu et al.) and for multicomponent flow of incompressible fluids with adsorption (Johansen and Winther, Dahl et al., Shapiro et al.).
In this paper, we analyze the simultaneous flow of water and gas containing multiple adsorbing components by the method of characteristics. Mixtures of N2, CH4, CO2 and H2O are used to represent the ECBMflue gas process. The displacement behavior is demonstrated to be strongly dependent on the relative adsorption strength of the gas components.
N2 and CO2 recover CH4 through different mechanisms: CO2 preferentially adsorbs onto the coal surface, resulting in a shock solution; while N2 displaces CH4 by reducing the partial pressure of CH4, resulting in a rarefaction solution. When mixtures of N2 and CO2 are injected, the displacement exhibits both shock and rarefaction features. For CO2rich flue gas, a path that includes a switch between branches of nontieline paths is observed, a feature not previously reported for gasliquid displacements. In these solutions, an additional key tieline, at which the switch between nontieline paths occurs, is required. In the shock along this tieline, the nontieline eigenvalues of injection and initial segments and the tieline shock velocity are equal.
Analytical solutions to ECBM processes provide insight into the complex interplay among adsorption, phase behavior and convection. Improved understanding of the physics of these displacements will aid in developing more efficient and physically accurate techniques for predicting the fate of injected CO2 in the subsurface.



Measures of Efficiency for Assisted History Matching
Authors G. J. Walker and S. PettigrewThe recent advances in computer assisted history matching have enabled the asset team to investigate multiple alternative reservoir descriptions (SPE89974, Williams et al). The systems are being asked to work in a high number of dimensions, and yet we know that the problem is tractable as we are able to find models that satisfy history match criterion. What we need is a measure of efficiency or elegance to the finding of the alternative solutions, to then allow optimization of the search and an objective discussion of the way in which different strategies interact with the task. This paper covers two case studies, at different stages in their lifecycle, and how different choices of genetic algorithm parameters modify the efficiency.



Rapid, Stochastic Updating of Reservoir Models to Dynamic Data – An Evaluation of PField Approach
Authors S. T. Reinlie and S. SrinivasanContinuous dynamic data such as well flowing bottomhole pressure carry information that characterizes reservoir heterogeneity. A novel approach to analyze continuous monitoring pressure data and to update reservoir models based on incremental information is presented. First, the pressure transient data is analyzed to identify the size and shape of permeability heterogeneity in the presence of fluctuations in rate and pressure. Unlike the complicated pressure or rate deconvolution algorithms presented in the literature, a simple semianalytical approach is presented here that attempts to reconstruct the bottom hole pressure response after removing the effect of rate fluctuation. Using the reconstructed pressure profile, estimates of the radius to the boundary of the heterogeneous region and the effective average permeability are obtained by applying a simple optimization procedure for fitting the pressure and pressurederivative plots.
Once the configuration of the reservoir heterogeneity in the vicinity of wells has been identified that information is used to condition highresolution reservoir models. The conditional probability distribution that characterises the uncertainty in permeability value at any location is perturbed using the dynamic pressure response as conditioning information. The gradual deformation of the conditional probability distribution is carried out within a pfield simulation framework. In the pfield approach, permeability values are sampled from the conditional distributions using a correlated field of random numbers (or uniform probability values). This approach retains the computational efficiency of the traditional gradual deformation algorithm, while at the same time is amenable to modelling nonGaussian permeability fields that exhibit severe discontinuities such as facies/indicator type distributions. Moreover, since the updated conditional distribution is available in the pfield approach, uncertainty assessment is possible by sampling several realizations from the updated distribution. The application of the proposed method is demonstrated on a realistic 3D example.



Measuring the Value of TimeLapse (4D) Seismic as Part of History Matching in the Schiehallion UKCS Field
More LessIn seismic history matching we use production data from wells and timelapse (4D) seismic to constrain simulation models so that they better represent reservoir properties and behaviour. Together, these data types reduce the nonuniqueness of the problem, and therefore reduce the uncertainty of both the reservoir description and also the estimation of future behaviour. The more constraints we have, however, the harder it is to find the best models and more simulations may be required to search the parameter space. This leads to increasing computing costs, which must be balanced against the reduction in model uncertainty.
We have developed a method of performing a cost:benefit analysis of including extra data and simulations in the history matching process. We use a Neighbourhood Algorithm to sample the parameter space and work in a Bayesian framework to determine model probabilities. After history matching, we then resample the posterior probability density to estimate parameter uncertainties. In addition, the parameter sampling has a density roughly in proportion to the probability distribution of the models. With this property of our method and with sufficient models, we then determine the most likely model outcome and its uncertainty. This enables calculation of expected saturation and pressure distributions at the time the data was measured and into the future. This is beneficial for reservoir management, particularly for identifying unswept areas.
We apply our method to a UKCS field and analyse how the uncertainty changes in response to adding the seismic data to the history match. We also analyse the change in uncertainty as a function of the number of simulations carried out. We identify an optimum number of models that are required before we enter the domain of diminishing returns. We confirm that seismic is important if we wish to describe the reservoir some distance from production wells. We also find that some parameters may be determined more quickly than others, depending on their location relative to the data being used.



Wavelet Based Regularization of the Well Test Deconvolution Problem
By O. SaevareidPermanent downhole measurements provide "well test" data in abundance, but their behaviour often reflects the erratic environment of everyday well operations rather than the more sterile conditions typical of a traditional well test. Consequently, processing and interpretation demand an increased level of sophistication.
In reference [1] the authors make a strong case for total least squares (TLS) applied to deconvolution of well test data. Their approach includes regularization based on penalizing the total curvature of the response function.
The present effort combines a TLS approach with regularization based on a discrete wavelet transform. As indicated in reference [2], this allows a systematic approach where the regularization have a concrete interpretation in terms of resolved details of the response function. Highly contaminated data only allows the most prominent details to be determined within relevant accuracy, while improved data quality allows correspondingly more details to be revealed.
Preliminary results demonstrating feasibility, accuracy and robustness will be presented.
[1] von Schroeter,T., Hollaender,F and Gringarten,A.C.  Deconvolution of Well Test Data as a Nonlinear Total Least Squares Problem  SPEJ 9 no 4, 375390 2004
[2] Nikolaou,M. and Vuthandam,P.  FIR Model Identification: Achieving Parsimony through Kernel Compression with Wavelets  AIChE J. 44 no 1, 141150 1998



Ensemble Kalman Filter for Field Estimation – Investigations on the Effect of the Ensemble Size
Authors G. Naevdal and K. ThulinThe ensemble Kalman filter has obtained popularity for field estimation, in particular for updating of reservoir simulation models. Here we will present a through study on the effect of the ensemble size. To be able to run a large number of simulations, we have simplified the modeling by focusing on the heat equation, but will also include a discussion on what to expect for reservoir simulation models. We find that a modest ensemble size is enough to match the measurements, but the ensemble size has to be increased significantly to get a proper picture of the model uncertainty.



Estimation of Distribution Algorithms for History Matching
Authors I. Petrovska and J. N. CarterIn previous history matching studies it has been shown that there may be multiple local optima to the response surface, even when the inverse problem is well defined in a mathematical sense. Practical algorithms that allow the identification of these local optima, such as Genetic Algorithms, have been demonstrated to work. Whilst it is useful to know where in parameter space multiple solutions exist, it is not every thing we would wish to know. At each local optimum we would like to know the range of uncertainty for each parameter and how the parameters are correlated. This will allow us to make more useful predictions, including better estimates of the uncertainty in those predictions.
In this paper we demonstrate the use of a simple Estimation of Distribution Algorithm, Probability Based Incremental Learning, on a simple reservoir cross sectional model with three parameters which is known to have multiple high quality local optima. The probability distribution function, for each parameter, can be approximated by a histogram which is adjusted using the results of the search. The sampling of the parameter space is guided by the current pdfs. We show that this algorithm can evolve steadystate pdfs which would allow us to sample the parameter space more efficiently when estimating uncertainty.
We have also introduced a modified version of the sum of square objective function. This allows a better treatment of water break through as part of the objective function. The result is that some of the optima are wider and this allows optimisers to find them more easily.



History Matching and Uncertainty Quantification Assisted by Global Optimization Techniques
Authors A. Castellini, B. Yeten, U. Singh, A. Vahedi and R. SawirisRanges in production forecast provide critical information for reservoir management decisions. Well developed methodologies exist for handling subsurface uncertainties for new field developments. The task is more challenging for fields that have been produced for several years as all models need to be conditioned to available production data in order to deliver reliable predictions. The computing cost associated with the exhaustive search of models that reproduce historical data is in general prohibitive. This paper describes an efficient method that combines the strength of various techniques, including optimization algorithms, experimental design and nonlinear response surfaces. It is applied to an offshore field in the Philippines for which multiple historymatched models are obtained in a reasonable timeframe.



Direct Inversion by detrending stochastic permeability Fields
More LessIn this paper a new method for direct inversion of production data to a permeability field is proposed. With this direct inversion methodology, history matching is physics based, honouring both the measurements and retaining the full variability in geology.
The method involves adjustment of the geological permeability realizations using the measurements. First, the realizations are split in a trend and in a residual. The trend is a smooth function satisfying the same boundary conditions as the original realization; the residual then contains the geological stochastic information and satisfies homogenous boundary conditions. The split is made both for the potential and for the stream function. The trends from a realization are in serious conflict with the measurements. Therefore we replace the wrong trends in the realization by trends calculated with the measurements. The result is a model that obeys all measurements as well as the geology, and the objective function in terms of indirect inversion is zero. Darcy’s law is explicitly included, continuity in flow is guaranteed and the boundary conditions are fulfilled. The prize that has to be paid is anisotropy in the resulting permeability field. However, in general this anisotropy is not serious in most of the model cells. We are still working on solutions to suppress this anisotropy.
A possible extension of our current method is to incorporate the water cut measurements by a multiplication factor applied directly to the stream function. Another extension is from 2D to 3D: the velocity can then replace the stream function and the permeability will be derived from the velocity vectors and the potential gradients.



Comparing Facies Realizations – Defining Metrices on Realization Space
Authors H. H. Soleng, A. R. Syversveen and O. KolbjørnsenA typical workflow for generating flow simulation grids goes through a facies modelling step. This typically involves setting up a stochastic model that is supposed to capture the important properties of the facies bodies in the reservoir volume in question and their uncertainties. This may be difficult or impossible within a particular modelling framework. Either we end up with a model too simple to be able to reproduce the characteristics of the reservoir, or the model parametres become too many and too difficult to specify. Hence users ask for methods able to reproduce the properties of a training image automatically. In any case one would like objective measures of similarity of facies realizations so that one is able to determine if a set of realizations have the properties that one wants.
Here we discuss possible components in a metric on a space of facies realizations and present an implementation of a facies realization analyzer program. The algorithm simply scans a number of realizations and computes the global volume fractions of each facies and the number of facies bodies of each type. Then it computes the surface areas, volumes, and extensions in each directions for the bodies and performs simple statistical analysis of the realizations and compare it with properties of a training image.
We present results of applying this software on facies realizations produced with variogram based methods, multipoint methods, and sequential Markov random fields.
The analyzer algorithm is fast, applicable in 2D and 3D, and the results are in excellent agreement with the subjective impression of similarity or dissimilarity obtained through visual inspection.



Facies Modelling in Fault Zones
Authors A. R. Syversveen, A. Skorstad, H. H. Soleng, P. Røe and J. TverangerTraditionally fault impact on fluid flow is included by assigning transmissibility multipliers to flow simulation grid cell faces colocated with the fault plane (Manzocchi et al. 1999). A new method, called Fault Facies modelling (Tveranger et al. 2004, 2005), captures fault impact by considering faults as deformed rock volumes rather than simple planes. Architectures and petrophysical properties of these deformed volumes (i.e. fault zones) are linked to a range of factors such as lithology, host rock petrophysical properties, tectonic regime, orientation, magnitude, and distribution of stress, as well as the burial depth at the time of faulting. By understanding these links and identifying bounding values for distributions and parameters, fault zone architectures and properties, as well as uncertainties attached to these, can be forecasted.
The fault facies approach allows 3D features such as anisotropic permeability fields, capillarity effects and tortuosity of flow paths inside the fault zone to be explicitly represented in the reservoir models. Furthermore, on the simulation grid scale, flow between cells on opposite sides of faults, as well as any uncertainty attached to this, can be estimated a priori rather than set deterministically a posteriori using history matching.
The paper compares fluid flow behaviour of conventional transmissibility multipliertype fault property models and fault facies type models through a series of simple tests. The study demonstrates that the fault facies concept is a technically feasible methodology that represents an alternative or supplement to standard industrial fault modelling methods.



MPG Simulation of Facies Thicknesses Interpreted through Sequence Stratigraphy – Application on a Carbonate Outcrop
Authors L. Dovera, J. Caers and J. BorgomanoOutcrop models can provide important information about reservoir architecture and heterogeneity. However, geological information from outcrop, because of its great variety and complexity, cannot be integrated to its fullest extent using traditional geostatistics. The conventional variogrambased modeling techniques typically fail to capture complex geological structures.
Multiplepoint geostatistics encompasses a set of an innovative modeling technique that performs the simulation starting from a training image, a 3D conceptual visual representation of how heterogeneities are believed to be distributed in the actual reservoir. The training image forms a gateway for geological expertise and interpretation to be quantitatively used in reservoir modelling. Training images can be constructed using objectbased simulation or using outcrop data. The latter will be the topic of investigation in this paper.
This paper presents an application of simpat, a multiplepoint stochastic simulation method for generating reservoir facies models by capturing the complex facies series of a carbonate outcrop and anchoring them to subsurface well data. The outcrop consists of a complex depositional sequence of mound and lobe bodies with complex spatial relationships. The idea is to generate reservoir models that reflect the observed geological complexity, yet at the same time are constrained to any available reservoir data. In order to reproduce this complex architecture in such models and to make explicit use of the sequence stratigraphic depositional information, the depositional bodies were not simulated directly but with a new approach relying on the simulation of facies thicknesses interpreted through sequence stratigraphy.
The main idea of this new approach is to separate thickness and facies information into two different yet coupled 3D training images, then to rely on the capabilities of simpat to jointly simulate these two properties. This approach can easily integrate complex information regarding bodies’ geometry and sequence stratigraphy and it has the potential to be applied to several different geological depositional system.



Direct Spatiotemporal Interpolation of Reservoir Flow Responses
Authors S. Srinivasan and S. SrinivasanThe traditional reservoir modeling workflow consists of first developing a reservoir model, performing flow simulation on that model, validating the model by performing history matching and finally using the history matched model to make predictions of future performance. In contrast, this paper presents two approaches to directly analyze the spatiotemporal variations of dynamic responses such as pressure and well flowrates and perform interpolation. Both these techniques are anchored to the data at the wells. Therefore, the resultant spatiotemporal predictions of dynamic response are history matched by construction. Interpolation or extrapolation of dynamic response to locations away from wells is possible using both the approaches. Therefore, the proposed approaches can be used to quickly determine optimal location to drill additional wells and to gauge the influence of reservoir management decisions.
In the first approach, dynamic responses such as pressure transients are treated as time series data. They are analyzed using wavelets that facilitate multiscale decomposition of pressure signals. Using a waveletlifting scheme, the transient signal is decomposed into averages and residuals. The corresponding filter coefficients defining the wavelets are treated as spatial random variables and estimated using geostatistics at locations away from wells. Pressure response is reconstructed at unsampled locations by employing the inverse wavelet transform.
In the alternate approach, direct spatiotemporal extrapolation of pressure is performed. The transient pressure data at the wells are first analyzed using correlation measures such as semivariograms. For simplicity, time is taken as another spatial dimension and variogram values corresponding to the resultant lagvectors are inferred and subsequently modeled. Spatiotemporal extrapolation is then performed to obtain the response at any location in space and at any instant in time. The robustness of both these approaches is verified on a number of case examples.



A Probabilistic Approach to Integration of Well Log, Geological Information, 3D/4D Seismic and Production Data
More LessThe ultimate goal of reservoir modeling is to obtain a model of the reservoir that is able to predict future flow performance. Achieving this challenging goal requires the model to honor all available static (well log, geological information and 3D seismic) and dynamic (4D seismic and production) data. This paper introduces a general methodology and workflow for reservoir modeling that integrates data from multiple and diverse sources, using a probabilistic approach addressing the possible inconsistency and/or redundancy between various data sources. The goal of the workflow is to model an unknown A (facies/petrophysical property) using data from different sources D1, D2, …, Dn. The workflow requires modeling the information content of each data source as a spatial distribution model termed P(ADi). Next, all P(ADi) are combined into a joint conditional P(AD1,D2,D3) from which reservoir models are drawn using sequential simulation. The procedure followed to modeling the information content of each data source as a spatial probability distribution model depends on the data source. Geological information about A is made quantitative through a training image, a 3D geological analog representation of the subsurface heterogeneity. Multiplepoint geostatistics translates training image information into local conditional probabilities of the unknown variable A. 3D seismic data is translated to a spatial probability distribution using a calibration between well logs and the 3D seismic data itself. Production and 4D seismic data are translated using and iterative procedure termed the probability perturbation method. Using Journel's tau model, a model for data redundancy, the spatial probability distributions are combined into a joint probability model used for sequential simulation. This paper shows that the reservoir model obtained using this approach honors both static and dynamic data simultaneously while explicitly accounting for data redundancy.



Conditioning Sedimentary Models to WellLog Data – An Application of Ensemble Kalman Filter
Authors A. Barrera, R. A. Rmaileh, S. Srinivasan and C. HuhGeological models consistent with the physics of sediment transport and deposition are increasingly becoming popular for the evaluation and development of water, oil and gas reservoirs. Several statistical methods have been used for the geologic modeling but they are based on very scarce information which does not consider the underlying physical principles that control the process of sediment deposition. Flowbased sedimentary models that consider depositional and physical principles for the distribution of sediments in space can be used to help generate geologically realistic models. However, these models have an important limitation, in that, it is impossible to assimilate reservoir specific information in these models. In order to overcome this limitation, an approach utilizing the Ensemble Kalman Filter (EnKF) technique for sequentially conditioning the sedimentary process model to chronostratigraphic information recorded at wells is presented in this paper. The EnKF technique allows the conditioning of sedimentary models to petrophysical and poroelastic acoustic measurements from seismic and welllogs. Starting from initial probability distributions representing the uncertainty in state variables, an iterative updating process is implemented within the EnKF framework. The ensemble of state variables is adjusted so that the updated model parameters, when used in the sediment deposition model will yield a match to the conditioning data. This method was initially tested with a 1Dimensional sedimentary deposit model for quasisteady state turbidity currents. The geologic model was then extended to a 2Dimensional space and the probabilistic method of Lattice Boltzmann Simulation was implemented using the threespeed D2Q9 lattice.
Besides presenting a novel application of EnKF for arriving at dataconditioned, sedimentary process models, the paper also sheds light on the capabilities and limitations of the EnKF approach for data conditioning.



Data Assimilation in Reservoir Management Using the Representer Method and the Ensemble Kalman Filter
Authors J. R. Rommelse, O. Kleptsova, J. D. Jansen and A. W. HeeminkThe aim of data assimilation is to improve numerical models by adding measurement information. In case of petroleum engineering, the model might be the combination of a reservoir simulator, a rockphysics model and a wave propagation package. Measurements can originate from geology, seismics, petrophysics, downhole sensors and surface facilities. In theory, one tries to maximize the likelihood of the parameters given the measurements, where the numerical model is used as a weak constraint. In practice the problem is often reduced to a least squares problem by assuming Gaussian error statistics, resulting in a variety of related data assimilation algorithms. In this paper the Ensemble Kalman Filter (EnKF) and the Representer Method (RM) are compared. For linear systems they solve the same least squares problem; for nonlinear systems, like multiphase flow in porous media, they have their own peculiarities and utilization. A variational method like the RM might get stuck in a local minimum of the squared data misfit objective function, whereas this objective does not even have a physical or probabilistic interpretation for nonlinear models or nonGaussian probability distributions. The measurement update of a filter overestimates the jump of the forecasted reservoir states towards the observed reservoir states. These errors accumulate and cannot be corrected in an iterative process, unlike what can be done in a variational method. The RM is computationally more demanding than an ENKF, especially when the number of measurements increases. Unlike the ENKF, the RM not only calculates estimates of state variables and model parameters, but it also quantifies what the isolated effect of every measurement in space and time is on the final estimate. It is therefore a promising method to quantify the “value of information” of a specific measurement.



Characterizing Data Measurement Errors with the EM Algorithm
Authors A. C. Reynolds, Y. Zhao and G. LiThe characterization of errors in measured data is important if one
wishes to condition reservoir models to diverse data sets because measurement/processing errors determine the proper relative weights of the data. In the literature, the measurement error for each data type is often estimated by some smoothing technique in the whole data domain, which often oversmoothes the data particularly around
points where the underlying true data change sharply and results in
over estimation of the measurement error. Here, we apply a modified
EM (ExpectationMaximization) algorithm to separate measured data
into groups based on the value of the measurement and the spatial
location. By applying a moving polynomial fit within each group, we generate estimates of the mean and covariance of measurement errors. The algorithm is applied to both synthetic and field time lapse seismic data as well as production data sets. For synthetic cases, the covariance functions estimated with the EM algorithm are superior to those obtained using smoothing with a moving window. The EM based procedure also appears to give more reasonable results for the field data.



Improved Modeling of 4D Seismic Response Using FlowBased Downscaling of Coarse Grid Saturations
Authors S. Castro, J. Caers and L. Durlofsky4D seismic data is used to monitor the movement of fluids in the reservoir with time and can be incorporated into the history matching process by minimizing the difference between the 4D seismic observed in the field and the 4D seismic computed from the reservoir model. Modeling 4D seismic data involves the computation of a “base” 3D seismic, which is straightforward if the reservoir has a known initial distribution of fluids, and the “time lapse” 3D seismic, which entails use of the simulated fluid distribution. The flow simulations that provide these fluid distributions are typically performed at relatively coarse scales.
It has been observed that fine details in the saturation distribution (although below seismic resolution) can impact the seismic response (Mavko and Mukerji, 1998; Sengupta and Mavko, 1998; Sengupta, 2000). Downscaling saturation outputs from the flow simulator may therefore be required to correctly model the 4D seismic response. In this paper we propose an approximate method for downscaling saturations where local fine scale flows are simulated to reconstruct the fine scale saturation using local boundary conditions determined from the global coarse scale twophase flow solution. This reconstruction does not require any global fine scale computations, guarantees flux continuity across fine scale cells in neighboring coarse blocks, and accounts for subgrid heterogeneity. Using a 2D synthetic example, we demonstrate how ignoring the fine scale effects can produce erroneous 4D seismic responses. We demonstrate that our flowbased downscaling procedure improves these results significantly.



Using Time Domain Seismic Attributes for History Matching
Authors A. Fornel, F. Roggero and B. NoetingerThe goal of reservoir characterization is to enhance the description of oil and gas reservoirs to get a reliable prediction of their dynamic behaviour using numerical simulations. In optimized management, geological models need to be updated with each new available data (well tests, production data, seismic data, ...). The value of each individual data does not lie in its isolated use, but rather in the value it adds when integrated with other data.
We present an integrated methodology for reservoir characterization, based on a nonlinear optimization loop, which uses both production (rather local information) and 4D seismic (spatial information) informations for an improved geological modelling. The main difficulty is to manage the depth/time conversion for the integration of 4D attributes and to keep it consistent with fluid flow simulation results. To do this, we propose a new methodology to match seismic attributes in time by updating the depth/time conversion in the inversion workflow.
The proposed methodology is based on a simulation workflow integrating geological modelling, upscaling, multiphase fluid flow simulation coupled with a rock physics modelling, depth/time conversion and frequency filtering. The reliability of the geological model is improved through the minimisation of a weighted, leastsquares objective function, which measures the mismatch between the simulated results and the data (both production and 4D seismic).
The huge number of data brought by the 4D seismic increases the difficulty in solving the numerical problem during optimization process. Thus, we propose an innovative approach to compute gradients of the objective function with respect to optimization parameters.
Our methodology was successfully applied for a multiscale reservoir characterization process. The 3D model is constrained by compressional/shear impedances and the twoway travel times for compressional waves at the base seismic survey.



Using the Ensemble Kalman Filter with 4D Data to Estimate Properties and Lithology of Reservoir Rocks
Authors J. A. Skjervheim, B. O. Ruud, S. I. Aanonsen, G. Evensen and T. A. JohansenImprovement of the seismic reservoir characterisation is performed by introducing an uncertainty connected to parameters in the rock physics model. Traditionally the Ensemble Kalman filter (EnKF) method has been used to estimate permeability and porosity.
In this paper we show that when seismic difference data are available also the lithology can be estimated, which is coupled to the effective bulk modulus via the rock physics model. Incorporation of inverted seismic difference data in the EnKF introduces a large amount of data in the assimilation step. Thus, to improve the results a methodology based on a combination of a global and a local analysis scheme is proposed. The global and the local analysis are used to assimilate respectively the production data and the inverted seismic difference data, where the local scheme assume that only seismic data within a certain distance from a state variable will impact the analysis in this state variable. The technique is applied to synthetic 2D and 3D reservoir models, where effects of using local versus global analysis schemes on different inverted seismic difference data, such as acoustic impedance and Poisson's ratio, are investigated. Other evaluated factors are the effects of using an incorrect seismic data error model in the analysis schemes.



A Multiscale and Metamodel SimulationBased Method for History Matching
Authors A. A. Rodriguez, H. M. Klie, S. G. Thomas and M. F. WheelerThis paper presents a novel framework for history matching using the concept of simulationbased optimization with guided search sampling, multiscale resolution and incremental metamodel (surrogate model) generation, aimed to mitigate the computational burden of largescale history matching.
The initial stage of the framework consists of a multiscale treatment of the permeability field through successive wavelet transformations. The coarsest grid which represents a highly constrained parameter space, is sampled with the aid of a derivativefree stochastic optimization algorithm that detects the most promising search regions. Due to the size of the coarse grid, thousands of simulation runs are possible at a low computational cost. Next, a sequence of intermediate metamodels is built iteratively by gradually increasing the number of sampling points in the decision space and using these temporary models to guide an incremental sampling. This incremental sampling is dictated by the use of an optimization method that finds a local optimum solution in a few iteration steps. The iterative refinement process is terminated when the metamodel solution is capable of reproducing (within a predefined tolerance) the reservoir simulator response. These metamodels are constructed using a support vector machine approach that captures the causal relations embedded in reservoir simulation by discriminating the true signal from the noise without overfitting the simulation results. Finally, the coarse grid optimal solution is used as an initial point for the next finer grid level with the use of the inverse wavelet transform. The procedure is repeated with a decreasing number of function evaluations as the grid resolution level is increased. The objective function includes well production data and sensor measurements. Numerical experiments on realistic data reveal that the proposed framework improves the history matching process, not only in terms of computing savings and the accuracy of the estimated permeability field.



Iterative Forms of the Ensemble Kalman Filter
Authors A. C. Reynolds, M. Zafari and G. LiRecently, the Ensemble Kalman filter (EnKF) has been proposed as an alternative to traditional history matching. Because of its computational efficiency, ease of implementation into a reservoir simulator and ability to provide an evaluation of uncertainty in the reservoir model and in production forecasts, EnKF appears to be a more attractive method for integrating essentially continuous streams of dynamic data to update reservoir models and characterize uncertainty than automatic history matching based on methods such as randomized maximum likelihood using the adjoint/LBFGS approach for optimization.
Although the standard theoretical underpinnings of EnKF rest on Bayesian updating with Gaussian priors, we show that the EnKF update equations can also be derived as an approximation to the GaussNewton method, which uses an ``average'' sensitivity matrix.
This suggests that for highly nonlinear, nonGaussian problems, EnKF may not provide an appropriate characterization of uncertainty and that some form of iteration is required. By viewing EnKF through the lens of optimization, instead of Monte Carlo sampling, we derive an iterative EnKF procedure for nonlinear problems. Although the iterative scheme incorporates some of the main features of EnKF, the computational efficiency of the basic EnKF method is not preserved.
We show, however, that for some problems, iteration can provide an improved characterization of uncertainty.



Decision Making Methodology in a Risk Prone Environment  Application to Production Scheme Management
Authors I. ZabalzaMezghani, M. Feraille, B. Guard and E. ManceauPetroleum field understanding and management is always associated to uncertainties, whose importance varies during all the production period. However, even if MonteCarlo technique is well known to investigate uncertainty assessment problems, in case of reservoir simulation it is not relevant anymore due to both the high number of simulations required, and the high computation time per simulation.
In recent years, probabilistic forecasting has gained popularity and has become the preferred approach when assessing the value of a project, given the uncertainty of many input variables. Reservoir understanding and production forecasting may involve as uncertain parameters both continuous uncertain parameters, which vary inside a range of possible values, and discrete parameters to model physical or geological possible scenarios, or options to be evaluated and tested to optimize the field management.
Experimental design technique has proven its efficiency to assess risk on reservoir performances related to uncertainties on continuous uncertain parameters. On the other hand, decision tree analysis is a widely validated technique, used when a problem involves subsequent decisions, to model nested scenarios and configurations as well as decision options. It is extremely useful for quantifying the impact of each scenario and configuration and to set up the value of each possible decision, finally leading to discriminate the optimal decision among all possible options.
We present here an integrated approach, which by taking the best of both experimental design techniques and decision tree analysis will allow to manage decision making in an uncertain framework, taking into account risk associated with both technical reservoir uncertain parameters and economical uncertain parameters.
This methodology has been applied to a reservoir synthetic case to highlight the need to integrate all available sources of uncertainty, both on technical and economical parameters, to avoid suboptimal decisions in terms of maximizing the economical value of the field.



Adaptive Experimental Design for NonLinear Modeling – Application to Quantification of Risk for Real Field Production
Authors I. ZabalzaMezghani, C. Scheidt, M. Feraille, B. Guard and D. CollombierOne of the most important aspect in reservoir engineering is to quantify uncertainty in reservoir behavior. Because of the large number of parameters and the physical complexity of the reservoir, fluid flow models are complex and time consuming. In order to control the cost of an uncertainty study, traditional uncertainty management is routinely performed using proxy models of the production, advocated by experimental design methodology. This problem is complex since the impact of variables in reservoir performances is often nonregular. By selecting optimally the simulation to perform, experimental design technique allows to fit polynomial models of the response but often ignores nonregularity. In this paper, we propose an original methodology to construct irregular proxy models of the fluid flow simulator. Contrary to classical experimental designs which assume a polynomial behavior of the response, we propose to built evolutive experimental design to fit gradually the potentially irregular shape of uncertainty. This methodology benefits from the advantage of experimental design, which allows to control the number of fluid flow simulations, combined with the flexibility to study nonregular behavior. We propose here an original way to increase the prior predictivity of the approximation in the nonexplored areas of the experimental domain. Based on the pilot points methodology, the research of the more predictive estimator is realized by constraining the approximation to fictitious data. These data, which are not simulated, are calibrated to ensure a better robustness and quality of the approximation. The proxy model obtained with the evolutive methodology can be considered as a good representation of the fluid flow simulator, whose evaluation is not expensive and therefore allows better risk analysis using Monte Carlo sampling. This innovative approach has been applied to model production behavior for an offshore Brazilian field, and thus to quantify the risk associated with the main reservoir uncertainties.



The Impact of Data Errors on Uncertainty Analysis
Authors G. E. Pickup, M. A. Christie and M. SambridgeSince there is much uncertainty in reservoir modelling, it makes sense to start with coarsescale models, so that a wide range of scenarios can be assessed rapidly, before focussing on fewer, more detailed models. The simplest model for reservoir analysis is the material balance equation, and this forms a good starting point for uncertainty appraisal. Although there are drawbacks with this method, such as the assumption of pressure equilibration throughout the reservoir (or compartment), there is the advantage that a minimum number of a priori assumptions are made regarding the reservoir volume and drive mechanism.
As the first stage in a topdown reservoir evaluation procedure, we have applied stochastic history matching and uncertainty analysis to a material balance problem, using a synthetic reservoir model which had aquifer influx and high rock compressibility. A truth case simulation was run and noise was added to the resulting fluid production and pressure values to generate synthetic data sets. The parameters adjusted were the volume of oil (STOIIP), the initial aquifer size and the rock compressibility. A thorough analysis of the errors was performed, including propagation of errors in the pressure data to determine their effect on the modelled production. The Neighbourhood Approximation (NA) method was used to home in on models with low misfit. Then the posterior probability distributions and their correlations were assessed using a Bayesian approach.
Results showed that the shape of the posterior probability distributions (PPDs) depended on the assumed level of the noise. In particular, they indicated that, if the amount of noise is not assessed correctly, the position of the maximum likelihood value may be estimated incorrectly.



Quantification of Uncertainty in CoarseScale Relative Permeabilities Due to SubGrid Heterogeneity
Authors H. Okano, G. E. Pickup, M. A. Christie, S. Subbey and M. SambridgeReservoir production forecasts are essentially uncertain due to the lack of data. Specifically, it is impossible to estimate detailed heterogeneity in a reservoir. In order to mitigate the ambiguity of a model, production data is incorporated into a historymatching process. However, there is insufficient data to constrain the subsurface properties all over the field.
We investigated the parameterisation of coarsescale relative permeabilities for historymatching and uncertainty quantification. Coarsescale models are often employed in history matching, because of computational cost. The results of the investigation provided us with guidelines for historymatching a coarsescale model to the observed data by adjusting relative permeabilities.
This paper addresses two issues. Firstly, because the coarsescale model inevitably misses out subgrid heterogeneity, physical dispersion is ignored in the simulation. Secondly, the smallscale heterogeneity is not explicitly known and can only be inferred by historymatching. To solve these problems, local features in the coarsescale relative permeability curves were adjusted in historymatching to capture the effect of physical dispersion and to compensate for the effect of numerical dispersion.
The success of historymatching relative permeabilities depends on the flexibility of the saturation function. We applied the flexible Bspline function as well as a conventional power or exponential function, namely the Corey or Chierici functions, respectively. We compared these parameterisations in terms of the resulting relative permeabilities during historymatching and uncertainty appraisal.
The historymatched relative permeabilities and their uncertainty envelopes were examined in comparison with the twophase upscaling results. We used a synthetic data set for which the true solution is known. The twophase upscaling was conducted using the truth model to give a reference set of coarsescale relative permeability curves. We also compared the truth production profiles with the uncertainty envelopes which were quantified in Bayesian framework. Our results highlight the fact that the parameterisation affects the width of uncertainty envelope.



Uncertainty Quantification in Producing Fields using the Neighbourhood Algorithm
Authors M. Christie, G. Nicotra, M. Rotondi and A. GodiNowadays, the majority of the world’s most important oil and gas provinces have reached a mature stage. A proper exploitation and recovery maximisation primarily rely on the ability to foresee the consequences of different reservoir management decisions and production scenarios.
Because of the inherent nonuniqueness of the history match process, using only one dataconditioned reservoir model to forecast hydrocarbon production may lead to erroneous interpretations and discrepancies with reality. A possible solution to this inverse problem and the related uncertainty quantification has been recently introduced by a novel methodology based on a stochastic search technique, called the Neighbourhood Algorithm (NA). The above methodology consists of two main steps: an optimization phase and an uncertainty assessment phase carried out in a Bayesian framework (NABayes). After sampling acceptable datafit regions of the parameter space, the posterior probabilities are calculated and a quantitative inference is performed.
In this paper the NANAB approach was firstly evaluated by means of four analytical test functions (Branin, SixHump Camel Back, Goldstein and Price, Levy) with the aim of assessing its efficiency in searching and sampling the parameters space, verifying its stability and robustness and examining the effects of the algorithm control parameters.
Two case studies are reported in order to investigate the suitability of the methodology for real fields. The hydrocarbon production, the pressure and the water cut were selected as match variables for the considered oil and gas reservoirs. The approach was promising and the results obtained were similar to those of the manually history matched model though achieved with a significant time reduction.
In addition to the increased speed of history matching, the procedure allowed a proper uncertainty quantification by means of multiple production forecasts that better quantify risk and uncertainty in reservoir performance, which is crucial for economical evaluation and decision making.



Uncertainty Assessment of Transport in Porous Media Based on a Probability Density Function Method
Authors P. Jenny, H. A. Tchelepi and D. W. MeyerA new approach for uncertainty assessment in porous media is devised. The goal is to study tracer or phase transport while assuming that the multipoint velocity statistics is known. The method depends on solving a transport equation for the joint probability density function (PDF) of velocity and concentration (or phase saturation). Similar PDF methods have been developed for turbulent flow simulations and proved extremely successful in providing joint statistics of species concentration and velocity required for turbulent combustion modeling. As in the case of turbulent flows, the joint PDF equation for uncertainty assessment of transport in porous media is defined over a high dimensional space (physical, velocity, concentration or phase, and time). This results in severe computational limitations, if a finitevolume, finitedifference or finiteelement method is employed, which is the motivation for the use of a particle method. The crucial advantage of the PDF modeling approach is its flexibility and high level of closure. For example, no modeling is required for macrodispersion (ensemble mean of the product of velocityconcentration fluctuations). Moreover, no assumptions are made regarding correlation length scales and unlike perturbationbased Statistical Moment Equation (SME) methods, the PDF method is not restricted to small input variance. Here, the PDF methodology is presented for incompressible singlephase tracer transport in heterogeneous porous media, but more general applications involving multiphase transport are discussed. A detailed description of the PDF method is given and its accuracy is demonstrated by comparison with published Monte Carlo results.



ProjectionBased Approximation Methods for the Optimal Control of Smart Oil Fields
Authors E. Gildin, H. Klie, A. Rodriguez, M. F. Wheeler and R. H. BishopComputational improvements of instrumented largescale reservoir simulation are becoming one of the main research topics in the oil industry. In particular, the problem of closedloop control is capturing a great deal of interest for reliable reservoir management. One of the main difficulties in designing controllers for largescale reservoir systems has to do with the high dimensional statespace and parameter uncertainties. Hence, lower dimensional models, linear or nonlinear, that approximate the full order system are desirable to either mitigate the cost of largescale reservoir simulation or design efficient closedloop control systems. This work aims to compare recent advances in model order reduction techniques applied to reservoir simulation. In general, the problem of reducing the order of a largescale model is known as approximation of dynamical systems. Several techniques have been developed in the linear dynamical systems framework, namely, the Balanced Truncation, Moment Matching by Krylov techniques, among others and in the nonlinear setting, namely the use of the Proper Orthogonal Decomposition (POD) and its variants. They all share a common approach: they are based on projection techniques. This work provides a comparative analysis of these techniques with particular emphasis on Krylov approaches since they are becoming one of the most active areas of research in largescale optimal control but yet, they has not been broadly reported within the reservoir community. Preliminary computational experiments reveal that these methods offer promising opportunities to design closedloop loworder controllers for the management of largescale smart fields.



Integrated, MultipleComponent Proxy Model for Optimizing Production
Authors D. Johnson, A. S. Cullick and G. ShiThis paper addresses integrated asset optimization models that are enabled by multiplecomponent, multilayer neural network proxy models. The proxy models are used to represent nonlinear production behavior for realtime production optimization.
The paper introduces a methodology for training neural networks that are robust proxy models for nonlinear physicsbased simulators. The neural network architecture, the training, and model combinations will be presented. The training of a proxy initiates with a design of experiments. If multiple simulators, e.g. reservoir, well nodal analysis, and gathering network process analysis represent different components of the physical system, then each can be executed individually. Each results in a trained proxy which are then combined mathematically into a single proxy for the integrated system. Two example applications are presented.
The first example develops a nodal analysis model, which is critical for realtime management in fields with hundreds or even thousands of wells. The model represents the system analysis for the Inflow Performance Relationship (IPR) of the reservoir and the Vertical Lift Performance (VLP) for the tubing. The example goes through a multidimension table generation for training the neural network models. An optimization solver is used to find the oil rate where the bottomhole pressures equal tubinghead pressures for a range of conditions. The complete IPR/VLP “nodal” solution is represented by the proxy. We expand this to include not just a single IPR/VLP pair of equations but a series of subsystems.
Another example builds a proxy from flow simulation of a reservoir operating with water injection and wells with downhole control valves. A multicomponent proxy is used both to optimize the valve settings proactively in response to reservoir performance and to quickly update the reservoir model history match.
Significant contributions:
1. Demonstrate that a properly trained multilayer neural network can be a robust proxy for complex, nonlinear multiple physicsbased simulations of surface, reservoir, and well performance.
2. Demonstrate use of gradient optimization of a proxy model, for proactive optimization of field and well operations in response to highfrequency data input.
3. Show how multiple neural network proxies that represent different physical components of a field operating system can be combined mathematically rather easily to form a fully integrated model.



BangBang Control in Reservoir Flooding
Various studies have shown that dynamic optimization of waterflooding using optimal control theory has a significant potential to increase Net Present Value (NPV). In these studies, gradientbased optimization methods are used, where the gradients are usually obtained with an adjoint formulation. However, the shape of the optimal injection and production settings is generally not known beforehand. The main contribution of this paper is to show that a whole variety of reservoir flooding problems can be formulated as optimal control problems that are linear in the control and that, if the only constraints are upper and lower bounds on the control, these problems will sometimes have bangbang (onoff) optimal solutions. This is supported by a waterflooding example of a 3dimensional reservoir in a fluvial depositional environment, modeled with 18.553 grid blocks. The valve settings of 8 injection and 4 production wells are optimized over the life of the reservoir, with the objective to maximize NPV. For various situations, the optimal solution is either bangbang, or a bangbang solution exists that is only slightly suboptimal. This has obvious practical implications, since bangbang solutions can be implemented with simple onoff control valves.



Production Optimization under Constraints Using Adjoint Gradients
Authors P. de Montleau, A. Cominelli, K. Neylon, D. Rowan, I. Pallister, O. Tesaker and I. NygardThe introduction of controllable downhole devices has greatly improved the ability of the reservoir engineer to implement complex well control strategies to optimize hydrocarbon recovery. The determination of these optimal control strategies, subject to limitations imposed by production and injection constraints, is an area of much active research and generally involves coupling some form of control logic to a reservoir simulator. Some of these strategies are reactive: interventions are made when conditions are met at particular wells or valves, with no account taken for the effect on the future lifetime of the reservoir. Moreover, it may be too late to prevent unwanted breakthrough when the intervention is applied. Alternative proactive strategies may be applied to the lifetime of the field and fluid flow controlled early enough to delay breakthrough. This paper presents a proactive, gradientbased method to optimize production throughout the field life. This method requires the formulation of a constrained optimization problem, where bottomhole pressure or target flow rates of wells, or flow rates of groups, represent the controllable parameters. To control a large number of wells or groups at a reasonably high frequency, efficient calculation of accurate well sensitivities (gradients) is required. Hence, the adjoint method has been implemented in a commercial reservoir simulator to compute these gradients. Once these have been calculated, the simulator can be run in optimization mode to find a locally optimal objective function (e.g., cumulative production). This optimization procedure usually involves progressively activating constraints, with each new constraint representing a significant improvement in the objective. Proper management of degrees of freedom of the parameters is essential when calculating the constrained optimization search direction. Adjoint methods have already been used for production optimization within reservoir simulation; however, an accurate analysis of optimal management of active and inactive constraints for different type of recovery processes in fieldlike cases has not been discussed to our knowledge.



Efficient Optimization of Production from Smart Wells Based on the Augmented Lagrangian Method
Authors D. C. Doublet, R. Martinsen, S. I. Aanonsen and X. C. TaiA new method for dynamic optimization of water flooding with smart wells is developed. The algorithm finds optimal injection and production well or well segment rates. In the new method, we solve a constrained optimization problem where the net present value is maximized and the reservoir flow equations are considered as constraints. The problem is formulated as finding the saddle point of the associated augmented Lagrangian functional, and solved efficiently. The method is compared with a more traditional optimalcontrol method, based on solving the adjoint system of equations. In the examples tested the new method obtains the same maximum profit as the adjoint method using approximately the same number of iterations. An advantage of the new method is that we do not solve the flow equations exactly at each iteration. As the optimization proceeds, the flow equations will be fulfilled at convergence. Thus, each iteration of the minimization algorithm is much cheaper than for the adjoint method. The method is tested on a small 2D model, but the results should be valid also for larger, 3D models.



Reservoir Flow Uncertainty Management in Presence of Stochastic Parameters
By E. FetelThis paper focuses on the management of uncertainty associated with production variables in presence of stochastic uncertain input parameters. In particular, it aims at dealing with ndimensional nonlinear response surfaces. A stochastic parameter is defined when the relationship between its variations and flow response variations is purely random. A typical example is the seed for geostatistical simulations. Alternatively, if the relationship is not random the parameter is said continuous. Here, the key idea is to model not a single response surface but a probability density function varying in the ndimensional space of the continuous parameters. In this framework, this paper develops (1) a response surface building approach, (2) a variance based sensitivity analysis scheme for identifying influential parameters and (3) a bayesian inversion technique for integrating a given production history. The proposed techniques do not require any prior regression model and are based on Monte Carlo sampling. Thus, the developed approach is suitable for ndimensional and nonlinear problems. Finally, the approach is validated on a fluviatilelike reservoir model.



Computational Techniques for ClosedLoop Reservoir Modeling with Application to a Realistic Reservoir
Authors P. Sarma, L. J. Durlofsky and K. AzizThis paper extends and applies novel computational procedures for the efficient closedloop optimal control of petroleum reservoirs under uncertainty. It addresses two important issues that were present in our earlier implementation [2] that limited the application of the procedure to practical problems.
Specifically, the previous approach encountered difficulties in handling nonlinear path constraints (constraints that must be satisfied at every timestep of the forward model) during optimization. Such constraints (e.g., maximum liquid production rate) are frequently present in practical problems. To address this issue, an approximate feasible direction optimization algorithm was proposed. The algorithm uses the objective function gradient and a combined gradient of the active constraints [3], both of which can be obtained efficiently with adjoint models. The second limitation of the implementation in [2] was the use of the standard KarhunenLoeve (KL) expansion for parameterization of the input random fields of the simulation model. This parameterization is computationally expensive and preserves only twopoint statistics of the random field. It is thus not suitable for large simulation models or for complex geological scenarios, such as channelized systems. In another paper [4], a nonlinear form of the KL expansion, referred to as kernel PCA, is applied for parameterizing arbitrary random fields. Kernel PCA successfully addresses the limitations of the KL expansion, and is differentiable, meaning that gradientbased methods can be utilized in conjunction with this parameterization within the closedloop.
An example based on a Gulf of Mexico reservoir model is considered. For this case it is demonstrated that the proposed algorithms indeed provide a viable realtime closedloop optimization framework. Application of the closedloop methodology is shown to result in a 25% increase in NPV over the base case. This is almost the same improvement achieved using an openloop approach, which is an idealized formulation in which the geological model is assumed to be known.



Stochastic Subspace Projection Methods for Efficient Multiphase Flow Uncertainty Assessment
Authors H. M. Klie, M. F. Wheeler, G. Liu and D. ZhangThis work introduces an efficient Krylov subspace strategy for the implementation of the KarhunenLoève moment equation (KLME) method. The KLME method has recently emerged as a competitive alternative for subsurface uncertainty assessment since it involves simulations at a lower resolution level than Monte Carlo simulations. Algebraically, the KLME method reduces to the solution of a sequence of linear systems with multiple righthand sides. We propose a Krylov subspace projection method to efficiently compute different stochastic orders and moments of the primary variable response from the zeroorder solution. The Krylov basis is recycled to deflate and improve the initial guess for the block and seed treatment of righthand sides. Numerical results are encouraging to extend the capabilities of the proposed stochastic framework to address more complex simulation models.



Some performance investigations of a compositional reservoir flow simulator
Authors B. O. Heimsund, E. Øian and G. E. FladmarkA reservoir fluid flow simulator for the study of advanced subsurface processes is described. It is based on the volume balance formulation for implicitly determining the fluid pressure. A set of component conservation equations are used to determine the molar masses of each component, and an energy conservation law is used to determine the system emperature. Fluid property calculations are clearly separate from flow computations, and this is illustrated by a brief discussion on a compositional blackoil formulation. As wellflow is not of primary interest, only some very simple source and sink treatment has been implemented. Different degrees of implicitness have been tested. To try to keep the formulation simple, we have preferred to use a sequential approach in which the different equations are solved in sequence. Then for each equation, and implicit or explicit solver may be used. For the pressure, only implicit solvers have been used, while for the molar masses, explicit and different types of implicit discretization have been tried. However, we were not able to develop an implicit molar mass solver which yielded higher performance than a simple explicit solver. This appeared to be due to large decoupling errors and the increased time to solve the additional implicit equations. Consequently, the simulator is using an implicit pressure solver followed by a temperature and explicit molar mass solvers. Another aspect of the simulator is its use of general unstructured grids. On those grids, it is often necessary to apply a multipoint flux approximation scheme (MPFA). As many commercial simulators experience performance degradation using MPFA rather than a twopoint flux approximation (TPFA), we provide some examples which show that in our case, MPFA is not considerably slower than TPFA on unstructured meshes. This is partly due to the choice of not using a high degree of implicitness, which reduces the size of the linear systems. Finally, the simulator parallelization strategy is outlined. A domain decomposition preconditioner is used as a solver, and the its performance is ascertained using a varying number of subdomains.



A Sequential Splitting Strategy for CO2 Storage Modelling
Authors L. Trenty, A. Michel, E. Tillier and Y. Le GalloResearch and development methodologies for the storage of CO2 in geological formation are in developing over the last 10 years. In this context, numerical simulators are the practical tools to understand the physical processes involved by acid gas injection and evaluate the long term stability of the storage.
CO2 storage models can be seen as a mix between two types of models: a reservoir model coupling multiphase flow in porous media with local phase equilibria and an hydrogeochemical model coupling transport in aqueous phase with local chemical equilibria and kinetic reaction laws.
In a recent paper, Nghiem et al. [1] proposed a fullycoupled method to solve this problem extending the local PVTunknowns elimination strategy to geochemical problems. This paper presents a different approach which combine the two models in a single computation using a sequential splitting method. The main advantage of this approach is the ability to choose a specific numerical scheme for each model.
From a mathematical point of view, this method can be seen as a two stage time integration scheme. At each stage we compute an approximate solution of the main problem using some predicted terms estimated from previous stages. This scheme is a sequential splitting strategy so it may induce some numerical errors relevant for the physical coherence of the model. To correct this error without using any iterative method, we propose to use a correction stage.
This method has been implemented in a research computer code. The code is applied to model a fieldscale CO2 storage in an heterogeneous saline aquifer.
[1] Long Nghiem, Peter Sammon, Jim Grabenstetter and Hiroshi Ohkuma, ”Modeling CO2 Storage in Aquifers with a FullyCoupled Geochemical EOS Compositional Simulator”, SPE 89474, April 004.



Unstructured Higher Resolution Convective Schemes for Flow in Porous Media
Authors S. Lamine and M. G. EdwardsNovel high resolution schemes for convective flow approximation are developed and coupled with general continuous fulltensor Darcy flux approximations. This development leads to new locally conservative formulations for multiphase flow in porous media.
The higher order schemes are designed for fluid transport on unstructured grids and are constructed using slope limiters such that they are stable with a maximum principle that ensures solutions are free of spurious oscillations. The schemes are developed to handle unstructured meshes with variable grid spacing and different formulations are compared.
Benefits of the resulting schemes are demonstrated for classical test problems in reservoir simulation including cases with full tensor permeability fields. The test cases involve a range of unstructured grids with variations in grid spacing, orientation and permeability that lead to flow fields that are poorly resolved by standard simulation methods. The new formulation is compared with standard reservoir simulation schemes and controlvolume finite element methods. The new schemes are shown to effectively reduce numerical diffusion while resolving flow induced by rapid changes in permeability, leading to superior resolution of concentration and saturation fronts compared to other schemes.



Symmetric Positive Definite Subcell CVD Schemes for CellCentred Quadrilateral Grids
Authors M. G. Edwards and M. PalA new family of locally conservative subcell fluxcontinuous schemes is presented for cellcentred approximation of the generaltensor pressure equation on quadrilateral grids. The local position of flux continuity quadrature point defines the scheme. This work continues the development of symmetric positive definite subcell ControlVolume Distributed (CVD) schemes first introduced in [1, 2], where a piecewise constant general geometrypermeability tensor approximation is introduced over each subcell of a controlvolume. Physicalspace fluxcontinuous schemes possess nonsymmetric matrices for general quadrilateral cells, or indeed any general cell type. The subcell tensor approximation ensures that a fluxcontinuous finite volume scheme is obtained with a symmetric positive definite (SPD) discretization matrix on any grid, structured or unstructured [1, 2, 3].
By definition, tensor approximation at the subcell level leads to a finer scale tensor approximation compared to the cell level, and consequently can be expected to yield a superior SPD approximation to that of earlier cellcentred cellwise constant SPD tensor schemes. A numerical convergence study confirms that the SPD subcell tensor approximation reduces solution errors when compared to the SPD cellwise constant tensor schemes. In addition, constructing the subcell tensor approximation using controlvolume face geometry yields the best results consistent with design of the schemes. Comparisons are also made with the physical space schemes. A particular quadrature point is found to give the best numerical convergence of the subcell schemes for the cases tested.



Algebraic Multigrid Based Preconditionners for Oil Reservoir Simulation
Authors P. Bonneau, Y. Achdou, R. Masson and P. QuandalleAlgebraic MultiGrid (AMG) methods have become attractive methods when used as preconditioners for Krylov subspace linear solvers.
We present here two preconditioning strategies using AMG methods for the solving of large ill conditioned systems arising from the discretization and linearization of fluid mass conservation laws.
The first strategy is a "combinative" preconditioner combining an ILU(0) step on all types of unknowns(pressure, saturation) and a more specific step on pressure unknowns with an AMG method.
The second strategy is performed in a single step with a Block Aggregation AMG. The Block Aggregation AMG method has been developed to cope with matrices araising from the discretization and the linearization of PDE systems. It is also known for its low computational cost and its simple implementation in comparison to classical AMG. This simplicity has to be paid: Aggregation AMG is less efficient than classical AMG, which led us to enhance it with an ILU(0) smoother.
One difficulty of PDE systems is the coupling existing between the different types of unknowns. A widespread idea, in reservoir simulation, proposes to perform a "local decoupling" by scaling the original matrix of the system by the inverse of the block diagonal matrix built with the diagonal blocs of the original matrix.
We will see that the scaling/decoupling of the original system is not well suited for AMG based preconditioners, resulting in bad performances of the linear solver. As a consequence, we will abandon the scaling/decoupling of the original matrix when using an AMG preconditioner.
The two AMG based preconditioners and a simple ILU(0) preconditioner will be compared on "black oil" simulations on fields with highly heterogeneous permeabilities.



Multiscale Mortar Mixed Finite Elements for Multiphase Flow in Porous Media
Authors M. F. Wheeler, G. Pencheva, S. Thomas and I. YotovWe present a multiscale mortar mixed finite element method for multiphase flow in porous media. The method is based on a domain decomposition (or coarse grid). Mass balance equations in matching fine grids of scale h, while continuity of fluxes is imposed via mortar finite elements on a coarse scale H. Higher order mortar spaces on appropriately chosen coarse interface grids are used to provide optimal fine scale convergence. For example, for the lowest order RaviartThomas mixed method or cellcentered finite differences, a choice of H = O(sqrt h) and quadratic mortars gives O(h) convergence for both the pressure and the velocity. The nonlinear algebraic system in a fully implicit discretization is solved via a non overlapping domain decomposition algorithm, which reduces the global problem to an interface problem for one pressure and one saturation. Computational experiments for oilwater displacement in highly heterogeneous media illustrate the efficiency of the multiscale mortar formulation versus the singlescale
approach.



Control Volume Discretisation on NonMatching Meshes in 3D
Authors E. Øian, B. O. Heimsund, G. T. Eigestad and I. AavatsmarkComplex geometric and geological features of realistic reservoirs motivate the need for flexible and robust flux discretisation for forecasting of fluid flow in porous media. Geomodelling tools allow faults in the stratigraphy and local grid refinements which typically result in nonmatching meshes with hanging nodes. Flow simulations based on cell centred control volume discretisations can be performed when connectivities are established and transmissibilities are computed across the nonmatching interfaces. Ideally, the interface couplings should be in a form that allows the use of robust flux discretisation techniques. Existing approaches for such problems include generating prismatic ghost cells along the nonmatching interfaces or deriving specialised stencils for possible mesh configurations.
We study an alternative technique that deals with nonmatching meshes by combining a conversion to a topologically conforming mesh with a general Omethod multi point flux approximation. Hangingnodes are removed by adding them to the mesh interfaces, creating general polyhedral mesh cells. Then the general pointbased algorithm for setting up Omethod interaction regions is used.
Since, in three dimensions, this approach can lead to a mismatch between the number of degrees of freedom and the number of continuity conditions, the mesh must be created with care. Due to the advantages from both an implementation perspective and potential good quality in the solution, it is nevertheless worthwhile to examine the properties of this approach. We present guidelines to ensure that the resulting conforming mesh is consistent with the Omethod and also how to apply a reduced flux stencil for arbitrary polyhedral meshes.



Nonlinear TwoPoint Flux Approximations for Simulating Subsurface Flows with FullTensor Anisotropy
Authors B. T. Mallison, Y. Chen and L. J. DurlofskyFulltensor anisotropy effects are often encountered in subsurface flow due to either grid nonorthogonality or permeability anisotropy. A multipoint flux approximation (MPFA) is generally needed to accurately simulate flow for such systems, though the resulting discrete system is more complex and may be less robust than that resulting from a twopoint flux approximation (TPFA). In this paper, we present and apply a different approach, nonlinear twopoint flux approximation (NTPFA), for modeling systems with fulltensor effects and grid nonorthogonality. NTPFA incorporates global flow into the determination of the twopoint flux transmissibilities, which is analogous to the use of twopoint flux transmissibility in global and localglobal upscaling procedures. For a given flow scenario, the global MPFA solution can be used to compute the twopoint flux transmissibility, leading to a global NTPFA scheme. To avoid solving the full global MPFA system, we have also developed a localglobal NTPFA procedure, in which the global flow is approximated by local MPFA solutions iteratively coupled with a global TPFA solution. The NTPFA methods described here are applied to 2D parallelogram grids, heterogeneous fulltensor permeability fields, and a 3D multiblock model with grid nonorthogonality. Results from both NTPFA schemes are generally in close agreement with the MPFA solution and provide considerable improvement over the standard TPFA scheme.



Fast Sequential Implicit Porous Media Flow Simulations Using Multiscale Finite Elements and Reordering of Cells for Solution of Nonlinear Transport Equation
Authors J. E. Aarnes, S. Krogstad, K. A. Lie and J. R. NatvigIt is demonstrated previously in the literature that multiscale methods can used to provide accurate highresolution velocity fields at a low computational cost. However, to achieve enhanced accuracy in flow simulations compared with a standard approach, the multiscale method must be accompanied by a transport solver that can account for the finescale structures of the velocity fields. In this paper, we use the standard implicit singlepoint upwind (SPU) finitevolume method for computing transport. This method requires that a nonlinear system is solved at each timestep. However, if we assume (as in streamline methods) that capillary forces can be disregarded and that gravity can be treated by operator splitting, and reordering the cells in an optimal way, the nonlinear systems in each implicit advective step can be solved on a cellbycell (or blockbyblock) basis. This approach makes the standard SPU method at least as fast as a streamline method, even on geocellular models with multimillion cells, and alleviates many limitations that streamline methods have. In particular, the method is mass conservative, compressibility can be handled in a straightforward manner, and pressure can be updated frequently without severely influencing the computational efficiency. By combining this transport solver with a multiscale pressure solver, we obtain a very efficient solution method capable of direct simulation of geocellular models with multimillion cells within an acceptable timeframe on a single desktop computer.



A Compact MPFA Method with Improved Robustness
Authors I. Aavatsmark, G. T. Eigestad and J. M. NordbottenMPFA methods were introduced to solve control volume formulations on general grids. While these methods are general in the sense that they may be applied to any grid, their convergence properties vary.
An important property for multiphase flow is the monotonicity of the numerical elliptic operator. In a recent paper, conditions for monotonicity on quadrilateral grids have been developed. These conditions indicate that MPFA formulations which lead to smaller flux stencils, are desirable for grids with high aspect ratio or severe skewness and for media with strong anisotropy or strong heterogeneity.
We introduce a new MPFA method for quadrilateral grids termed the Lmethod. The methodology is valid for general media. For homogeneous media and uniform grids, this method has fourpoint flux stencils and sevenpoint cell stencils in two dimensions. The reduced stencil appears as a consequence of adapting the method to the closest
neighboring cells.
We have tested the convergence and monotonicity properties for this method, and compared it with the Omethod. For moderate grids the convergence rates are the same, but for rough grids with large aspect ratios, the convergence of the Omethods is lost, while the Lmethod converges with a reduced convergence rate.
The Lmethod has a somewhat larger monotonicity range than the Omethods, but the dominant difference is that when monotonicity is lost, the Omethods may give large oscillations, while the oscillations with the Lmethod are small or absent.



A New Stability Criterion for the IMPSAT Formulation of Compositional Fluid Flow
Authors J. Haukås, I. Aavatsmark, M. Espedal and E. ReisoA new stability criterion for the IMPSAT formulation of compositional fluid flow has been developed. The new criterion may allow for significantly larger stable timesteps than the conventional IMPSAT stability criterion.
The IMPSAT formulation implies that interblock flow terms are evaluated with pressure and saturations from the current time level (implicitly), but with additional variables from the previous
time level (explicitly). The new criterion is associated with explicit treatment of variables that may be interpreted as complementary to volumes, referred to as isochoric variables. Analysis shows that the isochoric variables represent the part of the system that may change even though pressure and volumes are kept fixed. The interpretation leads to the identification of two mutually orthogonal subspaces of the computational space, a volume space and an isochoric space. The stability criterion is based on projections onto the isochoric space, and requires that the hyperbolic, isochoric part of the flow must not exceed the corresponding part of the mass present in a gridblock.
Whereas the conventional IMPSAT stability criterion puts restrictions on the entire flow of each component, the new criterion only limits the isochoric part of the flow. In many cases, larger timesteps are therefore allowed. The relative improvement may be in the range of 50150 %, as shown in numerical examples, and the predicted maximum stable timesteps are reasonably close to the observed maximum stable timesteps. However, for cases where there is little or no saturation change between the hydrocarbon phases, e.g., for retrograde gas condensate cases or single hydrocarbon phase cases, the improvement is insignificant.
