# Water-rich sublithospheric melt channel in the equatorial Atlantic Ocean

Fares Mehouachi<fnref n="1" /> and Satish C. Singh<fnref n="1" />,²*

The lithosphere-asthenosphere boundary is the most extensive boundary on Earth, separating the mobile plate above from the convecting mantle below, but its nature remains a matter of debate. Using an ultra-deep seismic reflection technique, here we show a systematic seismic image of two deep reflectors that we interpret as the upper and lower limits of the lithosphere-asthenosphere boundary beneath a 40-70-million-year-old oceanic lithosphere in the Atlantic Ocean. These two reflections correspond to 1,260 °C and 1,355 °C isotherms and bound a low-velocity channel, suggesting that the lithosphere-asthenosphere boundary is thermally controlled. We observe a clear age dependency of this sublithospheric channel: its depth increases with age from 72 km where it is 40-Myr-old to 88 km where it is 70-Myr-old, whereas its thickness decreases with age from 18 km to 12 km. We suggest that partial melting, facilitated by water, is the main mechanism responsible for the low-velocity channel. The required water concentration for melting increases with age; nevertheless, its corresponding total mass remains relatively constant, suggesting that most of the volatiles in the oceanic sublithospheric channel originate from a horizontal flux near the ridge axis.

The plate tectonic theory requires the presence of a rigid lithosphere floating over a ductile asthenosphere, separated by a lithosphere-asthenosphere boundary (LAB). The LAB corresponds to a change in the rheological properties of the upper mantle, which is generally associated with a seismic low-velocity zone (LVZ) whose origin has been extensively debated and is still a matter of controversy<fnref n="1" />. The smooth shear velocity drop imaged by surface wave tomographic studies²,³ is inconsistent with the sharp S-wave contrast reported for receiver function observations⁴⁻⁷. A recent seismic reflection study in a subduction setting⁸ indicated the presence of two P-wave reflections 10 km apart at a depth 100 km for a 120-Myr-old subducting oceanic lithosphere, further enhancing the debate. Partial melting and grain boundary sliding (GBS) are the two mechanisms that could potentially explain the seismic velocity variation at the base of the lithosphere⁹⁻¹². The GBS model predicts a single age-independent boundary for the LAB⁶. On the other hand, melt existence is not guaranteed since the temperatures at the imaged LAB depths are too low to sustain melt. To better understand the nature of the LAB, we acquired seismic reflection data over normal oceanic lithosphere covering an extended lithospheric age range with a fine sampling usually used for oil and gas exploration.

## Ultradeep seismic reflection survey

Our study area lies in the equatorial Atlantic Ocean across the St Paul Fracture Zone (FZ) over 40 to 70-Myr-old oceanic lithosphere (Fig. 1). The equatorial Atlantic Ocean opened during the Late Jurassic epoch, breaking the West Gondwanaland and separating northern South America from central Africa initially in a trans- extensional regime creating large fracture zones (Chain, Romanche and St Paul)¹³. The initial half-spreading rate was 27 mm yr⁻<fnref n="1" /> 110 million years ago (Ma), decreasing to 16 mm yr⁻<fnref n="1" /> around 50 Ma, increasing to 20 mm yr⁻<fnref n="1" /> 40 Ma and then remaining about 18 mm yr⁻<fnref n="1" /> over the past 20 Myr (ref. ¹⁴).

We acquired a 215-km-long ultradeep seismic reflection profile across the St Paul FZ. Although the total offset across the St Paul FZ

is ~580 km, corresponding to an age contrast of 30 Myr, the detailed analysis of the bathymetry and magnetic anomaly data suggests that the fracture zone is about 60-km wide containing small spreading segments with discrete age variations (Fig. 1). Furthermore, our seismic profile is south-north and crosses the fracture zone at 86°, making the age variation more gradual along the profile. The main advantage of acquiring data across the St Paul FZ is that we can cover an age variation of 30 Myr with just a 215-km-long profile. Moreover, this area has been affected neither by subduction or plume and hence provides information about normal oceanic lithosphere. Moreover, the presence of thin sediments over the older lithosphere reduces the scattering effect from rough oceanic seafloor, particularly for the seafloor formed at slow spreading centres.

The seismic reflection data were acquired using a 12-km-long multi-sensor streamer towed at a depth of 30 m (see Methods). The source was an array of air guns with a total volume of 1,671 deployed at 15 m water depth to generate low-frequency energy for deep penetration. The shot interval was 75 m and the record length was 30.7 s to observe reflected energy down to 100-km depth. Only frequencies between 1.5 Hz and 6 Hz were used for deep seismic imaging. The final seismic image in the time domain is shown in Fig. 2.

## Deep mantle seismic reflections across the St Paul FZ

The seismic image shows the presence of a veneer of sediments covering the oceanic crust along the profile. The main branch of the St Paul FZ crosses the profile between the 100- and 120-km distance range and is characterized by a basement low, filled with sediments. The Moho is intermittently imaged at 10 s two-way travel time (Fig. 2). On the southern end of the profile, we observe a reflection at 23.5 s that can be followed all along the profile deepening to 27 s at the northern end. This reflection could be a residual multiple, a sideswipe or a real reflection in the mantle. The overlapping plot of the first three multiples (Supplementary Fig. 1) indicates that this reflection lies between the second and third multiples from the 0- to 120-km distance range, south of the St Paul FZ and below the third

<footnote n="1">¹Equipe de Géosciences Marines, Institut de Physique du Globe de Paris (CNRS, Paris Diderot, Sorbonne Paris Cité), Paris, France. ²Earth Observatory of Singapore, Nanyang Technological University, Singapore, Singapore. *e-mail: singh@ipgp.fr</footnote>
<figure><img src="image_1.png" /><figcaption>**Fig. 1** | Study area of the seismic survey. White thick line: seismic profile [HI]. Red line: Mid-Atlantic Ridge (MAR). The numbers indicate the age of the oceanic lithosphere in Myr. Yellow lines: the major fracture zones in the equatorial Atlantic. The bathymetry around the seismic line is flat with a water depth of around 5 km, which is favourable for deep-target imaging.</figcaption></figure>

