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

HAL
open science

# On the role of density and attenuation in three-dimensional multiparameter viscoacoustic VTI frequency-domain FWI: an OBC case study from the North Sea

S. Operto, A. Miniussi

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

## To cite this version:

S. Operto, A. Miniussi. On the role of density and attenuation in three-dimensional multiparameter viscoacoustic VTI frequency-domain FWI: an OBC case study from the North Sea. Geophysical Journal International, 2018, 213 (3), pp.2037 - 2059. 10.1093/gji/ggy103. hal-01852750

HAL Id: hal-01852750

https://hal.science/hal-01852750

Submitted on 13 Oct 2021

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_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_4.png" /></figure>

# On the role of density and attenuation in three-dimensional multiparameter viscoacoustic VTI frequency-domain FWI: an OBC case study from the North Sea

S. Operto<fnref n="1" /> and A. Miniussi<fnref n="2" />

Received 2018 February 12; in original form 2017 October 30

## SUMMARY

3-D frequency-domain full waveform inversion (FWI) is applied on North Sea wide-azimuth ocean-bottom cable data at low frequencies (≤10 Hz) to jointly update vertical wave speed, density and quality factor Q in the viscoacoustic VTI approximation. We assess whether density and Q should be viewed as proxy to absorb artefacts resulting from approximate wave physics or are valuable for interpretation in the presence of soft sediments and gas cloud. FWI is performed in the frequency domain to account for attenuation easily. Multiparameter frequency-domain FWI is efficiently performed with a few discrete frequencies following a multiscale frequency continuation. However, grouping a few frequencies during each multiscale step is necessary to mitigate acquisition footprint and match dispersive shallow guided waves. Q and density absorb a significant part of the acquisition footprint hence cleaning the velocity model from this pollution. Low Q perturbations correlate with low-velocity zones associated with soft sediments and gas cloud. However, the amplitudes of the Q perturbations show significant variations when the inversion tuning is modified. This dispersion in the Q reconstructions is however not passed on the velocity parameter suggesting that cross-talks between first-order kinematic and second-order dynamic parameters are limited. The density model shows a good match with a well log at shallow depths. Moreover, the impedance built a posteriori from the FWI velocity and density models shows a well-focused image with however local differences with the velocity model near the sea bed where density might have absorbed elastic effects. The FWI models are finally assessed against time-domain synthetic seismogram modelling performed with the same frequency-domain modelling engine used for FWI.

**Key words:** Body waves; Computational seismology; Controlled source seismology; Inverse theory; Seismic anisotropy; Wave propagation.

## INTRODUCTION

Besides compressional wave speed, imaging different earth properties (Poisson ratio, density, attenuation, anisotropy, etc.) from broadband long-offset seismic data is a topical issue in full waveform inversion (FWI) as these seismic properties are helpful to decrease the ambiguities about the rocks and fluid types and draw inferences on petrophysical attributes through downscaling approaches at meso- and microscales (Dupuy et al. 2016). However, multiparameter FWI is challenging either because of the inherent trade-off between parameters or the different imprint of the parameters in the seismic response in terms of strength or nature (kinematic versus dynamic parameters; Operto et al. 2013).

In FWI, the sensitivity of the seismic response to a local parameter perturbation is represented by the wavefield emitted by the seismic source and scattered by this perturbation, referred to as the

partial derivative wavefield (Pratt et al. 1998). The secondary source formed by the perturbation has a radiation pattern, which depends on the parameter type and controls the amplitude versus scattering angle variations of the scattered wavefield. The radiation pattern is the only term in the gradient of the FWI misfit function that allows the inversion to discriminate the contribution of different parameter types in the seismic response. Moreover, the spectral component (i.e. wavenumber vector) of the model perturbation constrained by a seismic event is related to the local wavelength and the scattering angle by

$$ \mathbf{k} = \frac{2\omega}{c} \cos(\theta/2) (\cos \phi, \sin \phi), \qquad (1) $$

where ω is the angular frequency, c the local wave speed, θ the scattering angle and φ the dip angle (Miller et al. 1987; Wu & Toksöz

<footnote n="1">¹Université Côte d'Azur, CNRS, IRD, OCA, Géoazur, F 06905 SOPHIA ANTIPOLIS Cedex, France. E-mail: operto@geoazur.unice.fr</footnote>

<footnote n="2">²Université Côte d'Azur, Observatoire de la Côte d'Azur (OCA), F 06394 NICE Cedex 4, France</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Resolution and trade-off assessment. (a) From left to right, true model contains a *V*₀, *ρ*, *Q* mono-parameter inclusion. Three mono-parameter FWIs are performed to update the parameter in question with 80 iterations. (b) Same as (a) for 10 FWI iterations. (c) True model contains a multiparameter inclusion. Multiparameter FWI is performed to jointly update *V*₀ (left), *ρ* (middle) and *Q* (right) with 80 iterations. Maps show one slice across the reconstructed inclusion. Right panel shows the direct comparison between the *x*, *y*, *z* true (black line) and reconstructed profiles intersecting the centre of the inclusion. Bottom panel shows the true (black) and updated (grey) wavenumber spectrum of a profile intersecting the inclusion.</figcaption></figure>

1987; Lambaré *et al.* 2003). This expression shows that when broadband sources and wide-aperture long-offset surveys sample a wide range of θ angles, then one **k** is sampled in a redundant way by many (ω, θ) pairs. While this redundancy is often decimated in efficient mono-parameter frequency-domain FWI for the sake of computational efficiency (Sirgue & Pratt 2004), it may be useful to preserve it during multiparameter FWI to mitigate parameter cross-talks (Podgornova *et al.* 2015). The reason is that the amplitude versus θ variations of the partial derivative wavefields will be more fully exploited during each **k** reconstruction. In the opposite extreme

case, when decimation is fully removed during a mono-frequency inversion step, one **k** is imaged with one θ, hence making FWI blind to the radiation pattern.

Furthermore, FWI is generally applied on narrow frequency band (cut-off frequency of the order of 5–10 Hz) either for computational or nonlinearity reasons or because it is used as a velocity model building tool for migration rather than for structural and lithological interpretation. Moreover, frequency continuation strategies to mitigate cycle skipping further contribute to narrow the frequency band used during each multiscale step (e.g. Warner *et al.* 2013). All
<figure><img src="image_6.png" /><figcaption>**Figure 2.** OBC acquisition. The black dots and black circles represent the shots and the hydrophones, respectively. The white and grey circles point two receivers, whose gathers are shown in Fig. 3. A depth slice across the low-velocity gas cloud (black area) is superimposed in transparency. The parallel (Y-oriented) shot profile above the first receiver (white circle) intersects the gas cloud, while the one above the second receiver (grey circle) is away from the gas cloud. The white and grey stars point the positions of two sonic well logs for vertical wave speed and density, respectively (see also Figs 9 and 14).</figcaption></figure>

these strategies contribute to decrease the wavenumber redundancy and hence the well-posedness of the multiparameter inversion.

In practice, multiparameter FWI may not be only driven by radiation patterns to discriminate between different parameter classes because they control the amplitudes of the partial derivative wavefields, not their traveltimes. At wide apertures and low frequencies, FWI primarily updates the long wavelengths of the waves peeds to fit the phase of the wavefields before updating parameters more closely related to amplitudes such as attenuation and density. When moving at higher frequencies and shorter scattering angles, parameter cross-talks may become a more significant issue when the inversion attempts to fit the amplitudes of the short-spread reflections to update the short wavelengths of the subsurface. In short, it is likely that the relative influence of each parameter on traveltime and amplitudes during wave propagation establish a natural hierarchy among them that reduces the risk of cross-talk at low frequencies.

This study presents a real 3-D data case study of multiparameter FWI where we jointly update the vertical wave speed ($V_0$), the density ($\rho$) and the quality factor $Q$ in the viscoacoustic vertical transverse isotropic (VTI) approximation and the 3.5–10 Hz frequency band, while the Thomsen's parameters $\delta$ and $\epsilon$ are used in a passive way. Data have been collected during a well-documented wide-azimuth ocean-bottom cable (OBC) survey in the North Sea (Sirgue *et al.* 2010; Barkved *et al.* 2010). The OBC geometry allows for the recording of diving waves and reflections up to the critical incidence. The geology is characterized by soft sediments in the near surface and a gas cloud above a top hard chalk reflector. The gas cloud as well as the top hard chalk reflector generates E106 03 significant impedance contrasts that allow for the recording of energetic short spread reflections amenable to density update. Moreover, the soft sediments and the gas cloud should make the imprint of the attenuation significant in the seismic response. We perform FWI in the frequency domain where accurate attenuation mechanisms can be easily implemented without computational overheads. 3-D frequency-domain FWI is computationally tractable as long as the inversion is limited to a few discrete frequencies (Operto *et al.* 2015; Amestoy *et al.* 2016). As explained above, this frequency decimation reduces the redundancy with which the spectral components of the subsurface are sampled by different ($\omega, \theta$) pairs and hence

questions the relevance of multiparameter frequency-domain FWI to manage parameter cross-talks. We apply FWI using different frequency samplings and groupings and different inversion tunings to confront the experimental results with the theoretical framework above reviewed. Our results are validated against well logs and full waveform seismic modelling in the frequency and time domains.

There have been several investigations to assess either the feasibility of density or impedance (Jeong *et al.* 2012; Przebindowska *et al.* 2012; Bai & Yingst 2014; Qin & Lambaré 2016; Yang *et al.* 2016b) or attenuation (e.g. Hicks & Pratt 2001; Askan *et al.* 2007; Malinowski *et al.* 2011; Kamei & Pratt 2013; Kurzmann *et al.* 2013; Takam Takougang & Calvert 2013; Stopin *et al.* 2016) imaging by FWI. However, there are very few synthetic or real-data case studies where these two parameters are jointly updated at the notable exception of Prieux *et al.* (2013) and Métivier *et al.* (2015). This gap might sound surprising as these two parameters control amplitudes.

