Next-generation seismic experiments – II: wide-angle, multi-azimuth, 3-D, full-waveform inversion of sparse ﬁeld data

3-D full-waveform inversion (FWI) is an advanced seismic imaging technique that has been widely adopted by the oil and gas industry to obtain high-ﬁdelity models of P -wave velocity that lead to improvements in migrated images of the reservoir. Most industrial applications of 3-D FWI model the acoustic waveﬁeld, often account for the kinematic effect of anisotropy, and focus on matching the low-frequency component of the early arriving refractions that are most sensitive to P -wave velocity structure. Here, we have adopted the same approach in an application of 3-D acoustic, anisotropic FWI to an ocean-bottom-seismometer (OBS) ﬁeld data set acquired across the Endeavour oceanic spreading centre in the northeastern Paciﬁc. Starting models for P -wave velocity and anisotropy were obtained from traveltime tomography; during FWI, velocity is updated whereas anisotropy is kept ﬁxed. We demonstrate that, for the Endeavour ﬁeld data set, 3-D FWI is able to recover ﬁne-scale velocity structure with a resolution that is 2–4 times better than conventional traveltime tomography. Quality assurance procedures have been employed to monitor each step of the workﬂow; these are time consuming but critical to the development of a successful inversion strategy. Finally, a suite of checkerboard tests has been performed which shows that the full potential resolution of FWI can be obtained if we acquire a 3-D survey with a slightly denser shot and receiver spacing than is usual for an academic experiment. We anticipate that this exciting development will encourage future seismic investigations of earth science targets that would beneﬁt from the superior resolution offered by 3-D FWI.


I N T RO D U C T I O N
Practice within the petroleum sector has been transformed since the development of 3-D full-waveform inversion (3-D FWI). This technology is now commonly used to obtain fine-scale models of P-wave velocity which, when used in pre-stack or reverse-time depth migrations, lead to a significant improvement in reflection images (e.g. Plessix & Perkins 2010;Sirgue et al. 2010;Ratcliffe et al. 2011;Kapoor et al. 2012;Houbiers et al. 2013;Jones et al. 2013;Selwood et al. 2013;Warner et al. 2013). To date, 3-D FWI has been predominantly applied to industrial streamer , ocean-bottom-cable (Sirgue et al. 2010) and dense oceanbottom-node data (Bansal et al. 2013) which more-than-adequately sample the subsurface for FWI applications, meaning there is significant data redundancy and only a subset of the data set is utilized when running inversions . The relatively short shot-receiver offsets (<10 km) in such experiments, as well as the dense coverage, means that noise is not a significant issue and that these inversions can often be run with only minimal regularization (Plessix & Perkins 2010).
In a previous paper, synthetic tests were used to assess the potential of 3-D, wide-angle, long-offset, low-frequency FWI as a tool to determine P-wave velocity structure to address scientific targets (Morgan et al. 2013). Here, we apply this technology to a marine seismic data set acquired across an active ocean-spreading centre, using airgun shots and ocean-bottom seismometers (OBS). In petroleum applications the focus is on improving the migrated reflection image but, here, the interest is in recovering and then interpreting the FWI velocity model directly-in this case to better understand the accretion of oceanic crust (Arnoux et al. 2014). Our starting models for P-wave velocity and anisotropy were obtained using traveltime tomography (Weekly et al. 2014). Prior to initiating this project we recognized there would be additional challenges associated with applying 3-D FWI to an academic data set, most notably dealing with data sparsity, due to typical shot spacings of between several 10s m and a few 100 m, and ocean-bottom receiver spacings of between 1 and 20 km. In addition, we anticipated that any inversion of longer-offset data, which are required to obtain models of whole-crustal velocity, could suffer from problems with noise, in particular at the low-frequency range (<4.5 Hz). Given these issues, it was unclear as to whether FWI would be able to resolve velocity anomalies of the order of half a seismic wavelength, the theoretical resolution of FWI (Pratt et al. 1996) or even be successful at all.
With industrial data sets, we are able to check whether our inversions are successful through a suite of different quality-control procedures, including verifying that the inverted velocity model better matches sonic well logs and checking whether the migrated image is improved (e.g. Sirgue et al. 2010;Warner et al. 2013). Here, without drill holes or 3-D seismic reflection data, we have to rely on (1) verifying the improvement in match between observed and predicted data, (2) using phase plots to check for cycle skipping, and (3) using phase plots to check whether the inversion gradually converges to a global minimum . We have performed a suite of inversions on the field data to seek a successful workflow that is able to deal effectively with data sparsity and noise, and have run a suite of synthetic tests to verify that our chosen workflow is able to recover the anomalies that we see in our inverted velocity models. Finally, we ran additional checkerboard tests to explore how acquisition parameters affect resolution, and how we could improve future experiments in order to exploit the full potential resolution of FWI.
The purpose of this paper is to demonstrate that 3-D FWI, a technology developed for the petroleum industry, can be applied to sparse marine seismic data sets that are typically acquired in scientific investigations of crustal targets. In a separate paper, a synthesis of the resulting FWI velocity model is presented together with other existing data and used to investigate ridge processes at the Endeavour spreading centre (Arnoux et al. 2016), including placing constraints on the reaction zone that links magmatic and hydrothermal systems, and controls the pattern of heat transfer.