multiple, north of the St Paul FZ, and is not related to multiples. Furthermore, the deep reflection dips towards north from 23.5 s to 27 s, whereas the seafloor reflection gently dips towards south. The bathymetry data indicate that there are no features either east or west of the profile that could produce such a continuous 200-km-long reflection between 23 s and 27 s, and therefore, we argue that these reflections are real within the oceanic mantle. We also observe a second reflection at 27.5 s at the southern end of the profile that can also be followed along the whole profile more or less parallel to the first reflection deepening towards the north. Although it is coincident with the third multiple from the 0- to 35-km distance range, it is very distinct from the multiple afterwards. Further north, the second reflection could be followed until the end of the profile at 30 s, but the reflectivity is not that continuous. The bright reflectivity between the 35- and 120-km distance range suggests that this reflection is also real, not an artefact. The time difference between the two reflections decreases from ~4 s in the south of the St Paul FZ to ~2.5 s in the north.

## Low-velocity sublithospheric channel

The polarities of these reflections provide insight into the velocity contrast across the interface. In the zone where these reflections are bright, we find that the upper reflection has a negative polarity as compared with the seafloor reflection whereas the lower reflection has a positive polarity (Supplementary Fig. 1) indicating that the upper reflection is likely to be due to a decrease in velocity contrast whereas the lower reflection is due to an increase in velocity contrast. We found that the upper reflection requires a P-wave velocity decrease of 8.5 ± 4% and the lower reflection a velocity increase of 8.5 ± 4% (see the Methods) (Supplementary Fig. 2).

We depth-converted the time-domain image using one-dimensional simple seismic velocity (Supplementary Fig. 3) and the depth-converted seismic image is shown in Supplementary Fig. 4. The velocity values in the one-dimensional model are based on typical values for oceanic crust from ref. ¹⁵ and for oceanic mantle from LITHO1.0 (refs ⁸,¹⁶). We find that the upper reflection lies at 72 ± 2 km depth for a 40-Myr-old oceanic lithosphere increasing to 88 ± 2 km for a 70-million-year-old oceanic lithosphere. The second reflection lies at a depth of 90 ± 3 km where it is 40 Myr old and at a depth of 100 ± 3 km where it is 70 Myr old (Fig. 3 and Supplementary Fig. 5). The depth separation between the two reflections varies from 18 ± 0.8 km at 40 Myr old to 12 ± 0.5 km at 70 Myr old. We computed the isotherm for the slow-spreading

<figure><img src="image_2.png" /><figcaption>**Fig. 2** | Time-domain seismic image. The seismic image along profile [HI]. Green arrows: Moho reflections; yellow arrows: top LAB reflections; red arrows: bottom LAB reflections. The St Paul FZ (white arrow and dashed line) lies between 100- and 120-km distance, separating 70-Myr-old oceanic lithosphere in the north from 40-Myr-old lithosphere in the south. The base of the LAB is difficult to pick beyond the 180-km distance range, as the seismic image reaches its recording limit.</figcaption></figure>

Mid-Atlantic Ridge using a cooling plate model¹⁷ with a plate thickness of 106 km (ref. ¹⁸) and found that the upper reflection corresponds to a 1,260 ± 20°C isotherm, whereas the lower reflection corresponds to a 1,355 ± 15°C isotherm (Fig. 3b and Supplementary Fig. 5). Here, we have assumed a mantle melting temperature of 1,410°C (ref. ¹⁹) and an average spreading rate of 20 mm yr⁻¹.

Assuming a dominant frequency of 4–5 Hz and an average P-wave velocity of 8–8.5 km s⁻¹, the wavelength will be 1,600–2,000 m, and the vertical resolution (a measure of how thick an object needs to be in order to be seen in seismic) will be 400–500 m. These two reflections could be due to thin low- and high-velocity layers, or a sharp decrease and then increase in velocity from the two interfaces. An 8.5% decrease in the velocity for the upper reflection could be explained by a thin low-velocity layer at the base of the lithosphere, but there is no explanation for a thin high-velocity layer for the second reflection, and therefore, we suggest that these pairs of reflections represent the top and bottom of a LVZ.

Although our super common mid-point spacing is 50 m, the Fresnel zone at 70–100 km depth will be 5–8 km, and hence we can resolve 5–10-km-wide structures laterally. The blow-up of the seismic image (Fig. 3b and Supplementary Fig. 5) shows that the upper reflection suddenly deepens by 2 km at the 40-km and 80-km distance ranges, possibly related to the age contrast created by small ridge segments south of the St Paul transform fault (Fig. 3a). Across the St Paul FZ, the depth of the upper reflection increases from 77 km to 88 km linearly in a 30-km-wide zone, suggesting a smooth variation of temperature across the fracture zone. A similar depth variation is observed for the lower reflection, but it is less pronounced.

Both the depths and the isotherms suggest that these pairs of reflections represent the LAB. A pair of reflections, 10 km apart, were also observed for the subducting Pacific plate lithosphere of 120 Ma beneath New Zealand at 100 km depth⁸. The authors of ref. ²⁰ have observed a 25-km-thick conductive layer for the subducting Nazca Plate at 45 km depth for a 23-Myr-old lithosphere, which they interpret as the LAB. Although the thermal regime might get modified near the subduction zones, these observations are coherent with our new results. This constitutes strong evidence for the age dependency of the LAB both in depth and thickness.

## Nature of the LAB

The correspondence between the isotherms and the reflections suggests that the depth of the LAB and the thickness of the channel
<figure><img src="image_3.png" /><figcaption>**Fig. 3** | The LAB image. **a**, Age variation along the profile. **b**, Blow-up of the seismic image between the 65 and 100 km depth range. The white dashed line highlights the top and bottom reflectors and the red lines indicate the estimated positions of 1,260 °C and 1,355 °C isotherms. Inset: bathymetry image showing small segments along the St Paul FZ, affecting the depth of the LAB. The red arrows indicate small spreading segments on the MAR (inset) and associated age variations (main panel).</figcaption></figure>

