 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
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.



Grid Optimization for Improved Monotonicity of MPFA Solutions on Unstructured Grids
Authors M. J. Mlacnik and L. J. DurlofskyUnstructured grids are useful for resolving key geological or flow features in reservoir simulation. Accurate finite volume discretization of the reservoir flow equations on such grids generally requires the use of multipoint flux approximation (MPFA). MPFA can be applied to heterogeneous, anisotropic systems on generally unstructured grids, though it may suffer from loss of monotonicity of the inverse of the resulting linear operator at moderate to high permeability anisotropy ratios. As a consequence, the resulting pressure solution can show errors in the form of spurious oscillations. In recent work, we developed new algorithms for 2D and 3D systems that address this loss of monotonicity from a grid optimization perspective. Given an underlying finescale heterogeneous permeability field, the approaches can be generalized (via iteration) to couple permeability upscaling with the unstructured grid optimization.
In our previous work, solutions of the singlephase incompressible pressure equation were considered for 2D and 3D models. In this paper we describe our grid optimization procedures, consider new example cases, and apply the approach to a twophase flow problem. We demonstrate that the overall procedure is capable of providing accurate solutions (for both singlephase and twophase flows) that are free of spurious pressure oscillations. An upscaling example demonstrates results in reasonably close agreement with the corresponding finegrid solution. Although not presented here, the method can be coupled with a flowbased grid generation procedure to provide an overall gridding capability that resolves key flow regions while minimizing or eliminating spurious pressure oscillations.



A SemiAnalytical Model for Productivity Testing of Complex Well Configurations
Authors P. A. Fokker, G. K. Brouwer, F. Verga and D. FerreroThis paper presents a semianalytical method for the modeling of productivity testing of vertical, horizontal or multilateral wells. The method, which is applicable to both oil and gas reservoirs, automatically accounts for well interference. The use of analytical expressions ensures proper handling of transient shorttime behavior and semisteadystate longtime behavior, both close to the well and further into the reservoir. Calculation times are still very limited, in the order of a few minutes down to a few seconds when there are vertical wells only. This makes the tool suitable for well testing evaluation.
The approach is based on an earlier derived productivity prediction tool, in which the steadystate equations were solved. It has now been extended to solve the timedependent diffusion equation and it is thus more rigorous than the extension to timedependent behavior using solutions to the Laplace equation and moving pressure boundaries, which was presented recently. In our current method, the equations have first been transformed using the Laplace transformation. The expressions for the producing wells are combined with auxiliary sources outside the reservoir. The core of the semianalytic method involves an adjustment of the positions and strengths of these sources in order to approximate the boundary conditions at the reservoir boundaries. The solution that is obtained is transformed back into the time domain using a Stehfest algorithm.
The new approach has been validated with numerical tools, including both reservoir simulators and welltest interpretation software. Validations were performed with artificial cases using both singlewell and multiplewell production tests. The results of these tests were excellent.



Application of SFunction for Well Test Analysis
By T. ShimamotoOnedimensional radial flow equation, of which permeability and porosity are in inversely proportional to the radius, is equivalent to the equation of linear flow. In other words, if we take the radial variations of permeability and porosity into consideration, we can represent the linear flow by onedimensional radial flow equation. To expand this idea and to apply it for various types of inner and outer boundary shapes, we defined the Sfunction as "The functions of reservoir properties in radial direction which should be used for a radial composite model in order to represent the pressure transient behaviors at the well location for given inner and/or outer boundary shapes". We found that the approximate Sfunction can be easily estimated from complex velocity potentials for simple boundary shapes(SPE100174).
In this paper, the Sfunctional analyses are expanded to freeform boundary shapes. For this purpose, we have to derive the complex velocity potential for arbitrary boundary shapes. The solutions for these kinds of problems are known as “SchwarzChristoffel conformal mapping”. In decades, the development of the computer technologies enables us to calculate the conformal mapping easily.
The typical SchwarzChristoffel conformal mapping is a formula that will transform the complex velocity potential in a canonical domain into a polygon of the physical domain. This theory can be extended to the problem that sinks and /or sources (singular points) exist. Thus, we can obtain the complex velocity potential and at last we can calculate the pressure behavior for arbitrary shaped boundary problems. Here we will mainly discuss a “strip type transform” for the bent channel connected to fan system and about a “disk type transform” for the several hydraulic plane fractures which have different directions each other.



Deflation Accelerated Preconditioned Conjugate Gradient Method in Finite Element Methods in Oil Reservoirs
Authors F. J. Vermolen, P. L. J. Zitha and C. VuikOil reservoirs generally contain several layers with highly varying permeabilities. Hence modeling oil and water flows often involves solving partial differential equations with very large contrasts in the coefficients. Since the injection and production wells are very small compared to dimensions of the reservoir, the injection wells and production wells are modelled by the use of deltafunctions appearing in the right hand side of the partial differential equation for the pressure. In the presentation we consider a reservoir that consists of several layers with extreme contrasts of the permeability at the interfaces between the adjacent layers. Further, the finite element mesh is refined in the vicinity of the production and injection wells.
The finite element discretization of the above equation gives a stiffness matrix with extremely varying coefficients and hence the spectrum consists of large eigenvalues and eigenvalues that are almost zero, which gives a very high condition number. Hence a very bad convergence behavior for an iterative solver such as the conjugate
gradient method results. A preconditioner, like ILU, removes almost all the small eigenvalues, however, some small eigenvalues due to the large ratio of the coefficients at the interfaces persist. These small eigenvalues are removed by deflation based on a set of vectors that approximate the span of the corresponding eigenvectors.
Herewith the speed of convergence is successfully enhanced and the computational cost are reduced significantly. By the use of a proper choice for the deflation vectors, we show that the speed of convergence of our method does not depend on either the value of the contrasts in the coefficients or on the number of layers with varying coefficients. Further, the method is scalable in a parallel computing environment.



