# CO2 monitoring at Sleipner field using reflection oriented full waveform inversion: part 2 - 4D investigation

Fandy Adji Fachtony, Romain Brossier, Bhargav Boddupalli, Bastien Dupuy, Ludovic Métivier, Anouar Romdhane

<figure><img src="image_1.png" /></figure>

## To cite this version:

Fandy Adji Fachtony, Romain Brossier, Bhargav Boddupalli, Bastien Dupuy, Ludovic Métivier, et al.. CO2 monitoring at Sleipner field using reflection oriented full waveform inversion: part 2 - 4D investigation. Geophysical Journal International, 2025, 242 (1), pp.ggaf185. 10.1093/gji/ggaf185. hal-05336518

HAL Id: hal-05336518

https://hal.science/hal-05336518v1

Submitted on 29 Oct 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# CO₂ monitoring at Sleipner field using reflection oriented full waveform inversion: part 2—4D investigation

F. A. Fachtony,<fnref n="1" /> R. Brossier,<fnref n="1" /> B. Boddupalli,<fnref n="2" /> B. Dupuy,<fnref n="2" /> L. Métivier<fnref n="1" />,<fnref n="3" />
and A. Romdhane<fnref n="2" />,*

Accepted 2025 May 17. Received 2025 May 12; in original form 2024 October 21

## SUMMARY

In this work, we apply a dedicated 4-D full waveform inversion workflow to short offset streamer data from the Sleipner CO₂ storage field in the North Sea. We consider a baseline data set acquired in 1994 and a monitor data set acquired in 2008. Accessing to only short offset data raises significant difficulties for full waveform inversion. In this case the penetration of diving waves, which controls the depth where quantitative updates of the velocity can be expected, do not reach the zone of interest where the CO₂ is injected. For this reason, we propose to combine an efficient time-lapse full waveform inversion strategy, which we call simultaneous, with a reflection oriented full waveform inversion workflow. The latter has been introduced in the literature as a way to circumvent short-offset limitation and increase the ability of full waveform inversion to update the macro-velocity model at depth by exploiting the reflection paths, using a prior step of impedance reconstruction. We first illustrate the interest of this combined strategy on a 2-D synthetic model inspired from the Sleipner area. Then we apply it to the Sleipner field data, using as baseline model the one we present in a companion paper, where our reflection oriented workflow is presented. Our combined approach yields reliable estimates of the changes due to the CO₂ injection, characterized by velocity reductions of up to 400 m s⁻<fnref n="1" /> and strong impedance contrasts at depths of 800–1000 m, which consistent with previous full waveform inversion studies. Furthermore, the spatial distribution of CO₂ changes aligns with conventional seismic time-migration results from earlier studies, following a north–south migration trend.

**Key words:** Tomography; Computational seismology; Controlled source seismology; Crustal imaging; Waveform inversion.

## 1 INTRODUCTION

Carbon capture storage is considered by the Intergovernmental Panel on Climate Change (IPCC) as an effective means to reduce the concentration of CO₂ in the atmosphere and to mitigate the effects of the global climate changes. The process includes separation of CO₂ from industrial sources followed by the injection inside geological formation for long-term isolation. Despite its recent development, the current existing Carbon Capture and Storage (CCS) facilities around the world are still not sufficient to meet the target of 1000 Mt injection CO₂ per year for the Net Zero Emission scenario by 2030. In 2024, the deployment of CCS facilities are increasing all over the world, with the United States and

Europe being the leading region, while Asia-Pacific region showing their initial developments. Among existing CCS facilities, for almost three decades of operation, Sleipner storage facility in the north sea of Norway has been injecting 1 Mt CO₂ yr⁻<fnref n="1" /> since 1996. It is also one of the pioneering CO₂ offshore storage at the industrial scale. At the Sleipner field, CO₂ has been transported from the producing natural gas field from Sleipner west field, and injected into the Utsira formation (Eiken 2019) at the Sleipner east field.

Several 3-D seismic data sets have been acquired over the years to monitor CO₂ changes in the Sleipner field. For this purpose, conventional seismic monitoring methods have been adopted (Furre & Eiken 2014), where the quantification of time-lapse changes is performed in the image domain, involving pre-stack time migrated images of the baseline and monitor (Lumley 2001). Both

<footnote n="1">¹ISTerre, Univ. Grenoble Alpes, F-38000 Grenoble, France. E-mail: fandy.fachtony@univ-grenoble-alpes.fr</footnote>

<footnote n="2">²SINTEF Industry, Applied Geosciences Department, S.P. Andersens veg 15A, 7031 Trondheim, Norway</footnote>

<footnote n="3">³Univ. Grenoble Alpes, CNRS, LJK, F-38000 Grenoble, France</footnote>

<footnote>*Now at: Equinor, Norway.
© The Author(s) 2025. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.</footnote>
baseline and monitor data sets are migrated with the same velocity model, assuming a linear relationship between them. Prior to the 4-D analysis, a time-shift correction is applied to the monitor image to match the interfaces position on the baseline data. However, in case of strong changes of the rock properties between the monitor and baseline surveys, this linear assumption is violated. In addition, this workflow operates in the time-domain, and additional time-to-depth conversion on the data set could introduce errors due to conversion effects. Therefore, more robust, quantitative, monitoring analysis in the depth-domain is necessary to give more reliable time-lapse changes model.

Full waveform inversion (FWI) is a powerful method to estimate quantitatively subsurface parameters (Lailly 1983; Tarantola 1984, 1986), it has now become a standard practice both in industry and academia, with applications at various scales, from global (Thras tarson et al. 2024) and regional seismology (Fichtner & Villaseñor 2015; Davy et al. 2017; Górszczyk et al. 2021), to exploration (Operto et al. 2015; Hicks et al. 2016), near surface characteriza- tion (Provenzano et al. 2018; Irnaka et al. 2022), and even medical imaging (Guasch et al. 2020). The capability of FWI to take into ac- count complex propagation effects makes it a suitable technique to estimate the subsurface in a comprehensive and effective manner. However, it also results in FWI being a highly nonlinear inver- sion method, therefore a careful treatment of the data is necessary to produce reliable models. For time-lapse application, the use of FWI could be advantageous as it should be able to estimate di- rectly the mechanical properties from the baseline and monitor data set, therefore making the quantification of the time-lapse changes more straightforward, from the differences between the monitor and baseline models. Nonetheless, the inversion of two different data sets can lead to convergence issues due to the nonlinearity of the inversion, which may be responsible for time-lapse artifacts. The non-repeatability of the acquisition geometry between baseline and monitor data is another source of such artifacts. For this rea- son, a proper time-lapse strategy needs to be considered to handle potential drawbacks on the 4-D FWI applications.

In recent years, the application of time lapse FWI to field data has been rather limited. For most of them, highly repeatable surveys have been used. For instance, Hicks *et al.* (2016) applied time lapse FWI on repeated permanent-reservoir-monitoring system with long offset availability, good azimuth coverage and low frequencies, at Grane field. In this case, they were able to show that time-lapse FWI can reliably recover *P*-wave velocity changes related to gas replacing oil in hydrocarbon production. Similarly, Yang *et al.* (2016a) applied time lapse FWI on ocean bottom cable (OBC) data from Valhall field with mild non-repeatability problems between baseline and monitor. A double difference waveform inversion was applied and produced credible time lapse *P*-wave velocity changes.

At the Sleipner field, several studies on the application of time
lapse FWI have been already proposed in the literature, from 2-D
(Queißer & Singh 2013; Romdhane & Querendez 2014) to 3-D
applications (Raknes et al. 2015a; Mispel et al. 2019). Romdhane
& Querendez (2014) and Mispel et al. (2019) mainly work with
the monitor data to characterize the CO₂ anomaly changes, while
Queißer & Singh (2013) and Raknes et al. (2015a) are working
with both baseline and monitor data. Queißer & Singh (2013) and
Mispel et al. (2019) were able to get access to relatively long off-
set data, with maximum offset of 3600 and 6000 m, respectively,
which helps significantly to constrain the CO₂ anomaly at the deeper
part of the subsurface. Romdhane & Querendez (2014) and Rak-
nes et al. (2015a) work with a limited offset data (maximum offset

less than 1700 m), making the estimation of CO₂ related changes
challenging and less constrained. Because of this constraint in off-
sets, the information carried out by the diving waves sampling the
model at the depth of the reservoir is absent, making it difficult to
obtain reliable *P*-wave velocity changes. Instead, only reflections
sample this zone, making FWI behave similarly as a least-squares
migration algorithm, mainly adding high wavenumber information
to the model. One possibility to circumvent this difficulty while
relying on finite-frequency full waveform modelling approaches is
to use Reflection oriented Waveform Inversion, by exploiting the
tomographic contribution of the deep reflection wavepaths with a
prior reflectivity model. Building on this idea, Zhou *et al.* (2015)
propose to jointly invert the diving and reflection waves to pro-
mote tomographic updates both at the shallow and deep parts of the
model.