The density has a radiation pattern which shows a maximum sensitivity of the partial derivative wavefield at the zero scattering angle when it is associated with the wave speed in an acoustic parametrization, this sensitivity gradually decreasing with scattering angle to zero at the 180° scattering angle (e.g. Operto *et al.* 2013, their fig. 1). This radiation pattern raises the issue of the trade-off with velocities at short- to intermediate-scattering angles discussed above. This trade-off can be mitigated by moving to a velocity-impedance parametrization at the expense of the resolution with which velocity is updated (Prieux *et al.* 2013). Moreover, as the radiation pattern can be seen as a $\theta$-dependent weighting operator in the FWI misfit gradient, the high wavenumbers of the density model will tend to be updated before the smaller ones and hence will break down the scale continuation from the low wavenumbers to the higher ones fostered by the frequency continuation in the data space. This prompts Yang *et al.* (2016b) to design an FWI pre-conditioner which balances the effects of the density radiation pattern.

The radiation pattern of $Q$ is isotropic as the one of the wave speed in the viscoacoustic approximation (Malinowski *et al.* 2011). This can be easily understood as viscous effects are implemented in the imaginary part of the complex-valued wave speed. Accordingly, the partial derivative of the wavefields with respect to wave speed and $Q$ are related to a Hilbert transform, which prompts Mulder & Hak (2009) to conclude that the joint reconstruction of wave speed and attenuation is ambiguous in linear reflection waveform inversion. This is consistent with the conclusions of Ribodetti *et al.* (2000) who showed that the Hessian of the viscoacoustic ray + Born waveform inversion is singular when a reflector is sampled by one-side reflections (from above or beneath) as it might be the case in surface acquisitions. Only a two-side illumination as in medical imaging allows to remove the Hessian singularity. However, Hak & Mulder (2011) show that the velocity and attenuation can be jointly imaged in nonlinear FWI provided that the viscous model honours causality and a large number of iterations are performed. Furthermore, this conclusion was supported by several synthetic or real data applications (see the above-mentioned reference list). Kurzmann *et al.* (2013) performed a sensitivity analysis of mono-parameter FWI for velocity to attenuation effects. They show that reliable velocity models are built provided that a crude homogeneous or smooth $Q$ model can be used in a passive way. Note that this statement does not mean that attenuation belongs to the null space of the FWI (see the discussion on this issue in Prieux *et al.* (2013) and Operto *et al.* (2013)).

This study is not an attempt to investigate different subsurface parametrization for FWI and we will limit ourselves to a quite standard ($V_0, \rho, Q, \delta, \epsilon$) parametrization. Our motivation beyond
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Common-receiver gathers. The receiver positions are denoted by the white (a) and grey (b) circles in Fig. 2. (a) The parallel shot profile cross-cuts the gas cloud. (b) Same as (a) for a shot profile away from the gas cloud. The yellow arrows point the reflections from the top of the gas cloud and the top hard chalk reflector at around 1.8 and 2.7 s two-way traveltimes. This phase identification is checked in Fig. 4 where two inline common-receiver gathers are superimposed on the coincident sections of an FWI model (FWI4) after depth-to-time conversion. The bottom panels show the two gathers with a reduction velocity of 2.5 km s⁻¹. The dashed white lines highlight at long offsets the tangency relationship between the diving waves from above the gas cloud and the reflection from the top of this zone. The black arrow points a possible refracted arrival from the reservoir.</figcaption></figure>

this parametrization is the isotropic radiation pattern of $V_0$ which guarantees a broad-band reconstruction of this first-order parameter. Moreover, we believe that the long wavelengths contained in the available $\epsilon$ model are accurate enough to use this parameter in a passive way during inversion (Operto et al. 2015).

In the first part of this study, we first review the main ingredients of our frequency-domain FWI implementation. The second part describes the real data case study. After a description of the acquisition and data, we describe our FWI experimental setup in terms of frequency management and inversion tuning. Then we discuss the relevance of the $V_0$, $\rho$ and $Q$ models updated during each test. Basically, we show how we gradually improve the reconstruction of $V_0$ in terms of signal-to-noise ratio and focusing as we perform frequency grouping at each multiscale step and use $\rho$ and $Q$ as optimization parameters rather than passive parameters. We show how the acquisition footprint tends to be focused in the $\rho$ and $Q$ updates, hence cleaning the $V_0$ model from their imprint. We also show the dispersion with which $Q$ and $\rho$ are updated depending on the inversion tuning. However, this dispersion is not passed on $V_0$, which prompts us to conclude that $V_0-\rho$ and $V_0-Q$ trade-offs are quite limited for this case study. In the final part, we perform a comparative analysis of the data fit achieved by the FWI models developed in this study. This comparative analysis provides clear insights on the hierarchical role of $V_0$, $\rho$ and $Q$ to fit different parts of the wavefields as well as the importance of the frequency management during frequency-domain FWI to fit the dispersive part of the wavefield. This study supports the feasibility and computational efficiency of multiparameter frequency domain FWI of wide-azimuth

stationary-receiver surveys as we limit the inversion to six to ten discrete frequencies in the 3.5–10 Hz frequency band.

## MULTIPARAMETER FREQUENCY-DOMAIN FWI

### Forward problem

Frequency-domain seismic modelling is a boundary-value problem, which requires solving a large and sparse system of linear equation per frequency, whose right-hand sides are the seismic sources and the solution are the monochromatic wavefields (Marfurt 1984). We solve this system with a Gauss elimination technique suitable for sparse matrices to process efficiently the reciprocal sources of the 3-D OBC survey (Operto et al. 2015; Amestoy et al. 2016). However, the computational burden of the lower-upper (LU) decomposition, which is performed once and for all before the forward and backward substitution steps, requires to limit the inversion to a few discrete frequencies.

We discretize the VTI viscoacoustic wave equation with the finite-different method of Operto et al. (2014), which has been designed to mitigate the memory burden associated with sparse direct solvers while providing accurate solutions for a discretization rule of four grid points per wavelength, which is optimal for the half-wavelength resolution of FWI.

The discrete VTI viscoacoustic wave equation can be written in matrix form as

$$A_h \mathbf{p}_h = \mathbf{s}', \qquad (2)$$
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Vertical sections of a $V_0$ FWI model (Table 1, FWI4) plotted as a function of two-way traveltimes. The high-wavenumber components of the $V_0$ model are superimposed to highlight the main reflectors (top of the gas cloud at around 1.6 s, top hard chalk reflector at 2.7–3 s and base cretaceous reflector at 3.5–4 s). (a) The section intersects the velocity log (Fig. 2, white star) and hence is away from the gas cloud. (b) The section intersects the gas cloud. In the bottom panels, two common-receiver gathers are superimposed in transparency to correlate the reflections identified in these gathers with the reflectivity of the FWI model. See the text for interpretation.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Initial models for FWI ($V_0$). Depth slices at (a) 175 m depth, (b) 500 m depth, (c) 1 km depth and (d) 3.35 km depth.</figcaption></figure>

$$ \mathbf{p}_v = \mathbf{A}_v \mathbf{p}_h + \mathbf{s}'' \quad (3) $$ where $\mathbf{p}_h, \mathbf{p}_v$ and $\mathbf{p}$ are the so-called horizontal, vertical and physical pressure vectors. The operators $\mathbf{A}_h$ and $\mathbf{A}_v$ discretized in matrices

$$ \mathbf{p} = \frac{1}{3} (2\mathbf{p}_h + \mathbf{p}_v), \quad (4) $$
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Initial models for FWI ($V_0$). Vertical section at (a) $x = 5.575$ km, (b) $x = 6.25$ km and (c) $x = 9.5$ km.</figcaption></figure>

$\mathbf{A}_h$ and $\mathbf{A}_v$ are given by

$$A_h = \omega^2 \left[ \frac{\omega^2}{\kappa_0} + (1 + 2\epsilon)(\mathcal{X} + \mathcal{Y}) + \sqrt{1 + 2\delta} \mathcal{Z} \frac{1}{\sqrt{1 + 2\delta}} \right] + 2\sqrt{1 + 2\delta} \mathcal{Z} \frac{\kappa_0(\epsilon - \delta)}{\sqrt{1 + 2\delta}} (\mathcal{X} + \mathcal{Y}) \quad (5)$$

$$A_v = \frac{1}{\sqrt{1 + 2\delta}} + \frac{2(\epsilon - \delta)\kappa_0}{\omega^2 \sqrt{1 + 2\delta}} (\mathcal{X} + \mathcal{Y}) \quad (6)$$

where $\kappa_0 = \rho V_0^2$, $\rho$ is density, $V_0$ is the vertical wave speed, $\omega$ is the angular frequency, $\delta$ and $\epsilon$ are Thomsen's parameters (Thomsen 1986). Differential operators $\mathcal{X}$, $\mathcal{Y}$ and $\mathcal{Z}$ are given by $\partial_{\tilde{x}} b \partial_{\tilde{x}}$, $\partial_{\tilde{y}} b \partial_{\tilde{y}}$ and $\partial_{\tilde{z}} b \partial_{\tilde{z}}$, respectively, where $b = 1/\rho$ is buoyancy and $(\tilde{x}, \tilde{y}, \tilde{z})$ define a complex-valued coordinate system in which absorbing boundary condition are implemented (Operto et al. 2007). The expression of the right-hand sides $\mathbf{s}'$ and $\mathbf{s}''$ are given in Operto et al. (2014).

We implement viscous effects assuming that $Q$ is constant over the range of frequencies considered in this study. Under this assumption, the Kolsky-Futterman model (Kolsky 1956; Futterman 1962) honours causality and leads to the following complex-valued vertical wave speed

$$\frac{1}{\tilde{V}_0(\omega)} = \frac{1}{V_0(\omega)} \left( 1 - \frac{1}{\pi Q} \ln \left| \frac{\omega}{\omega_r} \right| \right), \quad (7)$$

where $\omega_r$ is a reference angular frequency (Aki & Richards 2002, page 61).

## Inverse problem

The FWI is performed by a classical iterative least-squares minimization of the difference between recorded and modelled monochromatic pressure:

$$\min_{\mathbf{m}} C(\mathbf{m}) = \min_{\mathbf{s}} \sum_{\omega} \|\Delta\mathbf{d}_{s,\omega}\|^2, \quad (8)$$

