Elastic versus acoustic 3-D Full Waveform Inversion at the East Pacific Rise 9°50’N

a


I N T RO D U C T I O N
The theoretical background of Full Waveform Inversion (FWI), formulated as an optimization problem that aims at minimizing differences between observed and synthetically calculated seismic data for attaining high-resolution velocity models of the subsurface, was laid out several decades ago (Tarantola 1984). At first, it was focused on obtaining 1-D velocity functions (e.g. Collier & Singh 1997;Singh et al. 1998). With fast technological advancements in data acquisition and increase in computational power, the development and application of different waveform inversion techniques were extended to 2-D (e.g. Pratt et al. 1996;Shipp & Singh 2002;Operto et al. 2006;Arnulf et al. 2011), and more recently to 3-D (e.g. Sirgue et al. 2010;Plessix et al. 2013;Warner et al. 2013;Raknes et al. 2015). However, most of the available 3-D FWI case studies assume acoustic media. A number of synthetic studies have been conducted to evaluate the effect of the acoustic assumption when applying FWI (e.g. Barnes & Charara 2009;He et al. 2018). The results suggest that the elastic effect could potentially have a big impact on the final velocity models of the subsurface obtained from acoustic FWI, especially in the presence of high-velocity contrasts (Barnes & Charara 2009;Pérez Solano et al. 2013).
We apply multiparameter 3-D acoustic and elastic FWI using a real, narrow-azimuth, 3-D multichannel seismic (MCS) data set to examine the differences between the two approaches and quantify their effects on the final inversion results. The data used in this study were collected in a deep-marine environment (the shallowest depth ∼2.5 km), at the East Pacific Rise (EPR) 9 • 50 N (Fig. 1), a fast spreading centre along which oceanic crust forms.
The upper portion of the crust formed at the EPR (first 1.5-2 km below seafloor) is represented by two layers: (1) layer 2A, which mainly consists of pillow basalts, and (2) layer 2B, represented by intrusive dikes (see Appendix A). Due to the maximum 6 km source-receiver offset available in our study, the penetration depth of the transmitted/diving waves is limited to the first 600-800 m of the upper crust (Marjanović et al. 2017). Thus, we focus mainly on recovering and comparing velocities from layer 2A and the transition zone up to the layer 2A/2B boundary. The reflectivity within layer 2 (upper crust) is relatively weak, and normal moveout reflection times could not be picked (see Appendix D). The application of FWI in this context is adequate since it can exploit the transmitted/diving waves that turn due to the velocity gradient. At low frequencies, the elastic properties' variations are of the size of the wavelength, indicating that elastic FWI is required according to the scattering C The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society.

1497
Downloaded from https://academic.oup.com/gji/article/216/3/1497/5210091 by guest on 24 February 2022 Figure 1. Seismic lines collected during 3-D multichannel seismic survey at the East Pacific Rise 9 o 50 N used in 3-D Full Waveform Inversions (FWI). The area covered by the 3-D grid for FWI is indicated by the red rectangle, with indicated inlines and crosslines. The above is superimposed on a regional seafloor bathymetry map, which is a composite derived from the Global Multi-Resolution Topography database (Ryan et al. 2009) and multibeam sonar data acquired during the cruise, gridded at 50 m. The red star in the globe inset shows the location of the EPR. theory. Furthermore, as the topmost basaltic layer is not covered by the sediments, a high velocity contrast is present at the seafloor, specially for S-wave. Thus, the EPR is a well-suited real example for conducting a comparative study between acoustic and elastic FWI.

