# High-resolution multiparameter characterization of the subsurface using full waveform inversion on broad-band data: application to the oceanic crust in the North Sea using a dense ocean bottom cable data set

Ludovic Métivier, Romain Brossier, Andrzej Górszczyk, Giuseppe Provenzano, Serge Sambolian, Alizia Tarayoun

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

To cite this version:

Ludovic Métivier, Romain Brossier, Andrzej Górszczyk, Giuseppe Provenzano, Serge Sambolian, et al.. High-resolution multiparameter characterization of the subsurface using full waveform inversion on broad-band data: application to the oceanic crust in the North Sea using a dense ocean bottom cable data set. Geophysical Journal International, 2025, 243 (2), pp.ggaf352. 10.1093/gji/ggaf352 hal-05336536

HAL Id: hal-05336536

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

Submitted on 29 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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
# High-resolution multiparameter characterization of the subsurface using full waveform inversion on broad-band data: application to the oceanic crust in the North Sea using a dense ocean bottom cable data set

L. Métivier,<fnref n="1" />,<fnref n="2" /> R. Brossier,<fnref n="1" /> A. Górszczyk,<fnref n="1" />,<fnref n="3" /> G. Provenzano,<fnref n="1" /> S. Sambolian,<fnref n="1" />* and A. Tarayoun<fnref n="1" />

Accepted 2025 August 26. Received 2025 August 21; in original form 2025 April 14

## SUMMARY

This study focuses on the hydrophone component of a dense ocean bottom cable data set from the North Sea. This data had already been used in the past to illustrate the high resolution power of full waveform inversion based strategies. We have developed a highly scalable implementation of a visco-acoustic full waveform inversion engine making it possible to double the frequency content of the inverted data compared to previous studies, using simultaneously up to almost 50 000 CPU. This results in a multiparameter reconstruction of the subsurface, where the P-wave velocity, the density and the quality factor are reliably reconstructed down to 2 km depth, with a resolution of about 10 m.

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

## 1 INTRODUCTION

Since the industrial revolution at the turn of the 19th century, the development of our societies has been sustained by fossil energy sources found in the first kilometres of the crust. These are still the main sources for producing electricity and powering fret, public and private transport, with massive emission of greenhouse gases, such as carbon-dioxide (CO₂), in the atmosphere. Today, the transition to alternative sources of energy is key to mitigate the impact of human development on global climate change. To make this transition possible, a detailed knowledge and characterization of the properties and processes in the top first kilometres of the Earth’s crust plays a paramount role.

Indeed, harnessing wind and solar energy relies on rare-earth ore resources required to produce batteries. Exploration of uranium is a key to produce nuclear energy. Building the related nuclear plants requires a precise geo-technical knowledge of the top hundred metres of the subsurface. The crust itself offers alternative renewable energies. Earth’s natural heat can be exploited to power electrical

plants and feed urban heating systems. Enhanced Geothermal systems (EGS), relying on fluid injection and circulation, extend the scope of geothermal energy beyond natural hydrothermal sites, requiring a detailed structural characterization (Salaun *et al.* 2020). Naturally occurring white hydrogen has started to be explored recently, with promising advances in the knowledge of its production and location in the shallow crust (Malvoisin & Brunet 2023).

Along with a transition to renewable energies, the capture of greenhouse gases (especially CO₂) and their storage in viable geological sites will be crucial (Lee *et al.* 2023), as emphasized in the latest IPCC reports. Only the shallow crust, for instance in depleted gas reservoirs, or in saline aquifers, can provide sufficiently large confined volumes to host efficiently supercritical CO₂. Aside permanent CO₂ sequestration, the temporary storage of gas dedicated to energy production as CH₄ or H₂, with injection and pumping phases to accommodate supply and demand flow rates, can also take benefit from natural reservoirs hosted in the subsurface if coupled with the required monitoring technology. Similar challenges are posed in the context of nuclear waste disposals monitoring.

To answer the need of this detailed subsurface characterization, for mineral exploration prospecting, safety monitoring of CO₂ sequestration sites and geothermal energy systems, advances in geo-physical data processing is required. In particular, seismic imaging

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

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

<footnote n="3">³Institute of Geophysics, Polish Academy of Sciences, Warsaw, Poland</footnote>

<footnote>*Now at: Univ. Strasbourg, CNRS, Institut Terre et Environnement de Strasbourg, UMR 7063, F-67084 Strasbourg, France</footnote>
approaches play a key role. Conventional seismic imaging ‘migrates’ the mechanical waves recorded at the surface back to the reflective physical discontinuities in the subsurface (Claerbout 1985) corresponding, for example, to changes in lithology or pore-fluid content. The reliability of the obtained image depends on the accuracy of a low resolution wave velocity model determined in a prior step through traveltime tomography (Nolet 1987). This two-step workflow relies on a time-consuming and rather subjective series of data processing steps to isolate reflections from unwanted energy and ‘picked’ traveltimes. The output reflectivity image is interpreted in terms of subsurface structure and lithology. However, to reconstruct and identify time-lapse changes of subsurface parameters related to fluid migration pathways and accumulation, or discriminate between gas and liquid phases, a 3-D multiparameter medium characterization (i.e. seismic wave velocities, density and attenuation) is required with a resolution in the order of few tens of metres (Raknes et al. 2015).

Since 2010 (Sirgue et al. 2010), Full Waveform Inversion (FWI) (Lailly 1983; Tarantola 1984; Virieux & Operto 2009; Virieux et al. 2017) has progressively reached maturity as an alternative all-att Once imaging approach, ideally achieving the characterization of short- (reflective) and large-scale (kinematic) features of the subsurface mechanical properties from almost raw data. FWI iteratively updates 3-D maps of subsurface mechanical properties, by minimizing the difference between computed and recorded seismic data. This has the potential to reduce the turnaround time from acquisition to characterization all the while accounting for the complexity of the full recorded wavefield.

FWI draws its imaging power from a quantitative, physics-based interpretation of the full wavefield-structure interactions, that is, reflections, refractions, absorption and mode conversions. Nevertheless, due to the strongly nonlinear mapping from the observed waveforms to the subsurface mechanical parameters, exploiting this wealth of information comes at the price of solving an ill-posed nonlinear inverse problem, requiring dedicated regularization strategies. A second difficulty comes from the computational cost of solving the thousands of 3-D wave propagation problems for the iterative model update. In the theory of diffraction tomography, the upper limit of spatial resolution equals half the minimum propagated wavelength (Devaney 1984; Wu & Toksöz 1987). Therefore, accessing high-resolution features requires increasing the data bandwidth. However, the number of floating-point operations to perform to compute a wavefield depends on the power 4 of the data corner frequency, namely, doubling the corner frequency multiplies by 16 the computational cost.

For these reasons, FWI has been often restricted to low frequency data (below 10 Hz), in areas where the subsurface structure is mildly complex. The main output from FWI is in this case a refined velocity model, which is used to improve the quality of the reflectivity image. More recently, industry applications, leveraging on an increased available computing power, have shown how FWI could be pushed towards higher frequency bands, reaching frequencies in general used for the migration step (Huang et al. 2021). The directional derivatives of the reconstructed high resolution velocity models are then used as a proxy for reflectivity models, named FWI-image. The latter shows a broader wavenumber content than reflectivity models derived from more conventional (least-squares) reverse time migration. The advantage of such strategy is thus two-fold: an accelerated data-to-image workflow in an integrated single step, and a broader wavenumber high quality final reflectivity image. The latter result shall not be a surprise: FWI in high frequency

regime works as a nonlinear least-squares migration algorithm, accommodating mainly amplitude information, the kinematic of the wave propagation having been interpreted in the lower frequency bands (Mora 1989).

In this context, we present in this study a high resolution application of FWI to a dense Ocean Bottom Cable seismic data acquired in the North Sea, with the following two main contributions. First, we question the recurring practice of inverting for a single velocity parameter from which the reflectivity is derived in a second step. Multiparameter FWI is challenging, especially for the simultaneous reconstruction of the density and the quality factor governing the attenuation. However, from broad-band data, for subsurface zone illuminated by a sufficient diversity of waves, offering a continuity from large-angle transmission to narrow-angle reflection data, FWI has in theory the potential to provide reliable multiparameter information (Operto et al. 2013). Here we assess this strategy within the framework of visco-acoustic FWI, inverting an OBC hydrophone data set from the North Sea, aiming at reconstructing simultaneously the *P*-wave velocity, the density and quality factor *Q*<sub>P</sub> associated with the attenuation of *P* waves.

We reconstruct from this data 3-D maps of the *P*-wave velocity and density models, which exhibit high resolution features that we can interpret geologically. In particular, near sea-bed paleo-channels are exhibited with decametric resolution. A fine layering structures appears in the shallow part of the *P*-wave velocity model which is consistent with the known stratigraphy in this part of the North Sea. The pseudo-reflectivity model extracted from the reconstructed *P*-impedance makes it possible to highlight the main reconstructed structures. More strikingly, the fit to the *P*-wave velocity and density logs show that our FWI algorithm, taking benefit from the large frequency band and the diversity in illumination, is capable of reconstructing quantitatively these two parameters without trade-off down to 2 km depth. Finally, the reconstructed quality factor is shown to consistently attenuate the reflectivity in the central low-velocity region, while helping in fitting the observed data.

