Bayesian seismic 4D inversion for lithology and fluid prediction

Seismic inversion is a technique that aims to estimate the subsurface properties of the earth from seismic data. It can provide valuable information for reservoir characterization during exploration and production. Using repeated seismic surveys taken at different calendar times over the same area allows for 4D (also denoted time-lapse) seismic inversion. Dynamic changes in the subsurface due to fluid change, pressure depletion, or compaction caused by injection or production lead to changes in the seismic response. These changes can be detected by comparing the surveys (Landrø, 2015). Therefore, 4D seismic inversion has become a valuable tool for monitoring reservoirs during hydrocarbon production or for monitoring the flow of CO₂ during a storage operation (see the application overviews in Lumley [2001] and Johnston [2013]).

The inherent nonuniqueness of the seismic inversion problem, as well as the complexity and heterogeneity of geologic structures, means that there is large uncertainty in the results. This has led the way to probabilistic methods (e.g., Tarantola and Valette, 1982), wherein the nonuniqueness and uncertainty are addressed. The Bayesian approach, wherein a prior probability model is updated to a posterior probability distribution based on data, is common in solving geophysical inverse problems (Buland and Omre, 2003; Eidsvik et al., 2004; Larsen et al., 2006; Grana, 2016; Jullum and Kolbjørnsen, 2016; Kolbjørnsen et al., 2020). The Bayesian methodology captures the nonuniqueness well and constrains the results by carefully defining the prior model from geologic knowledge. However, the computational cost related to the inversion of large-scale spatially coupled models is challenging and even more so in a 4D setting.

In its simplest form, 4D seismic inversion is performed as a sequence of 3D inversions (e.g., Rosa et al., 2020). By comparing the inversion results of a repeat vintage with the inversion of a baseline survey, it is possible to infer the dynamic behavior of the reservoir and its impact on the elastic properties. However, inverting only single surveys does not take into account the correlations in the signal and noise between the vintages. In addition, having no coupling between the inversions leaves the inferred results vulnerable to inconsistencies.

An alternative is to do inversion directly on the difference in seismic data between two vintages. Because signal and noise are highly correlated between surveys, the difference between the two surveys highlights the changes between vintages. Inversions based on difference data were explored by Landrø (2001), Meadows (2001), and Trani et al. (2011) to extract changes in fluid saturation and pressure. Buland and Ouair (2006) apply Bayesian methodology to the difference in prestack amplitude data for two vintages to infer changes in elastic properties. In Grana and Mukerji (2014), the Bayesian method was applied to difference data to estimate static and dynamic time-lapse reservoir properties, whereas Côrte et al. (2023) invert for pressure and saturation changes in a thin reservoir zone using a Bayesian methodology on seismic difference data in a map setting.

Rosa et al. (2020) carry out a comparison of the outcomes of a Bayesian 4D difference inversion, four dimensions as a sequence of Bayesian 3D inversions, and four dimensions as a sequence of deterministic 3D inversions. It was shown that the Bayesian 4D difference inversion was the most accurate in terms of the misfit estimators chosen for the evaluation of algorithm performances. However, the inversion of difference data does not give direct access to the properties of the reservoir at the time of the baseline survey; only the 4D changes are addressed. Hence, a methodology that solves for the baseline reservoir properties as well as the changes between the vintages in a computationally feasible way is desired.

Forberg et al. (2021) perform the joint Bayesian inversion of two seismic amplitude-variation-with-offset (AVO) vintages, focusing on the evolution of reservoir properties (e.g., water saturation) in a dynamic setting of ongoing hydrocarbon production. They combine linear Gaussian seismic and rock-physics likelihood models with a selection Gaussian (Omre and Rimstad, 2021) prior model for the reservoir properties. This allowed the authors to use the Markov chain Monte Carlo methods to assess the selection of Gaussian posterior distribution. Using cross-vintage correlations in the multimodal prior model, consistent results for water saturation in both vintages were obtained.

In this paper, we perform simultaneous Bayesian inversion of two seismic vintages. This is similar to Forberg et al. (2021), but instead of inverting for reservoir properties, we invert for discrete lithology-fluid classes (LFCs) in both vintages. Gaussian elastic prior models per LFC, including cross-LFC and cross-vintage correlations, account for the multimodality in the elastic domain that connects the LFCs with the seismic response.

Our method is an extension of the approach by Kolbjørnsen et al. (2020). This paper presented an efficient methodology for the Bayesian prediction of lithology and pore fluid from prestack seismic data. The nonuniqueness was handled by evaluating all possible lithology and fluid combinations. Doing this for a full survey is computationally infeasible, so it was done trace by trace and only within a local neighborhood around each sample point. These local inversions were then combined into a consistent result for the complete trace. This approach makes it possible to avoid Monte Carlo sampling. In Ndingwan et al. (2022), the methodology was used to infer the changes in elastic properties from a 4D difference inversion.