Scalability and Load Balancing Problems in Parallel Reservoir Simulation
Authors T. Guignon, J. M. Gratien, J. F. Magras, P. Q. Quandalle and O. R. RicoisScalability and Load Balancing Problems in Parallel Reservoir Simulation
New parallel reservoir simulator software designed for Linux clusters enable to overcome hardware limitation and to simulate models with large amount of data. Reservoir Engineering industry is very interested in using ever growing dataset with more and more complex physics and detailed models. The key issue still remains running simulations in an acceptable CPU time. As, the trend in hardware technologies is not to improve drastically the performance of individual CPUs but to facilitate the aggregation of computation facilities (with high bandwidth network, multicore architectures, ...), the challenge is to improve the efficiency of reservoir simulation software on a large number of processors.
New numerical difficulties and performance problems appear when the number of cells and the number of processors are growing. As a matter of fact, the architecture of Linux clusters is very sensible to memory distribution and load balancing:
 the cost of parallel solver algorithm is usually sensible to the size of the reservoir model (lack of scalability) and the consequences on CPU performance can no more be neglected;
 the domain decomposition algorithms used to distribute data between processors have a great influence on the computing load balancing between processors;
 using adaptive numerical schemes with dynamic space criteria (AIM schemes, flash algorithms based on the thermodynamical state of each cell) is a source of unbalance that cannot statically be resolved;
 simulation result storages on irregular data structures, such as unstructured grids, multilateral smart wells and perforated cells, lead to store an important amount of information during the simulation. With the variety of IO subsystems found on Linux clusters the simulator must be able to adapt its IO strategy to the underlaying IO library/file system and hardware.
In this paper, we present different approaches to overcome these kinds of problems. We discuss technical choices such like :
 advanced scalable linear solver algorithm ;
 load balancing issue with different domain decomposition strategies ;
 mesh partitioner strategy and parallel solver performance management;
 flexible IO strategy from simple file system to more complex parallel file system or database.
We have developed and benchmarked these different solutions on published reference large scale problems and actual case studies with several tens millions of cells. We analyze the results and discuss the efficiency of each solution to overcome the scalability difficulties and performance limitations due to load unbalance.



NonEqulibrium TwoVelocity Effects in GasCondensate Flow through Porous Media
Authors I. Panfilova, M. Panfilov and S. OladyshkinThe flow of gascondensate mixture in porous media is characterized by three properties which determine a significant nonequilibrium in transfers between gas and liquid: 1) a capillary and gravity coagulation of small liquid drops with forming large aggregates, 2) a high difference in diffusion coefficients for liquid and gas, which leads to a delay in establishing the liquid aggregate composition, and 3) a high difference between the liquid and gas mobility, which causes the dependence of the nonequilibrium parameters on the relative phase velocity. The objective of this study is to construct a closed nonequilibrium compositional flow model and to reveal various regimes of the nonequilibrium behaviour.
The analysis was based on separating in time the capillarygravity coagulation and the phase exchange. The coagulation was studied using the method of pore network modelling. An isolated liquid drop limited by two meniscus inside a pore can move in direction of the resulting capillary or gravity force. The drop motion is limited by percolation conditions saying that the gas from the neighbouring pores can not be displaced if it has no cluster connection to the medium exit. As the result we obtained the relations between the correlation length of the pore radii field or the permeability field and the scale of the liquid aggregate.
At the next step we constructed a set of the averaged compositional flow models taking into account the phase nonequilibrium. The averaging was performed by twoscale asymptotic expansion method. Three regimes of the nonequilibrium were revealed. At a small relative phase velocity the macroscale exchange is diffusion limited, being described by a nonlocal integrodifferential operator. Its kernel is calculated as the result of solution to a cell problem. A growth of the relative phase velocity determines an increasing influence of the rotational flow inside liquid aggregates on mass transfer between the phases. Such a diffusionrotation regime is described in terms of the double relaxation model, presenting a secondorder nonlinear kinetic differential equation. The diffusion and rotation relaxation times are calculated as the result of solution to a problem of gas flow around a macroscale liquid aggregate in a porous medium with mass transfer. The fluid flow was described by the Brinkman equations which allow the rotations, in contrast to the Darcy law. The “slip regime” arises at very high relative phase velocities, when the nonequilibrium degree becomes to decrease, however the system tends to a new equilibrium state in which the time of contact between liquid and gas appears to be much lower than the time of mass exchanges. The generalized kinetic model of the slip regime is obtained by homogenization technique. A number of examples are simulated for gascondensate flow in the vicinity of a well, where the nonequilibrium degree appears to be the most significant.
The influence of the capillary number, the Forchheimer effect and the velocitydependent relative permeability on the nonequilibrium was analysed.
We develop generalized relations which describes uniformly all the three nonequilibrium regimes, which may be used as a plugin to the existing PVT or hybrid thermodynamichydrodynamic software, with providing a new option to simulate the nonequilibrium behaviour of multicomponent gascondensate mixtures.
The research is financed by the Schlumberger Moscow Research Center.