where $\mathbf{m}$ gathers the subsurface model parameters of different kinds and $\Delta\mathbf{d}_{s,\omega} = \mathbf{d}_{s,\omega}(\mathbf{m}) - \mathbf{d}_{s,\omega}^*$ denotes the monochromatic data residuals for source $s$ with $\mathbf{d}_{s,\omega}(\mathbf{m})$ and $\mathbf{d}_{s,\omega}^*$ the modelled and recorded data, respectively. In this study, we follow a frequency continuation approach by proceeding from the low frequencies to the higher ones to mitigate cycle skipping (e.g. Sirgue & Pratt 2004). Moreover, in the multiparameter framework, we may need to preserve some degree of redundancy in the wavenumber sampling during each multiscale step through the summation over a subset of frequencies in eq. (8). The subset of frequencies involved during one multiscale step will be referred to as frequency group in the following. The subsurface model updated at iteration $k + 1$ is given by

$$\mathbf{m}_{k+1} = \mathbf{m}_k - \gamma_k \mathbf{H}_k \nabla_m C_k, \quad (9)$$

where $\mathbf{H}_k$ and $\gamma_k$ denotes an approximate Hessian and the step length, respectively.

Operto et al. (2015) show that the gradient of the misfit function, $\nabla_m C_k$ can be approximated by

$$\nabla C_m \approx \sum_s \sum_{\omega} \Re \left\{ \left( \frac{\partial \mathbf{A}_h(\omega, \mathbf{m})}{\partial m} \mathbf{p}_h(s, \omega) \right)^{\dagger} \mathbf{a}_1(s, \omega) \right\}, \quad (10)$$

where the adjoint wavefields $\mathbf{a}_1$ satisfies

$$\mathbf{A}_h^{\dagger}(\omega, \mathbf{m}) \mathbf{a}_1(s, \omega) = \frac{1}{3} (\mathbf{A}_v^{\dagger}(\omega, \mathbf{m}) + 2\mathbf{I}) \mathbf{R}^{\mathrm{T}} \Delta \mathbf{d}_{s,\omega}, \quad (11)$$

where $\mathbf{R}^{\mathrm{T}}$ augment with 0 are the residuals at the receiver positions in the full computational domain and $\mathbf{I}$ is the identity matrix.

Using the SEISCORE toolbox (Métivier & Brossier 2016), we perform the optimization with either a pre-conditioned steepest-descent or l-BFGS methods depending which frequency is processed (see the next section for details). As pre-conditioner, we use the diagonal elements of the so-called pseudo-Hessian matrix (Shin et al. 2001), the aim of which is to balance the gradient with respect to depth by removing geometrical spreading effects. The expression of the gradient pre-conditioner is given by

$$\mathbf{H} = 1 / [\mathbf{P} + \beta \max(\mathbf{P})], \quad (12)$$

where

$$\mathbf{P} = \mathrm{diag} \sum_s \sum_{\omega} \Re \left\{ \left( \frac{\partial \mathbf{A}_h(\omega, \mathbf{m})}{\partial m} \mathbf{p}_h(s, \omega) \right)^{\dagger} \times \left( \frac{\partial \mathbf{A}_h(\omega \mathbf{m}, \cdot)}{\partial m} \mathbf{p}_h(s, \omega) \right) \right\}$$

and $(\partial \mathbf{A}_h(\omega, \mathbf{m})/\partial m)\mathbf{p}_h(s, \omega)$ are the so-called virtual sources (Pratt et al. 1998) whose radiation pattern are $\partial \mathbf{A}_h(\omega, \mathbf{m})/\partial m$. The pre-whitening factor $\beta$, which damps the balancing with respect to depth of the gradient, can be adapted to each parameter class. Furthermore, we scale each parameter class by a constant factor to balance the relative amplitudes of the multiparameter gradient during the early FWI iterations when the Hessian action is not yet accurately estimated by l-BFGS.

As in Operto et al. (2015), we do not apply regularization and we update the source signature at each FWI iteration by matching in the least-squares sense the monochromatic Green functions with the real data (Pratt 1999).
<figure><img src="image_11.png" /><figcaption>**Figure 7.** $V_0$ FWI models. Depth slices at (a) 175 m depth across glacial sand channel deposits, (b) 500 m depth across a wide low-velocity zone and scrapes left by drifting icebergs on the paleo-seafloor, (c) 1 km depth across the gas cloud and (d) 3.35 km depth across the base cretaceous reflector. In (a)–(d), each panel is labelled by the test number (Table 1).</figcaption></figure>

## Subsurface parametrization

We use the ($V_0$, $\rho$, $Q$, $\delta$, $\epsilon$) parametrization to perform FWI. For this parametrization, a $V_0$ and $\epsilon$ perturbations scatter energy at all scattering angles and wide angles, respectively (Gholami et al. 2013). This implies that a broad-band $V_0$ can be built with this parametrization, while only the long wavelengths of $\epsilon$ will be attainable. Therefore, $\epsilon$ can be used as a passive parameter during FWI as long as its long wavelengths have been accurately built during a former stage (traveltime tomography or others). In this study, we use the same VTI initial model as that used in Operto et al. (2015), who showed that the $\epsilon$ model is accurate enough to match first-arrival traveltimes and post-critical reflections from the reservoir during seismic modelling with a kinematic error small enough to prevent significant cycle skipping at a starting frequency of 3.5 Hz. This prompts us to

use $\epsilon$ as a passive parameter in this study and focus on the update of $V_0$, $\rho$ and $Q$. As in Operto et al. (2015), we build a starting $\rho$ model from the initial $V_0$ model using a Gardner law.

In order to have a rough idea of the FWI results that should be expected with this parametrization before the real data application, we perform multiparameter FWI for a toy model corresponding to a spherical inclusion embedded in a homogeneous background (Fig. 1). The sources and the receivers are deployed along the six faces of the computational domain leading to a complete illumination of the inclusion in terms of scattering angles. Five frequencies between 4 and 8 Hz are simultaneously inverted starting from the homogeneous background model. The $V_0$, $\rho$ and $Q$ values are 1600 m s⁻¹, 1000 kg m⁻³ and 150 in the background model and 1800, m s⁻¹, 1200 kg m⁻³ and 10 in the inclusion. The intrinsic
<figure><img src="image_12.png" /><figcaption>**Figure 8.** $V_0$ FWI models. Vertical sections (a) across gas cloud ($x = 5.575$ km), (b) near the periphery of the gas cloud ($x = 6.25$ km), (c) away from the gas cloud ($x = 9.5$ km). The reflectivity component of the velocity model, obtained by high-pass time-domain filtering after depth-to-time conversion, is superimposed in transparency on the velocity field to highlight the main structural discontinuities. In (a)–(c), each panel is labelled by the test number (Table 1).</figcaption></figure>

resolution with which each parameter can be reconstructed is assessed by performing three separate mono-parameter inversions for $V_0$, $\rho$ and $Q$ where the corresponding true models contains a $V_0$, $\rho$ and $Q$ inclusion, respectively (Fig. 1a). The results show as expected that the full spectrum of $V_0$ is reconstructed. In contrast, the inversion fails to update the small wavenumbers of $\rho$ according to its radiation pattern. Moreover, looking at the FWI results after a limited number of iterations shows that the inversion starts updating the high wavenumbers of $\rho$ from the high amplitudes of the partial derivative wavefields at short scattering angles before updating the intermediate wavenumbers from the weaker amplitudes of the partial derivative wavefields at wider scattering angles (Fig. 1b). $Q$ is updated with a resolution close to that of $V_0$ although the inversion has more difficulties to update its high wavenumbers. When moving to the joint reconstruction of the three parameter classes in a multiparameter inclusion model, $V_0$ and $Q$ are reconstructed in a very similar way without any cross-talk artefacts (Fig. 1c). These results are consistent with the 2-D ones of Malinowski *et al.* (2011). This is also the case when the five frequencies are processed sequentially rather than simultaneously. In contrast, the inversion failed to

update the intermediate wavenumbers of $\rho$ due to too limited sensitivity. This failure is even more severe in the case of successive mono-frequency inversions.

We would like to mention that the gradient of the frequency-domain FWI may be sensitive to the numerical scheme used to discretize the radiation pattern matrices $\partial A_h/\partial m$. In our code, the matrix $A_h$ used to generate the radiation pattern matrices is built with a basic seven-point 3-D stencil for implementation convenience, while the matrix $A_h$ built for seismic modelling is built with a 27-point stencil (Operto *et al.* 2014). This potentially makes the $\rho$ updates particularly sensitive to this discretization as it is embedded in the operators $X$, $Y$ and $Z$, eq. (5). In contrast, gradients of time-domain FWI can be made less sensitive to space discretization by using a velocity-stress state equation written in pseudo-conservative form (Vigh *et al.* 2014; Yang *et al.* 2016c). We check this statement by comparing the results of our frequency-domain FWI code with a time-domain implementation. In the latter case, a broader wavenumber spectrum of $\rho$ was built compared to Fig. 1(a) at the expense of a very large number of iterations. Ongoing work aims to build the $\partial A_h/\partial m$ from a 27-point stencil to regularize the $\rho$ update by decreasing the footprint of the discretization.

## APPLICATION TO OBC DATA

### Study area and data