are, to some extent, thermally controlled. Melt-free GBS can explain the smooth S-wave velocity drop imaged by surface wave tomography<sup>21</sup>. In an effort to reconcile the GBS model with receiver function observations, the presence of a steep water gradient at 70-km depth producing a sharp negative S-wave velocity contrast has been suggested<sup>9</sup>. The high temperature around the base of the lithosphere and the possible existence of water are likely to trigger GBS and weaken the shear resistance of rocks. However, GBS does not impact the bulk modulus<sup>22</sup> and hence the P-wave velocity is insensitive to its effects. As a consequence, a GBS-based LVZ would be transparent on a P-wave seismic section, suggesting that GBS is not sufficient to explain the velocity drop.

The second possible mechanism to reduce the P-wave velocity is the presence of melt<sup>10</sup>: the upper reflection could be interpreted as a freezing front separating partially molten material in the channel below from the impermeable lithosphere above. On the basis of the effective medium theory<sup>23</sup>, we estimated the melt concentration in the channel. To decrease the P-wave velocity by 8.5±4%, one would require about 1.4±0.8% of melt in the form of thin films with a thickness to length ratio of the order of 1:100 (Supplementary Fig. 6). Similar values were obtained in ref.<sup>8</sup>. The thin films of melt could be produced by localized shearing due to plate motion above and mantle convection below<sup>24</sup>. The stress-driven segregation of melt into thin films is crucial to produce the required velocity drop with a reasonable amount of melting (Supplementary Fig. 6). Above the channel, water concentration would be lower and the amount of frozen gabroic melt would be small and would not affect the velocity much. Below the channel, either the melt fraction diminishes abruptly or the melt inclusion shape changes from a thin film to a spherical form (Supplementary Fig. 6). Both mechanisms could coexist and explain the P-wave velocity increase at the base of the channel.

If we assume that the LAB is locally defined by the two isotherms, then the deepening of the LAB with age is consistent with the plate cooling model<sup>18</sup>, but brings out the issue of the melt stability. Water is considered as the main volatile in Earth that facilitates the melting of the mantle<sup>25</sup>. We estimated the water fraction needed to stabilize a 1.4% of melt throughout the channel<sup>26</sup> (see the Methods) and found that the required water concentration increases with age (Supplementary Fig. 7), varying from 945±345 ppm where it is 40-Myr-old to 1,515±380 ppm where it is 70-Myr-old for the upper freezing front. The water requirement for the lower freezing front is lower and varies from 610±285 ppm to 890±310 ppm, respectively.

Despite these apparent age variations in concentration, the total mass of water (normalized per surface unit section) in the channel is relatively constant: 13,960±5,700 ppm km where it is 40-Myr-old and 14,440±4,100 ppm km where it is 70-Myr-old. The average value varies by less than 3.4% over an age span of 30 Myr (~1% per 10 Myr) and is statistically meaningless given the relatively high uncertainties (Supplementary Fig. 8b). The asthenosphere is richer in water than the dry lithosphere<sup>9,27</sup> and water could potentially percolate from the asthenosphere to enrich the LAB. However, the small variation of the average total water content from 40 to 70 Myr makes such a mechanism improbable and suggests that most of the LAB's volatiles have been trapped in the channel at a younger age: the LAB could therefore be thicker near the ridge, with a lower water concentration that is focused and transported within the ever thinning LAB (Fig. 4 and Supplementary Fig. 8a).

Corner flow beneath the ridge axis was evoked as a possible source of volatiles in the LVZ<sup>10</sup>. However, its mechanism is not quantified. In addition, most of the existing melt transport models assume a vertical flow of melt that focuses at the ridge axis<sup>28,29</sup>.
<figure><img src="image_4.png" /></figure>

**Fig. 4** | The LAB as a function of age. A schematic diagram showing the evolution of LAB as a function of age between 40 and 70 Myr old. The white horizontal ellipsoids indicate the sheared melt bands responsible for the velocity drop in the LAB. Also marked are the depths and thickness variations of melt channels from refs. ⁸,²⁰. The correspondence of the isotherms with the LAB channel is local. The majority of volatiles in the LAB are trapped and focused near the ridge³¹.

Although it was believed that lateral flux is unlikely³⁰, recent modeling studies³¹,³² of melt production and focusing near the ridge support the existence of such as horizontal flux, indicating that 60% of deep volatile-rich melt could possibly be channelled along the LAB and therefore create the required feeding mechanism. If we extrapolate the constant total water hypothesis to the ridge and divide it by a 160 km base of the melting regime³⁰, then 90 ± 35 ppm of water fraction should be channelled in the LAB. This value is coherent with the 60% estimation from ref. ³¹ of the typical water content found under the ridge (50–200 ppm)³³–³⁵.

The thinning of the LAB could be due to freezing from the top and mechanical shearing. It has been suggested³⁶ that ductile mantle rock softening by water could lead to strain localization and grain size reduction that would weaken the rock fabric and create a micro-scale (grain size) pumping mechanism. The fluids (water/melt) would preferentially accumulate in the weakened layer of the LAB. This water accumulation would in turn enhance strain localization and therefore create a feedback mechanism that would lead to the thinning of the sheared fluid-rich layer.

Within this paradigm of the thinning water-rich melt channel, defining the LAB with two isotherms would be valid between 40 and 70 Myr old but would not hold for a younger lithosphere, since the isotherms converge near the ridge. In addition, the nature of the lower reflection might depend primarily on the melt focusing process along the LAB, and its seismic observation for younger lithospheric age might not be guaranteed given the potential interference of the established vertical melt focusing under the ridge.

The effect of melt on the S-wave velocity is more important (Supplementary Fig. 6b), which combined with the smearing effect and the large wavelengths of surface waves, might explain the large depth extent of the low S-wave velocity observed by the tomographic studies². The melt present in the LAB in the form of thin films will also produce radial anisotropy³⁷, and can explain the observed anisotropy³ below the oceanic lithosphere.