We do not explicitly invert the 4D difference in seismic data but rather use coupled seismic noise and correlations between the elastic properties of the baseline and repeat vintages. This still requires the alignment of seismic surveys, and we assume that time shifts are corrected for in the input data (for an overview of 4D data acquisition and processing, see, e.g., Johnston, 2013). By simultaneously inverting two vintages, taking correlations into account, we ensure that although the exact elastic properties at each location are uncertain, they stay consistent between the vintages. Furthermore, we are able to handle correlated noise between the vintages. Separate inversions of the two vintages cannot do this.

The Bayesian framework allows us to incorporate prior information on the stratigraphic layering, lithology and fluid distribution, and layer thicknesses. For the baseline vintage, we specify the possible LFCs within each stratigraphic layer. Some baseline LFCs, e.g., hydrocarbon sand, are expected to change in the repeat vintage, whereas others, e.g., shales, remain the same. The changes are modeled by new discrete LFCs for the repeat vintage. Each baseline LFC can transition into one or more LFCs in the repeat vintage, and each such pair of LFCs is conveniently called a 4D LFC.

Typically, the possible cross-vintage changes are the effects of fluid flow, pressure changes, compaction, or temperature changes. Thus, the elastic properties of the repeat vintage will be strongly correlated with the elastic properties of the baseline vintage. In addition, seismic noise is highly correlated across vintages. We include prior information on elastic correlations as well as the cross-vintage covariance of seismic noise. These covariances between vintages are the key components as we here make a 4D extension of the model in Kolbjørnsen et al. (2020).

The posterior probability distribution for 4D LFCs is computed. We obtain posterior probabilities that are consistent across the vintages, and 4D changes are easily inferred from this. An additional advantage of our proposed method is that the results for the baseline will benefit not only from the baseline seismic data but also from both data vintages. This helps to constrain the baseline results. Results using the main methodology were previously presented in Kjønsberg et al. (2020, 2022). In this paper, we thoroughly describe the method, illustrate it with a synthetic case, and present results for a real case.

The approach that we use here is an extension of the 3D inversion presented in Kolbjørnsen et al. (2020). We will follow the same notation and organization and stress the differences introduced by the move to four dimensions. Our concern is to do a simultaneous inversion of two seismic vintages, predicting joint probabilities for discrete LFCs.

We will use boldface small letters for vectors and boldface capital letters for matrices. If referring to a single vintage, an index will indicate this. Whenever we consider one specific vertical time sample, this is indicated with a second index. For instance, $dv$ is the vector over the time samples and angles for the seismic data of vintage $v$ and $dv,i$ is the vector over the angles for the seismic data of vintage $v$ in location $i$⁠. If the $i$ sample of a vector is scalar, we show this by using plain notation. For instance, $fv,i$ is the LFC of vintage $v$ in location $i$⁠. The different symbols are explained consecutively in the text.

The matrices $Σm1|f1$ and $Σm2|f2$ in equation 12 are the elastic variance terms for each survey. The cross-covariance of elastic parameters for the two vintages $Σm1,m2|f1,f2$ is a new feature in this work and needs to be discussed further. First, we consider the case when there is no change to the rock between the two surveys. In this case, the correlation is one. Next, we consider the case of a fluid and/or pressure change. There will still be a high degree of correlation between the two distinct time steps because the dry rock matrix is the same. The effect of this is that the resulting joint covariance matrix is singular. However, this is a strength rather than a weakness because it effectively reduces the number of parameters we need to estimate. For example, when there is no change, we only need to estimate three rather than six parameters in each position because $m1,i=m2,i$ (see equation 4).

We now consider how to assess the joint distribution of the two vintages when there is a change. This can be derived by generating dual point clouds, that is, by generating elastic properties at the time of vintage 1 and computing the effect of fluid and pressure changes to provide the corresponding rock under the conditions of vintage 2. The sensitivity to fluid and pressure changes is a property of the lithology, whereas the variability of fluid and pressure changes is specific to a lithology and fluid class. Thus, multiple fluid classes can belong to one lithology, and we can compute the joint distribution of all these fluid classes for the lithology. The required mean and covariances then can be derived from these point clouds (for details, see Appendix  A).

We have a double set of LFCs at each location, one value for each of the two vintages. We use the term 4D LFC for this combination of LFCs in the two vintages. The definition of classes for vintage 1 is just as before: one class for each possible lithology and fluid combination. For the vintage 2 classes, we need the same classes as are used in vintage 1 (to represent the areas of no transition). In addition, we need extra classes for 4D transitions, representing the new (e.g., fluid) state in vintage 2. Multiple classes can be used to represent the different strengths of one transition, such as partially and fully depleted oil.

By using a one-step Markov chain, we can take into account knowledge about the order of LFCs within each stratigraphic layer for each vintage. This makes it possible to model, for instance, strict fluid ordering in the baseline vintage while allowing for fingering in the repeat vintage.

The seismic noise covariance matrix, equation 13, consists of vintage-specific noise covariance matrices $Σe1$ and $Σe2$ and the cross-vintage covariance matrix $Σe1,e2$⁠. Note that the use of cross-vintage noise covariance distinguishes our approach from that of difference inversions, wherein one of the key features is that correlated noise is removed. The vintage specific noise models can be estimated as before (e.g., Buland and Omre, 2003). As our approach has separate covariance matrices for noise in the two surveys, we can easily handle different noise levels. To specify the cross-vintage noise covariance, we need to determine the size of the covariances, represented by the diagonal of the matrix, and the cross-vintage correlations, represented by the off-diagonal structure. The details of how we have done this are presented in Appendix  B. Because the cross-vintage correlation is a new feature, we here provide an illustration of the principle. The likelihood equations of the 4D model are valid however one chooses to estimate the noise correlations.