E N D E AV O U R F I E L D DATA S E T Geology
The study area is the Endeavour segment of the Juan de Fuca Ridge in the northeastern Pacific, which is bounded by two large-offset overlapping spreading centres (Fig. 1). The Endeavour segment is ∼90 km long and has a full-spreading rate of 52 mm yr −1 (DeMets et al. 2010). Multichannel seismic (MCS) studies have shown that an axial magma chamber (AMC) lies beneath the central portion of the Endeavour segment (Van Ark et al. 2007;Carbotte et al. 2008), where there are several high-temperature hydrothermal vent fields (Kelley et al. 2012). These vents have evolved geochemically with time (Lilley et al. 2003), and show along-axis gradients in temperature and chemistry (Butterfield et al. 1994). The area has an interesting recent history, with periods of northward and southward propagation at each end of the Endeavour segment (Davis & Lister 1977;Johnson et al. 1983;Shoberg et al. 1991;Davis & Villinger 1992). There have been several notable episodes of large volcanic earthquake swarms (Johnson et al. 2000;Bohnenstiehl et al. 2004;Hooft et al. 2010;Weekly et al. 2013), and microearthquakes from the central axial region are interpreted to represent ongoing magma inflation (Wilcock et al. 2009). There are some seamount chains, and a large elevated plateau that may be the site of enhanced crustal production related to a hot-spot anomaly that is associated with the Heckle Seamount chain (Karsten & Delaney 1989;Carbotte et al. 2008). The geological target for the FWI is the along-and acrossaxis variation in ridge structure and hydrothermal processes to the base of the crust ∼6 km below the seafloor.

Experiment
An OBS data set was acquired in 2009 across the Endeavour segment of the Juan de Fuca ridge to investigate crustal accretion and hydrothermal processes at an intermediate spreading segment (Weekly et al. 2014). The data were acquired with the R/V Marcus G. Langseth; a 36 element, 6600 in. 3 airgun array was used to fire ∼5500 shots which were recorded on 68 four-component OBS ( Fig. 1) from the Scripps Institution of Oceanography (SIO) and Woods Hole Oceanographic Institution. For this study, we have used only data acquired within the red box (Fig. 1). The OBS were all SIO instruments, and the average OBS spacing was ∼5 km. The shot spacing along a shot line was ∼450 m; the shot-line spacing in the centre of the study area was ∼450 m, and increased to ∼1 km near the edge of the red box. Either the vertical geophone or hydrophone data could have been used as input to FWI, but here we used the hydrophone data only. The amplitude response of these hydrophones is relatively flat between 6 and 80 Hz, and the response rolls off slowly below 6 Hz such that the amplitude at 3 Hz is about one third of that at 6 Hz. The shots used in this study were all fired at a 9-m tow depth. The seafloor bathymetry was mapped using a multibeam system, and water depth varies between 1750 and 3250 m across the survey area (see Fig. 1). Water-column velocities were determined using expendable bathythermograph profiles collected throughout the experiment (Weekly et al. 2014).

P-wave velocity and anisotropy
Crustal arrivals (P g ) are clear for the majority of stations and could be picked with uncertainties of between 10-15 ms; the RMS mean uncertainty for all the picks was 11.3 ms (Weekly et al. 2014). A tomographic regularized traveltime inversion (Toomey et al. 1994;Dunn et al. 2005) was used to invert first-arrival data and obtain models of isotropic slowness and anisotropy (Figs 2 and 3; Weekly et al. 2013). The chi-squared value was reduced from 13.19 to 1.17 after traveltime tomography (Weekly et al. 2014), which corresponded to an RMS misfit reduction from 39.3 to 11.7 ms. This meant that the majority of the predicted traveltimes matched the observed well and were within, or close to, the pick uncertainty. In the final velocity model, the uppermost oceanic crust (layer 2A) has a Pwave velocity of ∼2200-2400 ms −1 , and jumps to over 4000 ms −1 below layer 2A, and then gradually increases to ∼6900 ms −1 near the base of the crust. The recovered isotropic P-wave velocity model shows substantial lateral and vertical heterogeneity, with lower velocities along the ridge axis, surrounded by strips of high-and low-velocity zones at shallow depths . At depths greater than 2 km beneath the seafloor (Figs 2f-g), there is a relatively large zone of elevated velocities (blue) in the central zone away from the ridge, and reduced velocities (red) at both ends of the Endeavour segment. The nature of the P-wave anisotropy (Fig. 3) is as expected for oceanic crust at a spreading centre (Dunn & Toomey 2001), with the fast direction subparallel to the ridge and anisotropy decreasing both with increasing depth and distance from the ridge (Weekly et al. 2014). These models of P-wave velocity and anisotropy were used to construct starting models for FWI.