It has been suggested²⁵ that silicate and carbonatite melts would be stable at different depths. The continuity of the LAB across the St Paul FZ suggests that these two types of melt could potentially coexist to create a coherent boundary. Carbonatite melt has a conductance that is orders of magnitude higher than that of silicate

melt and this might explain the high conductivity observed by the authors of ref. ²⁰ at shallow depth. At greater depths, only silicate melt is stable in the oceanic LAB.

The small amount of melt in the LAB would reduce the viscosity significantly (three order of magnitude), and would decouple the floating lithosphere above with the convecting mantle underneath³⁸. The presence of molten LAB will have a significant impact on the plate tectonic process, requiring the plates to be driven mainly by slab pull and ridge push, not by a drag due to the convection in the mantle. The presence of the melt channel at the base of the lithosphere might also explain the long-distance transport of melt along leaky-transform faults³⁹ or from plume to ridges, such as the Rodrigues Ridge⁴⁰ or petit-spot¹¹.

## Methods

Methods, including statements of data availability and any associated accession codes and references, are available at https://doi.org/10.1038/s41561-017-0034-z.

Received: 2 August 2017; Accepted: 21 November 2017; Published online: 18 December 2017

## References

1. Eaton, D. W. et al. The elusive lithosphere-asthenosphere boundary (LAB) beneath cratons. *Lithos* **109**, 1–22 (2009).

2. Priestley, K. & McKenzie, D. The thermal structure of the lithosphere from shear wave velocities. *Earth Planet. Sci. Lett.* **244**, 285–301 (2006).

3. Burgos, G. et al. Oceanic lithosphere-asthenosphere boundary from surface wave dispersion data. *J. Geophys. Res. Solid Earth* **119**, 1079–1093 (2014).

4. Kawakatsu, H. et al. Seismic evidence for sharp lithosphere-asthenosphere boundaries of oceanic plates. *Science* **324**, 499–502 (2009).

5. Rychert, C. A. & Shearer, P. M. Imaging the lithosphere-asthenosphere boundary beneath the Pacific using SS waveform modeling. *J. Geophys. Res. Solid Earth* **116**, Q0AK10 (2011).

6. Olugboji, T. M., Park, J., Karato, S. I. & Shinohara, M. Nature of the seismic lithosphere-asthenosphere boundary within normal oceanic mantle from high-resolution receiver functions. *Geochem. Geophys. Geosyst.* **17**, 1265–1282 (2016).

7. Tharimena, S., Rychert, C. A. & Harmon, N. Seismic imaging of a mid-lithospheric discontinuity beneath Ontong Java Plateau. *Earth Planet. Sci. Lett.* **450**, 62–70 (2016).
8. Stern, T. A. et al. A seismic reflection image for the base of a tectonic plate. *Nature* **518**, 85 (2015).

9. Karato, S. I. On the origin of the asthenosphere. *Earth Planet. Sci. Lett.* **321**, 95–103 (2012).

10. Hirschmann, M. M. Partial melt in the oceanic low velocity zone. *Phys. Earth Planet. Inter.* **179**, 60–71 (2010).

11. Machida, S., Kogiso, T. & Hirano, N. Petit-spot as definitive evidence for partial melting in the asthenosphere caused by CO₂. *Nat. Commun.* **8**, 14302 (2017).

12. Tharimena, S., Rychert, C. & Harmon, N. A unified continental thickness from seismology and diamonds suggests a melt-defined plate. *Science* **357**, 580–583 (2017).

13. Darros De Matos, R. M. *Atlantic Riffs and Continental Margins* (American Geophysical Union, Washington, DC, 2000); 331–354.

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

15. White, R. S., McKenzie, D. & O'Nions, R. K. Oceanic crustal thickness from seismic measurements and rare earth element inversions. *J. Geophys. Res. Solid Earth* **97**, 19683–19715 (1992).

16. Pasyanos, M., Masters, T., Laske, G. & Ma, Z. LITHO1.0: An updated crust and lithospheric model of the Earth. *J. Geophys. Res. Solid Earth* **119**, 2153–2173 (2014).

17. Stein, C. A. & Stein, S. A model for the global variation in oceanic depth and heat flow with lithospheric age. *Nature* **359**, 123–129 (1992).

18. McKenzie, D., Jackson, J. & Priestley, K. Thermal structure of oceanic and continental lithosphere. *Earth Planet. Sci. Lett.* **233**, 337–349 (2005).

19. Sarafian, E., Gaetani, G. A., Hauri, E. H. & Sarafian, A. R. Experimental constraints on the damp peridotite solidus and oceanic mantle potential temperature. *Science* **355**, 942–945 (2017).

20. Naif, S., Key, K., Constable, S. & Evans, R. L. Melt-rich channel observed at the lithosphere–asthenosphere boundary. *Nature* **495**, 356 (2013).

21. Faul, U. H. & Jackson, I. The seismological signature of temperature and grain size variations in the upper mantle. *Earth Planet. Sci. Lett.* **234**, 119–134 (2005).

22. Ghahremani, F. Effect of grain boundary sliding on anelasticity of polycrystals. *Int. J. Solids Struct.* **16**, 825–845 (1980).

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

24. Holtzman, B. K. & Kohlstedt, D. L. Stress-driven melt segregation and strain partitioning in partially molten rocks: Effects of stress and strain. *J. Petrol.* **48**, 2379–2406 (2007).

25. Hirschmann, M. M., Tenner, T., Aubaud, C. & Withers, A. C. Dehydration melting of nominally anhydrous mantle: The primacy of partitioning. *Phys. Earth Planet. Int.* **176**, 54–68 (2009).

26. Katz, R. F., Spiegelman, M. & Langmuir, C. H. A new parameterization of hydrous mantle melting. *Geochem. Geophys. Geosyst.* **4**, 1073 (2003).

