# Nature of upper crust beneath the Lucky Strike volcano using elastic full waveform inversion of streamer data

A. F. Arnulf, A. J. Harding, S. C. Singh, G. M. Kent, W. C. Crawford

► To cite this version:

A. F. Arnulf, A. J. Harding, S. C. Singh, G. M. Kent, W. C. Crawford. Nature of upper crust beneath the Lucky Strike volcano using elastic full waveform inversion of streamer data. Geophysical Journal International, 2014, 196, pp.1471-1491. 10.1093/gji/ggt461 . insu-03581144

HAL Id: insu-03581144

https://insu.hal.science/insu-03581144

Submitted on 19 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Nature of upper crust beneath the Lucky Strike volcano using elastic full waveform inversion of streamer data

A.F. Arnulf,<fnref n="1" /> A.J. Harding,<fnref n="1" /> S.C. Singh,<fnref n="2" /> G.M. Kent<fnref n="3" /> and W.C. Crawford<fnref n="2" />

Accepted 2013 November 8. Received 2013 October 25; in original form 2013 July 9

## SUMMARY

Seismic full waveform is an emerging technique for determining the fine-scale velocity structure of the subsurface. Here, we present results of elastic full waveform inversion (FWI) along three multichannel seismic lines at the Lucky Strike volcano on the Mid-Atlantic ridge that provides a velocity image of the upper oceanic crust with unprecedented resolution (50–100 m). We have used a two-step process combining downward continuation with a time-domain, elastic FWI. The downward continuation procedure enhances the refracted arrivals and wide-angle reflections, and reduces the scattering noise due to rough seafloor. Since both sources and receivers are downward continued to the seafloor, the computational cost of FWI is reduced, as one does not need to model the thick water layer. Our results clearly demarcate two layers within seismic Layer 2A; a low-velocity, highly heterogeneous layer likely reflecting the complexity of accretion that is underlain by a more homogeneous high-velocity gradient layer. The base of Layer 2A is defined as a lithological boundary that can be offset by faulting. Thick (>400 m) units of anomalously low-velocity material (<2.5 km s⁻<fnref n="1" />) beneath different summital edifices on the central volcano indicate that a thick pile of high-porosity extrusive rocks can be supported without collapsing, suggesting that while in general there is pore closure with depth this is not the cause of high velocities we observe. Hydrothermal deposition sealing of small-scale porosity is shown to be a secondary process, which likely explains the upper crustal velocity increase with age, but is not responsible for the high-velocity gradient Layer 2A. Finally, the rapid thinning of the entire Layer 2A in the vicinity of the main normal faults suggests the tectonic thinning of a geologically defined layer, further confirming the lithological origin of the high-velocity gradient zone at the base of seismic Layer 2A.

**Key words:** Tomography; Composition of the oceanic crust; Controlled source seismology; Mid-ocean ridge processes; Submarine tectonics and volcanism.

## 1 INTRODUCTION

Oceanic crust is formed by a combination of magmatic and tectonic processes. Based on early seismic studies (Spudich & Orcutt 1980; Lewis & Garmany 1982; White & Purdy 1983) and observations from ophiolites (Nicolas et al. 1988), the igneous oceanic crust was divided into three layers: Layer 2A, Layer 2B and Layer 3. The uppermost (150–1000 m) Layer 2A has a low-velocity (~2.5 km s⁻<fnref n="1" />) and high-velocity gradient at its base that separates it from Layer 2B. Layer 2B has a P-wave velocity of 4.5–6 km s⁻<fnref n="1" />, and roughly encompasses the sheeted dyke complex. The lowest layer, Layer 3, has a small velocity gradient with velocity reaching up to 6.8 km s⁻<fnref n="1" /> at its base. Over the last few decades, the origin of the velocity transition at the base of Layer 2A has been interpreted as the lithological

boundary between extrusives and sheeted dyke complex (Toomey et al. 1990; Harding et al. 1993) or as a porosity limit within the extrusive section associated with an hydrothermal alteration front or a fracture front (McClain et al. 1985; Burnett et al. 1989; Wilcock et al. 1992; Christeson et al. 2007).

Layers 2B and 3 are well studied using ocean bottom seismometers (OBS; Cudrak & Clowes 1993; Van Avendonk et al. 1998; Canales et al. 2000; Seher et al. 2010a,b), but the nature of Layer 2A is poorly constrained by OBS arrays because the spacing of instruments is generally too large for detailed studies. Since there is a strong velocity gradient at the base of Layer 2A, Harding et al. (1993) stacked the turning rays present in seismic reflection data to determine the thickness of Layer 2A. By deploying both the sources and the receivers close to the seafloor during a NOBEL experiment,

<footnote n="1">¹Cecil H. and Ida M. Green Institute of Geophysics and Planetary Physics, Scripps Institution of Oceanography, University of California San Diego, La Jolla, CA 92093, USA. E-mail: aarnulf@ucsd.edu</footnote>

<footnote n="2">²Laboratoire de Géosciences Marines, Institut de Physique du Globe de Paris, Boîte 89, 4 place Jussieu, F-75252 Paris CEDEX 05, France</footnote>

<footnote n="3">³Nevada Seismological Laboratory/0174, University of Nevada, Reno, Reno, NV 89557, USA</footnote>
Christeson *et al.* (1994) provided the first well-constrained velocity structure of the uppermost several hundred metres of the crust. By using multiple components of the seismic wavefield Christeson *et al.* (1997) modelled the shear and compressional-wave structure of Layer 2A at the East Pacific Rise, but resolution remained a major issue. Collier & Singh (1998) and Hussenoeder *et al.* (2002a,b) performed 1-D waveform inversion of surface multichannel streamer data, obtaining better resolution but analysis was restricted to areas of smooth seafloor where the 1-D assumption could be justified.

Lucky Strike volcano lies on the Mid-Atlantic Ridge (MAR) south of the Azores triple junction at 37.3° N. Lucky Strike has several characteristics typical of a slow spreading ridge segment, including a full spreading rate of ~2.1 cm yr⁻¹ (Demets *et al.* 1994), a 15- to 20-km-wide central valley, an axial high at the centre and nodal basins near the segment ends (Fig. 1a). The Lucky Strike volcano, located at the centre of the segment, is ~7 km in diameter, hosts an active hydrothermal field (Langmuir *et al.* 1997) and presents a wide variety of recent volcanic deposits that we will present in the geological setting section. From a geochemical standpoint, basalts from the volcano suggest an enhanced melt supply from the Azores hotspot (Moreira & Allegre 2002) and a well-homogenized and long-lived magma chamber system (Moreira *et al.* 2011).

The multichannel seismic (MCS) data were collected in 2005 onboard the R/V l'Atalante as part of the SISMOMAR experiment (Singh *et al.* 2006). The data set consists of 39 lines perpendicular to the spreading axis and spaced 100 m apart (Fig. 1b). The data were acquired using a 4.5-km-long hydrophone digital streamer towed at a nominal depth of ~14 m. The streamer was composed of 360 channels with 12.5 m receiver intervals. The source was fired each 37.5 m and consisted of a 42.5 L (2594 cubic inch) 18-airgun array towed at a depth of ~14 m and tuned to provide a broad-band (5–50 Hz) energy with a strong first arrival (Fig. S1). A total of 19 942 shots were fired during the experiment and data were recorded in 11-s-long records sampled every 2 ms. Initial investigations of the reflection data concentrated on the imaging of a ~7-km-long, ~3-km-wide axial magma chamber (AMC) reflector at a depth of 3.4 km below the seafloor (bsf; Singh *et al.* 2006; Combier 2007), and on the imaging of the Layer 2A reflection using conventional methods (Combier 2007; Seher *et al.* 2010c). OBS-based tomographic studies have revealed two low-velocity zones (LVZs) within the axial valley: an upper crustal LVZ that has been linked to a recent high-porosity formation (Seher *et al.* 2010b) and a lower crustal zone beneath the AMC reflector that has been identified as a zone of partial melt (Seher *et al.* 2010a).

To determine the fine-scale structure of the upper crust, Arnulf *et al.* (2011, 2012, 2013) subsequently reprocessed the MCS data using a combination of downward continuation and 3-D traveltime tomography (Arnulf *et al.* 2011, 2013) and using full waveform inversion (FWI) of one selected seismic line (Arnulf *et al.* 2012). The combination of downward continuation with traveltime tomography and ultimately FWI provides a high-resolution velocity image of the upper ~1 km of crust beneath the volcano, with each different method providing a different insight into the data through its different resolution/accuracy. Based on the 3-D traveltime tomography results, Arnulf *et al.* (2011) proposed that the porosity boundary between the relatively low-porosity lava lake (LL) and the surrounding high-porosity volcanic cones could act as a preferential pathway for the upgoing hydrothermal flow. Based on full wave inversion of a single line north of the volcano summit, Arnulf *et al.* (2012) suggested that the young shallow crust beneath the Lucky Strike volcano is formed from a stack of

volcanic sequences emplaced within the axial valley. At last, from an extensive study of the 3-D traveltime tomography model, Arnulf *et al.* (2013) inferred that the entire upper part of the central volcano is young relative to the underlying median valley floor and suggested that hydrothermal sealing of small-scale porosity is progressing at normal to enhanced rates.

Here we extend the work of Arnulf *et al.* (2012) to a broader region of the Lucky Strike volcano and explore the fundamental links between magmatic, tectonic and hydrothermal processes. We invert three seismic profiles (Fig. 1b) that span the 3-D seismic reflection data box (Singh *et al.* 2006) using a FWI procedure. These profiles were chosen to explore the upper oceanic crust beneath several important features of the central volcano. Seismic line 14 spans the southern part of the central volcano, which appears to be slightly less tectonized (Fig. 1b) and crosses the southern summital volcanic edifice (S, Fig. 1). Seismic line 24 crosses the northern part of the LL, the hydrothermal field and the southern flanks of the two northern summital volcanic highs (NW and NE, Fig. 1). Seismic line 39, was already presented by Arnulf *et al.* (2012) and crosses the nascent central graben (CG) deforming the northern part of the volcano (CG, Fig. 1). At all these locations, the streamer recorded not only significant refraction energy from Layer 2A at far-offsets, but also strongly scattered energy from the rough igneous seafloor (Fig. S1b).