Splitting the Thermodynamlics and Hydrodynamics in Compositional GasLiquid Flow through Porous Reservoirs
Authors S. Oladyshkin and M. PanfilovFor twophase compositional flow with mass transfer in porous media, it is shown that the existing steadystate solutions are unstable for gasliquid systems, which are characterized by a high relative phase mobility. Instead of them the process model allows stable “semistationary” solutions which correspond to a limit model when the liquidgas mobility ratio tends to zero. In a semistationary model the pressure and phase concentration fields are quasistationary, while the liquid saturation, in conditions of a continuous condensation/evaporation, is strongly nonstationary and has no a stationary limit in time. Within the framework of the semistationary model we have obtained a full splitting between hydrodynamics and thermodynamics, thanks to the fact that N2 equations describing the phase concentration transport allow an explicit first integration along streamlines. This enables to transform them into differential equations of thermodynamic type (out of space and time) with differentiating over pressure. Being added to the usual equilibrium and EOS equations, such a system constitutes a new closed thermodynamic model which depends on the pressure only. This model describes the thermodynamic behaviour in an open system which corresponds to gasliquid flow along each streamline. The remaining two flow equations determine the pressure and the saturation fields, with coefficients dependent on thermodynamics. To determine them it is sufficient to calculate the thermodynamic system one time, in contrast to the full compositional model where the thermodynamics is calculated at each space and time point.
The obtained hydrodynamic model allows obtaining the analytical solutions along streamlines, by using the singular perturbation method with respect to the relative phase mobility parameter. The irregular character of the asymptotic expansions is determined by the fact that the formal limit solution corresponds to an immobile liquid, which means a full pore plugging by liquid and breaking of gas movement. The full analytical solutions to the non selfsimilar multicomponent flow problems towards a well were constructed with using the matching asymptotic expansion method. This result was used to propose a new method of gascondensate well representation in reservoir simulations. The comparison of the solutions obtained for the split thermodynamic and hydrodynamic models with numerical ECLIPSEbased fullcompositional simulations has shown a very good agreement. The suggested method was used to develop a new kind of the streamline compositional simulator. We illustrate some examples of its functioning.
The research was financed by the Schlumberger Abingdon Technology Center.



The distillation mechanism in steam displacement of oil
Authors D. Marchesin and J. BruiningDuring steam drive recovery of oil a thin zone of high volatile oil concentration is formed between the upstream hot steam region and the downstream cold liquid region containing water and the oil at its initial composition. Distillable components, coinjected or stripped from the oil remaining in the upstream steam region, are transported and build up this region. The existence of this zone is responsible for the high recovery of oil from the steam zone. This paper investigates the growth and the stability of this zone, for the case when relatively small amounts of volatile oil are present. We show that this zone is described by ODE's typical of traveling waves. The growth and stability are illustrated numerically. Our results suggest that steam drive recovery can be improved by the preinjection of a small slug of volatile oil.



Efficient Integration of Stiff Kinetics in Reactive Compositional and Thermal Porous Media Processes
Authors M. R. Kristensen, M. Gerritsen, P. G. Thomsen, M. L. Michelsen and E. H. StenbyWe propose the use of implicit onestep ESDIRK (Explicit Singly Diagonal Implicit RungeKutta) methods for integration of the stiff kinetics in reactive, compositional and thermal processes that are solved using operatorsplitting type approaches. To facilitate the algorithmic development we construct a virtual kinetic cell model. The model serves both as a tool for the development and testing of tailored solvers and as a testbed for studying the interactions between chemical kinetics and phase behavior. As case study, two chemical kinetics models with 6 and 14 components, respectively, are implemented for insitu combustion, a thermal oil recovery process. Through benchmark studies using the 14 component reaction model the new ESDIRK solvers
are shown to improve computational speed when compared to the widely used multistep BDF methods DASSL and LSODE.
Phase changes are known to cause convergence problems for the integration method. We propose an algorithm for detection and location of phase changes based on discrete event system theory. Experiments show that the algorithm improves the robustness of the integration process near phase oundaries by lowering the number convergence and error test failures by more than 50% compared to direct integration without the new algorithm.



Subsidence Simulations for Depleting Reservoirs – Considerations on the Role of Hydromechanical Nonlinearities
Authors G. Musso, D. Selvaggi and S. ManticaThe paper aims at investigating some issues related to the theoretical and numerical modelization of subsidence phenomena induced by hydrocarbon withdrawal. First a brief review of the implications of facing the problem as coupled / uncoupled is provided accordingly to the theory of mechanics of porous media as formulated by Biot. On the basis of these considerations, the classical nucleus of strain solution given by Geertsma is coupled to the hydrodinamic problem of a single phase fluid draining towards a well. The approach implemented has the advantage of providing a series of subsidence profiles evolving with time, together with the well production. The problem is handled both introducing or not the possibility of hydromechanical nonlinearities in the behaviour of the reservoir rock. When accounted, the mechanical nonlinearity is accomplished by assuming that the rock deforms accordingly to classical logarithmic relationships for oedometric stress paths, as foreseen by critical state models such as Cam Clay. On its turn the hydraulic nonlinearity is accounted for by implementing a literature relationship valid for Mexico Gulf sands, postulating that an exponential relationship exists between porosity and permeability.



