In this work, we investigate different approaches for history matching of imperfect reservoir models while accounting for model error. The first approach (base case scenario) relies on direct Bayesian inversion using iterative ensemble smoothing with annealing schedules without accounting for model error. In the second approach the residual, obtained after calibration, is used to iteratively update the covariance matrix of the total error, that is a combination of model error and data error. In the third approach, PCA-based error model is used to represent the model discrepancy during history matching. However, the prior for the PCA weights is quite subjective and is generally hard to define. Here the prior statistics of model error parameters are estimated using pairs of accurate and inaccurate models. The fourth approach, inspired from Köpke et al. (2017) , relies on building an orthogonal basis for the error model misfit component, which is obtained from difference between PCA-based error model and corresponding actual realizations of prior error. The fifth approach is similar to third approach, however the additional covariance matrix of error model misfit is also computed from the prior model error statistics and added into the covariance matrix of the measurement error. The sixth approach, inspired from Oliver and Alfonzo (2018) , is the combination of second and third approach, i.e. PCA-based error model is used along with the iterative update of the covariance matrix of the total error during history matching. Based on the results, we conclude that a good parameterization of the error model is needed in order to obtain good estimate of physical model parameters and to provide better predictions. In this study, the last three approaches (i.e. 4, 5, 6) outperform the others in terms of the quality of the estimated parameters and the prediction accuracy (reliability of the calibrated models).

The use of full-physics models in close-loop reservoir management can be computationally prohibitive as a large number of simulation runs are required for history matching and optimization. In this paper we propose the use of a physics-based data-driven model to accelerate reservoir management and we describe how it could be implemented with a commercial simulator. In the proposed model, the reservoir is modeled as a network of 1D flow paths connecting perforations at different wells. These flow paths are discretized and the properties at each grid block along each flow path are derived from history matching of production data. To simulate flow in this network model through a commercial simulator with all the physics, an equivalent 2D Cartesian model is set up in which each row corresponds to one of the 1D flow paths. Finally, the history matching is performed with ensemble smoother with multiple data assimilation (ESMDA). The proposed network model is tested on both waterflood and steamflood problems. It is demonstrated that the proposed model matches with well-level production history (including pressure and phase flow rate) well. The calibrated ensemble from ESMDA also provided a satisfactory probabilistic forecast of future production that almost always envelops the true solutions. This indicates that the proposed model, after calibrated with production data, is accurate enough for production forecast and optimization. In addition, the use of commercial simulator in the network model provided flexibility to account for complex physics, as demonstrated by the successfully application to the steamflood problem. Compared with traditional workflow that goes through the full cycle of geological modelling, history matching and probabilistic forecasting, the proposed network model only requires production data and can be built within hours. The resulted network model also runs much faster than a full-physics as it typically has much less grid blocks. We expect the proposed method to be most useful for mature fields when abundant of production data is available. As far as we know, this is first time a physics-based data-driven model is implemented with a commercial simulator. The use of commercial simulator makes it easy to extend the model for complex reservoir such as thermal or compositional reservoirs.

The need to deliver well-informed decisions within stringent production cycles in unconventional plays is motivating the quest for practical models that can assimilate increasing volumes of data and satisfactorily account for observed production trends. The present work introduces the extended Dynamic Mode Decomposition (EDMD) as a suitable data-driven framework for learning the reservoir dynamics entailed by flow/fracture interactions in unconventional shales. The proposed EDMD approach builds on the approximation of infinite dimensional linear operators combined with the power of deep learning autoencoder networks to extract salient transient features from pressure/stress fields and bulks of production data. The data-driven model is demonstrated on three illustrative examples involving single and two-phase coupled flow/geomechanics simulations and a real production dataset from Vaca Muerta unconventional shale formation in Argentina. Given relatively moderate data requirements, we show that it is possible to attain a high level of predictability from hidden field state variables and well production data. As the main conclusion of this work, EDMD stands as a promising data- driven choice for efficiently reconstructing flow/fracture dynamics that are either partially or entirely unknown, or that are too complex to formulate using known simulation tools on unconventional plays.

