Workflow for adding 4D seismic data in history matching

In this document we present a workﬂow for ensemble-based 4D seismic history matching. Ensemble-based history matching has become standard for production data, but 4D seismic data poses a number of additional challenges. One issue is that the amount of data is considerably larger, but another, probably more complicating factor is that for utilizing the seismic data, either the seismic data must be inverted to properties that is included in the reservoir simulation model, or a seismic response must be modeled, given the current estimate of the reservoir properties. This leads to a number of choices on how to utilize the information of the 4D seismic data. We will discuss this, as well as point to approaches for handling large amounts of data in ensemble-based history matching. The developed approach has been applied on the Norne ﬁeld and is currently being evaluated at the Ekoﬁsk ﬁeld. This document is primarily addressed to reservoir engineers and researchers that are working on history matching 4D seismic data, but it might also be of interest to those working with 4D seismic data from a geophysical perspective. After all, 4D seismic history matching should be viewed as an interdisciplinary subject. Although, our focus has been on ensemble-based history matching, some of the choices that have to be made in utilizing 4D seismic data is independent of the actual method used for history matching.


Introduction
Although history matching of production data is a common standard in the industry, and there exist a large amount of 4D seismic data, utilizing the 4D seismic data for history matching is not considered as a standard workflow [25].The 4D seismic data is mostly used for interpretation to find regions of the reservoir where changes in the pressure or fluid content have occurred.However, this means that data sets that would have a large degree of complementary information to the production data are not utilized to improve the reservoir models, and thereby contribute to improved reservoir forecasting.However, utilizing 4D seismic data for history matching has been a grand challenge, and developing methodology for solving this challenge has been a major undertaking at the IOR centre for its full lifetime.In this document, we will present major findings, and contributions in solving this challenge.
For testing the developed methodology, two fields have been selected.First we have worked on the Norne field.It was decided in communication with Equinor that starting out with the data set from Norne that was publicly available [22] was a suitable starting point.Most of the development has been done based on this data set, including some additional information from Equinor.As a second data set, a model of Ekofisk (provided by ConocoPhillips) has been utilized.Although, it might have been advantageous to test the methodology on more data sets, this has not been possible due to the overhead cost of getting started on working with such large data sets.
Developing a workflow for 4D seismic history matching consists of a large number of separate tasks, as we have illustrated in Figure 1.As we have considered history matching utilizing production data a matured topic, the focus has been on handling seismic data (the blocks "Rock physics model parameters", "Sim2seis", "Simulated seismic attributes" and "Adjust parameters").For the first three of the mentioned blocks an important decision that has to be made is on which format the seismic data should be used.For the block "Adjust parameters", a particular challenge is how to address large amount of data while doing the update part.Based on our previous experience, we have only considered ensemble-based methods for updating the reservoir simulation model.This seems natural, as we would like to change a large number of parameters in the model [24].

Methodological Approach
We have split the description of the methodology into several subsections in the document.First we discuss how to handle the seismic data and the choices that must be done in that respect in the two sections "Seismic data integration" and "Approach to integrate seismic data".For ensemble-based history matching generation of the initial ensemble is an important part, and the section "Generation of initial ensemble" is devoted to that.Thereafter we discuss how to handle the corresponding large data sets while adjusting the parameters (in the ensemble-based approach) in the sections "Sparse data representation" and "Correlationbased localization", before we discuss about some available choices for ensemble-based methods for adjusting the parameters in the section "Ensemble methods".

Seismic data types
When performing an assisted history matching using seismic data, we can integrate the seismic data in various levels, i.e. (i) the seismic level, (ii) the elastic parameter level and (iii) the simulation model level [25].The amount of work for either on forward modeling or on real data preparation depends on the level of data integration.For example, in the seismic domain level no inversion of seismic data is required, rather observed seismic data are compared directly with simulated seismic data.However, one needs to generate synthetic/simulated data which could be a complex and time consuming process.Further, the quality of the simulated data depends on various factors such as a good petro-elastic model (PEM) and/or sim2seis model as well as the quality of the underlying reservoir or geological models.In the elastic parameter level, the simulated elastic parameters from the reservoir model are compared with inverted impedance or elastic parameters.Although the need for forward modeling is reduced, data integration in this level is still challenging as the seismic inversion process can increase the uncertainty in the inverted seismic data.In the simulation model level, the seismic data are inverted into reservoir parameters/states such as changes in pressure and saturation fields and then directly compared with reservoir simulation outputs.This level does not require any seismic forward modeling.However, the biggest challenge is the large uncertainty in the inverted saturation and pressure fields from the seismic data.

