# Lower oceanic crust formed by in situ melt crystallization revealed by seismic layering

Peng Guo, Satish C. Singh, Venkata A. Vaddineni, Ingo Grevemeyer, Erdinc Saygin

► To cite this version:

Peng Guo, Satish C. Singh, Venkata A. Vaddineni, Ingo Grevemeyer, Erdinc Saygin. Lower oceanic crust formed by in situ melt crystallization revealed by seismic layering. Nature Geoscience, 2022, 15, pp.591-596. 10.1038/s41561-022-00963-w. insu-03748520

HAL Id: insu-03748520

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

Submitted on 9 Aug 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>

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

OPEN

# Lower oceanic crust formed by in situ melt crystallization revealed by seismic layering

Peng Guo<sup>1, ∗</sup>, Satish C. Singh<sup>1, †</sup>, Venkata A. Vaddineni<sup>2</sup>, Ingo Grevemeyer<sup>1, ‡</sup> and Erdinc Saygin<sup>1, §</sup>

Oceanic crust forms at mid-ocean spreading centres through a combination of magmatic and tectonic processes, with the magmatic processes creating two distinct layers: the upper and the lower crust. While the upper crust is known to form from lava flows and basaltic dykes based on geophysical and drilling results, the formation of the gabbroic lower crust is still debated. Here we perform a full waveform inversion of wide-angle seismic data from relatively young (7-12-Myr-old) crust formed at the slow-spreading Mid-Atlantic Ridge. The seismic velocity model reveals alternating, 400-500 m thick, high- and low-velocity layers with ±200 m s⁻<fnref n="1" /> velocity variations, below ~2 km from the oceanic basement. The uppermost low-velocity layer is consistent with hydrothermal alteration, defining the base of extensive hydrothermal circulation near the ridge axis. The underlying layering supports that the lower crust is formed through the intrusion of melt as sills at different depths, which cool and crystallize in situ. The layering extends up to 5-15 km distance along the seismic profile, covering 300,000-800,000 years, suggesting that this form of lower crustal accretion is a stable process.

Oceanic crust is formed at the ocean spreading centres with different accretionary mechanisms<fnref n="1" />⁻³, depending on the spreading rates. At fast- and intermediate-spreading centres, the crustal accretion is mainly dominated by magmatic processes, whereas at slow- and ultra-slow-spreading ridges, the crustal formation is influenced by active tectonic processes such as faulting, leading to heterogeneous crustal structures⁴. Additionally, hydrothermal circulations play important roles by releasing heat and altering rock properties⁵. Based on insights from ophiolite studies⁶,⁷, ocean drilling results⁸ and geophysical studies⁹,¹⁰, the crust formed by magmatic processes can be divided into two layers, an upper and a lower crust. The upper crust is primarily composed of pillow lavas and the underlying sheeted dykes with a high velocity gradient (1-2 s⁻<fnref n="1" />)⁹,¹¹, whereas the lower crust contains mainly gabbro where the velocity increases very slowly (~0.1 s⁻<fnref n="1" />)³,⁹,¹⁰. When the magma supply is low, the tectonic process dominates, and the exhumation of mantle rocks in the crust and on the seafloor can occur along detachment faults⁴. Although the upper crustal structure and accretion process have been well studied⁵,⁸⁻¹³, our knowledge about the lower crust remains poorly constrained¹⁴. The lower crust could represent two-thirds of the oceanic crust³, therefore it is important to understand how this part of the crust is formed and evolves.

The current debate on the accretion process of the lower magmatic oceanic crust centres on two end-member models. The 'gabbro glacier' model<fnref n="1" />,¹⁵ suggests that the lower crust is formed by cooling, crystallization and subsequent subsidence of the gabbro from the axial melt lens (AML) at the upper-lower crust boundary, while the 'sheeted sill' model⁶,⁷ suggests that the lower crustal gabbro formation could occur through in situ cooling and crystallization of melt sills. Recent seismic discoveries of melt lenses in the lower crust beneath the fast-¹⁶⁻¹⁸ and intermediate-spreading mid-ocean ridges¹⁹⁻²¹ provide evidence for the sheeted sill model for the lower crustal accretion. However, no such melt lens has been observed in the lower crust beneath slow-spreading ridges,

although the presence of partial melt has been suggested based on low velocity anomalies underneath the ridge axis²²⁻²⁴. Furthermore, how these structures beneath the ridge manifest themselves as they move away from the ridge axis remains unknown.

## Seismic velocity model of the oceanic crust

In 2017, wide-angle seismic data²⁵ were acquired aboard the German RV *Maria S. Marian* in the equatorial Atlantic Ocean (Fig. 1) over crust formed at the slow-spreading Mid-Atlantic Ridge (MAR) with a half spreading rate of ~16 mm yr⁻<fnref n="1" /> (ref. ²⁶). An air gun array with a total volume of 5,440 cubic inches was fired at a ~400 m interval, which was recorded on ocean bottom seismometers (OBSs; Supplementary Fig. 1) deployed at a spacing of 10–20 km on the seafloor. The crustal (Pg), mantle arrivals (Pn) and wide-angle reflections from the Moho (PmP) were used to obtain the P-wave velocity structures. Here, we use a part of these data (8 OBSs) where the OBS spacing was dominantly 10 km (Fig. 1), covering 7–12-Myr-old seafloor on the African plate. The velocity model obtained using the travel-time tomography¹⁴ (Methods) is shown in Extended Data Fig. 1. The crustal thickness in the study area is 5.6 ± 0.2 km (ref. ¹⁴). As expected, the tomographic velocity model contains only the large-scale velocity structures but clearly shows upper crust with high velocity gradient down to ~2.4 km below the basement underlain by a low-velocity-gradient lower crust¹⁴.

Starting from the travel-time tomographic model¹⁴ (Extended Data Fig. 1), we performed a full waveform inversion (FWI; Methods)²⁷⁻²⁹. The FWI is based on matching full wavefield of the real data and synthetically calculated data, and provides a detailed velocity information of the sub-surface²⁹. The FWI was performed in two steps: (1) a trace-normalized FWI inversion²⁷ was performed to obtain an intermediate-scale velocity model; and (2) then a true-amplitude FWI²⁸,²⁹ was performed to obtain fine-scale velocity information (Methods). The final P-wave velocity model from the true-amplitude FWI is shown in Fig. 2. The synthetically

<footnote n="1">¹Deep Earth Imaging Future Science Platform, The Commonwealth Scientific and Industrial Research Organisation (CSIRO), Kensington, Western Australia, Australia. ²Laboratoire de Géosciences Marines, Institut de Physique du Globe de Paris, Université de Paris Cité, Paris, France. ³GEOMAR Helmholtz Centre for Ocean Research, Kiel, Germany. ⁴Department of Physics, School of Physics, Mathematics and Computing, Faculty of Engineering and Mathematical Sciences, University of Western Australia, Crawley, Western Australia, Australia. ⁵e-mail: peng.guo@csiro.au</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1 | Map of the seismic survey area near the MAR in the equatorial Atlantic Ocean.** The dashed black line refers to the MAR and the triangles show the locations of the OBSs used for this study. The bold black line passing through the OBSs shows the part of shot profile used here, and the two pink bars indicate the boundaries of the seismic images shown in Fig. 2. The thin black line and numbers indicate the distance from the MAR. The contour labelling 10 indicates the 10 Myr lithospheric age²⁶. The bathymetry was taken from ref. ¹⁴. See the inset map for the location of the study area.</figcaption></figure>

calculated data after the FWI show a much better match to the crustal Pg arrivals of the field data than those from the tomographic model (Fig. 3 and Extended Data Figs. 2–5). To make sure the inverted results are real, not artefacts of inversion, extensive modelling and inversion tests were performed (Methods and Supplementary Figs. 2–23), which shows that the inverted velocity models are required by the data, and hence represent real structures of the sub-surface. Resolution analysis suggests that structures on the scale of 400 m vertically and 5 km horizontally can be resolved by the data (Methods and Supplementary Figs. 14–20).