A Predictor for Accelerated Coupled Rock Mechanics and Reservoir Simulation
By O. PettersenThe impact of the stress field on reservoir fluid flow and production can be significant for many kinds of reservoirs, and hence coupled Rock Mechanics and Reservoir Simulation has been seeing a growing popularity. A much used scheme is iterative coupling, where compaction is computed at each stress step by iteratively updating cell pore volumes in the reservoir simulator by values calculated from strain in the stress simulator.
Although the procedure works satisfactory, it may be slow, as often many iterations are needed. Further, the pore volume corrections will only be performed at selected stress time steps, such that pressure and compaction in the flow simulator are not continuous in time. Many reported schemes assume specific poroelastoplastic models, as e.g. linear elastic, and also require modification of code.
It is well known that compaction is a function of strain, while reservoir simulators use fluid pressure, the only compaction energy available. On this background few if any coupled procedures utilize the compaction vs. fluid relationship at all.
In this paper we show that the relationship can nevertheless be used as basis for constructing a predictor for the actual stress / strain computations, which leads to significant speedup. Many of the features of the predictor can be determined from the first stress time step only, and for later stress steps it can be improved with small effort. The scheme is valid irrespective of the poroelastoplastic model, and is based on information exchange, so no simulator code modification is necessary.
The compaction state is primarily dependent on the materials, boundary conditions, and the production process, with the geometry dependency as the governing. The predictor is constructed by modifying compaction vs. fluid pressure to take account of geometry variation. A good predictor will result in an improved pressure field as computed by the reservoir simulator, hence providing the stress simulator with a better pseudoinitialiser, such that it converges quicker, and in the pore volume iteration scheme fewer if any iterations are required.
In total we have experienced a reduction in total computer time of more than 90% in some cases, and as a bonus the fluid pressure field is continuous in time.



Modelling Induced Fracture Mechanics in Reservoir Simulations – A New Coupled Simulator
Authors D. Zwarts, B. Hustedt and P. J. van den HoekWaterinjection induced fractures are key factors influencing successful waterflooding projects. Controlling the dynamic fracture growth can lead to largely improved water management strategies and potentially to increased oil recovery and reduced operational costs (reduction in well count and water treatment facilities etc.).
The primary tool that a reservoir engineer requires to design an optimal waterflood is an appropriate reservoir simulator that is capable of handling the dynamic fracturing process in complex reservoir simulation grids.
We have developed a new modelling strategy that adds fracturegrowth to our standard fluidflow reservoir simulator. A first prototype simulator was successfully tested and applied to field cases.
The dynamic fractures are free to propagate asymmetrically in length, height, and widthdirection controlled by the pressures and poro and thermoelastic stresses acting on the fracture face. The stresses are calculated in the reservoir simulation grid. The simulator determines new fracture sizes by evaluating fracture propagation criteria, based on a Barenblatt condition on the four fracture tips: left, right, up, and down, and on a volumebalance on the width. This nonlinear fivedimensional coupled set of equations is solved every timestep using a Broyden approach moderated with LevenbergMarquardt techniques to handle nonlinearities. Since the changing fracture sizes affect the pressure and stress profiles in the reservoir, the equations are solved in an implicit scheme.



A New Model for Solution Gas Drive in Heavy Oils
Authors S. Zaleski, M. Chraibi and F. FrancoWe introduce a new Darcy scale model in which the permeability of the gas phase is replaced by a new expression for the gas phase velocity. When the bubbles are smaller that a critical radius, capillary forces anchor the bubbles and the gas phase does not move. Above the critical radius the bubbles move at a velocity of the same order as the liquid phase velocity.
The number of bubbles per unit volume may vary through coalescence between neighboring bubbles. The coalescence rate is an important parameter in the model. When coalescence is fast, the bubbles reach rapidly a large size at which they become mobile. The dependence of the coalescence rate on the fluid parameters and the bubble size is discussed.
The predictions of the model are then compared to several laboratory scale experiments.



Methodology of Fast Transmissibility Determination for MediumScale Fractures Systems
By S. VitelProductivity of fractured reservoirs is mainly affected by the fracture network configuration, which impact on the flow paths cannot be neglected. Today, two main approaches, opposite and complementary, are used to perform flow simulation on fractured reservoirs: (1) the continuum approximations, simple and efficient but approximate in the actual fracture geometry integration, and (2) the discrete model formulations, quite complex to put into practice and much slower, but respectful of the discrete fracture network and thus more accurate. Yet, fractured reservoir simulation still needs to balance the speed up of the calculation processes while keeping the fracture system paths accurate.
This paper presents a methodology to evaluate equivalent transmissibilities in naturally fractured reservoirs, based on a discrete fracture network. It is divided in two parts: (1) first the codiscretization of the fractures and the simulation grid, (2) then the determination of the equivalent transmissibility between two gridblocks of the simulation grid. The first phase consists in representing the reservoir with a connectivity list. Nodes represent control volumes, holding porosity, volume, pressure and saturation properties; pipes stand for connections between those control volumes, holding transmissibility property. First this list is extracted from the simulation grid and from the fractures that have been clipped by the gridblocks. Then both lists are connected together between each fracture piece and the surrounding gridblock. The second phase consists in computing the interblock transmissibilities from this connectivity list. These transmissibilities are evaluated for each couple of gridblocks in each direction by applying electrical simplification theorems.
The discretization tool and the simplification technique make the method really fast. Its performance is demonstrated on a 3D highly fractured reservoir. The main limitation is due to the electricity analogy limitations that does not take into account complex flow mechanisms related to fractures; even though, results show good agreement with those of a fine grid.