F U L L -WAV E F O R M I N V E R S I O N
FWI seeks to find a quantitative model of the subsurface that is capable of predicting the waveforms of seismic field data in detail (Tarantola 1984). This involves iteratively updating an initial starting model using a linearized least-squares local inversion (Pratt 1999;Plessix 2008). The principal benefit of FWI is that it has the potential to resolve subsurface properties to about half the seismic wavelength (Pratt et al. 1996;Virieux & Operto 2009), and this is a significant improvement on conventional traveltime tomography. Early FWI codes were 2-D (e.g. Pratt et al. 1996;Pratt & Shipp 1999;Shipp & Singh 2002;Brenders & Pratt 2007;Delescluse et al. 2011;Morgan et al. 2011), but FWI only became widely adopted by the petroleum industry when 3-D applications became practicable (e.g. Sirgue et al. 2010). The spectacular recent success of 3-D FWI has been achieved by the multi-azimuthal coverage, with many crossing wavefields providing multiple independent observations across each cell in a velocity model (Sirgue et al. 2007;Plessix & Perkins  2010; Li et al. 2011). FWI, as it has come to be widely applied across the petroleum industry, uses wide-angle, long-offset, lowfrequency data that are dominated by forward-scattered, refracted, transmitted arrivals (Sirgue 2006;Prieux et al. 2011;Vigh et al. 2011;Mothi et al. 2013;Vigh et al. 2013a;Yoon et al. 2014). With few exceptions Lu et al. 2013;Vigh et al. 2013b), commercial FWI uses an acoustic approximation to the wave equation (Virieux & Operto 2009;Kapoor et al. 2013) and commonly includes the kinematic (traveltime) effects of P-wave anisotropy (Bansal et al. 2013;Jones et al. 2013;Selwood et al. 2013;Warner et al. 2013). The clear verifiable success of acoustic inversions is somewhat surprising given that they cannot properly account for viscoelastic and other effects in field data.
For the inversions shown here, we use a time-domain, acoustic, anisotropic 3-D FWI code developed at Imperial-see Warner et al. (2013) for a detailed description of the code. To mitigate against issues with elasticity, the field data can be processed to remove S-wave arrivals and surface waves. These data are then bandpass filtered, and trace amplitudes normalized so that the rms value of each time-windowed trace is identical . The RMS amplitude of the predicted data is then scaled so that it matches that of the field data, trace-by-trace. In the time-domain FWI, used here, the misfit is defined as the sum of the squares of the difference between the normalized predicted and observed data, for every time sample. The misfit is then minimized, and this leads to an improvement in the match of both the amplitude and phase spectra of every trace over the bandwidth of the pre-processed data, and it matches the waveform of every trace over that bandwidth. Tests using both synthetic models and field data (for which inverted velocity models can be directly compared to sonic logs) indicate that this approach is robust and appears to avoid the introduction of artefacts into the recovered P-wave velocity caused by viscoelasticity and other hard-to-determine parameters that affect the dynamics (amplitudes) of the propagating wavefield (Morgan et al. 2013). It is rarely the case that good independent anisotropic models for attenuation, density and S-wave velocities exist, hence the use of viscoelastic codes is a challenging prospect as each of these parameters affects amplitudes. Although, ultimately it would be ideal to be able to propagate a viscoelastic wavefield in 3-D, and properly model absolute amplitudes, the extra computation time as well as crosstalk between parameters (Virieux & Operto 2009;Prieux et al. 2013) means that such inversions are currently not normally practicable in 3-D for field data and may give misleading results .
Note that, in the time-domain FWI used here, we do not invert at single frequencies and always invert finite-bandwidth data. Where particular frequencies are mentioned they refer to the cut-off frequency of a low-pass filter . As noted earlier, smoothing is not normally required when inverting densely sampled industry data, but is more likely to be necessary for sparser data. In the FWI algorithms used here, smoothing is applied as part of the pre-conditioning of the gradient rather than as a penalty within the objective function. Specifically, the raw gradient of the objective function with respect to the model parameters is computed together with a simple diagonal approximation to the Hessian matrix; the latter predominantly takes account of differences in illumination energy within the subsurface. This approximate Hessian is then used to pre-condition the gradient. Smoothing is then applied to this pre-conditioned gradient prior to the calculation of the step-length.

M E T H O D : W O R K F L O W
For the inversions shown here, we invert for P-wave velocity while keeping anisotropy fixed. The adopted workflow is the same as that described in Warner et al. (2013), which is summarized below: (1) build the source wavelet; (2) choose the starting frequency; (3) check adequacy of starting model; (4) pre-process the data; (5) devise a modelling and inversion strategy; (6) invert the data with continued quality assurance.

Build the source wavelet
For FWI, we require a source wavelet that is accurate at the frequency range we invert for. Fig. 4(a) shows a vertical-incidence source wavelet that has been modelled using the Nucleus program (H. Carton, private communication, 2015), which includes a source ghost (reflection from the sea-surface) for an airgun tow depth of 9 m that was used for acquisition. This wavelet is then deghosted as a free surface is used in the modelling which serves to reapply the ghost and thus properly accounts for directional changes in the propagating wavefield. A low-pass filter is applied to the source ( Fig. 4b) to pass only frequencies than can be accurately modelled given the chosen grid spacing-see modelling strategy below. In this case the filter was a minimum-phase Ormsby bandpass filter that rolled off from 6 to 9 Hz. The same filter is applied to the field data prior to input to the inversion. The adequacy of the source wavelet is verified later when we generate predicted data for the starting model and compare it with the field data-see below.
Note that, this method for deriving a source wavelet, automatically incorporates the response of the receiver. Although the amplitude response of hydrophones dips at low frequencies (<6 Hz), the adopted procedure of inverting narrow finite-bandwidth data in several steps with different frequency ranges, means that relative amplitude differences at different frequencies do not affect performance.