The subsurface target is a North Sea oil field in a 70 m depth shallow water environment (Barkved *et al.* 2010; Sirgue *et al.* 2010; Haller *et al.* 2016; Yang *et al.* 2016a). The overburden is characterized by soft sediments in the near surface and a low-velocity gas cloud above the chalk reservoir, which is located at around 2.5 km depth. The gas cloud, whose main zone of influence is shown on the survey map in Fig. 2, makes seismic imaging at the reservoir depths challenging. The parallel geometry of the OBC acquisition consists of 2302 hydrophones which record 49 954 explosive sources located 5 m below the sea surface (Fig. 2). The seismic sources cover an area of 145 km² leading to a maximum offset of 14.5 km. In this study, the maximum depth of investigation is 4.5 km. Two common-receiver gathers are shown in Fig. 3. In Fig. 4, we link the main reflections contained in these gathers to the main reflectors of a depth-to-time converted FWI velocity model developed in this study. This correlation shows that the arrivals recorded at around 1.6 and 2.8 s two-way traveltimes in Fig. 3, yellow arrows, are the reflections from the top of the gas cloud and the top hard chalk reflector. The more discontinuous and shingling pattern of the reflections recorded along a profile intersecting the centre of the gas cloud highlights the complex interaction of the wavefield with it (Fig. 3a). The diving waves mainly propagate above the gas cloud as shown by their tangency relationship at long offsets with the reflection from the top of the gas cloud (Fig. 3b, white dashed lines). The refracted waves from the top hard chalk reflector are not recorded as a first arrival but may be identified as a second arrival (Fig. 3a, black arrow). The wavefield also shows two categories of surface waves: Scholte waves (Fig. 3, S), that will be processed as noise by viscoacoustic FWI, and guided waves propagating as leaking modes with phase velocities higher than the water wave speed (Fig. 3, GW). In soft sedimentary environment, these leaking modes are mostly composed of reverberating P-wave reflections in the waveguide bounded by the free surface and the bottom of the weathering layer, and hence can be used to constrain the compressional wave speeds in the near surface by inversion of dispersion
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Validation against well log ($V_0$). (a) Comparison between the well log (black line) and the coincident profiles of the initial (red line) and final FWI (blue line) models for the six FWI tests. The FWI velocity model around the well log is shown in the background. Each panel is labelled by the test number (Table 1). (b) Parallel section of the FWI $V_0$ model (FWI4) intersecting the well log. The black line represents the well log (see Fig. 2, white star, for its position on the survey map).</figcaption></figure>

curves (Roth *et al.* 1998; Shtivelman 2004; Boiero *et al.* 2013, 2014; Wiarda *et al.* 2014). These high-amplitude dispersive wave trains have a shingling pattern that partially overprints deep post-critical reflections. Similar to the surface-related multiples, we leave these quasi-$P$ wave trains in the data and invert them by viscoacoustic FWI.

## FWI experimental setup

We perform six FWI tests, referred to as FWI1 to FWI6, during which we use different frequency groupings and samplings, and $Q$ scaling (Table 1).

For all tests, we use a multiscale frequency continuation approach and adapt the grid interval to the maximum frequency of the current frequency group to optimize the computational cost and regularize the inversion by re-parametrization. We use grid intervals of 70, 50 and 35 m for $f_{max} < 5$ Hz, $5 Hz < f_{max} \le 7$ Hz and $7 Hz < f_{max} \le 10$ Hz, respectively (Table 1).

Between 3.5 and 4.5 Hz when the signal-to-noise ratio is limited, we perform successive mono-frequency inversions for $V_0$ with a pre-conditioned steepest-descent algorithm. For a given sets of discrete frequencies, successive mono-frequency inversions, by opposition to joint inversion of multiple frequencies, are safer to prevent cycle skipping, while a steepest-descent algorithm is less prone to boosting high-wavenumber artefacts in the FWI gradients. For frequencies greater or equal to 5Hz, we switch to a pre-conditioned I-BFGS algorithm.

In the first five tests, we use 10–11 frequencies between 3.5 and 10 Hz as in Operto *et al.* (2015), while we use a coarser frequency sampling (six frequencies) for the last test as in Amestoy *et al.* (2016) (Table 1).

The first two tests (FWI1 and FWI2) are mono-parameter FWI for $V_0$ and the frequencies are processed sequentially without grouping until the 10 Hz frequency. In FWI1, no viscous effects are implemented during forward modelling ($Q = 1000$), while a homogeneous background with $Q = 200$ below the sea floor is used during FWI2. Comparing the results of these two tests will provide first insights on the footprint of attenuation and assess whether even crude $Q$ models improve the FWI results as concluded by Kurzmann *et al.* (2013).

In the last four tests (FWI3–FWI6), we perform ($V_0, \rho, Q$) multiparameter FWI from the 5 Hz frequency. In all these tests, we start inverting the 5 Hz frequency component, although mono-frequency inversion is not recommended for multiparameter update. At end of the 5 Hz inversion, we check that the $V_0$ model is consistent with the one obtained by Operto *et al.* (2015). During FWI3, we perform seven successive multiparameter mono-frequency inversions between 5 and 10 Hz. In theory, this test should not give reliable results since the wavenumber redundancy is fully removed during each multiscale step. During FWI4 and FWI5, we perform two successive multifrequency inversions with the frequency groups [5.2 Hz, 5.8 Hz, 6.4 Hz, 7Hz] and [8 Hz, 10Hz] after the 5 Hz inversion. Comparing the results of FWI3 with those of FWI4 and FWI5 will highlight the benefit of frequency grouping for multiparameter inversion. Compared to FWI4, we modified the scaling of $Q$ during FWI5 to favour its shallow update at the expense of the deeper one (Table 1). Comparing the results of FWI4 and FWI5 will give some qualitative insights on the dispersion with which $Q$ can be imaged and the impact on this variations on the $V_0$ and $\rho$ updates. FWI6 is similar to FWI5 test except that the two multifrequency inversions are performed with the frequency groups [5 Hz, 7 Hz] and [7 Hz, 10 Hz] with one frequency overlap. The aim of this last test is to
<figure><img src="image_14.png" /><figcaption>**Figure 10.** FWI perturbation models (Q). Depth slices at (a) 175 m depth, (b) 500 m depth, (c) 1 km depth and (d) 3.35 km depth. In (a)–(d), each panel is labelled by the test number (Table 1). The geological features intersected by each slice are reviewed in the caption of Fig. 7.</figcaption></figure>

gain some insights on the sensitivity of the multiparameter updates to the frequency sampling in a given band, keeping in mind that the computational cost of frequency-domain FWI scales to the number of frequencies involved in the inversion.

We choose the values of the constant scaling factors so that they have the order of magnitude of the physical parameters. This scaling balances the amplitudes of the multiparameter gradient by removing the imprint of these orders of magnitude (or physical units). Then, we tune the quantity of descent at first iteration through an absolute scaling factor $\alpha$. We set $\alpha$ so that the multiparameter gradient after multiplication with this factor has the order of magnitude of a small fraction of the scaled parameters. This is implemented through

$$ \alpha = \frac{\epsilon_g \bar{m}_{mg}}{\max(\nabla_m \bar{C})}, \quad (13) $$

where $\nabla_m \bar{C}$ denotes the scaled multiparameter gradient, $\epsilon_g$ a small pre-set value and $\bar{m}_{mg}$ the value of the scaled parameter at the subsurface position where the maximum of the scaled gradient was found. Moreover, we check that our first guess of the relative scaling factors makes the quantity of descent to be controlled by the first-order kinematic parameter $V_0$ (namely, the maximum value of the scaled multiparameter gradient must be reached for $V_0$). With these guidelines in mind, we refine if necessary our initial guess of the scaling factors by trial-and-error to better balance the relative updates of the multiple parameter classes during the early iterations.

We use a reference frequency of 50 Hz in the Kolsky–Futterman attenuation model, eq. (7). This value has been chosen according to the mono-parameter FWI results presented by Operto et al. (2015).

They found a 40 ms time-shift between source signatures estimated from short offset and long offset data, respectively. Among the possible interpretations, they propose that this traveltime shift was generated by the underestimation of attenuation effects during seismic modelling and inversion. By comparing seismograms computed in a simple sedimentary two-layer model with and without attenuation, they found that a reference frequency of the order of 50 Hz was generating a 40ms time shift between the two sets of seismograms.

We did not design a quantitative stopping criterion of iterations. We stopped the iterations when the decrease of the misfit function and the model update become negligible with the risk to over fit the data. Note that we never met a line search failure with the SEISCORE toolbox.

## FWI results

### $V_0$ models

The initial $V_0$ model (Figs 5 and 6) and the FWI models (Figs 7 and 8) are shown along depth slices at 175 m, 500 m, 1 km and 3.35 km depth and vertical sections across the centre of the gas cloud ($x = 5.575$ km), its periphery ($x = 6.25$ km) and away from the gas cloud ($x = 9.5$ km). From top to bottom, the depth slices intersect glacial sand channel deposits, a wide low-velocity zone intersected by scrapes left by drifting icebergs on the paleo-seafloor, the gas cloud and the base cretaceous reflector, respectively. These structures have been well documented by FWI in Sirgue et al. (2010), Barkved et al. (2010), Operto et al. (2015), Yang et al. (2016a) and
<figure><img src="image_15.png" /><figcaption>**Figure 11.** FWI perturbation models ($Q$). Vertical sections (a) across gas cloud ($x = 5.575$ km), (b) near the periphery of the gas cloud ($x = 6.25$ km) and (c) away from the gas cloud ($x = 9.5$ km). In (a)–(c), each panel is labelled by the test number (Table 1).</figcaption></figure>

Amestoy et al. (2016) and the reader is referred to these publications for the assessment of the results presented in this study. Note that, in all the figures showing the vertical sections of the FWI $V_0$ models, we superimpose on the $V_0$ models their reflectivity component, generated by high-pass time-domain filtering after depth-to-time conversion, for an easier identification of the structural discontinuities reconstructed by FWI.

The $V_0$ model built by mono-parameter FWI without attenuation (FWI1) shows a strong acquisition footprint at 175 m depth (Fig. 7a), significant artefacts at 500 m depth (Figs 7b and 8b) and overall noisy images. The results of the second mono-parameter inversion (FWI2), when a homogeneous background with $Q = 200$

is used, confirms that using a crude $Q$ model in a passive way significantly mitigates the previous artefacts. For example, the acquisition footprint is reduced at 175 m depth and the artefacts at 500 m depth are removed (Figs 7a and b). The setup of the FWI2 test is similar to the one used in Operto et al. (2015) except that we switch to 1-BFGS from the 5 Hz frequency instead of using a steepest-descent algorithm for all the frequencies. Comparing the FWI2 $V_0$ models with those of Operto et al. (2015) shows how the deconvolution action of the Hessian performed by 1-BFGS has injected high-frequency noise in the models. Regularization would have been necessary to build $V_0$ models of similar quality as those of Operto et al. (2015). This illustrates the usual compromise between
<figure><img src="image_16.png" /><figcaption>**Figure 12.** FWI perturbation models ($\rho$). Depth slices at (a) 175 m depth, (b) 500 m depth, (c) 1 km depth and (d) 3.35 km depth. In (a)–(d), each panel is labelled by the test number (Table 1). The geological features intersected by each slice are reviewed in the captions of Fig. 7.</figcaption></figure>