In a companion paper (Fachtony et al. 2025), we have applied reflection oriented FWI to the limited-offset baseline Sleipner data, showing that conventional FWI fails to yield significant model updates at the target depth due to a lack of diving waves penetration. We thus use the reflection oriented workflow FWI, proposed by Provenzano et al. (2023), which consists of impedance model building with an asymptotic preconditioner for fast and balanced reconstruction, reformulation of Joint FWI (JFWI) gradient in a pseudotime domain (Plessix 2013), ensuring velocity-reflectivity consistency with a velocity-impedance parametrization, and a graph space optimal transport objective function to mitigate cycle skipping (Métivier et al. 2019). We show that this reflection-oriented FWI provides more constraints within the target region, leading to a more accurate and reliable baseline model. In this study, this improved baseline model serves as the starting point for an investigation on time-lapse FWI strategies. We attempt to reconstruct time lapse changes for CO₂ monitoring from 3-D short offset data with this reflection oriented workflow combined with a time-lapse strategy which simultaneously invert the baseline and the monitor data set. We first perform a synthetic study on a 2-D Sleipner model to illustrate the benefits of this simultaneous 4-D FWI strategy over more conventional ones, and how it can be efficiently combine with JFWI. Within a long offset scenario we show that the simultaneous strategy provides more accurate time-lapse models and is more resilient to data noise. Within a short offset scenario, we illustrate how combining this simultaneous strategy and JFWI, we are able to better constrain the CO₂ distribution compared with conventional FWI. Besides, a subsequent FWI run starting from the reconstructed JFWI models effectively improves the reconstruction of the CO₂ changes. In a second part, we apply this workflow to the 3-D Sleipner monitor and baseline data. In this field data application, we show that our combined workflow is able to reconstruct a meaningful low-frequency Vp model and a high-resolution impedance model of the baseline and monitor, and provides reliable time-lapse model related to the CO₂ changes.

The paper is organized as follows: In Section 2 we briefly de-
scribe the time-lapse FWI, we review previously proposed time-
lapse strategies and we propose a time-lapse strategy applicable to
the reflection oriented workflow. In Section 3, we compare selected
time-lapse strategies on synthetic tests based on a 2-D Sleipner
model. In Section 4, we propose a workflow combining JFWI and
the simultaneous time-lapse strategy and illustrate its performance
on further synthetic tests on the 2-D Sleipner model. In Section 5,
we apply this workflow to the 3-D Sleipner monitor and baseline
data. The discussions and conclusions are given in Sections 6 and
7, respectively.
<figure><img src="image_3.png" /><figcaption>Figure 1. Workflow of time lapse strategies: (a) parallel, (b) sequential, (c) double difference, (d) central difference and (e) joint time lapse.</figcaption></figure>

## 2 REVIEW OF TIME LAPSE FULL WAVEFORM INVERSION STRATEGY

FWI is formulated as the minimization of a misfit function starting from an initial model $m_0$. Using the conventional $L^2$ misfit function (Tarantola 1984), FWI is formulated as

$$ \min_m = \frac{1}{2} \|(d_{\text{cal}}(m) - d_{\text{obs}})\|_{L^2}^2, \quad (1) $$

where $d_{\text{obs}}$ denotes the observed data, and $d_{\text{cal}}(m)$ is the solution of the forward modelling problem

$$ A(m)u = s, \quad (2) $$

$$ d_{\text{cal}} = Ru, \quad (3) $$

where $A(m)$ is a wave propagation modelling operator, based on the model parameter $m$, $u(m)$ is the simulated wavefield, $s$ is a seismic source and $R$ is an extraction operator extracting the wavefield values at the receiver positions based on the prior knowledge of the acquisition geometry. FWI applications in time-lapse problems seek to recover induced changes in the subsurface model using multiple data sets from different acquisition periods. For two surveys separated in time, those data sets (and the associated models) are called baseline and monitor, and represent the data prior and after induced changes in the subsurface, respectively. In this case, FWI can be applied to the data sets following various strategies. We provide a short review of the existing strategies below, as illustrated in Fig. 1.

### 2.1 Parallel

The parallel strategy considers two independent inversions, where the baseline and monitor data sets are inverted separately using the same starting model (Plessix & Perkins 2010) as illustrated in

Fig. 1(a). The time-lapse change is estimated by a simple subtraction between the reconstructed monitor and baseline models.

### 2.2 Sequential

Time-lapse FWI can be performed in a sequential manner (Routh et al. 2012; Asnaashari et al. 2015). It can be performed in forward fashion, from baseline to monitor, or in reverse fashion, from monitor to baseline. As illustrated in Fig. 1(b), the first inversion is performed on one selected data set to reconstruct a first subsurface model. In a following step, it is used as an initial model to invert the second data set. The underlying assumption is that once the background model is well estimated through the first inversion, the subsequent inversion will only update the missing perturbations related to the time lapse changes in the data. Starting from the baseline model should reduce the number of necessary iterations for the inversion of the monitor data set. At the end, the time-lapse model is obtained by a subtraction between the two reconstructed models.

### 2.3 Double difference

In principle, time-lapse changes between different seismic survey can be observed directly on the seismic data differences. In this spirit, inverting directly the data differences between baseline and monitor data instead of the full data set should help to focus the inversion on the related time-lapse model changes. This is the essence of the double-difference strategy (Watanabe et al. 2004; Denli & Huang 2009), which can be formulated as follows. We introduce the double-difference data as

$$ \Delta d = (d_{\text{obsm}} - d_{\text{obsb}}) - (d_{\text{calm}} - d_{\text{calb}}) \quad (4) $$

where $d_{\text{obsm}}$ and $d_{\text{obsb}}$ denote the monitor and baseline observed data, respectively, while $d_{\text{calm}}$ and $d_{\text{calb}}$ denote the monitor and observed
calculated data, respectively. A prior inversion of the baseline data set is necessary to be able to compute $d_{calb}$ based on the reconstructed baseline model similar to what is performed following the sequential strategy. The overall workflow of double difference strategy is illustrated in Fig. 1(c). We can then reorder eq. (4) to create a composite data set as

$$ \Delta d = (d_{\text{obsm}} - d_{\text{obsb}} + d_{\text{calb}}) - d_{\text{calm}} = d_{\text{comp}} - d_{\text{calm}}. \quad (5) $$

FWI problem is then reformulated as

$$ \min_m = \frac{1}{2} \|d_{\text{comp}} - d_{\text{cal}}(m)\|_{L^2}^2. \quad (6) $$

The aim of double-difference strategy is to emphasize the updates directly related to the perturbation on the data with an assumption that the background model is sufficiently accurate to predict the P-value of the data. As in sequential strategy, the time-lapse model can be estimated from the subtraction between the reconstructed monitor and baseline models.

## 2.4 Central difference

In the central difference strategy (Zhou & Lumley 2021a), two sequential inversions are performed independently starting from the same initial model in forward and reverse manner, with a total of four inversions runs as shown in Fig. 1(d). Each forward and reverse sequential inversion, yield a 4-D model difference we denote by $\Delta m^+$ and $\Delta m^-$, respectively. Zhou & Lumley (2021a) show that these two 4-D model differences have different signs at the locations where most of the artifacts are located, while the sought anomalies have a consistent sign. Based on this observation, and inspired by the finite-difference method, Zhou & Lumley (2021a) propose to estimate the 4-D changes as the average of the two 4-D model differences as

$$ \Delta m = \frac{\Delta m^+ + \Delta m^-}{2}. \quad (7) $$

## 2.5 Joint time lapse

Maharramov & Biondi (2015) propose a joint time-lapse FWI (JT-FWI) to invert baseline and monitor data simultaneously (Fig. 1e) solving the following problem

$$ \min_{m_b, m_m} = \frac{1}{2} \|d_{\text{cal}}(m_b) - d_{\text{obsb}}\|_2^2 + \frac{1}{2} \|d_{\text{cal}}(m_m) - d_{\text{obsm}}\|_2^2 \\ + \frac{\delta}{2} \|\nabla(m_m - m_b - \Delta m^{\text{Prior}})\|_1. \quad (8) $$

Maharramov & Biondi (2015) parametrize the misfit function into two independent baseline $m_b$ and monitor $m_m$ models. The first two terms correspond to the data term for the baseline and monitor respectively, while the third term corresponds to the model difference regularization which couple the baseline and monitor data set. In Maharramov & Biondi (2015), a total variation regularization penalization is chosen, as it promotes sparsity in the 4-D changes model.

## 2.6 Simultaneous

Inspired by Maharramov & Biondi (2015) and Ayeni & Biondi (2010), we propose in this work a simultaneous 4-D strategy with a parametrization using a baseline model $m_b$ and a perturbation

model $\Delta m$:

$$ \min_{m_b, \Delta m} = \frac{1}{2} \|d_{\text{cal}}(m_b) - d_{\text{obsb}}\|_2^2 + \frac{1}{2} \|d_{\text{cal}}(m_b + \Delta m) - d_{\text{obsm}}\|_2^2. \quad (9) $$

In this formulation, we propose an objective function applied to the FWI problem which can couple the model parameters implicitly within the data term, instead of relying almost exclusively on regularization to couple the baseline and monitor models. In this case, we favour direct contribution from the data sets represented by their respective gradients as follows:

$$ \begin{aligned} \nabla f(m_b) &= J(m_b)^T (d_{\text{cal}}(m_b) - d_{\text{obsb}}) \\ &\quad + J(m_b + \Delta m)^T (d_{\text{cal}}(m_b + \Delta m) - d_{\text{obsm}}), \\ \nabla f(\Delta m) &= J(m_b + \Delta m)^T (d_{\text{cal}}(m_b + \Delta m) - d_{\text{obsm}}), \end{aligned} \quad (10) $$

where $J(m) = \frac{\partial d_{\text{cal}}}{\partial m}$ is the Jacobian matrix. The sum of the baseline and monitor gradients is directly contributing to the common updates for the baseline model, while the monitor model relies only on the perturbation updates, which is present only on the monitor data, translated to its gradient contribution. In addition, regularization terms can be added easily if necessary to help the inversion.