Approach to integrate seismic data
In good practice, there is a need for ensuring the validity of both the modeling and the inversion schemes to be able to compare at all these levels.First of all we need to ensure that the modeling is consistently going from the reservoir model to the seismic data set taking all relevant aspects into account.On the other hand, we should ensure that the inversion scheme actually can go from seismic to the reservoir scale, which we are comparing on, in a consistent way.In the end we should choose the dataset which contributes with highest information content in our case.Selection of one or more attributes for history matching should be based on a presensitivity analysis to decide how dynamic changes affect the potential attributes.Interpretation and understanding of 4D seismic data further increases ones confidence on the selected attribute.Extraction of a proper attribute depends on other factors too, for example, increase of precision of the extracted attribute escalates both the demand for computational power as well as the complexity of the seismic analysis.Further, the potential of the extracted attributes is partially dependent on the goodness of corresponding simulated attributes, thus on seismic and rock physics forward modeling.Therefore, deciding on attributes is somewhat constrained by the availability of forward seismic modeling.For example, 2D map based attribute is the most robust attribute to select, though, it doesn't contain any vertical information.On the other hand, 3D attribute contains vertical information, however might have more challenges in up-scaling, interpolating between reservoir and seismic domain [25].
In most cases, even after performing proper inversion and forward modeling as good as you can, there still remains discrepancies between modeled and read seismic data/attributes.There could be many reasons for the deviations, such as structural discrepancy in the static model, simplification in the PEM, bias in dynamic parts in the PEM, etc., [23].Therefore, there is always a need to scale/normalize both modeled and real field data at the final stage.The scaling can be performed in various ways using mean, median or any percentile value.Further, the most important data can be extracted by using a mask before scaling.
The process of creating a mask to put different weight on different aspects of the data should be guided by domain understanding.Geophysicists should be able to provide information about which 4D seismic observations are significant and interpretable, and which observations are more questionable, e.g.observations in regions with a poor signal-to-noise ratio, seismic interference or processing issues.Reservoir engineers should be able to provide information about significant well events within each survey interval, e.g.well start-up, shutdown, significant changes in pressures and/or rates, breakthrough, etc.The mask approach gives an option to put more weight on 4D effects that are significant and interpretable, e.g.4D effects close to important wells.The cross-domain process resulting in a suitable mask/set of masks will increase our confidence on the assimilated data further.
Before calibrating the model, one should check the adequacy of the model by comparing actual observations with simulated data.This model diagnostic process would help us to identify the source of deficiency in the model, such as identifying the source of conflict between the observed and simulated data, outliers etc.Here, an important question is what level of simulated differences are required to produce a significant 4D effect, and whether the model tends to over-predict or under-predict 4D effects.By checking whether the remaining discrepancies between models and data are within the range explained by model parameter variations, we can decide whether to accept or reject the model instead of performing the full history matching and conclude about the inadequacy of the model [25].

Generation of initial ensemble
The initial ensemble should be generated based on the available information, which might be from a number of sources, as interpreted seismic, information from wells, like core samples and well testing, analog outcrops, and other information available.This information might be summarized into a geostatistical model that ideally should be used for generating an ensemble of models.For the two ensemble-based history matching studies presented later (Norne and Ekofisk) the geostatistical model was not available, and the prior ensemble had to be based on a provided model, that was used to give an initial mean of the ensemble.To create an ensemble of models, the ensemble members was created by perturbing the mean model with perturbations in the permeability and porosity fields using a variogram model.The generation procedure is available at https://github.com/rolfjl/Norne-Initial-Ensemble.
An important question is whether the ensemble approach captures and preserves the geological interpretation of the reservoir.An example of significant geological information is the location of faults, fracture corridors, channels, thief zones, etc. which may have a significant impact on the fluid flow, e.g. by reducing the permeability (e.g.sealing fault) or enhancing it (e.g.open faults / fractures).In many cases, 3D seismic attributes can be used to automatically locate these features, but additional information is needed to determine appropriate permeability levels for these features.If the initial ensemble is parametrized according to the geological interpretation (features populated with either high or low permeability locally), the assimilation process will typically not alter the location of the features.Hence, the geological interpretation can be preserved.Using well information and seismic data to characterize faults and fractures, which can be used as prior information in ensemble modeling, is discussed in [2].