Data acquisition and processing
The 3-D MCS data were collected during the 3-D seismic survey aboard R/V Marcus G. Langseth (MGL0812) in 2008 (Mutter et al. 2008;Fig. 1). The data were acquired in a racetrack geometry with 300 m separation between the sail lines at 82 • azimuth for the lines shot from west to east (i.e. 278 • azimuth for the east-west lines). A single sail line comprised four 6-km long streamers with 150 m separation, each containing 468 channels at 12.5 m group interval. Two gun arrays, containing nine guns each, with a total volume of 3300 in 3 (54.1 dm 3 ) in each array, were fired in alternating mode at 37.5 m interval. The nominal streamer and gun array tow depths were 10 and 7.5 m, respectively. The primary 3-D survey was composed of 115 sail lines and spanned over ∼27.5 km in the south-north direction (between 9å42 N and 9å57 N on the ridge axis). South of the primary 3-D box, a 14 line swath extended for ∼4 km (from 9å37.8 N to 9å39.9 N on the ridge axis). All of the collected lines were included in the inversion resulting in a total of ∼90 000 shots and covering an area of ∼1000 km 2 (excluding the gaps; Fig. 1).
The data were processed following a processing scheme developed for FWI application on narrow-azimuth, low-frequency seismic data (Milcik et al. 2014). The frequency analysis of the data set showed the presence of coherent energy at frequencies above 3 Hz, with low signal-to-noise ratio for frequencies below 5 Hz. Further processing was focused on removing noise from the data, notably swell and cable noise (see Appendix B). After de-noising we applied spectral boost to the data to flatten the spectrum and enhance low frequencies. We inverted the data within 2-6 km source-receiver offsets that mainly contain refraction/diving waves and the seafloor reflection. We did not include offsets smaller than 2 km to reduce the weight of the deeper reflection events (e.g. the presence of axial magma lenses at ∼1.5 km below the seafloor). An example of a processed shot gather for the above-mentioned offset range is shown in Fig. 2(a).
The source-wavelet estimate is an important component in FWI modelling. As a far field signature for the survey was not available, the source wavelet was obtained from stacking of near-offset traces from all available common mid-point gathers using water-velocity moveout and adjusting for variable topography. Both source and receiver ghosts were eliminated from the estimated wavelet, since our modelling includes free surface effect. Similarly to the observed data, we applied spectral boost to the source wavelet to enhance the low frequencies (Milcik et al. 2014).

3-D FWI
A comprehensive overview of the FWI technique, its theoretical background, different choices on minimization criteria, representation of the data in the misfit, linearization methods and different approximations in solving the inverse problem are presented by Virieux & Operto (2009), and references therein. Here, we focus exclusively on the comparison and evaluation of the results obtained by conducting both acoustic and elastic 3-D FWI (with least-squares waveform residual as the objective function) in the presence of high velocity contrast interface where strong elastic effects are expected. In both cases, we use the same starting velocity model (Fig. 3a). The size of the area covered by the velocity model is 2420 (44 × 55) km 2 ( Fig. 1), with 100 m horizontal and 25 m vertical grid point spacings. The inversion grid is bigger than the area covered by the data to account for the edge effects. To reduce the inversion artefacts (for instance, short wavelength variations caused by noise), the gradient is smoothed using a Gaussian filter over 200 m in horizontal, and 100 m in vertical directions.

Acoustic approach
The acoustic velocity model was obtained using a time-domain acoustic (isotropic) 3-D FWI formulation of the wave equation (see Appendix C). To establish an optimal inversion strategy for the EPR data environment, we conducted different sets of tests on limited 3-D data swath (∼2 km wide in latitudinal extent and centred ∼9 • 50 N; Fig. 1). All the tests were done for 10, 15 and 20 iterations and for three frequency bands: 2-3-5-7Hz, 2-3-7-9Hz, and 2-3-10-12Hz. However, we observed that by increasing the number of iterations and including higher frequencies (above 7 Hz) the  (c) show initial residual gathers obtained by subtracting calculated from observed data, using acoustic and elastic 3-D FWI inversion, respectively; in (d) and (e) we show residual gathers for acoustic and elastic 3-D FWI modelling, respectively, after 15 iterations. In addition to source-receiver offsets, for each panel we provide corresponding value of angle of incidence (θ ) valid for water bottom interface; with dashed line we indicate the value of the critical angle (∼38.5 o ). After reviewing the tests results, we carried out a simultaneous inversion of compressional velocity (Vp) and impedance (Ip) parameters using acoustic 3-D FWI for the entire data set within 2-3-5-7Hz frequency band, and performed 15 iterations (see Appendix C). Examples of initial and final residual signals obtained using acoustic inversion are shown in Figs 2(b) and (d), respectively;

