# CO2 monitoring at Sleipner field using reflection oriented full waveform inversion: Part 1 - baseline reconstruction

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 1 - baseline reconstruction. Geophysical Journal International, 2025, 241 (2), pp.1226-1242. 10.1093/gji/ggaf097. hal-05336515

HAL Id: hal-05336515

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

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 1: baseline reconstruction

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 February 28. Received 2025 February 20; in original form 2024 October 21

## SUMMARY

In this study, we illustrate the application of a 3-D reflection oriented workflow for full waveform inversion (FWI) to the offset data from the Sleipner field. The data set is having maximum offset of less than 2000 m, and has been pre-processed with a low-cut filter below 6 Hz, imposing strong challenges for the FWI application. To tackle these challenges, a reflection oriented FWI is applied to the data set, which utilize joint full waveform inversion (JFWI) to constrain the low wavenumber updates at the deepest part of the reconstructed model. It consists of two steps, an impedance model building serving as a prior reflector information followed by a velocity model building. In this case, JFWI workflow is taking advantage of the pseudo-time formulation to honour the zero offset traveltime, fast and robust asymptotic pre-conditioner for impedance model building, and graph space optimal transport misfit function to mitigate cycle skipping. To show the effects of limited offset, conventional FWI is performed. In this case, it is clear that the meaningful updates coming from the diving waves are restricted to the shallow part no deeper than 500 m of depth, while no meaningful perturbation are observed beyond the diving waves penetration. Taking advantage of the meaningful shallow updates, diving wave only inversion is performed prior to the impedance model building, and then followed by the JFWI workflow. The results of the field data application show that JFWI is able to produce meaningful velocity updates both in shallow part and the deeper part. The result is supported by satisfactory fit of the calculated data based on the JFWI model compared to the observed data. In addition, the velocity model fits the low wavenumber trend of the well log data. A subsequent run of conventional FWI is performed starting from JFWI model, in order to improve the resolution of the velocity model. The results is able to introduce higher wavenumber content to the velocity model, producing satisfactory fit with the observed data, and matching the well log data.

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

## 1 INTRODUCTION

Nowadays, the importance of carbon capture and storage (CCS) is rising. CCS by definition is a process consisting of the separation of CO₂ from industrial and energy-related sources, its transport to a storage location and long-term isolation from the atmosphere. It is considered by the Intergovernmental Panel on Climate Change (IPCC) as an effective means to reduce the concentration of CO₂ in the atmosphere and mitigate the effects of the global climate

changes. Based on the existing project around the world, dedicated CO₂ storage capacity could reach around 615 Mt CO₂ yr⁻<fnref n="1" /> by 2030. However, this remains insufficient to meet the target of 1000 Mt CO₂ yr⁻<fnref n="1" /> by 2030 for the Net Zero Emissions (NZE) scenario (CCS 2023). The development of CCS facilities in major regions has advanced significantly in 2024, with United States and Europe being at the forefront. The Asia-Pacific region is also rising their commitment towards CO₂ infrastructures, several projects are announced by China, Malaysia, Singapore and Thailand. Further expansion of the CO₂ management industry is also starting in Japan, Korea and Indonesia. The Sleipner storage facility located in the North Sea

<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</footnote>
of Norway is one of the pioneering CO₂ storage offshore, the first operation at industrial scale. Since 1994, the Sleipner facility has been piloted and the first CO₂ injection into the Utsira formation started at 1996. Approximately, 1 Mt of CO₂ are being injected every year. At the Sleipner field, CO₂ has been transported from the produced natural gas field and injected into the Utsira subsurface (Eiken 2019).

In order to monitor the spreading of CO₂ and to ensure storage safety at Sleipner facility, several 3-D surface seismic streamer have been acquired. A pre-injection survey was carried out in 1994 as the baseline data, followed by several post-injection survey in 1999, 2001, 2002, 2004, 2006, 2008 and 2010 as monitor data. Conventionally, seismic monitoring is performed in image domain involving pre-stack time migrated images of the baseline and monitor (Lumley 2001). This method assumes a linear relationship between baseline and monitor, both data sets being migrated with the same migration velocity model. It involves a time-shift correction on the monitor to match the interfaces of the baseline before 4-D analysis. In this case, some conditions where time lapse changes induce strong changes to the rock properties will violate the linear assumption. For instance, in gas breakout case can cause a sudden drop in seismic velocity and density, leading to significant time-shifts and amplitude variations. In addition, this conventional method is performed in the time-domain, which might not be straightforward and needs further time-to-depth conversion before interpretation, and it is prone to additional error due to conversion effects. Therefore, more straightforward quantitative monitoring analysis techniques in depth domain are necessary to give more robust time-lapse changes model.

Full waveform inversion (FWI) is a powerful method used to estimate subsurface parameters (Lailly 1983; Tarantola 1984, 1986). It involves formulating an optimization problem to minimize differences between field and synthetic data, starting with an initial model and iteratively refining it using linear approximations to the nonlinear problem. This method operates on the assumption that field data can be approximated by solving a wave equation, capturing the complex interactions between heterogeneous media and waves, including mode conversions, reflections, transmissions, refractions, multiple reflections and surface waves. FWI seeks to reconstruct subsurface physical properties by minimizing the misfit between predicted and recorded waveforms, offering higher resolution than traveltime tomography, down to half a propagated wavelength (Woodward 1992). Over the years, FWI has found applications at various scales, from global (Thrastarson et al. 2024) and regional seismology (Fichtner & Villaseñor 2015; Davy et al. 2017; Górszczyk et al. 2021; Nouibat et al. 2023) to exploration (Operto et al. 2015; Hicks et al. 2016) and shallow site characterization (Provenzano et al. 2018; Irnaka et al. 2022) both in marine and land settings. Recently, FWI has also been applied in medical imaging (Guasch et al. 2020). For time lapse applications, FWI inverts the baseline and monitor data to estimate their respective mechanical properties of the subsurface. In recent years, FWI has been applied to various 4-D problems from reservoir monitoring, leakage detection and CO₂ monitoring (Queißer & Singh 2013; Zhang & Huang 2013; Raknes & Arntsen 2014). The quantification of the 4-D changes can be directly estimated from the differences between monitor and baseline velocity models. In addition, reservoir properties such as saturation and porosity can be derived from the reconstructed FWI model (Dupuy et al. 2016).