To ensure consistency with the noise covariance matrices for the baseline and repeat vintages, we use the same approach to build the covariance structure for the cross-covariance, as detailed in Appendix  B.

The core assumption regarding the data that is new in this approach, compared with standard probabilistic approaches (e.g., Buland and Omre, 2003; Grana, 2016), is that the 4D vintages are aligned. We use a convolutional model with linearized weak contrast approximation (Aki and Richards, 1980) and colored Gaussian noise. Errors in the alignment between surveys are not directly included in the noise model as presented. The model allows for flexibility in noise levels and wavelets between vintages, unlike a standard difference model. Due to the use of a convolution model, our approach is best suited for time-domain data.

Data acquisition and processing (Johnston, 2013) may impact the results of our inversion algorithm in two ways. First, the processing can improve the correspondence between the data and our assumptions. Deviations from these assumptions will be interpreted as noise, giving a higher noise level and less precise results. Second, the amount of noise can be impacted, and, in particular, the amount of uncorrelated noise between the surveys. Reducing the level of noise in general and the uncorrelated noise in particular will improve the results.

We assume stationary rock-physics distributions. The main motivation for this is to ensure the computational feasibility of the model. However, the spatial variation in complex reservoirs is a relevant topic. In principle, our model is valid for nonstationary rock-physics distributions, such as those of Rimstad et al. (2012). Implementing the general nonstationarity of the mean and covariance within our framework would increase the run time of the inversion significantly, as this requires an inversion of the conditional seismic covariance matrix. Nonstationarity only in the mean values can still be accounted for by a simple shift in the mean response of the signal and is a useful extension of the model presented.

We have now established a framework for computing the posterior LFC distribution $p(f|d)$ as the product of a Gaussian likelihood $p(d|f)$ and a prior LFC probability $p(f)$⁠. Our likelihood is based on two key steps: First, we use the linearized forward model with Gaussian noise. We then combine this with a conditional Gaussian distribution for elastic parameters given LFC. This allows us to compute the posterior density $p(f|d)$ for any $f$ and $d$⁠, except for a normalizing constant.

Unfortunately, being able to compute the posterior theoretically is of limited use here. As $f$ is high dimensional and discrete, any integral-based inference is infeasible. Sampling is also very difficult, and finding the $f$ with the highest posterior density is difficult, even for the smaller problem of 3D inversion (see, e.g., Gunning and Sams, 2018). However, with the framework we have established, we can apply the local approach from Jullum and Kolbjørnsen (2016). The idea here is to map out the complete posterior probability for a local window of a few cells at a location. Because the number of cells is limited, we can map out all possible configurations. These local distributions are then combined to full trace predictions of LFC probabilities using a Markov chain. The details are given in Kolbjørnsen et al. (2020); the only modifications needed for use in the 4D setting are extending the likelihood computation to four dimensions and using 4D LFCs in the Markov chain.

To illustrate the methodology of this paper, we have created a synthetic test case for a gas injection scenario. The system consists of a sand wedge constrained above and below by two different shales. In vintage 1, the reservoir sand has a layer of gas above the oil. In vintage 2, some of the oil has been replaced by gas, thus moving the gas-oil contact (GOC) downward, and this is shown in Figure 2. The LFC true model is shown in the top row, and the seismic data generated from the LFC true model are shown in the subsequent rows. The GOC is at 2352 ms in vintage 1 and 2452 ms in vintage 2.

Given the LFC true model in Figure 2, seismic data are generated by using the forward model of equation 7. We sampled the elastic parameters of vintage 1 and then performed Gassmann fluid substitution to find elastic parameters for sand filled with injection gas for the vertical interval between 2352 and 2452 ms. Sampled elastic parameters for vintage 1, along with the fluid-substituted vintage 2, are shown for one trace in Figure 3. The elastic properties of the true model are identical across the vintages for the intervals without fluid substitution but different in the interval wherein the oil is replaced by injected gas. After convolving the weak contrast reflections of the elastic parameters with an appropriate wavelet, we add seismic noise that follows the stratigraphy and is correlated between the vintages.

Full details on the generation of the synthetic data are shown in Appendix  C.

A stratigraphic model with three layers is defined in Figure 4. The reservoir layer is separated from the layers above and below by two horizons, both with uncertainty. The layers above and below the reservoir consist of two different shales. In the reservoir, vintage 1 can have oil- or gas-filled sand. The LFC prior model for this vintage is modeled by specifying each LFC’s probability at the top of a given layer and the vertical transition probabilities among the possible LFCs. This lays the groundwork for a one-step Markov model, as detailed in Figure 5. We assume that gas and oil are ordered in vintage 1, meaning that the vertical transition probability from oil-sand to gas-sand is zero.