Waveform inversion consists of minimizing a misfit function between observed and synthetically calculated seismograms, sample by sample (Shipp & Singh 2002). A serious disadvantage of this technique is that it attempts to fit all parts of the waveform equally whether it is signal, noise or artefact. So, in areas of rough igneous seafloor, scattering noise can easily dominate the misfit and thus the inversion results. A common way to tackle this problem is to invert only energy that arrives ahead of the seafloor reflection on the MCS data (Canales 2010) but in deep-sea environments this part of the data does not contain information on shallowest part of the crust. Arnulf *et al.* (2011) proposed extrapolating the surface wavefield to the seafloor, mimicking a pure sea-bottom experiment (Christeson *et al.* 1994) and thus bringing the crustal refraction energy out ahead of the seafloor reflection energy. The details of this method are discussed in Arnulf *et al.* (2011, 2013), and will not be presented here.

In the following sections, we first present a description of the geological setting that is specific to the Lucky Strike area, then we present a brief outline of the FWI method that is used in this study.

## 2 GEOLOGICAL SETTING

Lucky Strike volcano is a prominent seamount at the centre of the Lucky Strike segment of the MAR (Fig. 1). Multibeam bathymetry shows the volcano is divided into two parts separated by a NNE-trending, and 50- to 100-m-deep graben (CG, Fig. 1b). The western part of the volcano features an elongated narrow ridge that Fouquet *et al.* (1994) have called the western volcanic ridge (WVR, Fig. 1b), while the eastern part has a semi-circular shape (eastern ridge: ER, Fig. 1b). The volcano summit rises 0.3–0.4 km above the rift valley floor and features an LL surrounded by three volcanic cones, namely: southern (S), north eastern (NE) and north western (NW) summital volcanic cones (Fig. 1c).

The LL has a circular shape, is about 300 m in diameter and exhibits very fresh low vesicularity massive lava flows. In contrast, most of the three summital volcanic cones as well as the WVR are composed of older, highly vesicular volcanic breccias and pillows
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Bathymetric maps of the Lucky Strike segment with important features of the seafloor. (a) Regional bathymetric map of the Lucky Strike segment showing the central volcano in the middle of the axial valley. The black rectangles mark the extent of the area shown in (b) and (c); the inset globe shows the location of the Lucky Strike segment on the MAR. (b) Selected portion of the ridge segment with illumination across the ridge axis. Black lines highlight faults, dashed green lines correspond to the unrifted backside of some rifted seamounts (RC). The cyan rectangle shows the extent of the SISMOMAR 3-D survey. The blue lines mark the extent of the seismic structures presented in this paper; from south to north: seismic lines 14, 24 and 39. (c) Selected portion of the volcano summit with illumination across the ridge axis. White lines outline the different summital constructional features of Lucky Strike volcano [S, southern; NE, northeastern; NW, northwestern volcanoes; WVR, western volcanic ridge on (b)]. The red dots mark the hydrothermal vents that surround the central lava lake.</figcaption></figure>
(Fouquet et al. 1994; Ondréas et al. 1997, 2009; Humphris et al. 2002). Active hydrothermal sites, including high-temperature black smoker chimneys venting fluids up to 333 °C and low temperature diffuse flow (Langmuir et al. 1997) surround the LL (Figs 1b and c). These hydrothermal sites consist of massive sulphide deposits and standing inactive and active chimneys.

Two major normal faults bound the rift valley, with vertical throws of >500–900 m (e.g. the eastern and western bounding faults—EBFs and WBFs, Fig. 1b). DSL-120 side scan data indicate that Lucky Strike volcano is highly tectonized and is cut by numerous faults (Humphris et al. 2002; Ondréas et al. 2009). Each of the three cones on the summit is cut by prominent faults that clearly post-date construction of the cones. Faults on the summit have predominantly ridge-parallel orientations, vary in length from 200 m to 3 km, and exhibit spacings from tens to hundreds of metres. In the depression between the cones, volcanic flows of the LL truncate faults that are presumed to underlie the most recent volcanism (Humphris et al. 2002; Ondréas et al. 2009).

Furthermore, from different reflection images of the central volcano, Arnulf et al. (2011) imaged a new shallow reflector lying at 100–250 m below the seafloor, distinct from the usual strong Layer 2A reflector at depth (marking the transition between Layers 2A and B), which they call Layer 2Aa. This reflector was most clearly seen beneath the volcano summital edifices, and Arnulf et al. (2011) suggested that it might correspond to the boundary between lava sequences of different ages. The second deeper reflector was renamed 2Ab and corresponds to conventional Layer 2A.

## 3 ELASTIC FWI OF DOWNWARD CONTINUED DATA

Elastic FWI is an attractive method that attempts to match the observed seismograms in order to discover the elastic parameters of the Earth. Instead of exhaustively or randomly exploring the whole model space, efficient 2-D waveform inversion methods utilizes an iterative method of local search to solve the non-linear problem. Convergence to the global minimum rather than a local minimum however requires an accurate starting model that contains the long and, possibly, the intermediate wavelengths of the true velocity model (Neves & Singh 1996). Given a good starting model, the key to waveform inversion is the use of an efficient method to formulate the gradient operator. Different approaches have been formulated either in the frequency domain (Pratt & Worthington 1990; Pratt 1999; Sirgue & Pratt 2004) or in the time domain (Tarantola 1984, 1986; Shipp & Singh 2002), here we use the method developed by Shipp & Singh (2002).

### 3.1 Theory

The elastic FWI algorithm is based on the iterative minimization of a misfit function between observed and synthetically calculated seismograms, point by point, in a least-squares sense, where the misfit function S is formulated (Shipp & Singh 2002) as

$$S = \sum_{\text{shots}} \int_{0}^{T} dt \sum_{\text{recs}} [d_{\text{syn}} - d_{\text{obs}}]^2, \quad (1)$$

where $d_{\text{syn}}$ and $d_{\text{obs}}$ are, respectively, the forward-modelled synthetic data and the observed data, with sample differences summed over the total recording time T.

To iteratively solve the inversion problem we have used a traditional approach where the steepest descent gradient is condi-

tioned using conjugate gradients. The conjugate gradient method uses knowledge gained from previous iterations to ensure that search directions are not repeated, resulting in a faster convergence. To ensure an efficient calculation of the steepest descent gradient direction, we used the formulation of Mora (1987), in which the gradient direction for each Lamé parameter, $\delta L$, is derived from the cross-correlation of the forward synthetic wavefield $\vec{u}$ and the backpropagated data residuals wavefield $\vec{\psi}$:

$$\delta L = \sum_{\text{shots}} \int_{0}^{T} dt (\vec{u}, \vec{\psi}). \quad (2)$$

The gradient g for the P- and S-wave velocity models is a linear combination of the Lamé parameter gradients. An optimal linearized step length $\alpha$ is computed in the manner of Pica et al. (1990), requiring an extra modelling step. The updated model m at the iteration n + 1 is then computed as follows:

$$m_{n+1} = m_n - \alpha_n g_n. \quad (3)$$

Our method allows us to compute the gradient of the misfit function using three forward modelling steps (one to create the synthetic seismograms and calculate the data residuals; one to backpropagate the data residuals and compute the gradient direction; and a final one to compute the step length). An efficient forward modelling algorithm is thus required to obtain an accurate solution of the inverse problem.