## 2.7 Discussion on time-lapse strategies

We discuss now the advantages and disadvantages of the previously mentioned strategies. Among all, the parallel strategy remains the simplest and most straightforward. It is not restricted by non-repeatability of the acquisition between baseline and monitor data, as the time-lapse model is directly estimated from subtraction between models. On the contrary, subtraction in the data domain need well-matched acquisition geometries between baseline and monitor, which is a severe drawback for the double difference strategy. However, due to nonlinearity of FWI, two independent inversions converge to different points, which may create 4-D artifacts in the time-lapse model. In addition, different levels and types of noise between data sets could create 4-D artifacts related to the different noise contribution.

Similar to the parallel strategy, the sequential strategy is able to work efficiently with different acquisition geometries for the baseline and monitor data sets. The potential drawback is whenever the reconstructed baseline model contain inversion artifacts. In this case, the sequential inversion can propagate further these artifacts, and potentially produce more artifacts in the final model.

As already mentioned, the double-difference strategy requires to have well-matched surveys (e.g. permanent sensors or OBC configurations) to perform the subtraction between data sets. In case of non-matching surveys, data difference assumption in the data domain would be violated, therefore making it prone to specific artifacts.

Regarding the central difference method, it is observed that if the baseline and monitor data set have different levels and types of noise, the efficiency of the approach is reduced. This method also requires twice as much computational effort than the other 4-D strategies.

The main advantage of the joint time-lapse FWI strategy is the capability to invert the data sets simultaneously and to provide a better coupling between data sets. However, the coupling in Maharramov & Biondi (2015) formulation relies heavily on the regularization. From eq. (8) we can see that without the regularization the data are only coupled by the competition between the two parts of the misfit within the optimization process. For additive regularization
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Synthetic 2-D Sleipner model: (a) true baseline model, (b) true monitor model, (c) true 4-D differences, (d) initial baseline model, (e) initial monitor model and (f) initial 4-D differences.</figcaption></figure>

in eq. (8), the choice of the hyperparameter $\delta$ is not trivial and directly affects the inversion convergence. Therefore, implementation of regularization that controls the coupling heavily could be suboptimal.

Our reparametrization in eq. (9) aims at providing a stronger coupling between the baseline and monitor inversion. The coupling is implicitly derived from the data itself, while maintaining the advantage of inverting the data sets simultaneously.

## 3 SYNTHETIC STUDY ON TIME-LAPSE STRATEGIES

In this section, we assess the robustness of our proposed simultaneous strategy in eq. (9) against other previously proposed strategies. We focus on a comparison between our strategy, the parallel strategy, and the joint-time lapse FWI strategy, as the latter has been shown to overcome most of the limitations associated with previously mentioned 4-D FWI methods. We present a synthetic case study based on a 2-D Sleipner acoustic model [Figs 2a–c ($V_p$ and $\rho$)], comprising a layered structure with a low velocity anomaly representing the CO₂ changes (Romdhane & Querendez 2014). We use a streamer-like acquisition for a total of 96 shots, with a shot spacing of 30 m. We first set the maximum offset to 3000 m, to accommodate relatively long offset, sufficient to sample the CO₂ anomaly target between 800–1000 m. We opt for such a long offset acquisition to focus this test on the efficiency of the 4-D FWI strategies before considering additional difficulties associated with short offset and lack of illumination at the reservoir depth.

The numerical tests are performed with TOYxDAC_TIME, a finite-difference full waveform modelling and inversion code developed in the SEISCOP project (Yang et al. 2016b). We model the data with a band-pass filtered source wavelet between 6–10 Hz. The modelling is performed with time-domain fourth-order finite-difference, within a computational area of 4.8 km x 1.2 km, for a total grid size of 481 x 121. The spatial discretization step $h$ is set to 10 m, and the time discretization step to 1 ms. Perfectly

matched layers (PML; Bérenger 1994) are applied to attenuate spurious reflections at the boundary, and a free surface boundary condition is applied at the surface to simulate and accommodate surface related multiples and ghosts. The initial model for the inversion is a Gaussian-smoothed version of the true model shown in Figs 2(d)–(f).

We first simulate noise-free data for the inversion. As illustrated in Fig. 3, the 4-D anomalies are represented by reflections at 1.5 s. In this synthetic study, we compare parallel, joint time-lapse FWI, and our simultaneous strategy under ideal conditions. For the joint time-lapse FWI strategy, the regularization is not enforced during the inversion, in order to focus the comparison on the coupling effect due to different parametrization. The comparison of the inversion results strategies is presented in Figs 4–6. The CO₂ anomaly from the parallel strategy in Fig. 4(c) appears less smooth compared with the two other strategies, some distinct artifacts around the CO₂ anomaly are also stronger than the others. The results from joint time-lapse FWI (Fig. 5) also shows more artifacts compared with our simultaneous strategy (Fig. 5). To better quantify the performance, we calculate the root-mean-square (RMS) error between the reconstructed 4-D model and the true 4-D model. In this case, the result from simultaneous strategy shows the lowest rms model error compared with the other strategies.

In the second part of this test, we add Gaussian noise to the data to test the robustness of different strategies. As illustrated in Fig. 7, the model reconstructed with the parallel strategy exhibits stronger artifacts outside the CO₂ anomaly, similar with the results from JTFWI in Fig. 8 while the simultaneous strategy in Fig. 9 exhibits cleaner results than the two others, and has the lowest rms model among other strategies. This can be explained because of the summation of the gradient from the two data set in the simultaneous strategy is able to improve the signal to noise ratio.

To compare the accuracy of the estimated CO₂ anomaly, we compare 1-D profiles extracted at $x = 1.5$ km in Fig. 10. In the noise-free case, the estimation of the CO₂ anomaly at 800–1000 m from different strategies is comparable. However, we can see that the
<figure><img src="image_5.png" /><figcaption>Figure 3. Synthetic shot gather generated based on 2-D Sleipner model: (a) baseline, (b) monitor and (c) differences.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4. Inversion results on noise-free data for parallel strategies: (a) baseline, (b) monitor and (c) 4-D differences.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5. Inversion results on noise-free data for JTFWI strategies: (a) baseline, (b) monitor and (c) 4-D differences.</figcaption></figure>

parallel strategy wrongly estimates the area below the CO₂, while both joint time-lapse FWI and our simultaneous strategy estimate it very well. In the noisy case, the simultaneous strategy has smoother features at the CO₂ location and estimates the area below CO₂ better than the two others. Based on these synthetic tests, our simultaneous strategy appears, in overall, to provide a more robust estimation of the time-lapse changes than the two other methods. To understand the trade-offs between the methods, we present Table 1, which summarizes various computational cost parameters. The main trade-off between JTFWI and the simultaneous strategy is that they double the computational cost compared with the parallel strategy.

# 4 COMBINATION OF REFLECTION ORIENTED WORKFLOW AND SIMULTANEOUS STRATEGY: METHODOLOGY AND VALIDATION

## 4.1 Methodology

We can generalize the $L^2$ formulation of the simultaneous strategy in eq. (9) to take into account a general misfit function $G(d_{obs}, d_{cal})$ other than the conventional least-squares misfit function.

$$f(m_b, \Delta m) = G(d_{obs_b}, d_{cal}(m_b)) + G(d_{obs_m}, d_{cal}(m_b + \Delta m)) \quad (11)$$
<figure><img src="image_8.png" /><figcaption>Figure 6. Inversion results on noise-free data for simultaneous strategies: (a) baseline, (b) monitor and (c) 4-D differences.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7. Inversion results on noisy data for parallel strategies: (a) baseline, (b) monitor and (c) 4-D differences.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8. Inversion results on noisy data for JTFWI strategies: (a) baseline, (b) monitor and (c) 4-D differences.</figcaption></figure>

For this general misfit function, we are going to take into account a JFWI misfit function, which consists of a weighted sum of two terms measuring the misfit between early arrivals and reflections of observed $d_{obs}$ and calculated data $d_{cal}$, respectively (Provenzano et al. 2023). This general misfit function $G(d_{obs}, d_{cal})$ can be expressed as

$$G(d_{obs}, d_{cal}) = G(W^e d_{obs}^e, W^e d_{cal}^e) + G(W^r d_{obs}^r, W^r d_{cal}^r) \quad (12)$$

where $W^e$ and $W^r$ are weighting matrices applied to early arrivals and reflections, respectively, isolating and scaling the diving waves

in $d^e$ and the reflection waves in $d^r$. This allows JFWI to simultaneously invert the diving waves and the reflected waves, thus enriching the low wavenumber component from the transmission at the shallow part and along the reflection wavepaths at the deep part (Zhou & Lumley 2021b). The spirit of combining JFWI with simultaneous time-lapse strategy is to reconstruct successfully 4-D changes from limited offset data set with a good parameter coupling, while minimizing the 4-D artifacts. In practice, we combine the simultaneous strategy and JFWI during velocity model building part of JFWI, while the impedance model estimation remains independent.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Inversion results on noisy data for simultaneous strategies: (a) baseline, (b) monitor and (c) 4-D differences.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Comparison of extracted 1-D profiles at *x* = 1.5 km: (a) noise-free case and (b) noisy case.</figcaption></figure>