The second main contribution of this study is to show that we were able to apply FWI in a range of frequencies and a final resolution comparable to the ones considered currently by industry research groups. The frequency content we reach, namely 22.5 Hz, remains lower to what is presented in industrial case studies where inversion of the data beyond 50 Hz is reported (see for instance Kalinicheva et al. 2021; Zhang et al. 2022; Espin et al. 2023; Hirst et al. 2024) but is already twice as high as previously reported FWI results from academic group on the same data set (Operto & Miniussi 2018). The results reported here have been obtained thanks to the access to the French national machines. To give few numbers, finalizing this study required approximately 28 millions Central Processing Units (CPU) hours on the new AMD-GENOA partition of the national machine ADASTRA, in CINES, France, with FWI jobs running on 49152 CPU for the largest frequency bands. The ability to develop efficient multiparameter FWI algorithms running on national machines is crucial because high quality seismic data from dense networks are now more and more accessible, to help in earthquakes detection (Wu et al. 2013), or study complex targets such as the Alps (Hetényi et al. 2018), active faults (Brenguier et al. 2023), subduction zones (Górszczyk et al. 2025), volcanic areas (Aubert et al. 2025), geothermal sites (Lehujeur et al. 2018) and even glacier (Gimbert et al. 2021) to give few examples. The quantitative, 3-D, multiparameter inversion of such data set by FWI will be key to improve our knowledge of the Earth's crust and the physical process in action within it.
The structure of the study is as follows. We give in the first section an presentation of our FWI methodology, relating the different ingredients we have developed and how they are implemented. In the second section, we present the North Sea data set we consider, the workflow we design to invert for it and the high resolution reconstruction of the P-wave velocity, the density and the quality factor we achieve. A discussion section follows, before we draw the conclusion from this study.

## 2 MULTIPARAMETER FWI METHODOLOGY

FWI has been conceptualized in the 80s by Lailly (1983) and Tarantola (1984). The specific focus of our FWI implementation is on the mitigation of coupling problems in multiparameter reconstruction (Operto et al. 2013) and the development of a computationally efficient parallel algorithm. For the multiparameter aspect, we combine a pre-conditioned I-BFGS optimization strategy (Nocedal 1980) with a unity-based normalization, and a parametrization where we invert for the vertical velocity, the density and the inverse of the quality factor (Yang et al. 2018). Related to the computational efficiency, we combine a shot subsampling strategy (Warner et al. 2013) with an implementation exploiting domain decomposition at all stages of the algorithm, whether it is for the computation of visco-acoustic wavefields, the regularization using anisotropic diffusion (Métivier & Brossier 2022) or the nonlinear optimization operations (Métivier & Brossier 2016). The gradient computation implementation is versatile and can be adapted depending on the available memory resources, from basic storing techniques to more evolved stabilized checkpointing-assisted recomputation techniques (Yang et al. 2016). In the following we provide a detailed view of this multiparameter inversion algorithm.

### 2.1 Mathematical formalism

The observed seismic data $d_{obs}$ is denoted as a function depending on the active seismic source position $x_s$, the receiver positions $x_r$ and time $t$ as $d_{obs} = d_{obs,s}(x_r, t)$ for $1 \le s \le N_s$, $1 \le r \le N_r$ and $t \in [0, T]$ with $N_s$ and $N_r$ the number of active seismic sources and receivers, respectively, while $T$ is the total recording time (usually to the order of few seconds for crustal imaging experiments).

FWI is based on the iterative matching between these observed data and the calculated data $d_{cal,s}(x_r, t)$. To compute this quantity we rely first on a mathematical operator representing the wave propagation within the subsurface. We adopt a general operator notation, where the wavefield $u_s$ is defined as the solution of the wave equation

$$A(m)u_s = b_s, \qquad (1)$$

where $A(m)$ is the wave equation operator, $b_s$ is a source term representing an active seismic source at position $x_s$ and $m$ is the subsurface model that we ultimately reconstruct. The modelling operator $A(m)$ is not specified here: it can represent a simple scalar acoustic wave equation depending on a single velocity parameter, up to a fully anisotropic visco-elastic approximation, where the subsurface mechanical properties would be described by the density and a set of 21 independent coefficients in the Hooke's tensor, plus 21 quality factors, one per each of these coefficients. In this study we consider the case of the vertical transverse isotropic (VTI) visco-acoustic approximation, where the subsurface is described by the vertical pressure wave (P wave) velocity, the density, a quality

factor and two parameters representing the VTI anisotropy of the subsurface (Thomsen 1986; Yang et al. 2018).

The solution of the wave eq. (1) is denoted by $u_s[m](x, t)$ where the dependence of $u_s$ on the subsurface model $m$ through the wave equation is denoted using brackets, and its dependence to the spatial coordinates $x \in \Omega \subset \mathbb{R}^3$ and time $t \in [0, T]$ is denoted using parenthesis. Note that the nonlinearity of the FWI problem originates from the nonlinearity of the mapping between the model $m$ and the wavefield $u[m]$.

Once the wavefield $u_s[m]$ is computed, the calculated data is obtained as the extraction of the wavefield values at the receiver locations following an extraction operator $R$, such that

$$d_{cal,s}[m](\mathbf{x}_r, t) = R u_s[m](\mathbf{x}, t), \qquad (2)$$

which can be interpreted as a convolution of the wavefield with a delta Dirac function $\delta(\mathbf{x} - \mathbf{x}_r)$ centred at the receiver location $\mathbf{x}_r$

$$R u_s[m](\mathbf{x}, t) = \int_{\Omega} \delta(\mathbf{x} - \mathbf{x}_r) u_s[m](\mathbf{x}, t) d\mathbf{x}. \qquad (3)$$

The standard FWI formulation is the nonlinear least-squares minimization problem

$$\begin{aligned} \min_m f(m), \quad \text{where} \quad f(m) &= \frac{1}{2} \|d_{cal}[m] - d_{obs}\|^2 \\ &= \frac{1}{2} \sum_{s=1}^{N_s} \sum_{r=1}^{N_r} \int_0^T |d_{cal,s}[m](x_r, t) - d_{obs,s}(x_r, t)|^2 dt. \end{aligned} \qquad (4)$$

FWI can thus be interpreted as the problem of finding the subsurface parameters $m$ which provide the best fit between the calculated and the observed data in the least-squares sense. More general comparison functions can be used to mitigate the nonlinearity of the problem (Luo & Schuster 1991; van Leeuwen & Mulder 2010; Bozdağ et al. 2011; Górszczyk et al. 2021, see for instance). In this study we use the least-squares misfit function. The reason is the initial tomography P-wave velocity model we consider is sufficiently accurate to predict the data within half a period time-shift and a conventional multiscale strategy interpreting the data from low- to high-frequency bands is enough to mitigate the risk to converge towards a non-informative local minimum of the misfit function (an issue often referred to as cycle skipping).

### 2.2 Local optimization method and multiparameter strategy

The numerical solution of the FWI problem in eq. (4) is computed using local optimization strategies. The initial guess of the subsurface model is denoted by $m_0$, and a minimizing sequence of models $m_k$ is computed following the recurrence

$$m_{k+1} = m_k + \alpha_k \Delta m_k, \qquad (5)$$

where $\alpha_k$, referred to as the step size, is a positive scalar number computed via a line search procedure (Nocedal & Wright 2006; Bonnans et al. 2006), and $\Delta m_k$ is a descent direction. First-order methods, also known as steepest descent methods, use the opposite of the gradient of the misfit function $f(m_k)$ as a descent direction

$$\Delta m_k^{\text{steepest-descent}} = -\nabla f(m_k). \qquad (6)$$

These methods are easy to implement, however, they can exhibit very slow convergence. Second-order methods are based on the Newton descent direction

$$\Delta m_k^{\text{Newton}} = -H(m_k)^{-1} \nabla f(m_k), \qquad (7)$$
where $H(m_k)$ is the operator corresponding to the second-order partial derivatives of the misfit function. Newton-based descent methods exhibit a second-order convergence, and are therefore much faster than the steepest-descent method. The Hessian operator also carries the imprint of the parameter coupling in the frame of multiparameter inversion. Deconvolving the gradient by the inverse Hessian operator therefore mitigate the parameter coupling (Pratt et al. 1998; Operto et al. 2013; Métivier et al. 2015). The solution of the linear system in the Newton equation is however computationally too expensive for these methods to be used in practice for large realistic size problems (Métivier et al. 2017).

For this reason we rely on a compromise between the two approaches, where the inverse Hessian operator $(H(m_k))^{-1}$ is approximated at iteration $k$ from the $l$ most recent values of the gradient of the misfit function computed along the course of iterations (Nocedal 1980). The descent direction $\Delta m_k$ is thus given by

$$ \Delta m_k = -Q_k \nabla f(m_k), \quad (8) $$

where $Q_k$ is an approximation of the inverse Hessian

$$ Q_k \approx (\nabla^2 f(m_k))^{-1}. \quad (9) $$

This method represents a very good compromise in terms of efficiency, convergence speed and parameter decoupling. We complement it with a diagonal pre-conditioner which compensates for the energy loss of the wavefield with depth in the scenario of a surface acquisition. This pre-conditioner serves as a basis for the approximation of the inverse Hessian operator (Nocedal & Wright 2006; Métivier & Brossier 2016).

We also implement a unity-based normalization: each of the three parameters we invert for are rescaled between 0 and 1, to avoid any artificial numerical instabilities which would be due to the inversion of different parameters with different physical units and amplitude ranges. This rescaling is possible because the parameters we consider are bounded.