When compared with the tomographic model, the velocity anomaly obtained using trace-normalized FWI²⁷ (Methods and Extended Data Figs. 6 and 7) indicates an overall increase in the velocity (200–400 m s⁻¹ positive anomaly) above 6–7 km depth and a decrease in velocity (100–300 m s⁻¹ negative anomaly) below this depth. Interestingly, there are alternate high- and low-velocity layers below ~2 km depth from the basement in the vertical velocity gradient (Fig. 2c and Extended Data Fig. 8c) from the true-amplitude FWI²⁸,²⁹, extending up to 5–15 km distance along the profile.

One-dimensional (1D) velocity profiles (Fig. 4a and Extended Data Fig. 9) show three distinct features: (1) a layer with a much higher vertical velocity gradient than that of the tomographic model from the basement down to ~1.2 km depth, where the velocity linearly increases from an average of ~4.8 km s⁻¹ at the top to 6.2 km s⁻¹ at its base; (2) an underlying medium velocity gradient layer of 600–800 m thickness, with the velocity increasing to ~6.6 km s⁻¹; and (3) the alternate high- and low-velocity layers (400–500 m thick) with a velocity variation of ±200 m s⁻¹ underlying.

## The geological structure of oceanic crust

The upper layer with high velocity gradient is probably associated with the lava flow-dyke (layer 2A/2B, Fig. 4a) boundary¹,¹²,¹³,³⁰, where drilling results³⁰,³¹ suggest that the rapid reduction of porosity is the main reason for this fast velocity change. A porosity decrease will cause an increase of velocity (Supplementary Fig. 24). An abrupt transition from low-temperature hydrous alteration above to high-temperature hydrothermal alteration below has been observed

at the lava-dyke transition zone⁸,³⁰, which corresponds to lithology changes including a sudden appearance of hydrothermal minerals, and an abrupt reduction of permeability and porosity³⁰.

A recent tomographic study³² using ultra-long streamer data coincident with the OBS data indicates this boundary at 600–800 m below basement for ~7–12 Ma, suggesting a shallower 2A/2B transition than our results. However, these authors found the lava flow-dyke boundary at ~900 m below the seafloor underneath the ridge axis³². While streamer data may have better sampled the uppermost crust, the relatively low resolution of the tomographic inversion might smear the 2A/2B transition with gentler velocity gradient and cause this discrepancy.

The 600–800 m thick, medium velocity gradient layer may represent the dyke sequence. The drilling results from the Ocean Drilling Program (ODP) Hole 1256D indicate a thin (350 m) dyke sequence³³, whereas a thick (>1,000 m) dyke sequence is observed in Hole 504B³⁴; our results lie between these two extreme values. If the base of this sequence marks the base of the upper crust, the upper crust would be 1.8–2.0 km thick (Fig. 4a).

The alternating high- and low-velocity layers below ~2 km depth from the basement could be explained by the presence of: (1) hydrothermal alteration⁸,³⁵,³⁶; (2) alternate gabbro sills and serpentinized peridotite within the lower crust, possibly associated with oceanic core complexes (OCCs)⁴,³⁷; and/or (3) layering introduced either by the gabbro glacier model¹,¹⁵ (crystallization and subsidence of the gabbro from the shallow AML) or the sheeted sill model (in situ crystallization of melt sills in the lower crust beneath the ridge axis)⁶,⁷.

The hydrothermal circulation is expected to reach down to the top of the AML beneath the ridge axis¹⁰,³⁵. A seismic reflection study at the Lucky Strike segment of the MAR has shown the presence of an AML at 3 km depth from the seafloor³⁸, much deeper than the top low-velocity layer in the FWI model at 2 km depth, therefore the top layer is unlikely to be the fossil melt lens from the ridge axis. On the other hand, previous seismic investigation³⁵ has imaged a 150–200 m thick, low-velocity (0.3–0.5 km s⁻¹ decrease) zone above the roof of the AML, which was interpreted to be formed from enhanced hydrothermal alteration³⁵. Furthermore, a velocity decrease of 0.4–0.6 km s⁻¹ in a 100–200 m thick layer has been observed in Hole 1256D at the dyke-gabbro transition³⁶, with noticeable changes in mineralogy and decreasing porosity associated with hydrothermal alteration and contact metamorphism. A recent drilling result³⁹ from the Oman ophiolite of the dyke-gabbro transition reveals that this transition is characterized by low-Mg and Zr-enriched gabbro and diorite, possibly related to the interaction between the magmatic and hydrothermal systems at the AML roof³⁹,⁴⁰. Therefore, the topmost low-velocity layer at ~2 km depth could be associated with an enhanced hydrothermal circulation, and the underlying high-velocity layer at ~3 km depth could be the roof of the AML³⁵, consistent with the observation at the Lucky Strike³⁸. Below this depth, the role of the hydrothermal circulation would be much reduced⁵.

## Lower crustal accretion processes