A local neighborhood approximation with a window length of three is used. The possible window configurations of our synthetic case are shown in Figure 6.

The possible changes from vintage 1 to vintage 2 are described by stationary cross-vintage transition matrices for discrete LFCs. In this example, we allow the oil in vintage 1 to be replaced by injected gas in vintage 2. For each cell in the trace, given that vintage 1 has oil, the probability that vintage 2 has instead injected gas is set to 0.5. The probability that no fluid change has happened is also set to 0.5. A priori, we thus assume no knowledge of whether fluid substitution has occurred. For the other LFCs of vintage 1, we assume that there is no change. Figure 7 summarizes the possible cross-vintage transition probabilities.

For vertical continuity in vintage 2, we set the free parameter in equation 17 to $r=0.95$⁠. Thus, we have specified all terms of the 4D one-step Markov chain in equation 18. The final 4D vertical transition matrices within each stratigraphic layer are shown in Figure 8. The strict ordering of fluids in vintage 1 prevents vintage 2 from having injected gas above native gas. However, the transition matrix from the 4D combination OilSand/OilSand to OilSand/GasInjectSand is nonzero. This means that the fingering of oil and injected gas is possible in vintage 2.

In the time intervals with horizon uncertainty, the 4D vertical transition matrices are modified relative to what is shown in Figure 8. For instance, in a small interval around the top reservoir horizon, there will be a nonzero probability of going from Shale 1/Shale 1 to any of the reservoir 4D LFCs. This is similar to the method for a 3D inversion (Kolbjørnsen et al., 2020) and is not covered in more detail here.

The Gaussian prior model for elastic parameters per LFC is shown in Figure 9a. The injected gas fluid represents a mixture of pure oil and gas, and it has properties similar to, but not identical to, pure gas. We use uniform (homogeneous) saturation changes. The rock-physics distributions shown in Figure 9 are in accordance with the Gassmann fluid substitution used to define the true model data.

Covariance between the elastic parameters in the two vintages is crucial to take advantage of the simultaneous inversion of both vintages (equation 12). The cross-LFC correlations for the possible fluid changes make important contributions to this covariance. The prior elastic correlations for oil-sand and sand with injected gas are shown in Figure 10. When there is no change in the LFC from vintage 1 to vintage 2, the elastic correlations across the vintages are identical to the correlations of the LFC itself. In particular, the diagonal elements of these correlation matrices are one.

Gaussian seismic noise is added. The standard deviations for the noise are identical to what was used for creating the seismic data. The intervintage noise correlation is 0.9 (equation 19). This scales the cross-vintage covariance matrix. No lateral correlations are used for the prior noise model for inversion. The seismic noise is colored with the same wavelet that was used for data generation.

The results of the simultaneous inversion are shown in Figure 11, with the prior probabilities for each 4D LFC in Figure 11a and the posterior probabilities in Figure 11b. The prior probabilities are determined by the 4D top probabilities and vertical transition probabilities, as summarized in Figures 8 and 12. As is standard in single-vintage inversion, we have included a prior probability of 0.05 for an unknown LFC, representing the possibility that none of the permissible LFC configurations can explain the seismic data. Hence, the priors of Figure 11 sum to 0.95 in each grid cell.

Focusing on the reservoir layer, we observe that the posterior results find that above approximately 2350 ms, both vintages have sand filled with native gas. Between 2350 and 2450 ms, vintage 1 has oil, whereas vintage 2 has sand filled with injected gas. Below 2350 ms, both vintages have oil-filled sand. Detection of the GOC in the true model in vintage 1 (2352 ms) and vintage 2 (2452 ms) is thus excellent in this simple synthetic case.

Figure 13 shows the results for three different traces, focusing on the reservoir sands. The chosen traces are shown in Figure 14. Inline 10 (Figure 13a) corresponds to a trace where both vintages’ GOCs are within the reservoir. For inline 40 (Figure 13b), the vintage 2 GOC at 2452 ms is not present in the reservoir; thus, the interval where we expect fluid substitution is located directly on top of the base layer consisting of shale 2. Inline 60 corresponds to a case wherein neither vintages’ GOCs are present in the reservoir layer; thus, the interval where no fluid substitution is expected rests directly on top of the base layer. In all three cases, the detection of fluids in both vintages is as expected for the simultaneous inversion (columns two and three), with posterior probabilities close to one.

In the “Methodology” section, we have stressed that the elastic correlations between the vintages (equation 12), as well as the correlated seismic noise (equations 13 and B-2), are essential for the simultaneous inversion of two vintages. To illustrate the effect of these correlations, we can compare the 4D posterior results to what is achieved when these correlations are set to zero. The results of zero-correlated inversions are shown in the two rightmost columns of Figure 13, and they are denoted as “independent.” We have set up two prior models, one for vintage 1 and one for vintage 2, and for each model, the LFC prior is identical to the corresponding vintage prior for the 4D model. The independent vintage 1 model allows for oil and native gas, whereas 2 model also includes the possibility of injected gas. Both independent models have a strict ordering of oil and native gas, but the independent vintage 2 prior model allows for oil above the injected gas. The covariances of the elastic models are identical to the diagonal submatrices of equation 12 and similarly for seismic noise.