Choose the starting frequency
FWI typically starts by inverting the lowest possible frequencies in the data, to recover the long-wavelength velocity structure, and then gradually includes increasingly higher frequencies to recover the finer-scale structure (Bunks et al. 1995;Sirgue & Pratt 2004). Fig. 5(a) shows a phase plot at 3 Hz for the observed hydrophone data from a single OBS; the red rectangle in Fig. 1 shows the location of this plot. Dots are plotted at shot locations and their colour is the phase of a single trace which is obtained from data windowed in time using a Gaussian window centred upon the early refracted arrivals ). If source-generated signal exists in data at the plotted frequency, we will see a concentric structure around the OBS with phase for the observed data gradually changing with offset. Fig. 5(a) shows that there is clearly a good ratio of sourcegenerated signal to ambient noise at 3 Hz out to ∼10 km, and thus that the inversions can start at 3 Hz if the shot-receiver offsets are restricted to <10 km.

Check adequacy of starting model
The isotropic P-wave velocity and anisotropy models shown in Figs 2 and 3 were used to generate starting models for the inversion. In these figures, velocity is the mean P-wave velocity (Vp mean ) and fractional anisotropy is defined as (Vp Fast − Vp slow )/Vp mean , where Vp Fast and Vp slow are the P-wave velocities in the fast and slow direc-tion, respectively. We have adapted these velocity models for input to FWI using horizontal transverse isotropy, which assumes a slow axis of symmetry in the horizontal direction, and have converted the fractional anisotropy to the Thomsen's parameters (Thomsen 1986), Epsilon and Delta, using the assumption that the anisotropy is elliptical. Although Epsilon and Delta are unlikely to be equal, we have no independent control on these two parameters and this is the simplest assumption. Note that this model of anisotropy leads to both azimuthal and vertical variation in wave speed. An example vertical slice through the velocity model is shown Fig. 6(a).
An additional reason as to why the inversion is started at the lowest possible inversion frequency is to avoid problems with cycle skipping, which occur if the starting model is unable to predict the majority of data to within half a cycle of the field Figure 5. Phase plots for a receiver gather, OBS 32, viewed in the frequency domain at a single frequency of 3 Hz (a-c) and 4.5 Hz (d-f). The red rectangles in Figs 1-3 show the location of this plot. Coloured dots are the phase of a particular trace and are plotted at the shot location. The OBS location is towards the bottom right and, for the observed (a, d) and predicted (b, e) data, is the centre of the concentric features; spatially coherent structure in these receiver gathers represents source-generated signal. (c, f) Plots of the phase residual (observed minus predicted phase) at 3 and 4.5 Hz, respectively. data (Sirgue et al. 2010). For a 3-Hz starting inversion frequency, this corresponds to ∼167 ms. Given that the difference between the observed and predicted first-arrival traveltimes was much smaller than this (mean RMS misfit of 11.7 ms, Weekly et al. 2014), we would not expect to have a problem with cycle skipping.
The source wavelet shown in Fig. 4(b) and starting models for P-wave velocity (Fig. 6a) and anisotropy were used to generate predicted data for the real experimental geometry. Phase plots are used to verify whether data predicted using the starting model are cycle skipped (Fig. 5). Figs 5(b) and (e) show phase plots for the predicted data at 3 and 4.5 Hz, respectively, and Figs 5(c) and (f) show the residual (observed-predicted phase). The dots in residual plots are green when the two phases match (phase difference = zero). The residual phase plot at 3 Hz (Fig. 5c) shows the data appear to be well matched-the colours are close to green. There is no evidence of a significant problem with cycle skipping for offsets out to 10 km, which would be indicated by a rapid change from blue to red (a phase switch from −180 • to +180 • ) or vice versa. At 3 Hz the longer offset data are particularly noisy, as shown by the rapid variation in colour between adjacent dots in the top left hand corner of the plot in Fig. 5(a). Plots for this and other OBS suggest the starting model is sufficiently accurate, with the majority of the data at offsets <10 km being predicted to within half a cycle of the observed at 3 Hz. Fig. 7 shows some representative examples of observed and predicted traces for the hydrophone channel of OBS 32 at offsets of around 8 km, for which the predicted data are obtained using the starting model. It is clear that individual peaks and troughs between about 3.3 and 4 s, which are all Pg in these plots, are closely matched for the majority of the data, giving us additional confidence that the starting models for velocity and anisotropy, as well as the source Downloaded from https://academic.oup.com/gji/article-abstract/204/2/1342/596388 by Imperial College London Library user on 08 January 2018 Figure 7. Examples of the match between field and modelled data for the hydrophone channel of OBS 32 (see Fig. 1 for location). Shot-receiver offsets are ∼8 km. Predicted data are generated using the starting model in Fig. 6(a) and the source wavelet in Fig. 4(b).
wavelet, are all adequate. The later arrivals, however, are not well matched. Given that the traces in Fig. 7 are for adjacent shots along a line, and that shots are only about 450 m apart, the secondary arrivals are quite variable in both the observed and predicted data, and the relative amplitudes are larger in the predicted data. Modelling suggests that these secondary arrivals are Pg waves that have been reflected back downwards from the sea-bottom interface before being refracted back towards the OBS. It is clear, by comparing the observed and predicted data at >4 s in Fig. 7, that changes in seafloor topography and/or impedance contrast have not been captured in the starting model. The amplitudes of the later arrivals in the predicted data may be too high because the seafloor density contrast is too large, or possibly because the acoustic code does not model leakage into S-wave energy.
Notwithstanding this mismatch for later arrivals, these tests demonstrate that the starting model and source wavelet do adequately predict the early arriving refracted arrivals.