Sparse data representation
One feature of spatially distributed seismic data (e.g., AVA, impedance) is their huge volumes.While seismic data provide valuable information for reservoir characterization, they also impose numerical challenges, in terms of computational resources that are required to handle such big datasets in ensemble-based history matching workflows.
To handle the issue of big seismic data, we adopt the concept of sparse data representation, of which the main idea is to represent the original seismic data more tightly in another domain.While in the literature there are different methods of implementing the idea of sparse data representation, in our work [16,17,12,13] we have focused on using discrete wavelet transform (DWT), as this is a well developed method and computationally very efficient in large-scale problems.
Figure 2 illustrates the procedure of sparse data representation through DWT.Given a seismic dataset, which can be 2D, 3D, or even 4D, one first conducts a DWT to obtain a set of wavelet coefficients, which are the representations of the original seismic data in the wavelet domain.For the purpose of sparse representation, we choose to only keep those wavelet coefficients whose magnitudes are above certain threshold values, while setting to zero the wavelet coefficients whose magnitudes are below the threshold values.During history matching, we only use the kept leading wavelet coefficients as the observations.By doing so, we can significantly reduce the size of observation data, which helps to substantially reduce the consumption of computational resources.
The threshold values for sparse data representation are determined based on the estimated noise levels of the wavelet coefficients.Based on the estimated noise levels, we can construct as a by-product the observation error covariance matrix [14], which is a required input in order to deploy an ensemble history matching algorithm.An additional remark is that the procedure of sparse data representation is generic and flexible, and can be applied to handle different types of seismic data (e.g., AVA, impedance), or even other types of spatially distributed field datasets (e.g., electromagnetic data).As such, we also foresee the possibility to extend the same procedure/history matching workflow to deal with geophysical inverse problems in general.

Correlation-based localization
When applying an ensemble history matching algorithm to large-scale problems, it is common to use a relatively small ensemble (of order 100) of reservoir models to reduce the computational cost.One consequence of a small ensemble size is that there would be substantial sampling errors, which often deteriorates the performance of an ensemble history matching algorithm.To tackle this problem, an auxiliary technique, called localization, is introduced.
In a conventional ensemble history matching workflow, a localization scheme uses the distances between the physical locations of both observations and model variables to determine whether an observation data point should be used to update a given model variable, or not.This type of localization schemes is referred to as distance-based localization.
While distance-based localization has been shown to work well in many history matching problems, there are also certain noticed issues, including, e.g., the need for physical locations of both observations and model variables to compute the distances, the difficulties to handle observations that may have non-local influences on model update, the inconvenience to deal with large field datasets, and the deficiency of ignoring the time-lapse effects of field data.Some of the aforementioned issues may affect the applicability of distance-based localization to seismic history matching problems.This is particularly true if we choose to represent seismic data in the wavelet domain.In this case, the observations consist of a set of leading wavelet coefficients, which do not possess clearly defined physical locations.
To circumvent the difficulties of distance-based localization, we have developed a new type of methods [18,15], in which, instead of relying on physical distances, it is the sample correlations between the simulated observations and model variables that are adopted to construct the localization scheme.For this reason, this new type of localization schemes is referred to as correlation-based localization.As elaborated in [18,15], correlation-based localization is able to overcome the aforementioned issues that may arise in distance-based localization, and its efficacy is demonstrated in a few field case studies [19,20,13].