TwoPhase Flow in Fractured Media – Homogenized Model with Mixing and Upscaling by a StreamConfiguration Method
Authors S. Skachkov and M. PanfilovThe paper deals with the displacement process of two fluids, immiscible at the molecular scale, through a porous medium. The macroscale mixing between the phases is caused by the medium heterogeneity determined by a fracture network whose scale is much greater than the pore size, but much smaller than the reservoir size. At the heterogeneity scale the flow is described by a twophase flow model with a generalized Darcy law, classical phase permeabilities and capillary pressure. The macroscale flow model is obtained by the twoscale asymptotic homogenization method. To capture the effect of dynamic mixing, the first order model is derived.
The mixing is described by the dispersion effect and the convective velocity renormalization. The dispersion tensor, the effective phase permeabilities and the velocity renormalization are defined through the cell problems as the functions of saturation, viscosity ratio and flow velocity.
To solve the cell problems we have developed the twophase version of the streamconfiguration method proposed earlier for singlephase flow. The method is based on the following features of a system of fine fractures: a) the limit flow is locally onedimensional; b) the saturation and heterogeneity can be factorized, so the cell problem becomes independent of saturation. At the same time, the limit solution is shown to be nonunique due to a loss of information about the stream configuration geometry in the nodes of fracture intersection. The regularization procedure is developed which proposes additional conditions describing the type of stream configuration in each node. The types of configuration are determined from the principle of local energy minimum or entropy maximum. For a periodic regular fracture network, the method allows obtaining analytical solution for the dispersion tensor. For disordered networks, the method reduces the cell problem to an algebraic system of a rank equal to the number of nodes in a unite cell. A version of method is developed, based on a combination between the streamconfiguration and the bordering techniques.
A singular dispersion regime is revealed in which the dispersion tensor becomes unbounded due to arising of stagnant zones in several fracture segments and fluid trapping. The data on comparison of our results with those obtained using other methods are illustrated.
The industrial application of this research consists of a new numerical algorithm of upscaling twophase flow in fractured media, which is much faster than other methods.
The research is supported by the Schlumberger Technology Center in Abingdon.



Mathematical Model for ThreePhase Compressible Dual Porosity Streamline Simulation
Authors A. Myasnikov, A. Kozlova, F. Bratvedt and K. BratvedtFlow simulation of fractured reservoirs usually is performed using a dual porosity model. The dual porosity system is modeled by using two coupled grids: one for matrix and one for fracture. These two continua communicate by transfer functions. Until now, there were no mathematical models of dual porosity, threephase, compressible flow for streamline simulators. To realize this model, it was necessary to reformulate the matrix and fracture pressure equations The conventional transfer function has been incorporated as a source/sink term, not only in the streamline saturation equations (as it was in incompressible case), but also in the pressure equation.
The dual porosity model has been implemented into a streamline simulator. This tool has its main application in the geological modeling domain for analyzing uncertainty, model ranking and screening of geologically detailed models, including fractures.
This paper describes the mathematical model for a threephase compressible dual porosity streamline simulator and compares the results and run times of the streamlinebased approach with a conventional dual porosity gridbased commercial simulator. The results from the streamline simulator for dual porosity show good agreement with those produced by a commercial finite difference simulator with order of magnitude improvement in simulation time.
Streamline methods as a reservoir simulation tool have generated much interest in petroleum engineering because of the capability to calculate fluid flow in multimillion cell geological models with reasonable CPU times. However, important physical properties of geoscale fluid flow models are still not properly modeled by streamline methods. Enhancing the range of physical properties that can be simulated accurately in a timely manner will enable improved workflows in the geological modeling domain.



Rigorous Derivation of the KazemiGilmanElsharkawy Generalized Dual Porosity Shape Factor
Authors Z. E. Heinemann and G. M. MittermeirKazemi et al.[7] suggested to use an empirical matrix/fracture transfer function, verified based on experimental data of Mattax and Kyte[9]. Kazemi et al. showed that for rectangles and cylinders the formula reduces to the well known forms of the shape factor. In the mean time many authors indicated the validity of the formula, but no theoretical proof was offered so far.
This paper derives the KazemiGilmanElsharkawy (KSE) shape factor using Control Volume Finite Difference discretization on the fracturematrix dual continuum.
It will be shown that the KSE shape factor is valid for multiphase fow of compressible fluids irrespective of matrix block shape. Additionally an extension to tensorial matrix permeability is given.