Pre-process the data
In industry data sets, there is often some data redundancy and it is usual to decimate field data prior to input to the inversion. Here, every hydrophone trace is input to the inversion. The source ghosts are left in the data and a free surface is used in the modelling; this procedure was found to give the most reliable and stable results as it properly accounts for directional changes in the source wavefield . For input to the inversion data are bandpassed using the same Ormsby filter that was applied to the source wavelet.

Devise a modelling and inversion strategy
As data sparsity was expected to be problematic, an area was selected for our initial inversions where both the shot and OBS spacings were most dense (red rectangle in Fig. 1), with a shot and shot line spacing of ∼450 m and OBS spacing of ∼5 km. This region includes 22 OBS and 1673 shots, but station 41 failed to record data, and station 45 had timing issues -see section on quality assurance presented after the results. Hence, data from a total of 21 OBS were used in the initial tests, but only data from 20 OBS were input into the final inversions. The maximum frequency we can invert for is determined by the requirement to have ∼4 grid points per wavelength. The minimum velocity is slightly less than 1500 ms −1 in the water column, and a grid spacing of 50 m was chosen to allow inversions up to a maximum frequency of about 6.5 Hz. The velocity model is 20 × 27 km wide and 8.95 km deep, corresponding to 401 × 541 × 180 cells. Data reciprocity was applied so that the 21 receivers were treated as shots, and the 1673 shots were treated as receivers. This arrangement is more efficient since all computations for an individual source are performed on a single node. All shots and receivers were inverted at every iteration.
Downloaded from https://academic.oup.com/gji/article-abstract/204/2/1342/596388 by Imperial College London Library user on 08 January 2018 Figure 8. Example of depth slice through recovered checkerboard (a) without and (b) with horizontal smoothing of half a seismic wavelength for each iteration. A checkerboard of 1 × 1 × 1 km in size was added at a depth of between 3.5 and 4.5 km bsl to a velocity model recovered using FWI. The velocity perturbation is ±100 ms −1 and the depth slice is at 3.8 km bsl. For the inversion, 10 iterations have been run at each of 3.0, 3.4 and 3.9 Hz. The synthetic data have been windowed prior to inversion using a window length of 1200 ms in order to replicate the procedure adopted when inverting field data. Windowing the input data is important in checkerboard inversions since, if all the data are included, the recovered check is surprisingly good due to the inversion of scattered energy from the checkerboard boundaries.
During the inversion P-wave velocity is updated and anisotropy is kept fixed. A number of different inversion strategies were tested, with different frequency ranges and data window lengths, as well as a variety of different vertical and horizontal smoothing parameters. Checkerboard tests were then used to verify whether the particular inversion strategy could resolve structures of the size recovered in the velocity models, and investigate whether the strategy led to the generation of any velocity artefacts. It was apparent from these tests that it was essential to use extra smoothing during each iteration, as the recovered checkerboard was otherwise quite noisy (see Fig. 8). This is in contrast to inversions performed on industry data, for which only minimal smoothing of the gradient appears to be required (Plessix & Perkins 2010), presumably because of data redundancy.
A top mute (trace value set to zero) was applied ahead of the first arrivals and data were windowed to pass only the early arrivals. Window lengths of between 700 and 1800 ms were investigated when devising the final inversion strategy. From these tests it was clear that we were unable to reproduce the complexity of the secondary arrivals and that, when we used long window lengths, the inversion failed to improve the match to the secondary arrivals while, at the same time, the fit between the early arrivals often worsened (Fig. 9a). We considered that the problem may be related to the abrupt 2A boundary at about 400 m below the seafloor in the starting model (Fig. 6a), which had been included in the original traveltime tomographic model by Weekly et al. (2014) using estimates of its thickness from MCS data (Van Ark et al. 2007). It is generally preferable to start FWI with a smooth model, unless intermediate velocities and depths to any sharp boundaries are well constrained, such as is often the case for the seafloor. The original starting P-wave velocity model (Fig. 6a) was smoothed vertically and horizontally below the seafloor, over a distance of about half a seismic wavelength at the lowest inversion frequency, 3 Hz, to remove structure below the theoretical resolution of FWI (Fig. 6b). Apart from smoothing layer 2A, this also led to the removal of the vertical striping that was present in the original velocity model, and caused by discretizing the model and bathymetry on a 50-m grid. Predicted data that are generated using the smoothed starting velocity model better match the early arrivals in the observed data (Fig. 9b), but not the later arrivals, and, inversions using long window lengths still led to a visible deterioration in overall fit. Only short window lengths (<800 ms) led to a consistent improvement in misfit during the inversion.
In the results shown in the next section, the input data are windowed using a window length of 750 ms, which was chosen to include only the early arrivals that were matched well by the Figure 9. (a) Example of field data recorded on the hydrophone channel of OBS 32 (see Fig. 1 for location) and modelled data generated using the starting model in Fig. 6(a), and for the inverted FWI velocity model after 6 and 22 iterations, with 10 iterations at 3.0 Hz, 10 iterations at 3.4 Hz and 2 iterations at 3.9 Hz. Shot-receiver offset is ∼6.5 km. (b) Observed trace is as in (a), and modelled data are generated using the smooth starting model in Fig. 6(b).
Downloaded from https://academic.oup.com/gji/article-abstract/204/2/1342/596388 by Imperial College London Library user on 08 January 2018 starting model (see Figs 7 and 9b). Smoothing of 1.2 and 0.8 seismic wavelengths is applied in the horizontal and vertical direction, respectively. These values represent the minimum smoothing that could be used and still recover a checkerboard reasonably well when using the relatively short window length-see section on quality assurance after the results. Data between offsets of 2.9-10 km were inverted with 10 iterations at 3.0, 3.4, and 3.9 Hz, and then offsets between 2.9 and 15 km were input into the inversion and 5 iterations performed at 3.3 and 3.8 Hz, and then 10 iterations at 4.4 and 5.1 Hz. Longer offsets were excluded from the inversion to avoid data that were adversely affected by noise (Figs 5a and d). At offsets closer than 2.9 km, the first arrival is a complex interference pattern between large amplitude direct arrivals through the water layer and turning rays in the uppermost crust. As there was independent data on water velocity and seafloor topography, these were both kept fixed and the arrivals at offsets <2.9 km were excluded in the inversion. Finally, in applications of the Imperial 3-D FWI codes to industry data, spatial preconditioning is used to boost the gradient (gradient of the functional with respect to model parameters) where the incident wavefield is small, and reduce it where the wavefield is large; a strategy that works well for highly redundant over-sampled data. However, when checkerboard tests were performed for the Endeavour experiment, spatial preconditioning was found to produce artefacts in deeper parts of the velocity model, presumably because the deeper model regions are under-sampled. Hence, spatial preconditioning was switched off for the final inversions of the field data.