Finally, we choose to invert simultaneously for the vertical P-wave velocity, the density and the inverse of the quality factor. The velocity-density choice is motivated by the analysis of the parameter coupling based on radiation patterns (Operto et al. 2013). In these settings, the vertical velocity is sensitive to seismic events propagating in a broad range of angles, from small to large illumination angles. The density is sensitive only to events propagating with small illumination angles. This ensures a broad-band reconstruction of the velocity, a high-resolution for the density and helps to decouple between these two parameters in regions sampled by both diving and reflected waves. The choice to reconstruct the inverse of the quality factor is made to increase the sensitivity to the attenuation, which is otherwise weaker, and strongly nonlinear, in the range of variations we are considering. Using a logarithm of the quality factor could have also been considered as an alternative.

## 2.3 Gradient computation: the adjoint state method

The implementation of our local optimization scheme relies essentially on the computation of the gradient of the misfit function. A direct computation of the gradient requires the computation of the sensitivity matrix $\frac{\partial u}{\partial m}$ which is not affordable both in terms of computation and storage for the size of discrete problems associated with the FWI applications we are interested in (typically up to several hundred of millions of discrete unknown parameters). The key for such large scale optimization problem is a computing technique invented by the French mathematician Jacques-Louis Lions

in his seminal work on optimal control of partial differential equations (Lions 1968), named as the adjoint state method. A review of its application in geosciences is proposed in Plessix Plessix (2006), and we have recently proposed a rigorous mathematical derivation of these equations for visco-elastic FWI (Métivier & Brossier 2024). Following this technique, the gradient can be computed in our context as the sum for each source of the zero-lag time-correlation between the incident wavefield $u_s[m]$, solution of (1), and an adjoint field $\lambda_s[m]$ solution of a similar wave equation with a different source term, backpropagating in time. We denote it as

$$ A(m)^T \lambda_s = R^T (d_{cal,s}[m] - d_{obs,s}), \quad (10) $$

and

$$ \nabla f(m) = \sum_{s=1}^{N_s} \int_{0}^{T} u_s[m](x, t) \lambda_s[m](x, t) dt. \quad (11) $$

In eq. (10) we see that the adjoint field $\lambda_s$ is the solution of an adjoint wave equation involving the adjoint operator $A(m)^T$ and a source term which is equal to the difference between the calculated and the observed data for the seismic source $s$, localized at the receiver position through the adjoint operator $R^T$. The total computation cost for one gradient is equivalent to the computation of $2N_s$ wavefields and their time correlation, which makes possible the implementation of FWI to solve the large scale discrete problems associated with 3-D high-resolution imaging.

## 2.4 Regularization strategy: anisotropic diffusion smoothing

As an ill-posed nonlinear inverse problem, FWI requires regularization to converge towards meaningful and exploitable solutions. We rely on a recently developed anisotropic diffusion-based regularization strategy (Métivier & Brossier 2022; Métivier et al. 2024). This method has the advantage over more conventional strategies (such as Gaussian smoothing) to preserve the main structures in the gradient, which are automatically identified by the design of an anisotropic diffusion operator based on the structure tensor of the gradient (Weickert 1998; Lewis et al. 2014; Lee & Pyun 2018). Recent application on OBS data acquired offshore Australian on the Gorgon field illustrate the interest of this approach when strong contrasts are expected (Sambolian et al. 2023).

Instead of using directly the gradient $\nabla f(m)$ in the $l$-BFGS scheme presented in eqs (5) to (9), we first apply a structure-preserving smoothing filter which removes spurious oscillations from the gradient. These spurious oscillations originate typically from data noise, cycle-skipping of late arrivals or mismatch in Amplitude Versus Offset (AVO) variations. This smoothed version of the gradient is obtained by the solution of the following diffusion partial differential equations