Elastic approach
The elastic velocity model is obtained using the time-domain elastic (isotropic) 3-D FWI formulation. Following the exercises done for the acoustic approach, for the elastic 3-D FWI, we used the frequency range 2-3-5-7Hz, and run 15 iterations. Also, we used simultaneous, multiparameter approach, inverting for Vp, Ip, and the ratio between shear and compressional impedances (which also corresponds to the ratio between the shear and compressional velocities). The examples of initial and final residuals are shown in Figs 2(c) and (e), respectively. The resulting profiles of Vp and Vs (extracted from the ratio) are shown in Figs 3(e) and (f) and 4, respectively.

R E S U LT S A N D D I S C U S S I O N
While the conventional waveform modelling focuses on recovering the kinematics of the P-wave velocity in the subsurface, which is well represented by assuming acoustic formulation of the wave propagation, the physics argues for the presence of more complex, elastic Earth. The comparison of the average energy in the residual signal obtained using acoustic and elastic approximation to the wave equation after 15 iterations (e.g. Figs 2d and e) shows that the elastic approach accounts for ∼10-15 per cent more of the data signal when compared to the acoustically obtained one. This argues that the choice of elastic formulation of the wave equation and elastic physical model more appropriately explain the observed data. However, we need to acknowledge that not all of the observed signal is matched by the elastic approach. A portion of the remaining signal may be due to the presence of P-wave anisotropy, estimated to be ∼4 per cent for <1 km of the oceanic crust (Dunn & Toomey 2001), or due to attenuation (e.g. Yang et al. 2016), which we have not taken into account. Furthermore, additional iterations within the elastic inversion could have also reduced the observed difference. In addition to the above, it is important to note that, while the residuals after the elastic inversion show higher reduction in the overall misfit, locally amplitudes of the reflection signal seem to be explained better by the acoustic approach (Figs 2d and e). Interestingly, the most prominent, apparent better fit in the acoustic residual shows up around the critical angle for the water bottom interface (θ C ∼38 • or ∼4 km source-receiver offset). To examine the cause of this behaviour we conduct simple synthetic tests. The theoretically calculated curves of reflected/transmitted energy partitioning at the seafloor show that for incident angles >20 • the reflection coefficients of the reflected P wave for purely acoustic versus elastic media start to differ (Fig. 5). This difference increases with increasing angle of incidence to reach close to an order of magnitude after the critical angle. Under the acoustic assumption we thus model reflections that have higher amplitude than the real data. To reduce the misfit, the inversion therefore needs to adjust the velocity in such a way as to reduce the modelled amplitude (e.g. lowering the velocity as we can see in the P-wave acoustic model), leading to artefacts in the velocity model that have been reported by previous studies (e.g. Raknes et al. 2015). The anomalous velocity perturbations are also seen in the Vp model obtained after acoustic 3-D FWI (Fig. 3) and are addressed below.
While, in both acoustic and elastic 3-D Vp models, we clearly identify the base of the high velocity gradient zone, delimited by the Vp = 4.2 km s −1 contour, which we interpret as a lithological contact between layer 2A and layer 2B, the shape of the contact and inferred thickness of the layer 2A is very different in the two models. These differences originate exclusively from different approximations used in the wave equation: acoustic and elastic. In the final acoustic model, this boundary surface is extremely rough, with the presence of random velocity perturbations. In contrast, the final elastic model displays a relatively smooth surface. This difference becomes more evident when the effect of ridge topography is removed (Fig. 6). The velocity peaks present on the flanks in the acoustic FWI model, supposedly random intrusions of dikes are not likely to occur, not at the steady-state spreading centres, such as the EPR 9 • 50 N. Moreover, all the available seismic images of the layer 2A/2B boundary assembled at fast spreading centre display behaviour that is very similar to the one we observe in the elastic FWI model (e.g. Kent et al. 1994;Carbotte et al. 2000). The differences are even more prominent within the axial region (∼5 km around the ridge axis): while the elastic model provides relatively gradual increase in thickness of accumulated extrusives, the acoustic one shows more prominent increase, with ∼300 m of step in Figure 5. (a) A schematic representation of the incident/reflected/transmitted waves in the case of liquid-solid interface, that is, seafloor, assuming elastic (left) and acoustic (right) earth.ÀP marks down-going P wave,ÀPP´indicates reflected P-wave energy (blue),ÀPPÀ indicates transmitted P-wave (green), whereas APSÀ indicates transmitted S-wave energy (red) from the P-to-S conversion at the seafloor (only for the elastic case). The reflected and transmitted P waves for the acoustic case are shown in dashed lines. (b) For each of the identified phases, we show variation of their relative elastic/acoustic energy with angle of incidence for the fluid-solid interface, using the colour-code, indicated in the legend. The calculations are done placing the seafloor at 2500 m depth; for the water layer we assume Vp = 1500 m s −1 , Vs = 0 m s −1 and density 1 kg m −3 ; for the media directly below seafloor for the elastic case we assume Vp = 2400 m s −1 , Vs = 1300 m s −1 and density 2.6 kg m −3 ; for the acoustic case in the subsurface we assume Vs = 0 m s −1 . In (c) and (d) we show synthetically calculated shot gathers for the liquid-solid interface, assuming elastic and acoustic medium, respectively (as shown in panel a). In (e) we show the difference between (c) and (d).
thickness of the intrusive layer thickness within a narrow (<1km) region around the ridge axis (Fig. 3). This abrupt subsidence of dikes is not impossible, but we would expect it to be reflected in the shape of the seafloor. However, we do not observe any correlation between the dike subsidence and seafloor morphology.
The preference for the elastic over acoustic inversion for midocean ridge (MOR) environment has been suggested earlier from a 2-D waveform inversion modelling, using single streamer, and including source-receiver offsets from 1 to 3.3 km of downwardcontinued data set (Arnulf et al. 2014). However, except for local differences in the Vp amplitude, the two approaches provided an almost identical image of the subsurface ( fig. 13 in Arnulf et al. 2014). It has to be emphasized that the above study was done in 2-D for highly 3-D nature of the Mid-Atlantic Ridge. By conducting 3-D acoustic and elastic FWI using a short offset (<2 km) 3-D MCS data set, Raknes et al. (2015) show that besides differences in absolute velocities, the final velocity model obtained using acoustic approach is highly contaminated by artefacts, and observe a downward shift in the depth of layers well constrained by sonic logs.
In addition to the Vp model of the subsurface, the elastic 3-D FWI provides a Vp (Fig. 4) model. Thus, this result may potentially represent an additional step forward in studying MOR and can be used to better constrain other crustal properties, such as porosity. However, due to limited information on the S wave, which is mainly represented by P-to-S converted phases the absolute value of this parameter may not be well constrained, and only relative variations should be taken into account in interpretations.