Ensemble methods
In this section we give an overview of the most common ensemble methods used for assimilating seismic data.The text in this section is based on a summary of the methodology presented in Oliver et al. [25].
Currently the industry standard is use of iterative ensemble smoothers that assimilate all available data for each iteration.The benefit compared to traditional Ensemble Kalman Filters (EnKF) is that time consuming restarts of the flow simulator are avoided.For detailed reviews of the traditional ensemble Kalman filter in reservoir engineering we refer to Aanonsen et al. [1] and Oliver and Chen [24].
The most common objective function used for history matching is obtained from Bayes rule for the computation of the posterior probability distribution for model parameters [32].If the prior uncertainty in model parameters can be modeled adequately as multivariate Gaussian with covariance C m , and if the observation errors are additive and Gaussian, then to compute the most probable vector of model parameters, m j for ensemble member j, one need only solve for the minimizer of In the above equation each m pr j represents an ensemble member from the prior distribution and d obs j = d obs + j are perturbed observations using samples j from the data error distribution N (0, C d ).The forward model g(m) maps model parameters to simulated data.An updated ensemble (m j ) of model realizations are computed by minimizing the objective function for each ensemble member, although it it is known that, for a general non-linear observation operator the distribution of samples obtained by minimization is only an approximation of the posteriori distribution.A common feature when deriving the ensemble smoothers is use of a first order Taylor approximation.An exception from this approach is a recently developed ensemble subspace formulation [26,9] of the Ensemble Randomized Maximum Likelihood method (EnRML) [3].In the subspace formulation approximations are justified using linear regression.
The Ensemble Smoother with Multiple Data Assimilation (ES-MDA) ES-MDA was introduced in Emerick and Reynolds [7] as an approach to improve the traditional Ensemble Smoother (ES) [29] when working with non-linear problems.Instead of performing a single update step, the algorithm assimilates the data N MDA times with inflated noise variance.In the linear Gaussian case, it is shown [7] that multiple data assimilation is equivalent to the single ES data assimilation step.The algorithm is based on the fact that the likelihood term in Bayes formula can, in the Gaussian case, be rewritten as a product of N MDA exponential terms with inflated data covariance.The resulting update formula for each ensemble member (j = 1, . . ., N ) is then given by: where i = 0, . . ., N MDA −1, and N MDA −1 In the above equation we use the definitions: and ḡi = N −1 N 1 g(m i j ).The ES-MDA method was utilized for 4D seismic history matching of a heavy-oil turbidite reservoir in Campos Basin [8] where results were compared to the results obtained using the standard EnKF.The performance of the methodology was further analyzed in Emerick [6], where the same field was investigated, but ES-MDA was compared to the standard ES.The conclusions from the papers are that ES-MDA outperformed both the EnKF and the ES in the joint assimilation of production and seismic data, but that the reduction in the ensemble variance was excessive.There exist numerous other applications of ES-MDA for 4D seismic history matching [e.g., 10,5,33].In general, the mismatches to both seismic and production data were reduced at reasonable computational expense.
The Ensemble Randomized Maximum Likelihood method (EnRML) EnRML [3] searches for the minimum of the objective functions (Eq. 1) using an ensemble-based approximation of the Gauss-Newton method.It does not rely on a predefined number of assimilation steps, and differs in that sense from the ES-MDA.The method does, however, require a convergence criteria and the step length parameter is determined by standard line search.The methodology was later improved using the Levenberg-Marquardt algorithm (LM-EnRML), for better selection of the step size and faster convergence [4].However, the original LM-EnRML involves the inverse of the state covariance matrix ( Ci mm ), which is found to be sensitive to the level of truncation in the singular value decomposition (see also the discussion below), which has lead to the most widely used form of the LM-EnRML ignoring the updates from the model mismatch term ( Ci mm ) −1 (m i j − m pr j ) (see Chen and Oliver [4] for details): Here S m and S d are the same as for ES-MDA and given by Eqs. 3 and 4. The Levenberg-Marquardt tuning parameter is denoted λ i .For details regarding adaptive updates of λ i we refer to Chen and Oliver [4].
A flavor of the EnRML is presented in Luo et al. [21].In that work a Regularized Levenberg-Marquardt algorithm is used to solve a Minimum Average Cost problem (RLM-MAC).The two methods are compared in Luo et al. [21] and reports better performance for the RLM-MAC algorithm, but as the authors clearly states, the conclusion may be different if the methodologies are investigated on a broader set of experiments.The RLM-MAC was used for sequential assimilation of first production data and then seismic data in a digital twin experiment based on the Norne field [12].That study showed that it was possible to obtain good data match for both data types.In a subsequent study using real data from the Norne field, Lorentzen et al. [13] showed that it is also possible to assimilate both data types simultaneously, with reduction in mismatch for both production and seismic data.

Validation
Here we report several case studies, ranging from simple 2D synthetic to more complicated 3D real-field case studies, in which our proposed workflow (or part of it) was developed and tested.

Synthetic studies in Norne and Brugge models
The concept of DWT based sparse data representation was first developed and tested in a synthetic 2D Norne model [16].In this study, amplitude versus angle (AVA) data was used as the seismic attribute, and was history-matched by an iterative ES [21].It was shown that through DWT based sparse data representation, one can significantly reduce the data size, while keeping the most important data features.In addition, adopting both AVA and production data led to better history-matching performance than using either type of data only.In a follow-up paper [17], the efficacy of DWT based sparse data representation procedure was further validated in a 3D Brugge benchmark case study, where the 4D seismic history matching workflow of [16] was adopted, and the same conclusions as those in [16] were obtained.
In [17], it was also reported that the iterative ES suffered from adverse effects of a relatively small ensemble size.This led to subsequent developments of correlation-based localization schemes [18,15,31].As the origin of these developments, the 3D Brugge model with the same settings as in [17] was among the selected case studies to demonstrate the performance of correlation-based localization schemes, while a 3D Norne model was also adopted in a few other studies [19,12].In all these case studies, correlation-based localization was observed to perform well, help mitigate the aforementioned adverse effects, and work very efficiently for field-scale reservoir models.