A Learning Computational Engine for History Matching
Authors R. E. Banchs, H. Klie, A. Rodríguez, S. G. Thomas and M. F. WheelerThe main objective of the present work is to propose and evaluate a learning computational engine for history matching, which is based on a hybrid multilevel search methodology. According to this methodology, the parameter space is globally explored and sampled by the simultaneous perturbation stochastic approximation (SPSA) algorithm at a given resolution level. This estimation is followed by further analysis by using a neural learning engine for evaluating the sensitiveness of the objective function with respect to variations of each individual model parameter in the vicinity of the promising optimal solution explored by the SPSA algorithm.
The proposed methodology is used to numerically determine how additional sources of information may aid in reducing the illposedness associated with permeability estimation via conventional history matching procedures. The additional sources of information considered in this work are related to pressures, concentrations and fluid velocities at given locations in a reliable fashion, which in practical scenarios might be estimated from high resolution seismic surveys, or directly obtained as in situ measurements provided by sensors. This additional information is incorporated, along with production data, into a multiobjective function that is mismatched between the observed and the predicted data. The preliminary results presented in this work shed light on future research avenues for optimizing the use of additional sources of information such as seismic or sensor data in history matching procedures.



Hybrid Mesh Generation for Local Grid Refinement in Reservoir Flow Simulation
Authors C. Bennis, H. Borouchaki and N. FlandrinMesh generation becomes a crucial step in reservoir flow simulation of new generation. The mesh must faithfully represent the architecture of the reservoir and its heterogeneity. For more accuracy, the mesh must also follow the flow directions in the well vicinity. The current industrial standard meshes based on Corner Point Geometry (CPG) grids have already shown their limits. They are very practical and easy to use, but they fail to represent complex objects due to their structured aspect. More recently, other approaches have been proposed, in particular using the PErpendicular BIssector (PEBI) grids, which are completely unstructured. These grids are very flexible and can model most complex shapes. But, they are often difficult to manage in 3D due to their lack of structure.
A three dimensional hybrid mesh model was proposed, in ECMOR 9, to capture the radial characteristics of the flows around the wells. In this hybrid mesh, the reservoir is described by a non uniform Cartesian structured mesh and the drainage areas around the wells are represented by structured radial circular meshes. Unstructured polyhedral meshes are used to connect these two kinds of structured grids. The construction of these transition meshes is based on 3D Power Diagrams to ensure finite volume properties : mesh conformity, dual orthogonality and cell convexity. In this paper, we propose an extension of this hybrid model to the case of real CPG reservoir grids. At first, the CPG grid is mapped, in a reference space, into a non uniform Cartesian grid by minimising the mapping deformation. Then, a hybrid mesh is generated in this reference space using the previous method. Finally, this mesh is mapped back into the real space. This mesh can also be optimised with respect to the deformation metric between the real space and the reference one.



Homogenization for Fe²+ deposition near drink water tube wells during arsenic remediation
Authors J. Bruining and M. I. M. DarwishHomogenization is a well established methodology that can be used for upscaling. The advantage of homogenization over other upscaling methods is that it only requires the choice of a periodic unit cell to overcome the closure problem. The choice of the unit cell is completely flexible and only limited by computer capacity. Our interest is in the upscaling of reactive diffusion convection flows in connection with arsenic contaminated water remediation. We will rederive the homogenized equations with emphasis on the physical aspects. By considering a simple twodimensional example we cover all computational aspects linked to the actual application of the method. However, the method can be easily extended to more complicated configurations in 3D. The computed dispersion coefficient shows a similar (qualitative) dependency on the Peclet number as literature data, albeit that the values for this 2D example are considerably smaller. Moreover in the case of nonlinear adsorption the dispersion coefficient depends on the concentration.



Optimization of the Oil Field Putting into Operation Stage
Authors S. A. Ermolaev and A. I. ErmolaevThe strategy of field putting into operation is to solve the following problem: “Which wells are to be put into operation at the given moment?” There are different strategies which differ greatly by the volumes of oil recovery for the period of the field putting into operation. Consequently there is a problem how to form and choose the rational strategy. By rational strategy is meant the best strategy for given criteria. In the present paper we propose some algorithms to solve this problem.
There are two aspects while solving the problem of creating and choosing the rational strategy. The first aspect is to fulfill detailed hydrodynamic simulation of the field development. The second aspect is to enumerate a great number of possible strategies. At present either nobody realizes detailed simulation or very few possible strategies are analyzed. Thus an irrational strategy may be selected.
The suggested method of creating and choosing rational strategy of putting into operation is a synthesis of simulation and optimization algorithms. This method allows one to take into account main reasons of strategies different efficiency (reservoir heterogeneity, nonuniform watering of field sections). Also we realize purposeful and shortcut enumeration of all possible strategies of field putting into operation.
The main stages of this method are:
1. Wells grouping into blocks (number of wells in blocks is to be approximately equal, not only neighboring wells can be included in one block).
2. Computation of potential volumes of oil recovery produced by wells blocks while the moments of their putting into operation are different (these parameters are input data for optimization; these computations use hydrodynamic simulation software; also the suggested technology allows one to estimate the injection wells deposit into total oil recovery).
3. The optimum moments of putting blocks into operation or their optimum order are found (these problems are formulated as transport models; standard linear programming algorithms may be used to solve).
In this paper we consider examples using suggested method for a real deposit.
The method can be expediently used to choose rational strategies of field putting into operation when it is impossible to follow specialists experience and intuition.