C O N C L U S I O N S
We present a comparative study of acoustic and elastic 3-D FWI conducted on a 3-D seismic data set collected within deep-water setting and with high velocity contrast at the seafloor interface. The differences we observe in the two resulting compressional velocity models are significant. Besides the differences in absolute velocity values on the order of a couple of hundreds of meters per second, we also observe prominent structural variations with important implications on our understanding of the crustal accretion processes. The differences are also readily visible in shot gathers with a clear tendency for 'overfitting' of the S-wave signal in the case of the acoustic inversion approach. Here, we suggest the resulting elastic model as a geologically more plausible solution for the topmost oceanic crust formed at the EPR. Finally, our general recommendation is to apply elastic FWI when high velocity contrasts are present, and especially when mid source-receiver offsets are included in the inversion.

A C K N O W L E D G E M E N T S
This work is funded by Marie S. Curie grant awarded to MM within the call: H2020-H2020-MSCA-IF-2014. We would like to thank Shell Global Solutions International B.V. for all the support and means provided to conduct the project. Also, we are very grateful to Rose Anne Weissel (Lamont-Doherty Earth Observatory) and Siebren Reker (Shell) for their help with the seismic data transfer. Jörg Renner,Óscar Calderon Agudo and an anonymous reviewer provided thoughtful and constructive comments and corrections that helped improve the manuscript. Finally, we would like to thank PIs, Suzanne M. Carbotte and John C. Mutter, and co-PIs, Mladen Nedimović and Juan Pablo Canales, of the MGL0812 experiment. Seismic data from this study are archived with the IEDA MGDS (Mutter et al. 2008).