Compared with the conventional workflow, FWI offers more direct quantitative estimation of the mechanical properties with high resolution. In general, only minimal pre-processing to the

data is needed for FWI. In contrary intense pre-processing steps must be applied to the data prior to the final time lapse processing in a conventional workflow. However, FWI is highly nonlinear, a careful workflow needs to be implemented to obtain optimal results. In time lapse problems, the inversion of two different data set independently can lead to convergence issues due to this nonlinearity, potentially creating 4-D artefacts. In addition, the non-repeatability of the acquisition between baseline and monitor is a challenging factor, as it is also a potential source of error and artefacts in 4-D FWI. Therefore, a careful time lapse strategy for FWI needs to be considered to handle potential 4-D problems.

In seismic monitoring, the cost of deployment of the surveys becomes a crucial factor, emphasizing the efficiency and cost effectiveness of highly repeatable monitor surveys. In general, contractors do not have the luxury to afford a monitor survey with the most advanced technology and the same large coverage as baseline survey, especially on carbon capture and storage projects where they operate more in a non-commercial manner. Therefore, in most cases it limits the monitor survey to be as cheap as possible, and more target oriented. In this case, the use of a survey with limited offset is more likely, as in a vintage data set acquired in previous years. This imposes another challenge to FWI, where the availability of long offset and low frequency is a key to an ideal imaging. Therefore, the workaround with the methodological aspect needs to be considered in time lapse application, especially since there is limited room for compromise on the acquisition side.

In this study, we consider a 3-D Sleipner vintage streamer data set, acquired over two decades ago, limited with short offset. This imposes a challenging condition for FWI to be able to produce meaningful model. In seismic data sets acquired with surface-deployed sources and receivers, limited source-receiver offsets and frequency bandwidth restrict FWI's retrieval of low-wavenumber information to shallow depths sampled by diving waves. Conversely, at deep targets accessible only by reflections, FWI exhibits less sensitivity to model kinematics, behaving similar to least-squares migration and reconstructing only the reflective component of the subsurface mechanical properties (Jannane et al. 1989; Mora 1989; Woodward 1992). This scenario highlights the necessity for an accurate starting model to cover the intermediate-wavenumber resolution gap.

Previously, Raknes et al. (2015) performed 3-D elastic FWI to the same data, demonstrating the difficulties to recover the deeper part of the model where it roughly matches with the well log profile despite meaningful shallow model reconstruction. They suggested the use of parametrization of P-wave velocity-impedance ($V_p-I_p$) instead of S-wave velocity density ($V_s-\rho$) could be more favourable for the short offset data. They also demonstrated the difficulties to reconstruct the S-wave velocity ($V_s$) model, whereas the inversion only provides minor updates in the upper part and quickly ran into local minimum, due to small amount of S-wave information in the field data set.

One way to circumvent the offset limitation while still relying on finite-frequency full waveform modelling approach is reflection waveform inversion (RWI). RWI uses prior reflectors information to create second-order scattering wavefield which interact with the first-order wavefield, producing low wavenumber contributions to the sensitivity kernel. In addition to RWI, Zhou et al. (2015) proposed joint FWI (JFWI) to enrich the RWI sensitivity kernel with the shallow contribution from diving waves, by
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Overview of the acquisition geometry. The blue dots are the source position, and the grey area are the receivers position. The red circle is the exploration well location.</figcaption></figure>

jointly fitting the early wave arrivals and reflections. This provides further robustness and stability of the overall updates on the model, where models kinematic at all depth are simultaneously improved.

In our work, we are interested to apply a reflection-oriented workflow in full time lapse mode to the limited offset Sleipner data for CO₂ monitoring. We split our works into two papers. The first paper focuses on the application on the 3-D baseline data set. The second paper focuses on the reflection-oriented workflow for time lapse FWI, its strategy and the application to the baseline and monitor data. In this first paper, we illustrate the application of this reflection oriented workflow to reconstruct the subsurface velocity and impedance model on the 3-D short offset data set as a prior step to construct a well-defined baseline model to be utilized for time lapse FWI. To do so, we rely on the methodology proposed by Provenzano et al. (2024) of velocity-impedance parametrization JFWI which consist of impedance model building with an asymptotic pre-conditioner (Li & Chauris 2018) for fast and balanced reconstruction, a pseudo-time reformulation of the JFWI gradient to ensure velocity-reflectivity consistency (Plessix 2013), and a graph space optimal transport (GSOT) objective function to mitigate cycle skipping (Métivier et al. 2016). We show that, this reflection oriented workflow is able to reconstruct a meaningful low-frequency model and high resolution impedance model which is consistent with the full-band migrated image. Subsequent run of standard FWI, starting from the JFWI model, is able to reconstruct a higher resolution model, providing better estimated broad-band model compared to the conventional FWI.

The paper is organized as follows: In Section 2 we describe the area where the field data set was acquired, in addition to details about the data set. In Section 3 we illustrate the application of the standard FWI on the field data, showing the limitation on the limited offset case. Section 4 is devoted to the inversion on the field data with reflection-oriented workflow, the result for JFWI approach and

the quality check. The discussions and conclusions are given in Sections 5 and 6, respectively.

## 2 OVERVIEW OF THE SLEIPNER DATA SET

The Sleipner area is located at the west of Norway in the North Sea. It consists of the Sleipner West gas field and the Sleipner East condensate field. The CO₂ is extracted from Sleipner West field and injected to the Sleipner East field into the 200–300 m thick Utsira formation. The Utsira formation lies between about 800–1100 m depth. The sandstone in Utsira formation is weakly cemented and has high porosity and permeability. On top of the Utsira formation, a thick and impermeable shale layer is acting as a seal.

The data set is a 3-D towed streamer survey from 1994 pre-injection. Originally, the data set was acquired with 3000 m cable length and a recording length of 5.5 s for each shot. The air gun source was towed at 6 m depth with 18.75 m shooting interval. The acquisition geometry is illustrated in Fig. 1. However, several processing steps have been applied to the data we have access to. They include: offset restriction to 1700 m, recording length reduction to 2.3 s, swell noise filtering, low cut filtering at 6.0 Hz, time-step resampling to 2.0 ms and $t^2$ amplitude scaling to compensate for geometrical spreading. These restrictions are applied to focus the investigation on the CO₂ injection area, avoiding the investigation on production area below. This means the diving waves are reaching only the shallow depth, while the target area at the deeper part only covered by the reflections. The absence of the low frequency due to pre-processing also imposes another challenge for FWI application in addition to the offset limitation.

