Lower oceanic crust formed by in situ melt crystallisation revealed by seismic layering.

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 dikes 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-million-year-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 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, that cool and crystallise 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[1-3], depending upon the spreading rates. At fast and intermediate spreading centres, the crustal accretion is mainly dominated by magmatic process, 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[4]. Additionally, hydrothermal circulations play important roles by releasing heat and altering rock properties[5]. Based on insights from ophiolite studies[6,7], ocean drilling results[8] and geophysical studies[3,9-12], the crust formed by magmatic process 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 dikes 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 faults[4]. Although the upper crustal structure and accretion process have been well studied[3,8-13], our knowledge about the lower crust remains poorly constrained[14]. The lower crust could represent two-thirds of the oceanic crust[3], 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[1,15] suggests that the lower crust is formed by cooling, crystallisation and subsequent subsidence of the gabbro from the axial melt lens (AML) at the upper-lower crust boundary, while the ‘sheeted sill’ model[6,7] suggests that the lower crustal gabbro formation could occur through in situ cooling and crystallisation of melt sills. Recent seismic discoveries of melt lenses in the lower crust beneath the fast[16-18] and intermediate spreading mid-ocean ridges [19-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 axis[22-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 data[25] were acquired aboard the German R/V Maria S. Marian in the equatorial Atlantic Ocean (Fig. 1) over crust formed at the slow spreading Mid-Atlantic Ridge with a half spreading rate of ~16 mm/year[26]. An air gun array with a total volume of 89.14 litre was fired at ~400 m interval, which was recorded on ocean bottom seismometers (OBSs, Fig. S1) deployed at a spacing of 10 to 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 to 12 Ma old seafloor on the African plate. The velocity model obtained using the travel-time tomography[14] (see Methods) is shown in Extended Data Fig. 1. The crustal thickness in the study area is 5.6 ± 0.2 km[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 crust[14].

Fig. 1

Map of the seismic survey area near the Mid Atlantic Ridge (MAR) in the equatorial Atlantic Ocean.

The dashed black line refers to the MAR, and the triangles show the locations of the ocean bottom seismometers (OBSs) used for this study. The bold black line passing through 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 Ma lithospheric age[26]. The bathymetry was taken from Vaddineni et al. (2021)[14]. See the inset map for the location of the study area.

Extended Data Fig. 1

The tomographic seismic velocity model[14] from travel-time tomography.

It serves as the starting model for the FWI.

Starting from the travel-time tomographic model[14] (Extended Data Fig. 1), we performed a full waveform inversion (FWI)[27-29] (FWI, see 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[29]. The FWI was performed in two steps: (1) a trace-normalised FWI inversion[27] was performed to obtain an intermediate-scale velocity model, and (2) then a true-amplitude FWI[28,29] was performed to obtain fine-scale velocity information (See Methods). The final 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, Extended Data Figs. 2-5). To make sure the inverted results are real, not artefacts of inversion, extensive modelling and inversion tests were performed (See Methods, Figs. S2-S23), 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 (See Methods, Figs. S14-S20).

Fig. 2

Seismic velocity models of the oceanic crust.

(a) FWI[28] velocity model starting from the tomographic seismic velocity model[14], (b) the velocity anomaly (the difference between the velocity models from FWI and tomography), and (c) the vertical velocity gradient (the derivative of velocity model from FWI with respect to depth). Black triangles in (a) mark the OBS locations. The velocity contours 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 lithospheric age in (c) is calculated using a spreading rate of 16 mm/year[26]. The horizontal distance starts with 0 km at the MAR.

Fig. 3

Observed and synthetic seismic waveform data from OBS57.

(a) The field data (black) and the synthetic waveform (red) using the tomographic model (Extended Data Fig. 1), and (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.

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. 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 Equation 2 in the Method section and then normalised.

When compared with the tomographic model, the velocity anomaly obtained using trace-normalised FWI[27] (see Methods, 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, Extended Data Fig. 8c) from the true-amplitude FWI[28,29], extending up to 5-15 km distance along the profile.

Extended Data Fig. 6

Seismic 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[26]. 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.

One-dimensional velocity profiles (Fig. 4a, 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.

Fig. 4

Modelled layering in the lower oceanic crust.

(a) One-dimensional velocity-depth profiles at the distance of 182 km from the MAR, derived from tomography (blue) and FWI (black), respectively. 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 9 km depth (~5.4 km from basement at distance 182 km), because of potential Moho reflection interference for greater depths (Fig. S22). (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 dikes. The blue ellipsoids refer to hydrothermal alteration above the roof of AML that could contribute to the top low-velocity layer, which is much reduced below AML. 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 be comprised of multiple thinner layers as indicated by the darker brown patches, but is beyond data resolution.

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.

The Geological Structure of Oceanic Crust

The upper layer with high velocity gradient is probably associated with the lava flow-dike (layer 2A/2B, Fig. 4a) boundary[12,13,30], where drilling results[30,31] suggest that the rapid reduction of porosity is the main reason for this fast velocity changes. A porosity decrease will cause an increase of velocity (Fig. S24). An abrupt transition from low-temperature hydrous alteration above to high-temperature hydrothermal alteration below has been observed at the lava-dike transition zone[8,30], which correspond to lithology changes including a sudden appearance of hydrothermal minerals, and an abrupt reduction of permeability and porosity[30].

A recent tomographic study[32] 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-dike boundary at ~ 900 m below the seafloor underneath the ridge axis[32]. 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 dike sequence. The drilling results from the Ocean Drilling Program (ODP) Hole 1256 D indicate a thin (350 m) dike sequence[33] whereas a thick (>1000 m) dike sequence is observed in Hole 504 B[34]; 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[8,35,36], (2) alternate gabbro sills and serpentinised peridotite within the lower crust, possibly associated with oceanic core complexes[4,37], and/or (3) layering introduced either by the gabbro glacier model[1,15] (crystallisation and subsidence of the gabbro from the shallow AML) or the sheeted sill model (in situ crystallisation 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 axis[10,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 seafloor[38], 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[35] has imaged a 150-200 m thick low-velocity (0.3 – 0.5 km/s decrease) zone above the roof of AML, which was interpreted to be formed from enhanced hydrothermal alteration[35]. Furthermore, a velocity decrease of 0.4 – 0.6 km/s in 100-200 m thick layer has been observed in Hole 1256 D at the dike/gabbro transition[36], with noticeable changes in minerology and decreasing porosity associated with hydrothermal alteration and contact metamorphism. A recent drilling result[39] from the Oman ophiolite of the dike/gabbro transition reveals that this transition is characterised by low-Mg and Zr-enriched gabbro and diorite, possibly related to the interaction between the magmatic and hydrothermal systems at the AML roof[39,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 AML[35], consistent with the observation at the Lucky Strike[38]. Below this depth, the role of the hydrothermal circulation would be much reduced[5].

Lower Crustal Accretion Processes

Oceanic crust formed in a slow spreading environment is suggested to be heterogeneous and may even contain mantle rocks[2,4]. Tectonic processes, such as long-lived detachment faults[4], can produce oceanic core complexes (OCCs), exhuming lower crustal and mantle rocks to the seafloor[4,37]. For example, FWI of multi-channel streamer data from the Kane OCC[37] of the MAR has revealed a ~300 m thick low-velocity layer at ~1.5 km depth from seafloor interpreted as serpentinised peridotite underlain by a thick gabbroic body, and a deeper ~200 m thick high-velocity layer interpreted as a gabbroic sill embedded in serpentinised peridotite. The P-wave velocity of serpentinised peridotite can have a wide range depending on the degree of serpentinisation, and the velocity in the layered structure may well be within this range[41]. However, the tomographic results[14] in our 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[14]. 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[1,15] or the sheeted sill model[6,7] 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[19-21] from seismic reflection studies support the lower crustal accretion by in situ melt intrusion and crystallisation. For slow-spreading ridges, melt lenses have been imaged at the upper/lower crust boundary for melt-rich spreading centre segments[38,42], while our study region is supposed to lie in the cold mantle temperature regime[14,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-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 (cm to 100 m) of alternating strata rich in mafic minerals (olivine and clinopyroxene) at the base and plagioclase at the top have been observed[6,7]. The geochemical analysis of these layers suggests the presence of secondary sills and that the lower crust could be formed by cooling and crystallisation 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 observed[44].

A more than 1500-m of the lower crustal gabbro was drilled during the ODP Hole 735 B[45] from an ultra-slow spreading environment containing several layered units of > 250 m thickness[45,46]. Analysis of a geochemical unit in the lower gabbro suggests that it constitutes of a single magma reservoir but is separated in two parts[46]. 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[46,47]. This process would lead to a layering in the lower crust, containing olivine-rich gabbro and troctolites at the base with a distinct boundary separating more evolved olivine and olivine-bearing gabbro with decreasing olivine and increasing plagioclase.

Observations from IODP/ODP drilling and ophiolites[39,45-47] indicate that the lower crustal gabbroic rocks are mainly composed of olivine (Ol), clinopyroxene (Cpx), and plagioclase (Pl), with seismic velocities VOl > VCpx > VPl (Extended Data Table 1a), where V represents the velocity. Rocks rich in Ol and Cpx, indicative of more primitive melt, would have higher velocities, whereas those of more evolved rocks rich in Pl would have relatively lower velocities[31,45]. Using the Voigt-Reuss-Hill averaging method[48] (See Methods), we computed the P-wave velocities for different gabbro compositions from the Hole 735 B and found that the velocity of the gabbro varies from 6.6 km/s to 7.1 km/s (Extended Data Table 1b), 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 layer may be rich in rocks with low Ol/Cpx (<10%) and high Pl contents, such as Ol-bearing gabbro.

Extended Data Table 1

Rock physics modelling using Voigt-Reuss-Hill (VRH) averaging.

(a) Physical properties of Olivine (Ol), Clinopyroxene (Cpx) and Plagioclase(Pl) from [67]. (b) The effective seismic P-wave velocities from Voigt-Reuss-Hill (VRH) averaging using different fraction volumes of mineral components (Ol, Cpx and Pl).

(a)
Physical properties	Olivine (01)	Clinopyroxene (Cpx)	Plagioclase (Pl)
Bulk moduli (GPa)	129.2	107.8	55.1
Shear moduli (GPa)	81.2	75.7	29.7
Density (g/cm3)	3.23	3.20	2.61

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 to 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[4], 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 (Figs. S15, S17), suggesting a more variable magma supply, but is beyond data resolution. Studies at the Kane[37] and Rainbow[49] 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 Crystallisation in the Lower Crust

Fig. 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-years geophysical[16-24] and ODP/IODP drilling[45-47] studies, the discovery of lower crust stacked-layering off axis shows that in our study area the lower oceanic crust is more likely formed by in situ melt sill intrusion from mantle and crystallisation 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 likely 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 crystallisation requires extensive seawater circulation down to the Moho depth along the sides of crystal-rich mush zone beneath the magma chambers for cooling[5], which can be provided through well-developed faults/fractures in the slow-spreading environment[46]. Moreover, the hydrothermal activities and the magmatic reaction with host rocks could lead to remelting and assimilation[40,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 likely formed by in situ cooling and crystallisation 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 full waveform inversion 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 upon the formation and evolution of oceanic lithosphere.

Methods

Data acquisition

The field OBS data was acquired during the LITHOS-iLAB cruise in 2017 onboard the Germain R/V Maria S. Merian to study the upper lithosphere from 0-50 Ma[14,51]. A total of 71 instruments consisting of 55 ocean bottom seismometers (OBSs) and 16 ocean bottom hydrophones (OBHs) were deployed along a 1100-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 was sampled at 250 Hz. The active seismic source used in the survey comprised of 6 G-gun clusters (12 guns) configured as two sub-arrays with a total volume of 89.14 litre, 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[28,29] (multiplying the amplitudes of field data by , where t is the two-way travel-time, and convolving with ) 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[52], 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) ± 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[14] (Fig. S1).

Travel time tomography

Travel time tomography[14,53] was first applied for estimating a P wave velocity model (Extended Data Fig. 1). Tomography was carried out through a linearised 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[53]; the model update was obtained by least-squares penalties on the data misfit, together with smoothing and damping for regularisation[53].

The velocity model from travel time tomography (Extended Data Fig. 1) gives very good travel-time fit for the first arrivals[14]. 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[14]. Details of the travel-time tomography can be found in Vaddineni et al (2021)[14].

Full waveform inversion (FWI)

FWI is the current state-of-the-art technique for high-resolution subsurface quantitative imaging[27-29,51]. Unlike tomography only using travel time[53], the FWI is based on minimising 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[28]. For the numerical implementation, a gradient-based linearised 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[28,29]. We used time-domain staggered-grid finite-difference[54] 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 (CPML) absorbing boundary condition[55] was applied at the model boundaries.

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

Although the tomographic inversion has converged with good travel-time fit[14], 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-normalised FWI[27] was applied, with the misfit function defined as: where s and d and are synthetic and field data, respectively, || || is the l-2 norm, Ns and Nr are the number of sources and receivers. 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 normalisation but the amplitude variation within each trace remains. An important point of the trace-normalised inversion is to determine the amplitude normalisation factors between the synthetic and real data. In the second stage, we performed a true-amplitude FWI[28,29], with the misfit function as

Source wavelet is important for seismic wave modelling in FWI and was estimated by stacking near-offset free-surface multiples of water waves[51] (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 data has the most linear behaviour with subsurface properties and the wide-aperture data is sensitive to both the upper- and lower-crustal structures[58]. We didn’t apply FWI to the PmP arrivals, because of its strong nonlinearity around critical angles[51]. A time-window of 0.5 s (Fig. S2) 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[61] regression fit. Density was linked to the P-wave velocity based on the empirical relation in Hamilton (1978)[62] for velocity smaller than 2.2 km/s, and the Gardner et al. (1974)’s relation[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 regularisation.

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 (Fig. S3), 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. Fig. S4 shows the seismograms are better aligned after trace-normalised 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).

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.

The velocity models from the trace-normalised FWI are shown in Extended Data Figs. 6 and 7, and those from the true-amplitude FWI are shown in Figs. 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-normalised FWI highlights the large-scale velocity variations (Extended Data Figs. 6b and 7b), which sharpens the upper/lower crust transition, whereas the vertical velocity gradient from true-amplitude FWI enhances the layered anomalies in the lower crust (Figs. 2c and Extended Data Fig. 8c). Models in Extended Data Figs. 7 and 8 have better horizontal continuity than those of 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-normalised FWI step (Fig. S5). The recovered models are generally consistent with those in Fig. 2, with relatively larger data misfit (Fig. S6). Considering the data frequency is 3.5-10 Hz, the dominant wavelength is 700-1600 m for a velocity of 6.8 km/s, indicating that the vertical resolution in the lower crust is ~ 400 m.

Extended Data Fig. 9 shows one-dimensional (1-D) 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 to the tomographic results. Extended Data Figs. 9c and 9d show a subset of two (1-D) velocity-depth profiles, with the interpreted boundaries of Layer 2A/2B and 2/3.

We didn’t consider velocity anisotropy. Studies[23,64] at MAR indicate velocity anisotropy 2-4% in the upper crust, and becomes nearly isotropic from 3 km below the seafloor. The observed anisotropy decreases fast off-axis over 5-10 km[3]. Therefore, we consider the influence of P-wave anisotropy 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[65]. Using surface-wave data, Russell et al (2019)[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[66]; however we can only invert layering of a quarter of wavelength using limited offset range. We consider intrinsic attenuation at the lithospheric age 7 -12 Ma weak, and thus not included for FWI. With significant attenuation and viscoelastic modelling, we expect the layered anomaly 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 (Fig. S7), by replacing the layering structures with the corresponding trace-normalised 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. Fig. S8 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 lower crust also appear at near offsets but at 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 (Figs. S9, S10) and quantified (Fig. S11a, S11b) especially for the far (larger than 15 km, Fig. S11b) offsets for models without layering. We used the modified model as the starting model to run FWI. The inverted model is shown in Fig. S12 with remarkable similarity to Fig. 2. In addition, we also ran FWI for the field data in which the velocity structure greater than depths 2.5 km from basement was fixed, so the inversion was forced to search for update in the upper crust to explain the data. The final models are shown in Fig. S13, with the main difference of the upper crust with Fig. 2 being a high velocity anomaly around 180 km distance. Waveform misfit is larger than updating the whole model, especially for large offsets (Fig. S11c, S11d). 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 (Fig. S14) 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 Fig. S15, we found that inversion results with 0.5 km and 1 km horizontal smoothing contain significant ‘smiling’ artifacts. Inversion discontinuities can also be observed with small horizontal smoothing and extensive anomalies. Conversely, large horizontal smoothing length may cause anomaly size 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 Fig. S16a and S16b 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 (Fig. S16c). With the vertical smoothing length being the wavelength, the layered anomalies cannot be reconstructed (Fig. S16d). 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. Fig. S17 shows that thin layers will be inverted as larger thickness (Fig. S17e and Fig. S17h). If there are multiple layers with separation smaller than FWI resolution, they may not be recovered properly (Fig. S17f and Fig. S17g).

Checkerboard 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, Fig. S18a; 10 km long and 0.5 km thick, Fig. S19a) to the tomographic velocity model. Another test with random anomalies of thickness 400 m but variable horizontal extension is shown in Fig. S20. Most of the anomalies are well resolved with good match to the true models (Fig. S18-S20). Different anomalies with distance < 2 km were imaged as one (Fig. S20b). 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 artifacts 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 (Fig. S21a)[14] for generating the ‘observed’ data. For FWI, the velocity model after removing the Moho and upper mantle (Fig. S21b) was used as the starting model. We observed weaker anomalies around the crustal base in the inverted model (Fig. S22), but their influence for the observed crustal layering is limited. These tests confirm that the layered structure in the lower crust, especially those away from the crustal base, represents real features of the sub-surface; the Moho reflection is unlikely attributed to the layered structures in the lower crust. We also test the influence of muting time window size for FWI. Fig S23 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 rocks 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 (VRH) averaging[48] is a common and simple approach for computing the effective elastic moduli of rocks from their mineral constituents. Evidence[39,45-47] from IODP/ODP drilling and ophiolite show that the gabbroic rocks of the lower crust are mainly composed of olivine (Ol), clinopyroxene (Cpx), and plagioclase (Pl). Therefore, we calculated the P-wave velocities using different combinations of Ol, Cpx and Pl volume fractions. The elastic moduli and densities[67] of the three mineral components are shown in Extended Data Table 1a. 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 Table 1b), 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[31].