A P P E N D I X A : G E O L O G I C A L C O N T E X T
Most of our knowledge about lithological properties of oceanic crust formed at the EPR is based on: (1) a handful of 1-D velocity models obtained from seismic studies (e.g. Vera & Diebold 1994;Christenson et al. 1994), (2) results from ocean drilling programs (e.g. Alt et al. 1996) and (3) observations from ophiolites, exposed oceanic crust (e.g. Oman; Christensen & Smewing 1981). The results from 1-D seismic studies showed that the topmost part of layer 2A is characterized by a gradual increase in Vp, from 2.2 km s −1 at the seafloor to ∼3 km s −1 within the first hundreds of metres of the oceanic crust (Singh & Nicolas 2015). The lower portion of layer 2A (∼100 m) is dominated by the presence of a high vertical velocity gradient, within which the velocity increases for an additional 1-1.5 km s −1 . This portion of layer 2A is often referred to as the transition zone, the origin of which is poorly understood. At the boundary between layers 2A and 2B, the velocity reaches 4-4.5 km s −1 and increases to ∼5.5-6 km s −1 further down within layer 2B. In Fig. A1 we illustrate a simplified 1-D velocity function and corresponding lithological composition. Figure A1. Illustration of oceanic crustal structure formed at the East Pacific Rise 9 • N as seen in 1-D Vp model (modified from Vera & Diebold 1994 Plea), and its petrological interpretation.

A P P E N D I X B : DATA P RO C E S S I N G
We applied a standard marine data processing scheme that consists of removing swell noise and bad traces, band-pass filtering, and enhancing the low frequencies (Fig. B1). Regarding the latter, we first need to mention that the streamer data contain source and receiver ghosts, which degrade data frequency content. To compensate for the damping of the low frequencies by the presence of receiver and source ghosts, we derive a filter, equivalent to an inverse ghost filter, which we apply to the observed data. The same filter is applied to the source wavelet. With this approach, we could start the inversion at around 3 Hz.

A P P E N D I X C : I N V E R S I O N S T R AT E G Y
Multiple-parameter inversion is challenging because of the potential trade-off between parameters. In this work, the optimization is a quasi-Newton algorithm (L-BGFS-B; Nocedal & Wright 2006), and an initial inverse of the Hessian is estimated to compensate for the geometrical spreading. Moreover, we determine a gradient scaling per inverted parameter. This is done by trial and error by carrying out few inversions on a subset of the data and then analysing the results.
In terms of parameters for which we do the inversions, we choose to invert for compressional velocity (Vp), compressional impedance (Ip), and in the elastic case for the ratio between shear and compressional impedances. We decided to invert for compressional impedance contrast instead of density to have a better decoupling between the parameters, since the data contain pre-and pos-tcritical reflections, and diving waves (Fig. 2). The results of our tests show that the Vp model obtained from acoustic approach when we invert for Vp and density displays the presence of random, relatively small anomalies that cannot be explained by geology, which we define as artefacts. In contrast, when inverting for Vp and Ip these artefacts are not observed (Fig. C1).
In addition, in the elastic case we invert for impedance ratio because the data do not contain pure shear waves, and the main information on the shear parameters comes from the amplitudes Downloaded from https://academic.oup.com/gji/article/216/3/1497/5210091 by guest on 24 February 2022 Figure B1. An example of a shot gather recorded along four streamers: (a) raw data and (b) data after processing. In panel (c), we show the frequency content of the raw data (panel a) in blue line, and in red we show the frequency content of the processed data (panel b).
of the compressional waves. We should also mention that at low frequencies, the elastic parameter variations are of the order of the wavelengths. Therefore, we have significant interferences between the different converted phases.
Finally, we would like to mention that in this work, the wave equation is solved with a finite-difference scheme. In the elastic approach, we do not treat specifically the water layer, and thus the minimum shear velocity is around 1000 m s −1 , two-thirds of the water velocity, which leads to an elastic inversion roughly 10 times more expensive than an acoustic one.