27. Masuti, S., Barbot, S. D., Karato, S. I., Feng, L. & Banerjee, P. Upper-mantle water stratification inferred from observations of the 2012 Indian Ocean earthquake. *Nature* **538**, 373–377 (2016).

28. McKenzie, D. The generation and compaction of partially molten rock. *J. Petrol.* **25**, 713–765 (1984).

29. Asimow, P. D. & Langmuir, C. H. The importance of water to oceanic mantle melting regimes. *Nature* **421**, 815 (2003).

30. Keller, T., Katz, R. F. & Hirschmann, M. M. Volatiles beneath mid-ocean ridges: Deep melting, channelised transport, focusing, and metasomatism. *Earth Planet. Sci. Lett.* **464**, 55–68 (2017).

31. Keller, T. & Katz, R. F. The role of volatiles in reactive melt transport in the asthenosphere. *J. Petrology* **57**, 1073–1108 (2016).

32. Plank, T. & Langmuir, C. H. Effects of the melting regime on the composition of the oceanic crust. *J. Geophys. Res. Solid Earth* **97**, 19749–19770 (1992).

33. Michael, P. Regionally distinctive sources of depleted MORB: Evidence from trace elements and H₂O. *Earth Planet. Sci. Lett.* **131**, 301–320 (1995).

34. Danyushevsky, L. V., Eggins, S. M., Falloon, T. J. & Christie, D. M. H₂O abundance in depleted to moderately enriched mid-ocean ridge magmas; part I: incompatible behaviour, implications for mantle storage, and origin of regional variations. *J. Petrol.* **41**, 1329–1364 (2000).

35. Sobolev, A. V. & Chaussidon, M. H₂O concentrations in primary melts from supra-subduction zones and mid-ocean ridges: implications for H₂O storage and recycling in the mantle. *Earth Planet. Sci. Lett.* **137**, 45–55 (1996).

36. Précigout, J., Prigent, C., Palasse, L. & Pochon, A. Water pumping in mantle shear zones. *Nat. Commun.* **8**, 15736 (2017).

37. Schmerr, N. The Gutenberg discontinuity: Melt at the lithosphere-asthenosphere boundary. *Science* **335**, 1480–1483 (2012).

38. Takei, Y. & Holtzman, B. K. Viscous constitutive relations of solid-liquid composites in terms of grain boundary contiguity: 1. Grain boundary diffusion control model. *J. Geophys. Res. Solid Earth* **114**, B06205 (2009).

39. Murton, B. J. Anomalous oceanic lithosphere formed in a leaky transform fault: evidence from the Western Limassol Forest Complex, Cyprus. *J. Geol. Soc.* **143**, 845–854 (1986).

40. Machida, S. et al. Regional mantle heterogeneity regulates melt production along the Réunion hotspot-influenced Central Indian Ridge. *Geochem. J.* **48**, 433–449 (2014).

## Acknowledgements

The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007–2013)/ERC Advance Grant agreement no. 339442\_TransAtlanticILAB. The discussions with C. Langmuir on petrology were extremely useful. This manuscript was substantially improved in response to the thoughtful suggestions of B. Romanowicz and C. Langmuir. T. Stern provided a very constructive review. This is an Intitut de Physique du Globe de Paris contribution number 3899.

## Author contributions

S.C.S. designed the project, led the data acquisition, supervised F.M. and wrote the paper. F.M. processed the seismic data, carried out all the calculations, produced the figures and participated in the writing of the paper.

## Competing interests

The authors declare no competing financial interests.

## Additional information

Supplementary information is available for this paper at https://doi.org/10.1038/s41561-017-0034-z.

Reprints and permissions information is available at www.nature.com/reprints.

Correspondence and requests for materials should be addressed to S.C.S.

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

**Data acquisition.** The WesternGeco marine seismic vessel Western Trident acquired the 215-km-long offset deep seismic reflection profile across the equatorial Atlantic Ocean in March-April 2015. The profile was shot in the south-southeast-north-northwest direction, nearly parallel to the Mid-Atlantic Ridge at an angle of about 4° and crossing the St Paul FZ. To enhance the low-frequency signals, the seismic data were acquired using a 12-km-long multi-sensor Isométrix streamer¹¹, towed at 30 m depth. The streamer contained two hydrophones and three-component accelerometers (a<sub>u</sub>, a<sub>p</sub>, a<sub>z</sub>) spaced at 0.75 m. The pressure data were pre-processed in real time to produce a receiver spacing of 3.125 m. The nearest receiver was at 135 m from the source. The seismic source was composed of 6 sub-arrays with 8 air-guns each, providing a total volume of compressed air of 10,170 cubic inch (1671) at an operating pressure of 2,000 ± 100 psi. The seismic source was towed at 15 m and the shot interval was 75 m, allowing for a recording length of 30,720 ms. Special effort during the acquisition was made to avoid barnacles, bend and connector noises.

**Data processing.** The initial processing was carried out on board. Anomalous amplitude attenuations were corrected, swell and linear noises were removed and random noise was attenuated. The source variation effects from shot to shot were corrected. Deghosting was carried out using both the pressure and the vertical component to provide an optimal reconstruction of the missing frequencies notches and enhancing the bandwidth, particularly in the low frequencies range. Multiples were attenuated using three proprietary techniques from Schlumberger: two-dimensional generalized surface multiple prediction, weighted least-squares demultiplex and post-radon isolating multiple algorithm demultiplex. These data were resampled to 6.25 m receiver interval and 8 ms time interval. The second phase of the processing was carried at the Institut de Physique du Globe de Paris. Super CMP (common mid-point) gathers were used for processing instead of normal CMP gathers to enhance the signal-to-noise ratio. Eight adjacent CMPs were combined to form a super CMP. Fourier filtering was carried out to keep only frequencies in the 1.5 to 6 Hz band. A smoothing two-dimensional median filter was also applied to enhance the lateral coherency and remove any spikes. Constant-velocity stack analysis was carried out to build a two-dimensional velocity model that yielded the clearest image of a time window between 20 and 30 s. Given that our primary imaging target is located between 70 to 100 km depths and the streamer was only 12 km long, no significant move-out should exist for real late arrival events. As such, non-flat events were filtered from the normal-move-out-corrected super CMP gathers before applying a near and far offset mute and stacking. A time linear gain was applied to enhance the deeper events. The post stack image was not migrated but instead was dip filtered to attenuate the tails of diffractions. A windowed frequency-wave number (FK) spectrum powering and a soft threshold were used to further enhance the seismic image iteratively until a satisfactory result was obtained. Empirical mode decomposition filtering in the frequency-distance (FX) domain was applied to further enhance the quality of the image and attenuate noise introduced by the previous steps¹². The final seismic image is shown in Fig. 2 and Supplementary Fig. 1.