The results clearly illustrate that uncorrelated inversions, wherein the individual vintage posterior probabilities are conditioned only on the seismic data of the given vintage, are less able to detect the correct GOCs. Because the rock physics of gas sand and sand with injected gas are quite similar (Figure 9), the independent vintage 2 inversion has problems distinguishing them. In the rightmost column of Figure 13, we have therefore included the sum of the pointwise probabilities for these LFCs. The results are then comparable to the 4D vintage 2 results for the traces without a vintage 2 GOC present in the reservoir (inlines 40 and 60). For trace 10, the vintage 2 GOC at 2452 ms is suggested but not clear.

Comparing the posterior results for the 4D and independent data, it is clear that in the 4D inversion, each vintage also receives information from the other vintage’s seismic data. In particular, the posterior results for vintage 1 are also affected by vintage 2 seismic data. This is, of course, not the case for the independent results, wherein, for instance, only vintage 1 seismic data have an impact on the posterior results for this vintage.

To further show the behavior of the 4D inversion algorithm, it has been applied to seismic data from the Edvard Grieg oil field. This field is located in the Utsira High area in the central North Sea, approximately 35 km south of the Grane and Balder fields and to the west of the Johan Sverdrup field. Its location is shown in Figure 15. The field is operated by Aker BP, with partners OMV and Wintershall DEA. Following production startup in 2015, the field is produced by pressure support from water injection. It produces undersaturated oil from the shallow marine sandstones of the Cretaceous period, as well as the aeolian and fluvial sandstones and conglomerates of the Permian or Triassic period. The reservoir quality varies from moderate to very good in the sandstones, whereas the quality is relatively poor in the conglomerates. Oil is also proven in the altered basement in the north of the field. The top of the reservoir is encountered at a depth of approximately 1900 m below sea level. Details on the sedimentation of Edvard Grieg and reservoir characterization by means of Bayesian 3D inversion can be found in Ndingwan et al. (2018).

The 4D inversion presented here uses ocean-bottom seismic surveys acquired in 2016 (vintage 1) and 2018 (vintage 2). The seismic data are given as near- (0°–10°), mid- (10°–20°), and far- (20°–30°) angle stacks to capture AVO effects (Ndingwan et al., 2022).

Ten different lithology classes are defined to use in the inversion. There are two different shales in the overburden (Lista and Base Lista) and two types of chalk (acoustically fast and slow) between the overburden and the reservoir. In the reservoir, there are aeolian sandstones (Aeolian) and fluvial/aeolian reworked sandstones (FAR), a sandy conglomerate (Sandy-CGL), and a shale (FloodPlain-Shales). The underburden consists of an altered or solid basement. The reservoir sands and conglomerates may have different fluid fillings, represented by categorization into three distinct fluid filling classes. These represent saturated oil filling (Aeolian-85OIL, FAR-85OIL, and Sandy-CGL-60-OIL), depleted oil filling (Aeolian-15OIL, FAR-15OIL, and Sandy-CGL-20-OIL), and 100% water filling (Aeolian-100WTR, FAR-100WTR, and Sandy-CGL-100WTR). We use homogeneous mixing in Aeolian and FAR, whereas in the conglomerate model, we use a patchy model (Brie et al., 1995). Altogether, we have 16 LFCs.

The rock-physics model consists of Gaussian models for the P- and S-wave velocities, $VP$ and V_S, and for density $ρ$⁠. Figure 16 shows acoustic impedance (AI) and the $VP$/$VS$ ratio for each LFC. The Aeolian and FAR sandstones are known to have excellent reservoir properties. They have comparable properties for $VP$/$VS$ and AI but separate well from the sandy conglomerate (Sandy-CGL). The conglomerate’s elastic properties overlap some with the chalks; however, they occur in different stratigraphic zones.

Each of the sandstones and conglomerate is represented by three different LFCs. For a given lithology, the LFCs differ only in fluid filling, whereas the rock frame is unchanged. This means the elastic properties are strongly correlated. We use pairwise cross-LFC correlations for the following groups:

• Aeolian-85OIL, Aeolian-15OIL, and Aeolian-100WTR.

• FAR-85OIL, FAR-15OIL, and FAR-100WTR

• Sandy-CGL-60-OIL, Sandy-CGL-20-OIL, and Sandy-CGL-100WTR.

There is an almost perfect correlation between $VS$ and density across the fluid filling. Crosscorrelations across the three groups are zero, reflecting that different geologic histories for three lithologies make their elastic properties uncorrelated. An almost complete overlap between the elastic properties of the depleted Aeolian (Aeolian-15OIL) and the oil-saturated FAR (FAR-85OIL) makes them difficult to distinguish.

We divide the inversion interval into five different zones, which is schematically shown in Figure 17. Top Chalk and Top Basement are seismically interpreted horizons, whereas Top Cromer Knoll has been derived from Top Chalk because it is in tuning. We also explicitly use the oil-water contact (OWC) because its location is well known in the Edvard Grieg field. We allow for some vertical uncertainty for each horizon defining a zone boundary.