<table><caption>**Table 1.** Computational cost for various time-lapse strategies.</caption><thead><tr><td>No.</td><td>Strategy</td><td>N<sub>run</sub></td><td>Memory</td><td>N<sub>forward modeling</sub></td><td>N<sub>gradient</sub></td></tr></thead><tbody><tr><td>1</td><td>Parallel</td><td>2</td><td>1x</td><td>200</td><td>1</td></tr><tr><td>2</td><td>JTFWI</td><td>1</td><td>2x</td><td>400</td><td>2</td></tr><tr><td>3</td><td>Simultaneous</td><td>1</td><td>2x</td><td>400</td><td>2</td></tr></tbody></table>

Note: The memory requirements are relative to the memory needed for a single inversion with the same inversion setup. The number of gradient is the number of gradient computation on each iteration.

Therefore, baseline and monitor data set are inverted simultaneously following JFWI workflow to produce smooth baseline model $m_b$ and smooth 4-D perturbation model $\Delta m$ at the same time. A full description of the JFWI methodology we employ can be found in our companion paper (Fachtony et al. 2025) where we apply it on the Sleipner baseline data set we consider in this study.

## 4.2 Synthetic validation

In this test, we want to assess the robustness of the combination of JFWI and our simultaneous time-lapse strategy on limited offset data, based on the 2-D synthetic Sleipner model. We reduce the maximum offset of the acquisition from 3000 to 1700 m, which corresponds to the maximum offset for the field data we had access
<figure><img src="image_13.png" /><figcaption>**Figure 11.** FWI results: (a) baseline, (b) monitor and (c) 4-D differences.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** JFWI results: (a) baseline, (b) monitor and (c) 4-D differences.</figcaption></figure>

to. In addition, we start the inversion with a 1-D linearly increasing model, to illustrate the capability of our combined JFWI + simultaneous time-lapse strategy to provide reliable 4-D model changes in a realistic configuration.

As a first step, and to better illustrate the importance of using JFWI with such short offset data, we apply a simultaneous time-lapse strategy with a conventional FWI. Both baseline and monitor data are inverted starting from the same 1-D initial model. In Figs 11(a)–(c), we show the reconstructed baseline and monitor, and the estimated 4-D model. The latter is poorly resolved. Only the outline of the interface of the 4-D change is reconstructed. The quantitative $V_p$ values are not properly estimated. This shows that the conventional FWI is not able to constrain the kinematic updates on the CO₂ anomaly due to the absence of low wavenumber updates from the diving waves.

To tackle the offset limitation of the synthetic data, both baseline and monitor data are inverted with our combined workflow. Thanks to JFWI, we are able to reconstruct a smooth high velocity layer at 800 m in both the baseline and monitor model as shown in Fig. 12. In addition, JFWI is able to predict small portions of the low velocity anomaly at 800–1000 m which are related to the CO₂ anomaly. In contrast to the previous FWI result, the effect of the low wavenumber contribution from reflections in JFWI is able to provide kinematic updates to constrain the CO₂ anomaly better, even if the 4-D changes are not perfectly reconstructed.

We perform a subsequent run of conventional FWI starting from the reconstructed JFWI model to test if we can improve the resolution of the JFWI model. In Figs 13(a)–(c), we show that this

subsequent run of FWI is able to focus further on the high velocity layer both in the baseline and monitor models, and to provide a better estimation of the CO₂ anomaly with more accurate $V_p$ value and its interface outline. The comparison of the extracted 1-D profiles from different inversions at $x = 1.5$ km (Fig. 14) vouch for the benefit of the JFWI which provides a good estimation of the CO₂ anomaly. The subsequent run of conventional FWI starting from the JFWI model further improves the CO₂ anomaly.

Based on these synthetic tests, we show that the simultaneous strategy is able to provide more robust time-lapse model compared with other 4-D strategies. We show the capability of our combined JFWI + simultaneous inversion workflow to tackle the offset limitation on the synthetic data set and to provide a better time-lapse model compared with the conventional FWI workflow. The next step is to apply it on the 3-D Sleipner data.

# 5 APPLICATION TO THE SLEIPNER DATA SET

## 5.1 Overview of the Sleipner data set

The Sleipner field located in the North Sea is one of the pioneering site of carbon and capture storage. The CO₂ is transported from the Sleipner West gas field, and injected in the Utsira formation at Sleipner East. This is a 200–300 m thick sandstone formation lying in between 800 and 1100 m depth. It is characterized by a weak cementation layer with high porosity and permeability, and a thick and impermeable shale layer sealing its top.
<figure><img src="image_15.png" /><figcaption>**Figure 13.** JFWI + FWI results: (a) baseline, (b) monitor and (c) 4-D differences.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Comparison of the extracted 1-D profiles at 1.5 km for different inversion on short offset data.</figcaption></figure>

The two data set considered in this study are the baseline data, acquired in 1994, and the monitor data, acquired in 2008. Originally, both of the data set were acquired with 3000 m cable length streamer data and a recording duration of 5.5 s for each shot. The air gun source was towed at 6 m depth with 18.75 m shooting interval. As illustrated in Fig. 15, the baseline and monitor data have different acquisition geometry layouts, implying a specific non-repeatability challenge for time-lapse analysis. Several processing steps have been applied to the data, including swell noise filtering, low-cut filtering at 6.0 Hz, time step resampling to 2.0 ms, and $t^2$ amplitude scaling for geometrical spreading compensation. More importantly, an offset restriction from 3000 to 1700 m, and a recording length reduction from 5.5 to 2.3 s have been imposed. This time/offset restriction, combined with the absence of low frequencies, pose a significant challenge to conventional FWI applications.

We show in Fig. 16(a) comparison between shot gathers from the baseline and monitor, located near the expected CO₂ change area. The CO₂ signature in the monitor data is quite strong. On

full band data, the CO₂ signature can be easily identified as the strong amplitude reflections around 1 s in the monitor data, which are absent in the baseline data. On low-frequency data, these strong amplitude reflections are still visible but shifted to around 1.5 s of time due to the effect of minimum-phase filtering. In general, the arrivals between the baseline and monitor are quite consistent except for arrivals interfering with the CO₂ anomaly. This implies that no major changes in terms of background model are expected when we invert the data. In this case, our efforts to reconstruct a sufficiently accurate baseline model in the companion paper (Fachtony et al. 2025) should help us to explain the arrivals in the monitor data quite well, except for the CO₂ signature. Therefore, we can expect that an inversion of the monitor data set to provide the CO₂ changes to the model within relatively few iterations.

## 5.2 Inversion setup

The inversion is performed within a computational area of 5 km x 8 km x 2.3 km. The inversion is performed under the acoustic approximation with the same time-domain finite-difference modelling and inversion engine as the one used for the synthetic study (TOYx-DAC_TIME, Yang et al. 2016b). The finite-difference discretization is performed using a fourth-order in space and second-order in time staggered grid for the wave equation modelling. The models are discretized in space with a grid spacing of $h = 20$ m for the starting 6–10 Hz band, to match the sea floor depth at 80 m and ensure at least five grid points per wavelength for the lowest velocity. A time step of $\Delta t = 2$ ms is chosen to satisfy the Courant–Friedrichs–Lewy stability condition for this explicit scheme. Hicks (2002) interpolation is used to interpolate source and receivers positions in between grid points. A free surface condition is applied on the top of the model to account for ghosts and multiples in the data. PML are applied to attenuate spurious reflections to the other boundaries (Bérenger 1994).

A shot decimation strategy is applied to both baseline and monitor data to reduce the computational cost, selecting one shot out of three. After shot decimation, the final data set consists of a total of 2512 shots for the baseline and 2622 shots for monitor. An additional random subsampling strategy with a batch size of 256 sources is applied to further decrease the computational cost (Warner et al. 2013; Kamath et al. 2021). Prior to inversion, a minimal processing is applied to the data, which includes reverse $t^2$ amplitude scaling and minimum phase band-pass frequency filtering. The inversion is carried out in a multiscale manner (Bunks et al. 1995), from 6–10 Hz to 6–14 Hz with 2 Hz frequency increment.
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Overview of the acquisition geometry. The blue dots represent the source positions of the baseline, the green dots represent the source positions of the monitor, and the grey area is the receivers positions. The red circle is the injection well location.</figcaption></figure>

For JFWI, two independent windowing functions are applied to the baseline and monitor to separate transmission and reflection arrivals. The design of an appropriate window for this separation is assisted with synthetic data computation using the initial model to identify a proper time-offset separation range between transmission and reflection. A taper zone of 0.15 s is applied to avoid contamination between transmission and reflection arrivals and to accommodate traveltime changes during inversion. The windowing function is applied to each shot and is re-estimated for each frequency band. The taper zone is also reduced as the frequency increases. For each data, a shot-averaged source wavelet, shown in Fig. 17, is estimated using a frequency domain method, selecting near-offset part of the direct arrivals.

## 5.3 Baseline inversion

In this study, we follow the workflow outlined in Fig. 18, which consists of two phases of model reconstruction: baseline reconstruction and 4-D reconstruction. The baseline reconstruction phase is described in more detail in the companion paper (Fachtony et al. 2025), but we will also provide a summary of it in this section. The baseline model from the first phase, obtained through JFWI, serves as the starting point for our 4-D strategy. To summarize the processes, we apply a reflection oriented FWI to the baseline data starting from a 1-D linearly increasing $V_p$ model and a density model derived from the following Gardner's relationship

$$ \rho = \begin{cases} 1000 \text{ kg/m}^3 & \text{if depth} \le 80 \text{ m} \\ 320 V_p^{0.25} \text{ kg/m}^3 & \text{if depth} > 80 \text{ m.} \end{cases} \qquad (13) $$