$$ \left\{ \begin{array}{l} \frac{\partial \varphi}{\partial \tau} - \operatorname{div}(D \nabla \varphi) = 0 \\ \varphi(x, 0) = \nabla f(m), \end{array} \right. \quad (12) $$

with homogeneous Neumann boundary conditions. The smoothed gradient corresponds to the solution $\varphi$ after a certain diffusion time $\tau_{\text{diff}}$. We observe that the non-smoothed gradient $\nabla f(m)$ is the initial condition for the diffusion process (12). This strategy presents the advantage to smooth the gradient along structures present in the gradient but not across them. These structures are detected through the diffusion matrix $D$, built upon a singular value decomposition of the structure tensor associated with the gradient. The procedure originates from the work of Weickert (1998).
<figure><img src="image_3.png" /><figcaption>**Figure 1.** 2-D illustration of the anisotropic diffusion smoothing. (a) gradient before smoothing, (b) with isocontours superimposed, (c) and corresponding wavenumber spectrum, (d) gradient after conventional Gaussian smoothing, (e) with isocontours superimposed, (f) and corresponding wavenumber spectrum, (g) gradient after anisotropic diffusion smoothing, (h) with isocontours superimposed, (i) and corresponding wavenumber spectrum. One can appreciate how small scale oscillations are removed while the deep structure of the gradient is preserved using the anisotropic diffusion smoothing. Using a conventional Gaussian smoothing the small scale oscillations are efficiently removed but the deep structure is also destroyed. This loss of information is visible in the $(k_x, k_z)$ spectrum as the reduction of the spread of in the $k_z$ direction when using the Gaussian smoothing.</figcaption></figure>

A 2-D illustration is presented in Fig. 1. The initial gradient, without smoothing, is compared with a smoothed gradient obtained with a conventional Gaussian smoothing technique, and the smoothed gradient obtained with our anisotropic diffusion filter. In each case, we present a colour view of the (smoothed) gradient, the same view with iscontours overlaid and the corresponding wavenumber spectrum in $(k_x, k_z)$ domain. One can appreciate how small scale oscillations are removed while the deep structure of the gradient is preserved using the anisotropic diffusion smoothing. Using a conventional Gaussian smoothing the small scale oscillations are efficiently removed but the deep structure is also destroyed. The latter is visible in the $(k_x, k_z)$ spectrum where we can see how the spread of the energy is reduced along the $k_z$ dimension when using the Gaussian filter, which results from the loss of resolution in the vertical direction.

This structure preserving feature is especially key when applying FWI to broad-band data with the perspective to achieve high-resolution multiparameter reconstruction. This smoothing also automatically adapts to the scales of the structures to be preserved in the gradient, which is very interesting when the latter are diverse, which is often the case for practical applications. Finally, this type of regularization accelerates the convergence of the whole inversion process. Conventional smoothing such as the Gaussian smoothing removes partially the structural information from the gradient. This structural information is however associated with the fit to the data. This situation is typical of a slowly convergent process, where two sources of information contradict and lead the inversion towards orthogonal directions. A structure-preserving smoothing mitigates this shortcoming.

## 2.5 Fully scalable implementation

We summarize the FWI implementation as the flowchart presented in Fig. 2. Given a data set $d_{\text{obs}}$, and an initial model $m_0$, a minimizing sequence $m_k$ is built. When the misfit function has decreased below

a user-defined criterion $\epsilon$, the last computed model $m_k$ is output. The main building block of this flowchart is composed of five steps:

(i) Computing the misfit function gradient thanks to the adjoint state method (eqs 1, 10 and 11).

(ii) Smoothing the misfit function gradient using the diffusion-based approach (eq. 12).

(iii) Computing the descent direction using the *l*-BFGS approach (eq. 8).

(iv) Computing the line search parameter $\alpha_k$.

(v) Computing the new model $m_{k+1}$ (eq. 5).

For a FWI implementation to be fully scalable and able to exploit efficiently a HPC platform, each of these five steps needs to exploit computing resources in parallel, so as to avoid any sequential bottleneck.

The main step in terms of computation cost is the step 1, as it requires to solve one forward and one adjoint 3-D wave equation problem per source. The implementation we propose for step 1 relies on three nested levels of parallelism. The two first levels are based upon Message Passing Interface (MPI) parallelism. The third one exploits OpenMP multithreading parallelism. The MPI parallelism works both intra- and internodes, and does not rely on shared memory. The OpenMP parallelism exploits only intranode parallelism, and relies on the shared memory of the node. The combination of both parallelisms makes it possible to exploit at best the performances of HPC machines. The first MPI level is a parallelization over the source positions: each incident and adjoint wavefield can be computed independently for each source. The second MPI level involves a domain decomposition algorithm. The implementation we adopt relies on a finite-difference discretization of the VTI visco-ascoustic wave equation, and the representation of the subsurface model and discrete wavefields on a 3-D Cartesian mesh. We use a staggered-grid finite-difference discretization which is fourth-order in space and second-order in time, with an explicit time discretization (Yang et al. 2018).
<figure><img src="image_4.png" /><figcaption>**Figure 2.** General FWI workflow. Given input data $d_{\text{obs}}$, an initial model $m_0$ and a stopping criterion $\epsilon$, the model is iteratively updated to minimize the misfit function $f(m)$. Each step requires the computation of the misfit gradient $\nabla f(m)$, the smoothing of this gradient and the computation of the model update thanks to the l-BFGS formula. The model providing the best fit is output.</figcaption></figure>

The 3-D mesh is decomposed in subdomains of smaller size. The solution of the incident and adjoint wave propagation problems makes easy the use of this domain decomposition algorithm: at each time-step, the solution in each domain is computed by a set of independent computing units. To compute the solution at the boundaries of the subdomain, information from its neighbours is required. The computation of the wavefield inside the subdomain is then followed by a communication step with the neighbouring domains. An overlapping partition of the whole domain is performed, where the size of the overlap is adapted to the fourth-order scheme: the higher the order of the scheme, the larger the size of the overlap, as high-order scheme requires larger number of neighbouring points to evaluate the partial derivatives. The inner level of parallelization corresponds to multithreading through OpenMP directives to accelerates the 3-D loop required to compute the incident and adjoint field inside each domain at each time-step.

The subsequent steps (2 to 5) are less computationally intensive as they work on a single model or a single gradient. They need however to be treated with care so as to optimize the usage of HPC resources. In our implementation, we exploit the domain decomposition paradigm for all of these steps, complemented with OpenMP acceleration. For the smoothing operation (step 2), this is possible due to the fact that the anisotropic diffusion is formulated as the solution of a partial differential diffusion equation. We design an explicit finite-difference scheme to solve this diffusion equation which then exploits the domain decomposition algorithm initially set up for the solution of the wave propagation problem. Note that this is a specificity of the smoothing operator we consider: more conventional ones based on convolution with smooth functions such as Gaussian functions are less easy to make parallel using a domain decomposition algorithm. This is another interesting feature of the partial-differential-equations based smoothing we consider here.

For the linear algebra operations involved in the local optimization steps from 3 to 5 the parallelization is straightforward. In practice, all the scalar products and matrix-vector products required in

the model space are performed in parallel on each subdomain: the model or gradient in the whole domain are never assembled on a single computing unit, they remain split on different computing units thanks to the domain decomposition algorithm. This is facilitated by the reverse communication paradigm on which relies the SEISCOPE optimization toolbox which we use for local optimization (Métivier & Brossier 2016).

Finally, let us note that the time-correlation between the incident and adjoint fields required to compute the gradient (eq. 11) is complicated by the fact that the adjoint field backpropagates in time from a final state, while the incident field propagates from an initial state. For this reason it is difficult to access both fields at the same time-step, which is mandatory to evaluate the time-correlation between the two fields. Different strategies exist, either based on storage or recomputation techniques. In this study we store in memory the incident field whenever it is possible (a strategy referred to as the RAM technique in the next section), and switch to the Checkpointing Assisted Reverse Forward Simulation (CARFS) technique when the memory requirement becomes too important. The latter technique proposes a smart combination of memory storing and backward recomputation of the incident from its final state and decimated time history on the boundaries, to exploit at best the available memory and perform the least possible recomputation (Yang et al. 2016).

A scheme summarizing this three-levels parallelization is presented in Fig. 3. In this scheme, three nodes of 8 CPUs are available, for a total of 24 CPUs. Two source positions are considered. The whole domain is decomposed in 6 subdomains (we consider a 2-D illustration for the sake of simplicity, but this extends straightforwardly to the 3-D configuration used in this study). Two CPUs are assigned to each subdomain. This partition is consistent with the fact that the inner parallelism level is performed through OpenMP directive, which requires that each pair of CPU assigned to a subdomain is on the same node and share the same memory.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Typical computing resources allocation. Three computing nodes, each equipped with 8 CPUs, are available. The data set contains 2 seismic sources. For each source, the domain is partitioned into 6 subdomains. On each subdomain, 2 CPUs are used in parallel using OpenMP multithreading parallelization. The scheme presents how the available CPU are partitioned on each subdomain. As OpenMP uses shared memory, it is important that the processors assigned to each subdomain are on the same node. The domain-decomposition based MPI parallelism is not subjected to this constraint: for source 1, neighbouring domains $\Omega_4$ and $\Omega_6$ have been assigned CPUs from two different nodes (node 1 and node 2 respectively).</figcaption></figure>

# 3 APPLICATION TO A NORTH SEA OCEAN BOTTOM CABLE DATA SET

Now that we have reviewed our FWI methodology, we present its application to a dense Ocean Bottom Cable (OBC) data set recorded in the North Sea at Valhall, starting with a presentation of this data set.

## 3.1 Data presentation

The Valhall oil & gas field is an Upper Cretaceous chalk reservoir located in the Norwegian sector of the North Sea. The reservoir (consisting of the Tor and Hod formations) is located 2400 m below the seabed (Olofsson et al. 2003). Since the production started in Valhall in 1982, the originally highly porous chalk (more than 50 per cent) has continuously lost its porosity because of the effective stress increase induced by the hydrocarbon production. The pore collapse causes reservoir compaction (locally exceeding 10 m) and an associated seafloor subsidence reaching 6 m at the platform complex (Kristiansen & Flatebø 2010).

A permanent array of approximately 2400 multicomponent geophones and hydrophones has been deployed on the sea-floor in 2003, on a 300 m × 50 m nominal grid using OBC, for a total length of 120 km buried at 1.5 m to 2 m below the seabed. This permanent array has been used on 17 surveys, for monitoring purpose, from 2003 through 2014. It constitutes for the exploration community an exceptionally dense data set, considered as a reference for exploration and monitoring of a reservoir. For each survey, about 50 000 positions of airgun sources towed at 5 m depth in the water column were used, covering an area of 9 × 16 = 144 km² (Barkved et al. 2003). We use in this study the hydrophone component of the data acquired in 2011, the same which was used in previous FWI studies (Operto et al. 2015; Operto & Miniussi 2018; Kamath et al. 2021; Pladys et al. 2022), including the seminal paper (Sirgue et al. 2010). Fig. 4(a) displays the layout of the corresponding acquisition, where receivers along OBC appear in blue and the shot coverage appears in red. The density in terms of sources and receivers yields an outstanding quality data set with useful signal in a broad 2.5 – 100 Hz frequency band.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** (a) Layout of the Valhall acquisition overlapped on an horizontal *P*-wave velocity slice at 1 km depth obtained by FWI. Sources coverage (red cloud) and receivers (blue diamonds). One receiver position (A) is located with a black star. Cables A ($x$ = 2950 m), B ($x$ = 5530 m) and C ($x$ = 3080 m) are identified. Black dots denote the position of three available *P*-wave velocity sonic-logs and density measurement displayed in Figs 9, 10 and 11. (b) Cross-section through the 3-D volume of the common-receiver gather recorded at location marked by a black star in (a). The vertical dimension is the recording time while the horizontal directions correspond to the different shot positions. A spherical wave propagation pattern is visible on the horizontal layers. The solid white line (TR) delineates transmitted arrivals. The dashed white line (RF) delineates reflections. The dashed black line (ST) denotes Scholte waves (waves propagating at the water/subsurface interface). The dotted white line (MP) highlights multiple reflections propagating resulting from the interaction between the wavefield and the free surface at the water/air interface.</figcaption></figure>

An example of a 3-D common receiver gather is presented in Fig. 4(b). The seismograms are filtered in the 2.5–22.5 Hz frequency band we consider in this study. It exhibits complex transmitted (TR, solid white line) and reflected (RF, dashed white line) waveforms, indicating an extensive interaction with the subsurface heterogeneities. Moreover, we observe an imprint of the Scholte waves (ST, dashed black line) propagating at slower velocity along the interface between the water and the seabed, and

multiple reflections (MP, dotted white line) resulting from the interaction between the wavefield and the free surface at the water/air interface.

## 3.2 Practical application and workflow

The workflow we set up here relies on several years of investigation on this data set. It has been tackled first in the frequency domain
<table><caption>**Table 1.** Multiscale approach to FWI. Characteristics of the 8 frequency bands and corresponding meshes. Number of FWI performed per frequency band.</caption><thead><tr><th></th><th>Frequency band</th><th>Mesh size</th><th>Discrete points</th><th>Number of FWI run</th></tr></thead><tbody><tr><td>band 1</td><td>2.5–5 Hz</td><td>70 m</td><td>66 × 130 × 230</td><td>2</td></tr><tr><td>band 2</td><td>2.5–7 Hz</td><td>35 m</td><td>131 × 259 × 459</td><td>2</td></tr><tr><td>band 3</td><td>2.5–10 Hz</td><td>35 m</td><td>131 × 259 × 459</td><td>3</td></tr><tr><td>band 4</td><td>2.5–12.5 Hz</td><td>17.5 m</td><td>261 × 517 × 917</td><td>2</td></tr><tr><td>band 5</td><td>2.5–15 Hz</td><td>17.5 m</td><td>261 × 517 × 917</td><td>3</td></tr><tr><td>band 6</td><td>2.5–17.5 Hz</td><td>17.5 m</td><td>261 × 517 × 917</td><td>3</td></tr><tr><td>band 7</td><td>2.5–20 Hz</td><td>17.5 m</td><td>261 × 517 × 917</td><td>2</td></tr><tr><td>band 8</td><td>2.5–22.5 Hz</td><td>17.5 m</td><td>261 × 517 × 917</td><td>3</td></tr></tbody></table>

(Operto et al. 2015), then in the time-domain, focusing first on viscosity effects (Kamath et al. 2021), then on the mitigation of the inversion nonlinearity by the use of dedicated misfit functions (Plaudys et al. 2022), the exploitation of reflections (Provenzano et al. 2024), and the use of an anisotropic diffusion filter to enhance the resolution of the estimates in a multiparameter velocity-density inversion (Métivier et al. 2024). It relies on a frequency continuation strategy, the design of suitable initial models, the use of reciprocity between sources and hydrophones, a random shot subsampling algorithm and a source-wavelet estimation based on windowed arrivals. We present these ingredients in the following. Note that even if the specific settings of these ingredients depend on the Valhall data set, the proposed workflow in itself is generic and can be translated to other data sets.

### 3.2.1 Frequency continuation

The full frequency band we consider in this study, from 3 to 22.5 Hz, is split into 8 frequency bands. Each frequency band is obtained by applying a minimum-phase filter. For each frequency band we perform two or three runs of FWI depending on the evolution of the misfit function and the model. The output model from frequency band *i* serves as an input model for the frequency band *i* + 1. Such a strategy is documented as a multiscale approach in the literature (Bunks et al. 1995). The rationale is to invert first for the low frequency content of the data, to reconstruct first the large scale variation of the subsurface model parameters and progressively introduce higher frequency content in the inversion, to reconstruct smaller scale variations and increase the subsurface model resolution. This strategy acts as a regularization of the inverse problem, as the nonlinearity of the FWI problem is milder at lower frequency, and it makes it possible to start each FWI run from an accurate initial guess favoring the convergence.

Following this strategy, we consider three different spatial discretization steps: for the first frequency band, the step is set to 70 m (coarse grid). For the frequency band 2 and 3, the mesh is refined by a factor 2, the step being set to 35 m (medium grid). For frequency bands 4 to 8, an refinement by a factor 2 is used and the step size reaches 17.5 m (fine grid). The choice of the mesh is governed by the need to minimize the numerical dispersion associated with the finite-difference discretization of the wave equation, and the (crucial) importance of representing accurately the sea bottom interface. The water column being constant and equal to 70 m, we choose mesh sizes which are multiples of this value so as to have an integer number of discretization point in the water column. This information is summarized in Table 1.

### 3.2.2 Initial models

We parametrize the subsurface with the *P*-wave velocity $v_P$, the density $\rho$, a quality factor $Q_P$ and two Thomsen's parameters $\epsilon$ and $\delta$ accounting for the VTI anisotropy. We consider an initial *P*-wave velocity model $v_{P,0}$ obtained by reflection tomography. The initial density model $\rho_0$ is derived from $v_{P,0}$ through the following Gardner's law

$$ \rho = 309.6 V_p^{0.25}. \qquad (13) $$

The initial $Q_P$ model is taken constant, equal to 200 in the whole volume except in the water column where it is set to 1000 (no attenuation in the water layer). The initial $\epsilon$ and $\delta$ models are obtained by tomography and well log information respectively. The two latter are kept fixed while we simultaneously update $v_P$, $\rho$ and the inverse of the quality factor $Q_P^{-1}$ in the inversion.

### 3.2.3 Reciprocity and random shot subsampling

The acquisition is composed of more than 50 000 thousands shot positions and 2048 receivers on the seabed. The computation cost of FWI depending linearly on the number of seismic sources we consider, we resort to the reciprocity principle to transform our data into a set of 2048 seismic gathers which correspond to the data which would have been recorded by hydrophones at the shot positions, for explosive sources at the receiver positions. This simple data manipulation reduces by almost a factor 25 the total computation cost.

To decrease further the computational cost, it is usual to exploit a random source subsampling scheme (Warner et al. 2013). Rather than taking into account the full 2048 seismic gathers simultaneously during one FWI run, the data set is split into random subgroups of shot-gathers, with the constraint that the overlap between two subgroups is empty (a source belongs to only one subgroup). A sequence of FWI runs is then performed on each subgroup for few iterations, the subsurface model estimated for one subgroup being the initial guess for FWI on the next subgroup. In this study, we split the data, for each frequency band, into 16 subgroups each containing 128 shot-gathers. Three iterations of the *l*-BFGS algorithm are performed on each subgroup, for a total of 48 iterations per FWI run. Note that the subgroups are re-generated for each frequency band, to mitigate any imprint from this random selection in the final results. The random subsampling exploits the data redundancy, and finally decreases the total computation cost by a factor 16. The quality of the obtained model after 48 iterations is similar to what would have been obtained with the same number of iterations using the full data set (Kamath et al. 2021).
<table><caption>**Table 2.** Partitioning of the ADASTRA resources and computation time depending on the considered frequency band and the gradient computation technique. In the column titles, freq. band stands for 'frequency band', d.d. stands for 'domain decomposition', elaps. time for 'elapsed time', s.h. for 'scalar hours', the latter being the number of equivalent running time on a single CPU. The increase in computational cost from bands 4–5 to bands 6–8 is due to the gradient computation strategy, from storing the incident wavefield in memory (RAM) to recomputing it from a final condition and stored values at the boundaries (CARFS). The latter is less demanding in memory to the expense of additional computational efforts to recompute the incident field.</caption><thead><tr><th>freq. band</th><th>mesh size</th><th>gradient mode</th><th>d.d.</th><th># OpenMP</th><th># CPU</th><th>elaps. time</th><th>s.h.</th></tr></thead><tbody><tr><td>band 1</td><td>70 m</td><td>RAM</td><td>1 × 2 × 6</td><td>2</td><td>3072</td><td>6096 s</td><td>5200</td></tr><tr><td>band 2–3</td><td>35 m</td><td>RAM</td><td>1 × 2 × 8</td><td>3</td><td>6144</td><td>17047 s</td><td>29000</td></tr><tr><td>band 4–5</td><td>17.5 m</td><td>RAM</td><td>2 × 8 × 12</td><td>2</td><td>49152</td><td>22892 s</td><td>312000</td></tr><tr><td>band 6–8</td><td>17.5 m</td><td>CARFS</td><td>2 × 8 × 12</td><td>2</td><td>49152</td><td>82847 s</td><td>1131000</td></tr></tbody></table>

### 3.2.4 Source wavelet estimation

As a prior step before FWI on each frequency band, the time-signature of the seismic sources needs to be estimated. The sea bottom being flat and laterally homogeneous, we make the assumption that the virtual explosive sources (resorting on reciprocity) share the same time signature. We estimate this time signature using a linear deconvolution approach formulated in the frequency domain (Pratt 1999). This source signature remains constant in the FWI workflow and for subsequent FWI runs on the same frequency-band. It is re-estimated only when the frequency band changes. Note that the estimation of the source time signature is based on the early arrival of short offset data. This is achieved by a specific time-windowing strategy.

### 3.2.5 High-performance computing machine

We use in this study the Central Processing Unit (CPU) partition of the ADASTRA HPC platform, currently the largest in France, operated by CINES, and opened for academic research projects since November 2023. We had a special access to this machine earlier in June 2023, as our project was selected along few others to run on the machine. In total, the CPU partition consists of 536 nodes equipped with Bi-socket AMD Genoa CPU containing 2 × 96 computing cores at 2.4 GHz clock frequency. On each node, these computing units share 768 Gigabyte of DDR5 Random Access Memory (RAM). The peak performance of this CPU partition is evaluated at 3.9 Petaflops, that is, $3.9 \times 10^{15}$ floating-point operations per second.

We report in Table 2 the partitioning of the ADASTRA resources depending on the size of the mesh, as well as the elapsed time and scalar time for one FWI run. We were able to exploit simultaneously up to 256 nodes (half the Adastra CPU partition), making a simultaneous use of 49 152 computing units to invert our seismic data. Scalar time is a measure of computational cost equal to the elapsed time multiplied by the number of CPU used during one run. On the fine mesh, starting from frequency band 6, the computation of the gradient using the RAM technique is not possible anymore, because the correlation of the incident and adjoint wavefield requires to be done at a finer time sampling according to Nyquist criterion. There is no sufficient memory available to store the wavefield in RAM, even relying on the domain decomposition strategy. We thus switch to the CARFS technique, which requires less memory, but, as can be seen from the Table 2, is three times more expensive in terms of computational cost (Yang et al. 2016).

## 3.3 Results: high-resolution crust characterization

Among the three parameters we reconstruct, the P-wave velocity has a first-order impact on the kinematics and the dynamics

of the waves. It is therefore considered as the primary parameter. In Figs 5 and 6 we present the initial (tomography) and final 3-D FWI P-wave velocity models (after FWI up to 22.5 Hz). One can observe the tremendous resolution improvement in the final FWI model (Fig. 5b) compared with the initial model (Fig. 5a). Features to the order of few tens of metres become accessible. Vertical sections in Fig. 5(b) are crossing through the low-velocity region in the central part of the 3-D volume we further refer to as 'gas-cloud' region (GC in Fig. 5b). In this zone the seismic imaging is severely degraded, presumably due to sediments charged with hydrocarbons to varying extents, resulting strong lateral and vertical variations in P-wave velocities, strong attenuation of the high frequency part of the data and associated phase-velocity dispersion.

The resolution increase is not only pronounced within the volume comprising this region. On the sides of this anomaly, the final FWI model exhibits a clear layering structure, which conforms with the available stratigraphic information (see the stratigraphic table in Fig. 5b). In the shallow seabed, the youngest Quaternary layers (made of clay, silts and sands) exhibit P-wave velocity values up to ~1900 m s⁻¹. Below, the Neogene/Paleogene package consists of claystone, siltstone and shale in its upper part and claystone with limestone stingers and limestone in its lower part. The Late Miocene and Pliocene sections is mainly claystone. They act as a caprock to the underlying hydrocarbon-charged sequence of early to late Miocene low-velocity siltstone and silty sandstone at the base, grading upwards into claystone (Munns 1985; Kristiansen et al. 2000).

Disturbed layering is also clearly reconstructed in the deeper part of the model, that is, below the boundary around ~2.6 km. Around this depth, the rapid increase of velocity is consistent with a change of lithology marking the Upper Cretaceous chalk with Tor and Hod formations. A strong deeper reflector, below the reservoir (R in Fig. 5b) appears clearly, with a remarkable continuity. Moreover, continuous structures, probably due to dipping strata, also appear in the horizontal plane at 4.5 km depth, which are absent from the initial tomography model (Fig. 5a).

In Fig. 6 two depth-slice sections at 200 and 1000 m depth are presented. At 200 m the top of the low-velocity anomaly emerges together with a detailed network of paleo-channels at scales ranging from few tens to few hundreds metres. At about 1 km depth, we can see the lateral extent of the low-velocity anomaly related to presence of hydrocarbons. Primary and secondary faults and fractures show low-velocities indicating possible migration pathways for these hydrocarbon-charge sediments.

An impressive amount of geological information can be already inferred from the vertical velocity model in Fig. 5(b). Moreover, thanks to our multiparameter FWI approach, we also reconstruct 3-D density models. In Fig. 7, we present an intermediate density
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Cross-sections through the 3-D vertical velocity (a) initial and (b) final FWI models. The stratigraphic table shows the changes in lithology derived from the drilling in the central Graben structure (Munns 1985; Kristiansen et al. 2000). The reservoir R is located below the gas-cloud GC in the oil-bearing upper-cretaceous Tor and Hod formations T&H.</figcaption></figure>

model obtained at 10 Hz and the final density model obtained at 22.5 Hz where the resolution improvement at 22.5 Hz is already visible.

We also build 3-D reflectivity models, taking the vertical derivative of the 3-D impedance model, that is, the product of the *P*-wave velocity and the density. These reflectivity models at 10 and 22.5 Hz are presented in Fig. 8. The resolution enhancement is remarkable in Fig. 8(b), which displays a 3-D view of the reflectivity in the GC region. The image obtained at 22.5 Hz shows a clear delineation of sedimentary layers. This is true both for the shallower sedimentary part, near the seafloor, and the deeper part of the model, where the previously identified continuous structures, associated with dipping strata, invisible at 10 Hz, appear in the horizontal plane at 4.5 km depth in the 22.5 Hz image.

To control the quality of the *P*-wave velocity and density reconstruction, we have access to three wells inside which sonic logging (Fig. 9) and density measurement (Fig. 10) have been performed. Their locations are indicated in Fig. 4. The comparison between the reconstructed models and the well-log quantities is presented in Figs 9 and 10. The reconstructed velocity fits correctly the sonic

logging down to 3 km depth, with a clear improvement from 10 to 22.5 Hz. Interestingly, the reconstructed density also matches the measurement, down to 2 km depth, with an improvement from 10 to 22.5 Hz, especially from 1.5 to 2 km depth. Deeper, the density variations appear overestimated, especially on log 3, while the measured impedance is still fitted, as expected from radiation pattern interpretation (Operto et al. 2013). In regions where the illumination is rich enough (reflections and refraction), our algorithm appears as capable of exploiting the increase in frequency band to decouple efficiently density and velocity parameters. Deeper, the decoupling is less effective, while the overall reconstructed impedance is still quantitatively estimated, as can be seen on the impedance measurement logs presented in Fig. 11.

We reconstruct also the quality factor model, which characterizes the wavefield intrinsic attenuation, namely the energy loss per wavelength. In the overburden, the attenuation effect of the fluid-filled sediment and shale is well acknowledged (Munns 1985), although there is no possibility to directly measure this parameter from wells. It is also difficult to reconstruct the spatial variations of this parameter, especially on narrow frequency band data (Kamath et al. 2021).
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Depth-slice sections extracted from the 3-D vertical velocity model in Fig. 5(b) at (a) 200 m and (b) 1000 m.</figcaption></figure>

In this study, as we increase the frequency content of the inverted data, we increase the sensitivity to this parameter. To assess the reliability of our 3-D reconstruction, we compare it with the 3-D structure of the reconstructed reflectivity model (Fig. 12) in the gas cloud region. This comparison (Fig. 12b), shows that the reflectivity within the GC is decreasing rapidly with depth, as soon as the value of the quality factor starts to increase—indicating a higher attenuation. The 3-D view of the fluid-filled region in Fig. 12(b) also illustrates how the different layers are interconnected through

a network of faults and fractures up to the top surface of this region at 1 km depth.

An additional quality control in the reconstruction of the quality factor $Q_P$ is to compare modelled data in the final P-wave velocity and density models with a homogeneous $Q_P$ model set at $Q_P = 2000$ and with the reconstructed heterogeneous $Q_P$ model. This comparison is presented in Figs 13 and 14. In Fig. 13, we focus on residual shot gathers. The effect of the heterogeneous $Q_P$ model remains subtle, however it is clear that it helps reducing the
<figure><img src="image_9.png" /><figcaption>Figure 7. Cross-sections through the 3-D FWI density models at (a) 10 Hz and (b) 22.5 Hz.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Cross-sections through the 3-D FWI pseudo-reflectivity, that is, vertical derivative of impedance models at (a) 10 Hz and (b) 22.5 Hz.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Sonic logging compared to the initial model (dotted red line), the reconstructed model at 10 Hz (solid purple line), and the reconstructed model at 22.5 Hz (solid orange line) along the three wells whose location are indicated in Fig. 4.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Density measurement compared to the initial model (dotted red line), the reconstructed model at 10 Hz (solid purple line) and the reconstructed model at 22.5 Hz (solid orange line) along the three wells whose location are indicated in Fig. 4.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Impedance measurement compared to the initial model (dotted red line), the reconstructed model at 10 Hz (solid purple line) and the reconstructed model at 22.5 Hz (solid orange line) along the three wells whose location are indicated in Fig. 4.</figcaption></figure>

misfit on the transmitted arrivals (diving waves) pointed by the red arrows A1 to A4. It is anticipated that $Q_p$ would have a stronger effect on transmitted waves as it has a cumulative effect for phases propagating in this regime.

To further analyse the effect of the heterogeneous $Q_p$ on the data fit, we present in Fig. 14 traces extracted at the receiver positions corresponding to the arrows A1 to A4 in Fig. 13. For each receiver, we plot the observed trace (in black), and the corresponding synthetic trace using the heterogeneous $Q_p$ model (in blue) and the constant $Q_p = 200$ model (in red). We see that the heterogeneous $Q_p$ model helps in better fitting the amplitude of the transmitted data. Its effect in terms of phase dispersion remains almost invisible, despite the relatively broad frequency band considered here.

## 4 DISCUSSION

The results achieved in this study, in terms of multiparameter reconstruction of the subsurface mechanical properties from dense seismic data, open interesting perspectives for the quantification and monitoring of subsurface properties at specific targets such as carbon capture storage reservoir, going potentially toward data-driven prediction of loss of integrity, migration pathways and fluid escape. To achieve such characterization requires the ability to correctly predict and integrate in the inversion the subtle (nonlinear)

effects the subsurface heterogeneities have on the recorded seismogram. This can only be done at the expense of modelling seismic waves in a broad frequency band, which requires to use efficiently high performance computing resources thanks to dedicated, fully parallel algorithms, such as the one we present in this study.

Nonetheless, it is worth pointing out that so far we are only partially exploiting the wealth of information contained in the data. The frequency content could still be increased, at the expense of a higher computational cost. This is a viable goal. However, at higher frequency, it is expected that elastic effects become predominant, which in turn will require to adopt an elastic wave propagation modelling rather than the acoustic one we use in this study. Secondly, only the hydrophone data is fully exploited in this study. We do not use the ground motion recorded by the three-components geophones. These recordings bear a strong imprint of the visco-elasticity of the subsurface: the Scholte waves and converted shear-waves dominate the data, which cannot be predicted in the acoustic approximation.

We use the visco-acoustic approximation for three reasons. First, in media where wave velocity contrasts are relatively mild, as is the case here, this approximation remains relevant in the considered frequency band. This would not be the case in the presence of salt structures for instance (Rivera *et al.* 2019). Secondly, the nonlinearity between the wavefield and the subsurface parameters increases
<figure><img src="image_14.png" /><figcaption>**Figure 12.** (a) Cross-section through the 3-D pseudo-reflectivity model with vanishing reflectivity in the gas-cloud area. (b) Same as (a) but with superimposed model of the quality factor showing high attenuation zone within the area of low reflectivity in (a).</figcaption></figure>

significantly in the visco-elastic approximation, making the FWI problem even more challenging to solve. Thirdly, the computation cost of visco-elastic FWI is at least two orders of magnitudes higher than its visco-acoustic counterpart.

On the other hand, a fully visco-elastic characterization of the crust is desirable, giving the possibility to improve further the detection of fluids and discriminate between gas and liquids through P-wave velocity to shear wave (S-wave) velocity ratio. Please note
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Comparison of residuals for one common receiver gather along cable A (see Fig. 4). (a) residuals computed in the final *P*-wave velocity and density models and the final reconstructed heterogeneous *Q*<sub>P</sub> model. (b) same but using a homogeneous *Q*<sub>P</sub> model taken equal to 200 in the whole model. Red arrows point diving waves where the residuals are notably reduced thanks to the heterogeneous *Q*<sub>P</sub> model.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Traces comparison for the four receivers corresponding to arrows A1, A2, A3, A4 in the residual shot gathers presented in Fig. 13. (a) Comparison at receiver A1, (b) comparison at receiver A2, (c) comparison at receiver A3, (d) comparison at receiver A4. For each receiver, the observed trace, in black, is compared to the synthetic trace after inversion, using the heterogeneous reconstructed $Q_p$ model (in blue) or the homogeneous $Q_p$ model (in red).</figcaption></figure>

that current elastic FWI implementations promoted in exploration seismic mostly rely on an automatic conversion from P-wave velocity to S-wave velocity model or focus on the interpretation of P wave and converted P-to-S waves rather than trying to invert for the full elastic waveform. We believe however that only elastic FWI of the whole data will make possible to achieve the fine-scale quantitative characterization of the crust we are targeting.

This will be possible only by fully exploiting the multicomponent geophone data, that is horizontal, transverse and vertical geophone data, while only hydrophone or vertical geophone data are usually considered. There is however a wealth of information in the multicomponent data, notably the polarization state of each of the recorded seismic phases. Recent work show the interest to accurately take into account this information using FWI with a time-dependent polarization attribute (Sambolian et al. 2022). The latter study puts forward the idea that schemes dedicated to a proper exploitation of the multicomponent data should be developed further in the future.

Towards elastic FWI, we also have implemented parallel numerical tools for visco-elastic full waveform modelling and inversion in our code SEM46 (Trinh et al. 2019; Cao et al. 2022b). The progress to make now is mostly in mitigating the nonlinearity of visco-elastic FWI and the development of a suitable HPC strategy

able to exploit the newest generation of supercomputers, reaching exascale performances (10¹⁸ FLOPS). Regarding the nonlinearity of the FWI problem, it is important to develop adequate hierarchical schemes integrating progressively the multicomponent data, combined with suitable parametrization. Misfit function design will be also key to mitigate nonlinearities such as cycle skipping (Luo & Schuster 1991; van Leeuwen & Mulder 2010; Bozdağ et al. 2011; Górszczyk et al. 2021). Finally, regularization techniques which could exploit homogenization tools (Cupillard & Capdeville 2018) appear promising for elastic FWI. This theory assesses that sub-wavelength heterogeneities are sensed by the wavefield as a locally smooth medium, potentially anisotropic and provides numerical tools to compute equivalent smooth anisotropic media from fine scale media, in a given frequency band. This in turn has the consequence that FWI should be able to reconstruct at best an homogenized version of the exact subsurface mechanical properties in the frequency band of the inverted data. This has been highlighted on schematic 2-D examples (Capdeville & Métivier 2018) and we are now working in the generalization of these results in more realistic settings.

Regarding the HPC strategy, versatile implementations able to exploit simultaneously CPU and graphic cards accelerated computing units (GPU) are required (Cao et al. 2022a). We have also
performed first steps in this direction, by designing an orchestrator overarching the general FWI algorithm to distribute dynamically and automatically the computation on a given cluster. This makes it possible to run FWI in a much more flexible way, adapting to the available computational resources at a given time, and making the FWI code much more resilient to localized failures, with the capability to automatically detect failures and relaunch affected parts of the computation without losing the results of the whole computation (Schuster et al. 2024).

## 5 CONCLUSION

The fine-scale characterization of the first kilometres of the crust is a particularly challenging problem which will be essential in the implementation of green energy transition. This task requires going beyond conventional seismic imaging, towards 3-D multiparametric quantitative reconstructions of the subsurface mechanical parameters, bridging the gap between seismic imaging and reservoir/site characterization.

In this study we show how such characterization can be achieved by inverting high quality seismic data with a suitable FWI algorithm. We reconstruct a 3-D multiparameter model of P-wave velocity, density and quality factor, revealing 3-D fine-scale geological structures such as paleo-channels nearby the seafloor and complex networks of fluid-filled faults and fractures. Comparison with direct measurements in wells show that not only the reconstructed P-wave velocity model but also the density model can be quantitatively interpreted down to approximately 2 km depth, showing the possibility for FWI to reliably estimate simultaneously multiple subsurface mechanical parameters provided the data is rich enough.

## ACKNOWLEDGMENTS

The authors would like to express their gratitude to Espen Raknes, Ross Milne and Einar Jarle Kjos from Aker BP for sharing the data and their help and advice in assembling this study. This study was partially funded by the SEISCOPE consortium (http://seiscope2.osug.fr), sponsored by Aker BP, DUG, EXXONMOBIL, GEOLINKS, JGI, PETROBRAS, SHEARWATER, 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.

## DATA AVAILABILITY

The data used in this study is an industrial data set acquired and made available to us by the AkerBP company. Besides technical difficulties to make it available openly (the data set represents several tens of terabytes) we do not own this data and cannot make it accessible openly. It would however be available upon request to the Aker BP company.

The full waveform inversion methodology used in this study is implemented in the TOYxDAC_TIME FORTRAN90 code we have developed within the framework of the SEISCOPE project, funded by industrial partners. Under the contracts gathering our partners in the SEISCOPE project, it has been agreed that this code is not open-source. However the TOYxDAC_TIME code can be (and is

currently) shared with academic researchers. The interested readers are encouraged to contact the corresponding author.

## REFERENCES

Aubert, C. et al., 2025. The MACIV multiscale seismic experiments in the French Massif Central (2023-2027): deployment, data quality and availability, *Ann. Geophys.*, **68**, doi:10.4401/ag-9264

Barkved, O., Bærheim, A., Howe, D., Kommedal, J. & Nicol, G., 2003. Life of Field Seismic Implementation - Another “First at Valhall”, in *65th EAGE Workshop, Stavanger*.

Bonnans, J.F., Gilbert, J.C., Lemaréchal, C. & Sagastizábal, C.A., 2006. *Numerical Optimization, Theoretical and Practical Aspects*, Springer.

Bozdağ, E., Trampert, J. & Tromp, J., 2011. Misfit functions for full waveform inversion based on instantaneous phase and envelope measurements, *Geophys. J. Int.*, **185**(2), 845–870.

Brenguier, F. et al., 2023. A long-term, 300 nodes seismic array to monitor seismic velocities across the San Jacinto fault (South. Cal.), in *AGU Fall Meeting Abstracts*, Vol. **2023**, pp. S23D–0396.

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

Cao, J., Brossier, R., Cabrera, E., de la Puente, J., Métivier, L. & Tarayoun, A., 2022a. GPU accelerated computing towards a fast and scalable seismic wave modelling in SEISCOPE SEM46 code, in *Sixth EAGE Workshop on High Performance Computing*, European Association of Geoscientists & Engineers.

Cao, J., Brossier, R., Górszczyk, A., Métivier, L. & Virieux, J., 202b. 3D multi-parameter full-waveform inversion for ocean-bottom seismic data using an efficient fluid-solid coupled spectral-element solver, *Geophys. J. Int.*, **229**(1), 671–703.

Capdeville, Y. & Métivier, L., 2018. Elastic full waveform inversion based on the homogenization method: theoretical framework and 2-D numerical illustrations, *Geophys. J. Int.*, **213**(2), 1093–1112.

Claerbout, J., 1985. *Imaging the Earth's interior*, Blackwell Scientific Publication.

Cupillard, P. & Capdeville, Y., 2018. Non-periodic homogenization of 3-D elastic media for the seismic wave equation, *Geophys. J. Int.*, **213**(2), 983–1001.

Devaney, A., 1984. Geophysical diffraction tomography, *IEEE Trans. Geosci. Remote Sens.*, **GE-22**(1), 3–13.

Espin, I., Salaun, N., Jiang, H. & Reinier, M., 2023. From FWI to ultra-high-resolution imaging, *Leading Edge*, **42**(1), 16–23.

Gimbert, F. et al., 2021. A multi-physics experiment with a temporary dense seismic array on the Argentière glacier, french alps: The RESOLVE project, *Seismol. Res. Lett.*, **92**(2A), 1185–1201.

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.

Górszczyk, A., Almeida, R., Brossier, R., Métivier, L. & Operto, S., 2025. Constraining the effect of bathymetric high subduction in an accretionary wedge: Evidence from the Tokai area of the Nankai trough, Japan, *Geochem. Geophys. Geosyst.*, **26**(3), e2024GC011754.

Hetényi, G., Molinari, I., Clinton, J. &, 2018. The AlpArray Seismic Network: a large-scale european experiment to image the Alpine Orogen, *Surv. Geophys.*, **39**.

Hirst, S., Taylor, S., Jones, O., Dancer, K. & Dheile, P., 2024. High-frequency MP-FWI imaging and AVA for shallow targets: Barents sea, in *84th EAGE Annual Conference & Exhibition*, Vol. **2024**, pp. 1–5, European Association of Geoscientists & Engineers.

Huang, R., Zhang, Z., Wu, Z., Wei, Z., Mei, J. & Wang, P., 2021. Full-waveform inversion for full-wavefield imaging: decades in the making, *Leading Edge*, **40**(5), 324–334.

Kalinicheva, T., Warner, M. & Mancini, F., 2021. Full-bandwidth FWI, in *84th Annual EAGE Meeting (Madrid)*, European Association of Geoscientists & Engineers.
Kamath, N., Brossier, R., Métivier, L., Pladys, A. & Yang, P., 2021. Multi-parameter full-waveform inversion of 3D ocean-bottom cable data from the Valhall field, *Geophysics*, **86**(1), B15–B35.

Kristiansen, T.G. & Flatebø, R.E., 2010. Sixty days ahead of schedule: reducing drilling risk at Valhall using computational geomechanics, *SPE Drilling Complet.*, **25**(04), 544–554.

Kristiansen, T.G., Barkved, O. & Pattillo, P.D., 2000. Use of passive seismic monitoring in well and casing design in the compacting and subsiding Valhall field, North Sea, SPE Europec, in *EAGE Conference and Exhibition*, pp. SPE-65134-MS, EAGE.

Lailly, P., 1983. The seismic inverse problem as a sequence of before stack migrations, in *Conference on Inverse Scattering, Theory and application*, Society for Industrial and Applied Mathematics, Philadelphia, pp. 206–220.

Lee, D. & Pyun, S., 2018. Adaptive preconditioning of full-waveform inversion based on structure-oriented smoothing filter, in *SEG Technical Program Expanded Abstracts 2018*, pp. 1048–1052.

Lee, H. et al., 2023. *Climate change 2023: Synthesis report. contribution of working groups I, II and III to the sixth assessment report of the Intergovernmental Panel on Climate Change*, IPCC, Geneva, pp. 35115.

Lehujeur, M., Vergne, J., Schmittbuhl, J., Zigone, D., Le Chenadec, A. & Team E., 2018. Reservoir imaging using ambient noise correlation from a dense seismic network, *J. geophys. Res.: Solid Earth*, **123**(8), 6671–6686.

Lewis, W., Amazonas, D., Vigh, D. & Coates, R., 2014. Geologically constrained full-waveform inversion unsing an anisotropic diffusion based regularization scheme: application to a 3D offshore Brazil data set, in *SEG Expanded Abstracts*, pp. 1083–1088, Society of Exploration Geophysicists.

Lions, J.L., 1968. *Contrôle optimal de systèmes gouvernés par des équations aux dérivées partielles*, Dunod.

Luo, Y. & Schuster, G.T., 1991. Wave-equation traveltime inversion, *Geophysics*, **56**(5), 645–653.

Malvoisin, B. & Brunet, F., 2023. Barren ground depressions, natural H2 and orogenic gold deposits: Spatial link and geochemical model, *Sci. Total Environ.*, **856**, 158969.

Métivier, L. & Brossier, R., 2016. The SEISCOPE optimization toolbox: a large-scale nonlinear optimization library based on reverse communication, *Geophysics*, **81**(2), F11–F25.

Métivier, L. & Brossier, R., 2022. On the use of nonlinear anisotropic diffusion filters for seismic imaging using the full waveform, *Inverse Probl.*, **38**(11), 115001.

Métivier, L., Brossier, R., Operto, S. & Virieux, J., 2015. Acoustic multi-parameter FWI for the reconstruction of *P*-wave velocity, density and attenuation: preconditioned truncated Newton approach, in *SEG Technical Program Expanded Abstracts*, pp. 1198–1203, SEG.

Métivier, L., Brossier, R., Operto, S. &, 2017. Full waveform inversion and the truncated Newton method, *SIAM Rev.*, **59**(1), 153–195.

Métivier, L., Brossier, R., Hoffmann, A., Mirebeau, J.-M., Provenzano, G., Tarayoun, A. & Yong, P., 2024. Coherence-enhancing anisotropic diffusion filter for 3D high resolution reconstruction of *P*-wave velocity and density using full waveform inversion: application to a North Sea Ocean Bottom Cable data set, *Geophysics*, **89**(1), R33–R58.

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

Munns, J.W., 1985. The Valhall field: a geological overview, *Mar. Petrol. Geol.*, **2**, 23–43.

Métivier, L. & Brossier, R., 2024. On the adjoint state method for the gradient computation in full waveform inversion: a complete mathematical derivation for the (visco-)elastodynamics approximation, *Geophys. J. Int.*, **240**(2), 942–966.

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

Nocedal, J. & Wright, S.J., 2006. *Numerical Optimization*, 2nd edn, Springer.

Nolet, G., 1987. *Seismic Tomography with Applications in Global Seismology and Exploration Geophysics*, D. Reidel publishing Company.

Olofsson, B., Probert, T., Kommedal, J. & Barkved, O., 2003. Azimuthal anisotropy from the valhall 4C 3D survey, *Leading Edge*, **22**, 1228–1235.

Operto, S. & Miniussi, A., 2018. On the role of density and attenuation in 3D multi-parameter visco-acoustic VTI frequency-domain FWI: an OBC case study from the North Sea, *Geophys. J. Int.*, **213**, 2037–2059.

Operto, S., Brossier, R., Gholami, Y., Métivier, L., Prieux, V., Ribodetti, A. & Virieux, J., 2013. A guided tour of multiparameter full waveform inversion for multicomponent data: from theory to practice, *Leading Edge*, **32**(9), 1040–1054.

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.

Pladys, A., Brossier, R., Kamath, N. & Métivier, L., 2022. Robust FWI with graph space optimal transport: application to 3D OBC Valhall data, *Geophysics*, **87**(3), 1–76.

Plessix, R.E., 2006. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications, *Geophys. J. Int.*, **167**(2), 495–503.

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.

Pratt, R.G., Shin, C. & Hicks, G.J., 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion, *Geophys. J. Int.*, **133**, 341–362.

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

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

Rivera, C.A., Trinh, P.-T., Bergounioux, E. & Duquet, B., 2019. Elastic multiparameter FWI in sharp contrast medium, in *SEG Technical Program Expanded Abstracts 2019*, pp. 1410–1414, Society of Exploration Geophysicists.

Salaun, N., Toubiana, H., Mitschler, J.-B., Gigou, G., Carriere, X., Maurer, V. & Richard, A., 2020. High-resolution 3D seismic imaging and refined velocity model building improve the image of a deep geothermal reservoir in the Upper Rhine Graben, *Leading Edge*, **39**(12), 857–863.

Sambolian, S., Brossier, R. & Métivier, L., 2022. Exploiting the richness of multi-component data: a time-dependent polarization-based fwi approach, in *83rd Annual EAGE Meeting (Madrid)*, European Association of Geoscientists & Engineers.

Sambolian, S., Brossier, R. & Métivier, L., 2023. Time-domain FWI of full azimuth OBN data from offshore North West Australia, in *84th Annual EAGE Meeting (Vienna)*, Vol. **2023**, European Association of Geoscientists & Engineers.

Schuster, D., Brossier, R., Métivier, L. & Tarayoun, A., 2024. On the efficient use of HPC resources: disaggregation of the full waveform inversion algorithm, in *85th EAGE Annual Conference & Exhibition*, Vol. **2024**, pp. 1–5, European Association of Geoscientists & Engineers.

Sirgue, L., Barkved, O.I., Dellinger, J., Etgen, J., Albertin, U. & Kommedal, J.H., 2010. Full waveform inversion: the next leap forward in imaging at Valhall, *First Break*, **28**, 65–70.

Tarantola, A., 1984. Linearized inversion of seismic reflection data, *Geophys. Prospect.*, **32**, 998–1015.

Thomsen, L.A., 1986. Weak elastic anisotropy, *Geophysics*, **51**, 1954–1966.

Trinh, P.T., Brossier, R., Métivier, L., Tavard, L. & Virieux, J., 2019. Efficient 3D time-domain elastic and viscoelastic Full Waveform Inversion using a spectral-element method on flexible Cartesian-based mesh, *Geophysics*, **84**(1), R75–R97.

van Leeuwen, T. & Mulder, W.A., 2010. A correlation-based misfit criterion for wave-equation traveltime tomography, *Geophys. J. Int.*, **182**(3), 1383–1394.

Virieux, J. & Operto, S., 2009. An overview of full waveform inversion in exploration geophysics, *Geophysics*, **74**(6), WCC1–WCC26.
Virieux, J., Asnaashari, A., Brossier, R., Métivier, L., Ribodetti, A. & Zhou, W., 2017. An introduction to Full Waveform Inversion, in *Encyclopedia of Exploration Geophysics*, pp. R1–1-R1-40, eds Grechka, V. & Wapenaar, K., Society of Exploration Geophysics.

Warner, M. et al., 2013. Anisotropic 3D full-waveform inversion, *Geophysics*, **78**(2), R59–R80.

Weickert, J., 1998. *Anisotropic Diffusion in Image Processing*, Treubner Verlag.

Wu, R.S. & Toksöz, M.N., 1987. Diffraction tomography and multisource holography applied to seismic imaging, *Geophysics*, **52**, 11–25.

Wu, Y., Chen, D., Lin, T., Hsieh, C., Chin, T., Chang, W., Li, W. & Ker, S., 2013. A high-density seismic network for earthquake early warning in

taiwan based on low cost sensors, *Seismol. Res. Lett.*, **84**(6), 1048–1054.

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.

Yang, P., Brossier, R., Métivier, L., Virieux, J. & Zhou, W., 2018. A time-domain preconditioned truncated newton approach to multiparameter visco-acoustic full waveform inversion, *SIAM J. Sci. Comput.*, **40**(4), B1101–B1130.

Zhang, L., Sengupta, S., Parekh, C. & Eliott-Lockhart, R., 2022. High-frequency FWI velocity model building for hydrocarbon delineation, in *Second International Meeting for Applied Geoscience & Energy*, pp. 817–821, Society of Exploration Geophysicists.