The overburden, chalk, and basement zones all consist of LFCs that remain unchanged from vintage 1 to vintage 2.

The reservoir is divided into two zones, separated by the OWC. It consists of sandstones, conglomerates, and shale. Above the OWC, the only fluid filling allowed in vintage 1 is saturated oil, with Aeolian-85OIL dominating. In vintage 2, the sands and conglomerates may remain oil saturated, or they may be substituted by their depleted version. The only possible LFC change from vintage 1 to vintage 2 is that an oil-filled LFC can be replaced by its depleted counterpart, e.g., Aeolian-85OIL with Aeolian-15OIL. We set the prior transition probability of change to 50% (see Figure 18). This is used as the reference level when we evaluate the posterior to look for evidence for or against change. Below the OWC, we find the water-filled sands and conglomerate as well as the shale called FloodPlain-Shales. The conglomerate is assumed to be dominating here. No LFC transitions between the vintages are allowed in the lower part of the reservoir.

The cross-vintage dependency factor $ρs$ in equation A-3 is set to 0.9. This increases the possibility for change in the repeat vintage slightly compared with a default assumption of 1. To ensure vertical continuity in vintage 2, we set the free parameter $r$ in equation 17 to $r=0.95$⁠, representing a high probability that an LFC cross-vintage transition at one location will continue as long as the same LFC is present in vintage 1. A high value of $r$ prevents unwanted LFC oscillations in the repeat vintage.

We focus on fluid substitution in the reservoir by investigating depleted LFC posterior probabilities in vintage 2. These are equivalent to the difference between vintages 1 and 2 of the oil-saturated LFC posterior probabilities. For a more concise presentation, we will treat the three oil-saturated LFCs (Aeolian-85OIL, FAR-85OIL, and Sandy-CGL-60-OIL) as a single unit in the sense that we only investigate their summed probabilities. The same is done for the three depleted LFCs (Aeolian-15OIL, FAR-15OIL, and Sandy-CGL-20-OIL).

In Figure 19, we look at inversion results for a single cross section (inline 460). The posterior probabilities for the oil-saturated LFCs in vintage 1 are shown in Figure 19a, the posterior probabilities for the depleted LFCs in vintage 2 are shown in Figure 19b, and the difference between the posterior and prior probabilities for the depleted LFCs in vintage 2 given oil-filled LFCs in vintage 1, scaled by the posterior probabilities for oil-filled LFCs in vintage 1, are shown in Figure 19c. Plots of the full seismic stack of vintage 1 and the corresponding seismic difference between the vintages are shown in Figure 19d and 19e, respectively. Figure 19 indicates that significant posterior probabilities for the depleted LFCs in vintage 2 coincide with large seismic differences.

The same arrangement of figures from a cross section (inline 888) in an area without any significant seismic differences is shown in Figure 20 (inline 888). Here, we can see some indication of the structures in the seismic difference, possibly due to misalignment, but nothing as distinctive as shown in Figure 19. The amount of depletion is much more uniform across the section.

A survey-wide summary of the inversion results is shown as a map in Figure 21. Figure 21a shows the normalized root mean square (nrms) of the seismic difference in the upper reservoir zone. Note that this zone is not strictly limited by the Top Cromer Knoll/OWC horizons because they are defined with some degree of uncertainty. Figure 21b shows the maximum difference between the conditional posterior and prior probabilities of depleted LFCs in vintage 2, given the oil-filled LFCs in vintage 1 scaled by the posterior probabilities of the oil-filled LFCs in vintage 1. The areas that indicate the largest increase in posterior probabilities compared with the prior probabilities in Figure 21b are clearly colocated with the areas that have significant seismic differences in Figure 21a.

Figure 22 shows a more detailed presentation of a single trace from the cross section in Figure 19, taken at a crossline, which shows a significant seismic difference. The seismic difference of the three angle stacks is shown in Figure 22a, and the posterior probability of oil-saturated LFCs in both vintages is shown in Figure 22b, together with the difference between the conditional posterior and prior probabilities of the depleted LFCs in vintage 2, given the oil-filled LFCs in vintage 1 scaled by the posterior probabilities of the oil-filled LFCs in vintage 1. We can see that a higher degree of substitution than the prior takes place within the interval indicated by the dashed black lines, corresponding to the interval with a significant seismic difference.

The Bayesian method presented in this paper simultaneously solves for the LFC probabilities in two vintages. Here, we discuss four main aspects of the method: the inherent intervintage consistency, the use of discrete LFCs, the role of the prior model, and the computational complexity.

By simultaneously inverting for vintage 1 and vintage 2 posterior LFC probabilities, we ensure that the results are consistent across the vintages. If, for instance, the only permitted change is to replace pure oil with depleted oil, vintage 1 must have pure oil if vintage 2 has depleted oil. Likewise, a sample with shale in vintage 1 must also have shale for vintage 2. This applies to any given vertical time sample. Consistency is an inherent part of our approach, and it cannot be ensured by independent inversions for the two vintages.