Oceanic crust formed in a slow-spreading environment is suggested to be heterogeneous and may even contain mantle rocks²,⁴. Tectonic processes, such as long-lived detachment faults⁴, can produce OCCs, exhuming lower crustal and mantle rocks to the seafloor⁴,³⁷. For example, FWI of multi-channel streamer data from the Kane OCC³⁷ of the MAR has revealed a ~300 m thick, low-velocity layer at ~1.5 km depth from the seafloor interpreted as serpentinized peridotite underlain by a thick gabbroic body, and a deeper ~200 m thick, high-velocity layer interpreted as a gabbroic sill embedded in serpentinized peridotite. The P-wave velocity of serpentinized peridotite can have a wide range depending on the degree of serpentinization, and the velocity in the layered structure may well be within this range⁴¹. However, the tomographic results¹⁴ in our
<figure><img src="image_4.png" /><figcaption>**Fig. 2** | Seismic P-wave velocity models of the oceanic crust. **a–c**, FWI<sup>28</sup> velocity model starting from the tomographic seismic velocity model<sup>14</sup> (**a**), the velocity anomaly (the difference between the velocity models from FWI and tomography; **b** and the vertical velocity gradient (the derivative of velocity model from FWI with respect to depth; **c**. Black triangles in **a** mark the OBS locations. The velocity contours are from 5 to 7 km s<sup>-1</sup> with an increment of 0.2 km s<sup>-1</sup>. The coloured parts in **b** and **c** start from the basement (the top of layer 2). The lithospheric age in **c** is calculated using a spreading rate of 16 mm yr<sup>-1</sup> (ref. <sup>26</sup>). The horizontal distance starts with 0 km at the MAR.</figcaption></figure>

area suggest that lower crustal velocity increases from ~6.6 km s⁻¹ at 2 km below the basement to ~7 km s⁻¹ just above the Moho, suggestive of a gabbroic origin. The mantle velocity just below the Moho is ~8 km s⁻¹, confirming the presence of mantle peridotite¹⁴. The presence of a sharp Moho rules out the presence of gabbro sills in uplifted mantle peridotite for explaining the observed lower crustal seismic layering.

It would be difficult to distinguish if the lower crustal layering is caused by the gabbro glacier model¹,¹⁵ or the sheeted sill model⁶,⁷ from the velocity models alone; however, the recent discoveries of stacked melt lenses in the lower crust at the fast-16–18 and intermediate-spreading ridges¹⁹–²¹ from seismic reflection studies support the lower crustal accretion by in situ melt intrusion and crystallization. For slow-spreading ridges, melt lenses have been imaged at the upper-lower crust boundary for melt-rich spreading centre segments³⁸,⁴², while our study region is supposed to lie in the cold mantle temperature regime¹⁴,⁴³. Although no secondary melt sills were observed in the lower crust, the existence of low velocity anomalies has been reported underneath slow-spreading ridges, suggesting the presence of partial melt (mush)²²–²⁴. The rough sea-floor bathymetry and the weak P-wave velocity contrasts between the host mush and the injected melt sills could cause the absence of lower crustal melt sill evidence from seismic reflection study.

The layered structure in the lower crust has also been observed in exposed ophiolites. In the Oman ophiolite, which is a representative of fast-spreading mid-ocean ridge, thin layers (centimetres to 100 m) of alternating strata rich in mafic minerals (olivine (Ol) and clinopyroxene (Cpx)) at the base and plagioclase (Pl) at the top have been observed⁶,⁷. The geochemical analysis of these layers suggests the presence of secondary sills and that the lower crust could be formed by cooling and crystallization of melt sills in situ (the sheeted sill model)⁶,⁷. In the Bay of Islands ophiolite, a representative of slow-spreading environment, layers with thickness of several hundred metres to a kilometre have been observed⁴⁴.

A more than 1,500 m of the lower crustal gabbro was drilled during ODP Hole 735B⁴⁵ from an ultra-slow-spreading environment containing several layered units of >250 m thickness⁴⁵,⁴⁶. Analysis of a geochemical unit in the lower gabbro suggests that it constitutes a single magma reservoir but is separated in two parts⁴⁶. The lower part is formed by a stack of repeated recharge of primitive thin melt sills, whereas the upper part contains a homogeneous evolved magma mush formed by upward reactive porous flow, progressive differentiation and accumulation⁴⁶,⁴⁷. This process would lead to a layering in the lower crust, containing Ol-rich gabbro and troctolites at the base with a distinct boundary separating more evolved Ol and Ol-bearing gabbro with decreasing Ol and increasing Pl.
<figure><img src="image_5.png" /><figcaption>**Fig. 3** | Observed and synthetic seismic waveform data from OBS 57. **a,** The field data (black) and the synthetic waveform (red) using the tomographic model (Extended Data Fig. 1). **b,** The field data (black) and the synthetic waveform (red) using the true-amplitude FWI model (Fig. 2). Seismic waveforms from more OBSs can be found in Extended Data Figs. 3 and 4. A reduced travel time of 7 km s⁻¹ velocity was applied to both the field and synthetic data.</figcaption></figure>

Observations from Integrated Ocean Drilling Program (IODP)/ ODP drilling and ophiolites³⁹,⁴⁵⁻⁴⁷ indicate that the lower crustal gabbroic rocks are mainly composed of Ol, Cpx and Pl, with seismic velocities $V_{Ol} > V_{Cpx} > V_{Pl}$ (Extended Data Fig. 10a), where V represents the velocity. Rocks rich in Ol and Cpx, indicative of more primitive melt, would have higher velocities, whereas more evolved rocks rich in Pl would have relatively lower velocities³¹,⁴⁵. Using the Voigt-Reuss-Hill averaging method⁴⁸ (Methods), we computed the P-wave velocities for different gabbro compositions from Hole 735B and found that the velocity of the gabbro varies from 6.6 to 7.1 km s⁻¹ (Extended Data Fig. 10b) for 10–20% changes in Ol/Cpx and Pl contents. Therefore, high-velocity layers may contain gabbroic rocks with relatively high Ol/Cpx concentration (usually >10%) such as Ol and Ol-rich gabbro; low-velocity layers may be rich in rocks with low Ol/Cpx (<10%) and high Pl contents, such as Ol-bearing gabbro.

The observed lower crustal layering of frozen melt extends 5–15 km horizontal distance along the profile (Fig. 2c), suggesting that the crustal accretion process is stable for 300,000–800,000 years. A sustained melt supply and a relatively stable crustal accretion process can indeed be a general phenomenon at slow-spreading ridges, especially when the process is not interrupted by the presence of detachment faults⁴, which seems to be the case in our study region. Furthermore, the bathymetric data show simple abyssal hill fabrics indicating a normal mode of seafloor spreading. A robust and stable magma supply is also consistent with the high-velocity lower crust. It is also possible that these interpreted frozen sills are thinner and shorter (Supplementary Figs. 15 and 17), suggesting a more variable magma supply, but this is beyond data resolution. Studies at the Kane³⁷ and Rainbow⁴⁹,⁵⁰ Massif OCCs of the MAR, both of low-magma-supply environment, have reported sill intrusions with shorter extension in the uplifted ultramafic rocks. Large normal faults (with vertical offsets >400 m) at the ridge axis might disrupt the continuity of layers in the lower crust.

**Model for melt sill injection and crystallization in the lower crust**

Figure 4 shows a schematic diagram about the geometry of the melt lens and the oceanic crust. Together with the overwhelming evidence for the sheeted sill model from recent geophysical¹⁶⁻²⁴ and ODP/IODP drilling⁴⁵⁻⁴⁷ studies, the discovery of lower crust stacked layering off-axis shows that in our study area the lower oceanic crust is more probably formed by in situ melt sill intrusion from mantle and crystallization than the gabbro glacier model, and not dominated by tectonic extension with emplaced mantle-derived peridotite. The alternate high and low velocities indicate the progressive extraction and assimilation of cyclically replenished melts. The high-velocity layers are probably formed from primitive melt intrusions at the base of a magma reservoir, whereas the low-velocity layers can be associated with the more evolved, upward-migrated melt residue (mush). In situ crystallization requires extensive seawater circulation down to the Moho depth along the sides of the crystal-rich mush zone beneath the magma chambers for cooling⁵, which can be provided through well-developed faults/fractures in the slow-spreading environment¹⁶. Moreover, the hydrothermal activities and the magmatic reaction with host rocks could lead to remelting and assimilation⁴⁰,⁴⁶,⁴⁷, altering the petrological constituents of igneous rocks, hence changing the rock properties further.

Our results provide the quantitative seismic velocity model for a layered lower oceanic crust away from the spreading centre. We suggest that the oceanic lower crust is generally composed of high- and low-velocity layers, and that the magmatic oceanic lower crust at slow-spreading ridges is probably formed by in situ cooling and crystallization of cyclic magma upwelling from the upper mantle. The uppermost low-velocity layer could be associated with hydrothermal alteration. These results demonstrate the capability of the advanced FWI technique for determining high-resolution quantitative velocity models of oceanic crust and could be applied to other
<figure><img src="image_6.png" /></figure>

**Fig. 4** | Modelled layering in the lower oceanic crust. **a**, 1D velocity-depth profiles at the distance of 182 km from the MAR, derived from tomography (blue) and FWI (black). Vp is the P-wave velocity of seismic waves. The thin dashed black line is the interpreted boundary of layer 2A/2B from the FWI model; the bold dashed blue and black lines indicate the boundaries of layer 2/3 from the tomographic and the FWI models, respectively. We invert the velocity model from FWI down to a maximum of 9 km depth (-5.4 km from basement at distance 182 km), because of potential Moho reflection interference for greater depths (Supplementary Fig. 22). **b**, A schematic diagram illustrating the oceanic crustal structure at the MAR and away from the ridge axis. The black balls in layer 2A and stripes in layer 2B refer to pillow lava and basaltic dykes, respectively. The blue ellipsoids refer to hydrothermal alteration above the roof of the AML that could contribute to the top low-velocity layer. The red vertically elongated ellipsoids suggest magma from mantle upwelling. The light brown layers in the lower crust refer to the low-velocity layers from FWI. The imaged layers in the lower crust may comprise multiple thinner layers as indicated by the darker brown patches, but are beyond data resolution.

similar data from different oceanic settings, to shed light on the formation and evolution of oceanic lithosphere.

## Online content

Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41561-022-00963-w.

Received: 26 January 2021; Accepted: 6 May 2022;
Published online: 13 June 2022

## References

1. Morgan, J. P. & Chen, Y. J. The genesis of oceanic crust: magma injection, hydrothermal circulation, and crustal flow. *J. Geophys. Res. Solid Earth* **98**, 6283-6297 (1993).

2. Cannat, M. et al. Thin crust, ultramafic exposures, and rugged faulting patterns at the Mid-Atlantic Ridge (22–24 N). *Geology* **23**, 49–52 (1995).

3. Christeson, G. L., Goff, J. A. & Reece, R. S. Synthesis of oceanic crustal structure from two-dimensional seismic profiles. *Rev. Geophys.* **57**, 504–529 (2019).

4. Cannat, M. Emplacement of mantle rocks in the seafloor at mid-ocean ridges. *J. Geophys. Res. Solid Earth* **98**, 4163–4172 (1993).

5. MacIennan, J., Hulme, T. & Singh, S. C. Thermal models of oceanic crustal accretion: linking geophysical, geological and petrological observations. *Geochem. Geophys. Geosyst.* **5**, Q02F25 (2004).

6. Boudier, F., Nicolas, A. & Ildefonse, B. Magma chambers in the Oman ophiolite: fed from the top and the bottom. *Earth Planet. Sci. Lett.* **144**, 239–250 (1996).

7. Kelemen, P. B., Koga, K. & Shimizu, N. Geochemistry of gabbro sills in the crust-mantle transition zone of the Oman ophiolite: implications for the origin of the oceanic lower crust. *Earth Planet. Sci. Lett.* **146**, 475–488 (1997).

8. Wilson, D. S. et al. Drilling to gabbro in intact ocean crust. *Science* **312**, 1016–1020 (2006).

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

10. Toomey, D. R., Solomon, S. C. & Purdy, G. M. Tomographic imaging of the shallow crustal structure of the East Pacific Rise at 9°30′ N. *J. Geophys. Res. Solid Earth* **99**, 24135–24157 (1994).

11. Grevemeyer, I., Ranero, C. R. & Ivandic, M. Structure of oceanic crust and serpentinization at subduction trenches. *Geosphere* **14**, 395–418 (2018).

12. Detrick, R. S. et al. Seismic structure of the southern East Pacific Rise. *Science* **259**, 499–503 (1993).

13. Harding, A. J., Kent, G. M. & Orcutt, J. A. A multichannel seismic investigation of upper crustal structure at 9° N on the East Pacific Rise: implications for crustal accretion. *J. Geophys. Res. Solid Earth* **98**, 13925–13944 (1993).

14. Vaddineni, V. A., Singh, S. C., Grevemeyer, I., Audhkhasi, P. & Papenberg, C. Evolution of the crustal and upper mantle seismic structure from 0–27 Ma in the equatorial Atlantic Ocean at 2° 43'S. *J. Geophys. Res. Solid Earth* **126**, e2020JB021390 (2021).

15. Quick, J. E. & Denlinger, R. P. Ductile deformation and the origin of layered gabbro in ophiolites. *J. Geophys. Res. Solid Earth* **98**, 14015–14027 (1993).

16. Marjanović, M. et al. A multi-sill magma plumbing system beneath the axis of the East Pacific Rise. *Nat. Geosci.* **7**, 825–829 (2014).

17. Arnulf, A. F., Singh, S. C. & Pye, J. W. Seismic evidence of a complex multi-lens melt reservoir beneath the 9°N overlapping spreading center at the East Pacific Rise. *Geophys. Res. Lett.* **41**, 6109–6115 (2014).

18. Canales, J. P. et al. Network of off-axis melt bodies at the East Pacific Rise. *Nat. Geosci.* **5**, 279–283 (2012).

19. Canales, J. P., Nedimović, M. R., Kent, G. M., Carbotte, S. M. & Detrick, R. S. Seismic reflection images of a near-axis melt sill within the lower crust at the Juan de Fuca ridge. *Nature* **460**, 89–93 (2009).

20. Carbotte, S. M. et al. Stacked sills forming a deep melt-mush feeder conduit beneath Axial Seamount. *Geology* **48**, 693–697 (2020).

21. Carbotte, S. M. et al. Stacked magma lenses beneath mid-ocean ridges: insights from new seismic observations and synthesis with prior geophysical and geologic findings. *J. Geophys. Res. Solid Earth* **126**, e2020JB021434 (2021).

22. Canales, J. P., Collins, J. A., Escartin, J. & Detrick, R. S. Seismic structure across the rift valley of the Mid-Atlantic Ridge at 23 20' (MARK area): implications for crustal accretion processes at slow spreading ridges. *J. Geophys. Res. Solid Earth* **105**, 28411–28425 (2000).

23. Seher, T. et al. Crustal velocity structure of the Lucky Strike segment of the Mid-Atlantic Ridge at 37°N from seismic refraction measurements. *J. Geophys. Res. Solid Earth* **115**, B03103 (2010).

24. Dunn, R. A., Lekić, V., Detrick, R. S. & Toomey, D. R. Three-dimensional seismic structure of the Mid-Atlantic Ridge (35° N): evidence for focused melt supply and lower crustal dike injection. *J. Geophys. Res. Solid Earth* **110**, B09101 (2005).

25. Grevemeyer, I., Singh, S. C. & Papenberg, C. Ocean Bottom Seismometer and Ocean Bottom Hydrophone Seismic Refraction and Wide Angle Data from Profile P01 of Maria S. Merian Cruise MSM69 with Links to sgy Data Files (PANGAEA, 2020); https://doi.pangaea.de/10.1594/PANGAEA.914912
26. Müller, R. D., Sdrolias, M., Gaina, C. & Roest, W. R. Age, spreading rates, and spreading asymmetry of the world's ocean crust. *Geochem. Geophys. Geosyst.* **9**, Q04006 (2008).

27. Wang, Z., Singh, S. C. & Noble, M. True-amplitude versus trace-normalized full waveform inversion. *Geophys. J. Int.* **220**, 1421–1435 (2020).

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

29. Tarantola, A. Inversion of seismic reflection data in the acoustic approximation. *Geophys* **49**, 1259–1266 (1984).

30. Carlson, R. L. The effect of hydrothermal alteration on the seismic structure of the upper oceanic crust: evidence from Holes 504B and 1256D. *Geochem. Geophys. Geosyst.* **12**, Q09013 (2011).

31. Carlson, R. L. & Jay Miller, D. Influence of pressure and mineralogy on seismic velocities in oceanic gabbros: implications for the composition and state of the lower oceanic crust. *J. Geophys. Res. Solid Earth* **109**, B09205 (2004).

32. Audhkhasi, P. & Singh, S. C. Seismic structure of the upper crust from 0–75 Ma in the equatorial Atlantic Ocean on the African plate using ultralong offset seismic data. *Geochem. Geophys. Geosyst.* **20**, 6140–6162 (2019).

33. Teagle, D. A. H. et al. Superfast spreading rate crust 3. *Proc. Integr. Ocean. Drill. Program* Vol. 312, 1–168 (2006).

34. Alt, J. C. et al. Hydrothermal alteration of a section of upper oceanic crust in the eastern equatorial Pacific: a synthesis of results from Site 504 (DSDP LEGS 69, 70, and 83, and ODP LEGS 111, 137,140, and 148). *Proc. Ocean Drill. Program* **148**, 417–434 (1996).

35. Singh, S. C., Collier, J. S., Harding, A. J., Kent, G. M. & Orcutt, J. A. Seismic evidence for a hydrothermal layer above the solid roof of the axial magma chamber at the southern East Pacific Rise. *Geology* **27**, 219–222 (1999).

36. Swift, S., Reichow, M., Tikku, A., Tominaga, M. & Gilbert, L. Velocity structure of upper ocean crust at Ocean Drilling Program Site 1256. *Geochem. Geophys. Geosyst.* **9**, Q10O13 (2008).

37. Canales, J. P. 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.* **37**, e2010GL044412 (2010).

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

39. Kelemen, P. B., Matter, J. M., Teagle, D. A. H., Coggon, J. A. & The Oman Drilling Project Science Team. Site GT3: sheeted dike to gabbro transition. In *Proc. of the Oman Drilling Project* (ed. The Oman Drilling Project Science Team) 1-152 (International Ocean Discovery Program, 2020); https://doi.org/10.14379/OmanDP.proc

40. France, L. et al. Quantifying the axial magma lens dynamics at the roof of oceanic magma reservoirs (dike/gabbro transition): Oman Drilling Project GT3 site survey. *J. Geophys. Res. Solid Earth* **126**, 2020JB021496 (2021).

41. Miller, D. J. & Christensen, N. I. Seismic velocities of lower crustal and upper mantle rocks from the slow-spreading Mid-Atlantic Ridge, south of the Kane Transform Zone (MARK). *Proc. Ocean Drill. Program* **153**, 437–456 (1997).

42. Sinha, M. C. et al. Magmatic processes at slow spreading ridges: implications of the RAMESSES experiment at 57 45'N on the Mid-Atlantic Ridge. *Geophys. J. Int.* **135**, 731–745 (1998).

43. Bonatti, E. et al. Transform migration and vertical tectonics at the Romanche fracture zone, equatorial Atlantic. *J. Geophys. Res. Solid Earth* **99**, 21779–21802 (1994).

44. Casey, J. F. & Karson, J. A. Magma chamber profiles from the Bay of Islands ophiolite complex. *Nature* **292**, 295–301 (1981).

45. Iturrino, G. J., Ildefonse, B. & Boitnott, G. Velocity structure of the lower oceanic crust: results from Hole 735B, Atlantic II Fracture Zone. In *Proc. ODP Sci. Results* (eds Natland, J. H., Dick H. J. B., Miller, D. J. & Von Herzen, R. P.) Vol. 176, 1–71 (Ocean Drilling Program, 2001).

46. Boulanger, M. et al. Magma reservoir formation and evolution at a slow-spreading center (Atlantis Bank, Southwest Indian Ridge). *Front. Earth Sci.* **8**, 554598 (2020).

47. Lissenberg, C. J., MacLeod, C. J., Howard, K. A. & Godard, M. Pervasive reactive melt migration through fast-spreading lower oceanic crust (Hess Deep, equatorial Pacific Ocean). *Earth Planet. Sci. Lett.* **361**, 436–447 (2013).

48. Hill, R. The elastic behaviour of a crystalline aggregate. *Proc. Phys. Soc., A.* **65**, 349 (1952).

49. Canales, J. P., Dunn, R. A., Arai, R. & Sohn, R. A. Seismic imaging of magma sills beneath an ultramafic-hosted hydrothermal system. *Geology* **45**, 451–454 (2017).

50. Dunn, R. A., Arai, R., Eason, D. E., Canales, J. P. & Sohn, R. A. Three-dimensional seismic structure of the Mid-Atlantic Ridge: an investigation of tectonic, magmatic, and hydrothermal processes in the Rainbow area. *J. Geophys. Res. Solid Earth* **122**, 9580–9602 (2017).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long

as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© Crown 2022
Methods

**Data acquisition.** The field OBS data were acquired during the LITHOS-iLAB cruise in 2017 onboard the Germain RV *Marie S. Merian* to study the upper lithosphere from 0 to 50 Ma<sup>14,51</sup>. A total of 71 instruments consisting of 55 OBSs and 16 ocean bottom hydrophones were deployed along a 1,100-km-long transect with a variable spacing of 10–20 km to record wide-angle refractions and reflection arrivals. All OBSs were equipped with a hydrophone (measuring pressure) and three geophones (measuring vertical and horizontal displacements). The data were sampled at 250 Hz. The active seismic source used in the survey comprised 6 G-gun clusters (12 guns) configured as two sub-arrays with a total volume of 5,440 cubic inches, which was towed at 7.5 m in depth and fired every ~400 m along the profile. We used a part of these data (8 OBSs, hydrophone component) where the OBS spacing was dominantly ~10 km covering a young (7–12 Ma) oceanic crust (Fig. 1).

**Data pre-processing.** We limit the data pre-processing to a minimum to keep the waveform information, including a zero-phase bandpass filtering of 3.5–10 Hz and 3D to 2D transformation<sup>28,29</sup> (multiplying the amplitudes of field data by $\sqrt{t}$, where $t$ is the two-way travel time, and convolving with $1/\sqrt{t}$) because 2D elastic wave equation modelling was used for simulating seismic data recorded in a 3D Earth. A predictive gapped deconvolution was applied to suppress the bubble effects from the air gun sources; we used the supef function from processing software Seismic Unix<sup>52</sup>, with value 0.2 s for parameter minlag and 0.6 s for maxlag. Clear crustal turning waves (Pg) with high signal-to-noise ratio can be observed starting from offsets (source–receiver distance) of ±6–8 km, up to offsets of ±20–26 km, followed by wide-aperture Moho reflection (PmP) and upper mantle refraction (Pn) at the far offsets<sup>14</sup> (Supplementary Fig. 1).

**Travel-time tomography.** Travel-time tomography<sup>14,53</sup> was first applied for estimating a P-wave velocity model (Extended Data Fig. 1). Tomography was carried out through a linearized iterative approach. At each iteration, a ray-tracing algorithm with a hybrid of graph (shortest path) method and ray bending was used for forward modelling of travel times<sup>53</sup>; the model update was obtained by least-squares penalties on the data misfit, together with smoothing and damping for regularization<sup>53</sup>.

The velocity model from travel-time tomography (Extended Data Fig. 1) gives very good travel-time fit for the first arrivals<sup>14</sup>. Since the travel time is mainly sensitive to large-scale velocity structures, the tomographic velocity model contains few details for the oceanic crust. The first arrival, Pg, rays penetrate down to ~3 km below the seafloor, thus the velocity in the lower crust is mainly determined using the PmP arrivals, and can be poorly constrained because of the trade-off between the lower crustal velocity and the Moho depth<sup>14</sup>. Details of the travel-time tomography can be found in ref. <sup>14</sup>.

**FWI.** FWI is the current state-of-the-art technique for high-resolution sub-surface quantitative imaging<sup>28,54,55</sup>. Unlike tomography only using travel time<sup>53</sup>, the FWI is based on minimizing the difference between the observed and synthetic seismic waveforms, with a numerical solution of the elastic wave equation for realistic simulation of seismic wave propagation<sup>28</sup>. For the numerical implementation, a gradient-based linearized inversion approach is used for updating the velocity model iteratively, with the gradients of the data misfit to model parameters efficiently calculated by the adjoint method from cross-correlation of the forward and adjoint wavefields<sup>28,29</sup>. We used time-domain staggered-grid finite-difference<sup>54</sup> with fourth-order spatial and second-order temporal accuracy for solving the elastic wave equation in the stress and particle-velocity formulation. The convolutionally perfectly matched layer absorbing boundary condition<sup>55</sup> was applied at the model boundaries.

FWI can provide high-resolution quantitative Earth internal models for physical properties and their vertical differentiations. Sensitivity studies<sup>56,57</sup> suggest that the travel-time data provide information on a scale greater than approximately five times the dominant wavelength, whereas the FWI provides updates of between a quarter of wavelength to a wavelength. Compared with ray tracing for the Pg arrivals mainly within the upper crust<sup>14</sup>, FWI can update the lower crustal model with wide-angle reflections/diffractions from lower crustal anomalies being accurately modelled<sup>58</sup>. The application of FWI is mainly in the upper oceanic crust<sup>17,28,35,37,59</sup>, with very few exceptions down to 3–5 km depth<sup>60</sup>.

Although the tomographic inversion has converged with good travel-time fit<sup>14</sup>, large waveform difference can be observed between the field data and the synthetic seismograms from wave equation modelling (Fig. 3a). We employed a multi-stage strategy for obtaining the final model. In the first stage of FWI, the tomographic model was used as the starting model. A trace-normalized FWI<sup>27</sup> was applied, with the misfit function (J) defined as:

$$J = \sum_{i=1}^{\text{Ns}} \sum_{j=1}^{\text{Nr}} s_{i,j} \left( \frac{s_{i,j}}{|s_{i,j}|} - \frac{d_{i,j}}{|d_{i,j}|} \right), \quad (1)$$

where s and d are synthetic and field data, respectively, $||||$ is the l-2 norm of a seismic trace (a time-series vector), Ns and Nr are the number of sources and receivers, and i and j are the indexes for the sources and receivers, respectively. This intermediate step was used to bridge the gap between tomography and the classic

FWI using 'true amplitude' waveform: the influence of amplitude-versus-offset is removed at the adjoint source by trace-by-trace normalization but the amplitude variation within each trace remains. An important point of the trace-normalized inversion is to determine the amplitude normalization factors between the synthetic and real data. In the second stage, we performed a true-amplitude FWI<sup>28,29</sup>, with the misfit function as:

$$J = \sum_{i=1}^{\text{Ns}} \sum_{j=1}^{\text{Nr}} \left| \frac{s_{i,j} - d_{i,j}}{\|s_{i,j}\|} \right|^2. \quad (2)$$

Source wavelet is important for seismic wave modelling in FWI and was estimated by stacking near-offset free-surface multiples of water waves<sup>51</sup> (Extended Data Fig. 2), thanks to a good separation between water waves and seafloor-related scatterings in their free-surface multiples. An amplitude scaling factor was then estimated by comparing the field and synthetic waveforms at near offsets.

For the inversion, we used the crustal Pg arrivals, because this part of the data has the most linear behaviour with sub-surface properties and the wide-aperture data is sensitive to both the upper and lower crustal structures<sup>38</sup>. We didn't apply FWI to the PmP arrivals, because of its strong nonlinearity around critical angles<sup>51</sup>. A time window of 0.5 s (Supplementary Fig. 2) was applied to the OBS gathers, by muting the data before 0.2 s and after 0.3 s of the picked Pg travel times, to reduce the influence of noise and to isolate Pg arrivals from the other seismic events including the PmP at far offsets. We updated the P-wave velocity model only, with the S-wave velocity model derived from the P-wave velocity using Brocher's<sup>61</sup> regression fit. Density was linked to the P-wave velocity based on the empirical relation in ref. 62 for velocity smaller than 2.2 km s<sup>-1</sup>, and from ref. 63 for higher velocities. We used a grid spacing of 20 m and the time step is 0.0012 s. Considering the sparse distribution of OBS, a Gaussian smoothing operator with 2 km horizontal and 0.4 km vertical lengths was applied to the gradient for regularization.

The inversion was carried out using a top-down approach: we inverted first for the relatively shallow structures, which were constrained by near-to-intermediate source–receiver offset ranges (from ±6–8 km up to ±15 km offsets) of crustal Pg arrivals (Supplementary Fig. 3), then far-offset arrivals were included for estimating the deeper crustal model. The updated velocity model from the prior stage was used as the starting model for the next stage of inversion. Supplementary Fig. 4 shows the seismograms are better aligned after trace-normalized inversion. Extended Data Figs. 3 and 4 show seismic waveform match between synthetic and field data has been much improved after true-amplitude FWI. Extended Data Fig. 5 shows the data misfits for all the 80 iterations. The data misfits after FWI were reduced to 30–60% for different OBS gathers (Extended Data Fig. 5b).

The velocity models from the trace-normalized FWI are shown in Extended Data Figs. 6 and 7, and those from the true-amplitude FWI are shown in Fig. 2 and Extended Data Fig. 8. To heighten the salient features of the results from FWI, besides the velocity models, we also show the velocity anomalies (the difference between the FWI and the tomographic models) as well as the vertical velocity gradients. For example, the velocity anomaly from the trace-normalized FWI highlights the large-scale velocity variations (Extended Data Figs. 6b and 7), which sharpens the upper–lower crust transition, whereas the vertical velocity gradient from true-amplitude FWI enhances the layered anomalies in the lower crust (Fig. 2c and Extended Data Fig. 8c). Models in Extended Data Figs. 7 and 8 have better horizontal continuity than those in Fig. 2 and Extended Data Fig. 6, because of a larger Gaussian horizontal smoothing (4 km compared with 2 km). We also performed FWI without the intermediate trace-normalized FWI step (Supplementary Fig. 5). The recovered models are generally consistent with those in Fig. 2, with relatively larger data misfit (Supplementary Fig. 6). Considering the data frequency is 3.5–10 Hz, the dominant wavelength is 700–1600 m for a velocity of 6.8 km s<sup>-1</sup>, indicating that the vertical resolution in the lower crust is ~400 m.

Extended Data Fig. 9 shows 1D velocity–depth profiles below basement between 112 km and 192 km distance along the profile from the tomographic (Extended Data Fig. 1) and FWI velocity models (Fig. 2). One can observe that there is a wider variation in velocity both in the upper and lower crust as compared with the tomographic results. Extended Data Fig. 9c,d shows a subset of two (1D) velocity–depth profiles, with the interpreted boundaries of Layers 2A/2B and 2/3.

We didn't consider velocity anisotropy. Studies<sup>23,64</sup> at the MAR indicate velocity anisotropy of 2–4% in the upper crust, which becomes nearly isotropic from 3 km below the seafloor. The observed anisotropy decreases fast off-axis over 5–10 km<sup>-1</sup>. Therefore, we consider the influence of P-wave anisotropy to be weak. If P-wave anisotropy is strong, the velocity model from isotropic FWI constrained by turning waves and wide-angle reflections will approximate to the horizontal velocity<sup>65</sup>. Using surface-wave data, ref. 66 observed a S-wave velocity anisotropy (4–5% radial anisotropy) in the lower crust, but its influence on the first-arrival P-waves should be negligible. Lower crustal anisotropy could be produced by thin layering<sup>66</sup>; however, we can only invert layering of a quarter of wavelength using limited offset range. We consider intrinsic attenuation at the lithospheric age of 7–12 Myr to be weak, and thus was not included for FWI. With notable attenuation and viscoelastic modelling, we expect the layered anomalies in the lower crust larger than those in Fig. 2 to generate stronger reflection accounting for dissipation loss.
**Resolution and uncertainty synthetic studies.** To gain confidence in our results, especially for the layered structures in the lower crust, we performed extensive numerical tests using synthetic and field data. When performing synthetic forward modelling and FWI, we used the same source and receiver positions and frequency ranges of the data as in the actual observation, and the same inversion parameters. The same data windowing was applied using Pg travel-time picks from the field data. Except for testing smoothing operators, a Gaussian smoothing with 2 km horizontal and 0.4 km vertical lengths was applied for model update.

We modified the velocity model from FWI starting from 6 km depth to the bottom of the model (Supplementary Fig. 7), by replacing the layering structures with the corresponding trace-normalized FWI model to see if the layered structures in the lower crust are required by the data and not due to smearing of the structures in the upper crust.

We first performed synthetic modelling by adding alternate ±200 m s⁻¹ velocity anomalies in the lower crust, to understand the footprint of lower crustal layering in the wide-angle OBS data. Supplementary Fig. 8 shows the added velocity perturbations and the corresponding synthetic seismograms. We found that reflections from the layered anomalies will interfere and merge with the crustal Pg arrivals, causing amplitude and waveform variations at large offsets. Reflections from the lower crust also appear at near offsets but at a later time, tailing the Pg arrivals with very weak amplitudes. This test suggests that crustal Pg arrivals contain information about the lower crust from interference of reflections originating in the lower crust.

Then we compared the field data and synthetic seismograms using models with and without layering in the lower crust; larger waveform misfit can be observed (Supplementary Figs. 9 and 10) and quantified (Supplementary Fig. 11a,b), especially for the far (larger than 15 km; Supplementary Fig. 11b) offsets for models without layering. We used the modified model as the starting model to run FWI. The inverted model is shown in Supplementary Fig. 12 with remarkable similarity to Fig. 2. In addition, we also ran FWI for the field data in which the velocity structure at depths greater than 2.5 km from basement was fixed, so the inversion was forced to search for updates in the upper crust to explain the data. The final models are shown in Supplementary Fig. 13, with the main difference of the upper crust in Fig. 2 being a high velocity anomaly around 180 km distance. Waveform misfit is larger than updating the whole model, especially for large offsets (Supplementary Fig. 11c,d). These tests suggest that the layered structure in the lower crust is required for explaining the data.

We tested the influence of the smoothing operators for FWI. First, we tested the horizontal Gaussian smoothing lengths. We added layered anomalies of velocity ±200 m s⁻¹ with 400 m thickness and different horizontal extensions (Supplementary Fig. 14) to the tomographic model (Extended Data Fig. 1). The perturbed velocity model was used for generating the 'observed' data, and the unperturbed tomographic model was used as starting model for FWI to see how well these perturbations can be recovered. From Supplementary Fig. 15, we found that inversion results with 0.5 km and 1 km horizontal smoothing contain significant 'smiling' artefacts. Inversion discontinuities can also be observed with small horizontal smoothing and extensive anomalies. Conversely, large horizontal smoothing length may cause anomaly size to be overestimated. With these results and considering that the Fresnel zone for a dominant frequency of 5–6 Hz and 6.8 km s⁻¹ velocity at 4 km depth from basement is ~1.8 km, we chose 2 km as the optimal horizontal smoothing length. Next we ran inversion with different vertical smoothing lengths of 200 m, 400 m, 800 m and 1.6 km, respectively. The results shown in Supplementary Fig. 16a,b are very similar, indicating that vertical smoothing length smaller than FWI vertical resolution (~400 m) will not improve the results. Smoothing length larger than a quarter wavelength will smear the inverted anomalies (Supplementary Fig. 16c). With the vertical smoothing length being the wavelength, the layered anomalies cannot be reconstructed (Supplementary Fig. 16d). Therefore we chose 0.4 km for vertical smoothing.

We also evaluated the performance of FWI when the anomaly thickness in the lower crust is 100 m, much smaller than data resolution. Supplementary Fig. 17 shows that thin layers will be inverted as larger thickness (Supplementary Fig. 17e,h). If there are multiple layers with separation smaller than FWI resolution, they may not be recovered properly (Supplementary Fig. 17f,g).

Chequerboard tests were performed for estimating the size of the anomaly that could be recovered using our FWI settings. We generated synthetic models by adding 8% positive and negative Gaussian-shaped velocity anomalies (10 km long and 1 km thick, Supplementary Fig. 18a; 10 km long and 0.5 km thick, Supplementary Fig. 19a) to the tomographic velocity model. Another test with random anomalies of thickness 400 m but variable horizontal extension is shown in Supplementary Fig. 20. Most of the anomalies are well resolved with good match to the true models (Supplementary Figs. 18–20). Different anomalies with distance <2 km were imaged as one (Supplementary Fig. 20b). Results suggest that the OBS data we used have constraints for anomalies at different depths, including the uppermost and lower crust. We also observe inversion artefacts accompanying real structures but with opposite values arising from sidelobes of seismic wavefront, therefore some of the anomalies in the inverted model could be artefacts.

We tested if the layering in the lower crust could be introduced by the potential presence of the Moho (PmP) reflections. We used the velocity model with Moho (Supplementary Fig. 21a)¹⁴ for generating the 'observed' data. For FWI, the velocity

model after removing the Moho and upper mantle (Supplementary Fig. 21b) was used as the starting model. We observed weaker anomalies around the crustal base in the inverted model (Supplementary Fig. 22), but their influence for the observed crustal layering is limited. These tests confirm that the layered structure in the lower crust, especially that away from the crustal base, represents real features of the sub-surface; the Moho reflection is unlikely to be attributed to the layered structures in the lower crust. We also test the influence of muting time window size for FWI. Supplementary Fig. 23 contains the inverted velocity models for time windows of 0.5 s (Fig. 2), 0.8 s, 1.2 s and 2.0 s, respectively; the velocity models are similar, except for some small differences mainly in the lower portion of the lower crust.

**Voigt–Reuss–Hill averaging for rock property modelling.** To shed light on the types of rock that can produce velocities comparable to those estimated from the FWI, we estimated velocities for rocks commonly present in the lower crust. Voigt–Reuss–Hill averaging⁴⁸ is a common and simple approach for computing the effective elastic moduli of rocks from their mineral constituents. Evidence¹⁹,⁴⁵⁻⁴⁷ from IODP/ODP drilling and ophiolite show that the gabbroic rocks of the lower crust are mainly composed of Ol, Cpx and Pl. Therefore, we calculated the P-wave velocities using different combinations of Ol, Cpx and Pl volume fractions. The elastic moduli and densities⁴⁹ of the three mineral components are shown in Extended Data Fig. 10a. We expect the velocities to increase with increasing Ol or Cpx, and decrease with increasing Pl. By varying the volume fractions of different components, we found a velocity variation of 200–400 m s⁻¹ (±100–200 m s⁻¹; Extended Data Fig. 10b), similar to those observed from the FWI. Pyroxene may play a more important role in increasing the velocities of gabbro than Ol, because the properties of Ol can get more easily altered, leading to low effective densities and elastic moduli³¹.

## Data availability

The raw OBS data used for this study is stored at the PANGAEA data centre https://doi.pangaea.de/10.1594/PANGAEA.914912. The derived velocity models and travel-time pickups can be accessed at https://doi.org/10.5281/zenodo.4390552.

## Code availability

The code for travel-time tomography and FWI can be accessed upon reasonable request from S.C.S. (singh@ipgp.fr) or P.G. (peng.guo@csiro.au).

## References

51. Guo, P. et al. Nonlinear full waveform inversion of wide-aperture OBS data for Moho structure using a trans-dimensional Bayesian method. *Geophys. J.* **224**, 1056–1078 (2021).

52. Stockwell, J. W. & Cohen, J. K. *The New SU User's Manual* (Center for Wave Phenomena, Colorado School of Mines, 2002).

53. Van Avendonk, H. J., Harding, A. J., Orcutt, J. A. & McClain, J. S. A two-dimensional tomographic study of the Clipperton transform fault. *J. Geophys. Res. Solid Earth* **103**, 17885–17899 (1998).

54. Virieux, J. P-SV wave propagation in heterogeneous media: velocity-stress finite-difference method. *Geophys* **51**, 889–901 (1986).

55. Komatitsch, D. & Martin, R. An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation. *Geophysics* **72**, SM155–SM167 (2007).

56. Jannane, M. et al. Wavelengths of Earth structures that can be resolved from seismic reflection data. *Geophysics* **54**, 906–910 (1989).

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

58. Ingale, V. & Singh, S. C. Insights from synthetic seismogram modelling study of oceanic lower crust and Moho transition zone. *Tectonophysics* **816**, 229030 (2021).

59. Xu, M., Zhao, X. & Canales, J. P. Structural variability within the Kane oceanic core complex from full waveform inversion and reverse time migration of streamer data. *Geophys. Res. Lett.* **47**, e2020GL087405 (2020).

60. Davy, R. G. et al. Resolving the fine-scale velocity structure of continental hyperextension at the Deep Galicia Margin using full-waveform inversion. *Geophys. J. Int.* **212**, 244–263 (2018).

61. Brocher, T. M. Empirical relations between elastic wavespeeds and density in the Earth's crust. *Bull. Seismol. Soc. Am.* **95**, 2081–2092 (2005).

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

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

64. Barclay, A. H., Toomey, D. R. & Solomon, S. C. Seismic structure and crustal magmatism at the Mid-Atlantic Ridge, 35° N. *J. Geophys. Res. Solid Earth* **103**, 17827–17844 (1998).

65. Prieux, V. et al. On the footprint of anisotropy on isotropic full waveform inversion: the Valhall case study. *Geophys. J. Int.* **187**, 1495–1515 (2011).

66. Russell, J. B. et al. High-resolution constraints on Pacific upper mantle petrofabric inferred from surface-wave anisotropy. *J. Geophys. Res. Solid Earth* **124**, 631–657 (2019).
67. Sobolev, S. V. & Babeyko, A. Y. Modeling of mineralogical composition, density and elastic wave velocities in anhydrous magmatic rocks. *Surv. Geophys.* **15**, 515–544 (1994).

## Acknowledgements

The OBS data were acquired during the LITHOS-iLAB experiment onboard the German RV *Maria S. Merian* in 2017. We thank L. France for the discussion on ODP Hole 735B. This project was funded by the European Research Council Advanced Grant agreement no. 339442 TransAtlanticiLAB. The cruise with the ID MSM69 was funded by the German Science Foundation (DFG). This research was funded by the Deep Earth Imaging Future Science Platform, CSIRO. This work was supported by resources provided by the Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia. We also thank the computational resources supported by S-CAPAD at IPGP. This work is IPGP contribution no. 4264.

## Author contributions

P.G. performed data processing and FWI of the OBS data, analysed the results, and wrote the paper; S.C.S. participated in the OBS data acquisition, designed and supervised the

project, interpreted the results, and wrote the paper; V.A.V. participated in the OBS data acquisition, and performed data processing and travel-time tomography; I.G. led the OBS data acquisition and discussed the results; E.S. discussed the results.

## Competing interests

The authors declare no competing interests.

## Additional information

**Extended data** is available for this paper at https://doi.org/10.1038/s41561-022-00963-w.

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41561-022-00963-w.

**Correspondence and requests for materials** should be addressed to Peng Guo.

**Peer review information** *Nature Geoscience* thanks Timothy Minshull, Gaye Bayrakci, Ryuta Arai and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editors: Simon Harold and Louise Hawkins, in collaboration with the *Nature Geoscience* team.

**Reprints and permissions information** is available at www.nature.com/reprints.
<figure><img src="image_8.png" /><figcaption>**Extended Data Fig. 1** | The tomographic seismic P-wave velocity model from travel-time tomography. It serves as the starting model for the FWI.</figcaption></figure>
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Extended Data Fig. 2 | See next page for caption.</figcaption></figure>
**Extended Data Fig. 2 | Source wavelet and comparison of synthetic and field data.** In (a), the left five black wiggles show the aligned free-surface multiples of water waves; the dashed lines show water waves before alignment; the green line is the stacked results of the black wiggles; the red line is the stacked results filtered into frequency of 3.5–10 Hz. Because of the reflection coefficient is -1 at free surface, the red wiggle needs to be multiplied with -1 for source wavelet. (b) shows comparison of field and synthetic data, which demonstrates good waveform match, suggesting a successful estimation of the source wavelet. For more details, refers to Guo et al. 2021.
<figure><img src="image_11.png" /></figure>

**Extended Data Fig. 3** | Observed and synthetic seismic waveform data from the tomographic model. The recorded field data (black) and the synthetic waveform (red) using the tomographic model (Extended Data Fig. 1) for OBS 55–59. A reduced travel time of 7 km/s velocity was applied to both the field and synthetic data. A scalar weighting factor (Offset/6) was multiplied for each trace to boost amplitude at large offsets for visualisation.
<figure><img src="image_12.png" /></figure>

**Extended Data Fig. 4** | Observed and synthetic seismic waveform data from the FWI model. The recorded field data (black) and the synthetic waveform (red) using the FWI model (Fig. 2) for OBS 55–59. A reduced travel time of 7 km/s velocity was applied to both the field and synthetic data. A scalar weighting factor (Offset/6) was multiplied for each trace to boost amplitude at large offsets for visualisation.
<figure><img src="image_13.png" /></figure>

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

**Extended Data Fig. 5** | Data misfits from FWI. (a) Sum of data misfit of all the 8 OBS gathers for the 80 iterations of FWI. Inversion starts with trace-normalised FWI, followed by true-amplitude of FWI of the first arrivals at the near offsets, the far offsets, and the full offsets. Note that the misfits at the beginning of each stage were normalised to 1. (b) Data misfit for each of the 8 OBS gathers from the tomographic model and the final FWI models. FWI models (2 km) and (4 km) are FWI with Gaussian smoothing of 2 km (Fig. 2) and 4 km (Extended Data Fig. 8) along the horizontal direction for velocity update, respectively. The waveform data misfit was calculated using Eq. 2 in the Method section and has been normalised.
<figure><img src="image_15.png" /></figure>

**Extended Data Fig. 6** | Seismic P-wave velocity models of the oceanic crust from the traced-normalised FWI. A Gaussian smoothing operator (smoothing lengths of 2 km and 0.4 km in the horizontal and vertical directions, respectively) was applied for regularising the velocity update. (a) The velocity model from the trace-normalised FWI, (b) the velocity anomaly (the difference between the velocity models from the FWI and the tomography), and (c) the vertical velocity gradient (the derivative of velocity with respect to depth). The lithospheric age is calculated using a spreading rate of 16 mm/ year²⁶. Black triangles in (a) mark the locations of OBS. The velocity contours in (a) are from 5 to 7 km/s with an increment of 0.2 km/s. The coloured parts in (b) and (c) start from the basement (the top of Layer 2). The horizontal distance starts with 0 km from the MAR.
<figure><img src="image_16.png" /></figure>

**Extended Data Fig. 7** | The velocity models similar to that in Extended Data Fig. 6 but with a larger Gaussian smoothing. A Gaussian smoothing operator with 4 km in the horizontal and 0.4 km in the vertical directions was used for regularising velocity update, and the rest of the parameters are the same as in Extended Data Fig. 6.
<figure><img src="image_17.png" /></figure>

**Extended Data Fig. 8** | The velocity models similar to that in Fig. 2 but with a larger Gaussian smoothing. A Gaussian smoothing operator with 4 km in the horizontal and 0.4 km in the vertical directions was used for regularising velocity update, and the rest of the parameters are the same as in Fig. 2.
<figure><img src="image_18.png" /></figure>

**Extended Data Fig. 9 | Velocity-depth profiles.** (a) One-dimensional velocity-depth profiles from the tomographic velocity model in Extended Data Fig. 1 and (b) from the FWI model in Fig. 2a, between the distances from 112 km to 192 km at an increment of 2 km. (c) and (d) show two velocity profiles from the tomographic model in Extended Data Fig. 1 and the FWI model in Fig. 2a, at distance 146 km and 182 km from the MAR, respectively. The thin dashed line indicates the interpreted Layer 2A/2B boundary from FWI; the bold dashed blue and black lines indicate the boundaries of Layer 2/3 from the tomographic model and the FWI model, respectively. The top low-velocity layer (indicated by letter 'H') could be related to hydrothermal alteration.
(a)

<table><thead><tr><th>Physical properties</th><th>Olivine (Ol)</th><th>Clinopyroxene (Cpx)</th><th>Plagioclase (Pl)</th></tr></thead><tbody><tr><td>Bulk moduli (GPa)</td><td>129.2</td><td>107.8</td><td>55.1</td></tr><tr><td>Shear moduli (GPa)</td><td>81.2</td><td>75.7</td><td>29.7</td></tr><tr><td>Density (g/cm³)</td><td>3.23</td><td>3.20</td><td>2.61</td></tr></tbody></table>

(b)

<table><thead><tr><th></th><th>Olivine<br>(Ol)</th><th>Clinopyroxene<br>(Cpx)</th><th>Plagioclase (Pl)</th><th>Vp (km/s)</th></tr></thead><tbody><tr><td>Olivine-rich gabbro: I</td><td>0.25</td><td>0.25</td><td>0.5</td><td>7.072</td></tr><tr><td>Olivine-rich gabbro: II</td><td>0.2</td><td>0.3</td><td>0.5</td><td>7.050</td></tr><tr><td>Olivine gabbro: I</td><td>0.15</td><td>0.3</td><td>0.55</td><td>6.934</td></tr><tr><td>Olivine gabbro: II</td><td>0.1</td><td>0.3</td><td>0.6</td><td>6.818</td></tr><tr><td>Olivine gabbro: III</td><td>0.1</td><td>0.4</td><td>0.5</td><td>7.005</td></tr><tr><td>Olivine-bearing<br>gabbro: I</td><td>0.05</td><td>0.3</td><td>0.65</td><td>6.703</td></tr><tr><td>Olivine-bearing<br>gabbro: II</td><td>0.05</td><td>0.25</td><td>0.7</td><td>6.611</td></tr></tbody></table>

**Extended Data Fig. 10 | Rock physics modelling using the Voigt-Reuss-Hill (VRH) averaging.** (a) Physical properties of Olivine (Ol), Clinopyroxene (Cpx) and Plagioclase (Pl) from [Sobolev and Babeyko 1994]. (b) The effective seismic P-wave velocities from Voigt-Reuss-Hill (VRH) averaging using different fraction volumes of mineral components (Ol, Cpx and Pl).