Initially, we reconstruct an $I_p$ model through impedance wave-form inversion (IPWI), using short offset reflections up to 600 m of offset. An asymptotic preconditioner makes it possible a fast and

robust $I_p$ reconstruction (Provenzano et al. 2023). This impedance inversion is performed in three iterations for each batch of 256 sources, for a total of 36 iterations, covering all available shot gathers within one cycle. In Fig. 19, we present the reconstructed $I_p$ model, which shows dominant layered interfaces, as expected from the Sleipner area geology which exhibit fairly weak lateral variations (Queißer & Singh 2013). Additionally, in Fig. 19 we display a depth slice of the $I_p$ model at 900 m depth, representing structures within the Utsira formation before CO₂ injection.

A subsequent run of JFWI is performed using the initial 1D $V_p$ model and the reconstructed $I_p$ model as prior reflectors. JFWI is performed on full offset up to 1700 m, for 72 preconditioned steepest descent (ST) iterations, which is equivalent to two cycles of three iterations of 256-shot random subsampling batch. To balance the contribution between reflections and diving waves, the weight of the diving waves term in JFWI is set to half of the reflections term. To compensate the contribution of the deeper part, a linear depth dependent preconditioner is applied to the gradient. Utilizing the pseudo-time domain formulation of JFWI, the positions of the $I_p$ model interface are updated simultaneously in depth along with the $V_p$, thus avoiding the need to re-estimate the $I_p$ during the JFWI iterations, and maintain consistency between kinematics and reflectivity at short offset.

The final JFWI model is shown in Fig. 20 revealing meaningful tomography updates. In the companion paper (Fachtony et al. 2025), we illustrate that JFWI provides a tomographic model which is consistent with the well log data and fit the observed data quite well. This is crucial, as a well-reconstructed background model is a pre-requisite to focus on 4-D updates rather than simultaneously updating both the background model and the 4-D changes. In this approach, we effectively decouple the 3-D problem from the 4-D problem, addressing them in a sequential manner.
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Common shot gather for different frequency band on baseline and monitor data: (a) the full frequency band for baseline data, (b) 6–14 Hz filtered baseline data, (c) the full frequency band for monitor data and (d) 6–14 Hz filtered monitor data.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 17.** Estimated source wavelet at 6–10 Hz frequency band: (a) baseline and (b) monitor.</figcaption></figure>

## 5.4 Simultaneous time-lapse FWI

Starting from the final baseline JFWI $V_p$ model, we estimate the monitor $I_p$ model through IPWI. In this case, we expect a better reconstruction of the impedance discontinuities thanks to a better estimated background velocity model, at exception of the area where the $CO_2$ changes are expected. As the velocity of the $CO_2$ is not yet estimated, a possible depth shift of the $CO_2$ impedance contrast is expected in this first impedance model, which should be progressively updated during the $V_p$ inversion as the velocity is reconstructed, aided by the pseudo-time domain formulation. The reconstructed $I_p$ model is shown in Fig. 21. We can observe that the monitor $I_p$ model exhibits a different background compared with the baseline $I_p$. This difference is due to the improved estimation

of the background velocity model, which is used as a starting point. Strong impedance contrasts can be clearly seen in *x* and *y* directions at 800–1000 m depth. The depth slice at 900 m depth shows the distribution of strong impedance contrast related to the $CO_2$ changes.

Following the last part of the workflow in Fig. 18, we perform a simultaneous time-lapse JFWI, jointly inverting the baseline data and the monitor data. To do so, we use the baseline $V_p$ model as a common initial model for both baseline and monitor, while their respective $I_p$ models are used as a prior reflector during the simultaneous inversion. In this sense, we directly estimate the tomographic 4-D changes $\Delta m$. For the simultaneous inversion, we rely solely on reflection arrivals while the contribution from diving waves is
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Reflection oriented workflow for time lapse FWI. The workflow consists of two phases, baseline and 4-D reconstruction. On each phase, it consists of impedance reconstruction (IPWI) and velocity reconstruction (JFWI). During velocity reconstruction, the impedance model is also updated in terms of interface re-positioning following the velocity model update, through pseudo-time formulation. Simultaneous strategy (ST-JFWI) is performed at the final step to estimate baseline Vp and perturbation model within single inversion.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 19.** Reconstructed $I_p$ for the baseline data after IPWI. Upper-left panel is a depth slices cutting through the volume at $z = 0.9$ km, bottom-left panel is a vertical slice at $x = 2$ km, and bottom right is a vertical slice at $y = 3$ km. Slices locations are indicated by the black dashed-lines.</figcaption></figure>

muted. In the companion paper (Fachtony *et al.* 2025), we show that the diving waves mainly contribute to the overburden area. Given that the overburden is well estimated and no significant changes in the overburden are expected during CO₂ injection, we exclude the diving waves from the inversion. The reason is also that potential additional complexities for the reconstruction of this the near surface part of the model (e.g. changes in sea water velocity) are not considered.

In Fig. 22, we show the estimated 4-D difference model from the simultaneous inversion. The CO₂ changes are represented by a

400 m s⁻¹ velocity decrease. The delineation of the CO₂ anomaly in the $V_p$ model difference is also consistent with the $I_p$ model, the distribution of the CO₂ being more smeared, while the $I_p$ model shows a distribution staying within the expected 800–1100 m depth range. This difference in resolution is expected, as the scale separation within JFWI provides a high resolution the $I_p$ model and a tomographic resolution smooth $V_p$ model.

To introduce higher wavenumber contributions to the $V_p$ model, we perform a conventional mono-parameter FWI starting from these results. The parametrization is reset to velocity–density ($V_p-\rho$),
<figure><img src="image_22.png" /><figcaption>**Figure 20.** Reconstructed $V_p$ for the baseline data after JFWI workflow. Upper-left panel is a depth slices cutting through the volume at $z = 0.9$ km, bottom-left panel is a vertical slice at $x = 2$ km, and bottom right is a vertical slice at $y = 3$ km. Slices locations are indicated by the black dashed-lines.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 21.** Reconstructed $I_p$ for the monitor data after IPWI. Upper-left panel is a depth slices cutting through the volume at $z = 0.9$ km, bottom-left panel is a vertical slice at $x = 2$ km, and bottom right is a vertical slice at $y = 3$ km. Slices locations are indicated by the black dashed-lines.</figcaption></figure>

with the density parameter being kept fixed, derived from Gardner's relationship based on the final $V_p$ model from JFWI. We invert both the baseline and monitor data simultaneously, starting from the previously reconstructed JFWI $V_p$ models, to add higher wavenumber

contributions to these models. The result of the conventional FWI on the baseline data (Fig. 23) is satisfactory, as presented in the companion paper (Fachtony et al. 2025), the reconstructed model accurately follows the events observed in the well log data. After
<figure><img src="image_24.png" /><figcaption>**Figure 22.** Reconstructed 4-D differences after simultaneous inversion. Upper-left panel is a depth slices cutting through the volume at *z* = 0.9 km, bottom-left panel is a vertical slice at *x* = 1.9 km, and bottom right is a vertical slice at *y* = 3 km. Slices locations are indicated by the black dashed-lines.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 23.** Reconstructed Vp after simultaneous FWI starting from the reconstructed JFWI model. The upper-left panel is a depth slice cutting through the volume at *z* = 0.9 km. The bottom-left panel is a vertical slice at *x* = 1.9 km. The bottom right is a vertical slice at *y* = 3 km. Slices locations are indicated by the black dashed-lines.</figcaption></figure>

conventional FWI, additional resolution is added to the CO₂ distribution as shown in Fig. 24, although it only marginally helps with the depth smearing. In addition, the 1-D velocity profiles in Fig. 25 illustrate that the CO₂ changes are spanning between 750–1200 m, both in JFWI model and the conventional FWI. In this case, performing FWI on higher frequency could help to better localize the anomaly.

## 5.5 Quality control

We perform a data fit as a first quality control of the inversion result, as shown in Fig. 26 for a shot gather located near the CO₂ anomaly. The synthetic data in Fig. 26(a), computed using the reconstructed baseline model, shows that the diving waves of the baseline data are correctly predicted. In addition, we see that the synthetic data
<figure><img src="image_26.png" /><figcaption>**Figure 24.** Reconstructed 4-D differences after simultaneous FWI starting from the reconstructed JFWI model. The upper-left panel is a depth slice cutting through the volume at *z* = 0.9 km. The bottom-left panel is a vertical slice at *x* = 1.9 km. The bottom right is a vertical slice at *y* = 3 km. Slices locations are indicated by the black dashed-lines.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure 25.** 1-D profiles extracted from 4-D differences model at *x* = 2 km and *y* = 3 km location: (a) after JFWI and (b) after JFWI + FWI.</figcaption></figure>

computed using the monitor model in Fig. 26(b) fit the arrival times of the strong reflections around 1.5 s associated with the CO₂ injection in the monitor data. This indicates a satisfactory kinematic match of these reflections. The amplitude of these events is however not perfectly predicted. This suggests that elastic effects might need to be compensated for in the future with elastic FWI. We also illustrate the misfit function evolution along iteration in Fig. 27.