resolution and signal-to-noise ratio that needs to be found during deconvolution processes.

When moving to the first multiparameter inversion by successive mono-frequency inversions (FWI3), the $V_0$ model keeps on slightly improving in terms of acquisition footprint and noise mitigation, although comparison with the results of Operto et al. (2015) shows that they clearly remain perfectible in terms of deep reflector focusing and signal-to-noise ratio (Figs 7 and 8). Moving from mono-frequency to multifrequency inversion during FWI4 significantly improves the $V_0$ model, which is similar to that of Operto et al. (2015): the acquisition footprint has been mitigated at all depths and the noise in the FWI model is weak although no regularization was applied during the I-BFGS optimization (Figs 7 and 8). Although the $V_0$ models built during mono-frequency FWI3 and multifrequency FWI4 inversions are consistent from the structural viewpoint, the values of the velocities show some differences, for example, across the gas cloud and the base cretaceous reflector (Figs 7c and d). Overall, the velocity perturbations built by mono-frequency inversions are sharper than those built by multifrequency inversions because mono-frequency inversions are more prone to over fit less-redundant monochromatic data at the expense of the signal-to-noise ratio. Another possible reason is related to parameter cross-talks between $V_0$, $\rho$ and $Q$ which can have a different strength depending to which extent the data redundancy is decimated during each multiscale step.

The $V_0$ model obtained during FWI5 when the $Q$ scaling was modified compared to the previous tests (Table 1), is similar to the one of FWI4, except that the acquisition footprint has been less efficiently mitigated at depths (e.g. compare the depth slices across

the gas cloud in Fig. 7(c) and the vertical sections above the gas cloud between 0.5 and 1.5 km depth in Fig. 8). The subsequent analysis of the $Q$ FWI models will justify these differences. Finally, the $V_0$ model obtained with a coarse frequency sampling during FWI6 is very close to the one of FWI5. This suggests that the fine frequency sampling of the 5.2–7 Hz frequency group during tests FWI3–FWI5 is not mandatory.

To identify more precisely differences between the $V_0$ models of the six tests, we assess the $V_0$ models against a well log in Fig. 9. Interestingly, the $V_0$ models built during the attenuation-free FWI1 shows a slight deepening of the reservoir reflector at 2.5 km depth that might result because dispersion effects generated by attenuation have not been taken into account during seismic modelling and inversion. Another striking feature is the simpler velocity structure between the reservoir at 2.5 km depth and the base cretaceous reflector at 3.5 km depth of the FWI4 test. Hypothetically, this might indicate a different weight of extrinsic versus intrinsic attenuation, a more significant $Q$ update at the reservoir depth during FWI4 translating into a less heterogeneous structure. With these two differences, the velocity profiles of the six tests remain quite consistent providing a first evidence of limited $V_0-\rho$ and $V_0-Q$ cross-talks.

### $Q$ models

The final FWI $Q$ models of tests FWI3–FWI6 are shown in Figs 10 (depth slices) and 11 (vertical slices). The $Q$ results are shown as perturbation models with respect to the initial homogeneous
<figure><img src="image_17.png" /><figcaption>**Figure 13.** FWI perturbation models ($\rho$). Vertical sections (a) across gas cloud ($x = 5.575$ km), (b) near the periphery of the gas cloud ($x = 6.25$ km) and (c) away from the gas cloud ($x = 9.5$ km). In (a)–(c), each panel is labelled by the test number (Table 1).</figcaption></figure>

model ($Q = 200$). The $Q$ parameter absorbed a significant part of the acquisition footprint due to its intrinsic link with amplitudes. This acquisition footprint is however less strong in the $Q$ model obtained with FWI4 due to frequency grouping and small damping in the pre-conditioner. For this test, the geometry of the gas cloud as well as small-scale low-$Q$ anomalies radiating from the gas cloud previously visible in the $V_0$ models are well delineated in Fig. 10(c). For all of the tests, we show a good correlation between low-velocity zones in the shallow sediments and the gas cloud and low-$Q$ anomalies, which indeed does not allow us to conclude that ($V_0, Q$) cross-talks did not occur. We also check that the attenuation imprint decreases as we move southward from the gas cloud (from

panels a to c in Fig. 11). Different pre-conditionings of the $Q$ update used for tests FWI4 and FWI5 manifest as expected by a different distribution of the $Q$ perturbations with depth. The $Q$ model inferred from FWI5 shows more focused and lower values of $Q$ in the shallow sediments (with values as low as 60), while the $Q$ model of test FWI4 shows milder $Q$ perturbations ($Q$ of the order of 110 in the shallow sediments) that are more uniformly distributed in depth (due to smaller damping $\beta$ of the pre-conditioner). Interestingly, we also show significant negative $Q$ perturbations at the reservoir depths below the gas cloud (Figs 11a and b).
<figure><img src="image_18.png" /><figcaption>**Figure 14.** Validation of $\rho$ against well log (see the position of the log in Fig. 2, grey star). (a) Comparison between density well log (black) and corresponding profiles from the initial (blue) and final (red) FWI density models. Each panel is labelled by the test number (Table 1). (b) Vertical section of the final (top) $V_0$ and (bottom) $\rho$ perturbation models (FWI4) across the well log ($x = 6.23$ km). The arrows point two low $V_0/\rho$ anomalies that are well matched on the density log.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 15.** Impedance model inferred from the FWI4 $V_0$ and $\rho$ models. As for the $V_0$ models, the reflectivity of the impedance model is superimposed in transparency to delineate the main discontinuities. Vertical sections at (a) $x = 5.75$ km, (b) $x = 6.25$ km and (c) $x = 9.5$ km.</figcaption></figure>

### $\rho$ models

As $Q$, the final FWI $\rho$ models built during tests FWI3–FWI6 are shown in the form of model perturbations with respect to the initial model [Figs 12 (depth slices) and 13 (vertical slices)]. We had to strongly damp the density update at reservoir depths where anisotropy is the most significant due to instabilities. Aggressive weighting of the density update was also performed in Goh *et al.* (2016, page SU20). Therefore, density values may be only reliable at shallow and intermediate depths where they can however have a significant imprint on the wavefield amplitudes in particular at the sea bed. As $Q$, the density absorbs a significant part of the acquisition footprint. We however delineate fairly well the geometry of the glacial sand channel deposits in the $\rho$ slices of tests FWI4–FWI6 (Fig. 12a). The vertical sections of the $\rho$ perturbations (Fig. 13) show the main reflectors (500 m depth reflector, top of the gas cloud, top hard chalk) as well as near-surface small-scale vertical features (Fig. 13b, $y = 11$ km), formally interpreted in Sirgue *et al.* (2010) and Barkved *et al.* (2010). This migration behaviour is consistent with the radiation pattern of $\rho$ in a ($V_0$, $\rho$)-type parametrization. The $\rho$ section inferred from FWI3 shows sharper perturbations relative to those inferred from the multifrequency inversions. Again, this can result either because mono-frequency inversions have more freedom to over fit the data or different parameter cross-talk footprint. It is also worth noting how the different $Q$ scaling used during the FWI4 and FWI5 tests manifest in the $\rho$ update: relative to FWI4,
<figure><img src="image_20.png" /><figcaption>**Figure 16.** (a–c) Mean (top panels) and standard deviation (bottom panels) of the multiparameter models (tests FWI2 to FWI7). Left and right panels show sections across the gas cloud and its periphery, respectively. (a) *V*₀ model. (b) *Q* model. (c) *ρ* model. The reflectivity of the mean *V*₀ model is superimposed on the mean *V*₀ model in (a) and the mean *Q* perturbation model in (b).</figcaption></figure>

perturbations built by FWI5 are more focused in the shallow sub-surface.

The density models built during FWI3–FWI6 are further validated against a well log (Fig. 14). The well log intersects two shallow low-velocity/low-density anomalies at around 500 and 700 m depths (Fig. 14b, black arrows). These two anomalies are fairly well reconstructed by the four FWI tests. Surprisingly, FWI3 seems even to show the best agreement with the well log.

A last indirect validation of the ρ models can be performed by building impedance models by nonlinear recombination (i.e. multiplication) of the V₀ and ρ FWI models. If ρ was a proxy parameter to absorb elastic effects (Borisov *et al.* 2014), it is likely that this recombination would produce unfocused impedance models. The impedance models inferred from the FWI4 V₀ and ρ models are

shown in Fig. 15 and show fairly well focused structures consistent with those shown in the *V*₀ models (Fig. 8). Structural differences between the *V*₀ and impedance models are mainly visible near the sea bed where density might have been useful to absorb elastic effects and match guided waves.

### Synthesis

As a synthesis, we generate the mean and standard deviation of the (*V*₀, *ρ*, *Q*) FWI models of tests FWI3 to FWI7 where FWI7 is the coarse-frequency analogous of test FWI4 (Table 1). Our aim is not, of course, to draw statistical inferences from this showing but merely to assess the dispersion with which the three parameters are reconstructed depending on the experimental setup. The standard
<figure><img src="image_21.png" /><figcaption>**Figure 17.** Frequency-domain data fit. 7 Hz common-receiver gather for the receiver near the gas cloud (grey star in Fig. 2). (a) Recorded gather. (b) Synthetic gather computed in the FWI model after the 7 Hz inversion. From left to right, each panel is labelled by the test number (Table 1). (c) Residuals. The dashed line represents the profile along which a direct comparison between recorded and modelled data is shown in Fig. 18(a).</figcaption></figure>

deviation of $V_0$ is negligible in the central part of the vertical section, nearby the gas cloud and reach a maximum value of 100 m s$^{-1}$ at reservoir depths near the ends of the section (Fig. 16a). In contrast, the standard deviation of $Q$ is more significant and can reach values as high as 80 in the shallow sediments and at the reservoir depths (Fig. 16b). Since the high wavenumbers of the $\rho$ and $Q$ updates are probably coupled to some extent, the significant standard deviations of $Q$ translate into significant standard deviations of $\rho$ along the three main reflectors at 500, 1500 and 2500 m depths (Fig. 16c). The mean $V_0$ sections across the gas cloud and its periphery show fairly well focused structures (Fig. 16a). We superimpose in transparency the reflectivity of the mean FWI $V_0$ model on the mean $Q$ perturbations in Fig. 16(b) to correlate the attenuation distribution with the geology.