Integrating Dynamic Data into Reservoir Models – A Multiple Point Perspective
Authors S. Srinivasan and K. EskandariConstraining reservoir models such as permeability distribution to dynamic well data has been traditionally accomplished using inverse theory. The resulting illconditioned system of equations needs to be solved for the required model parameters. Gradientbased optimization schemes require efficient computation of sensitivity coefficients (Chu, Reynolds and Oliver, 1992). Markov chain Monte Carlo (Omre and Tjelmeland, 1996), frequently used in earth sciences, is another iterative approach in which a prior distribution of permeability is perturbed into a posterior distribution using a probabilistic perturbation acceptance criteria. However, that technique like the pervious methods is expensive in term of CPU usage.
This paper investigates an approach to alleviate the computational expense associated with the iterative conditioning of reservoir models to dynamic information. A multiple point proxy is proposed that accounts for the configuration and orientation of geological patterns in the reservoir and their relationship to flow. The nonlinear relationship between the multiple point connectivity and the flow response is established by calibration. Once calibrated, the proxy expression acts as a surrogate to the full physics based flow simulators and can be used within an iterative framework such as Markov chain Monte Carlo algorithm to build reservoir models conditioned to well test responses. The mathematical formulation and calibration of such a proxy function for matching well pressure characteristics is presented in this paper.
Building reservoir models based on just static data using multiple point statistics instead of just two point statistics has received a lot of attention recently. The SNESIM algorithm (Strebelle, 2002) is such a multiple point (mp) simulation algorithm. This paper presents a multipoint simulation algorithm that is distinctly different from all other mpbased simulation algorithms available in the literature. The presented approach is extended to integrate dynamic data using an efficient probability perturbation based approach. The proposed method is computationally fast for retrieving mp statistics (scanning) and simulation using a unique pattern growthbased methodology. Dynamic data integration is achieved by gradually perturbing the multiple point probability distribution using a deformation parameter. The perturbed multiple point probability distribution is then merged with the prior multiple point distribution inferred from a training image. The permanence ratio of hypothesis is used to perform the merge. The resultant reservoir model is thus consistent with both the available production data as well as the prior geology.



A Numerical Solution to Groundwater Flow in Full Tensorial Permeability Fields
Authors M. Le RavalecDupin and L. RicardThe description of flow through heterogeneous porous media is of primary interest to many fields such as reservoir engineering or hydrology. Because of porosity and permeability heterogeneity, flow equations are usually solved on the basis of numerical techniques. The discretization technique most commonly used is the finitedifference method. Thus, a conventional discretization scheme involves a 5 point stencil for twodimensional media and a 7 point stencil for threedimensional media. Continuity of flux and pressure is readily incorporated into the discretization by deriving the coefficients at the interface between adjacent gridblocks from the harmonic average of the permeabilities of the two contiguous gridblocks. This approximation holds as far as permeability tensor is diagonal.
Today, very detailed geological models are built on grids containing millions of gridblocks. Although computers are growing ever more powerful, fluid flow simulations on such grids are tremendously CPUtime consuming. To make them tractable, the number of gridblocks must be reduced to about 105. This motivated the development of techniques for upscaling the geological model to a manageable level of detail. Upscaling consists in determining the equivalent permeabilities of the coarse gridblocks. These equivalent permeabilities are usually full tensors. Discretization techniques in which interface fluxes are determined from more than two points have been developed to solve full tensorial flow equations. In this case, the discretization stencil requires 9 points for two dimensional media and even more for threedimensional media. These techniques, known as multipoint flux approximations (MPFA) have recently enjoyed increased popularity. To date, to our knowledge, there is no straightforward generalization of the harmonic average to estimate fluxes at interfaces with MPFAs.
We develop an alternative to MPFAs based upon numerical spectral methods so that the full tensorial pressure equation is solved on regular grids. We also consider source or sink terms as well as gravity. Contrary to MPFAs, the proposed spectral technique does not call for any kind of approximation to estimate fluxes at the interface between adjacent gridblocks. Specific applications are presented and discussed.



Nonlinear Analysis of Fluid Flow Behaviour and Flow Control Factors in Fractured Systems Using Support Vector Machine
Authors J. Ma and G. CouplesPredicting the fluid flow behaviour of systems containing fractures or other flow conduits (which are typically subseismic in scale) is an important element of flow modelling in petroleum exploration and production. Since the fluid flow can be strongly influenced by multiple flow control factors, including the connectivity of the fracture network, the spatial distribution of fractures (e.g. fracture intensity), and the contrast of flow properties between fractures and the matrix, the fluid flow behaviour of a fractured system, with respect to flow control factors, may show multiple distinctive modes. Clearly the existence of different modes indicates a need for a multimodal approach to incorporating effective flow properties in the coarserscale flow simulation. However, the intrinsic nonlinearity of such relationships and the high dimensionality of the factors make it difficult to identify such relationships, let alone to distinguish the distinctive modes among them.
In this work, a nonlinear analysis, based on a machine learning method  Support Vector Machine (SVM), was considered for identifying the relationships. It was applied to two simple fractured systems in 2D, each of which was characterised by a set of distributional parameters and a stochastic fracture modelling procedure. For each system, equallyplausible fracture models were generated, and the singlephase steadystate fluid flow was simulated for each model. The relationships between the simulated fluid flow, in a form of upscaled permeability, and a number of flow control factors, were then analysed. The results showed the existence of complex and multimodal relationships between the upscaled permeability and the control factors in each system, with distinctive features between the systems. The implications of not employing a multimodal approach in a coarserscale simulation are obvious: the upscaled flow properties can be significantly misestimated.



