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, Denver, pp. 925–929, Society of Exploration Geophysics. Boiero, D., Wiarda, E. & Vermeer, P., 2013. Surfaceand 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 modelling 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., Ribodetti, 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–


I N T RO D U C T I O N
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 mesoand 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 where ω is the angular frequency, c the local wave speed, θ the scattering angle and φ the dip angle (Miller et al. 1987;Wu & Toksöz  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 crosstalks (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 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 lowvelocity 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).
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 (ρ) 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 δ and 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 generatesE106 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 (ω, θ) 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.
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 θ -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 , ρ, Q, δ, ) parametrization.Our motivation beyond  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 −1 .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.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 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 , ρ 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 ρ and Q as optimization parameters rather than passive parameters.We show how the acquisition footprint tends to be focused in the ρ and Q updates, hence cleaning the V 0 model from their imprint.We also show the dispersion with which Q and ρ are updated depending on the inversion tuning.However, this dispersion is not passed on V 0 , which prompts us to conclude that V 0 -ρ 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 , ρ 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.

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 finitedifferent 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 where p h , p v and p are the so-called horizontal, vertical and physical pressure vectors.The operators A h and A v discretized in matrices A h and A v are given by where κ 0 = ρV 2 0 , ρ is density, V 0 is the vertical wave speed, ω is the angular frequency, δ and are Thomsen's parameters (Thomsen 1986).Differential operators X , Y and Z are given by ∂ x b∂ x , ∂ ỹ b∂ ỹ and ∂ z b∂ z , respectively, where b = 1/ρ is buoyancy and ( x, ỹ, 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 s and 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 where ω 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 : where m gathers the subsurface model parameters of different kinds and d s,ω = d s,ω (m) − d * s,ω denotes the monochromatic data residuals for source s with d s, ω (m) and d * s,ω 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 where H k and γ k denotes an approximate Hessian and the step length, respectively.Operto et al. (2015) show that the gradient of the misfit function, ∇ m C k can be approximated by where the adjoint wavefields a 1 satisfies where R t augment with 0 are the residuals at the receiver positions in the full computational domain and I is the identity matrix.
Using the SEISCOPE toolbox (Métivier & Brossier 2016), we perform the optimization with either a pre-conditioned steepestdescent 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 where and (∂A h (ω, m)/∂m)p h (s, ω) are the so-called virtual sources (Pratt et al. 1998) whose radiation pattern are ∂A h (ω, m)/∂m.The prewhitening factor β, 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).1).

Subsurface parametrization
We use the (V 0 , ρ, Q, δ, ) parametrization to perform FWI.For this parametrization, a V 0 and 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 will be attainable.Therefore, 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 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 as a passive parameter in this study and focus on the update of V 0 , ρ and Q.As in Operto et al. (2015), we build a starting ρ 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.1).
resolution with which each parameter can be reconstructed is assessed by performing three separate mono-parameter inversions for V 0 , ρ and Q where the corresponding true models contains a V 0 , ρ 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 ρ 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 ρ 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 ρ 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 frequencydomain FWI may be sensitive to the numerical scheme used to discretize the radiation pattern matrices ∂A h /∂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 ρ updates particularly sensitive to this discretization as it is embedded in the operators X , Y and Z, eq. ( 5).In contrast, gradients of timedomain 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 ρ was built compared to Fig. 1(a) at the expense of a very large number of iterations.Ongoing work aims to build the ∂A h /∂m from a 27-point stencil to regularize the ρ update by decreasing the footprint of the discretization.

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 2 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 curves (Roth et al. 1998;Shtivelman 2004;Boiero et al. 2013Boiero et al. , 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 ≤ 7 Hz and 7 Hz < f max ≤ 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 preconditioned 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 l-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 , ρ, 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 ρ 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  1).The geological features intersected by each slice are reviewed in the caption of Fig. 7.
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 α. We set α 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 where ∇ m C denotes the scaled multiparameter gradient, g a small pre-set value and mmg 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 SEISCOPE toolbox.

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  1).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 l-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 l-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  1).The geological features intersected by each slice are reviewed in the captions of Fig. 7. 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 l-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 monofrequency 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 , ρ 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 attenuationfree 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 -ρ and V 0 -Q cross-talks.