Real-field study using the Norne model
In [13] both production data and seismic data were assimilated simultaneously for the Norne model.Angle-versus-offset data were inverted for acoustic impedance using a Bayesian approach.The data were further processed by taking the difference between the monitor surveys, and the base survey.The resulting differences were averaged over reservoir formation layers.Sparse representation based on truncation of wavelet coefficients were used, as described above and in [16].Further, issues related to limited ensemble size and erroneous covariance matrices are handled using the localization technique based on correlations between model parameters and measurements, see [18].This approach is flexible and can be applied to any data type, and for global parameters.That is crucial for our application as the wavelet coefficients do not have a spatial position and parameters such as those describing relative permeability curves are global.It is not possible to apply traditional distance based localization for these data and parameters.In addition, the correlation based localization method is easy to apply, as very little manual work and user input are required.The method is also adaptive, and an updated localization domain is computed every time model parameters are updated.
The results from this work show that we reduce the data mismatch for both production and seismic data, see Figure 3. Inspection of the data match for individual wells also reveal that the updated models outperform the manually history matched model.In addition, the updated static petrophysical fields are geologically credible, and ensemble collapse of model parameters is avoided.In addition to production forecasts, the updated models can be used to simulate flow in unswept parts of the reservoir, and have potential for providing useful information when planning infill wells or EOR strategies.Real-field study using the Ekofisk model For our second real field study, the Ekofisk reservoir model is considered.The field is a compacting reservoir which further adds complexity in the history matching.Here, we have considered time-lapse acoustic impedance data as the observation/ measurements in addition to production data.Here, both the seismic and production data are assimilated simultaneously.Two life-of-field seismic (LOFS) data difference of 3 years interval is considered here for the data assimilation.As the simulated acoustic impedance data does not fit well with the inverted real field acoustic impedance in each grid block, we only consider data that are close to the well location.Here, we define a well-mask that only extracts data near to the well location masking the rest away from the wells.Further, our newly developed correlation based localization technique [18] is used to mitigate ensemble collapse.The results from this work show that the data mismatch for both production and seismic data reduce if we con-sider well-mask to extract seismic data rather using seismic data coming from all the reservoir grids.

Conclusions and recommendations
A recommended practice has been presented for adding 4D seismic data in the history matching workflow.The workflow has been illustrated on some synthetic cases as well as on a data set from the Norne field.In addition, tests has been done on a data set from the Ekofisk field.The workflow consists of a number of different steps, and techniques has been developed for handling all of them.The Norne filed was successfully history matched, and it was possible to improve the fit to both seismic and production data simultaneously.However, the ongoing study on the Ekofisk field reveals that there are additional challenges that needs to be addressed.
Handling the complex problem of history matching 4D seismic data has inspired the developments of several tools described in this document.The concept of sparse data representation inspired a few works later, including the work of [30] (who used dictionary learning for sparse data representation), and that of [11] (who used convolutional neural network for sparse representation).For the concept of correlation-based localization, it has also been used to deal with a few other problems like CO2 monitoring, model error characterization through a machine learning model [20] and correlation-based local analysis [31].
The history matched model will be used for further planning of the reservoir development, including for important decisions as planning of new wells, decisions on using EOR based methods for improved reservoir drainage.A report on the effect of the research funded by Research Council of Norway (see [27,28]) reported about better decision support by utilizing ensemble-based methods which in certain cases could lead to faster reservoir drainage, or a reduction in the number of drilled wells.It is tempting to believe that adding more information into the process in terms of using 4D seismic data should give benefits, in terms of higher economic return, and some of the decisions might also reduce the emissions of greenhouse gases.
The presented workflow for 4D seismic history matching might not have reached its final form, and further testing on a range of different fields is required.This testing would be better done by the operators having access to all the data and information about the fields.There is reason to believe that this testing will bring up new research questions that have to be tackled.We find that handling seismic data is challenging, and there might be large systematic deviations between the actual seismic data and those obtained by simulations with the initial model.To handle this problem, we suggest to focus on matching key aspects of e.g.acoustic impedance, time-shift and / or time strain.
Part of the developed methodology is now planned to be or being implemented in different places.For instance, correlation-based localization was implemented in a software of USGS (https://pubs.er.usgs.gov/publication/tm7C26),and Equinor will implement correlation-based local analysis in their ERT (ensemble reservoir tool) system.

Figure 2 :
Figure 2: Procedure of sparse data representation through the discrete wavelet transform (DWT), and the related noise estimation problem in the wavelet domain.

Figure 3 :
Figure 3: Data mismatch for production data (top) and seismic data (bottom).The Y-axis is logarithmic.The red horizontal lines indicate the median, and the blue horizontal lines indicate the 25th and 75th percentiles.The whiskers indicate the most extreme values.[13]