The streamer data shown in Fig. 2 are dominated by unconverted PP arrivals (Queißer & Singh 2013). At near offset, approximately less than 400 m offset, the data are dominated by pre-critical reflections and multiples, while at far offsets it is dominated by refracted
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Common shot gather for different frequency band: (a) The full frequency band, (b) 6–10 Hz, (c) 6–12 Hz and (d) 6–14 Hz.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Amplitude spectrum of the full band data.</figcaption></figure>

arrivals, combination of head and turning waves, post-critical reflections and their ghosts and surface related multiple. The corresponding frequency spectrum is displayed in Fig. 3. The data are mainly distributed between 6 and 90 Hz, with some notches indicating the ghosts arrivals. Direct arrivals is hardly distinguished from other arrivals even at full frequency band. At lower frequency, the refracted and the direct arrivals are smeared together. As the frequency band is increased, refracted arrivals are easily seen and less smeared with the direct arrivals, showing more resolvable events.

Within Sleipner area, an exploration well was drilled. The corresponding sonic log data are shown in Fig. 4. However, the vertical profile of the P-wave velocity sonic log has no information above

250 m. Based on the sonic log, three high velocity layers are visible at approximately 600, 800 and 1100 m. In this case, the Utsira formation lies between the two high velocity layers at approximately 800 and 1100 m depth. On top of it, high velocity at 600 m correspond to the shale layers are sealing the Utsira formation. Missing from the sonic log data, on top of the shale layer, the overburden layers with a velocity range of 1600–1900 m s⁻¹ (Queißer & Singh 2013).

# 3 STANDARD FWI

Before we illustrate the application of our reflection workflow to the data set, we show the application of the conventional FWI, to illustrate the challenges and limitation of the short offset data application. We also describe the common inversion parameters and setup for the inversion workflows in this study.

## 3.1 Inversion Workflow

We perform inversion within a computational area of 5 km x 8 km x 2.3 km. The inversion is performed under the acoustic approximation with a time-domain finite-difference modelling engine to
<figure><img src="image_6.png" /><figcaption>Figure 4. Sonic well log from the inversion area.</figcaption></figure>

predict the calculated pressure data set. The finite-difference discretization is performed using a fourth-order in space and second-order in time staggered grid for the wave equation modelling (Yang et al. 2016). The models are discretized in space with h = 20 m for the starting 6–10 Hz band, to match the sea floor depth of 80 m and ensure at least five grid points per wavelength for the lowest velocity, and a time-step Δt = 2 ms. This value is selected to meet the Courant–Friedrichs–Lewy stability (CFL) conditions for our explicit finite-difference scheme. Hicks interpolation is used to interpolate source and receivers positions located in between grid points (Hicks 2002). A free surface condition is applied to the top boundary, to model ghosts and multiples in the data as they are still preserved in the data. Perfectly matched layer (PML) (Bérenger 1994) are applied to the other model boundaries to avoid spurious reflections and mimic an semi-infinite domain.

To reduce the computational cost, we decimate the field data by selecting one shot over three from the original number of shot, with a total of 2512 shots selected for the final data set. This increases the source interval from 18.75 to 37.5 m, which is still well below the Nyquist criteria for the reference wavelength. To further decrease the computational cost, we use a random subsampling strategy, with 256 sources per group of sources during the inversion (Warner et al. 2013). In order to recover the true amplitude of the data, we apply a reverse $l^2$ amplitude correction. In addition we also perform a minimum phase band-pass frequency filtering to the data. The

<figure><img src="image_7.png" /><figcaption>Figure 5. Estimated source wavelet at 6–10 Hz frequency band.</figcaption></figure>

inversion is performed in a multiscale manner (Bunks et al. 1995), from 6–10 to 6–14 Hz with 2 Hz frequency increment. At each inversion stage, a shot-averaged source wavelet is estimated using the frequency domain method proposed by Pratt (1999) (Fig. 5). The source wavelet is estimated using only the near offset part of the direct arrivals time-window.

The starting model is a 1-D linearly increasing velocity (Vp) as shown in Fig. 6. We use the Gardner’s relationship (Gardner et al. 1974) to derive the density model from the starting Vp model. This relationship is given as

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

The reason of utilizing a crude 1-D model in this study is to promote the reflection oriented workflow as an initial model building method for a challenging case related to the limited offset, compared to a more traditional reflection tomography.

## 3.2 FWI result

To investigate the effect of the short offset data restriction to the inversion result, we perform a standard $L^2$ mono-parameter FWI with a $V_p$ - $\rho$ parametrization, with the $\rho$ parameter kept fixed. No arrival separation through time-offset windowing is applied to the data to see the contribution from both transmission and reflections across the depth in the reconstructed model. We perform three iterations of 1-BFGS for each subgroup of sources, for a total of 36 iterations such that every source is used at least twice. The inversion is performed only on the starting frequency 6–10 Hz just to analyse the overall behaviour of the diving waves penetration limit.

We illustrate the limit of the diving waves penetration by showing the gradient in Fig. 7. The main low wavenumber perturbation is localized within a very shallow part, no deeper than 500 m. This indicates the diving waves penetration is unable to reach the expected target at 800 m depth. Consequently, the reconstructed $V_p$ model in Fig. 8 shows velocity updates which are localized only at the shallow part. The deeper part (below 500 m) is dominated by high wavenumber updates shown by the interfaces created. Comparison of 1-D profiles extracted at x = 1.2 km and y = 3 km location between initial $V_p$ and reconstructed $V_p$ in Fig. 9 show that there are no significant velocity update at the deeper part compared with the initial model. Nonetheless, we observe the results at shallow part are able to match the well log data quite well.