A P P E N D I X D : A D D I T I O N A L I N V E R S I O N R E S U LT S
To illustrate the update in the velocity models resulting from the acoustic approach (update in Vp) and elastic approach (updates in Vp and Vs), in Fig. D1 we show the differences between the final and initial models that are presented in Figs 3 and 4.
As it can be observed the Vp and Vs models (Figs 3 and 4) obtained from the elastic 3-D FWI approach display very similar subsurface structures. This observation can be explained by the fact that the shear energy present in the data arises entirely from converted energy, and that the information is extracted mainly from the amplitude variation embedded in the compressional signal (as mentioned in Appendix C). However, we cannot completely discard the possibility that at some extent we have a crosstalk between these two parameters.
As we mention in Section 2.2, using the acoustic approach on a subset of the seismic data volume we conducted inversions for three frequency ranges: 2-3-5-7Hz, 2-3-7-9Hz and 2-3-10-12Hz (Fig. D2). The results of the tests show that by increasing the frequency we enhance the presence of small, random anomalies, which cannot be explained by geology. These anomalies most probably originate from the fact that the ill-posed character of the inversion problem becomes more prominent at higher frequencies. The elastic effects on the pre-and post-critical reflection amplitudes may become more dominant when we increase frequency. We did not Figure C1. Compressional velocity models obtained from 3-D FWI using acoustic approach and inverting for: (a) Vp and Ip (final velocity model) and (b) Vp and density. The arrows indicate the presence of random anomalies that we term artefacts. Figure D1. Updates in velocity models (i.e. differences between final and initial models) are shown for two inlines, 380 (left) and 420 (right). In (a) and (b) we show the updates for Vp obtained from acoustic 3-D FWI approach, whereas in (c) and (d) we show the updates for Vp obtained from the elastic 3-D FWI approach. In panels (e) and (f) we show updates in shear velocity obtained from the elastic 3-D FWI approach. carry out elastic inversions with the two last frequency bands (due to the computational cost).
The evaluation of our results is based on observations in data (shot gathers) and model domain, and comparison with the background knowledge on the MORs (in particular the EPR), which has been acquired over decades (e.g. Fornari et al. 2012). As it is described in Appendix A, the boundary between the layers 2A and 2B is represented by a high velocity gradient zone. The presence of a prominent velocity gradient zone, instead of a sharp velocity contrast, results in a complicated seismic event dominated by transmitted energy (note the absence of reflected energy in the stacks shown in Fig. D3). While this transmitted energy is crucial for obtaining velocity from seismic data, it does not allow us to conduct tests to further evaluate the results form 3-D inversion (e.g. produce common image gathers to examine flatness of reflectors, or compare migrated sections). Figure D2. Compressional velocity models obtained using acoustic 3-D FWI approach (Vp and Ip are the inverted parameters) for frequency ranges: (a) 2-3-5-7 Hz, (b) 2-3-7-9 Hz and c) 2-3-10-12 Hz. The arrows indicate prominent anomalies and random anomalies, which we define as artefacts. Figure D3. Seismic sections extracted from the 3-D pre-stack depth migrated data volume obtained using: (a) final Vp model from the acoustic approach and (b) final Vp model from the elastic approach. Note the absence of reflection energy for the entire upper crust. Some improvements in the seafloor reflection can be observed in the stacked seismic section migrated using Vp from the elastic approach (b).