An implication of the consistency is that the simultaneous inversion can provide improved results for vintage 1 because both vintages’ seismic data affect the results. We see this clearly in Figure 13. A difference inversion, wherein the difference in seismic data is inverted to infer changes in the earth profile, does not provide explicit information on vintage 1 but does give internally consistent results for the 4D changes.

Another implication is that there may be a competition between the individual vintages’ seismic fits and the difference fit. If the vintage 2 seismic data strongly favor the depleted oil, this pushes vintage 1 toward having pure oil, even if the vintage 1 data alone might fit better with a different LFC. Careful specification of the prior model is therefore advisable. The permitted cross-vintage LFC changes, the rock-physics distributions per LFC, and the cross-vintage covariances are all essential as they enable physical constraints. Geophysical knowledge, including LFC and the rock-physics extraction from the wells and production data, is useful. The beneficial use of this kind of information is not unique to our approach, but sensitivity to the information content of the prior model is a cost the presented methodology pays to ensure consistency. The reduction of this sensitivity is a relevant research topic for the future.

Our model uses discrete LFCs for both vintages. Multiple classes can be used to represent the different strengths of the transitions from the baseline to the repeat vintage. Examples of repeat vintage LFCs are partially or fully depleted oil or an LFC for which the pressure effects have changed the elastic properties. Our model choice allows us to explore different fluids and/or pressure change regimes, but the correct prior specification and interpretation of the posterior probabilities can be challenging. A key reason for this is that in the areas affected by production, fluid states no longer fall nicely into discrete classes. The discretization becomes an approximation, and the probability for transition LFCs will represent the degree of reservoir changes. Using several transition LFCs with different degrees of fluid substitution gives a less coarse discretization. There are two challenges related to adding more transition LFCs to the model. Increasing the number of classes will increase the amount of overlap between them, making LFC separation more difficult. In addition, it leads to increased computing time because more LFC configurations need to be considered. Both of these aspects will limit the practical number of transition LFCs, but there is no simple way of determining the optimal number.

The dynamic 4D changes are functions of the static reservoir properties, i.e., a stiff reservoir sand will have smaller fluid and pressure sensitivity than a soft reservoir. Hence, the LFC changes and natural variability within each LFC can affect the 4D signatures. Thus, stiff reservoir sand has little to gain from using more than one transition LFC, whereas this might be beneficial for soft sand.

In our results, there is a tendency to predict the probability of change even in areas where there are no clear 4D effects. The drivers for this can be inferred from equation 5. The prior probabilities $p(f)$ for 4D LFCs with transitions are positive everywhere and are the root cause. This is a reasonable prior, as we know that transitions will occur, and we do not want to use the prior to rule out a transition at a location where we have predicted oil. We had hoped that the posterior would still be able to predict close to zero transition probabilities in areas with no 4D effects. However, we see that with our settings, there is a positive likelihood for transitions even in these areas. There are two possible reasons for this: either the evidence against the transition is not sufficiently strong, or our cross-vintage correlations are not calibrated precisely enough to discern this evidence. We managed to reduce the transition probability relative to the prior in these areas, but we did not manage to push it to zero, as we would have wanted.

When doing difference inversions in four dimensions, this problem is often handled by defining a prior, which is symmetric around zero. This way, no signal corresponds to no change, which is intuitive. However, note that our result is stronger; we claim that there is evidence against a change because our posterior probability is reduced compared with the prior distribution. In a Bayesian approach, we should always quantify our learnings by comparing them with the prior distribution. A key advantage of the pure difference seismic data is that the covariance structures are easier to quantify. The benefit of using baseline and repeat vintage seismic data, which is our approach, is that we have the opportunity to gain information not only about change but also about the initial and final states at the cost of specifying a more complex model.

The computational complexity of the Bayesian seismic inversion algorithm presented in Kolbjørnsen et al. (2020) is given as $O((Nnl2+nl3)nm)$⁠, where $nl$ is the wavelet length, $nm$ is the number of permissible local-neighborhood LFC configurations, and $N$ is the number of locations in which the result is computed. The computational complexity is increased in a joint 4D inversion for two reasons:

• Here, $nm$ increases because the number of permissible 4D LFCs is typically larger than the number of basic LFCs. The magnitude of the increase depends on the fraction of LFCs with nonzero cross-vintage transition probabilities to two or more LFCs. A worst-case approximation for a local neighborhood of $wl$ samples with no vertical ordering is $Mwl$ if we assume that any LFC in vintage 1 can transition to $M$ different LFCs in vintage 2.

• The joint likelihood calculations are performed on data of twice the length because of the two vintages involved, so the effect of the wavelet length $nl$ is doubled.

This leads to a computational complexity of $O((Nnl′2+nl′3)nm′)≤O(4(Nnl2+2nl3)Mwlnm)$ for $nl′=2nl$ and $nm′≤Mwlnm$⁠. The 4D inversion algorithm has an approximate worst-case increase in complexity by a factor $4Mwl$ compared with its independent-vintage counterpart, but only a limited amount of LFCs will experience 4D transitions in practice, so the actual increase could be significantly lower.