In Fig. 27(a), we can observe that JFWI significantly reduce the misfit during baseline data inversion. The fluctuations in Fig. 27(a) partially can be associated to the effect of source subsampling strategy, where different batch of sources are used during the inversion. Meanwhile, in Fig. 27(b), we can observe that the decrease of the misfit is less drastic and less fluctuative, thanks to prior reconstruction of baseline model.
<figure><img src="image_28.png" /><figcaption>**Figure 26.** Data fit QC on monitor shot gather near CO₂ location: (a) computed with the reconstructed baseline model and (b) computed with the reconstructed monitor model. We display the calculated data in two modes. The left panel shows a mirrored view, the left- and right-hand side represent observed data, while the middle part represent the synthetic data. The right panel shows an overlay view with a red and blue colour scale represent the synthetic data, while the observed data is shown in greyscale. If solely blue and black are apparent on the wiggles, a good fit is obtained, while red shows indicates unsatisfactory fit.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Figure 27.** (a) Misfit function evolution along iteration for baseline reconstruction. (b) Misfit function evolution along iteration for simultaneous baseline and 4-D model reconstruction.</figcaption></figure>

## 5.6 CO₂ saturation estimation

To relate the inverted P-wave velocity from inversion and CO₂ saturation, we utilize fluid substitution method (Gassmann 1951) which relates the dry rock bulk modulus $K_{dry}$ of a porous medium to the bulk modulus $K_f$ of the pore fluid and the bulk modulus $K_{grain}$ of the constitutive grain. First, we analyse the relationship between $V_p$ and saturation with the velocity-saturation curve shown in Fig. 28. To do so, we calculate an effective fluid bulk modulus using Brie's empirical law (Brie et al. 1995), given as

$$K_f = (K_w - K_{CO_2})S_w^e + K_{CO_2}, \quad (14)$$

where $K_f$ is the effective fluid bulk modulus, $K_w$ and $K_{CO_2}$ are the bulk modulus of brine and CO₂, respectively, and $S_w$ is the brine saturation. Finally, $e$ is the Brie's exponent, varying between 1 (patchy mixing) and 40 (uniform mixing). For rock physics calculation, we use the parameters described in Dupuy et al. (2017).

In our case, we also calculate a P-impedance-saturation curve shown in Fig. 28(b) to provide a better understanding of the relationship between $I_p$ and the CO₂ saturation. The curves confirm that both velocity and impedance change due to the CO₂ are variable depending on which Brie exponent is used. This shows that different exponent values representing different scenarios lead to large disparity of the $V_p$, suggesting ultimately high uncertainties of the CO₂ saturation in the Utsira formation.
<figure><img src="image_30.png" /><figcaption>**Figure 28.** a) P-velocity-saturation curve and b) P-impedance-saturation curve.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>**Figure 29.** CO₂ saturation derived from different parameters and Brie's exponent: (a) derived from *V*<sub>p</sub> JFWI and *e* = 1, (b) derived from *V*<sub>p</sub> FWI and *e* = 1, (c) derived from *I*<sub>p</sub> and *e* = 1, (d–f) the same as (a–c) but calculated with *e* = 3, and similarly (g–i) calculated with *e* = 5.</figcaption></figure>

Using the rock physics model, we calculate the CO₂ saturation at each grid points from the inverted velocity and impedance, using different value of Brie exponent. As the parameters described by Dupuy *et al.* (2017) are mainly related to Utsira sandstone formation, we limit the calculation between 800 and 1100 m of depth where the Utsira formation is expected. To estimate the CO₂ saturation, at each grid point, we compare the reconstructed *V*<sub>p</sub> or *I*<sub>p</sub> value from the inversion with the values obtained from rock physics calculation. If they do not match, the CO₂ saturation is increased or decreased and *V*<sub>p</sub> or *I*<sub>p</sub> is calculated again, until a match is achieved. First, we compare vertical slices of the saturation model at *y* = 3 km derived from *V*<sub>p</sub> model from JFWI, *I*<sub>p</sub> model, and *V*<sub>p</sub> model with different Brie's exponent as shown in Fig. 29. In this case, the CO₂ saturation derived from the *I*<sub>p</sub> model has a higher resolution compared with the others. This is consistent with the scale separation applied within the JFWI workflow. Among different Brie's exponent, the CO₂ saturation varies quite significantly. It ranges between 40–50 per cent for *e* = 1, around 20 per cent for *e* = 3, and less than 10 per cent for *e* = 5. In Fig. 30, we compare the CO₂ saturation derived from the *I*<sub>p</sub> model along different depths. We show that the highest saturation is located around 900–1000 m depth.

## 6 DISCUSSION

In our study, we use a reflection oriented workflow under the acoustic approximation, which remains a limitation for JFWI to accommodate the elastic effects in the field data. However, it has been shown in several previous studies, that the S-wave velocity is hard to be constrained by the Sleipner data (Raknes *et al.* 2015a). In this case, the use of multicomponent data might be necessary to successfully invert the S waves. Queißer & Singh (2013) shows that the Sleipner data are dominated by unconverted PP arrivals, although the mode conversions (e.g P–S–P waves) are still expected to be present in the data. However, these mode conversions occur in a longer propagation path and usually recorded at later time in seismic recording. With a recording length reduced to 2.3 s from the original 5.5 s, it seems difficult to record such converted waves. Nevertheless, some residual energy due to elastic effects probably remain in the data. On the other hand, the elastic reflection waveform inversion is still at a research stage, with some notable development by Masmoudi *et al.* (2021) and Aaker *et al.* (2024). Masmoudi *et al.* (2021) utilizes PS Born modelling to construct the PS sensitivity kernel while Aaker *et al.* (2024) suggests a PS mode-decomposition integrated into the misfit functional used for elastic reflection waveform inversion.
<figure><img src="image_32.png" /><figcaption>**Figure 30.** CO₂ saturation derived from *I*<sub>p</sub> model and different Brie's exponent at different depth: (a) at 800 m and *e* = 1, (b) at 900 m and *e* = 1, (c) at 1000 m and *e* = 1, (d–f) the same as (a–c) but calculated with *e* = 3, and similarly (g–i) calculated with *e* = 5.</figcaption></figure>

The impact of attenuation is also a factor that introduces uncer-
tainty in the inversion results, which is not accounted for in our
study. In general, attenuation mainly affect transmitted waves, mod-
ifying their amplitudes and distorting their phase, which in turns has
an effect on the reconstruction of the long wavelength variations of
the velocity. For near offset data, attenuation mainly affects the am-
plitudes of the reflections, potentially causing an underestimation of
the velocity. In our case, we see the introduction of the attenuation in
the modelling and inversion as a straightforward path towards better
velocity and impedance reconstruction. In addition, using multipa-
parameter FWI to estimate attenuation could be considered. This is
however a difficult problem, requiring a good prior knowledge on
the velocity model (see Yong *et al.* 2021 for instance). In this case,
our estimated velocity model could serve as an initial model for
subsequent attenuation model estimation.

It is interesting to compare the results in these studies with some
previous study on the Sleipner data. Queißer & Singh (2013) per-
formed 2-D elastic FWI on a single line from Sleipner data with
a maximum offset of 3500 m, reaching up to 38 Hz frequency. In
their time-lapse model (Queißer & Singh 2013, fig. 10), several
low velocity layers are visible at 800–1100 m, associated with CO2
changes. In these layers a velocity reduction up to 800 m s⁻¹ can be
observed. In our results, the low velocity perturbation on the time-
lapse model is at a similar depth range, despite the depth smearing
of the anomaly due to resolution limits of the JFWI Vp model.
Besides, our reconstructed Ip model complement our Vp model,
showing similar layering with strong impedance contrast associated
with the CO₂ response in a more confined zone.

In our results, the velocity reduction associated with the CO₂
anomaly is compared with the results of (Queißer & Singh 2013),
reaching only up to 400 m s⁻¹. From previous studies, the expected
minimum velocity of the area affected by the presence of CO₂
depend on the rock physical parameters used. Arts et al. (2002)

and Chadwick *et al.* (2005) expect a minimum velocity due to CO₂
injection to be around 1400 m s⁻¹, while Carcione *et al.* (2006)
expected it to be around 1200 m s⁻¹. In this case, the result of
(Queißer & Singh 2013) favours the scenario suggested by Carcione
*et al.* (2006), while the scenario suggested by Arts *et al.* (2002) and
Chadwick *et al.* (2005) could be more favourable to our result,
which is around 1600 m s⁻¹. The difference on the magnitude of
the velocity after the CO₂ injection could also be associated to the
fact that the frequency content used in our inversion is still lower
than the maximum frequency band considered in Queißer & Singh
(2013).

However, it is expected that a 3-D inversion is more beneficial
to derive better overall distribution of the CO₂ changes compared
with a 2-D inversion. In another study, Raknes *et al.* (2015b) es-
timated the CO₂ changes on limited offset data with 3-D elastic
FWI. In this case, the estimated time-lapse model is also associated
with a velocity reduction around 800–1100 m of depth, but widely
spread along the model. In addition, the time-lapse model shows
strong artifacts related to high velocity perturbations which are not
expected, while this study and Queißer & Singh (2013) do not show
such strong high-velocity perturbation. In this case, the artifacts
might be associated to the offset limitation on the Sleipner data,
such that the expected target at the deeper part are not well con-
strained. In our study, we show the benefit of the JFWI to constrain
better the deeper part of the model, thereby stabilizing the velocity
updates associated with the CO₂ changes.