Q models
The   1). 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 β of the pre-conditioner).Interestingly, we also show significant negative Q perturbations at the reservoir depths below the gas cloud (Figs 11a and b).  .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 ρ slices of tests FWI4-FWI6 (Fig. 12a).The vertical sections of the ρ 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 ρ in a (V 0 , ρ)-type parametrization.The ρ 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 ρ update: relative to FWI4, perturbations built by FWI5 are more focused in the shallow subsurface.
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 0 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 0 and ρ models are shown in Fig. 15 and show fairly well focused structures consistent with those shown in the V 0 models (Fig. 8).Structural differences between the V 0 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 0 , ρ, 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  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).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 ρ and Q updates are probably coupled to some extent, the significant standard deviations of Q translate into significant standard deviations of ρ 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 ρ 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 crosstalks, 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 ρ 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 ρ updates are mainly related to amplitudes variations.However, each Q and ρ 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 dash 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   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.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 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 Table 1.FWI experimental setup.h(m) and Optim: grid interval and optimization algorithm used during each multiscale step.pSTD: pre-conditioned steepestdescent.The superscript * means that V 0 , ρ and Q are jointly updated.Otherwise, only V 0 is updated.Q b : Q in the initial/background homogeneous Q model.Q 0 , β 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 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.

D I S C U S S I O N
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 frequencies shows reliable and consistent velocity updates provided that frequency grouping is used.We have shown how frequency grouping 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 test FWI6, although it is likely that increasing the number of frequencies 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 according to the available computational resources.The V 0 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 0 updates inferred from multiparameter reconstruction have not been significantly contaminated by ρ and Q cross-talks artefacts.A notable feature is that the ρ and Q updates absorb a significant part of the acquisition footprint hence cleaning the V 0 updates from this imprint.This allows us to use a regularization-free l-BFGS optimization, hence taking advantage of an improved convergence rate without boosting high frequency noise in the V 0 FWI models.The fact that the V 0 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 0 , ρ and Q).At low frequencies (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 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).

C O N C L U S I O N
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 frequency 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 contaminated by density and attenuation cross-talks.On the contrary, these amplitude-related parameters have been useful to absorb acquisition 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 correlate well with low-velocity anomalies associated with soft sediments and gas.Density seems more challenging to update.However, comparison 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 l-BFGS to truncated Newton optimization to improve the accounting of the Hessian operator during the early FWI iterations.

A C K N O W L E D G E M E N T S
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 research 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 allocation 046091 made by GENCI.We thank R. Brossier (IS-Terre), L. Combe (Geoazur), L. Métivier (ISTerre/LJK), A. Ribodetti (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.

Figure 1 .
Figure 1.Resolution and trade-off assessment.(a) From left to right, true model contains a V 0 , ρ, 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 0 (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.

Figure 2 .
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 lowvelocity 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).

Figure 3 .
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.4where 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 −1 .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.

Figure 4 .
Figure 4. Vertical sections of a V 0 FWI model (Table1, 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.

Figure 7 .
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 (Table1).
Five frequencies between 4 and 8 Hz are simultaneously inverted starting from the homogeneous background model.The V 0 , ρ and Q values are 1600 m s −1 , 1000 kg m −3 and 150 in the background model and 1800, m s −1 1200 kg m −3 and 10 in the inclusion.The intrinsic

Figure 8 .
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 (Table1).

Figure 9 .
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 (Table1).(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).

Figure 10 .
Figure 10.FWI perturbation models (Q).Depth slices at (a) 175 m depth, (b) 500 m depth, (c) 1km depth and (d) 3.35 km depth.In (a)-(d), each panel is labelled by the test number (Table1).The geological features intersected by each slice are reviewed in the caption of Fig.7.

Figure 12 .
Figure 12.FWI perturbation models (ρ).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 (Table1).The geological features intersected by each slice are reviewed in the captions of Fig.7.
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 Downloaded from https://academic.oup.com/gji/article-abstract/213/3/2037/4950495 by CNRS user on 17 April 2018

Figure 14 .
Figure 14.Validation of ρ 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 (Table1).(b) Vertical section of the final (top) V 0 and (bottom) ρ perturbation models (FWI4) across the well log (x = 6.23 km).The arrows point two low V 0 /ρ anomalies that are well matched on the density log.

Figure 15 .
Figure 15.Impedance model inferred from the FWI4 V 0 and ρ 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.

Figure 16 .
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 0 model.(b) Q model.(c) ρ model.The reflectivity of the mean V 0 model is superimposed on the mean V 0 model in (a) and the mean Q perturbation model in (b).

Figure 17 .
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).

Figure 18 .
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 (Table1).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.
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 (Table1).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.

Figure 19 .
Figure19.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 impulsional 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.

Figure 21 .
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).

Figure 23 .
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). .