EOS-based phase equilibrium calculations are usually used in compositional simulation to have accurate phase behaviour. Phase equilibrium calculations include two parts: phase stability tests and phase splitting calculations. Since the conventional methods for phase equilibrium calculations need to iteratively solve strongly nonlinear equations, the computational cost spent on the phase equilibrium calculations is huge, especially for the phase stability tests. In this work, we propose artificial neural network (ANN) models to accelerate the phase flash calculations in compositional simulations. For the phase stability tests, an ANN model is built to predict the saturation pressures at given temperature and compositions, and consequently the stability can be obtained by comparing the saturation pressure with the system pressure. The prediction accuracy is more than 99% according to our numerical results. For the phase splitting calculations, another ANN model is trained to provide initial guesses for the conventional methods. With these initial guesses, the nonlinear iterations can converge much faster. The numerical results show that 90% of the computation time spent on the phase flash calculations can be saved with the application of the ANN models.

Important decisions in the oil industry rely on reservoir simulation predictions. Unfortunately, most of the information available to build the necessary reservoir simulation models are uncertain, and one must quantify how this uncertainty propagates to the reservoir predictions. Recently, ensemble methods based on the Kalman filter have become very popular due to its relatively easy implementation and computational efficiency. However, ensemble methods based on the Kalman filter are developed based on an assumption of a linear relationship between reservoir parameters and reservoir simulation predictions as well as the assumption that the reservoir parameters follows a Gaussian distribution, and these assumptions do not hold for most practical applications. When these assumptions do not hold, ensemble methods only provide a rough approximation of the posterior probability density functions (pdf 's) for model parameters and predictions of future reservoir performance. However, in cases where the posterior pdf for the reservoir model parameters conditioned to dynamic observed data can be constructed from Bayes’ theorem, uncertainty quantification can be accomplished by sampling the posterior pdf. The Markov chain Monte Carlos (MCMC) method provides the means to sample the posterior pdf, although with an extremely high computational cost because, for each new state proposed in the Markov chain, the evaluation of the acceptance probability requires one reservoir simulation run. The primary objective of this work is to obtain a reliable least-squares support vector regression (LS-SVR) proxy to replace the reservoir simulator as the forward model when MCMC is used for sampling the posterior pdf of reservoir model parameters in order to characterize the uncertainty in reservoir parameters and future reservoir performance predictions using a practically feasible number of reservoir simulation runs. Application of LS-SVR to history-matching is also investigated.

The objective of this work is to design novel multi-layer neural network architectures for simulations of multi-phase flow taking into account the observed data (e.g., production data) and physical modeling concepts. Our approaches use deep learning concepts combined with model reduction methodologies to predict multi-phase flow dynamics. The use of reduced-order model concepts is important for constructing robust deep learning architectures. The reduced-order models provide fewer degrees of freedom and allow handling the cases relevant to reservoir engineering that is limited to production and near-well data. Multi-phase flow dynamics can be thought as multi-layer networks. More precisely, the solution, pressures and saturations, at the time instant n+1 depends on the solution at the time instant n and input parameters, such as permeability, well rates, and so on. Thus, one can regard the solution as a multi-layer network, where each layer is a nonlinear forward map. The number of time steps is user-defined quantity, which will be treated as an unknown within our deep learning algorithms. We will rely on rigorous model reduction concepts to define unknowns and connections for each layer. Novel proper orthogonal basis functions will be constructed such that the degrees of freedom have physical meanings (e.g., represent the solution values at selected locations) and basis functions have limited support, which will allow localizing the forward dynamics. This will allow writing the forward map for the solution values at selected locations with pre-computed neighborhood structure that will be used in deep learning algorithms. In each layer, our reduced-order models will provide a forward map, which will be modified (trained) using available data. It is critical to use reduced-order models for this purpose, which will identify the regions of influence and the appropriate number of variables. Because of the lack of available data, the training will be supplemented with computational data as needed and the interpolation between data-rich and data-deficient models. We will also use deep learning algorithms to train the elements of the reduced model discrete system. In this case, deep learning architectures will be employed to approximate the elements of the discrete system and reduced-order model basis functions. The numerical results will use deep learning architectures to predict the solution and reduced-order model variables. Trained basis functions will allow interpolating the solution between the observation points. We show how network architecture, which includes the neighborhood connection, number of layers, and neurons, affect the approximation. Our results show that with a fewer number of layers, the multi-phase flow dynamics can be approximated. The proposed approach uses physical model concepts and deep learning methods to design a novel forward map, which combines the available data and physical models. This will benefit to develop a fast and data-based algorithms for reservoir simulations.

We develop a new data-driven model for the assisted history matching of production data from a reservoir under waterflood. Although the model is developed from production data and requires no prior knowledge of rock property fields, it incorporates far more fundamental physics than that of the popular capacitance-resistance model (CRM). The new model also represents a substantial improvement on an interwell numerical simulation model (INSIM) which was presented previously in a paper co-authored by the last two authors of the current paper. The new model, which is referred to as INSIM-FT, eliminates the three deficiencies of the original INSIM data-driven model. (1) For some complex cases, e.g., when a producer is converted to an injector or when injected water from more than one injector passes through an intermediate well node, the INSIM procedure for calculation of water saturation degrades to an ad hoc calculation which introduces inaccuracies. Our new model uses an accurate front-tracking procedure to calculate water saturation, hence the name INSIM-FT. (2) The original INSIM formulation assumes relative permeabilities are known a priori which defeats the objective of finding a model without requiring knowledge of petrophysical properties; INSIM-FT estimates relative permeabilities by historymatching. (3) Unlike CRM, the original INSIM model does not provide a reasonable characterization of how water from an injector is allocated among producers and thus does not reliably identify large-scale geological features such as faults. INSIM-FT remedies this INSIM deficiency. The reliability of INSIM-FT for history-matching, future reservoir performance prediction and reservoir characterization is validated with two synthetic models, and its performance is compared with that of CRM. Finally, INSIM-FT is applied to a field case.

Objectives/Scope History matching large, complex fields has remained time consuming, and valid probabilistic uncertainty quantification has been a distant goal. Significant advances have been made recently in the development of improved sampling algorithms, efficient and accurate proxy models, and extremely high-performance forward simulation using GPUs and scalable algorithms. Combining these, we present a highly efficient workflow for robust probabilistic forecasts of a real field model with detailed reservoir characterization and over thirty years of production history. Methods, Procedures, Process The model used here describes a highly faulted offshore field produced from two platforms with 11 distinct sand sub-intervals and 22 major faults. Geostatistical methods were used to estimate each reservoir's static property distribution based on well logs. This static 3D model was upscaled and combined with reservoir engineering data and historical-interpretations to provide initial fluid distribution for use in flow simulation history matching and forecasts. Hamiltonian Markov Chain Monte Carlo techniques coupled with an efficient implementation of proxy models are used for uncertainty quantifications. A fast, highly efficient and scalable GPU-based reservoir simulator is employed for all forward simulations. Results, Observations, Conclusions The complexity of the geological properties of the reservoir requires a large number of cells to adequately describe the dynamics of production, and many adjustable parameters are needed to characterize the uncertainties in the model. While the use of advanced sampling and proxies vastly reduces the number of detailed flow simulations needed, hundreds of forward simulations are still required to calibrate the proxies and production forecasts. The bandwidth and computational throughput provided by GPUs allow such simulations to be performed extremely fast using only a modest amount of hardware. Through this case study, we show that by coupling an advanced history matching, prediction and optimization tool with a fast GPU reservoir simulator, accurate probabilistic uncertainty quantification can be made practical even for large models. Novel/Additive Information The workflow presented here can potentially transform the way probabilistic uncertainty quantifications are performed on large real field models. It will lay the foundation for the next generation of fast and accurate uncertainty tools and workflows.

When many daily measurements of a thermal EOR field are taken throughout production, it is not cost effective to manually interpret trends to update reservoir models, so we developed an automated data-driven approach for production prediction using machine learning techniques. This is a two-step scheme that first predicts auxilary field measurements from directly-controlled field settings, then uses these predicted field measurements to predict production. The full two-step prediction process needs further refinement, but the second step alone shows promise for aiding in automated interpretation of data. Time shifts from daily seismic surveys improved production predictions.

History matching is commonly performed in reservoir simulations to calibrate model parameters and to improve prediction accuracy. History matching problems often have non-unique solutions, i.e., there exist different combinations of parameter values that all yield the simulation results matching the measurements. In such a situation, finding a single solution matching the observations does not guarantee a correct prediction for future production. Alternatively, a more reliable prediction should be made with an uncertainty quantification based on all different possible scenarios of the model parameters. Bayesian theorem provides a theoretical foundation to represent different solutions and to quantify the uncertainty with the posterior probability density function (PDF). Lacking an analytical expression, the posterior PDF are often shown with a sample of realizations, each representing a possible scenario. This paper presents a novel sampling algorithm aiming to deal with two commonly encountered difficulties in the sampling process. First, a typical sampling method requires intensive model evaluations and hence may cause unaffordable computational burden. To alleviate this burden, our algorithm uses a Gaussian process (GP)-based surrogate as an approximation of the computationally expensive reservoir model to speed up the sampling process. The GP surrogate is adaptively refined locally such that the necessary approximation accuracy is achieved with a minimum level of computational cost. Secondly, when the dependent relationship between observation variables and input parameters is nonlinear, the posterior PDF could be in a complex form, such as multimodal, which is difficult to sample from. To tackle this difficulty, a Gaussian mixture model (GMM) is used as the proposal PDF to explore the parameter space. The GMM is flexible to approximate different distributions and is particularly efficient when the posterior is multimodal. The developed approach is tested with an illustrative history matching problem and shows its capability in handling the above-mentioned issues. Multimodal posterior of the testing problem is captured and are used to give a reliable production prediction with uncertainty quantification. The new algorithm reveals a great improvement in terms of computational efficiency comparing previously studied approaches for the sample problem.

Oil production forecasting methods can be categorized mainly into reservoir simulation, empirical and analytical models. Simulation requires laborious reservoir modeling and history matching. Even so, the forecast might be unreliable owing to various uncertainties. Empirical models (including data-driven models) find correlations between oil cut and production, etc., but they lack of explicit knowledge of the physical behavior. Classic analytical models are loyal to reservoir physics. Nevertheless, they often require estimation of water saturation as a function of time, which is difficult to obtain for multi-well reservoirs. To tackle these deficiencies, we propose to infer fractional flow models from field data using the Koval (1963) theory. We inversely solved the Koval fractional flow equation to obtain a relationship between water cut and dimensionless time. By history matching field water cut data, two model parameters, the Koval factor and the producer drainage volume, are estimated. Nevertheless, it is challenging to use the Koval approach as a predictive model directly since the injection contribution into each producer in a future time horizon must be evaluated first. To address this issue, we combine the Koval approach with the Capacitance Resistance Model (CRM), which characterizes the injector-producer connectivities and response time. The material balance of fluids is established in a producer-based drainage volume to consider the contributions from nearby injectors and the time lag in production caused by reservoir/fluids compressibility. A regression approach is simultaneously advanced to minimize the model error. Because of robustly integrating the reservoir physical behavior and the data driven approach, the combination of the Koval theory and the CRM model can result in a synergy that leads to accurate oil rate predictions. We validated this integrated method in a synthetic heterogeneous reservoir to test its reliability and further applied it to a field case in western Venezuela. Both case studies demonstrate that this integrated model can be used as a real-time tool to characterize inter-well connection and predict future oil production accurately.

We derive and implement an interwell numerical simulation model (INSIM) which can be used as a calculation tool to approximate the performance of a reservoir under waterflooding. In INSIM, the reservoir is characterized as a coarse model consisting of a number of interwell control units, where each unit has two specific parameters: transmissibility and control pore volume. By solving the mass material balance and front tracking equations for the control units, the interwell fluid rates and saturations are obtained so that phase producing rates can be predicted. INSIM is applied to perform history matching for parameter estimation and to infer the interwell connectivity and geological characteristics. INSIM has the following advantages: (1) the model parameters estimated from history matching provide a relative characterization of interwell formation properties. The model can handle changes in the flow directions caused by changing well rates, including shutting in wells or converting producers to injectors, whereas with the common correlation-based interwell connectivity method, the well interactions are assumed to be fixed; (2) the previous methods that bear computational similarity to INSIM can only provide the total liquid production rate whereas, with our procedure, we can calculate the oil and water flow rates and hence history-match water-cut data; (3) because we can calculate the oil and water flow rates, our method can be used for waterflooding optimization but with far less computational effort than with the traditional method based on the use of a reservoir simulator.

Real-time reservoir management, optimization and uncertainty assessment requires an effective assimilation of all available data into multiple reservoir models to be able to generate accurate forecasts and in turn, fast and reliable decisions. This process entails thousands of complex simulations that generally require several days and weeks to be completed even with the most powerful supercomputer facilities available today. Consequently, there has been an increasing attempt to develop physics-based surrogate models that could eventually replace the job of high-fidelity simulations in an accurate and consistent way within optimization and uncertainty quantification workflows. The present work proposes a novel combination of physics-sound analytical growth models with data driven techniques to generate surrogate models that can replace or mitigate the computational burden entailed by intensive numerical simulations and thus, unlock fast and reliable forecasts to decrease the turn-around time within the field development decision-making process. The analytical model relies on the association of cumulative production curves with growth functions and diffusion phenomena that have been extensively used in the description of population and cell growth models in several branches of Life, Social and Economic Sciences. The use of growth function models enables description of production profiles from the complex interaction of main flow drivers across a connectivity network or continuous heterogeneous system. The data-driven component allows for discovering secondary physical flow or production trends that reside hidden in the data that in turn, may aid at complementing and extending the predictability of the whole surrogate model. The data-driven component relies on machine learning techniques to construct universal interpolators via radial basis function and/or artificial neural networks. Both analytical and data-driven approaches are conceived in a non-intrusive fashion; that is, neither an explicit knowledge of the equations governing the physics nor modification of flow simulation code is required. Hence, the proposed surrogate workflow can be realized from either field data or data generated from black-box commercial simulation software. The proposed class of non-intrusive physics-based surrogate models is promising for generation of moderate to long range forecasts for SAGD operations in oil sands as well as improving the reliability of forecasting and reserves estimation in unconventional resources. Demonstration of the proposed methodology is illustrated with a variety of forecasting scenarios and for predicting inflow performance relationships.

The present paper proposes a novel non-intrusive model reduction approach based on Proper Orthogonal Decomposition (POD), the Discrete Empirical Interpolation Method (DEIM) and Radial Basis Function (RBF) networks to efficiently predict production of oil and gas reservoirs. Provided a representative set of training reservoir scenarios, either POD or DEIM allows for effectively projecting input parameters (e.g., permeability, porosity), states (e.g., pressure, saturations) and outputs (e.g., well production curves) into a much lower dimension that retains the main features contained in the simulation system. In this work, these projections are applied across multiple levels to be able to collapse a large number of spatio-temporal correlations. It is observed that these projections can be effectively performed at a large extent regardless of the underlying geological complexity and operational constraints associated with the reservoir model. The RBF network provides a powerful means for developing learning functions from input-output relationships described by the reservoir dynamics entailed by multiple combinations of inputs and controls. In order to achieve a high degree of predictability from the resulting reduced model, the RBF network exploits locality by a means of Gaussian basis functions that are maximal at the sampled point and decrease monotonically with distance. Compared to multilayer perceptron networks (i.e., traditional artificial neural networks) RBF networks require less training and are less sensitive to the presence of noise in the data. In this regard, POD or DEIM acts as a data filter that additionally aids at designing a more compact RBF network representation suitable for targeting fast reservoir predictions. Numerical results show significant accelerations with respect to running the original simulation model on a set of field-motivated cases.

Multi-phase flow in porous media is governed by a non-linear parabolic equation for the pressure field and degenerate parabolic equations for the saturation fields. Although it is well understood that the pressure and saturation fields vary on different spatial and temporal scales, most conventional reservoir simulators apply the same discretization to both. Substantial computational advantages can be realized by using a coarse mesh, which can effectively capture the pressure variation, but does not provide sufficient accuracy for calculating variations in saturation. Hence, in this paper we investigate the use of higher-order methods to capture the sub-grid fine scale variations in the saturation profiles. Specifically, the moving least squares (MLS) technique [ 1 ], which is conceptually similar to multi-point-flux-approximation (MPFA) [ 2 ] techniques, is used to approximate higher-order hyperbolic and viscous fluxes in a multi-point fashion for finite-volume methods applied to porous media flow. The primary focus of this study is to examine the feasibility of higher-order schemes for the purposes of reservoir simulation. To this end, two test cases are considered; (1) 1D strongly hyperbolic tracer flow (2) 2D standard quarter-five spot problem to illustrate the advantages of the MLS scheme over traditional two-point flux approximation (TPFA) [ 2 ] schemes. Predictions show that the MLS scheme allows for computation of high-order derivatives of field variables for a Godunov-type approach to hyperbolic problems in an efficient manner. Unlike standard piecewise continuous reconstruction, interpolation of field variables in MLS is continuous across interfaces. This facilitates direct and accurate reconstruction of fluxes of viscous nature. In addition, because the MLS scheme has the desirable attributes of a MPFA scheme, K-orthogonality of the underlying finite-volume mesh is enforced by default.

In a recent paper, we introduced a data assimilation procedure based on the application of the ensemble smoother (ES) with multiple data assimilation (MDA). The method, which has the acronym ES-MDA, was theoretically motivated by the equivalence between single and multiple data assimilation for the linear-Gaussian case and an approximate relationship between ES and one iteration of the Gauss-Newton method. ES-MDA is easy to combine with any commercial simulator, in fact, much easier than it is to couple the ensemble Kalman filter (EnKF) with a commercial simulator. Here, we apply EnKF and ES-MDA to generate multiple realizations of the porosity, net-to-gross ratio (NTG) and permeability fields by history matching production and seismic impedance data (3D and 4D) in a turbidite reservoir in the Campos Basin. 3D seismic data are integrated first and then 4D and production data are assimilated together, either simultaneously with ES-MDA or sequentially in time with EnKF.

Reservoir simulation is used to optimize the design of multilateral wells, where the parameters include the well path, well position, the completions design (e.g. ICD configuration), the production constraints and uncertainty from multiple realizations of a geological model. The target application is the design of new multi-lateral complex wells within an established large study. The optimized well is included in the prediction scenario of a history matched model. The target full field models are very large (millions of cells) and have large amounts of wells (hundreds). The parameterization of the well path and position is performed by constraining the well path to surfaces, and geometric translation. Laterals are independently parameterized allowing the tie-in depth, angle and length to be parameterized. Completion design is parameterized by compartment length, number of valves per compartment and valve inflow diameter. The complex completion design relies on advanced well modeling features of the reservoir simulator. Production constraints, such as target oil rates and water limits are considered within the optimization. Finally, the uncertainty from multiple geological realizations is accounted for; this is achieved by introduction of a risk aversion factor which penalizes the objective function by a scaled standard deviation. Examples of optimization demonstrate long multi-lateral wells being optimized within a Full field study. The sensitivity of each parameter is studied, i.e. the effect of oil production on the well position, effect of modifying ICD inflow area, the ICD compartment length, the water handling limits, and the geological uncertainty. The major significance of this work is the provision of a system to enable engineers to produce an optimal well configuration given geological, production, well placement and completion (ICD) uncertainties in a timely fashion. The use of optimization assists the engineer in designing an optimal solution given many variables.

In this work, an approach for estimating non-Gaussian permeability field is developed, through a probabilistic collocation based Kalman filter (PCKF). In this approach, the polynomial chaos expansion is used to parameterize the non-Gaussian permeability field. The state variables are expressed by the polynomial chaos expansion whose coefficients are sequentially updated when observations are available. The probabilistic collocation method is employed to solve for the coefficients of the polynomial chaos expansion. The probabilistic collocation method is non-intrusive to the reservoir model, thus it allows the forward simulations to be performed independently with existing reservoir simulators, as in the Monte Carlo simulation. The applicability of the approach is demonstrated with black oil problems in spatially heterogeneous reservoirs. Continuous production data are used as observations for sequentially updating the permeability field. The PCKF approach is also compared with the ensemble Kalman filter (EnKF) for investigating the accuracy, efficiency and applicability. It reveals that the PCKF other than the EnKF can be used to efficiently estimate the non-Gaussian permeability field. While keeping the advantages of EnKF, such as the sequential updating and parallelism, the PCKF is more suitable for permeability fields characterized by non-Gaussian random fields.

Because of its ease of implementation, computational efficiency and the fact that it generates multiple history-matched models, which conceptually allows one to characterize the uncertainty in reservoir description and future performance predictions, the ensemble Kalman filter (EnKF) provides a highly attractive technique for history matching production data. In this work, we apply EnKF with a recently proposed method of covariance localization to history-match production data from a real field to generate multiple realizations of the permeability field. A single manually history-matched model is available for comparisons. Only 7.6 years of the 10 years of history were matched with the remaining 2.4 years of history used to assess the predictive capability of the history-matched models. For this field case, covariance localization was necessary to avoid the propagation of spurious correlations and loss of variance and also resulted in better data matches and predictions than were obtained with EnKF without localization. EnKF with covariance localization also gave better data matches, more accurate "future" predictions and far more geologically realistic models than were obtained by manually matching production data. We also present results obtained using half-iteration EnKF (HI-EnKF) with covariance localization. For this field case, HI-EnKF gave a significant further improvement in the data match and predictions. However, as HI-EnKF requires rerunning the ensemble from time zero at every data assimilation step, it leads to a considerable increase in the computational time. The results for this field case indicate that we can reduce the computational cost of HI-EnKF, without compromising the quality of the results, by rerunning the ensemble from time zero only when "large" changes in the state vector occur.

Although there is intense interest in employing the ensemble Kalman filter (EnKF) for assisted history matching, it is well known that the approximation of covariance matrices from a finite ensemble of states and a finite ensemble of predicted data vectors can lead to a large underestimation of uncertainty in the posterior probability density (pdf) for the state vector. Here, we demonstrate that, regardless of whether covariance localization is used, EnKF can lead to an overestimation of the uncertainty in future predictions of reservoir performance. This overestimation occurs because, even though the data matches obtained with EnKF tend to appear reasonable, these matches are much worse than those that can be obtained when history matching dynamic data with a gradient-based method. The relatively poor data match obtained with EnKF means the ensemble of states generated by assimilating data with EnKF are of relatively low probability. Under reasonable assumptions, a Markov chain Monte Carlo (MCMC) algorithm will theoretically generate a correct sampling of the posterior pdf conditional to dynamic data. However, standard implementations of MCMC are not computationally feasible as one must typically generate thousands of states in the chain to complete the transitionary burn-in period before one begins to generate states that represent samples from the target pdf, and each proposed state in the chain requires a run of the forward model (reservoir simulator) to evaluate the probability of accepting the proposal. In this work, the EnKF and MCMC methodologies are combined to obtain a relatively efficient algorithm for sampling the posterior pdf for the model parameters. In this algorithm, a symmetric square root of the state vector's posterior covariance matrix is calculated from an ensemble of state vectors obtained from EnKF. It is shown that this square root can be used both to propose new states in the Markov chain and evaluate the probability of accepting the state without running the simulator, which results in an efficient algorithm. Because the posterior covariance obtained after assimilating data with EnKF depends strongly on the initial ensemble, to obtain a good estimate of the posterior pdf, we generate several ensembles with EnKF and several Markov chains which are then combined and resampled with importance sampling. We demonstrate the potential utility of the proposed EnKF-MCMC algorithm on a small three-dimensional two-phase-flow reservoir example. The problem is sufficiently small that using only MCMC, we are able to generate a long-chain based on two million proposed states. The posterior pdf calculated from this long-chain is assumed to be true. The posterior pdf's for predicted cumulative oil and water production generated with the new EnKF-MCMC algorithm gives results that are in reasonable agreement with the posterior pdf's for the cumulative productions obtained from the long Markov chain, whereas the corresponding pdf's obtained using EnKF with and without covariance localization exhibit a significantly higher variance, i.e., overestimate uncertainty in the predicted cumulative oil and water production.