The fact that significant standard deviations of the second-order parameters $Q$ and $\rho$ do not translate into significant standard deviation of the $V_0$ updates prompts us to conclude that the FWI $V_0$ models have not been significantly hampered by parameter cross-talks, although the overlapping radiation patterns of these three parameters would have suggested the contrary. The obvious reason is that, at such low frequencies, the $V_0$ updates are to the first-order tied by the need to fit the kinematic properties of the wavefield rather than their amplitudes. The significant standard deviations of the $Q$ and $\rho$ updates reflect the ill-posedness with which such second-order parameters, which are mainly sensitive to amplitudes, are reconstructed from narrow low frequency bands. Moreover, the significant standard deviation of $Q$ and $\rho$ updates are mainly related to amplitudes variations. However, each $Q$ and $\rho$ updates show consistent structural trend (with more or less smearing in depth) as

this trend is mainly controlled by the kinematic accuracy of the $V_0$ model.

## Quality control

In order to gain additional insights on the reliability of the multiparameter subsurface models, we compare the recorded data with the synthetic ones computed in each FWI models both in the frequency and time domains.

### Frequency-domain data fit

Operto et al. (2015) hypothesized a significant footprint of attenuation in this dataset after noting that the amplitudes of the modelled 7 Hz monochromatic gathers were overestimated compared to the recorded ones at long offsets along the parallel profile intersecting the gas cloud. We check now this hypothesis by comparing the recorded and modelled amplitude for the different FWI tests (Figs 17 and 18a).

Fig. 17 shows the recorded 7 Hz common-receiver gather 1 (Fig. 2, white circle) as well as the modelled gathers computed in the FWI models of tests FWI1, FWI3 and FWI4 obtained after the 7 Hz frequency group inversion (Table 1). The dashed line represents the profile that intersects the gas cloud along which recorded and modelled amplitudes are directly compared in Fig. 18(a). We show that the data fit achieved by FWI3 designed as successive mono-frequency inversions is superior than the one of FWI4 designed as successive multifrequency inversions. Again, this results because the later inversion processes more redundant data. We will
<figure><img src="image_22.png" /><figcaption>**Figure 18.** Frequency-domain data fit. (a and b) Direct comparison between recorded (black) and modelled (grey) monochromatic data extracted from the 7Hz common-receiver gathers associated with receivers 1 and 2 in Fig. 2. The recorded and modelled data are compared along a parallel line cross-cutting the receiver position. (a) The parallel line (x = 5.5 km) across receiver 1 cross-cuts the gas cloud (Fig. 17). (b) The parallel line (x = 6.65 km) across receiver 2 is away from the gas cloud [see the monochromatic gather in Operto *et al.* (2015, their fig. 15)]. In (a) and (b), each panel is labelled by the test number (Table 1). In (a), the shaded area highlights recorded data, whose amplitudes might have been damped by intrinsic attenuation in the gas cloud. Such attenuation is not shown away from the gas cloud in (b). A linear gain with offset was applied to amplitudes for geometrical spreading correction.</figcaption></figure>

discuss in the following section whether the same conclusion holds in time-domain synthetics that jointly include the effects of all of the frequencies.

The direct comparison of the amplitudes along the parallel profile intersecting the gas cloud confirms that the amplitude fit at long offsets improves from tests FWI1 to FWI3–FWI5 as attenuation effects are more accurately taken into account for (Fig. 18a). The amplitude fit is similar for FWI3 and FWI4 and seems slightly superior than the one achieved by FWI5 and FWI6. This is consistent with the fact that more significant attenuation updates have been performed in the deep part during FWI3 and FWI4 relative to FWI5 and FWI6.

We show also the direct comparison between recorded and modelled amplitudes for the receiver 2 along the parallel profile intersecting it (Fig. 2, grey circle) to illustrate that the amplitude fit is much more consistent from one test to the next when the wavefield propagates outside the gas cloud's area of influence.

### *Time-domain data fit*

We now show the data fit in the time domain to assess more precisely which parts of the wavefields have been matched during the

different FWI tests. We compute time-domain seismograms with the same frequency-domain forward modelling engine as the one used to perform FWI to prevent any bias in this assessment. For a simulation length of 8 s and a cut-off frequency of 10 Hz, this required to perform 81 frequency-domain modelling (LU decomposition + solve phases) on the 35 m grid before inverse Fourier transform. This number of modelling could have been, however, decreased by a factor two if a reduced time scale implemented with complex-valued frequencies would have been used as illustrated in Fig. 3, bottom panels, with the additional difficulty to control wraparound effects (Mallick & Frazer 1987).

As in Operto *et al.* (2015), we first compute Green functions with Dirac delta source signatures before estimating the source signature with the approach of Pratt (1999). These source signatures computed in the initial model and the final FWI models of tests FWI1 and FWI4 are shown in Figs 19 and 20 for the two receivers pointed in Fig. 2. These source signatures are consistent from one test to the next because their estimation is mainly driven by the high-amplitude early arrivals at short offsets when no weighting is applied to the data. In this case, the source estimation is poorly sensitive to the accuracy of the subsurface model (Operto *et al.* 2015). Comparison between recorded and modelled seismograms
<figure><img src="image_23.png" /><figcaption>**Figure 19.** Time-domain data fit. Receiver gather 1 (Fig. 2, white circle). In the top, middle and bottom rows the seismograms are computed in the initial model, in the final FWI model of test FWI1 and in the final model of the test FWI4, respectively. In each row, the first and third panels show the real data, while the second panel shows the modelled data with a mirror representation. The inset in the second panel shows the source wavelet estimated by matching the modelled impulsive seismograms with the real ones. The right panel shows a direct comparison between the real data (red/white/blue colour scale) and the modelled data (black/grey/white scale). The two sets of seismograms are in phase if the black/white pixels on the synthetic section overlay the red/blue ones of the recorded one.</figcaption></figure>

shows that, compared to the initial model, FWI contributes to absorb some time shifts of the diving waves at long offsets and match the short-spread reflections (shaded ellipse and triangle areas in Figs 19 and 20, respectively). Moreover, comparing seismograms computed in the FWI1 and FWI4 models confirms that FWI4 improves the match of the dispersive wavefield at long offsets for receiver 1 where attenuation effects are expected to be significant, while the fit achieved by the two FWI tests look more similar for receiver 2 where attenuation effects are expected to be smaller (Figs 19 and 20, ellipse).

For a more quantitative assessment of the data fit, Figs 21 and 22 show a direct comparison between recorded seismograms and modelled seismograms computed in the initial model and FWI models of tests FWI1 to FWI6 for receivers 1 and 2. We focus on three different part of the wavefields: the early arrivals in box 1, the short-spread reflections in box 2 and the dispersive guided waves obscuring the deep post-critical reflections in box 3. Comparing the

seismograms computed in the starting model and the FWI1 model shows how FWI aligns the two sets of seismograms in box 1 to a large extent and generate the short-spread reflections. Moving from tests FWI1 to FWI3 shows that the fit of the amplitude versus offset variations and the phase of the early arrivals are progressively improved. However, none of these mono-frequency inversions managed to fit the dispersive part of the wavefield (box3). Only FWI4 and FWI5 tests, when multiple frequencies are jointly inverted, allow to match the dispersive part of the wavefield. FWI4 performs slightly better than FWI5 to match the deep short-spread reflections probably due to a better accounting for attenuation at the reservoir depths. In the box 3, the overall match achieved by FWI4 and FWI5 is satisfactory. Depending which part of the waveforms is analysed, one or the other test fits better the data. Comparing the spectral amplitudes of the seismograms computed in the FWI model of tests FWI3, FWI4 and FWI5 with those of the real data at different offsets suggest that the FWI5 test slightly outperforms FWI4, while
<figure><img src="image_24.png" /><figcaption>**Figure 20.** Time-domain data fit. Same as Fig. 19 for receiver 2 (Fig. 2, grey circle).</figcaption></figure>

<table><caption>**Table 1.** FWI experimental setup. *h(m)* and *Optim*: grid interval and optimization algorithm used during each multiscale step. *pSTD*: pre-conditioned steepest-descent. The superscript * means that $V_0$, $\rho$ and $Q$ are jointly updated. Otherwise, only $V_0$ is updated. $Q_b$: $Q$ in the initial/background homogeneous $Q$ model. $Q_0$, $\beta_Q$: $Q$ constant scaling factor and damping factor of the $Q$ gradient pre-conditioner. Compared to tests 1–3, we skip the 9 Hz frequency during tests 4 and 5 because no enough iterations can be performed in 24 hr when three frequencies are jointly inverted on the 35 m grid (the operating rule of your cluster sets a maximum wall time of 24 hr). The results of test FWI7 are not shown here because they lead to similar conclusions than test FWI6 in terms of frequency sampling. The results of this test are just used to augment the population of FWI models used to compute a mean and standard deviation models in Fig. 16.</caption><thead><tr><th rowspan="3">Tests</th><th colspan="3">70</th><th colspan="3">50</th><th colspan="3">35</th><th rowspan="3">Q<sub>b</sub></th><th rowspan="3">Q<sub>0</sub>, β<sub>Q</sub></th></tr><tr><th colspan="3">pSTD</th><th colspan="3">pL-BFGS</th><th colspan="3"></th></tr><tr><th colspan="10">Frequency groups (Hz)</th></tr></thead><tbody><tr><td>FWI 1</td><td>3.5</td><td>4</td><td>4.5</td><td>5</td><td>5.2</td><td>5.8</td><td>6.4</td><td>7</td><td>8</td><td>9</td><td>10</td><td>1000</td><td>—</td></tr><tr><td>FWI 2</td><td>3.5</td><td>4</td><td>4.5</td><td>5</td><td>5.2</td><td>5.8</td><td>6.4</td><td>7</td><td>8</td><td>9</td><td>10</td><td>200</td><td>—</td></tr><tr><td>FWI 3</td><td>3.5</td><td>4</td><td>4.5</td><td>5*</td><td>5.2*</td><td>5.8*</td><td>6.4*</td><td>7*</td><td>8*</td><td>9*</td><td>10*</td><td>200</td><td>200, 2e<sup>-3</sup></td></tr><tr><td>FWI 4</td><td>3.5</td><td>4</td><td>4.5</td><td>5*</td><td></td><td>(5.2, 5.8, 6.4, 7)*</td><td></td><td></td><td>(8, 10)*</td><td></td><td></td><td>200</td><td>200, 2e<sup>-3</sup></td></tr><tr><td>FWI 5</td><td>3.5</td><td>4</td><td>4.5</td><td>5*</td><td></td><td>(5.2, 5.8, 6.4, 7)*</td><td></td><td></td><td>(8, 10)*</td><td></td><td></td><td>200</td><td>400, 4e<sup>-3</sup></td></tr><tr><td>FWI 6</td><td>3.5</td><td>4</td><td>4.5</td><td>5*</td><td></td><td>(5, 7)*</td><td></td><td></td><td>(7, 10)*</td><td></td><td></td><td>200</td><td>400, 4e<sup>-3</sup></td></tr><tr><td>FWI 7</td><td>3.5</td><td>4</td><td>4.5</td><td>5*</td><td></td><td>(5, 7)*</td><td></td><td></td><td>(7, 10)*</td><td></td><td></td><td>200</td><td>200, 2e<sup>-3</sup></td></tr></tbody></table>