We compare the computed synthetic data based on the reconstructed FWI model with the observed data in mirrored view in Fig. 10. We show the data fit from two shot gathers located at the left part of the model and the right side of the model. We can observe that after FWI, the fit for the diving waves is quite satisfactory, indicating the velocity at the shallow part is accurate enough to explain the diving waves. However, we can see clearly the synthetic data are not able to predict most of the reflections on the observed data, only two strong reflections at 1.0 and 1.5 s being visible, and the fit is
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Starting 1-D *V*<sub>p</sub> model. Density is estimated using a Gardner's law relationship. Upper-left panel is a depth slices cutting through the volume at 200 m depth, bottom-left panel is a vertical slices at x = 1.9 km, and bottom right is a vertical slices at y = 3 km. Slices locations are indicated by the black dashed-lines.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** The first gradient of the conventional FWI starting from a 1-D model. The upper-left panel is a depth slice cutting through the volume at 200 m depth. The bottom-left panel is a vertical slice at x = 2 km. The bottom right panel is a vertical slice at y = 3 km. Slices locations are indicated by the black dashed lines.</figcaption></figure>

not satisfactory as well. This experiment illustrates the limitation of
conventional FWI on this short offset data, only the reconstruction
of the shallow part being satisfactory. Thus, we are going to use this
FWI as an additional part of a more general workflow through an
early waveform inversion, selecting only the diving waves part of
the data.

4 JFWI

4.1 Review of the method

Reflection oriented FWI aims at circumventing the limitation of
the lack of diving waves penetration due to narrow source–receiver
offset. In this reflection oriented workflow, we utilize the joint full
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Reconstructed $V_p$ after conventional FWI starting from 1-D model. The upper-left panel is a depth slice cutting through the volume at 200 m depth. The bottom-left panel is a vertical slice at $x = 2$ km. The bottom right panel is a vertical slice at $y = 3$ km. Slices locations are indicated by the black dashed lines.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** 1-D profiles extracted at well location. Comparison between initial model and after conventional FWI.</figcaption></figure>

waveform inversion (JFWI) to simultaneously invert the transmission waves at the shallow part and the reflected waves (Zhou et al. 2015), with the aim of enriching the low wavenumber component from the transmission at the shallow part with those from the reflection wave paths at the deeper part. This can be achieved by suppressing the migration isochrones through scale separation between the tomographic updates and its high wavenumber component from reflection, in combination with a velocity-impedance ($V_p-I_p$) parametrization.

We recall here the basic principle of JFWI. The objective function is computed as a weighted sum of diving waves misfit and reflections misfit compared with the calculated data. We have

$$C_{JFWI}(V_p, I_p) = \frac{1}{2} \|W^{\text{div}}(d_{\text{cal}}(V_p, I_0) - d^{\text{div}})\|_{L_2}^2 \\ + \frac{1}{2} \|W^{\text{refl}}(d_{\text{cal}}(V_p, I_p) - d^{\text{refl}})\|_{L_2}^2, \quad (2)$$