Our forward modelling algorithm uses the fourth order in space second order in time finite difference solution of the 2-D isotropic elastic wave equation (Levander 1988). The main advantage of the time-domain approach is that a large frequency band can be modelled and inverted simultaneously, unlike frequency domain implementations that invert for a range of discrete frequencies one at a time (Pratt & Worthington 1990; Pratt 1999; Sirgue & Pratt 2004). On the other hand, a drawback of the time-domain implementation is the large computational requirement of the different modelling steps, which can be significantly reduced using a downward continuation technique (shots and receivers closer from the seafloor (Arnulf et al. 2012) or an injection grid method (Robertsson et al. 2000; Robertsson & Chapman 2000; Royle 2010).

### 3.2 Practical considerations

The strength of waveform inversion is that it can use the full seismic wavefield to create a single self-consistent, high-resolution velocity image of the crust that can be used for geological interpretation. However, limitations of the method can arise from the data (noise and scattering effects), from the forward-modelling engine and from the local approach by which the gradient is estimated; thus one has to satisfy several conditions to achieve a satisfactory convergence. In this section, we discuss the required inputs for the inversion, including pre-processed downward continued data (Fig. S1) and the source wavelet estimation (Fig. 2). We also present a rigorous convergence analysis of the iterative method. The creation of the long-wavelength background velocity model, used as the starting model for waveform inversion (S1a), was described previously by Arnulf et al. (2011, 2013).

#### 3.2.1 Data processing during the Synthetic Ocean Bottom Experiment (SOBE) method

The main difficulty in applying FWI to mid-ocean ridge data is the ubiquity of high-amplitude seafloor scattering that obscures later
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Initial source estimate and updates. (a) Trace obtained from stacking the nearest three traces of several shots over a region of relatively flat seafloor (b). Only the signal present between 0.13 and 0.18 s was selected and the secondary reflections were removed before processing the initial source wavelet. (c) Initial bandpass filtered source wavelet (pink line), first (green line) and final updated source wavelet (black line). (d) Normalized amplitude spectra for the source wavelets. For comparison the dashed black line shows the frequency spectrum of a ricker wavelet with a dominant frequency of 9 Hz.</figcaption></figure>

arrivals and is hard to predict deterministically. One way to tackle this problem is to carefully window the data and to invert only the energy ahead of the seafloor reflection (Canales 2010; Arnulf *et al.* 2012). The advantage of using the downward continued data for waveform inversion is that the wide-angle reflections and refractions are moved in front of the water-wave reflections and scattering and thus can be isolated and inverted, independently of the details of the basement interface (Fig. S1). In fact, since the scattering is mainly from the seafloor and is intrinsically 3-D, the downward continuation acts as a 2-D filter to improve the signal-to-noise ratio.

In order to avoid the generation of surface waves at the seafloor during modelling, the surface seismic data were downward continued to a horizon 75 m above the seafloor using a Kirchhoff implementation of the downward continuation operator (Berryhill 1984; Shtivelman & Canning 1988; Larkin & Levander 1996) as surface waves could dominate the misfit calculation and mislead the inversion. Pre-processing of the field data was accomplished during the downward continuation stage. A fourth-order Butterworth bandpass filter with corner frequencies of 3 and 10 Hz was first applied to constrain the frequencies of the downward continued data to ones that satisfy the stability criterion of the 12.5 m finite-difference grid (Fig. S1). Second, the observed signal was convolved by $H(t)/\sqrt{t}$ (being the Heaviside step function and $t$ being the time) to boost low frequencies with respect to high frequencies. Then, the wavefield of each shot gather was spatially resampled on the finite difference grid in such a way that the location of each shot corresponds to the closest grid node, and the location of the attributed 360 receivers extend from 12.5 m to 4.5 km behind the shot location. To avoid numerical dispersion within the 12.5 m finite-difference grid, the downward continued shot gather were also temporally resampled within a 2.5 s time window using a time increment of 1 ms. As a valuable by-product, the new geometry deployed in an SOBE experiment using the SISMOMAR data improves by a factor of 2 the efficiency of the time-domain FWI by deploying shots and receivers ~1.5 km deeper within the water column.

The waveform inversion concentrated on the wide-angle between the first refraction arrival and the seafloor reflection/direct water arrival that is present in the modelled output but not in the downward

continued data, Fig. 3. To this end, first arrival and seafloor reflection times generated from the starting model were used to define two surgical mutes, top and bottom, and carefully window the wide-angle data, Fig. 3. The final step in the data pre-processing was the computation of a single data-scaling factor by least-squares minimization of the amplitude difference between a windowed set of observed and synthetically calculated shot gathers. This scaling was adjusted during inversion as part of the source wavelet update (see later).

For the upper oceanic crust, the frequency of the refraction arrivals is notably lower than the seafloor reflection. In 1-D waveform inversions, this results in low estimated attenuation ($Q$) particularly for Layer 2A (Collier & Singh 1998), although, in part, the lower frequencies could be attributable to scattering by the heterogeneous crust. With our elastic FWI inversion scheme we deliberately restricted attention to arrivals, refractions and wide-angle reflections, whose spectra exhibit only limited offset dependence (Fig. S2). As a result, we assumed for our elastic FWI that the principal attenuation effects can be incorporated into an effective source wavelet that is updated as part of the inversion process. Although the initial wavelet was broader band, the final wavelet captures most of the spectral content of the filtered data.

### 3.2.2 Source estimation

The synthetically calculated wavefield depends on the source wavelet injected into the finite-difference grid. An error in the input source wavelet thus affects the synthetically calculated data and the misfit. Ideally the source signature should be measured during the MCS survey (Ziolkowski 1991), however if no source wavelet was recorded, a method for estimating an appropriate source is required. In this study, we followed the method presented by Collier & Singh (1997) to recover an appropriate starting source wavelet from stacks of the nearest three traces of several shots over a region of relatively flat seafloor. Fig. 2(a) shows the resultant trace obtained from the summing of multiple adjacent near offset traces in the region of interest. Several successive reflectors are observed. The initial
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Initial forward modelling for shot gather 38662 (location shown on Fig. S1). (a) Observed downward extrapolated shot gather with its corresponding top and bottom mutes (blue/green and blue lines, respectively). (b) Synthetic shot gather modelled from the starting velocity model but the first updated source (green line in Fig. 2c). Pink stars on the mute window mark the extent of the wiggle plots presented in (f, g and h). (c) Frequency spectrum of the observed data computed within the mute window. (d) Corresponding frequency spectrum of the synthetic shot gather. (e) Difference between the two frequency spectrums. The black line delimits regions where the differences are greater than 25 per cent. (f, g and h) Wiggle plots for three selected traces. Blue line and green line correspond, respectively, to the observed and synthetic data.</figcaption></figure>

pulse of the seafloor reflection present between 0.13 and 0.18 s in the stacked trace was selected as the unfiltered source wavelet and then subjected to the same Butterworth bandpass filters used during downward continuation. We thus started the waveform inversion with an approximately minimum phase wavelet with a correct dominant frequency (Fig. 2). The source wavelet was updated during the inversion by running special iterations of the waveform inversion that used the residual wavefields backpropagated to the source positions to calculate a source gradient (Mora 1987) rather than a model gradient (eq. 2; Figs 2c and d). The updated source wavelet reduced the average misfit and better incorporated any attenuation effects. Finally, no correction of the source or receiver ghosts were

applied to the source wavelet because in the SOBE geometry the seismic sources are located close to the seafloor, far from the sea-level and the free-surface ghosts are not removed during downward continuation.

### 3.2.3 Inversion strategy and misfit reduction

Waveform inversion principally targeted the energy related to P-wave refraction events and some wide-angle reflections present in the downward continued MCS data for offsets <3.45 km (Fig. 3). We followed a multistage FWI strategy alternating between model
updates, where we inverted simultaneously for $V_p$ and $V_s$, and source updates. We did not invert for density, which was updated using empirical relationships with P-wave velocities (Gardner et al. 1974; Hamilton 1978). We typically started the inversion with a run of eight model updates where all the frequencies present in the downward continued data were inverted simultaneously (Fig. S1). This was followed by a pair of source updates (Fig. 2) following the method of Mora (1987). The updated source was further adjusted to minimize the amplitude difference between the observed and synthetically calculated shot gathers. A set of eight model updates using the updated source was then run prior to a further source update. The inversions used every second shot from a line, in other words one shot every 75 m and approximately 20 model updates were run for the three different seismic lines. A 150 m (horizontal) by 37.5 m (vertical) Gaussian smoothing was applied to the gradient updates during inversion in order to remove short wavelength artefacts.

We estimated a starting S-wave velocity model using the empirical relationships found in Castagna et al. (1985) for low values and Fliedner et al. (1998) for higher P-wave velocities greater than 3.5 km s⁻¹. While the $V_p/V_s$ ratio is reasonable (1.9) for $V_p$ greater than 3500 m s⁻¹, we introduced two optional parameters to restrict the S-wave starting velocity model in the upper crust. The first parameter was used to restrict the maximum $V_p/V_s$ ratio allowed within the low-velocity basaltic crust and the second to control the minimum S-wave velocity in Layer 2A. After several tests, we observed that inversions with higher S-wave velocities in the shallow crust gave a slightly better misfit (a few per cent) and were better numerically because the dispersion limit, the requirement of a minimum number of grid points per wavelength, meant the maximum frequency could be pushed higher. Therefore, following the tests, we decided to restrict the $V_p/V_s$ ratio to 3, limit the minimum S-wave velocity to 800 m s⁻¹ and increase the bandwidth for subsequent inversions.

Waveform forward modelling using the first updated source wavelet (Fig. 2, green line) through the starting tomographic velocity model predicted first arrival traveltimes within ~10 ms (Fig. 3), consistent with tomographic inversion estimates (Arnulf et al. 2011, 2013). However, some significant phase misalignments, signal amplitude differences and amplitude spectrum inconsistencies were still present between the observed and modelled data (Fig. 3), which the FWI is expected to reduce. In fact, by definition the goodness of the fit between the observed and modelled data is given by the misfit value, which includes the amplitude and phase data residuals (Shipp & Singh 2002; Brenders & Pratt 2007; Sears et al. 2008, 2010; Canales 2010). In this study of three lines, the misfit was reduced, on average, by 75–80 per cent in approximately 20 iterations. The resulting final velocity models and the corresponding updated source wavelets greatly improved the amplitude, phase and frequency content of the final modelled data (Fig. 4), reducing the size of the data residuals (Fig. 5) and thus the average misfit value.

One could use the much smaller misfit to argue that the elastic FWI has achieved a satisfactory convergence to a minimum of the misfit function. However, reliance solely on the misfit reduction may not be sufficient to rigorously assess the convergence of the iterative method. Indeed, it is quite easy for the initial misfit to be 200 per cent of the observed data if the model predictions are the same size as the data and the average correlation between observations and model data is near zero, a not uncommon situation due to phase misalignments between the two data sets. In these circumstances a 50 per cent misfit reduction can be achieved by the FWI by reducing

the amplitude of the predicted data to almost zero; the FWI is then doing no harm but is performing no better than predicting nothing at all.

A more rigorous convergence analysis is required to assess the goodness of the fit between the observed and modelled data and validate the inverted model. Apart from assessing the reduction in the misfit value, we also tracked the size of the relative misfit, ratio of misfit to signal size. In this study, the size of the relative misfit was reduced from 130–150 to ~30 per cent in ~20 iterations. A successful inversion should have a low relative misfit at the end irrespective of the starting misfit. The other parameter we calculate to estimate the robustness of the inversion is a correlation coefficient value between the observed and modelled data (Minshull & Singh 1993). The correlation coefficient provides a direct estimation of the phase residuals and is amplitude invariant, and therefore supplies complementary information about the fitness. In this study, the average cross-correlation coefficient was increased from 0.25–0.35 to ~0.80 in ~20 iterations, testifying to the good similarity between the observed and modelled data (Fig. 5).

Fig. 6 presents a convergence analysis for seismic line 39. After 20 iterations, the initial absolute misfit value is reduced by ~76 per cent from 2830 to 667. The total data size is 2105, and thus the relative misfit is reduced from 134 to 31 per cent. Tracking of the misfit reduction shows that after 20 iterations the convergence has plateaued. Moreover, the final residuals are less than the data size for most of the shots, and the correlation coefficient is now positive for all shots and shows a high average value of ~0.8. For the shots with highest correlation coefficients, some secondary arrivals are correctly reproduced (Figs 4 and 5). In contrast, the shots with the lowest correlation coefficients appear to be located above the bounding walls of a nascent graben and thus correspond to the most complex regions of the velocity structure. A comparison of the elastic and acoustic FWI results will be covered more fully in a dedicated section.

## 4 RESULTS

Along the three profiles, the region covered by the waveform results includes the Lucky Strike volcano and part of the surrounding median valley from seafloor depths down to ~1 km bsf (Fig. 1). In order to highlight different aspects of the inversion results, we show the P-wave velocity structure, the P-wave velocity anomaly and the vertical velocity gradient for each line. We show the inverted S-wave velocity structure along seismic line 39 in Fig. S3, but we do not discuss extensively the results because we believe there is not enough information in the data to fully constrain S-wave velocity, although S-wave velocity information is necessary for proper convergence; we will discuss this later in the paper. Velocity anomaly sections were computed by subtracting a 1-D crustal velocity profile hung from a smoothed seafloor bathymetry. The 1-D model was the average profile of the 3-D tomographic cube (Arnulf et al. 2013). The vertical velocity gradient sections bring out the slope in the velocity models. At the base of Layer 2A a 300- to 400-m-thick high-velocity gradient region (>3 s⁻¹) marks the transition with Layer 2B. We delimited, respectively, the top and bottom of this transition region by a dashed black line and a plain black line on Figs 7–9. Similarly, we tried to delimit regions of anomalously low Layer 2A velocities (<2.5 km s⁻¹) from the more common Layer 2A velocities (>2.5 km s⁻¹) by a dotted line. At last, we inferred the extent of the different major faults recognized on the bathymetry using some yellow arrows (Figs 1 and 7–9).
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Final forward modelling for shot gather 38662. Same as Fig. 2 for the final synthetic data modelled within the best-fit velocity structure using the last updated source wavelet (black line in Fig. 2c).</figcaption></figure>

## 4.1 MCS line 14: the southern volcanic edifice

The MCS line 14 is the southernmost of the inverted seismic profiles and crosses the Lucky Strike volcano perpendicularly to the ridge axis, 700 m south of the centre of the LL. This line crosses the less tectonized southern half of the central volcano as well as the southernmost summital volcanic edifice (Fig. 1b). Overall our FWI model reveals velocities ranging between 1.7 and 5.4 km s⁻¹ with the lowest velocities (<2.5 km s⁻¹) mainly observed within the upper 200–450 m beneath the central volcano (Fig. 7a). The inverted model shows the presence of a persistent 50- to 100-m-thick low-velocity layer within the median valley. The western side of the central volcano, which also has a more tectonized seafloor (see Fig. 1), has a rougher velocity structure with several inclined features that may be related to fault zones (yellow arrows, e.g. −2 km distance range, Fig. 7). The vertical velocity gradient section (Fig. 7c) brings out the fine scale structure as well as the upper and lower boundary of the velocity gradient zone, that is, base of Layer 2A.

The Layer 2A gradient zone is less heterogeneous than the shallower structure and is about 400–500 m thick. The top of Layer 2B is recognized by the transition to a lower velocity gradient (<1.5 s⁻¹). It has a relatively homogeneous structure, and an absence of layered sequences. In addition, the velocity anomaly section highlights several low-velocity regions in the upper crust, including an ~800-m-thick, 1-km wide vertical region located underneath the southernmost summital volcanic edifice (Fig. 7b). Fig. 7(d) shows that Layer 2A velocities along seismic line 14 are on average increasing with age but at a rate slower than a reference aging curve (Carlson 2004; Arnulf et al. 2013). It is also worth noting that across axis, velocity horizons <3.5 km s⁻¹ appear to migrate upwards from the high gradient region but leave a persistent low-velocity surface layer in the first ~100 m of the median valley. Finally, velocities within the upper Layer 2B range between 4.5 and 5.3 km s⁻¹ with the fastest velocities being found in the axial valley floor in the vicinity of the bounding normal faults at −5 and 4.5 km distance (see Figs 1 and 7d).
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Evolution of data misfit for shot gather 38662 during waveform inversion. (a) Downward continued data at a datum 75 m above seafloor. Only refractions and reflections observed in front of the water wave were inverted. (b) Initial synthetic. (c) Final synthetic. (d) Final source wavelet. (e) Initial residuals (~153 per cent of the data size). (f) Final residuals (~11 per cent of the data size). The correlation coefficient between the synthetics and data was initially 0.07; it increased to 0.945 for the final model.</figcaption></figure>

The upper part of Layer 2A can be subdivided into a low-velocity anomaly region with velocity <2.5 km s⁻¹ and an underlying region with more common Layer 2A velocities (above and below the dotted black line, Fig. 7). Within the axial valley floor, this horizon lifts away from the top of the main gradient zone at the base of Layer 2A, and could mark the upward progression of hydrothermal mineralization processes through the basaltic crust (Peterson *et al.* 1986). Within the central volcano, this horizon coincides with the top of the gradient region, suggesting that the entire upper part of the central volcano is young and perhaps the top of the gradient region represents a palaeoseafloor below which hydrothermal alteration has increased velocities above 2.5 km s⁻¹. This horizon could correspond to the Layer 2Aa reflector described by Arnulf *et al.* (2011) (discussed further later in relation to line 39).

The main low-velocity anomaly observed beneath the southern submittal volcanic cone appears to be correlated with the most recent constructional activity that significantly disrupts the fabric of the lower, gradient region part of Layer 2A resulting in an ~1 km s⁻¹ anomaly. The root of this pipe-like feature of strong low-velocity anomaly could be explained by the presence of highly fractured and porous material (Arnulf *et al.* 2011).

## 4.2 MCS line 24: the LL and the surrounding hydrothermal vents

Seismic line 24 has a slightly more complex seismic structure. It crosses the central volcano 300 m to the north of the LL centre and includes the southern part of the WVR, the two northern summital

volcanic cones (NW and NE), the northern part of the LL and the surrounding hydrothermal field as well as an old rifted cone (RC) that lies on the eastern bounding wall (Figs 1 and 8).

The P-wave velocities (Fig. 8a) range between 1.7 and 5.5 km s⁻¹, with the lowest velocities (<2.5 km s⁻¹) observed in the topmost 200–250 m underneath the central volcano and the old RC. Usually low velocities (<2.5 km s⁻¹) are present beneath the two summital volcanic cones, the WVR, the eastern flank of the central volcano between 1.9 and 2.5 km distance and beneath the eastern valley footwall. Interestingly, the velocity structure along line 24 shows a relatively thinner layered formation (200–250 m on average) as compared to line 14 (>300 m), possibly indicating enhanced magmatic eruptions on the southern half of the central volcano. A lack of horizontal continuity is observed in the upper crustal units, which can be explained by a significantly more tectonized upper crust. Indeed, the layered units may be tilted due to normal faulting. The vertical velocity gradient plot (Fig. 8c) reveals the fine-scale structures of the upper crust and highlights the complexity along this seismic section. The Layer 2A high-velocity gradient region is 400–500 m thick but more heterogeneous than that along profile 14 further south. The base of the high-velocity gradient zone appears to be highly tectonized, with both heterogeneity and rapid changes in depth, especially beneath the LL and the surrounding volcanic cones between −2.5 and 1.5 km distance across axis. Interestingly, between −1.5 and 1 km, the base of Layer 2A is offset and slightly pulled down suggesting graben subsidence. Beneath the LL, a narrow, vertical positive velocity anomaly, ~400 m wide and +0.3 to 0.5 km s⁻¹, extends from the bottom of the model to the base of the layered section of Layer 2A (Fig. 8b). This pipe-like feature
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Misfit evolution along seismic line 39 for elastic and acoustic full waveform inversions. (a) Bathymetry along the seismic line. (b) Misfit for two selected model iterations and signal size. (c) Correlation coefficient between the observed and synthetic data. WVR, western volcanic ridge; ER, eastern ridge; CG, central graben.</figcaption></figure>

could correspond to a zone of intense dyke injection that penetrates the base of Layer 2A and feeds the upper crust with fresh magma from the underlying magma chamber. A zone of intense dyke injection could explain the existence of a porosity boundary between the relatively low-porosity LL and the surrounding high-porosity volcanic cones as suggested by Arnulf et al. (2011). Fig. 8(d) shows that, unlike line 14, upper Layer 2A velocities along seismic line 24 are on average close to the reference aging curve, except beneath the

various constructional features and the footwall of the EBF where lower velocities are present. As along line 14, the upper part of Layer 2B is characterized by a smaller velocity gradient (\<1.5 s⁻¹) and by velocities that vary between 4.5 and 5.5 km s⁻¹. These velocities for the upper part of Layer 2B are similar to those of Seher et al. (2010b) from an OBS tomography study and are, once again, slowest beneath the central volcano and the fastest in the vicinity of the main normal faults bounding the axial valley (Fig. 8d).
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Upper crustal velocity structure along line 14. (a) Velocity section. (b) Velocity anomaly section created by subtracting an average 1-D velocity profile. (c) Vertical velocity gradient. The thick black line and the thick dashed lines mark, respectively, the bottom and the top of the Layer 2A high-velocity gradient zone determined from Fig. 7(c). The dotted black line marks the limit between typical Layer 2A velocities (>2.5 km s⁻¹) and anomalously low Layer 2A velocities, which might be related to the Layer 2Aa reflector observed by Arnulf et al. (2011). S, southern summital volcanic edifice. (d) Across axis average Layer 2A (red line) and Layer 2B (blue line) velocity profiles. The Layer 2A profile was computed between 50 and 250 m bsf, while the Layer 2B profile was computed over the upper first 150 m of Layer 2B. The dashed black line is a reference velocity/aging rate curve (Arnulf et al. 2013). Yellow arrows show the extent of some inferred faults observed on the bathymetry (Fig. 1).</figcaption></figure>

## 4.3 MCS line 39: the nascent graben

Line 39 is the northernmost seismic line of the SISMOMAR MCS experiment. It crosses the ridge axis ~1.8 km north of the LL centre. This line crosses the shallowest part of the WVR and the nascent graben present in this region (Fig. 1). As line 39 was previously discussed by Arnulf et al. (2012), we will only describe it briefly here for comparison purposes.

As with the two previous lines, the entire upper crust beneath the central volcano is relatively low velocity and is separated from the underlying median valley crust by the Layer 2Aa reflector (Fig. 9), suggesting that it is relatively young (Arnulf et al. 2011, 2013). The velocity structure has unusually low velocities (<2.5 km s⁻¹) in the upper 300–450 m beneath the WVR and in the upper 200–250 m beneath the rest of the central volcano. The Layer 2Aa reflector (dotted line, Fig. 9) has sudden downward shifts of ~250 m beneath the nascent graben identified from surface faulting. The offsets are larger than those seen in the surface faulting, indicating that an earlier graben has been filled by later volcanic sequences. Similar to line 24, there is no anomalously low-velocity (<2.5 km s⁻¹) layer in the upper part of the median valley and upper Layer 2A velocities seem to follow the reference aging curve (Fig. 9d), which suggests that the northern half of the central volcano is less magmatically active than the southern one. As magmatic activity wanes, tectonic process seems to play an increasingly important role moving

northwards. The vertical velocity gradient section shows several 50- to 150-m-thick layered units in upper Layer 2A (Fig. 9c). The nature of these sequences beneath the WVR is different from those in the rest of the central volcano: conical rather than flat layering. Thus, it can be inferred from the velocity structure that the WVR is young, anomalous and is of constructional origin while the subsiding graben and the ER are older and have been faulted. The ER still has anomalously low velocities (Fig. 9b), but as with the old RC on line 24, this could persist even after it is rifted away from the axis. The faults mapped on the bathymetry image (Fig. 1) seem to be imaged at depth on the vertical velocity gradient section, particularly in the vicinity of the nascent graben walls. Finally, the high-velocity gradient zone at the base of Layer 2A is ~400 m thick and appears more homogeneous than that beneath profile 24.

# 5 DISCUSSION

## 5.1 Comparison with other seismic methods

As mentioned, the nature of Layer 2A at Lucky Strike has been studied using three different methods: traveltime tomography of OBS data (Seher et al. 2010a,b), stacking of turning rays using move-out velocity (Combier 2007; Seher et al. 2010c) and more recently
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Upper crustal velocity structure along line 24. Same as Fig. 7. NE, northeastern; NW, northwestern summit volcanic edifices; LL, lava lake; RC, rifted volcanic cone.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Upper crustal velocity structure along line 39. Same as Fig. 7. WVR, western volcanic ridge; ER, eastern ridge; CG, central graben.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Comparison of velocity structure along seismic line 39, from FWI, OBS and SOBE tomography. (a) OBS traveltime tomography model (Seher et al. 2010b). (b) SOBE traveltime tomography model (Arnulf et al. 2011, 2013). (c) Full waveform inversion model (this study). The corresponding vertical velocity gradient structures are presented in (d), (e) and (f), respectively. WVR, western volcanic ridge; ER, eastern ridge; CG, central graben.</figcaption></figure>

using the SOBE method (Arnulf *et al.* 2011, 2013). Below we discuss the results obtained using these three methods and compare them with the new FWI inversion results.

### 5.1.1 FWI versus traveltime tomography

Seher *et al.* (2010b) presented an OBS-based traveltime tomography model of Lucky Strike volcano (Figs 10a and d) derived from 3-D MCS shots and 18 OBSs spaced 5 km apart. The model resolved an off-axis velocity increase (~1 km s⁻¹), a low-velocity region within the median valley and a low-velocity anomaly underneath the central volcano (Fig. 10a). The vertical velocity gradient within Layer 2A was roughly constant at 2–2.5 s⁻¹ (Fig. 10d) and the lateral resolution of the velocity structure is ~5 km, the OBS spacing. The OBS tomographic model was shown to be moderately well constrained within the upper central volcano, with likely vertical smearing of structure to greater depths. Arnulf *et al.* (2011, 2013) presented a second tomography model of Lucky Strike volcano using an SOBE-based multichannel streamer tomography technique to obtain higher resolution (Figs 10b and e). The SOBE results revealed the strong lateral heterogeneity of the upper crust and indicated that the upper low-velocity anomaly is confined in the topmost ~500 m of the crust (Fig. 10b). An ~300-m-thick region of high vertical velocity gradient (>3.5 s⁻¹) was observed at the base of Layer 2A, as well as secondary shallower zones of high gradient (Fig. 10e). The SOBE-based tomography produced velocity models with horizontal and vertical resolutions of hundreds of metres, which are beyond the reach of conventional OBSs experiments with deployments on a 5 km grid. While the SOBE method improves the velocity information on the upper crust, the FWI procedure improves the gradient information (Figs 10c and f). The high vertical velocity gradient

zone that lies at the base of Layer 2A is sharpened and several ~50- to 150-m-thick layered features that may correspond to piled-up lava sequences are unambiguously resolved in the upper part of Layer 2A (Figs 10c and f). The FWI inversion of downward continued streamer data thus provides velocity structures with spatial resolutions of 50–100 m.

The explanation for the different resolutions comes from the intrinsic nature of these techniques and from the acquisition geometry employed in these experiments. The OBS-based traveltime tomography overestimates velocities at shallow depths due to a combination of the lack of turning rays in this region and the smoothness regularization operator used to stabilize traveltime tomography, which minimizes the gradient and the curvature of the velocity model. On the other hand, application of the SOBE method to MCS data takes advantage of the fine spatial sampling of the MCS profiling while at the same time ensuring dense sampling of shallowly turning energy. With this new source receiver geometry, the resolution of the upper crust is greatly improved. Nevertheless, the smoothness constraints included in the regularization operator and the limitations of traveltime picking still limits the resolution of the final tomography model. In contrast, the FWI procedure uses the full seismic wavefield in the generation of a high-resolution velocity model and does not require the application of a regularization operator. To conclude, the FWI provides crustal structure at spatial scales approaching 50–100 m, approximately one order of magnitude higher resolution than the SOBE-based tomography method and roughly two orders of magnitude better resolution than the OBS-based tomography method. The resolution of the waveform inversion could be further improved by inverting near offset data (Shipp & Singh 2002), which can provide resolution on the scale of the quarter of wavelength and spatial resolution <50–100 m. Therefore the FWI procedure is an attractive method that can
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Comparison of full waveform inversion model for seismic line 39 with stacked sections. (a) Velocity structure and (b) vertical velocity gradient section in the time/distance domain. (c) Stacked section created from a set of synthetically calculated seismograms modelled at a flat datum 1450 m below sea surface using the FWI model. (d) Stacked section created using the data from line 39 downwards continued to the same datum. Black (a, b) and red dashed lines (c, d) mark the top and the bottom of the Layer 2A high gradient region (b). Labels as for previous figures.</figcaption></figure>

image the velocity structure of the oceanic crust with unprecedented
resolution.

5.1.2 FWI versus stacking of turnings rays versus reverse-time migration (RTM)

Layer 2A at Lucky Strike volcano has been studied by stacking wide-
angle reflections/refractions (Singh et al. 2006; Combier 2007;
Seher et al. 2010c; Arnulf et al. 2011). In order to compare the
FWI results with the stacked image, we converted the FWI model
along seismic line 39 from depth to the time domain (Figs 11a and
b). The first stacked section was created from a set of synthetically
calculated seismograms modelled at a flat datum 1450 m below sea
surface using our best-fit FWI velocity model (Fig. 11c), while the
second stack section used the observed seismic data downwards
continued to the same datum (Fig. 11d). We used a constant normal
moveout velocity of 1900 m s⁻¹ to stack the Layer 2A event in
both sections. The observed data were bandpass filtered during the
downward extrapolation stage to match the frequency spectrum of
the modelled data, allowing a direct comparison between the two
stack sections. The initial impression is that the two sections are
very similar (Figs 11c and d). The bright seismic reflector at the
base of Layer 2A in the data stack is accurately reproduced in the
synthetic stack along with some of the shallow crustal reflectors
observed in the upper Layer 2A. The superimposed top and bot-
tom of the high-velocity gradient zone indicates that the Layer 2A
reflection in the stacked seismic image corresponds, in general, to
the bottom of the high gradient zone. The excellent correspondence
between these two stacks is a demonstration of the quality of the
inversion result. However, to appreciate the strength of the FWI

procedure, one can note that the FWI velocity model contains more
fine-scale structures in the upper crust than are visible in the stacked
section (Figs 11b and d). A closer comparison between the Layer
2A reflector in the stack section with the Layer 2A high gradient
region reveals some significant discrepancies, especially in places
where the base of Layer 2A has a significant dip. The main cause
of these differences is the large lateral velocity gradients, which vi-
olate the assumptions underpinning the seismic stacking approach.
In order to address these discrepancies, we performed a finite-
difference pre-stack RTM (Fig. 12). RTM is a sophisticated pre-
stack wave-equation migration for accurate imaging in and below
areas with great structural and velocity complexity (e.g. Baysal
et al. 1983). RTM discretizes the image volume into a mesh grid
and computes full two-way numerical solutions to the wave equa-
tion at each mesh point of the finite difference grid. As a valuable
by-product, it has no dip limitation and it handles extreme lateral
velocity variations using all possible arrivals. However, there is a
price for its superior imaging capabilities: RTM can be orders of
magnitude more computationally and memory intensive than more
routine imaging methods.

Our RTM migration algorithm implements the 2-D isotropic elas-
tic wave equation described by Levander (1988). To keep the com-
putational requirement to manageable level the maximum migration
frequency was set to ~40 Hz and the velocity model was discretized
on a 6.25 m by 6.25 m grid. A total of 356 surface shot gathers,
spaced every 12.5 m, were migrated along seismic line 39. For each
shot, a virtual source wavefield was propagated forward in time
and the recorded receiver wavefield was propagated backwards in
time. To produce the migrated section these two wavefields were
correlated at each time step then summed over all times and shots
(Fig. S4). We chose to mute the energy associated to the Layer
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Reverse-time migration (RTM) imaging. (a) Migrated image, (b) composite migrated/velocity image and (c) composite migrated/gradient image along line 39. The solid red (or white) lines and the dashed red (or white) lines mark, respectively, the bottom and the top of the Layer 2A high-velocity gradient zone. The dotted red (or white) lines mark the limit between typical Layer 2A velocities (>2.5 km s⁻¹) and anomalously low Layer 2A velocities, which we suggest is related to the Layer 2Aa reflector observed by Arnulf *et al.* (2011). Black lines highlight faults. Labels as for Fig. 9.</figcaption></figure>

2A/2B turning rays prior to migration to ensure the best execution of the RTM imaging (Fig. S4b).

Fig. 12(a), shows the migrated image along seismic line 39 with resolution of ~50 m. The upper part of Layer 2A appears to be composed of multitude of volcanic sequences that, in places, appear to be truncated and offset by normal faults (black lines, Fig. 12). Within the high-velocity gradient region at the base of Layer 2A, these shallow reflectors gradually disappear, suggesting some form of geological transition, perhaps between pillow basalts and the sheeted dyke complex, although the apparent reduction in velocity heterogeneity could also be a consequence of hydrothermal sealing of cracks. The composite migrated/velocity images (Figs 12b and c) show an excellent match between the Layer 2Aa reflector (Arnulf *et al.* 2011) and the base of the anomalously low-velocity region present within the upper section of the central volcano. One interpretation of this horizon is as a palaeoseafloor juxtaposing younger, less altered volcanics above with older, more altered volcanics below (Arnulf *et al.* 2013).

Fine-scale imaging of the upper oceanic crust requires imaging techniques that can take into account the great structural and velocity complexity of this region. We have demonstrated that FWI in conjunction with RTM can produce velocity and reflection images superior to more conventional methods. The principal advantage of FWI combined with downward continuation is that it unifies in a single velocity model the information from both refractions and wide-angle reflections. As illustrated for line 39, Fig. 11, a wide-angle stack of the Layer 2A reflection can be displaced relative to the base of Layer 2A defined via velocity gradient, perhaps leading

to discrepancies in different interpretations of Layer 2A (Christenson *et al.* 2007). Another advantage of FWI is that it yields absolute values of crustal elastic parameters (Tarantola 1984) and, provided that the errors due to the neglect of attenuation and anisotropy are not too severe, these can be used to better quantify and compare geological processes between ridge settings.

### 5.1.3 Elastic versus acoustic FWI

The first arrivals in the downward continued data are mostly transmitted compressional energy, which are mainly sensitive to P-wave speed. Therefore, one can question the need for fully elastic wave-form inversion in preference to an acoustic method (Canales 2010). To assess the difference, we have carried out FWI assuming that S-wave velocity is zero, that is, the medium is acoustic and compare the results (Fig. 13) and misfit functions (Fig. 6). Although the velocity structures from the elastic and acoustic inversions are similar (Fig. 13), the difference in the inverted velocities can be significant (±400 m s⁻¹), particularly near the seafloor (Fig. 13e), which would have a significant effect on porosity estimates. The overall misfit function using the elastic inversion is slightly better than that using acoustic inversion, but for certain shot points the improvements can be 10–20 per cent (Fig. 6). Even the cross-correlation between synthetic and real data shows improvement using the elastic inversion (0.8 compared to 0.72; Fig. 6c). The reason for this improvement could be the correct partitioning of seismic energy at the seafloor, a factor even for turning rays, and amplitude versus offset variations
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Elastic versus acoustic full waveform inversion. (a) Elastic FWI velocity model and its corresponding vertical velocity gradient structure (c). (b) Acoustic FWI velocity model and its corresponding gradient structure (d). (e) Difference between the elastic and acoustic FWI velocity models. Labels as for Fig. 9.</figcaption></figure>

in wide-angle reflections, included in the inversion window after the first arrivals. If the velocity contrast on the seafloor were much larger, then the improvement using elastic inversion would be much more drastic (Sears *et al.* 2008). Since the velocity increase at the seafloor is low (from 1.5 to 1.7 km s⁻¹), the effect of elasticity is rather moderate here.

## 5.2 Crustal aging and origin of Layer 2A

### 5.2.1 Implications of the velocity distribution along seismic line 39

In Fig. 14, we present the 1-D profiles of the waveform inversion results for five specific regions along line 39. The 1-D velocity models were created by picking the most common velocity from velocity histograms computed at each depth within the bounded areas (Fig. 14a). While three of these 1-D velocity profiles were extracted from regions of moderate faulting (regions bracketed by black lines, Fig. 14), two 1-D velocity profiles were computed in regions of extensive faulting (regions bracketed by pink lines, Fig. 14). The 1-D velocity profiles from regions of moderate faulting as well as the overall probability density function, exhibit a fair degree of heterogeneity in the upper part of Layer 2A but velocities at the base of Layer 2A are more homogeneous and the velocity at the top of Layer 2B is tightly clustered around 4.8 ± 0.15 km s⁻¹. Nevertheless, the distribution of Layer 2A velocities is not randomly distributed within the upper crust. Indeed, the lowest P-wave velocities are observed beneath the recently constructed WVR (Fig. 14, solid black line), while the highest P-wave velocities are observed beneath the median valley on the western side of the central volcano

(Fig. 14, dashed dotted black line). Within the central volcano but away from the WVR, the velocities follow an intermediate trend lying between the two endmembers (Fig. 14, dashed black line). It is widely accepted that velocity changes within Layer 2A are caused by changes in porosities (Becker *et al.* 1982; Carlson & Herrick 1990; Jacobson 1992). The velocity variations observed along seismic line 39 are consistent with the idea that hydrothermal mineralization seals the inherent small-scale porosity of the upper oceanic crust with age. However, the velocity increase is not uniform, being much greater within the upper part of Layer 2A, top ~300 m, than at depth within the transition zone, and indeed velocities at the base of Layer 2A appear anchored around 4.8 km s⁻¹. This behaviour would be consistent with hydrothermal mineralization preferentially sealing thin cracks within upper part of Layer 2A producing a large velocity increase (Wilkens *et al.* 1991). There is evidence that the velocity increase just below the seafloor, top 50–100 m, lags behind the one at depth, perhaps indicating a hydrothermal front progressing upwards through Layer 2A (Peterson *et al.* 1986). On these timescales, hydrothermal mineralization appears to be a secondary process that can explain the increase in upper crustal velocity with age but is not responsible for the high-velocity gradient at the base of Layer 2A.

At Lucky Strike, our results provide new evidence that an alteration front cannot be the sole cause of the velocity transition at the base of Layer 2A. In order to investigate further the nature of Layer 2A transition, we analysed the velocities at locations straddling two major fault zones (see RTM image, Fig. 12; regions bracketed by pink lines on Fig. 14a). At both, the top of the high gradient zone starts only 200 m bsf and the base of Layer 2A is anomalously shallow lying only 0.4–0.5 km bsf, Fig. 14(b). In addition, this pull up of the high gradient zone is abrupt at the fault zone, for example,
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Velocity distribution along line 39. (a) Full waveform inversion velocity structure. The five-bracketed areas delimit the regions where the corresponding 1-D velocity profiles in (b) were extracted. (b) Velocity probability density function for line 39. The white diamonds mark the base of Layer 2A on the individual 1-D profiles. The grey region bounds the velocities $4.8 \pm 0.15 \text{ km s}^{-1}$, which to first order indicates the base of Layer 2A. The inset cartoon illustrates how tectonic thinning may pull up the high-velocity gradient region at the base of Layer 2A. Labels as for Fig. 9.</figcaption></figure>

Fig. 12. This observation strongly suggests the tectonic thinning of a geologically defined layer (see sketch, Fig. 14b), and that the base of Layer 2A can be treated as a lithological boundary to infer fault offsets, even if the nature of the boundary is unresolved. Therefore, the nature of this lithological boundary is still open to debate. The Layer 2A high-velocity gradient zone could correspond to the upper part of the sheeted dyke transition rather than the whole transition region. If so the upper part of 2A would be almost exclusively pillow lavas, alternatively the Layer 2A transition could just span the interval of the first appearance of dykes and not the whole dyke transition (Karson *et al.* 2002).

We estimated the porosity in Layer 2A from its seismic velocity structure using a differential effective medium theory (DEMT; Berryman 1995; Taylor & Singh 2002; Fig. 15). A two-phase effective medium consisting of basalt and sea water was assumed (Seher *et al.* 2010b; Arnulf *et al.* 2011). Subsurface porosities were computed as function of velocities for different vertically aligned crack distributions. Cracks at mid-ocean ridges include region between individual lava flows, pore spaces within pillows, space between sheeted dykes, cooling cracks within individual sequences, faults and microcracks (Moos & Marion 1994; Collier & Singh 1998). Using 1-D FWI from data from southern East Pacific Rise, Collier & Singh (1998) have shown that porosity at the base of Layer 2A is ~6 per cent. Our results show that the presence of exclusively spherical inclusions (aspect ratio 1, Fig. 15e) requires unrealistically high porosity (>50 per cent) in the whole Layer 2A. Cracks with aspect ratios of 5–10 better explain the observed velocities but still have some unrealistically high porosities either in the

upper 300–400 m of crust or at the base of Layer 2A (Figs 15c and d). Cracks with aspect ratio >20 (Figs 15a and b) give a porosity of ~5 per cent at the base of Layer 2A, consistent with previous estimates (Collier & Singh 1998). The maximum porosity within the young constructional WVR is between 30 and 35 per cent, while typical porosity in the upper 200–300 m is around 10–20 per cent, consistent with upper-crustal porosity measurements from on-bottom gravimeters (25–30 per cent; Cochran *et al.* 1999; Gilbert & Johnson 1999; Gilbert *et al.* 2007). A small amount of porosity, mainly due to cracks with aspect ratios >20, is thus sufficient to explain the velocities observed across the Lucky Strike volcano. Other models, which consider a mixture of cracks, have further shown that much of the increase in velocity seen in upper levels of the oceanic crust is likely controlled by the sealing of thin cracks with aspects ratios >100 (Wilkens *et al.* 1991). Although a significant portion of this porosity may be explained by magmatic processes (boundaries between individual lava flows and pillows and their broken fragments, boundaries between the sheeted dykes, cooling cracks within individual sequences), a part of porosity is likely attributable to tectonic extension (faults, fissures).

### 5.2.2 Across axis velocity variations at the centre of the Lucky Strike segment

Fig. 16 shows the velocity models along lines 14, 24 and 39 plotted in terms of depth below the seafloor, highlighting a number of different features. Young constructional features such as the southern
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Estimates of subsurface porosity along line 39, using a vertically aligned crack distribution and differential effective medium theory. Subsurface porosity distribution using vertically aligned cracks with aspect ratios of 1/50 (a), 1/20 (b), 1/10 (c), 1/5 (d) and 1 = spherical inclusions (e) are represented. The thin black lines mark the porosity contours 10, 20, 30 and 40 per cent, while the thin white line marks the 50 per cent porosity contour. Labels as for Fig. 9.</figcaption></figure>

summital volcanic edifice on line 14 (Fig. 16a) or the WVR on line 39 (Fig. 16c) show up as notable depressions in the velocity with deep roots, suggestive of constructional origin. The ER along line 39 also seems to be constructional in origin but has higher velocity and is thus likely a relatively older feature (Fig. 16c). In contrast to the southern cone, the northwestern and northeastern volcanic cones (Fig. 16b) seem to be considerably older features. Unusually low velocities (<2.5 km s⁻¹) are observed at 400–450 m beneath the southern summital volcanic cone and the WVR, suggesting that constructional features can support a thick pile of high-porosity extrusives without collapsing. However, pore closure due to overburden pressure is one of the various suggestions for the cause of high-velocity gradient at the base of Layer 2A. So, while in general there is pore closure with depth, velocity contours are not horizontal with respect to the seafloor, and thus the lateral variations in velocity must have another cause. As discussed previously, the most obvious signature of crustal aging are the lateral increase of velocity in the

upper crust. Along the seismic lines 14 and 39 (see Fig. 16) velocities in the upper part of Layer 2A above 450 m increase rapidly from 2.5 km s⁻¹ beneath the young constructional seamounts to 3.5 km s⁻¹ within the median valley, suggesting hydrothermal filling of pore spaces with alteration products. Moreover, velocities in the upper crust along lines 24 and 39 are globally faster than those along line 14 on the southern half of the Lucky Strike volcano, Fig. 16(d). This observation could signify either an enhanced magma supply within the southern half of the Lucky Strike volcano, different types of lava units or a more intense hydrothermal circulation within the northern part of the Lucky Strike volcano, accompanied by enhanced hydrothermal mineralization. These inferences accord well with quantitative seafloor observations over the Lucky Strike segment, which indicate that the volcanic terrains to the south of the central volcano were erupted and subsequently faulted more recently (Le Saout *et al.* 2011). While the southern half of the Lucky Strike volcano was shown to be more magmatic, its northern counterpart is clearly more tectonized, which would favour hydrothermal mineralization. In other words, hydrothermal circulation might cause the upper crust to age faster within the northern part of the Lucky Strike volcano than in its southern counterpart.

## 6 CONCLUSION

We have presented three 2-D high-resolution velocity profiles of the upper crust beneath the Lucky Strike volcano using a new methodology that combines a downward extrapolation method and a 2-D elastic FWI procedure. We were able to clearly define two layers within Layer 2A, upper Layer 2A, mainly consisting of recent heterogeneous thin lava flows, and a homogeneous high-velocity gradient zone at the base of Layer 2A. These results helped us to understand the role of magmatic, tectonic and hydrothermal processes in defining the upper crustal structure of the Lucky Strike segment. We favour a model, where the base of Layer 2A is defined as a lithological boundary that can be offset by faulting. However, the nature of this lithological boundary is still open to debate, the Layer 2A high-velocity gradient zone could correspond to the upper part of the sheeted dyke transition rather than the whole transition region. We show that young volcanic cones can support a thick pile of high-porosity extrusive without collapsing, which suggest that pore closure due to overburden pressure with depth is not the cause for the high-velocity gradient at the base of Layer 2A. We demonstrate that hydrothermal sealing of small-scale porosity is a secondary process, which may explain the upper crustal velocity increase with age but in no instance is responsible for the high-velocity gradient layer present at the base of Layer 2A, in which little variations are noted across axis. Finally, we show that the main normal faults can fracture the entire lithologically defined Layer 2A.

## ACKNOWLEDGEMENTS

This research was carried out at the Institut de Physique du Globe de Paris and at the Scripps Institution of Oceanography. The National Science Foundation supported this work via grant OCE-0826481. The Cecil H. and Ida M. Green Foundation for Earth Sciences also supported part of this research. We are grateful to the captain, crew and seismic team from GENAVIR of cruise SISMOMAR of the RV L'Atalante for their assistance. The INSU MOMAR program funded the acquisition of SISMOMAR project.
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Velocity structure along seismic lines 14 (a), 24 (b) and 39 (c) plotted in terms of depth below seafloor. The corresponding bathymetry profile is included in each subplot. The thick dashed black lines mark the bottom and the top of the Layer 2A high gradient zone. Velocities are contoured at 0.5 km s⁻¹ intervals between 2.0 and 5.5 km s⁻¹ (grey lines). (d) Comparison of the average velocity structure along seismic lines 14, 24 and 39. Labels as for Fig. 9.</figcaption></figure>

## REFERENCES

Arnulf, A.F., Singh, S.C., Harding, A.J., Kent, G.M. & Crawford, W., 2011. Strong seismic heterogeneity in Layer 2A near hydrothermal vents at the Mid-Atlantic Ridge, *Geophys. Res. Lett.*, **38**, L13320, doi:10.1029/2011GL047753.

Arnulf, A.F., Harding, A.J., Singh, S.C., Kent, G.M. & Crawford, W., 2012. Fine-scale velocity structure of upper oceanic crust from full waveform inversion of downward continued seismic reflection data at the Lucky Strike Volcano, Mid-Atlantic Ridge, *Geophys. Res. Lett.*, **39**, L08303, doi:10.1029/2012GL051064.

Arnulf, A.F., Harding, A.J., Singh, S.C., Kent, G.M. & Crawford, W., 2013. Constraints on the shallow velocity structure of the Lucky Strike Volcano, Mid-Atlantic Ridge, from downward continued multi-channel streamer data, *J. geophys. Res.*, in press.

Baysal, E., Kosloff, D.D. & Sherwood, J.W.C., 1983. Reverse time migration, *Geophysics*, **48**, 1514–1524.

Becker, K. *et al.*, 1982. In situ electrical resistivity and bulk porosity of the oceanic crust Costa Rica Rift, *Nature*, **300**, 594–598.

Berryhill, J.R., 1984. Wave-equation datuming before stack, *Geophysics*, **49**, 2064–2066.

Berryman, J.G., 1995. Mixture theories for rock properties, in *AGU Handbook of Physical Constants*, pp. 205–228, ed. Ahrens, T.J., American Geophysical Union.

Brenders, A.J. & Pratt, R.G., 2007. Full waveform tomography for lithospheric imaging: results from a blind test in a realistic crustal model, *Geophys. J. Int.*, **168**, 133–151.

Burnett, M.S., Caress, D.W. & Orcutt, J.A., 1989. Tomographic image of the magma chamber at 12°50′N on the East Pacific Rise, *Nature*, **339**, 206–208.
Canales, J.P., 2010. Small-scale structure of the Kane oceanic core complex, Mid-Atlantic Ridge 23°30′N, from waveform tomography of multichannel seismic data, *Geophys. Res. Lett.*, **372**, L21305, doi:10.1029/2010GL044412.

Canales, J.P., Detrick, R.S., Lin, J. & Collins, J.A., 2000. Crustal and upper mantle seismic structure beneath the rift mountains and across a non-transform offset at the Mid-Atlantic Ridge (35°N), *J. geophys. Res.*, **105**, 2699–2719.

Carlson, R.L., 2004. Seismic properties of Layer 2A at 11 Ma: results of a vertical seismic profile at Ocean Program Site 1243, *Geophys. Res. Lett.*, **31**, L17601, doi:10.1029/2004GL020598.

Carlson, R.L. & Herrick, C.N., 1990. Densities and porosities in the oceanic crust and their variations with depth and age, *J. geophys. Res.*, **95**, 9153–9170.

Castagna, J.P., Batzle, M.L. & Eastwood, R.L., 1985. Relationships between compressional-wave and shear-wave velocities in clastic silicate rocks, *Geophysics*, **50**, 571–581.

Christeson, G., Purdy, G. & Fryer, G., 1994. Seismic constraints on shallow crustal emplacement processes at the fast spreading East Pacific Rise, *J. geophys. Res.*, **99**, 17 957–17 973.

Christeson, G., Shaw, P.R. & Garmany, J.D., 1997. Shear and compressional wave structure of the East Pacific Rise, 9°-10°N, *J. geophys. Res.*, **102**, 7821–7836.

Christeson, G.L., McIntosh, K.D. & Karson, J.A., 2007. Inconsistent correlation of seismic layer 2a and lava layer thickness in oceanic crust, *Nature*, **445**, 418–421.

Cochran, R., Fornari, D.J., Coakley, B.J., Herr, R. & Tivey, M.A., 1999. Continuous near-bottom gravity measurements made with a BGM-3 gravimeter in DSV Alvin on the East Pacific Rise crest near 9°40′N and 9°50′N, *J. geophys. Res.*, **104**, 10 841–10 861.

Collier, J.S. & Singh, S.C., 1997. Detailed structure of the top of the melt body beneath the East Pacific Rise at 9°40′N from waveform inversion of seismic reflection data, *J. geophys. Res.*, **1022**, 20 287–20 304.

Collier, J.S. & Singh, S.C., 1998. Poisson's ratio structure of young oceanic crust, *J. geophys. Res.*, **1032**, 20 981–20 996.

Combier, V., 2007. Mid-ocean ridge processes. Insights from 3D refraction seismics at 9°N OSC on the East Pacific Rise, and the Lucky Strike Volcano on the Mid-Atlantic Ridges, *PhD thesis, Laboratoire de Geosciences Marines*, Institut de Physique du Globe de Paris, Paris, France.

Cudrak, C.F. & Clowes, R.M., 1993. Crustal structure of Endeavour Ridge segment, Juan de Fuca Ridge, from a detailed seismic refraction survey, *J. geophys. Res.*, **98**(B4), 6329–6349.

Demets, C., Gordon, R.G., Argus, D.F. & Stein, S., 1994. Effect of recent revisions to the geomagnetic reversal timescale on estimates of current plate motions, *Geophys. Res. Lett.*, **21**, 2191–2194.

Fliedner, M.M., White, R.S. & Smallwood, J.R., 1998. Seismic velocity structure of basalt flows, *SEG Expanded Abstracts*, SP4.5, 1178–1181, http://dx.doi.org/10.1190/SEGEAB.17.

Fouquet, Y. et al., 1994. A detailed study of the Lucky Strike hydrothermal site and discovery of a new hydrothermal site: Menez Gwen: preliminary results of the DIVA1 cruise (5–29 May, 1994), *InterRidge News*, **3**(2), 14–17.

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

Gilbert, L.A. & Johnson, H.P., 1999. Direct measurements of oceanic crustal density at the northern Juan de Fuca Ridge, *Geophys. Res. Lett.*, **26**, 3633–3636.

Gilbert, L.A., McDuff, R.E. & Johnson, H.P., 2007. Porosity of the upper edifice axial seamount, *Geology*, **35**, 49–52.

Hamilton, E.L., 1978. Sound velocity-density relations in sea-floor sediments and rocks, *Acoust. Soc. Am. J.*, **63**, 366–377.

Harding, A.J., Kent, G.M. & Orcutt, J.A., 1993. A multichannel seismic investigation of the upper crustal structure at 9°N on the East Pacific Rise: implications for crustal accretion, *J. geophys. Res.*, **98**, 13 925–13 944.

Humphris, S.E., Fornari, D.J., Scheirer, D.S., German, C.R. & Parson, L.M., 2002. Geotectonic setting of hydrothermal activity on the summit of

Lucky Strike Seamount (37°17′N, Mid-Atlantic Ridge), *Geochem. Geo-phys. Geosyst.*, **3**, doi:10.1029/2001GC000284.

Hussenoeder, S.A., Kent, G.M. & Detrick, R.S., 2002a. Upper crustal seismic structure of the slow spreading Mid-Atlantic Ridge, 35°N: constraints on volcanic emplacement processes, *J. geophys. Res.*, **107**(B8), doi:10.1029/2001JB001691.

Hussenoeder, S.A., Detrick, R.S., Kent, G.M., Schouten, H. & Harding, A.J., 2002b. Fine-scale seismic structure of young upper crust at 17°20′S on the fast spreading East Pacific Rise, *J. geophys. Res.*, **107**(B8), doi:10.1029/2001JB001688.

Jacobson, R.S., 1992. The impact of crustal evolution on changes of the seismic properties of the uppermost ocean crust, *Rev. Geophys.*, **30**, 23–42.

Karson, J.A., Tivey, M.A. & Delaney, J.R., 2002. Internal structure of uppermost oceanic crust along the Western Blanco Transform Scarp: implications for subaxial accretion and deformation at the Juan de Fuca Ridge, *J. geophys. Res.*, **107**(B9), doi:10.1029/2000JB000051.

Langmuir, C. et al., 1997. Hydrothermal vents near a mantle hot spot: the Lucky Strike vent ffield at 37°N on the Mid-Atlantic Ridge, *Earth planet. Sci. Lett.*, **148**, 69–91.

Larkin, S.P. & Levander, A., 1996. Wave-equation dating for improving deep crustal seismic images, *Tectonophysics*, **264**, 371–379.

Le Saout, M., Deschamps, A., Escartin, J., Cannat, M. & Carlut, J., 2011. Along-axis variations of tectono-magmatic processes at Lucky-Strike segment, Mid-Atlantic Ridge, in *Abstract T23D-2427 Presented at 2011 Fall Meeting, AGU, San Francisco, CA*, 5–9 December 2011.

Levander, A.R., 1988. Fourth-order finite-difference P-SV seismograms, *Geophysics*, **53**, 1425–1436.

Lewis, B.T.R. & Garmany, J.D., 1982. Constraints on the structure of the East Pacific Rise from seismic refraction data, *J. geophys. Res.*, **87**(B10), 8417–8425.

McClain, J.S., Orcutt, J.A. & Burnett, M., 1985. The East Pacific Rise in cross section: a seismic model, *J. geophys. Res.*, **90**, 8627–8640.

Minshull, T.A. & Singh, S.C., 1993. Shallow structure of the oceanic crust in the western North Atlantic from seismic waveform inversion and modeling, *J. geophys. Res.*, **98**, 1777–1792.

Moos, D. & Marion, D., 1994. Morphology of extrusive basalts and its relationship to seismic velocities in the shallow oceanic crust, *J. geophys. Res.*, **99**, 2985–2994.

Mora, P., 1987. Nonlinear two-dimensional elastic inversion of multiooffset seismic data, *Geophysics*, **52**, 1211–1228.

Moreira, M. & Allegre, C.J., 2002. Rare gas systematics on Mid-Atlantic Ridge (37–40°N), *Earth planet. Sci. Lett.*, **198**, 401–416.

Moreira, M., Escartin, J., Gayer, E., Hamelin, C., Bezos, A., Guillon, F. & Cannat, M., 2011. Rare gas systematics on Lucky Strike basalts (37°N, North Atlantic): evidence for efficient homogenization in a long-lived magma chamber system? *Geophys. Res. Lett.*, **380**, L08304, doi:10.1029/2011GL046794.

Neves, F.A. & Singh, S.C., 1996. Sensitivity study of seismic reflection/refraction data, *Geophys. J. Int.*, **126**, 470–476.

Nicolas, A., Reuber, I. & Benn, K., 1988. A new magma chamber model based on structural studies in the Oman ophiolite, *Tectonophysics*, **151**, 87–105.

Ondréas, H., Fouquet, Y., Voisset, M. & Radford-Knoery, J., 1997. Detailed study of three contiguous segments of the Mid-Atlantic Ridge, south of the Azores (37°N to 38°30′N), using acoustic imaging coupled with submersible observations, *Mar. geophys. Res.*, **19**, 231–255.

Ondréas, H., Cannat, M., Fouquet, Y., Normand, A., Sarradin, P.M. & Sarrazin, J., 2009. Recent volcanic events and the distribution of hydrothermal venting at the Lucky Strike hydrothermal ffield, Mid-Atlantic Ridge, *Geochem. Geophys. Geosyst.*, **10**, doi:10.1029/2008GC002171.

Peterson, C., Duncan, R. & Scheidegger, K.F., 1986. Sequence and longevity of basalt alteration at deep sea drilling project site 597, in *Initial Reports of the Deep Sea Drilling Project Covering Leg 92 of the Cruises of the Drilling Vessel Glomar Challenger, Papeete, Tahiti, to Balboa, Panama, February-April, 1983*, pp. 505–515, available at: http://www.deepadrilling.org/92/volume/dsup92_32.pdf.
Pica, A., Diet, J.P. & Tarantola, A., 1990. Nonlinear inversion of seismic reflection data in a laterally invariant medium, *Geophysics*, **55**, 284–292.

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

Pratt, R.G. & Worthington, M.H., 1990. Inverse theory applied to multi-source cross-hole tomography. Part 1: acoustic wave-equation METHOD1, *Geophys. Prospect.*, **38**, 287–310.

Robertsson, J.O.A. & Chapman, C.H., 2000. An efficient method for calculating finite-difference seismograms after model alterations, *Geophysics*, **65**, 907–918.

Robertsson, J.O.A., Ryan-Grigor, S., Sayers, C.M. & Chapman, C.H., 2000. A finite-difference injection approach to modeling seismic fluid flow monitoring, *Geophysics*, **65**, 896–906.

Royle, G., 2010. Viscoelastic full waveform inversion for ocean-bottom cable seismic data, *PhD thesis*, Laboratoire de Geosciences Marines, Institut de Physique du Globe de Paris, Paris, France.

Sears, T.J., Singh, S.C. & Barton, P.J., 2008. Elastic full waveform inversion of multi-component OBC seismic data, *Geophys. Prospect.*, **56**, 843–862.

Sears, T.J., Barton, P.J. & Singh, S.C., 2010. Elastic full waveform inversion of multicomponent ocean-bottom cable seismic data: application to Alba Field, U. K. North Sea, *Geophysics*, **75**(6), 109–119.

Seher, T., Crawford, W.C., Singh, S.C., Cannat, M., Combier, V. & Dusunur, D., 2010a. Crustal velocity structure of the Lucky Strike segment of the Mid-Atlantic Ridge (37°N) from seismic refraction measurements, *J. geophys. Res.*, **115**, doi:10.1029/2009JB006650.

Seher, T., Singh, S.C., Crawford, W. & Escartin, J., 2010b. Upper crustal velocity structure beneath the central Lucky Strike Segment from seismic refraction measurements, *Geochem. Geophys. Geosyst.*, **11**, doi:10.1029/2009GC002894.

Seher, T., Crawford, W.C., Singh, S.C. & Cannat, M., 2010c. Seismic layer 2A variations in the Lucky Strike segment at the Mid-Atlantic Ridge from reflection measurements, *J. geophys. Res.*, **115**, B07107, doi:10.1029/2009JB006783.

Shipp, R.M. & Singh, S.C., 2002. Two-dimensional full wavefield inversion of wide-aperture marine seismic streamer data, *Geophys. J. Int.*, **151**(2), 325–344.

Shtivelman, V. & Canning, A., 1988. Datum correction by wave-equation extrapolation, *Geophysics*, **53**, 1311–1322.

Singh, S.C. et al., 2006. Discovery of a magma chamber and faults beneath a Mid-Atlantic Ridge hydrothermal feld, *Nature*, **442**, 1029–1032.

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

Spudich, P. & Orcutt, J., 1980. A new look at the seismic velocity structure of the oceanic crust, *Rev. geophys. Space Phys.*, **18**(3), 627–645.

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

Tarantola, A., 1986. A strategy for nonlinear elastic inversion of seismic reflection data, *Geophysics*, **51**, 1893–1903.

Taylor, M. & Singh, S., 2002. Composition and microstructure of magma bodies from effective medium theory, *Geophys. J. Int.*, **149**, 15–21.

Toomey, D.R., Purdy, G.M., Solomon, S.C. & Wilcock, W.S.D., 1990. The three-dimensional seismic velocity structure of the East Pacific Rise near latitude 9°30' N, *Nature*, **347**, 639–645.

Van Avendonk, H.J.A., Harding, A.J. & Orcutt, J.A., 1998. A two-dimensional tomographic study of the Clipperton transform fault, *J. geophys. Res.*, **103**, 17 885–17 899.

White, R.S. & Purdy, G.M., 1983. Crustal velocity structure on the flanks of the Mid-Atlantic Ridge at 24°N, *Geophys. J. R. astr. Soc.*, **75**, 361–385.

Wilcock, W.S.D., Solomon, S.C., Purdy, G.M. & Toomey, D.R., 1992. The seismic attenuation structure of a fast-spreading mid-ocean ridge, *Science*, **258**, 1470–1474.

Wilkens, R.H., Fryer, G.J. & Karsten, J., 1991. Evolution of porosity and

seismic structure of upper oceanic crust: importance of aspect ratios, *J. geophys. Res.*, **96**, 17 981–17 995.

Ziolkowski, A., 1991. Why don't we measure seismic signatures? *Geo-physics*, **56**, 190–201.

## SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this article:

**Figure S1.** Data processing via the Synthetic Ocean Bottom Experiment (SOBE) method. (a) Starting velocity model along seismic line 39 and synthetic seafloor geometry for the shot gather 38662. The black lines indicate ray paths for the first arrival crustal refraction events, the green lines mark the extent of the finite difference modelling window for this shot. (b) Raw shot gather recorded at the sea surface. (c) Same shot gather after downward extrapolation of the source and receivers to 75 m above the seafloor interface. During extrapolation the data was band pass filtered to satisfy limitations imposed by the FWI finite difference algorithm. The traces were spatially resampled to fall on finite difference grid nodes. Only refractions and reflections present within the blue triangle were inverted. (d) Normalized amplitude spectrum of the surface data (cyan line), of the downward extrapolated data (blue line) and of the initial and final synthetics, respectively the green and black lines (see Figs 3 and 4).

**Figure S2.** Offset dependence (1.75 km to 3.25 km) of the downward continued data amplitude spectrum (gray lines). With our elastic FWI inversion scheme we deliberately restricted attention to arrivals, refractions and wide-angle reflections, whose spectra exhibit only limited offset dependence (gray lines). As a result, we assumed for our elastic FWI that the principal attenuation effects can be incorporated into an effective source wavelet that is updated as part of the inversion process. Although the initial wavelet was broader band (green line), the final wavelet captures most of the spectral content of the filtered data (red line).

**Figure S3.** S-wave velocity structure recovered by our 2D elastic full waveform inversion along seismic line 39.

**Figure S4.** Constructing the reverse time migrated image. (a) Final FWI velocity model used for reverse time migration. Black lines mark the rays corresponding to the refracted wave and the seafloor reflection in a sea surface geometry for shot number 38730, located above the WVR along line 39. Green lines mark the limits of the finite difference-modelling window for this shot. (b) Surface shot gather band pass filtered to limit the frequency range for the finite-difference grid. The traces were spatially resampled to fall on a finite difference grid node. Only reflections present within the blue region were migrated. (c) Pre-stack depth migration of shot gather 38730. (d) Summation of all the migrated shot gathers. The black line marks the seafloor horizon (http://gji.oxfordjournals.org/lookup/suppl/doi:10.1093/gji/ggt461/-/DC1).

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.