significant improvement relative to FWI3 is confirmed (Fig. 23). The improved spectral match achieved by FWI5 may result because the high-amplitude guided waves have a dominant imprint in the

wavefield and the more significant attenuation in the shallow part of the FWI5 model has been useful to fit these guided waves.
<figure><img src="image_25.png" /><figcaption>**Figure 21.** Time-domain data fit. Direct comparison between recorded (black) and modelled (red) seismograms, these later being computed in the initial, FWI1 and FWI2 models. A reduced timescale is used as in Fig. 3. True amplitudes are shown with a gain with offset to allow for the data fit assessment at long offsets. The right and left panels correspond to receivers 1 and 2, respectively (see also Figs 19 and 20).</figcaption></figure>

DISCUSSION

We have presented a real data case study of multiparameter
frequency-domain FWI for the reconstruction of the vertical wave
speed, density and *Q* in the 3.5–10 Hz frequency band.

Does efficient multiparameter frequency-domain FWI
reliable?

From this study, we conclude that efficient multiparameter
frequency-domain FWI with a limited number of discrete frequen-
cies shows reliable and consistent velocity updates provided that
frequency grouping is used. We have shown how frequency group-
ing is helpful in frequency-domain FWI to match the dispersive
part of the wavefield associated with guided waves propagating in
the near surface. The match of these shallow events has probably
removed artefacts in the shallow part of the FWI models, which
prevents error accumulation at depth. The frequency band covered
by the frequency groups in this study is of the order of 2 Hz, which
is consistent with the experimental setup described in Stopin *et al.*
(2016). A fine sampling of these frequency groups as that used in
time-domain FWI seems not mandatory as shown by the results of

<figure><img src="image_26.png" /><figcaption>**Figure 22.** Following of Fig. 21 for tests FWI3 to FWI5.</figcaption></figure>

test FWI6, although it is likely that increasing the number of fre-
quencies in the group will make the multiparameter reconstruction
better posed. Clearly a trade-off between computational efficiency
and reliability of the multiparameter reconstruction should be found
to define the number of frequencies involved in the inversion ac-
cording to the available computational resources.

The V₀ models inferred from multiparameter reconstruction are
quite consistent with former models built by mono-parameter FWI
using a homogeneous Q model in a passive way (Operto et al.
2015). This strongly suggests that the V₀ updates inferred from
multiparameter reconstruction have not been significantly contam-
inated by ρ and Q cross-talks artefacts. A notable feature is that the
ρ and Q updates absorb a significant part of the acquisition foot-
print hence cleaning the V₀ updates from this imprint. This allows
us to use a regularization-free I-BFGS optimization, hence taking
advantage of an improved convergence rate without boosting high
frequency noise in the V₀ FWI models. The fact that the V₀ was
not contaminated by cross-talks artefacts shows that the reliability
of multiparameter FWI does not only rely on the radiation patterns
(since significant cross-talks would have been expected from the
strongly overlapping radiation patterns of V₀, ρ and Q). At low fre-
quencies (f < 10 Hz), FWI is firstly driven by the need to match
phase through the wave-speed update. If this update is not hampered
by cycle skipping artefacts, it is likely that FWI can update during
late FWI iterations second-order parameters such as attenuation and
density which are more closely related to amplitudes. However, the
<figure><img src="image_27.png" /><figcaption>**Figure 23.** Spectral analysis. Comparison between spectral amplitudes of recorded and modelled common-receiver gather 1 (Fig. 2, white circle). Modelling is performed in the final FWI model of test (a) FWI3, (b) FWI4 and (c) FWI5. The negative offsets for which the spectrum are computed are labelled above each panel (decreasing offsets from left to right according to the trace sorting in Fig. 3a).</figcaption></figure>

update is not mandatory to perform reliable update of wave speed as shown by our former mono-parameter application (Operto et al. 2015) and the sensitivity analysis of Prieux et al. (2013) and Kurzmann et al. (2013). The fact that smooth background models of second-order dynamic parameters are enough to perform reliable updates of the first-order kinematic parameters does not mean however that these second-order parameters belong to the null space of the FWI.

A more difficult question to answer is related to the reliability of the density and *Q* updates since their imprint in the wavefields remain small. We have first shown the positive correlation between low velocities associated with shallow soft sediments and gas and low-*Q* anomalies. We have shown how the attenuation updates allow to improve the fit of the amplitude versus offset variations of the diving waves. By performing frequency grouping, we nicely match the dispersive part of the wavefield. Since subsurface models inferred from multifrequency inversions tend to be less shingling than those inferred from mono-frequency inversions, this improved fit of the dispersive wavefield likely results from an improved estimation of the intrinsic attenuation rather than from artificial extrinsic attenuation.

Density is probably more challenging to update than attenuation as short-spread reflections are mostly sensitive to short-scale density perturbations. We had to strongly damp the density update at the reservoir depth to prevent instabilities where anisotropy is significant. These instabilities might have been partly generated by the basic seven-point stencil that we currently use to build the radiation pattern matrices. It is worth mentioning that the elliptic and anelliptic component of the anisotropy have been separated in the wave operator $A_h$ from which radiation pattern are computed, eq. (5). Both components impact upon the radiation pattern of the different parameter classes. Our practical experience is that the anelliptic component contributes to make the FWI updates more unstable. This issue requires further investigations, a pragmatical approach consisting in computing the radiation patterns assuming an elliptic anisotropy only.

To assess the relevance of the density update, we validate the FWI density model against a well log. From this comparison, we conclude that we successfully recover two low density anomalies that are

well correlated with two low-velocity anomalies at around 500 m depth. Moreover, we build a *posteriori* an impedance model by multiplying the vertical wave speed and the density. This nonlinear recombination didn't hamper the focusing of the structure in the impedance model. This would have been the case if density would have been used as a proxy to absorb elastic effects in most part of the subsurface. Structural differences are however shown between the shallow part of the $V_0$ and impedance models. This suggests that density might have been used as a proxy to absorb elastic effects near the sea bed although the soft sedimentary environment should prevent a significant imprint of shear waves.

From the optimization viewpoint, some improvements are possible. Instead of using a purely diagonal pre-conditioner, we could have used a block diagonal pre-conditioner where the diagonal coefficients of the off-diagonal blocks partially account for parameter cross-talk (Korta et al. 2013). This pre-conditioner can be coupled with a truncated Newton optimization whose implementation should be tractable in 3-D frequency domain FWI (Métivier et al. 2015). However, in order to avoid the re-computation of the LU factors associated with each frequencies of a group during the iterative resolution of the Newton system, each frequency of a group should be processed in parallel rather than sequentially. This implies that the computational resources should be augmented by a factor equal to the number of frequencies in the group. Decimation or blending of the reciprocal sources should also be viewed to mitigate the computational burden of the substitution steps during the resolution of the Newton system to make the truncated Newton approach tractable (Castellanos et al. 2015).

## CONCLUSION

We have shown with a real-data case study the feasibility of multiparameter frequency-domain FWI of stationary-receiver wideazimuth OBC data for the reconstruction of the vertical wave speed, density and *Q* in the 3.5–10 Hz frequency band. Although the simultaneous inversion of multiple frequencies is necessary to fit the dispersive part of the wavefield during each multiscale step, the number of frequencies involved in the inversion can remain
limited hence preserving the computational efficiency of the fre-
quency domain formulation in a stationary-receiver setting. Since
FWI is driven to the first order by the need to fit the phase of the
wavefield, the wave-speed update has not been significantly con-
taminated by density and attenuation cross-talks. On the contrary,
these amplitude-related parameters have been useful to absorb ac-
quisition footprint and hence remove this footprint from the velocity
update. The relevance of the attenuation and density updates is more
difficult to assess as they show a significant dispersion in terms of
amplitudes when the tuning of the inversion is modified. However,
the spatial distribution of the density and Q anomalies show a good
correlation with the geology. In particular, low-Q anomalies corre-
late well with low-velocity anomalies associated with soft sediments
and gas. Density seems more challenging to update. However, com-
parison with a well log as well as the reliable a posteriori building of
the impedance from the velocity and density FWI models suggests
that this parameter has been reliably reconstructed in the shallow
part of the subsurface. Our future work aims to move from I-BFGS
to truncated Newton optimization to improve the accounting of the
Hessian operator during the early FWI iterations.

ACKNOWLEDGEMENTS