Invert the data with continued quality assurance
The phase plot in Fig. 5 shows that data were not cycle-skipped at offsets less than 10 km at the starting frequency for the inversion, 3 Hz. Phase plots are also used to ascertain whether the inversion is iterating towards a global minimum. In addition, all investigated inversion strategies were validated with checkerboard tests to confirm that the chosen smoothing, number of iterations at each frequency, and window length could recover the anomalies seen in the inverted velocity models. These and some additional tests are presented after the results section.

R E S U LT S
Figs 10(a)-(f) show three vertical slices though the recovered velocity model using the inversion strategy outlined above, plotted next to the starting model. Figs 10(g-i) are three depth slices which, for display purposes, show the velocity perturbation-the difference between the inverted and starting velocity model. The velocity models are plotted down to 6 km below sea level (bsl) as this is considered to be the deepest depth where recovered velocity anomalies are likely to be reliable-see section on quality assurance. The vertical slices show sections perpendicular to the ridge (see Fig. 2 for location), and have been chosen as they illustrate that fine-scale velocity structure can be recovered using FWI, and that velocity anomalies beneath the ridge evolve along axis. Low-velocity zones beneath the ridge are likely to be related to regions with a higher temperature and melt content (Arnoux et al. 2014); the slice in Fig. 10(f) appears to show some off axis magma storage at about 4.5 km bsl, possibly fed by a conduit that slopes back towards the ridge axis. This low-velocity zone lies just below a topographic high formed by extrusive volcanism (see Fig. 2 for location). Evidence for off-axis magmatism has been reported elsewhere on the East Pacific Rise (Durant & Toomey 2009;Canales et al. 2012). Also apparent in the vertical slices are fine-scale changes in layer 2A, which are consistent with velocity variations in the uppermost oceanic crust in the original traveltime inversions (Weekly et al. 2014), as well as with MCS data: Van Ark et al. (2007) noted a variable thickness in layer 2A. The depth slices (Figs 10g-i) also show some interesting features. The slice at 2700 m bsl shows bands of reduced and elevated velocities beneath the ridge, whereas at 3800 m bsl (Fig. 10h) there is a zone of reduced velocity directly beneath the ridge that gradually becomes wider and more centralized by 5500 m bsl (Fig. 10i). A detailed analysis of the velocity model recovered with FWI is presented in Arnoux et al. (2016), together with its interpretation with respect to hydrothermal processes and the accretion of oceanic crust.