Determining the Filtration and Permeability Reduction Functions for Flow of Water with Particles in Porous Media
Authors D. Marchesin, A. C. Alvarez, G. Hime and P. G. BedrikovetskyDeep bed filtration of particle suspensions in porous media occurs during water injection into oil reservoirs, drilling fluid invasion of reservoir production zones, fines migration in oil fields, bacteria, viruses or contaminant transport in groundwater, industrial filtering, etc. The basic features of the process are particle capture by the porous medium and consequent permeability reduction.
Models for deep bed filtration contain two coefficients that represent rock and fluid properties: the filtration function, which is the fraction of captured particles per unit of particle path length, and formation damage function, which is the ratio between reduced and initial permeabilities.
The coefficients cannot be measured directly in the laboratory or in the field; therefore, they must be calculated indirectly by solving inverse problems. The practical petroleum and environmental engineering purpose is to predict injectivity loss and particle penetration depth around wells. Reliable prediction requires precise knowledge of these two coefficients.
In this work we determine these coefficients from pressure drop and effluent concentration histories, measured in onedimensional laboratory experiments.
The filtration function is recovered by optimizing a nonlinear functional with box constraints. The permeability reduction is recovered likewise, taking into account the filtration function already found.
The recovery method consists of optimizing Tikhonov's functionals in appropriate subdomains.
In both cases, the functionals are derived from least square formulations of the deviation between experimental data and quantities predicted by the model.



Numerical Approach to Optimising Oil and Gas Production in Thin Oil Rim Reservoirs
Authors G. G. Cassarà, C. Monico, D. Castano and S. DresdaThe development of a thin oil rim can involve complications in the oil production due to severe gas and/or water coning events.
This paper deals with a Nigerian oil and gas multilevel field where some of the levels bear a thin oil rim associated with a large gas cap. After a few years of production from three oil levels, it was decided to develop all gas and oil levels, to increase the oil production rate and satisfy the increasing demand for gas. Therefore, the thin oil rim reservoirs were considered as potential sources of oil and gas reserves.
A 3D black oil simulation model was built to study the development of the whole field. Gas PVT data were characterised through the vaporised oil gas ratio to simulate the effects of depletion on the condensate stream during production from gas levels.
Sensitivity runs were performed to investigate the development of the oil rim reservoirs. The model indicated that the wells located in the oil rim were affected by high water cut, whilst the wells in the gas cap area and far away from the oil rim, could produce large gas rates without draining the oil in the rim.
An intermediate location was considered in the gas cap flank close to the oil rim. The wells with this location showed early high gas production with condensate decline, followed by liquid hydrocarbon production increase with gas rate reduction. During the early time, the condensate production fell due to depletion. Afterwards, the liquid hydrocarbon rate started to increase since wells drained the oil rim through coning.
The results of the simulation encouraged the development of the thin oil rim reservoirs mainly through these intermediate locations to optimise both oil and gas production.



New Method to Calculate Water Saturation from LogDerived JFunction in Carbonte Reservoir
Authors T. A. Obeida and Y. S. Al MehaoriCalculation of initial fluid saturations is a critical step in any 3D reservoir modeling studies. The initial water saturation (Swi) distribution will dictate the original oil in place (STOIIP) estimation and will influence the subsequent steps in dynamic modeling (history match and predictions). Complex carbonate reservoirs always represent a quit a challenge to geologist and reservoir engineers to calculate the initial water saturation with limited or no SCAL data available. The proposed method in this study combines core data (permeability) from 32 cored wells with identifiable reservoir rock types (RRTs) and log data (porosity and Swi) to develop drainage logderived capillary pressure (Pc) based on rock quality index (RQI) and then calculate Jfunction for each RRT which was used in the simulator to calculate the initial water saturation in the reservoir.
The initialization results of the dynamic model indicate good Swi profile match between the calculated Swi and the logSwi for 70 wells across the field. The calculation of STOIP indicates a good agreement (within 3% difference) between the geological 3D model (31 million cells fine scale) and the upscaled dynamic model (1 million cells). The proposed method can be used in any heterogeneous media to calculate initial water saturation as a function of height as well as rock property (RQI) and thus estimate proper original oil in place. This method is also applicable to initialize huge fine scale static models.



Multiscale Averaging Algorithms for Flow Modeling in Heterogeneous Reservoir
Authors A. K. Pergament, V. A. Semiletov and M. Y. ZaslavskyThe averaging algorithms for Lebedev’s grids and for nonorthogonal grids by the Samarskii support operator method are considered in the article. The averaging process allows integrating the grid cells in case the approximate solutions in the integrated cell may be determined. The functions describing the features of solutions have been obtained for integrated cells of the Lebedev and arbitrary nonorthogonal grids. It results from approximating fluxes in each cell for linear span elements of the functions mentioned above that the symmetric permeability tensor may be constructed and the flux approximation of the first order may be obtained. The flux approximation leads to the strong convergence of the algorithm. It is essential the system grid cells may be multiscale, i.e. sizes of cells and the values of permeabilities may be very different.
The flux approximation distinguishes this article results from those of other authors of averaging algorithms. Then the methods of flow modeling in anisotropic media have been developed.