Seismic modelling was performed with the MUMPS package,
available on http://graal.ens-lyon.fr/MUMPS/index.html and http:
//mumps.enseeiht.fr. We thank the developers of the MUMPS re-
search team (P. Amestoy, A. Buttari, T. Mary, J.-Y. L'Excellent, C.
Puglisi) for useful discussions and recommendations concerning the
interfacing of MUMPS in our imaging code. This study was granted
access to the HPC resources of SIGAMM (http://crimson.oca.eu)
and CIMENT computer centres and CINES/IDRIS under the al-
location 046091 made by GENCI. We thank R. Brossier (IS-
Terre), L. Combe (Geoazur), L. Métivier (ISTerre/LJK), A. Ri-
bodetti (Geoazur) and J. Virieux (ISTerre) for useful discussions.
We also thank two anonymous reviewers for their useful comments
and the editor R.-E. Plessix. We thank AkerBP and their Valhall
partner Pandion Energy for allowing access to the OBC data set,
the well-log velocities and initial FWI models and for permission
to publish this work.

REFERENCES

Aki, K. & Richards, PG., 2002. *Quantitative Seismology, Theory and Methods*, 2nd edn, University Science Books.

Amestoy, P., Brossier, R., Buttari, A., L'Excellent, J.-Y., Mary, T., Métivier, L., Miniussi, A. & Operto, S., 2016. Fast 3D frequency-domain full waveform inversion with a parallel block low-rank multifrontal direct solver: application to OBC data from the North Sea, *Geophysics*, **81**, R363–R383.

Askan, A., Akcelik, V., Bielak, J. & Ghattas, O., 2007. Full waveform inversion for seismic velocity and anelastic losses in heterogeneous structures, *Bull. seism. Soc. Am.*, **97**, 1990–2008.

Bai, J. & Yingst, D., 2014. Simultaneous inversion of velocity and den-
sity in time-domain full waveform inversion, in *SEG Technical Program
Expanded Abstracts*, pp. 922–927, SEG.

Barkved, O., Albertin, U., Heavey, P., Kommedal, J., van Gestel, J., Synnove, R., Pettersen, H. & Kent, C., 2010. Business impact of full waveform inversion at Valhall, in *91 Annual SEG Meeting and Exposition, Expanded Abstracts*, Denver, pp. 925–929, Society of Exploration Geophysics.

Boiero, D., Wiarda, E. & Vermeer, P., 2013. Surface- and guided-wave inversion for near-surface modeling in land and shallow marine seismic data, *Leading Edge*, **32**, 638–646.

Boiero, D., Wiarda, E., West, L.R. & Vermeer, P., 2014. Near-surface mod-
elling in shallow marine environments using surface and guided waves,
in *76th EAGE Conference and Exhibition*, Amsterdam, We D201 10.

Borisov, D., Stopin, A. & Plessix, R.-E., 2014, Acoustic pseudo-density full waveform inversion in the presence of hard thin beds, *Presented at the 76th Annual EAGE Meeting*, Amsterdam.

Castellanos, C., Métivier, L., Operto, S., Brossier, R. & Virieux, J., 2015. Fast full waveform inversion with source encoding and second-order optimization methods, *Geophys. J. Int.*, **200**(2), 720–744.

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

Futterman, W., 1962. Dispersive body waves, *J. geophys. Res.*, **67**, 5279–5291.

Gholami, Y., Brossier, R., Operto, S., Ribodetti, A. & Virieux, J., 2013. Which parametrization is suitable for acoustic VTI full waveform inversion? - Part 1: sensitivity and trade-off analysis, *Geophysics*, **78**, R81–R105.

Goh, V., Halleland, K., Plessix, R.-E. & Stopin, A., 2016, Application of mul-
tiparameter full-waveform inversion in Central Luconia Basin, Sarawak,
*Interpretation*, **4**, SU17–SU24.

Hak, B. & Mulder, W.A., 2011. Seismic attenuation imaging with causality, *Geophys. J. Int.*, **184**, 439–451.

Haller, N., Flateboe, R., Twallin, C., Dahl-Eriksen, V., Heavey, P., Kjos, E., Milne, R. & Rietveld, W., 2016, Valhall case study - value of seismic technology for reducing risks in a reactive overburden, *Presented at the Expanded Abstracts, 78th Annual EAGE Meeting*, Vienna.

Hicks, G.J. & Pratt, R.G., 2001. Reflection waveform inversion using local descent methods: estimating attenuation and velocity over a gas-sand deposit, *Geophysics*, **66**, 598–612.

Jeong, W., Lee, H. & Min, D., 2012. Full waveform inversion strategy for
density in the frequency domain, *Geophys. J. Int.*, **188**, 1221–1242.

Kamei, R. & Pratt, R.G., 2013. Inversion strategies for visco-acoustic wave-
form inversion, *Geophys. J. Int.*, **194**, 859–894.

Kolsky, H., 1956. The propagation of stress pulses in viscoelastic solids,
*Phil. Mag.*, **1**, 693–710.

Korta, N., Fichtner, A. & Sallarès, V., 2013. Block-diagonal approximate Hessian for preconditioning in full waveform inversion, in *Presented at the Expanged Abstracts, 75th EAGE Conference & Exhibition incorporating SPE EUROPEC 2013*.

Kurzmann, A., Przebindowska, A., Kohn, D. & Bohlen, T., 2013. Acoustic full waveform tomography in the presence of attenuation: a sensitivity analysis, *Geophys. J. Int.*, **195**(2), 985–1000.

Lambaré, G., Operto, S., Podvin, P., Thierry, P. & Noble, M., 2003. 3-D
ray+Born migration/inversion—part 1: theory, *Geophysics*, **68**, 1348–
1356.

Malinowski, M., Operto, S. & Ribodetti, A., 2011. High-resolution seismic attenuation imaging from wide-aperture onshore data by visco-acoustic frequency-domain full waveform inversion, *Geophys. J. Int.*, **186**, 1179–1204.

Mallick, S. & Frazer, L.N., 1987. Practical aspects of reflectivity modeling,
*Geophysics*, **52**, 1355–1364.

Marfurt, K., 1984. Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations, *Geophysics*, **49**, 533–549.

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

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 Techni-
cal Program Expanded Abstracts 2015*, 1198–1203.

Miller, D., Oristaglio, M. & Beylkin, G., 1987. A new slant on seismic imaging: Migration and integral geometry, *Geophysics*, **52**, 943–964.

Mulder, W.A. & Hak, B., 2009. An ambiguity in attenuation scattering imaging, *Geophys. J. Int.*, **178**, 1614–1624.

Operto, S., Virieux, J., Amestoy, P., L'ÉExcellent, J.-Y., Giraud, L. & Ben Hadj Ali, H., 2007. 3D finite-difference frequency-domain modeling of
visco-acoustic wave propagation using a massively parallel direct solver: A feasibility study, Geophysics, **72**, SM195–SM211.

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**, 1040–1054.

Operto, S., Brossier, R., Combe, L., Métivier, L., Ribodetti, A. & Virieux, J., 2014. Computationally-efficient three-dimensional visco-acoustic finite-difference frequency-domain seismic modeling in vertical transversely isotropic media with sparse direct solver, *Geophysics*, **79**(5), T257–T275.

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**, 1362–1391.

Podgornova, O., Leaney, S. & Liang, L., 2015, Analysis of resolution limits of VTI anisotropy with full waveform inversion, in *SEG Technical Program Expanded Abstracts*, 1188–1192.

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.

Prieux, V., Brossier, R., Operto, S. & Virieux, J., 2013. Multiparameter full waveform inversion of multicomponent OBC data from Valhall. Part 1: imaging compressional wavespeed, density and attenuation, *Geophys. J. Int.*, **194**, 1640–1664.

Przebindowska, A., Kurzmann, A., Köhn, D. & Bohlen, T., 2012, The role of density in acoustic full waveform inversion of marine reflection seismics, in *Presented at the 74th EAGE Conference and Exhibition incorporating EUROPEC 2012*.

Qin, B. & Lambaré, G., 2016, Joint inversion of velocity and density in preserved-amplitude full waveform inversion, in *SEG Technical Program Expanded Abstracts 2016*, 1325–1330.

Ribodetti, A., Operto, S., Virieux, J., Lambaré, G., Valéro, H.-P. & Gibert, D., 2000. Asymptotic viscoacoustic diffraction tomography of ultrasonic laboratory data: a tool for rock properties analysis, *Geophys. J. Int.*, **140**, 324–340.

Roth, M., Holliger, K. & Green, A.G., 1998. Guided waves in near-surface

seismic surveys, Geophys. Res. Lett., **25**, 1071–1074.

Shin, C., Jang, S. & Min, D.J., 2001. Improved amplitude preservation for prestack depth migration by inverse scattering theory, *Geophys. Prospect.*, **49**, 592–606.

Shtivelman, V., 2004. Estimating shear wave velocities below the sea bed using surface waves, *Near Surf. Geophys.*, **2**, 241–247.

Sirgue, L. & Pratt, R.G., 2004. Efficient waveform inversion and imaging: a strategy for selecting temporal frequencies, *Geophysics*, **69**, 231–248.

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.

Stopin, A., Plessix, R.-E., Kuehl, H., Goh, V. & Overgaag, K., 2016, Application of visco-acoustic full waveform inversion for gas cloud imaging and velocity model building, in *Presented at the 78th EAGE Conference and Exhibition 2016*, EAGE.

Takam Takougang, E.M. & Calvert, A.J., 2013. Seismic waveform tomography across the Seattle fault zone in Puget Sound: resolution analysis and effectiveness of visco-acoustic inversion of viscoelastic data, *Geophys. J. Int.*, **193**, 763.

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

Vigh, D., Jiao, K., Watts, D. & Sun, D., 2014. Elastic full-waveform inversion application using multicomponent measurements of seismic data collection, *Geophysics*, **79**, R63–R77.

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

Wiarda, E., Shaw, S.A., Boiero, D., Gundersen, A. & West, L.R., 2014, Integrating FWI with surface-wave inversion to enhance near-surface modelling in a shallow water setting in Eilfisk, in *76th EAGE Conference and Exhibition*, Amsterdam, Th E106 03.

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

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

Yang, J., Liu, Y. & Dong, L., 2016b. Simultaneous estimation of velocity and density in acoustic multiparameter full-waveform inversion using an improved scattering-integral approach, *Geophysics*, **81**, R399–R415.

Yang, P., Brossier, R., Métivier, L. & Virieux, J., 2016c. A review on the systematic formulation of 3D multiparameter full waveform inversion in viscoelastic medium, *Geophys. J. Int.*, **207**, 129–149.