Q UA L I T Y A S S U R A N C E
The progress of the inversion was monitored using phase plots. Fig. 11 shows the residual phase (observed minus predicted phase) for the hydrophone channel of OBS 32 at 3 Hz for the: (a) starting model, (b) model after one iteration and (c) final model. The residual phase plots for the starting and first iteration were plotted for all 21 OBS to test for cycle skipping and check whether the inversion stepped in the correct direction. Synthetic tests show that, if cycle-skipping is an issue, the region of cycle-skipped data expands and moves inwards to shorter offsets ). There are no clear signs of cycle skipping in Fig. 11. Phase plots for OBS 45 revealed a higher noise content than other OBS, and the residual phase plot indicates that the observed and predicted phases were nearly 180 degrees apart (Fig. 12), hence data from this OBS were not included in the final inversions. All the other OBS had a good signal content for offsets out to ∼10 km, and the phase plots showed only minor changes after one iteration, as seen in Fig. 11(b). Subsequently, phase plots for a small selection of OBS across the model area were checked every 5 and then 10 iterations. In each case, the phase plots indicated that the inversions produced small or negligible improvements in the velocity model, and the phase residuals gradually moved towards 0 • (green). Fig. 11(c) is a plot of the phase residual for the final model at 3 Hz, which shows a clear improvement over the starting model. Figs 11(d and e) show the residual phase for OBS 32 at 4.5 Hz for the starting and final   Fig. 7 and (b) the same trace as shown in Fig. 9. Predicted data are generated using the smooth starting model in Fig. 6(b) and the final FWI model shown in Fig. 10. The location of the peaks and troughs for the data predicted by the FWI model are closer to the observed data overall, and the relative amplitudes of the peaks and troughs are much better matched for traces 1-3 in panel (a).
Downloaded from https://academic.oup.com/gji/article-abstract/204/2/1342/596388 by Imperial College London Library user on 08 January 2018 velocity model, respectively. Data from longer offsets are included in this plot as data out to 15 km offset were input to inversions at this frequency. The match between predicted and observed phase is clearly significantly better for the final velocity model. In summary, the phase plots show no evidence of cycle-skipping and suggest that FWI is moving the velocity model towards a global not a local minimum. Fig. 13(a) shows the same five traces as in Fig. 7. The observed traces are examples of data input to the inversion with the 750 ms window applied. The predicted data are generated from the starting and final FWI velocity models. The FWI model produces data that better match the observed: the travel times to individual peaks and troughs are closer overall, and there is significant improvement in the waveform match. For example, the amplitudes of the peaks and troughs in the later arrivals in traces 1-3 are all closer to the observed. Fig. 13(b) shows the same trace as in Fig. 9; here the traveltime fit is much better for the FWI model and the amplitude match is about the same. Note that the starting model is the smoothed velocity model shown in Fig. 6(b), for which the predicted data is shown in Fig. 9(b). Fig. 14 shows the change in misfit for each iteration for the final inversion strategy. The misfit is reduced by ∼14 per cent for the first 10 iterations at 3 Hz and then, thereafter, reduces by ∼5 per cent every 10 iterations. This relatively small reduction in misfit is in accordance with the observed improvement in match for individual traces (Fig. 13): the initial fit is quite good and the improvement is clear but often quite small.
While testing different inversion strategies for the field data, the same strategies (smoothing, window lengths, iteration sequence) were also applied to recover checkerboards of various sizes and depths within the velocity model. These checkerboard tests were run to verify that the adopted inversion strategy could recover the anomalies seen in the inverted velocity models. Fig. 15 shows depth slices (where depth is bsl) through checkerboards that were all recovered using the inversion strategy adopted in the final FWI models shown in Fig. 10. Fig. 15(a) shows the successful recovery of a 1.5 × 1.5 × 1.0-km sized checkerboard that was placed between 2.8 and 3.8 km bsl. Below this depth only larger checkerboards could be recovered. Fig. 15(b) shows the recovery of a checkerboard that was 2.5-km wide and 2-km deep and placed between 3 and 5 km bsl; Figs 15(c)-(d) show the recovery of the same checkerboard located between 4.5 and 6.5 km bsl. The smaller checks are well recovered down to about 3.5 km bsl, and the larger checks are well recovered down to about 4.5 km bsl. Below this, the checks are partially recovered down to about 6 km bsl, by which we mean that the velocity perturbations are at the correct location and have the correct sign, but the absolute size of the perturbation is not recovered. Hence, it can be concluded that, away from the model edges, the velocity structure in the final inversions is likely to be real, but the magnitude of the recovered anomalies is probably too small in deeper parts of the model.
To explore further whether the final inversions are robust and produce consistent results, two additional tests were performed. In the first, the final velocity model shown in Fig. 10 was used to generate synthetic data, which was then treated as the observed data and inverted using the original inversion strategy and starting model. Fig. 16 shows the results of this test above the original FWI model. Each of the depth slices have the same structural features and all the areas of high and low velocity are matched. The re-recovered velocity model is slightly smoother, and only about 80 per cent of the original velocity perturbation has been recovered. In this test, unlike for the field data, the observed data are noise free. Hence, the results shown in Fig. 16 suggest that the velocity anomalies obtained from inverting the field data are recoverable using the selected inversion strategy, and are not artefacts produced by noise within the field data. This test also suggests that our final inversion strategy will not fully recover the magnitude of the real velocity anomalies.
In the second test, a new starting model was generated by taking the velocity perturbations recovered by FWI, multiplying them by 1.5 and adding them to the original starting velocity model (Figs 17a-c). The inversion was performed again using the same input field data and with the original inversion strategy, but substituting this model as the starting model. Such an approach is useful in FWI since the inversion tends to head in the right direction but may take many iterations to get there. The process of multiplying the recovered perturbations by 1.5 can move the inversion more rapidly towards the global minimum. Figs 17(d)-(f) show that the inversion moves the recovered anomalies back towards the original FWI model (Figs 17g-i). The results of this test, together with the other quality assurance procedures that have been applied, lead us to the conclude that the velocity structure in our final velocity model is robust. The recovered features in the FWI model shown in Fig. 10 consistently appear in the inverted velocity models.

