It is well-known that time-lapse seismic data can indicate changes in the rock elastic properties caused by reservoir depletion when the water-front evolves inside a producing reservoir. Calibration of subsurface structures is an essential step in the modelling of fluid flow processes in a plethora of geoenvironmental applications such as petroleum and geothermal reservoir monitoring. History matching is known to be an ill-posed inverse process. Its primary goal is to find suitable reservoir model parameters honouring observations (and possibly other constraints) by integration of dynamic data (e.g., gas and oil rate, bottom hole pressure, water cut and subsidence/uplift) and static data (e.g., core, logging, and seismic) (Aanonsen et al., 2009). Ensemble-based data assimilation methods have been used for Seismic History Matching (SHM) for calibrating subsurface models that are tuned and conditioned to dynamic observations (Luo et al., 2018).
Seismic is a vital tool used for reservoir monitoring, exploration and characterisation and management. When compared to dynamic production data which is normally utilised in history matching, seismic data is obtained in lesser frequency during production time but has a much denser sample space. Hence, seismic data can be used to provide vital information for reservoir model calibration to complement production data (Luo et al. 2018).
In recently published work, the ensemble Kalman filter (EnKF) (Evensen 1994, Law et al 2012, Stuart 2010), the Ensemble Smoother (ES) (van Leeuwen & Evensen 1996) and Ensemble Smoother with Multiple Data Assimilation (ES-MDA) (Emerick & Reynolds 2013), as well as iterative variants of ES (Chen & Oliver 2011, Iglesias 2013), have been investigated extensively for history matching/data assimilation problems. However, it has been observed that standard ensemble-based data assimilation approaches have limited capability in preserving non-Gaussian distributions of model parameters such as facies (Aanonsen et al.,2009, Villegas et al., 2018, Lorentzen et al., 2013). Instead, in these classical ensemble-based data assimilation techniques model parameters tend to lose any initially assumed nonGaussianity which might have been encoded in the starting distributions based on structural constraints observed in hard data. This affects that parameters become increasingly Gaussian during these evolution schemes (Evensen & Eikrem 2018, Kim et al. 2018).
Several approaches have been recently proposed in the literature to maintain structural non-Gaussian features in these techniques. For example, Zhou et al. (2012) use normal score transforms in the ensemble-based data assimilation techniques to preserve non-Gaussian features of model parameters. In that approach, non-Gaussian model parameters are converted first into Gaussian model-states parameters using a normal score transform. EnKF is then applied to these Gaussian parameters. When completed, the final Gaussian profiles are transformed back to the corresponding non-Gaussian parameter representations by back-transformation. In this logic, transformation can take advantage of parameterisation for the cases where the number of required transformed parameters is smaller than the number of original parameters, which helps to reduce computational complexity and to emphasise certain dominant features in the parameters.
A particularly popular model used in such parametrisation techniques is the discrete cosine transform (DCT), which has been employed in the framework of reservoir characterisation by (Jafarpour & McLaughlin, 2007) (Liu & Jafarpour, 2013). An alternative approach is provided by the level set technique (Moreno & Aanonsen 2007, Dorn & Villegas 2008, Villegas et al. 2018, Chang et al. 2010, Lorentzen et al., 2013) which is particularly useful for the modelling of reservoirs where sharp interfaces between different geological regions are present. Moreover, sparse geologic dictionaries (Etienam et al., 2018, Kim et al., 2018, Khaninezhad et al., 2012) including Fourier transform-based methods and specific DCT based techniques are capable of capturing essential features such as patterns and main shapes of a facies channel reservoir (Khaninezhad, et al., 2012). However, they often reveal a deficiency in describing a clear contrast among different facies because of potential data loss during the corresponding inverse transformations (Kim et al. 2018, Khaninezhad et al., 2012, Tarrahi & Afra, 2016).
In the current paper, I will employ a sparse coding approach which is closely related to parametrisation techniques described above. In this approach, sparsity related concepts are applied to both, petrophysical (in particular permeability) and seismic parameters. Sparse coding entails the process of constructing representation coefficients based on the target signal and suitably chosen dictionaries (Candes, et al., 2006). These dictionaries contain groups of features that are capable of compactly representing various pre-specified phenomena typically observed in the environment. Accordingly, in geological modelling, sparse geologic dictionaries are utilised to represent spatial models with a sparse linear combination of basis vectors, which constitute essential geologic features of a reservoir (Khaninezhad, et al., 2012). Such sparse dictionaries can then be employed inside an ensemble-based history matching reducing overall computational complexity whereby maintaining structural consistency with prior information. However, obtaining suitable dictionaries for such a purpose is not straightforward. Aharon et al. (2006) demonstrated the efficacy of an unsupervised learning approach, in particular, K-singular value decomposition (K-SVD), for this purpose. Later Sana et al. (2016) built an archive of essential geologic features (called sparse geologic dictionary) from several static reservoir models using K-SVD and used it for calibrating reservoir models with these dictionaries and EnKF.
A similar approach has been taken by Etienam et al. (2018) for calibrating a synthetic channelised reservoir, however using the ES-MDA instead as the data assimilation technique. In this work, the authors construct sparse geologic dictionaries of the permeability field using K-SVD. They incorporate a greedy algorithm, in particular, orthogonal matching pursuit (OMP), for the sparse representations of the static properties (permeability field) which then is updated by ES-MDA. In a similar but slightly different approach, Kim et al. (2018) couple an iterative form of sparse coding in a transformed space with ES-MDA to construct a set of geologically plausible models that preserve the non-Gaussian distribution of lithofacies in a synthetic channelised reservoir. DCT of sand-shale facies is then followed by repeated K-SVD steps to construct sparse geologic dictionaries. This preserves geologic features of the channelised reservoir such as pattern and continuity.
Let me point out that none of the above mentioned two approaches has integrated seismic attributes into the inverse problem, in contrast to the study presented here. I propose a hybridised ES-MDA algorithm that implements sparse coding to outperform previous seismic history-matching methods in particular for highly non-Gaussian model parameters. The fundamental contribution of this paper is the simultaneous use of DCT analysis for the chosen seismic data combined with K-SVD and OMP based sparse analysis for petrophysical model parameters, and both integrated into an iterative ensemble history matching algorithm, here the ES-MDA. The main ingredient of the idea proposed in this paper SHM framework resides in two aspects. First I reduce the size of big data (which is the acoustic impedance) while retaining sufficient information required for data assimilation; second I parameterise the permeability field with the over-complete dictionary constructed using the unsupervised learning algorithm K-SVD.
The remaining sections of this paper are arranged as follows. Section 2 first introduces the various methodologies being used, in particular, the forward problem, ES-MDA, K-SVD, OMP and DCT. Then the novel SHM technique incorporating unsupervised learning is described. Section 3 presents and discusses the numerical experiments for evaluating this new technique. Section 4 provides some conclusions and indicates possible directions for future work.
This section aims to discuss the various methods that would be used to model the novel SHM technique. Aspects of the forward and inverse problem being used, sparse coding, 4D seismic and DCT for domain compression of the seismic attributes will all be explained in this section.
2.1 Forward Problem
combination of Darcy’s law and continuity equation. In the following, the index notation is used for indicating the quantities relating to oil, water and gas respectively. The continuity equation is then written as (Craft & Hawkins 1991),
gravity. are the source and sink terms representing the phase production and injection,
and
are the phase saturation and phase pressure,
is the phase relative permeability,
is the phase density,
is the phase viscosity.
is the phase formation volume factor,
the solution gas-oil ratio,
the gas-oil capillary pressure.
is the vertical coordinate (height) and
denotes the magnitude of the gravitational acceleration. Next, I will discuss the ES-MDA algorithm which will be used in solving the history matching inverse problem.
2.2 ES-MDA
The history matching problem considered here can be formulated mathematically as
Where denotes the objective function of the history matching problem and
denotes the state vector composed of reservoir variables (e.g., permeability and facies) to be estimated. The typical expression of
for ensemble-based history matching problems is presented as (Emerick & Reynolds, 2013, Oliver et al 2008, Tarantola 2005)
In this notation denotes a known estimate for the state vector and the superscript b indicates the background (sometimes denoted prior) distribution;
denotes the covariance matrix of
denotes the observed responses;
is the dynamic vector composed of simulated responses constructed by running a reservoir simulator
for the state vector
; and
denotes the covariance matrix of observation error. The right-hand side of Eqn.9 is the addition of background and observation error terms (Emerick & Reynolds, 2013).
in principle can contain any unknown variables such as facies indexes, permeability/porosity field, coefficients of discrete cosine functions or sparse coefficients.
can be used to derive the minimum of the cost function and update equation for
as (Emerick and Reynolds 2013, Nocedal 1999, Sherman & Morrison 1950)
where the subscript denotes the
ensemble member;
denotes the cross-covariance matrix of
and
denotes the auto-covariance matrix of
is the coefficient to inflate
, which denotes the covariance matrix of the observed data measurement error (Emerick and Reynolds, 2013);
denotes the observation data perturbed by the inflated observed data measurement error; and
is the ensemble size (i.e., number of reservoir models as columns of the matrix in the ensemble). Conventionally, ensemble-based history matching updates
reservoir models simultaneously. In Eqn.10,
denotes the Kalman gain K, which is normally computed with
regularization by SVD using 99.9% of the total energy in singular values (Emerick and Reynolds, 2013, Oliver et al 2008, Law et al 2012, Hanke 1997). ES-MDA assimilates every state vector times using an inflated covariance matrix of measurement error (Emerick and Reynolds, 2013). In this case,
is the number of assimilations in ES-MDA.
and
are defined as follows:
Eqn 12 where denotes the mean of state vectors and
denotes the mean of dynamic vectors.
In ES-MDA, is normally constrained to
The second quantity on the right-hand side of Eqn 14 is known as the perturbation term. It reflects the uncertainty associated with data measurement and processing. The stochastic characteristics of are reflected by
is the random error matrix to observations, which is constructed with a mean of zero and a standard deviation of
, where
is the number of time steps found in the observations. For more details see (Emerick and Reynolds 2013).
The goal in this approach is to use sparse representations inside the ES-MDA workflow. Sparsity, concerning petrophysical parameters requires efficient tools which will be described in the next paragraph. In particular, I will explain the unsupervised learning algorithm K-SVD and sparse coding with OMP.
2.3 Sparse Coding
Sparse coding can be described as the process of computing the representation coefficients, X, based on the given signal and a dictionary
. This process, commonly referred to as atom “decomposition", is proven to be an NP-hard problem (Aharon, et al., 2006), which is very difficult to solve. Therefore, it is often replaced by a “pursuit algorithm" that finds an approximate solution instead. Among the most straightforward pursuit algorithms are the Matching Pursuit (MP) and Orthogonal Matching Pursuit (OMP) (Tropp & Gilbert, 2007). These are known as greedy algorithms which select the dictionary atoms sequentially. These techniques are straightforward, involving the computation of inner products between the signal and the dictionary atoms, and possibly deploying some least squares solvers or projections. In this work, I employ the OMP algorithm which is described further below.
2.3.1 K-SVD Algorithm for obtaining a dictionary
A key criterion of sparse coding is the identification of a basis (domain) in which the field or signal under consideration can be modelled as sparse (Sana et al., 2016) A sparse domain describes one whereby any signal located in that domain can be adequately represented by a linear combination requiring a few basis elements (Candes, et al., 2006). The Wavelets (Jafarpour, 2011) and the Discrete Cosine Transform (DCT) (Jafarpour & McLaughlin, 2007) are examples of such sparse domains. In this work, I construct a sparse domain that is suited for inverse problems in reservoir engineering (history matching). Such algorithms modified for history matching are known as dictionary learning algorithms (Elsheikh, et al., 2013).
The term ‘unsupervised learning’ is a part of machine learning that learns from test data which has not been labelled or classified. Unsupervised learning identifies patterns or clusters in the training data and predicts the trends in each new piece of data (Bishop 2006). With this in mind, one drawback of K-SVD is its large size of sparse geologic dictionaries (Khaninezhad et al., 2012). References from the literature show that sparse coding implementation with a transformation of parameter space could significantly reduce both computational complexity and costs that are simultaneously required for model calibration (Khaninezhad et al., 2012, Etienam et al., 2018, Kim et al., 2018, Sana et al., 2016).
For a set of training signals of size
(which is also the dimension of the unknown permeability field) the dictionary learning algorithm finds the basis elements necessary to construct a dictionary
whereby its sparse linear combinations represents each of the vector space found in the training set of signals
. The corresponding weights of the signals are given as
and are calculated by solving the optimisation problem,
and
are the individual training signals and coefficient vectors in
and
respectively
are the atoms of the dictionary
selected such that the
of
satisfies Eqn 18. Finding
and
requires us to solve a non-convex problem, which is difficult even in an approximate manner (Elsheikh, et al., 2013).
Notice that Eqn 15 amounts of a simultaneous minimisation of . A computationally feasible dictionary learning algorithm that can solve this problem is the K-SVD algorithm. The K-SVD is NPhard (level of computational complexity, i.e. non-deterministic polynomial time), and it solves the
equivalent of the problem enumerated in equation Eqn 15 by using a cyclic three-step approach (Aharon, et al., 2006) indicated in Figure 1. The first of the three steps selects the cluster mean using K-means clustering. The second step is the sparse coding stage for updating the coefficient of the matrix
for a given iterate of
making use of any pursuit algorithm (a sparsity inducing method). Here OMP will be used which is explained in the next paragraph. In the third step, each column of the dictionary elements
will be updated in a sequential order by changing the values of its columns and also updating the sparse coefficients
. This is done using the signals
such that
, see again Figure 1. For further studies on the K-SVD, the reader may refer to (Aharon et al., 2006).
2.3.2 Orthogonal Matching Pursuit
In the process of establishing a suitable over-complete dictionary, we need to be able to represent vectors in a sparse form by atoms of this proposed dictionary. For this purpose, I employ OMP which works as follows. For an estimate of sparsity for a field that is represented by a vector of weights
in a d-dimensional space, the OMP algorithm attempts to recover the signal
by using a linear combination of
basis atoms (basis elements) from such an over-complete dictionary
. The OMP algorithm finds such a sparse representation by solving the following optimization problem
Eqn 17 Practically it solves this optimisation problem by projecting the data iteratively onto a dictionary and then finding a set of basis elements that are maximally correlated with the residuals. is the
-norm and it represents the signal sparsity level of
is the error tolerance.
The OMP algorithm requires as inputs an estimate of sparsity for the sparse signal
, the dictionary
of rank
with
being the same as the size of the row of the sparse coefficients . The idea in OMP is in constructing an approximation to the signal
(the spatial permeability field) by undergoing an iteration process. During each iteration, the local optimum of the solution is computed. This technique is calculated by finding the column vectors in
(the over-complete K-SVD learned dictionary) which closely matches a residual vector indicated as
. The reader may refer to Tropp & Gilbert, (2007) for further information.
As already mentioned before, in this work the K-SVD combined with OMP is used to generate the over-complete dictionary. This is done off-line and once before the commencement of history matching with the ES-MDA. A schematic showing the implementation of the K-SVD is shown in Figure 1.
Figure 1: Schematic showing the K-SVD algorithm
approach
2.4 4D Seismic
success of SHM I assume that the PEM is perfect and use the soft sand model as the PEM (Luo et al., 2018, Mavko et al., 2009). The PEM assumes that the cement is deposited away from the grain contacts. It also considers that the initial framework of the un-cemented sand rock is a densely random pack of spherically shaped grains having the porosity around 36%, which is the maximum porosity value that the rock could have before the suspension. I denote this parameter as the critical porosity (Luo et al, 2018). The shear modulus
and dry bulk modulus (
) at critical porosity is then calculated using the Hertz-Mindlin model as denoted in (Mindlin, 1949, Luo et al, 2018)
where are Poisson's ratio, grain shear modulus and effective stress, respectively.
denotes the coordination number and is the average number of contacts per sphere,
is the degree of the root in this case. Here, the
and
are 9 and 3, respectively. The modified Lower Hashin-Shtrikman (MLHS) bound is used to find the effective dry moduli for a porosity value less than
(Luo et al, 2018, Mavko et al, 2009). MLHS connects two end points in the elastic modulus-porosity plane. One end point corresponds to zero porosity, the other end point (
), corresponds to critical porosity
taking the moduli of the solid phase, i.e. quartz mineral (
). For porosity value
between zero and
, the lower bound for dry rock effective bulk (
) and shear (
) are;
is solid/mineral bulk modulus and
The saturation effect is infused using the Gassmann model (Gassmann,1951). The saturated bulk modulus (
) and shear modulus (
are defined as;
is the effective fluid bulk modulus. For a gas- oil-water mixture (as is the case in the PUNQ-S3 synthetic model),
is given by
where ,
and
are the bulk modulus of water/brine, bulk modulus of oil, bulk modulus of gas, saturation of water/brine, saturation of oil and saturation of gas, respectively. The saturated density (Gassmann, 1951) can be written as;
Since the volumetric scale of seismic data is large, SHM often can be referred to as a “big data” assimilation problem. Some numerical issues may arise for ensemble-based history matching algorithms when assimilating big data, e.g., ensemble collapse of the ensemble when tracking the true signal and large costs in storing and computing the Kalman gain matrices (Luo et al., 2018, Skjervheim 2007). Also, the history matching may become over-determined with big data (such that there is no solution that could match the observational data exactly) as opposed to the under-determined nature required for most successful implementation of history matching which may affect accuracy of the algorithm results (Luo et al. 2018, Tarantola 2005, Skjervheim 2007). Therefore, data reduction is a necessary approach in big data analytics, as seen in geosciences (Luo et al., 2018). I will use DCT for addressing this task concerning the seismic data.
2.5 Discrete Cosine Transform
While OMP is employed as described above for obtaining a sparse representation of the petrophysical parameters, a straightforward DCT will be more appropriate in our case for enforcing a sparse representation of the seismic attributes due to its different physical meaning and scale. Due to the low resolution of seismic in the vertical dimension (as opposed to a very fine resolution in the other two dimensions), it is furthermore sufficient to average seismic parameters along this direction on the scale of the reservoir which allows us to use DCT in two dimensions.
The two-dimensional forward DCT of an input image for an output image
can be described as (William & Mitchell 1993);
2.6 Novel SHM technique
In this section, I will state the novel history matching approach abbreviated as SHM-KED which stands for Seismic History Matching incorporating K-SVD with ES-MDA and DCT. The algorithm for the procedure is given next and explained below;
End (For) End (For) End program The initial realisations for this algorithm are generated using the Multiple-Point Statistics (MPS) algorithm. 200 realisations of static reservoir parameters, which are the permeability fields, are generated based on the well data. The changes of the dynamic parameters, pressure and saturation, and the simulated pressure-production data are generated using a three-phase flow simulator to obtain an ensemble of simulated states. The simulated impedance maps associated with the simulate states for one particular time step are calculated using the given petro-elastic model. The simulated changes in seismic
attribute, in particular, impedance maps, and the well production data, give rise to an ensemble of simulated measurements. The true or reference measurements are the seismic attributes and the historical production data, in particular, real impedance maps, computed from a prior elastic fullwaveform inversion scheme. The DCT algorithm parametrises this seismic attribute by selecting the leading cosine basis of these properties. The K-SVD learned dictionary enforces prior structural information on the permeability field during the ES-MDA inversion step. The updated state vector contains the corrected static and dynamic states of the reservoir after a single iteration of the ES-MDA.
3.1 The computational test case PUNQ-S3
In this section, I will test the novel SHM approach on a synthetic reservoir model. The PUNQ-S3 is a small-size synthetic 3D reservoir engineering model (Floris et al. 2001). The PUNQ-S3 model is made up of 19×28×5 grid blocks, where 1761 are active grid blocks. The grid-blocks have equal 180-meter sides in x- and y-directions. The heights of the grid blocks are varying. Figure 2 shows the top layer of the synthetic reservoir, including six production wells (PRO-1, PRO-4 PRO-5, PRO-11, PRO-12, and PRO-15). The reservoir is bounded by a strong aquifer in the west and north and a sealing fault in the east and south. The is no requirement for an injection well as the aquifer ensures high reservoir pressure. The red zone in Figure 2 below denotes a gas cap. There are two sections for the production schedule. The first is considered a history matching phase, and it lasts for 8 years. This section consists of four years of production, one year of extended well testing and a 3-year shut-in period. The second section comprises of 8.5 years of production. Measurements from all six wells were used during the assimilation period. The measurements includes bottom hole pressures, gas-oil ratios, and water cuts. The shut-in pressure has a noise value of 1 bar whereas flowing pressure has a noise level of 3 bar. The gas-oil ratios have a measurement uncertainty of 10 % before gas breakthrough and 25 % after the gas breakthrough. The water cut has 2 % uncertainty before water breakthrough and 5% after. The noise in the measurements is assumed uncorrelated, and I have used the same noise level generating the observations as when running the algorithm. A commercial reservoir simulator (Schlumberger GeoQuest, 2014) is used for reservoir simulation. Facies type is observed at well locations.
Figure 2: permeability field (Log K) of the true PUNQ-S3 synthetic model showing the 2D planar permeability field and
3.2 Computer specification
The computer specification for the desktop computer used for conducting numerical experiments is shown in Table 1.
Table 1: Computer specification for running the numerical experiments
Table 1 shows the computer specification used for running the numerical experiments in this paper.
3.3 Performance measure RMSE
To quantify our production data match to the true observed data, I use the root-mean-square-error (RMSE) function for each ensemble member (i), which is defined as
: Number of data assimilation time steps where measurements are assimilated (measurement times)
: Number of data collected at each time step
: Ensemble member index
: Time index
: Metric or response (history matched metric or response)
: Observed data metric for metrics j (Data equivalent in state space ensemble) at time step k.
: Simulated data from simulator for metrics j (Data equivalent in state space ensemble) at time
step k. : Observed data standard deviation for metrics j (Data equivalent in state space ensemble).
3.4 Presentation and Discussion of Results
In figure 3, I generate the initial ensemble of horizontal permeability field (using FILTERSIM (Wu et al, 2006). To generate the initial ensemble of porosity and permeability (in the z-direction) , I use the following statistical relationship derived from cores of the true synthetic model.
For Layer 2
For Layer 3
For Layer 4
For Layer 5
From this the state variable for ES-MDA or SHM-KED becomes;
Figure 3: an initial model of permeability realisations #21, #11, #31 and #67 generated with the FILTERSIM algorithm
Figure 4: Schematic showing how DCT can be used for image compression.
Figure 4 shows how DCT can be used to generate a sparse representation of the true impedance field. The leading coefficient (at the bottom lower corner) represents more than 98% of the information in the original impedance field.
Figure 5: Some initial permeability K-SVD basis atoms recovered after dictionary learning
Figure 5 shows some K-SVD basis recovered after activating K-SVD on an initial pool of realisations. Notice the patterns being maintained in this dictionary basis. With this idea, I will enforce prior structural information on the permeability updates during ES-MDA. In this numerical experiment, I use constantly for.
where
for the ES-MDA data assimilation algorithm.
Figure 6: Acoustic impedance reconstruction showing the two survey times (first survey on-top, second survey below):
Figure 6 shows the accuracy of our method in reconstructing the true impedance. SHM-KED recovered better the impedance signatures compared to standard ES-MDA with no sparse adaptation and no impedance data being assimilated. The proposed method was capable of recovering the streaks found in the true model.
Figure 7(a-c) shows the initial ensemble responses of the well bottom hole pressure, water cut and gas oil ratio. Notice the wide uncertainty in the ensemble to the true data. This is expected as we have yet to commence data assimilation. For the ES-MDA algorithm, Figures 8(a-c), we notice the uncertainty reduced but the algorithm fails to match the water cut at PRO-12. For the SHM-KED algorithm, Figures 9(a-c) we notice the algorithm does better and matches the water cut at PRO-12.
Figure 7(a): shows the bottom-hole-pressure for the six producer wells of the initial ensemble. The red curve represents the
Figure 7(b): shows the Gas-oil-ratio cut for the six producer wells of the initial ensemble. The red curve represents the true
Figure 7(c): shows the water cut for the six producer wells of the initial ensemble. The red curve represents the true data, and
Figure 8(a): shows the well bottom hole pressure for the six producer wells of the ES-MDA ensemble. The red curve
represents the true data, and the cyan overlay lines represent the realisations. The vertical dashed line represents the historical
Figure 8(b): shows the Gas-oil-ratio for the six producer wells of the ES-MDA ensemble. The red curve represents the true
Figure 8(c): shows the water cut for the six producer wells of the ES-MDA ensemble. The red curve represents the true data,
Figure 9(a): shows the well bottom hole pressure for the six producer wells of the ensemble recovered by SHM-KED. The
red curve represents the true data, and the cyan overlay lines represent the realisations. The vertical dashed line represents the
Figure 9(b): shows the gas-oil-ratio for the six producer wells of the ensemble recovered by SHM-KED. The red curve
represents the true data, and the cyan overlay lines represent the realisations. The vertical dashed line represents the historical
Figure 9(c): shows the water cut for the six producer wells of the ensemble recovered by SHM-KED. The red curve
represents the true data, and the cyan overlay lines represent the realisations. The vertical dashed line represents the historical
the SHM-KED and ES-MDA history matched reconstructed permeability realisation to the true model. SSIM is an image metric quantifier that analyses the visual impact of three identities of an image: structure, contrast, and luminance ( and
). Explaining SSIM further, a value of 1 means complete similarity whiles the value of
indicates complete dissimilarity.
Figure 10(a): shows 2D mean permeability field reconstruction
Figure 10(a) shows the 2D layer by layer mean permeability field and compares the performance of the SHM-KED with the true model directly. We observe that the algorithm is capable of reconstructing the main streak of the true model. Figure 10(b) shows the 3D analogues to Figure 10(a) with the 6 producer wells.
(a) (b) (c) : 3D mean permeability field reconstruction (a) mean initial permeability field, (b) mean permeability recovered
To calculate the performance of the SHM-KED with ES-MDA, I use 2 performance indexes, namely the RMSE from Eqn 30 and SSIM . I then derive the combined NORM value as described in Eqn 32.
Table 2: RMSE values of some selected models
Table 3: Combined NORM values of mean permeability fields from initial, ES-MDA and SHM-
To have a complete view of how the SHM-KED outperforms ES-MDA we look at the RMSE values for the whole 200 realisations. In Figure 11(a), we see the initial RMSE of all the realisations with high values. ES-MDA (shown in Figure 11(b)) reduces these RMSE values with SHM-KED (shown in Figure 11(c)) reducing it even further.
Figure 11: RMSE comparison (a) RMSE for initial ensemble, (b) RMSE for ensemble recovered by ES-MDA and (c)
The SHM-KED algorithm can calibrate a non-Gaussian synthetic reservoir model by updating the sparse representations of the permeability field. The first initial library of geological permeability models was generated using FILTERSIM which was then used for creating the K-SVD over-complete dictionary. Dimensionality reduction of the seismic attributes, in particular, acoustic impedance, constructed using DCT is effective in reducing the size of the acoustic impedance library by selecting leading coefficients concerning the DCT transform basis. The permeability field is converted to sparse coefficients using this K-SVD dictionary together with OMP. The weight matrix obtained by decomposing the permeability field is imported to ES-MDA as a state vector. In each assimilation of ES-MDA, an update of weights results in reservoir models that are well conditioned to both static and dynamic data.
The main novelty of this paper is the incorporation of sparse representations of permeability and acoustic impedance updates which is practically achieved by using a combination of machine learning and ensemble-based data assimilation methods. History-matching results of the PUNQ-S3 synthetic reservoirs indicate that the developed method is capable in capturing the permeability streak of the true reservoir model and the same time reconstructing the acoustic impedance of the true model. We also observe an improved matching accuracy in both history and forecast in terms of well production, the reduced dispersion of production behaviours and permeability distribution, and the well-connected channel body of reservoir models with geological plausibility. ES-MDA with the dictionary update and sparse representation of the impedance yield higher matching accuracy values and lower dispersion values than ES-MDA incorporated without any sparse parametrisation. The increase in computational costs invested during the initial construction of the dictionary (done once and off-line) is affordable compared to the assimilation algorithms not coupled with sparse coding. Improving the matching accuracy of some producer wells remains an outstanding task for the proposed technique despite the overall enhanced matching quality. Future works will be to adopt an adaptive inflation factor during the ES-MDA data assimilation and a surrogate based modelling for the forwarding of reservoir production data to reduce the expensive forwarding sequence using a full fidelity run. Also, a supervised learning technique involving deep neural networks for constructing these surrogate models or for permeability parametrisation is also planned for future research.
This project has been supported by the Niger Delta Development Commission of the Federal Republic of Nigeria (NDDC) and the University of Manchester. Thanks are due to Schlumberger for providing the reservoir simulator.
Aanonsen, S., Oliver, D., Reynolds, A. & Valles, B., 2009. The Ensemble Kalman Filter in Reservoir Engineering--a Review. SPE Journal, 14(3), pp. 393-412.
Aharon, M., Elad, M. & Bruckstein, A., 2006. K-SVD: An Algorithm for Designing Overcomplete Dictionaries for Sparse Representation. IEEE Transactions on Signal Processing, 54(11), pp. 4311-4322.
Bishop. Christopher M 2006. Pattern Recognition and Machine Learning (Information Science and Statistics). Springer-Verlag, Berlin, Heidelberg.
Candes, E. J., Romberg, J. & Tao, T., 2006. Robust uncertainty principles: Exact signal reconstruction from high incomplete frequency information. Information Theory, IEEE Transactions on, 52(2), pp. 489-509.
Chang, H., Zhang, D. & Lu, Z., 2010. History matching of facies distribution with the EnKF and level set parameterization. Journal of Computational Physics, 229(1), pp. 8011-8030.
Chen, Y. and Oliver,D., .2011 Ensemble randomized maximum likelihood method as an iterative ensemble smoother, Mathematical Geosciences, Online First.
Craft, B. C. & Hawkins, M., 1991. Applied Petroleum Reservoir Engineering. S .l.: Prentice Hall.
Crichlow, H. B., 1977. MODERN RESERVOIR ENGINEERING-A Simulation Approach. 1st ed. Oklahoma: Prentice-Hall, Inc.Englewood Cliffs, New Jersey 07632.
Dorn, O. & Villegas, R., 2008. History matching of petroleum reservoirs using a level set technique. Inverse problems, p. 24.
Elsheikh, A., Wheeler, M. & Hoteit, I., 2013. Sparse calibration of subsurface flow models using nonlinear orthogonal matching pursuit and an iterative stochastic ensemble method. Advances in Water Resources, 56(1), pp. 14-26.
Emerick, A. A. and Reynolds A. C., 2013 Ensemble smoother with multiple data assimilation, Computers & Geosciences.
Etienam, C., Villegas, R., & Dorn ,O 2018. An Improved Reservoir Model Calibration through Sparsity Promoting ES-MDA. ECMI conference proceeding. 20th European Conference on Mathematics for Industry (ECMI) 18-22 June 2018, Budapest, Hungary
Evensen, G., 2003. The Ensemble Kalman Filter: Theoretical formulation and practical implementation. Ocean Dynamics, 53(4), pp. 343-367.
Evensen, G., Eikrem, K, S 2018. Conditioning reservoir models on rate data using ensemble smoothers.Computational Geosciences,22(5), pp.1251-1270.
Floris, F.J.T., Bush, M.D., Cuypers, M, et al; 2001.Methods for Quantifying the Uncertainty of Production Forecasts: A Comparative study Petroleum Geoscience, 7(S) 87-96, 1, https://doi.org/10.1144/petgeo.7.S.S87
Gassmann, F., 1951. Uber die elastizitat poroser medien: Vierteljahrss-chrift der naturforschenden gesellschaft in zurich.
Hanke, M.,1997. A regularizing Levenberg-Marquardt scheme, with applications to inverse groundwater filtration problems. Inverse problems 13(1), pp. 79–95
Iglesias, M.A., Dawson, C.2013. The regularizing Levenberg-Marquardt scheme for history matching of petroleum reservoirs. Computational Geosciences 17, pp.1033–1053
Jafarpour, B., 2011. Wavelet reconstruction of geologic facies from nonlinear dynamic flow measurements. Geosciences and Remote sensing, IEEE Transactions, 49(5), pp. 1520-1535.
Jafarpour, B. & McLaughlin, D. B., 2007. History matching with an ensemble Kalman filter and discrete cosine parametrization. Anaheim, California, SPE.
Khaninezhad, M. M., Jafarpour, B. & Li, L., 2012. Sparse geologic dictionaries for subsurface flow model calibration: Part 1, Inversion formulation. Advances in Water Resources, 39(1), pp. 106-121.
Kim S, Min B, Lee K,& Jeong H, 2018 Integration of an Iterative Update of Sparse Geologic Dictionaries with ES-MDA for History Matching of Channelized Reservoirs. Geofluids, Volume 2018, Article ID 1532868, 21 pages.
Law K J H & Stuart A M, 2012 Evaluating Data Assimilation Algorithms Mon. Weather Rev 140 37-57
Liu, E. & Jafarpour, B., 2013. Learning sparse geologic dictionaries from low-rank representations of facies connectivity for flow model calibration. Water resources, Volume 49, pp. 7088-7101.
Lorentzen, R. J., Flornes, M. K. & Naevdal, G., 2012. History matching Channelized Reservoirs Using the Ensemble Kalman Filter. Society of Petroleum Engineers, 17(1).
Luo X, Bhakta T, Jakobsen M, Naevdal G, 2018. Efficient big data assimilation through sparse representation: A 3D benchmark case study in petroleum engineering. PLoS ONE 13(7): e0198586. https://doi.org/10.1371/journal.pone.0198586
Mavko G, Mukerji T, Dvorkin J.2009. The rock physics handbook: Tools for seismic analysis of porous media.Cambridge University Press.
Mindlin RD.1949 Compliance of elastic bodies in contact. Journal of Applied Mechanics.; 16:259-268
Moreno, D. L., 2011. Continuous Facies Updating Using the Ensemble Kalman Filter and the Level set method. s.l., Mathematical Geosciences.
Nocedal, J. and Wright, S.J, 1999. Numerical Optimization, Springer, New York.
Oliver, D. S. & Chen, Y., 2010. Recent progress on reservoir history matching: a review. s.l.: Computational Geoscience - Springer Science.
Oliver, D. S., Reynolds, A. C. & Liu, N., 2008. Inverse Theory for Petroleum Reservoir Characterization and History Matching. s.l.: Cambridge University Press.
Sana, F, Katterbauer, K , Al-Naffouri, T.Y and Hoteit, I.,2016 Orthogonal matching pursuit for enhanced recovery of sparse geological structures with the ensemble Kalman filter, IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 9, no. 4, pp. 1710–1724
Schlumberger GeoQuest, 2014. ECLIPSE 100 (Black Oil): Reference Manual and Technical Description. Houston: s.n.
Sherman, J., Morrison, W.J.1950. Adjustment of an inverse matrix corresponding to a change in one element of a given matrix. The Annals of Mathematical Statistics 21, 124–127.
Skjervheim J-A, 2007. Continuous Updating of a coupled reservoir-seismic model using ensemble Kalman filter technique, Bergen, Norway: PhD Thesis, University of Bergen.
Spanias AS, Jonsson SB, Stearns SD. 1991; Transform methods for seismic data compression. IEEE Transactions on Geoscience and Remote Sensing. ; 29(3):407±416. https://doi.org/10.1109/36.79431
Stuart A M ,2010 Inverse problems: A Bayesian perspective Acta Numerica 19 451-559.
Tarantola, 2005. Inverse problem theory and methods for model parameter estimation. 1 ed. Philadelphia: SIAM.
Tarrahi, M. & Afra, S., 2016. Improved Geological model calibration through Sparsity-promoting Ensemble Kalman Filter. Offshore Technology Conference.
Tropp, J. A. & Gilbert, C. A., 2007. Signal recovery from random measurements via orthogonal matching pursuit. Information Theory, IEEE Transactions on, 53(12), pp. 4655-4666.
Villegas, R., Etienam, C., Dorn, O., & Babaei, M.,2018. The shape and Distributed Parameter Estimation for History Matching using a Modified Ensemble Kalman Filter and Level Sets. ; Inverse problems Science and Engineering.
Wang, Z., Bovik, A. C., Sheikh, H. R. & Simoncelli, E. P., 2004. Image quality assessment: from error visibility to structural similarity. Image processing IEEE transactions, 13(4), pp. 600-612.
William P, Mitchell J,1993 JPEG: Still Image Data Compression Standard, Van Nostrand Reinhold,
Wu, J., Boucher, A. & Journel, G. A., 2006. A 3D code for mp simulation of continuous and categorical variables: FILTERSIM. SPE.
Zhou, H, Li, L and. Gómez-Hernández, J.,2012, “Characterizing curvilinear features using the localized normal-score ensemble Kalman filter,” Abstract and Applied Analysis, vol. 2012, Article ID 805707, 18 pages.