Skip to main content

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

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

Abstract

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−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.

Main

Oceanic crust is formed at the ocean spreading centres with different accretionary mechanisms1,2,3, 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 structures4. Additionally, hydrothermal circulations play important roles by releasing heat and altering rock properties5. Based on insights from ophiolite studies6,7, ocean drilling results8 and geophysical studies3,9,10,11,12, 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−1)9,11, whereas the lower crust contains mainly gabbro where the velocity increases very slowly (~0.1 s−1)3,9,10. 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 faults4. Although the upper crustal structure and accretion process have been well studied3,8,9,10,11,12,13, our knowledge about the lower crust remains poorly constrained14. The lower crust could represent two-thirds of the oceanic crust3, 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’ model1,15 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’ model6,7 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-16,17,18 and intermediate-spreading mid-ocean ridges19,20,21 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 axis22,23,24. 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 data25 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−1 (ref. 26). 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 tomography14 (Methods) is shown in Extended Data Fig. 1. The crustal thickness in the study area is 5.6 ± 0.2 km (ref. 14). 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 crust14.

Fig. 1: Map of the seismic survey area near the MAR in the equatorial Atlantic Ocean.
figure 1

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 age26. The bathymetry was taken from ref. 14. See the inset map for the location of the study area.

Starting from the travel-time tomographic model14 (Extended Data Fig. 1), we performed a full waveform inversion (FWI; Methods)27,28,29. 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-surface29. The FWI was performed in two steps: (1) a trace-normalized FWI inversion27 was performed to obtain an intermediate-scale velocity model; and (2) then a true-amplitude FWI28,29 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 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. 25). To make sure the inverted results are real, not artefacts of inversion, extensive modelling and inversion tests were performed (Methods and Supplementary Figs. 223), 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. 1420).

Fig. 2: Seismic P-wave velocity models of the oceanic crust.
figure 2

ac, FWI28 velocity model starting from the tomographic seismic velocity model14 (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−1 with an increment of 0.2 km s−1. 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−1 (ref. 26). The horizontal distance starts with 0 km at the MAR.

Fig. 3: Observed and synthetic seismic waveform data from OBS 57.
figure 3

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−1 velocity was applied to both the field and synthetic data.

When compared with the tomographic model, the velocity anomaly obtained using trace-normalized FWI27 (Methods and Extended Data Figs. 6 and 7) indicates an overall increase in the velocity (200–400 m s−1 positive anomaly) above 6–7 km depth and a decrease in velocity (100–300 m s−1 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 FWI28,29, 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−1 at the top to 6.2 km s−1 at its base; (2) an underlying medium velocity gradient layer of 600–800 m thickness, with the velocity increasing to ~6.6 km s−1; and (3) the alternate high- and low-velocity layers (400–500 m thick) with a velocity variation of ±200 m s−1 underlying.

Fig. 4: Modelled layering in the lower oceanic crust.
figure 4

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.

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) boundary12,13,30, where drilling results30,31 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 zone8,30, which corresponds to lithology changes including a sudden appearance of hydrothermal minerals, and an abrupt reduction of permeability and porosity30.

A recent tomographic study32 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 axis32. 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 sequence33, whereas a thick (>1,000 m) dyke sequence is observed in Hole 504B34; 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 alteration8,35,36; (2) alternate gabbro sills and serpentinized peridotite within the lower crust, possibly associated with oceanic core complexes (OCCs)4,37; and/or (3) layering introduced either by the gabbro glacier model1,15 (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)6,7.

The hydrothermal circulation is expected to reach down to the top of the AML beneath the ridge axis10,35. 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 seafloor38, 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 investigation35 has imaged a 150–200 m thick, low-velocity (0.3–0.5 km s−1 decrease) zone above the roof of the AML, which was interpreted to be formed from enhanced hydrothermal alteration35. Furthermore, a velocity decrease of 0.4–0.6 km s−1 in a 100–200 m thick layer has been observed in Hole 1256D at the dyke–gabbro transition36, with noticeable changes in minerology and decreasing porosity associated with hydrothermal alteration and contact metamorphism. A recent drilling result39 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 roof39,40. 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 AML35, consistent with the observation at the Lucky Strike38. Below this depth, the role of the hydrothermal circulation would be much reduced5.

Lower crustal accretion processes

Oceanic crust formed in a slow-spreading environment is suggested to be heterogeneous and may even contain mantle rocks2,4. Tectonic processes, such as long-lived detachment faults4, can produce OCCs, exhuming lower crustal and mantle rocks to the seafloor4,37. For example, FWI of multi-channel streamer data from the Kane OCC37 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 range41. However, the tomographic results14 in our area suggest that lower crustal velocity increases from ~6.6 km s−1 at 2 km below the basement to ~7 km s−1 just above the Moho, suggestive of a gabbroic origin. The mantle velocity just below the Moho is ~8 km s−1, confirming the presence of mantle peridotite14. 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 model1,15 or the sheeted sill model6,7 from the velocity models alone; however, the recent discoveries of stacked melt lenses in the lower crust at the fast-16,17,18 and intermediate-spreading ridges19,20,21 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 segments38,42, while our study region is supposed to lie in the cold mantle temperature regime14,43. 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)22,23,24. The rough seafloor 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 observed6,7. 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)6,7. In the Bay of Islands ophiolite, a representative of slow-spreading environment, layers with thickness of several hundred metres to a kilometre have been observed44.