D I S C U S S I O N
The checkerboard tests shown in Fig. 15 indicate that FWI can improve the model resolution by 2-4 times over traveltime tomography. A suite of checkerboard tests were performed by Weekly et al. (2014), and these showed that traveltime tomography could recover 2.5 × 2.5 × 1-km sized checks in the uppermost 2 km of the crust, and 5 × 5 × 2-km checkerboards down to 3 km below the seafloor. The recovery of finer-scale structure in the FWI model shown in Figs 10(b), (d) and (f) is encouraging and suggests that the approach adopted here could be successfully applied to other relatively sparse data sets. The theoretical resolution for FWI is half the seismic wavelength (Pratt et al. 1996) which, for the final inversion frequency of 5.1 Hz, corresponds to ∼300 m in the uppermost crust and ∼750 m in the lower crust at the Endeavour ridge. The resolution of the FWI velocity models presented here, which correspond to about 1.5 km in the upper crust and 2.5 km in the lower crust (Fig. 15), indicate that we have only obtained about 1/3-1/5 of the potential resolution of this technology. A suite of additional inversions were attempted to improve the resolution, including many Figure 15. Checkerboards at (a) 3200, (b) 3800, (c) 4800 and (d) 5400 m bsl, recovered using the same inversion strategy as in the final FWI model shown in Fig. 10. The checkerboard in panel (a) was 1.5 wide and 1.0 km deep and placed between 2.8 and 3.8 km bsl. In the other three plots, a checkerboard 2.5-km wide and 2-km deep was placed between (b) 3 and 5 km bsl and (c,d) 4.5 and 6.5 km bsl. The velocity perturbation was ±100 ms −1 .
Downloaded from https://academic.oup.com/gji/article-abstract/204/2/1342/596388 by Imperial College London Library user on 08 January 2018 In the re-recovered velocity model, the final velocity model shown in Fig. 10 was used to generate synthetic data, which were then treated as the observed data and inverted using the original inversion strategy and starting model. Each of the depth slices has the same structural features in the re-recovered and original FWI models, and all the areas of high and low velocity are matched. more iterations and the inversion of higher frequencies, but the recovered finer-scale features were more variable and the available quality assurance procedures were unable to distinguish whether one of these models was definitively better than another. Resolution for the Endeavour data set is almost certainly affected by the data sparsity, and may also be limited by noise and the fact that the original starting model was unable to accurately predict the secondary arrivals.
In order to investigate how the full potential resolution of FWI could be obtained, a suite of checkerboard tests were performed to explore the effect of denser experimental geometries. Since the expected resolution at 5 Hz is ∼300 m in the upper crust and ∼750 m in the lower crust, three checkerboards were placed in the model: a 300-m check was placed between 3000 and 3300 m, a 500-m check was placed between 4000 and 4500 m, and a 750-m check was placed between 5000 and 5750 m. Fig. 18 shows the recovered checkerboards for three experimental geometries: (ac) OBS spacing of 1 km and shot spacing of 250 m, (d-f) OBS spacing of 2 km and shot spacing of 250 m, and (g-i) OBS spacing of 2.5 km and shot spacing of 400 m. Not surprisingly, the checkerboards are best recovered using the closest OBS and shot spacing, but the checkerboard is still partially recovered using the  larger OBS and shot spacing. These tests suggest that we could obtain the full potential resolution of this technique using OBS spacings of between 2-2.5 km while inverting for frequencies up to 5 Hz.

C O N C L U S I O N S
The inverted velocity anomalies in the final FWI model are roughly 2-4 times finer than can be recovered using traveltime tomography. This improvement in resolution demonstrates that 3-D FWI can be used to recover fine-scale structure within the crust even when data are acquired using a relatively sparse shot and receiver spacing. In contrast to inversions performed on industrial data, additional regularization was required, with smoothing in both the horizontal and vertical directions. In addition, a laborious quality assurance procedure is needed to verify the validity of the adopted workflow, including checking the adequacy of the source wavelet and starting model, checking for problems with cycle skipping, investigating the effect of inputting different data bandwidths, window lengths and regularization, monitoring the progress of the inversion and performing additional model assessment. The full potential resolution of FWI could be obtained in future marine seismic surveys if data are collecting using a 3-D acquisition geometry with a denser OBS spacing than is typical for academic surveys.

A C K N O W L E D G E M E N T S
Two reviewers made constructive comments that helped improve the manuscript. Imperial College London wishes to thank the sponsors of the FULLWAVE consortia: BG, BP, Chevron, CGGVeritas, ConocoPhillips, DONG, DownUnder GeoSolutions, ENI, HESS, Maersk, Nexen, PGS, Rio Tinto, Statoil, Sub Salt Solutions, TGS, Tullow Oil and Woodside for support in developing the 3-D FWI software. The 3-D ProMAX/SeisSpace package, supplied by Halliburton Software and Services, a Halliburton Company, under a university software grant, was used to pre-process and analyse the field and synthetic data. The inversion software used in this study is available for application to academic problems through collaboration with Imperial College London, and to commercial partners through membership of the FULLWAVE consortium. The experiment and analysis were supported by the National Science Foundation (NSF) under grants numbered OCE-0454700 to the University of Washington and OCE-0454747 and OCE-0651123 to the University of Oregon. The Department of Geological Sciences, University of Oregon, supported the academic visits by Arnoux and Vander-Beek to Imperial College London. The data used in this research were provided by instruments from the Ocean Bottom Seismograph Instrument Pool (http://www.obsip.org) which is funded by the NSF. OBSIP data are archived at the IRIS Data Management Center (http://www.iris.edu).