**Depth conversion.** In the absence of a reliable two-dimensional P-wave velocity model, we decided to carry out the depth conversion using a simple one-dimensional model where the crust thickness is 6 km and the velocity is 6 km s⁻¹ (ref. ¹⁵), the water velocity is 1.5 km s⁻¹ and the mantle velocity varies linearly from 8 to 8.5 km s⁻¹ (ref. ¹⁶). As the mantle velocity could be between 8 and 8.5 km s⁻¹, we choose a linear increasing velocity averaging 8.25 km s⁻¹ velocity as a good compromise. As the velocity decrease in the LAB is about 8.5%, we used a P-wave velocity of 7.77 km s⁻¹. To avoid stretching artefacts at the boundaries, we used a smooth velocity variation in 0.5 s above and below the LAB. The depth conversion in this case consists of stretching and compressing the time axis along the corresponding velocity. As the velocity in the LAB velocity should be 8.5% less than that above (7.77 km s⁻¹), the uncertainty in the thickness of the LAB would be ±0.8 km where it is 40 Myr old and ±0.6 km where it is 70 Myr old (the resulting velocity is shown in Supplementary Fig. 3). Depth uncertainties are computed by changing the mantle velocity from a linear gradient to a constant value of 8 km s⁻¹ or 8.5 km s⁻¹. The LAB thickness depends on the velocity within the channel, and therefore the LAB velocity was reduced by 8.5 ± 4.3% from a mantle velocity of 8.5 km s⁻¹. We choose this maximal value to avoid unusually low velocities in the LAB⁹.

**Reflection coefficient estimation.** The reflection coefficient estimation was carried out by comparing the seismic amplitude of the target with a reference reflection whose impedance contrast is known¹⁴. In the absence of a strong continuous Moho, the seafloor was the only viable reflector we can use as a reference.

We choose to carry the reflectivity estimation in the stacked domain rather than the shot domain to tackle the issue of the low signal-to-noise ratio of deep reflections. Comparing (ratio) two stacked amplitudes is the same as comparing their average values. Averaging reduces the uncertainties (reduction by the square root of the sample size). The issue of true amplitude and move-out that arise in amplitude versus offset analysis could be considered here of a second order nature

given the depth of the imaging target (70 km) compared with the half-length of our streamer (6 km, full length is 12 km). As a consequence, the emergence angle of the rays will be small (less than 5° = tan⁻¹(6/70)), the move-out will be very small and stacking would not smear the amplitude significantly for the LAB. This approach is mainly valid for deep reflectors when the emergence angle is small.

The stacked seismic section was corrected for the spherical divergence and attenuation. The Earth intrinsic attenuation (A) is compensated using a quality factor (Q) of 1,000 as follows:

$$A = A_0 e^{-\frac{x\pi}{Q}}$$

where A₀ is the seafloor seismic amplitude, f is the dominant frequency and t is the two-way travel time.

For the upper reflector of the LAB, we picked the maximum negative amplitude within two wavelengths from the reflection (that is, ±160 ms) to account for manual picking errors. For the lower reflector, we did the same but with the positive amplitude. This process was also carried out for the seafloor.

We used a normal-incidence approximation since the maximum incidence angle for a target at about 70 km does not exceed 5° for a 12 km streamer. The normal-incidence reflection coefficient (R) in this case is:

$$R = \frac{\rho_2 v_2 - \rho_1 v_1}{\rho_2 v_2 + \rho_1 v_1}$$

The density in medium one and two is denoted by ρ₁ and ρ₂ respectively. The P-wave velocities are v₁ and v₂ in the first and the second media. One single average seafloor amplitude was used throughout the process, which is akin to considering that the variance of seafloor amplitude is much smaller than the variances of the LAB's amplitudes. The density variation between the lithosphere, the LAB and the asthenosphere is small. Therefore, we assume that the density is constant across the reflections, which simplifies the reflection coefficient to a ratio of velocities. We computed the reflection coefficient at each super CMP, produced a histogram and fitted a log-normal curve to estimate the velocity contrast and associated uncertainty (Supplementary Fig. 2). We use a maximum likelihood estimator (mode) and the uncertainty is fixed at one standard deviation.

**Melt estimation.** Melt estimation in the channel was carried out using the effective medium theory¹⁰. Melt inclusions are introduced incrementally into an olivine background matrix to produce a two-phase composite representative of the channel. The inclusions are introduced in the form of ellipsoids with semi-axes (a, c), which have a circular cross-section in the vertical plane (two of the three axes are equal). We consider ellipsoids with aspect ratios (a/c) varying from 1 (sphere) to 100 (thin film). To model the solid background material, we used the theoretically calculated elastic parameters of the olivine crystal from ref. ⁴⁴. We assumed that the inclusions are not connected to one another, which is valid for the low concentrations considered here. The properties of the melt are derived from ref. ⁴⁵. We made such a choice to be able to get the lowest estimation of the melt. The P-wave velocity is insensitive to the ellipsoid elongation above a a/c ratio of 100. On the other hand, the S-wave velocity remains sensitive to this ratio and decreases significantly for thinner inclusions²³. For the melt estimation we used the aspect ratio of 100, which yields a low estimation of the melt while keeping an acceptable shape of the melt structures. The P- and SV-wave velocities as a function of melt fraction are shown in Supplementary Fig. 6.