where $d_{\text{cal}}(V_p, I_0)$ is simulated based on the background velocity model $V_p$ and a smooth impedance model $I_0$ (typically derived from the initial model $V_p$ through a Gardner's law), while the second term $d_{\text{cal}}(V_p, I_p)$ is simulated using the background velocity model $V_p$ and the impedance model $I_p$. $W^{\text{div}}$ and $W^{\text{refl}}$ are data weighting operator to select and balance the contribution between early arrivals and reflections, respectively. During the inversion, diving and reflected waves contribute to the background velocity $V_p$ updates, while for $I_p$ updates only reflected waves are taken into account.

We minimize $C_{JFWI}$ using a variable projection approach. The gradient of $C_{JFWI}$ with respect to $V_p$ can be written in form as (Zhou et al. 2015)

$$\nabla_{V_p} C_{JFWI} = u_0(V_p, I_0) \star \lambda_0^{\text{div}}(V_p, I_0)(\text{banana} - \text{donuts}) \\ + \delta u(V_p, I_p) \star \lambda_0^{\text{refl}}(V_p, I_0) \\ + u_0(V_p, I_0) \star \delta \lambda^{\text{refl}}(V_p, I_p)(\text{rabbit ears}), \quad (3)$$

where $u_0(V_p, I_0)$ and $\lambda_0^{\text{div}}(V_p, I_0)$ denotes the incident wavefield and adjoint wavefield with the diving wave residuals as adjoint source propagating in the medium ($V_p, I_0$), and $\delta u(V_p, I_p)$ is the scattering component of the incident wavefield due to the presence of $I_p$. $\lambda_0^{\text{refl}}(V_p, I_0)$ denotes the adjoint wavefield with the reflection data residuals as the adjoint source inside the medium ($V_p, I_0$), and $\delta \lambda^{\text{refl}}(V_p, I_p)$ is the scattering component of the adjoint wavefield due to the presence of $I_p$. The $\star$ denotes the zero lag cross-correlation.
<figure><img src="image_12.png" /><figcaption>Figure 10. Data fit QC on common shot gather at (a) location x = 1 km and y = 3, and (b) x = 5 km and y = 4 km.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11. Construction of JFWI gradient extracted from Zhou et al. (2015). Gradient G1 is computed from the reflections part of the data, and G2 gradient is computed from diving waves and polarity-reversed reflections. The summation of G1 and G2 suppress the migration and isochrone and combine the low wavenumber contribution from diving and reflection waves.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 12. 3-D representation of the discrete graph of a reference Ricker function (red points) and a shifted in time Ricker function (blue points) scaled in amplitude by a factor 0.8. The grey arrows represent the assignment solution, which depends on the value of the parameter ΔT. (a) ΔT = 0.4 s, (b) tau = 4 s, (c) ΔT = 20 s [extracted from Métivier et al. (2019)].</figcaption></figure>

The gradient construction in (3) needs a reflectivity $I_p$ which we compute through Impedance Waveform Inversion (IPWI). This consists in minimizing $C_{JFWI}$ with respect to the impedance $I_p$. The gradient with respect to $I_p$ can be written as

$$ \nabla_l C_{JFWI} = u_0(V_p, I_p) \star \lambda_0^{\text{refl}}(V_p, I_p)(\text{migration isochrones}). \quad (4) $$

The gradient construction of JFWI can be illustrated in Fig. 11.

The original JFWI workflow needs cyclic re-estimation of the $I_p$ model through an $I_p$ waveform inversion along with the $V_p$ model updates, to ensure the structural consistency of the $I_p$ model. This is computationally demanding and therefore very time consuming. To overcome this issue, Provenzano et al. (2023) reformulate JFWI in the pseudo-time domain, following Plessix (2013) and Brossier et al. (2015). The gradient in eq. (3) is still computed in depth domain, but then reformulated from depth domain to pseudo-time

domain using the chain rule (Plessix 2013) as follows

$$ \nabla_{\tau(i)} C_{JFWI} = \nabla_{z(i)} C_{JFWI} - \int_{z_i}^{z_{\max}} \frac{dV_p^z}{dz} \frac{1}{V_p(z)} \nabla_z C_{JFWI} dz. \quad (5) $$

At each iteration, after $V_p$ is updated in the pseudo-time ($\tau$) domain, the $(V_p^\tau, I_p^\tau)$ model is transformed back to the depth ($z$) domain. By fixing $I_p$ in the pseudo-time rather than in the depth domain, we mitigate the velocity-depth ambiguity. In this case, the re-positioning of the $I_p$ model is implicitly done through time to depth domain conversion.

In our reflection workflow, we also replace the $L_2$ misfit function with a graph space optimal transport (GSOT) misfit function (Métivier et al. 2019; Métivier & Brossier 2021) to reduce the liability to cycle-skipping at intermediate-to-long offsets. In GSOT, seismic traces are transformed from oscillatory signals into 2-D distribution of unit-weight points in the time-amplitude graph-space domain illustrated as in Fig. 12. In this domain, the seismic
<figure><img src="image_15.png" /><figcaption>Figure 13. Data separation for inversion workflow. Area within red line is the arrivals used for source wavelet estimation, the area within blue line is used for impedance model reconstruction and green line separating the diving waves and reflections used in Vp reconstruction.</figcaption></figure>

traces are discretized into K samples, following Nyquist resampling. The goal is to find an optimal assignment between each observed ($t(j)$, $d_{\text{obs}}(j)$) and predicted ($t(i)$, $d_{\text{cal}}(i)$) sample. This assignment is achieved through an optimum permutation plan $\sigma^*$, which minimizes the effort required to transport one distribution ($d_{\text{cal}}$) into a target distribution $d_{\text{obs}}$ across all possible permutations. In the graph-space, the 2-Wasserstein distance between the two distributions is defined as follows:

$$h_{W^2} = \sum_{i=1}^{K} \frac{A^2}{\Delta T^2} |t_i - t_{\sigma^*(i)}|^2 + |d_i^{\text{cal}} - d_{\sigma^*(i)}^{\text{obs}}|^2, \quad (6)$$

$\Delta T$ is the maximum expected time-shift and $A$ is the maximum amplitude difference between observed and predicted data.

Convexity with respect to time-shifts is given by the monotonic dependence of the transport work on the terms $|t_i - t_j|^2$. Large $\Delta T$ values emphasize the convexity with respect to shifted patterns. As $\Delta t$ decreases, the mapping for each sample tends towards a local sample-by-sample mapping, and $h_{W^2}$ approximates the behaviour of the $L^2$ distance. In this case, $L^2$ can be viewed as a special case of GSOT distance when $\Delta T$ is close to zero.

In JFWI, the optimal transport distance of the reflection arrivals and diving wave arrivals can be computed independently, to honour the differences in amplitude range between the two propagation regimes. The objective function in eq. (2) can be re-written using eq. (6) as:

$$C_{\text{JFWI}}(V_p) = \sum_{\text{shot}} [W^e h_{W^2}^r(\sigma^{\text{div}*}) + W^r h_{W^2}^e(\sigma^{\text{refl}*})]. \quad (7)$$

In this case, each propagation regime is scaled by $A_{\text{div}}$ and $A_{\text{refl}}$ respectively, with the possibility of defining different maximum expected time-shifts $\Delta t$. For example, a larger one is used for diving waves than for reflections in Provenzano et al. (2024).

## 4.2 JFWI on Sleipner data set

As part of the JFWI workflow, separation between transmission arrivals and reflection arrivals is crucial. To design an appropriate window for the separation, synthetic data are computed with the starting model to identify the time-offset separation range between transmission and reflection. A taper zone of 0.15 s is applied to avoid contamination between arrivals. The windowing function is used every shot and re-estimated for each frequency band. The taper zone is also reduced as the frequency go higher. Separation of different portions of data for source wavelet estimation, impedance model building and velocity reconstruction are illustrated in Fig. 13.

Prior to JFWI, we perform an early wave inversion, a conventional FWI but only using the diving waves information. The reconstructed $V_p$ model is shown in Fig. 14. It shows that conventional FWI is able to provide a quite accurate model at the shallow part despite limitation at the deeper part, consistently with the conventional FWI results presented in Section 3. The inversion parameters for early wave inversion are identical with previous ones, except that we mute the reflected waves. The motivation of performing this early wave inversion is to better constrain the shallow part, in order to promote better accuracy and positioning of the $I_p$ model during the impedance model building in the next step.

We start the reflection oriented workflow by building the impedance through impedance waveform inversion (IPWI), windowing the short offset reflections up to 600 m of offset, and using an asymptotic pre-conditioner as in Provenzano et al. (2024). Three iterations are performed for each group of 256 sources batch, for a total of 36 iterations, so that each available shot gather is used at least once. The reconstruction of the $I_p$ model is crucial in the overall JFWI workflow as it provides prior reflectors information for velocity model building. The reconstructed $I_p$ model is presented in Fig. 15, showing dominance of layered interfaces, as expected from the Sleipner area which has relatively weak lateral impedance variations (Queißer & Singh 2013). Particularly interesting channel structure at 200 m is delineated quite well from the $I_p$ model as shown in the depth slice in Fig. 15.

We illustrate the benefit of JFWI on short offset data by presenting the first JFWI gradient in Fig. 16. It is a combination of the tomographic component from the diving waves and reflection waves. Clearly, the JFWI is able to promote a deeper reach of the gradient coverage on the deeper part, taking advantage of the tomographic component from the reflections. Compared with the conventional FWI, the JFWI gradient is able to constrain down to 2 km instead of 500 m with conventional FWI.

Utilizing the reconstructed $I_p$ model, a subsequent run of JFWI is performed with full offset for 72 pre-conditioned steepest descent (ST) iterations, three iterations each for 256 shot random resampling batch. The choice of ST instead of l-BFGS is to ensure more robust updates in reflection-based inversion, despite its slower convergence compared with the l-BFGS. This is because the presence of random and coherent (e.g. ghosts and multiples) noise may accumulate inaccurate information about the curvature of the objective function from the previously computed gradients.

The weight of the reflected waves term is two times as large as the one of diving waves, to balance the contribution between the diving and reflection waves. A depth dependent gain is also applied to the JFWI gradient to enhance deep velocity updates. The GSOT maximum expected shift ($\Delta T$) is set to 0.2 s for both diving and reflection waves by examining the phase difference between synthetic data generated from the initial model and the observed data. In pseudo-time domain formulation JFWI, the position of the
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Reconstructed $V_p$ after an early wave inversion. The upper-left panel is a depth slice cutting through the volume at 200 m depth. The bottom-left panel is a vertical slice at x = 2 km. The bottom right panel is a vertical slice at y = 3 km. Slices locations are indicated by the black dashed-lines.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Reconstructed $I_p$ after impedance waveform inversion. The upper-left panel is a depth slice cutting through the volume at 200 m depth. The bottom-left panel is a vertical slice at x = 2 km. The bottom right panel is a vertical slice at y = 3 km. Slices locations are indicated by the black dashed-lines.</figcaption></figure>

$I_p$ model interface is updated simultaneously in depth along with the $V_p$, therefore avoiding re-estimation of the $I_p$ through IPWI to maintain consistency between kinematics and reflectivity at short offset. However, we re-estimate the $I_p$ at the beginning of each frequency band starting from estimated $V_p$ from previous frequency band.

The final JFWI model is shown in Fig. 17. Compared with the model from the conventional FWI, the JFWI model is not only able to provide a meaningful update of the shallow part, but also of the deeper part greater than 500 m of depth. In addition, cor- related with the $I_p$ model, the reconstructed $V_p$ model is able to
<figure><img src="image_18.png" /><figcaption>**Figure 16.** The first gradient of JFWI starting from a 1-D model. The upper-left panel is a depth slice cutting through the volume at 200 m depth. The bottom-left panel is a vertical slice at x = 2 km. The bottom right panel is a vertical slice at y = 3 km. Slices locations are indicated by the black dashed-lines.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 17.** Reconstructed $V_p$ after JFWI workflow. The upper-left panel is a depth slice cutting through the volume at 200 m depth. The bottom-left panel is a vertical slice at x = 2 km. The bottom right panel is a vertical slice at y = 3 km. Slices locations are indicated by the black dashed-lines.</figcaption></figure>

delineate the channel structure which correspond to the lower velocity value compared to the surrounding. As expected, the $V_p$ model from JFWI is very smooth, comparable with the model one might obtain through tomography. In this case, the inversion will then proceed with a conventional FWI step starting from the final JFWI $V_p$ model to introduce a higher wavenumber content to the model.

### 4.3 FWI after JFWI

We perform a subsequent run of conventional FWI starting from the JFWI $V_p$ model with full offset up to 1.7 km with a linear depth pre-conditioning and three iterations per group of sources for a total of 36 iterations. The inversion is started from the lowest frequency band and continue in a multiscale manner. The parametrization is changed from [$V_p$, $I_p$] to [$V_p$, density ($\rho$)] (Operto 2013), where $\rho$
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Reconstructed $V_p$ after JFWI + FWI. The upper-left panel is a depth slice cutting through the volume at 200 m depth. The bottom-left panel is a vertical slice at x = 2 km. The bottom right panel is a vertical slice at y = 3 km. Slices locations are indicated by the black dashed-lines.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 19.** Data fit QC on common shot gather at location x = 1 km and y = 3 km: (a) Initial, (b) after JFWI workflow and (c) after JFWI + FWI workflow.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 20.** Data fit QC on common shot gather at location x = 5 km and y = 4 km: (a) Initial, and (b) after JFWI workflow and (c) after JFWI + FWI workflow.</figcaption></figure>

is reset to its initial value because density after JFWI in the [Vp, Ip] parametrization has no physical meaning. FWI is performed for 48 1-BFGS iterations (Nocedal 1980), three iterations each for 256 shot random subsampling batch. GSOT misfit function is still utilized to mitigate potential cycle skipping with the maximum expected time-shift set to 0.2 s.

The results of FWI following JFWI are shown in Fig. 18. The contribution of JFWI to FWI is apparent at all depths by introducing higher wavenumber contribution from the reconstructed JFWI model, in particular we can see more distinct delineation of the channel structure at 200 m compared to the reconstructed JFWI
<figure><img src="image_23.png" /><figcaption>**Figure 21.** 1-D profiles extracted at well location: (a) Comparison between initial model and after early wave FWI, (b) comparison between initial model and after JFWI and (c) Comparison between the well log profile with all reconstructed vp profile from different workflow.</figcaption></figure>

model. In contrast, the *V*<sub>p</sub> model from conventional FWI fails to update the model away from the 1-D starting model at those depths not sampled by diving waves.

## 4.4 Quality checks

To check the quality of reconstructed models from JFWI, we compare the data fit between synthetic and the observed data. Two shot gather located at the left part and the right part of the model are selected as shown in Figs 19 and 20. A mirrored view between synthetic and observed data is presented to showcase the data fit for: (a) The starting model, (b) The JFWI model and (c) the JFWI + FWI model. As shown in Figs 19(a) and 20(a), the initial model is unable to predict the diving waves accurately, in particular at long offset. In comparison, the results obtained with JFWI presented in Figs 19(b) and 20(b) are able to fit both diving waves and reflections waves quite nicely. Only marginal improvements in the data fit for the JFWI + FWI model presented in Figs 19(c) and 20(c) are observed, mostly for large offset reflections. In particular, reflection amplitude prediction is improved in synthetic data based on the JFWI + FWI model.

Fig. 21 shows a comparison between the smoothed sonic log from the exploration well and the reconstructed models. JFWI is able to estimate a *V*<sub>p</sub> model which matches the low wavenumber trend of the well log data. In addition, we can see meaningful contribution of reflection on JFWI to perturb the deeper part of the model, decreasing the velocity closer to the expected well log data. Compared with the JFWI result, the JFWI + FWI results shows more high wavenumber component which match quite well with the well log data.

To check the structural consistency, we compare the reconstructed structure from *V*<sub>p</sub> model and *I*<sub>p</sub> model with the full band post-stack time migrated (PSTM) image. The PSTM image is converted to depth with the smooth velocity from JFWI. In particular, we compare the channel structure at 200 m as shown in Fig. 22. We can observe clearly the channel structure on *I*<sub>p</sub> model (Fig. 22b) has higher resolution compared with the *V*<sub>p</sub> one (Fig. 22a) due to the scale separation. In comparison, delineation of the channel structure is consistent with the migrated image.

## 5 DISCUSSION

Application of conventional FWI on the limited offset Sleipner data set restricts the diving waves penetration at a maximum of 500 m depth, constraining the velocity updates only at shallow depth. Despite its limitation at the deeper part, the reconstructed *V*<sub>p</sub> at the shallow part is quite satisfactory, which motivates the use of early waves only inversion as and additional part in the workflow to better constrain the shallow part. The main motivation of this additional step in the workflow is to supply a better constrained shallow part model for the *I*<sub>p</sub> inversion, such that the interfaces at the shallow part of the *I*<sub>p</sub> model could be better estimated. Fig. 23 illustrates the Δ*I*<sub>p</sub> model. This is calculated by subtracting the final impedance model from the inversion and the starting impedance model to better emphasize the interfaces. It shows that the interfaces at the shallow part in the Δ*I*<sub>p</sub> model starting from *V*<sub>p</sub> after early wave inversion is more continuous compared with the one from initial 1-D model.

Although the shallow part of the model is well constrained by the diving waves through early wave inversion, the application of JFWI is still preferred instead of RWI. In areas of the subsurface sampled by reflections and refractions, the latter stabilize the inversion results in data containing surface-related multiples and elastic effects (Provenzano et al. 2024). This makes it possible to invert the data with minimal pre-processing, that is, no multiple suppression or dehosting, and within the acoustic approximation, despite expected elastic effects (e.g. amplitude variation with offset).

Ideally, FWI should be implemented with the most accurate physics described by the elastodynamic equation, to accommodate the complete wave propagation regimes. However, in Sleipner data the *S*-wave velocity is difficult to be constrained by the inversion. This can be explained by the lack of *S*-wave information in the data (Raknes et al. 2015). In this case, the use of multicomponent data might be necessary to successfully invert the *S* waves. Despite the data are dominated by unconverted PP arrivals, the mode conversions (e.g P-S-P waves) are still expected to be present in the data. The impact of converted P-S-P waves may challenge the acoustic approximation used in JFWI. However, usually these mode conversions occur in a longer propagation path, therefore it is expected to be recorded at late time in the seismic recording. In Sleipner
<figure><img src="image_24.png" /><figcaption>**Figure 22.** Comparison of the reconstructed channel at depth of 200 m for different model. (a) Reconstructed channel from the Vp model, (b) reconstructed impedance model and (c) post-stack migration image derived from full band data.</figcaption></figure>

data, the seismic data recording length has been reduced to 2.3 s from the original 5.5 s, implying less converted waves should be recorded within this time window. Nevertheless, there should be residual energy associated with the elastic effects which cannot be explained by the acoustic approximation. In the other hand, the elastic reflection waveform inversion is still at the research stage, with some notable development by Masmoudi et al. (2021) and Aaker et al. (2024). Masmoudi et al. (2021) takes advantage of the PS Born modelling to construct the PS sensitivity kernel while Aaker et al. (2024) suggests a PS mode-decomposition built into the misfit functional used for elastic reflection waveform inversion.

The results in this study are comparable with the one obtained in Queißer & Singh (2013), both in the overburden area and the

deeper part. In their study, they perform 2-D elastic FWI to the single line of Sleipner data up to 38 Hz frequency, and having access to the longer maximum offset of 3500 m compared with this study. Compared with their study, our results are not able to reconstruct quite well the high velocity at relatively shallow depth of 600 m, corresponds to the shale layer. This can be caused by the inaccurate density and unexplained elastic effects. Better-constrained shallow densities (e.g. aided by density logs) and additional pre-processing (e.g. surface-related multiple attenuation) may in principle reduce the shallow-model inaccuracies observed in this case study and thus their impact on the quality of the reconstructed model. Compared with the results obtained by Raknes et al. (2015) on the same data
<figure><img src="image_25.png" /><figcaption>**Figure 23.** Comparison of the reconstructed $\Delta I_p$ model starting from difference initial $V_p$ model. (a) Vertical profile extracted at $x = 2$ km, and (b) vertical profile extracted at $y = 3$ km are the $\Delta I_p$ model starting from 1-D $V_p$ model. (c) Vertical profile extracted at $x = 2$ km, and (d) vertical profile extracted at $y = 3$ km are the $\Delta I_p$ model starting from $V_p$ model after early wave inversion.</figcaption></figure>

set, our results shows clear improvement at the deeper part, where their model is unable to match the well log data below 800 m.

In our study, we use graph-space optimal transport (Métivier et al. 2019) misfit function in the framework of JFWI to mitigate cycle skipping on the data and produce robust inversion result against inaccurate starting models. Alternatively, we could also use another misfit function recently proposed by Yong et al. (2022) which implicitly measure the instantaneous time-shift between compared data through a local matching filter, called localized adaptive waveform inversion (LAWI). In combination with the JFWI, the phase-driven LAWI misfit function can balance the contributions of early arrivals and reflections to the JFWI gradient, taking benefit of the simultaneous reconstruction of shallow and deep velocity trends (Yong et al. 2024).

## 6 CONCLUSION

We have illustrated an application of the reflection oriented workflow full waveform inversion to the limited offset 3-D Sleipner data. We show the limitation of conventional FWI applied to the limited offset data, constraining only the shallow subsurface while the deeper part is not well updated. Taking advantage of the conventional FWI with diving waves only, we better constrain the velocity model at the shallow part with an early wave inversion, ensuring better input for the impedance estimation. Starting from the velocity model from the early wave inversion, we estimate the impedance model to serve as a prior reflectors information for JFWI. The result of JFWI shows clear improvement of the velocity model compared to the conventional FWI, especially at the deeper part, it is also able to match the low frequency trend of the well log data quite well and produces satisfactory fit with the observed data. A subsequent run of conventional FWI starting from the reconstructed velocity model from JFWI is able to introduce higher wavenumber content to the velocity model and match the events with the well log data, and show similar satisfactory fit with the

observed data compared with the JFWI model. The reconstructed baseline model in this study will be served as the starting point for an investigation on time-lapse FWI strategies in the second paper.

## ACKNOWLEDGMENTS

This study was partially funded by the SEISCOPE consortium (http://seiscope2.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. Bhargav Boddupalli, Bastien Dupuy and Anouar 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, Hafllund Oslo Celsio, KROHNE, Larvik Shipping, Norcem Heidelberg Cement, Offshore Norge, Quad Geometrics, Stratum Reservoir, TotalEnergies, 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 SEISCOPE consortium. Any request related to the materials availability needs to be addressed to Romain Brossier
(romain.brossier@univ-grenoble-alpes.fr) or Ludovic Metivier (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, 2024*(1), 1–5.

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

Brossier, R., Operto, S. & Virieux, J., 2015. Velocity model building from seismic reflection data by full waveform inversion, *Geophys. Prospect.*, **63**, 354–367.

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

CCS Institute, Global., 2023. The global status of ccs: 2023, Technical Report. Australia, Global CCS Institute.

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

Dupuy, B., Asnaashari, A., Brossier, R., Garambois, S., Métivier, L., Ribodetti, A. & Virieux, J., 2016. A downscaling strategy from FWI to microscale reservoir properties from high-resolution images, *Leading Edge*, **35**, 1146–150.

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

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.

Gardner, G. H.F., Gardner, L.W. & Gregory, A.R., 1974. Formation velocity and density—the diagnostic basics for stratigraphic traps, *Geophysics*, **39**, 770–780.

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 1d model, *J. geophys. Res.: Solid Earth*, **126**(5), e2020JB021504.

Guasch, L., Agudo, O.C., Tang, M.-X., Nachev, P. & Warner, M., 2020. Full-waveform inversion imaging of the human brain, *NPJ Digital Med.*, **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., Metivier, L., Bohlen, T. & Pan, Y., 2022. 3D Multi-component full waveform inversion for shallow-seismic target: Ettlingen line case study, *J. geophys. Int.*, **229**(2), 1017–1040.

Jannane, M. et al., 1989. Wavelengths of Earth structures that can be resolved from seismic reflection data, *Geophysics*, **54**(7), 906–910.

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

Li, Y. & Chauris, H., 2018. Coupling direct inversion to common-shot image-domain velocity analysis, *Geophysics*, **83**(5), R497–R514.

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

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., 2021. New insights on the graph space optimal transport distance for full waveform inversion, in *SEG Technical Program Expanded Abstracts 2021*. Society of Exploration Geophysicists, 812–816.

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 Probl.*, **35**(8), 085001.

Métivier, L., Brossier, R., Mérigot, Q., Oudet, E. & Virieux, J., 2016. Measuring the misfit between seismograms using an optimal transport distance: application to full waveform inversion, *J. geophys. Int.*, **205**, 345–377.

Mora, P.R., 1989. Inversion = migration + tomography, *Geophysics*, **54**(12), 1575–1586.

Nocedal, J., 1980. Updating quasi-Newton matrices with limited storage, *Math. Comput.*, **35**(151), 773–773.

Nouibat, A., Brossier, R., Stehly, L., Cao, J., Paul, A., Team, Cifalps. & AlpArray Working, Group., 2023. Ambient-noise wave-equation tomography of the alps and ligurian-provence basin, *J. geophys. Res.: Solid Earth*, **128**(10), e2023JB026776.

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, *J. geophys. Int.*, **202**(2), 1362–1391.

Operto, S., Gholami, Y., Prieux, V., Ribodetti, A., Brossier, R., Metivier, L. & Virieux, J., 2013. A guided tour of multiparameter full-waveform inversion with multicomponent data: From theory to practice, *The Leading Edge*, **32**(9), 1040.

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

Pratt, R.G., 1999. Seismic waveform inversion in the frequency domain, part I: theory and verification in a physical scale model, *Geophysics*, **64**, 888–901.

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, *J. geophys. 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.

Provenzano, G., Brossier, R. & Métivier, L., 2024. Robust and efficient waveform-based velocity-model building by optimal transport in the pseudotime domain: An ocean-bottom cable case study in the North Sea, *Geophysics*, **89**(2), B83–B103.

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

Raknes, E.B. & Arntsen, B., 2014. Time-lapse full-waveform inversion of limited-offset seismic data using a local migration regularization, *Geophysics*, **79**(3), WA117–WA128.

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

Tarantola, A., 1984. Inversion of seismic reflection data in the acoustic approximation, *Geophysics*, **49**(8), 1259–1266.

Tarantola, A., 1986. A strategy for non linear inversion of seismic reflection data, *Geophysics*, **51**(10), 1893–1903.

Thrustarson, 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. Full-waveform inversion of cycle-skipped seismic data by frequency down-shifting, chap. 176, pp. 903–907, Society of Exploration Geophysics.

Woodward, M.J., 1992. Wave-equation tomography, *Geophysics*, **57**, 15–26.

Yang, P., Brossier, R., Métivier, L. & Virieux, J., 2016. Wavefield reconstruction in attenuating media: A checkpointing-assisted reverse-forward simulation method, *Geophysics*, **81**(6), R349–R362.
Yong, P., Brossier, R., Métivier, L. & Virieux, J., 2022. Localized adaptive waveform inversion: theory and numerical verification, *J. geophys. Int.*, **233**(2), 1055–1080.

Yong, P., Provenzano, G., Brossier, R. & Métivier, L., 2024. Velocity model building by localized adaptive waveform inversion of early arrivals and reflections, *85th EAGE Annual Conference & Exhibition, 2024*(1), 1–5.

Zhang, Z. & Huang, L., 2013. Double-difference elastic-waveform inversion with prior information for time-lapse monitoring, *Geophysics*, **78**(6), R259–R273.

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