Previous studies at Sleipner also show that the CO₂ saturation is
difficult to estimate reliably. The lack of information on how the
fluids are mixed in this porous medium result in a strong uncer-
tainty (Queißer & Singh 2013; Ghosh *et al.* 2015; Dupuy *et al.*
2017; Yan *et al.* 2019). Two scenario represent two ends of a con-
tinuum of possibilities, from a patchy mixing to a uniform mixing.
The use of rock physics inversion such as in Dupuy *et al.* (2021)
could help quantify the uncertainties and provide more accurate estimations where lateral variations of the Brie's exponent could be estimated. In addition, using additional inputs from seismic inversion (attenuation, shear wave velocity, Amplitude Versus Offset (AVO) attributes) or from other geophysical techniques (gravity, Controlled Source Electro-Magnetic (CSEM)) could help to reduce the non-uniqueness of the inversion (Gao et al. 2012; Liang et al. 2016; Dupuy et al. 2021). Compared with previous studies, the highest CO₂ saturation in our study is located around 900–1000 m of depth, while the result from Dupuy et al. (2017) and Queißer & Singh (2013) suggest that it should be concentrated between 870 and 950 m. This difference might be attributed to the smearing effect from JFWI and the difference in frequency content, both previous study making use of higher frequency data.

Finally, the CO₂ changes at Sleipner have been estimated in previous studies using conventional seismic time-migration (Arts et al. 2004; Eiken et al. 2005; Chadwick et al. 2010). It is interesting to compare these conventional results with the reconstructed models from this study. The horizontal extent of the gas cloud is expected to be approximately 3.6 km × 1.0 km (Chadwick et al. 2009). It is also expected to follow a migration path toward a north-south trend (Bickle et al. 2007; Chadwick et al. 2010, fig. 1). Our reconstructed *I*<sub>p</sub> model shows a similar trend of the strong impedance contrast, with an extent of around 4 km × 1.5 km, and north-south orientation. Again, the difference in magnitude between the *I*<sub>p</sub> model and the migrated image could be attributed to the difference in frequency content, the latter being constructed from the full frequency band. The Utsira formation is expected to be located between 800 and 1100 m, where the CO₂ changes are expected. In our result, we can observe velocity changes below these depth which can be associated to vertical smearing due to JFWI effect. In this case, JFWI only provides low resolution velocity model as it is constructed only by low-wavenumber component gradient. In addition, As we work with relatively low frequency, the shape of prior reflectors from the reconstructed impedance model is ‘thicker’ due to finite frequency effect. In this case, ‘rabbits ears’ from reflection will sample deeper, therefore produce depth smearing.

## 7 CONCLUSION

We have illustrated the application of reflection oriented FWI in the context of time-lapse problems to the limited offset 3-D Sleipner data. To address the challenges associated with time-lapse FWI problem, we propose a simultaneous strategy which inverts both the baseline and monitor data set simultaneously with a baseline model *m* and a perturbation model *Δm* parametrization. To test the robustness of the proposed time-lapse strategy, we perform synthetic tests on a 2-D Sleipner model. The results of our synthetic tests with long offset scenario show that the simultaneous strategy is able to provide better time-lapse models compared with other 4-D strategies, and that it is also more robust against the presence of noise. In addition, the synthetic tests with short offset scenario show the benefit of combining this approach with a reflection-oriented workflow, JFWI, to better constrain the CO₂ at depth. A subsequent run of conventional FWI starting from the JFWI model is then able to provide a higher resolution perturbation model. In the field data application, this combined strategy is able to provide reliable CO₂ changes at an expected depth. To show the reliability of the results, we present data fit as a quality check and compare with the previous results from other studies. The estimated CO₂

changes are characterized by a velocity reduction up to 400 m s⁻¹, and strong impedance contrast between 800 and 1000 m depth, consistent with previous JFWI studies. In addition, the distribution of CO₂ changes aligns with conventional seismic time-migration from prior studies, following a north-south trend migration path.

## ACKNOWLEDGMENTS