A more than 1,500 m of the lower crustal gabbro was drilled during ODP Hole 735B45 from an ultra-slow-spreading environment containing several layered units of >250 m thickness45,46. Analysis of a geochemical unit in the lower gabbro suggests that it constitutes a single magma reservoir but is separated in two parts46. 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 accumulation46,47. 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.

Observations from Integrated Ocean Drilling Program (IODP)/ODP drilling and ophiolites39,45,46,47 indicate that the lower crustal gabbroic rocks are mainly composed of Ol, Cpx and Pl, with seismic velocities VOl > VCpx > VPl (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 velocities31,45. Using the Voigt–Reuss–Hill averaging method48 (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−1 (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 faults4, 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 Kane37 and Rainbow49,50 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 geophysical16,17,18,19,20,21,22,23,24 and ODP/IODP drilling45,46,47 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 cooling5, which can be provided through well-developed faults/fractures in the slow-spreading environment46. Moreover, the hydrothermal activities and the magmatic reaction with host rocks could lead to remelting and assimilation40,46,47, 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 similar data from different oceanic settings, to shed light on the formation and evolution of oceanic lithosphere.

Methods

Data acquisition

The field OBS data were acquired during the LITHOS-iLAB cruise in 2017 onboard the Germain RV Maria S. Merian to study the upper lithosphere from 0 to 50 Ma14,51. 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 transformation28,29 (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 Unix52, 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 offsets14 (Supplementary Fig. 1).

Travel-time tomography

Travel-time tomography14,53 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 times53; the model update was obtained by least-squares penalties on the data misfit, together with smoothing and damping for regularization53.

The velocity model from travel-time tomography (Extended Data Fig. 1) gives very good travel-time fit for the first arrivals14. 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 depth14. Details of the travel-time tomography can be found in ref. 14.

FWI

FWI is the current state-of-the-art technique for high-resolution sub-surface quantitative imaging27,28,29,51. Unlike tomography only using travel time53, 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 propagation28. 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 wavefields28,29. We used time-domain staggered-grid finite-difference54 with fourth-order spatial and second-order temporal accuracy for solving the elastic wave equation in the stress and particle-velocity formulation. The convolutional perfectly matched layer absorbing boundary condition55 was applied at the model boundaries.

FWI can provide high-resolution quantitative Earth internal models for physical properties and their vertical differentiations. Sensitivity studies56,57 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 crust14, FWI can update the lower crustal model with wide-angle reflections/diffractions from lower crustal anomalies being accurately modelled58. The application of FWI is mainly in the upper oceanic crust17,28,35,37,59, with very few exceptions down to 3–5 km depth60.

Although the tomographic inversion has converged with good travel-time fit14, 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 FWI27 was applied, with the misfit function (J) defined as:

$$J = \mathop {\sum }\limits_{i = 1}^{{\mathrm{Ns}}} \mathop {\sum }\limits_{j = 1}^{{\mathrm{Nr}}} s_{i,j}\left(\frac{{s_{i,j}}}{||{s_{i,j}}||} - \frac{{d_{i,j}}}{||{d_{i,j}}||}\right),$$
(1)

where s and d and 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 FWI28,29, with the misfit function as:

$$J = \mathop {\sum }\limits_{i = 1}^{{\mathrm{Ns}}} \mathop {\sum }\limits_{j = 1}^{{\mathrm{Nr}}} ||s_{i,j} - d_{i,j}||^2.$$
(2)

Source wavelet is important for seismic wave modelling in FWI and was estimated by stacking near-offset free-surface multiples of water waves51 (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 structures58. We didn’t apply FWI to the PmP arrivals, because of its strong nonlinearity around critical angles51. 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’s61 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−1, 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−1, 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. Studies23,64 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 km3. 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 velocity65. 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 layering66; 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−1 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−1 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−1 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. 1820). 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 waveform, 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)14 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 averaging48 is a common and simple approach for computing the effective elastic moduli of rocks from their mineral constituents. Evidence39,45,46,47 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 densities67 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−1 (±100–200 m s−1; 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 moduli31.

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 pickings 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

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

  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 2o 43’S. J. Geophys. Res. Solid Earth 126, e2020JB021390 (2021).

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

  22. Canales, J. P., Collins, J. A., Escartín, 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).

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

  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. Int. 224, 1056–1078 (2021).

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

Download references

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 information

Authors and Affiliations

Authors

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.

Corresponding author

Correspondence to Peng Guo.

Ethics declarations

Competing interests

The authors declare no competing interests.

Peer review

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.

Additional information

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

Extended data

Extended Data Fig. 1 The tomographic seismic P-wave velocity model from travel-time tomography.

It serves as the starting model for the FWI.

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.

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.

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.

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.

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/year26. 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.

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.

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.

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.

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

Supplementary information

Supplementary Information

Supplementary Figs. 1–24.

Rights and permissions

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/.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Guo, P., Singh, S.C., Vaddineni, V.A. et al. Lower oceanic crust formed by in situ melt crystallization revealed by seismic layering. Nat. Geosci. 15, 591–596 (2022). https://doi.org/10.1038/s41561-022-00963-w

Download citation

  • Received:

  • Accepted:

  • Published:

  • Issue Date:

  • DOI: https://doi.org/10.1038/s41561-022-00963-w

Search

Quick links

Nature Briefing

Sign up for the Nature Briefing newsletter — what matters in science, free to your inbox daily.

Get the most important science stories of the day, free in your inbox. Sign up for Nature Briefing