**Water fraction estimation for melt stabilization.** The isotherms defining the LAB are shallower than the dry solidus and therefore melt cannot exist unless volatiles such as water or carbon facilitate the melting. The weight fraction of water needed to stabilize the melt present in the LAB is calculated using the formulas from ref. ²⁶. A constant melt fraction of 1.4% wt is assumed throughout the channel. Carbon also facilitates the melting and reduces the water requirements⁶⁰; however, we do not consider the impact of carbon on melting since its concentration in the upper mantle is not very well constrained. The values computed here could be considered as upper limits for the required water fraction.

$$F(P, T, X_{\text{H}_2\text{O}}^{\text{bulk}}) = \left| \frac{T - \left( T_{\text{solidus}} - \Delta T \left( X_{\text{H}_2\text{O}} \left( X_{\text{H}_2\text{O}}^{\text{bulk}} P, F \right) \right) \right)^{1/\beta_1}}{T_{\text{herz}}^{\text{liquidus}} - T_{\text{solidus}}} \right|^{\beta_1}$$

T is the temperature (LAB's isotherms in this case), T<sub>solidus</sub> is the dry solidus temperature of olivine, T<sub>herz</sub><sup>liquidus</sup> is Iherzolite liquidus temperature, ΔT is the freezing point depreciation, P is the pressure, X<sub>H₂O</sub><sup>bulk</sup> is the bulk water fraction (wt %), X<sub>H₂O</sub><sup>solidus</sup> is the relative water fraction in the melt and β<sub>i</sub> is the power-law exponent of melt production. The values of these parameters are detailed in ref. ²⁶. The uncertainties are computed by applying the logarithm and differentiating the equations above. This method is chosen to avoid any assumption about uncertainty source correlation:
$$ \left| \frac{dX_{\text{H}_2\text{O}}^{\text{Bulk}}}{X_{\text{H}_2\text{O}}^{\text{Bulk}}} \right| \le \left| \frac{T-T_{\text{sol}}+\Delta T}{\gamma\beta_1\Delta T} - \frac{FD_{\text{H}_2\text{O}}}{D_{\text{H}_2\text{O}}+(1-D_{\text{H}_2\text{O}})F} \right| \left| \frac{dF}{F} \right| + \frac{1}{\gamma\Delta T} \left| F\frac{1}{\beta_1}\frac{\partial(T_{\text{liq}}^{\text{lherz}}-T_{\text{sol}})}{\partial P} + \frac{\partial T_{\text{sol}}}{\partial P} \right| \left| dP \right| + \frac{1}{\gamma\Delta T} \left| dT \right| $$

Supplementary Fig. 7b shows the water content as a function of pressure for different temperatures for a melt fraction of $F=1.4\%$.

The total mass of water for a given age (normalized per a surface unit section) is the surface delimited by the corresponding curve and the vertical axis. Numerical integration (along depth) of the water average value for a given age yields the average estimate of the total water mass. We assume that the uncertainties are perfectly correlated ($\rho_{corr}=1$) since the error sources (temperature, velocity, melt uncertainties) are the same; therefore, the standard deviation of the water integral (sum) is the integral (sum) of the standard deviations associated with each depth. Vertical integration along depth ($z$) or pressure in the LAB channel is the same since $dP=\rho dz$, and the density could be considered constant inside the channel.

**Water uncertainties in the LAB.** Two equations govern how much water is required to stabilize a given melt fraction$^{16}$.

$$ \left\{ \begin{aligned} F &= \left( \frac{T - T_{\text{sol}} + \Delta T}{T_{\text{liq}}^{\text{lherz}} - T_{\text{sol}}} \right)^{\beta_1} \\ \Delta T &= K \left( \frac{X_{\text{H}_2\text{O}}^{\text{Bulk}}}{D_{\text{H}_2\text{O}} + (1 - D_{\text{H}_2\text{O}})F} \right)^{\gamma} \end{aligned} \right. $$

Applying the logarithm we get:

$$ \left\{ \begin{aligned} \log F &= \beta_1 \log(T - T_{\text{sol}} + \Delta T) - \beta_1 \log(T_{\text{liq}}^{\text{lherz}} - T_{\text{sol}}) \\ \Delta T &= \log K + \gamma \log(X_{\text{H}_2\text{O}}^{\text{Bulk}}) - \gamma \log(D_{\text{H}_2\text{O}} + (1 - D_{\text{H}_2\text{O}})F) \end{aligned} \right. $$

Differentiation leads to:

$$ \left\{ \begin{aligned} \frac{dF}{F} &= \beta_1 \frac{d(T - T_{\text{sol}} + \Delta T)}{T - T_{\text{sol}} + \Delta T} - \beta_1 \frac{d(T_{\text{liq}}^{\text{lherz}} - T_{\text{sol}})}{T_{\text{liq}}^{\text{lherz}} - T_{\text{sol}}} \\ \frac{d\Delta T}{\Delta T} &= \gamma \frac{dX_{\text{H}_2\text{O}}^{\text{Bulk}}}{X_{\text{H}_2\text{O}}^{\text{Bulk}}} + \gamma \frac{D_{\text{H}_2\text{O}} dF}{D_{\text{H}_2\text{O}} + (1 - D_{\text{H}_2\text{O}})F} \end{aligned} \right. $$

These two equations can be reorganized as follows:

$$ \left\{ \begin{aligned} d\Delta T &= \frac{(T - T_{\text{sol}} + \Delta T)}{\beta_1 F} dF + \frac{T - T_{\text{sol}} + \Delta T}{T_{\text{liq}}^{\text{lherz}} - T_{\text{sol}}} d(T_{\text{liq}}^{\text{lherz}} - T_{\text{sol}}) - d(T - T_{\text{sol}}) \\ d\Delta T &= \gamma \frac{dX_{\text{H}_2\text{O}}^{\text{Bulk}}}{X_{\text{H}_2\text{O}}^{\text{Bulk}}} \Delta T + \gamma \frac{D_{\text{H}_2\text{O}} dF}{D_{\text{H}_2\text{O}} + (1 - D_{\text{H}_2\text{O}})F} \Delta T \end{aligned} \right. $$

Combining these two equations yields:

$$ \gamma \frac{dX_{\text{H}_2\text{O}}^{\text{Bulk}}}{X_{\text{H}_2\text{O}}^{\text{Bulk}}} \Delta T = \left( \frac{(T - T_{\text{sol}} + \Delta T)}{\beta_1} - \gamma \frac{FD_{\text{H}_2\text{O}} \Delta T}{D_{\text{H}_2\text{O}} + (1 - D_{\text{H}_2\text{O}})F} \right) \frac{dF}{F} + F \frac{1}{\beta_1} d(T_{\text{liq}}^{\text{lherz}} - T_{\text{sol}}) - d(T - T_{\text{sol}}) $$

Applying the norm and majoring gives the following inequality:

$$ \left| \frac{dX_{\text{H}_2\text{O}}^{\text{Bulk}}}{X_{\text{H}_2\text{O}}^{\text{Bulk}}} \right| \le \left| \frac{(T - T_{\text{sol}} + \Delta T)}{\gamma \beta_1 \Delta T} - \frac{FD_{\text{H}_2\text{O}}}{D_{\text{H}_2\text{O}} + (1 - D_{\text{H}_2\text{O}})F} \right| \left| \frac{dF}{F} \right| + \frac{1}{\gamma \Delta T} \left| F \frac{1}{\beta_1} \frac{\partial(T_{\text{liq}}^{\text{lherz}} - T_{\text{sol}})}{\partial P} + \frac{\partial T_{\text{sol}}}{\partial P} \right| \left| dP \right| + \frac{1}{\gamma \Delta T} \left| dT \right| $$

For the temperature $dT$ we have:

$$ \left\{ \begin{aligned} T &= 1,260 \,^{\circ}\text{C} \ dT \sim \pm 20 \,^{\circ}\text{C} \\ T &= 1,355 \,^{\circ}\text{C} \ dT \sim \pm 15 \,^{\circ}\text{C} \end{aligned} \right. $$

For the melt $dF$:

$$ F = 1.4\% \ dF \sim \pm 0.8\% $$

The pressure uncertainties $dP$ are related to the depth conversion:

$$ dP = \rho dz $$

The uncertainties for $dz$ are $\pm 2$ km for upper reflector and $\pm 3$ km for the lower reflector. To avoid any rounding issues, all depth uncertainties are fixed at $\pm 3.5$ km. $\rho$ is fixed at $3.3 \text{ kg m}^{-3}$.

**Viscosity estimation.** The effect of the melt on viscosity is quantified through the viscous contiguity model$^{38,47}$.

$$ \eta = 0.2 B \phi^2 \eta_0 $$

$$ \varphi = 1 - A \Phi^{1/2} \quad (A = 2.3) $$

$$ B = 0.12 $$

$\eta$ and $\eta_0$ are shear viscosities with and without melt. $\varphi$ is the contiguity parameter that depends on the melt fraction $\Phi$, which could be either weight- or volume-based, as the difference between weight and volume fractions of melt is negligible$^{47}$.

The melt-free viscosity can be estimated using the following equation:

$$ \eta = \frac{\sigma}{\epsilon} = \frac{A^{-1}}{\sigma^{n-1}} \exp\left(\frac{Q + PV}{RT}\right) $$

$\sigma$ is the differential stress, $\dot{c}$ is the strain rate, $A$ is multiplicative factor, $n$ is the stress exponent, $Q$ is the activation energy, $P$ is the pressure, $V$ is the activation volume, $R$ is the perfect gas constant and $T$ is the temperature in Kelvin. We used the parameters $(A, n, Q)$ from ref. $^{48}$.

The lithosphere viscosity is estimated to be around $10^{21}$ (Pa. s) at its base and drops sharply to $10^{18}$ (Pa. s) in the LAB due to the presence of melt.

**Data availability.** The authors declare that the data supporting the findings of this study are available within the article and its Supplementary Information. The processed seismic data used in this study will be publicly available at the end of the ERC Advance Grant agreement no. 339442\_TransAtlanticILAB, from 2020.

## References

41. Vassallo, M., Eggenberger, K., van Manen, D. J., Özbek, A. & Watterson, P. Broadband and beyond with marine towed streamers. *Leading Edge* **32**, 1356–1365 (2013).

42. Bekara, M. & Van der Baan, M. Random and coherent noise attenuation by empirical mode decomposition. *Geophysics* **74**, V89–V98 (2009).

43. Qin, Y. & Singh, S. C. Seismic evidence of a two-layer lithospheric deformation in the Indian Ocean. *Nat. Commun.* **6**, 8298 (2015).

44. Abers, G. A. & Hacker, B. R. A MATLAB toolbox and Excel workbook for calculating the densities, seismic wave speeds, and major element composition of minerals and rocks at pressure and temperature. *Geochem. Geophys. Geosyst.* **17**, 616 (2016).

45. Clark, A. N., Lesher, C. E., Jacobsen, S. D. & Wang, Y. Anomalous density and elastic properties of basalt at high pressure: Reevaluating the effect of melt fraction on seismic velocity in the Earth’s crust and upper mantle. *J. Geophys. Res. Solid Earth* **121**, 4232–4248 (2016).

46. Dasgupta, R. et al. Carbon-dioxide-rich silicate melt in the Earth’s upper mantle. *Nature* **493**, 211 (2013).

47. McCarthy, C. & Takei, Y. Anelasticity and viscosity of partially molten rock analogue: Toward seismic detection of small quantities of melt. *Geophys. Res. Lett.* **38**, L18306 (2011).

48. Hirth, G. & Kohlstedt, D. L. Water in the oceanic upper mantle: implications for rheology, melt extraction and the evolution of the lithosphere. *Earth Planet. Sci. Lett.* **144**, 93–108 (1996).