This study was partially funded by the SEISCPEO consortium (http://seiscpe2.osug.fr), sponsored by AKERBP, DUG, EXXON-MOBIL, GEOLINKS, JGI, PETROBRAS, SHELL, SINOPEC, TOTALENERGIES and VIRIDIEN. This study was granted access to the HPC resources provided by the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by Grenoble research communities, the HPC resources of Cray Marketing Partner Network (https://partners.cray.com), and those of CINES/IDRIS/TGCC under the allocation 046091 made by GENCI. We thank Equinor for sharing Sleipner raw data. B. Boddupalli, B. Dupuy and A. Romdhane contributions have been produced with support from the NCCS Research Centre, performed under the Norwegian research programme Centre for Environment-friendly Energy Research (FME) with the following partners: Aker BP, Aker Carbon Capture, Allton, Ansaldo Energia, Baker Hughes, CoorsTek Membrane Sciences, Elkem, Eramet, Equinor, Gassco, Hafslund Oslo Celsio, KROHNE, Larvik Shipping, Norcem Heidelberg Cement, Offshore Norge, Quad Geometrics, Stratum Reservoir, Total-Energies, Vår Energi, Wintershall DEA and the Research Council of Norway (257579/E20).

## DATA AVAILABILITY

Data and code associated with this research belong to the SEISCPEO consortium. Any request related to the materials availability needs to be addressed to Romain Brossier (romain.brossier@univ-grenoble-alpes.fr) or Ludovic Métivier (ludovic.metivier@univ-grenoble-alpes.fr).

## REFERENCES

* Aaker, O., Pedersen, Ø. & Raknes, E., 2024. Elastic reflection waveform inversion with a novel PS mode-decomposed misfit functional, *85th EAGE Annual Conference & Exhibition Vol. 2024*, (1), p. 1–5.

* Arts, R., Elsayed, R., van der Meer, L., Eiken, O., Ostmo, S., Chadwick, A., Kirby, G. & Zinszner, B., 2002. Estimation of the mass of injected co2 at Sleipner using time-lapse seismic data. *64th EAGE Conference & Exhibition*, European Association of Geoscientists & Engineers.

* Arts, R., Eiken, O., Chadwick, A., Zweigel, P., van der Meer, L. & Zinszner, B., 2004. Monitoring of CO2 injected at Sleipner using time-lapse seismic data, *Energy*, **29**(9), 1383–1392.

* Asnaashari, A., Brossier, R., Garambois, S., Audebert, F., Thore, P. & Virieux, J., 2015. Time-lapse seismic imaging using regularized full waveform inversion with prior model: which strategy?, *Geophys. Prospect.*, **63**(1), 78–98.

* Ayeni, G. & Biondi, B., 2010. Target-oriented joint least-squares migration/inversion of time-lapse seismic data sets, *Geophysics*, **75**(3), R61–R73.

* Bérénger, J.-P., 1994. A perfectly matched layer for absorption of electromagnetic waves, *J. Comput. Phys.*, **114**, 185–200.

* Bickle, M., Chadwick, A., Huppert, H. E., Hallworth, M. & Lyle, S., 2007. Modelling carbon dioxide accumulation at Sleipner: implications for underground carbon storage, *Earth planet. Sci. Lett.*, **255**(1), 164–176.
Brie, A., Pampuri, F., Marsala, A. & Meazza, O., 1995. Shear sonic interpretation in gas-bearing sands, *SPE Annual Technical Conf. 30595*, Society of Petroleum Engineers, pp. 701–710.

Bunks, C., Salek, F. M., Zaleski, S. & Chavent, G., 1995. Multiscale seismic waveform inversion, *Geophysics*, **60**(5), 1457–1473.

Carcione, J., Picotti, S., Gei, D. & Rossi, G., 2006. Physics and seismic modelling for monitoring CO₂ storage, *Pure appl. Geophys.*, **163**, 175–207.

Chadwick, A. et al., 2010. Quantitative analysis of time-lapse seismic monitoring data at the Sleipner CO₂ storage operation, *Leading Edge*, **29**(2), 170–177.

Chadwick, R., Noy, D., Arts, R. & Eiken, O., 2009. Latest time-lapse seismic data from Sleipner yield new insights into CO₂ plume development, *Energy Procedia*, **1**(1), 2103–2110.

Chadwick, R. A., Arts, R. & Eiken, O., 2005. 4-D seismic quantification of a growing CO₂ plume at Sleipner, North Sea, in *Petroleum Geology: North-West Europe and Global Perspectives—Proceedings of the 6th Petroleum Geology Conference*, Geological Society of London.

Davy, R.G. et al., 2017. Resolving the fine-scale velocity structure of continental hyperextension at the Deep Galicia Margin using full-waveform inversion, *Geophys. J. Int.*, **212**(1), 244–263.

Denli, H. & Huang, L., 2009. Double-difference elastic waveform tomography in the time domain, *SEG Tech. Program Expand. Abstr.*, **28**(1), 2302–2306.

Dupuy, B., Romdhane, A., Eliasson, P., Querendez, E., Yan, H., Torres, V.A. & Ghaderi, A., 2017. Quantitative seismic characterization of CO₂ at the Sleipner storage site, north sea, *Interpretation*, **5**(4), SS23–SS42.

Dupuy, B., Romdhane, A., Nordmann, P.-L., Eliasson, P. & Park, J., 2021. Bayesian rock-physics inversion: Application to CO₂ storage monitoring, *Geophysics*, **86**(4), M101–M122.

Eiken, O., 2019. Twenty years of monitoring CO₂ injection at Sleipner, in *Geophysics and Geosequestration*, pp. 209–234, eds, Davis, T.L., Landrø, M. & Wilson, M., Cambridge University Press.

Eiken, O., Brevik, I., Arts, R., Lindeberg, E. & Fagervik, K., 2005. Seismic monitoring of CO₂ injected into a marine aquifer, *SEG Technical Program Expanded Abstracts 2000*, Society of Exploration Geophysicists, pp. 1623–1626.

Fachtony, F., Brossier, R., Métivier, L., Boddupali, B., Dupuy, B. & Romdhane, A., CO₂ monitoring at Sleipner field using reflection oriented full waveform inversion: Part 1 - baseline reconstruction *Geophysical Journal International* 2025, **241**, 1226–1242.

Fichtner, A. & Villaseñor, A., 2015. Crust and upper mantle of the western mediterranean-constraints from full-waveform inversion, *Earth planet. Sci. Lett.*, **428**, 52–62.

Furre, A.-K. & Eiken, O., 2014. Dual sensor streamer technology used in Sleipner CO₂ injection monitoring, *Geophys. Prospect.*, **62**(5), 1075–1088.

Gao, G., Abubakar, A. & Habashy, T.M., 2012. Joint petrophysical inversion of electromagnetic and full-waveform seismic data, *Geophysics*, **77**(3), WA3–WA18.

Gassmann, F., 1951. Über die elastizität vorüber medien, *Vierteljahrsschrift der Naturforschenden Gesellschaft in Zürich*, **96**, 1–23.

Ghosh, R., Sen, M.K. & Vedanti, N., 2015. Quantitative interpretation of CO₂ plume from Sleipner (north sea), using post-stack inversion and rock physics modelling, *Int. J. Greenhouse Gas Control*, **32**, 147–158.

Górszczyk, A., Brossier, R. & Métivier, L., 2021. Graph-space optimal transport concept for time-domain full-waveform inversion of ocean-bottom seismometer data: Nankai trough velocity structure reconstructed from a 1-D model, *J. Geophys. Res. Solid Earth*, **126**(5), e2020JB021504, doi:10.1029/2020JB021504.

Guasch, L., Agudo, O.C., Tang, M.-X., Nachev, P. & Warner, M., 2020. Full-waveform inversion imaging of the human brain, *NPJ Digital Medicine*, **3**(1), 1–12.

Hicks, E., Hoeber, H., Houbiers, M., Lescoffit, S.P., Ratcliffe, A. & Vinje, V., 2016. Time-lapse full-waveform inversion as a reservoir-monitoring tool — A North Sea case study, *Leading Edge*, **35**(10), 850–858.

Hicks, G.J., 2002. Arbitrary source and receiver positioning in finite-difference schemes using Kaiser windowed sinc functions, *Geophysics*, **67**, 156–166.

Irnaka, T.M., Brossier, R., Métivier, L., Bohlen, T. & Pan, Y., 2022. 3-D Multi-component full waveform inversion for Shallow-seismic target: Ettlingen line case study, *Geophys. J. Int.*, **229**(2), 1017–1040.

Kamath, N., Brossier, R., Métivier, L., Pladys, A. & Yang, P., 2021. Multi-parameter full-waveform inversion of 3-D ocean-bottom cable data from the Valhall field, *Geophysics*, **86**(1), B15–B35.

Lailly, P., 1983. The seismic problem as a sequence of before-stack migrations, in *Conference on Inverse Scattering: Theory and Applications*, SIAM, Philadelphia.

Liang, L., Abubakar, A. & Habashy, T.M., 2016. Reservoir property mapping and monitoring from joint inversion of time-lapse seismic, electromagnetic, and production data, *Geophysics*, **81**(5), ID73–ID84.

Lumley, D.E., 2001. Time-lapse seismic reservoir monitoring, *Geophysics*, **66**(1), 50–53.

Maharramov, M. & Biondi, B., 2015. Robust simultaneous time-lapse full-waveform inversion with total-variation regularization of model difference, *77th EAGE Conference & Exhibition, Expanded Abstracts, We P3 09*, European Association of Geoscientists & Engineers.

Masmoudi, N., Ratcliffe, A., Wang, M., Xie, Y. & Wang, T., 2021. A practical implementation of converted-wave reflection full-waveform inversion, *82nd EAGE Annual Conference & Exhibition*, **2021**(1), 1–5.

Métivier, L., Brossier, R., Mérigot, Q. & Oudet, E., 2019. A graph space optimal transport distance as a generalization of $L^p$ distances: application to a seismic imaging inverse problem, *Inverse Problems*, **35**(8), 085001, doi:10.1088/1361-6420/ab206f.

Mispel, J., Furre, A., Sollid, A. & Maao, F.A., 2019. High frequency 3-D FWI at Sleipner: A closer look at the CO₂ plume, in *81st EAGE Annual Conference & Exhibition, Vol. 2019*, pp. 1–5, European Association of Geoscientists & Engineers.

Operto, S., Miniussi, A., Brossier, R., Combe, L., Métivier, L., Monteiller, V., Ribodetti, A. & Virieux, J., 2015. Efficient 3-D frequency-domain mono-parameter full-waveform inversion of ocean-bottom cable data: application to Valhall in the visco-acoustic vertical transverse isotropic approximation, *Geophys. J. Int.*, **202**(2), 1362–1391.

Plessix, R.E., 2013. A pseudo-time formulation for acoustic full waveform inversion, *Geophys. J. Int.*, **192**, 613–630.

Plessix, R.E. & Perkins, C., 2010. Full waveform inversion of a deep water ocean bottom seismometer data set, *First Break*, **28**, 71–78.

Provenzano, G., Vardy, M.E. & Henstock, T.J., 2018. Decimetric-resolution stochastic inversion of shallow marine seismic reflection data: dedicated strategy and application to a geohazard case study, *Geophys. J. Int.*, **214**, 1683–1700.

Provenzano, G., Brossier, R. & Métivier, L., 2023. Robust and efficient waveform-based velocity-model building by optimal transport in the pseudotime domain: methodology, *Geophysics*, **88**(2), U49–U70.

Queiüber, M. & Singh, S.C., 2013. Full waveform inversion in the time lapse mode applied to CO₂ storage at Sleipner, *Geophys. Prospect.*, **61**(3), 537–555.

Raknes, E.B., Arntsen, B. & Weibull, W., 2015a. Three-dimensional elastic full waveform inversion using seismic data from the sleipner area, *Geophys. J. Int.*, **202**(3), 1877–1894.

Raknes, E.B., Weibull, W. & Arntsen, B., 2015b. Seismic imaging of the carbon dioxide gas cloud at Sleipner using 3-D elastic time-lapse full waveform inversion, *Int. J. Greenhouse Gas Control*, **42**, 26–45.

Romdhane, A. & Querendez, E., 2014. CO₂ characterization at the Sleipner field with full waveform inversion: application to synthetic and real data, *Energy Procedia*, **63**, 4358–4365.

Routh, P., Palacharla, G., Chikichev, I. & Lazaratos, S., 2012. Full wavefield inversion of time-lapse data for improved imaging and reservoir characterization, *SEG Technical Program Expanded Abstracts 2012*, Society of Exploration Geophysicists, pp. 1–6.

Tarantola, A., 1984. Inversion of seismic reflection data in the acoustic approximation, *Geophysics*, **49**(8), 1259–1266.
Tarantola, A., 1986. A strategy for nonlinear inversion of seismic reflection data, *Geophysics*, **51**(10), 1893–1903.

Thrastarson, S., van Herwaarden, D., Noe, S., Josef Schiller, C. & Fichtner, A., 2024. REVEAL: a global full-waveform inversion model, *Bull. seism. Soc. Am.*, **114**(3), 1392–1406.

Warner, M., Nangoo, T., Shah, N., Umpleby, A. & Morgan, J., 2013. *SEG Technical Program Expanded Abstracts 2013*, Full-waveform inversion of cycle-skipped seismic data by frequency down-shifting, Society of Exploration Geophysics, pp. 903–907.

Watanabe, T., Shimizu, S., Asakawa, E. & Matsuoka, T., 2004. Differential waveform tomography for time-lapse crosswell seismic data with application to gas hydrate production monitoring, *SEG Tech. Program Expand. Abstr.*, **23**(1), 2323–2326.

Yan, H., Dupuy, B., Romdhane, A. & Arntsen, B., 2019. CO2 saturation estimates at Sleipner (north sea) from seismic tomography and rock physics inversion, *Geophys. Prospect.*, **67**(4), 1055–1071.

Yang, D., Liu, F., Morton, S., Malcolm, A. & Fehler, M., 2016a. Time-lapse full-waveform inversion with ocean-bottom-cable data: application on Valhall field, *Geophysics*, **81**(4), R225–R235.

Yang, P., Brossier, R., Métivier, L. & Virieux, J., 2016b. Wavefield reconstruction in attenuating media: A checkpoint-assisted reverse-forward simulation method, *Geophysics*, **81**(6), R349–R362.

Yong, P., Brossier, R., Métivier, L. & Virieux, J., 2021. Q estimation by full-waveform inversion: analysis and misfit functions comparison, *First International Meeting for Applied Geoscience & Energy*, pp. 827–831, Society of Exploration Geophysicists.

Zhou, W. & Lumley, D., 2021a. Central-difference time-lapse 4-D seismic full-waveform inversion, *Geophysics*, **86**(2), R161–R172.

Zhou, W. & Lumley, D., 2021b. Time-lapse joint full-waveform inversion of diving and reflected waves, *First International Meeting for Applied Geoscience & Energy Expanded Abstracts*, Society of Exploration Geophysicists, pp. 3495–3499.

Zhou, W., Brossier, R., Operto, S. & Virieux, J., 2015. Full waveform inversion of diving & reflected waves for velocity model building with impedance inversion based on scale separation, *Geophys. J. Int.*, **202**(3), 1535–1554.