We have presented a Bayesian seismic inversion method to find the joint posterior probabilities for discrete LFCs in two vintages of 4D seismic data. The method is based on the assumption that elastic parameters are different, yet highly correlated, between vintages due to 4D effects. The use of discrete LFCs and a prior model makes it possible to allow only certain 4D effects, such as particular forms of fluid substitution in lithologies. The increased computational complexity from doing a joint inversion is typically small enough to make it feasible for practical applications. The synthetic results demonstrate that the algorithm is able to use dual-vintage data together with a prior model specifying their correlations to calculate the joint LFC posterior probabilities for both vintages with a lower degree of uncertainty than independent single-vintage inversions. Analysis of the data from the Edvard Grieg field further indicates that the underlying model is sufficiently general to explain the 4D variations in seismic data using a reasonably simple prior model of 4D LFC changes. This is in line with previous work by others on 4D difference inversions but with the added value of being able to estimate the posterior LFC probabilities for each individual vintage and not just the probabilities of the relative changes between them.

We thank the sponsors of the GIG consortium: AkerBP, ConocoPhillips, Equinor, OMV, Sval, Total Energies, Vår Energi, and Wintershall DEA for financing the development of the algorithm and the work involved in the inversions presented here. We also thank AkerBP and their partners on the Edvard Grieg field, OMV and Wintershall DEA, for permission to present the data and publish the inversion results.

Data associated with this research are confidential and cannot be released.

The joint correlation of elastic parameters is a part of the prior in the proposed methodology, and it needs to be specified. Quantification of this property can be done under different sets of assumptions. Next, we show one way to define a consistent covariance structure based on pointwise covariance and a common underlying vertical time correlation. In addition to the time dependency and pointwise correlations, our method requires one additional parameter that can be interpreted as a cross-vintage dependency factor.

We focus on one single trace and let $A,X$ be an LFC vintage pair at time sample $i$⁠, while $B,Y$ is the pair at sample $j$ (see Figure A-1). For simplicity, the notation $Ak$ denotes the elastic parameter $k$ of LFC $A$⁠. For an isotropic medium, $k$ refers to $VP$⁠, $VS$⁠, or density. The pointwise correlations, valid for $j=i$⁠, between the elastic parameters $Ak$ and $Bl$ in vintage 1 are denoted by $ρ*(Ak,Bl)$⁠, and similarly for cross-vintage point-wise correlations $ρ*(Ak,Xl)$ and $ρ*(Bk,Yl)$⁠. The pointwise correlations and a vertical correlation structure $ρ(|ti−tj|)$ connecting time samples $i$ and $j$ are assumed to be known. In the following, we consider the case $j≠i$⁠. For simplicity, we denote vintage 1 as the baseline and vintage 2 as the repeat.

Here, $δ$ denotes the Kroenecker delta, and the scale factor $ρs∈[0,1]$ is fixed and can be interpreted as a cross-vintage dependency factor. It lowers the correlation when there is no LFC change in sample j. This allows for some flexibility in the repeat vintage. The known pointwise correlations are denoted as $ρ*$⁠, as described previously.

The middle section of equation A-5 represents the case in which there is an LFC change in one location but not the other. We then use the cross-vintage correlation in equation A-3 between the nonchanging LFC and the repeat LFC of the changed pair but lower it by scaling with the pointwise correlation for the LFC that does not change.

The last section of equation A-5 represents all the other possibilities, and we use an average of the pointwise correlations multiplied by the time-correlated cross-vintage correlations.

We emphasize that the quantification of elastic correlations can be done by different sets of assumptions. The expressions provided here show one possible way of setting up a consistent covariance structure and represent the method we have used.

An estimation of the seismic noise covariance matrix can be done in many different ways. This appendix illustrates the method used for the results in this paper.

The elastic properties of the synthetic test case’s true model are generated by first sampling the elastic properties of a system wherein the whole reservoir is oil filled, with shale 1 and shale 2 above and below, respectively (for the stratigraphic layout, see Figure 2). We sample the elastic parameters for this LFC configuration from a multinormal Gaussian distribution, with a pointwise mean and covariance structure according to the rock-physics prior model shown in Figure 10. The vertical and lateral correlations are included by means of a variogram with a lateral extent of 20 grid cells and a vertical extent of 20 ms. This gives elastic properties that are correlated vertically and laterally for the whole reservoir.

The reflection coefficients of the elastic properties shown in Figure C-1 are convolved with wavelets using the forward model defined in Kolbjørnsen et al. (2020). The same wavelet is used for all angles for both vintages (see Figure C-2).

We add Gaussian seismic noise. This is done by combining the 4D noise covariance matrix with a vintage-specific exponential variogram. The variogram has a lateral range of four grid cells and a negligible temporal range (0.4 ms). The 4D noise covariance matrix consists of vintage-specific noise matrices and a cross-vintage matrix. For each vintage-specific matrix, the colored noise correlation structure is computed from the wavelets, and the nugget ratio is 0.1 (the white noise component). The cross-covariance between the seismic noise in the two vintages is colored by the wavelet and scaled with a factor of 0.95 for numerical stability. The result is Gaussian seismic noise, correlated between vintage 1 and vintage 2 (see Figure C-3). The noise standard deviations are provided in Table C-1.

Biographies and photographs of the authors are not available.