Introduction

Earth is the only known planet featuring plate tectonics1, whose main driving force is the gravitational pull of the oceanic lithosphere sinking into the mantle, namely: slab subduction2,3. Water stored in the slab is fundamental to the onset of subduction, and is thus crucial to sustain plate tectonics4, by reducing solidus temperature5 and mechanical strength6 of the rocks. Numerical models7 suggest that ~2/3 of the water stored inside the slab is released before reaching 230 km of depth due to the decomposition of temperature sensitive hydrous minerals8. Dehydration processes trigger flux melting5, leading to arc volcanism9, and dehydration embrittlement6, which is responsible for intermediate depth seismicity10 (70−300 km). Some hydrous minerals, however, can survive longer inside the cold core of the slab11, potentially forming a deep water reservoir8 in the region extending from ~410 to ~660 km of depth, also known as Mantle Transition Zone (MTZ)12. Modelling slab thermal evolution and the depth of water release in the mantle is therefore essential for understanding the deep-Earth water cycle and plate tectonics.

Slab thermal evolution ∂T/∂t is usually computed with numerical models13, using Fourier’s Law14 of heat conduction \(Q=-\varLambda \nabla T\). This equation describes temperature changes as a function of the amount of thermal energy that enters or leaves an object i.e., the heat flux \(Q \ [W{m}^{-2}]\). The key parameter in this equation is the thermal conductivity, Λ  \([W{m}^{-1}{K}^{-1}]\), as it quantifies the amount of heat transported. In electrical insulators, such as silicate minerals, the heat propagates as acoustic vibrations of the crystal lattice (phonons)15, and as optical radiation (photons)16. The total thermal conductivity Λ is the sum of the lattice (\({\varLambda }_{{lat}}\)) and radiative (\({\varLambda }_{{rad}}\)) components: \(\varLambda={\varLambda }_{{lat}}+{\varLambda }_{{rad}}\). Importantly, the physical properties of subducting rocks change when the slab is subjected to high pressure \(P\left[{Pa}\right]\) and temperature \(T\left[K\right]\) conditions of the Earth’s interior17. The effects of increasing P, T on lattice heat transport of major upper mantle minerals are well known: P increases \({\varLambda }_{{lat}}\) by increasing mean phonon velocity18, whereas T reduces \({\varLambda }_{{lat}}\) by increasing phonons collisions19. In contrast, radiative heat transport increases rapidly with \({T}^{3}\)(Stephan-Boltzmann’s law)16, but decreases with increasing opacity of minerals, which itself is a function of P and T20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37. However, most geodynamic models do not consider the effects of pressure and temperature on lattice and radiative heat transport13, generally assuming a constant \({\varLambda }_{{lat}}\) and neglecting \({\varLambda }_{{rad}}\).

One of the reasons as to why radiative heat transport has been largely ignored in geodynamic models is because early studies23,24 reported that iron-bearing silicates minerals are likely opaque at the high P, T conditions of the upper mantle. In order to estimate the contribution of \({\varLambda }_{{rad}}\) of representative mantle minerals, it is crucial to probe the optical absorption coefficient in the wavelength range that includes the peak of thermal emission at mantle temperatures 1000 < T < 2500 K, i.e., the infrared (IR) range between 1.1 and 2.8 μm. Several studies have reported estimates of \({\varLambda }_{{rad}}\) for relevant iron-bearing silicate minerals: olivine22,23,24,25,26,27,28,29,30,31,32,33,34,35, pyroxene27,31,35, garnet27,31,35, wadsleyite36, ringwoodite20,36, bridgmanite21,37, and ferropericlase37. Most of these estimates, however, vary by several orders of magnitude, possibly because measurements were performed on samples of variable composition, optical quality, and at different P, T conditions. Quantitative estimates of the radiative component at high P, T conditions are very challenging because they require optical measurements on very small samples in the IR spectral range16,38. Olivine (Mg,Fe)2SiO4, the most abundant mineral in the slabs at upper mantle conditions, is a case in point. Even at room pressure, the radiative conductivity of olivine is poorly constrained with extant estimates ranging from 0.3127 to 2.3422 W m−1 K−1 (P = 1 atm, T ≤ 1300 K). At room pressure, crystal field bands (d-d excitation of Fe2+ atoms) are a prominent absorption feature in olivine in the near-IR range23,24. Measurements performed at high-T and room-P revealed that the crystal field bands broaden, and the lattice bands in the mid-IR intensify with increasing temperature22,25. Temperature-induced strengthening of the lattice bands is particularly important because these bands may effectively block radiative heat transport at upper mantle temperature. In addition, room-T, high-P (6−27 GPa) optical experiments on fayalite (i.e., the iron endmember of olivine, Fe2SiO4) revealed that the Fe-O absorption edge shifts into the visible and near-IR range, and becomes the dominant light absorption mechanism, potentially blocking radiative heat transport23,24. These studies, however, did not clarify which absorption mechanism(s) govern the optical properties of olivine at upper mantle P, T conditions due to the lack of experimental spectral coverage and/or limited P, T range27,28,29,30.

In this work, we probed, for the first time, the optical absorption coefficient of olivine (\({\alpha }^{{Ol}}\)) in the IR and visible spectral ranges at simultaneous high \(P\) and \(T\). This achievement was enabled by the use of a pulsed white laser probe synchronized to fast spectroscopic equipment in an optical setup for laser-heating diamond anvil cell (LH-DAC) experiments37,38,39,40. From our \({\alpha }^{{Ol}}\) measurements, we derived reliable estimates of the radiative thermal conductivity of olivine \(({\varLambda }_{{rad}}^{{Ol}})\) at upper mantle \(P\), \(T\) conditions. Our experiments reveal that olivine remains relatively transparent to IR light even at high \(P\),\(T\), which makes \({\varLambda }_{{rad}}^{{Ol}}\) a non-negligible component of diffusive heat transport in the mantle. Additionally, we derived a new formulation to compute \(T\)-dependent \({\varLambda }_{{rad}}^{{Ol}}\), which we introduced in a heat diffusion model to investigate the effects of significant radiative heat transport in the thermal evolution of subducting slabs. Recent thermal evolution models of subducting slabs13,41,42,43,44,45,46,47,48,49 included \(P\),\(T\)-dependent thermal conductivity formulations, and reported a large temperature difference (\(\Delta T=100\)–\(200 \ K\)) compared to the models with constant \(\varLambda\). Such large \(\Delta T\) has dramatic effects on subduction dynamics43,44 (i.e., slab buoyancy and sinking velocity), and affects the preservation of hydrous phases inside the slab49. In previous subduction models, however, the contribution of \({\varLambda }_{{rad}}\) was fixed to small values (<\(1.5\) \(W{m}^{-1}{K}^{-1}\)), which leads to negligible effects on subduction dynamics13,41,46,48. The inclusion of high \({\varLambda }_{{rad}}^{{Ol}}\), such as inferred in this work, might radically change the internal thermal conductivity distribution inside the slab. Our models show that the inclusion of both lattice and radiation heat diffusion mechanisms in geodynamic models is crucial for the correct computation of slab thermal evolution and the estimation of the amount of water released at \( < 230 \ {km}\)7,10 or delivered to the MTZ (\({{\boldsymbol{ \sim }}}410\)–\(660 \ {km}\))8.

Results

The absorption coefficient of olivine at high-P, T

The absorption coefficient of fayalite \({\alpha }^{{Fa}100}\) measured in the LH-DAC at ~1.1 \({GPa}\) is shown in Fig. 1. Across all frequencies probed, \({\alpha }^{{Fa}100}\) increases with heating. This increase is fully reversible upon quenching, as indicated by the room-temperature spectra measured after each heating cycle, suggesting no significant irreversible chemical alteration in fayalite at high temperatures. The agreement with \({\alpha }^{{Fa}100}\) measured at 1 \({atm}\) in a furnace30 is excellent, indicating a homogeneous temperature in the probed laser-heated fayalite. At pressures of ~5 and ~ 10 \({GPa}\), the measured absorption coefficients are very similar in magnitude to those at ~1.1 GPa, albeit the crystal field bands are blue-shifted, consistent with the expected pressure-induced shift24,26. The wavelength-dependent temperature-derivatives of \({\alpha }^{{Fa}100}\) measured in LH-DACs increase with temperature at all pressures studied here (Fig. S3). Similarly, the absorption coefficient, and its temperature-derivative for natural peridotitic olivine (Fa9.1) determined in the LINKAM experiment, also increase with temperature (Fig. S4). Using representative temperature-derivatives from the LH-DAC, the LINKAM experiments, and the data reported by ref. 25 we obtained a continuous in wavelength temperature-derivative of the absorption coefficient of peridotitic olivine (see Supplementary Materials, SM).

Fig. 1: The absorption coefficient of fayalite αFa100 at ~1.1  GPa, measured in the laser-heated diamond anvil cell (solid spectra).
figure 1

Eight overlapping black and blue spectra were measured before and after each of the high-temperature spectra (color-labeled), respectively. The dashed spectra are absorption coefficients of fayalite at room- and high temperature reported by ref. 30. The ICCD and InGaAs labels denote the spectral ranges recorded by these detectors. The vertical grey bar shows the region where our optical measurements are inaccurate due to the notch filter blocking the heating laser.

Full size image

Our results are in agreement with the extant studies of olivine opacity at frequencies > ~4000–5000 cm−1 25,28,29,30, in particular with the data reported by ref. 25 where they probed \({\alpha }^{{Ol}}\) at \(P=1\) \({atm}\) and \(T\le 1673\ K\) in the near- and mid- IR spectral range. The agreement with the optical data of ref. 25 indirectly testifies to the reliability of our measurements. In addition, we observe that the absorption coefficient at <4000 \({{cm}}^{-1}\) strongly increases with temperature due to the broadening of the lattice vibration bands, as also reported by ref. 25. This is an important observation that challenges the assumption made by ref. 28, that only the d-d crystal field bands contribute to the temperature-induced opacity of olivine at \(T > \) 673 \(K\). The major differences between our \({\alpha }^{{Ol}}\) measurements and most previous room-\(P\), high-\(T\) studies28,29,30 arise from their limited spectral coverage30, the investigated \(T\) range (e.g., \(T < 673 \ K\)28), or both29. At high temperatures, the optical anisotropy of olivine’s absorption coefficient \({\alpha }^{{Ol}}\) 50,51,52 results in only ~10% different \({\varLambda }_{{rad}}^{{Ol}}\)25 (Fig. 2, diamond markers). Therefore, the unpolarized optical measurements reported here enable an adequate assessment of the radiative thermal conductivity of olivine at upper mantle temperatures. We further discuss the anisotropy of \({\alpha }^{{Ol}}\) in the SM.

Fig. 2: Thermal conductivities of the Earth’s outer shell: 40 km-thick lithosphere, upper mantle (40–410 km), upper MTZ (410–520 km), and lower MTZ (520–660 km).
figure 2

Our calculated mantle adiabat53,54 assumes a potential temperature of \(1619 \ {K}\). The solid curves represent the \(\varLambda\) profiles used in our study, while the dashed and dashed-dotted curves indicate the \(\varLambda\) profiles reported in the literature. The \(P\), \(T\) dependence of \({\varLambda }_{{rad}}^{{Ol}}\) is shown in Fig. S5. 1. (solid black) total thermal conductivity computed as \(\varLambda={\varLambda }_{{lat}}+{\varLambda }_{{rad}}\). The grey shaded area represents the uncertainty in \(\varLambda\) estimates. 2. (solid red) radiative thermal conductivity \({\varLambda }_{{rad}}\) computed from: olivine[this study] (UM), wadsleyite36 (upper MTZ), and ringwoodite36 (lower MTZ). The red shaded area represents the \(\pm 30\%\) uncertainty in \({\varLambda }_{{rad}}^{{Ol}}\) estimates from LH-DAC measurements37,38,39. 3. (solid blue) lattice thermal conductivity \({\varLambda }_{{lat}}\) profile computed from: olivine55,68 (UM) and ringwoodite45 (MTZ). The blue shaded area represents \(\pm 15\%\) uncertainty in \({\varLambda }_{{lat}}^{{Ol}}\) estimates from Time-Domain Thermo-Reflectance (TDTR) measurements45,68. 4. (dashed black) total thermal conductivity \({\varLambda }_{{tot}}^{{Hof}1999}\), with \({\varLambda }_{{rad}}^{{Hof}1999} \sim 0.35\) \(W\) \({m}^{-1}\) \({K}^{-1}\) ref. 27. 5. (dashed-dotted crimson) radiative thermal conductivity \({\varLambda }_{{rad}}^{{Hof}2005}\), ref. 28. 6. (dashed crimson) radiative thermal conductivity \({\varLambda }_{{rad}}^{G\&A2019}\), ref. 35. 7. (coral diamonds) radiative thermal conductivity for the crystallographic directions a and c \({\varLambda }_{{rad}}^{Sha1979}\), ref. 25. Note that our estimated \({\varLambda }_{{rad}}\) (profile 2, solid red) rapidly increases in the lithosphere (linear \(T\) gradient \(35\) \(K/{km}\)) and flattens in the sub-lithospheric mantle (adiabatic \(T\) gradient \(\sim \!0.5\) \(K/{km}\)54).

Full size image

Olivine’s \({{\boldsymbol{P}}}\), \({{\boldsymbol{T}}}\)-dependent radiative thermal conductivity

From our absorption coefficient measurements, we computed \({\varLambda }_{{rad}}^{{Ol}}\) (SM) along the \(P\), \(T\) profile of the upper mantle53,54. We find that the radiative conductivity increases steeply below \(40 \ {km}\), but is \(\sim \!2.3 \ W{m}^{-1}{K}^{-1}\) throughout the upper mantle (Fig. 2). The relative error of our estimates is \({\varLambda }_{{rad}}^{{Ol}}\) \(\pm 30\%\). This is an empirical estimate that takes into account how the uncertainty in thickness and temperature of the sample under high P37,38,39 propagate to the error in the absorption coefficient. The uncertainty in temperature is particularly important because the absorbance of the sample is measured across a volume subject to a strong temperature gradient along the laser-heating direction. We find that the radiative component of the total thermal conductivity of olivine \({\varLambda }^{{Ol}}\) is almost insensitive to pressure but strongly dependent on temperature (Fig. S5). Moreover, our estimates demonstrate that heat radiation is the predominant mechanism for heat diffusion in the upper part of the sub-lithospheric mantle (\( < \!150 \ {km}\)), where \({\varLambda }_{{rad}}^{{Ol}} > \!0.5{\varLambda }^{{Ol}}\). Within the lithosphere (\( < \!40 \ {km}\)), the lattice component \({\varLambda }_{{lat}}^{{Ol}}\) (Fig. 2) is strongly damped55 due to the high-\(T\), low-\(P\) conditions. In the sub-lithospheric mantle, however, the pressure gradient exceeds the temperature gradient, hence \({\varLambda }_{{lat}}^{{Ol}}\) becomes progressively more relevant for \({\varLambda }^{{Ol}}\) with increasing depth (Fig. 2). Nevertheless, the high temperatures sustain intensive radiative heat transport down to a depth of \(410 \ {km}\), where \({\varLambda }_{{rad}}^{{Ol}} \sim 0.4{\varLambda }^{{Ol}}\).

Our radiative conductivity model for olivine is in reasonable agreement with that of ref. 25 but differs substantially from the estimates obtained from other formulations in the literature27,28,35 (Fig. 2). Please note that the formulation in ref. 35 is based on the optical absorption spectra reported by ref. 28, who assumed that only the d-d crystal field bands contribute to the temperature-induced opacity of olivine at \(T > \) 673 \(K\), an assumption that is not supported by our LINKAM experiments. As a result, the \(T\)-derivative of \({\varLambda }_{{rad}}^{{Ol}}\) from ref. 35 is ~2.5 times steeper than that from our work at \(T < 1200 \ K\). From our estimates, we conclude that olivine’s radiative and lattice conductivities are of approximately equal importance at upper mantle conditions.

Discussion

The effect of olivine’s radiative thermal conductivity on slab thermal evolution

The inclusion of \({\varLambda }_{{rad}}\) has a significant effect on the internal temperature distribution of the slab (Fig. 3). Each model shows a cold core (\(T\) < 1000 \(K\)) surrounded by hot outer edges56. However, accounting for radiative transport in heat diffusion (model set 02; \({\varLambda }_{{rad}} > 0\)) leads to slabs that are\(\sim 100\)–\(200 \ K\) warmer than in the simulations with only the lattice component (model set 01; \({\varLambda }_{{rad}}=0\)). This tendency can be observed in all models (Fig. S14). Using bound values for thermal conductivity (i.e., upper and lower limits of the grey shaded area in Fig. 2) will produce a \(\Delta T\) of\(\pm 50 \ K\) compared to the reference model in Fig. 3b (\({t}_{{slab}}=80\); \({v}_{{sink}}=5 \ {cm}\) \({{yr}}^{-1}\); \({\varLambda }_{{rad}} > 0\)), whereas the formulation from ref. 35 (eq. S12 in SM) produces \(\sim 20 \ K\) hotter slab than the reference in Fig. 3b.

Fig. 3: Examples of final 2D temperature distributions from a slab model with sinking velocity of vsink = 5 cm yr−1 and age of tslab = 80 Myrs.
figure 3

Panel (a) represents the model without radiative heat transport; panel (b) the model with radiative heat transport; panel (c) shows the temperature difference between (b) and (a). Each subplot plot shows only the upper \(40\) \({km}\) of the slab, starting from its surface at \(0\) \({km}\). For the plot we used the Scientific Color Maps56. The temperature profile across the whole slab thickness is reported in Fig. S15.

Full size image

The inclusion of high \({\varLambda }_{{rad}}\) formulation in heat diffusion models intensifies the heat flow leading to faster and deeper slab heating compared to \({\varLambda }_{{lat}}\)-only cases. The \(T\)-derivatives of \({\varLambda }_{{rad}}^{{Ol}}\) and \({\varLambda }_{{lat}}^{{Ol}}\) have opposite sign: \({\left(\frac{\partial {\varLambda }_{{rad}}}{\partial T}\right)}_{P} \sim\)1.61 × 10−3\(W{m}^{-1}{K}^{-2}\); \({\left(\frac{\partial {\varLambda }_{{lat}}}{\partial T}\right)}_{P} \sim\)−1.03 × 10−3 \(W{m}^{-1}{K}^{-2}\). Therefore, as temperature increases, radiative heat transport is enhanced and overcompensates for the reduction of lattice conductivity induced by \(T\) (Fig. 2), which gives a positive \(T\)-derivatives of total thermal conductivity \({\varLambda }_{{tot}}^{{Ol}}:{\left(\frac{\partial {\varLambda }_{{tot}}}{\partial T}\right)}_{P} \sim\) 0.58 × 10−3 \(W{m}^{-1}{K}^{-2}\). The positive feedback between temperature and \({\varLambda }_{{rad}}\) maintains high conductivity within the slab and smooths \(\varLambda\) variations across its thickness (Figs. S16, S17). The effect of pressure on \({\varLambda }_{{rad}}\), instead, is negligible \({\left(\frac{\partial {\varLambda }_{{rad}}}{\partial P}\right)}_{T} \sim\)0 \(W{m}^{-1}{K}^{-1}{{GPa}}^{-1}\), while the P-derivative of lattice conductivity is significant \({\left(\frac{\partial {\varLambda }_{{lat}}}{\partial P}\right)}_{T} \sim\)0.55 \(W{m}^{-1}{K}^{-1}{{GPa}}^{-1}\). Therefore, the variation in \({\varLambda }_{{tot}}^{{Ol}}\) along the slab dip is almost exclusively due to the increasing \({\varLambda }_{{lat}}(P)\), hence \({\left(\frac{\partial {\varLambda }_{{tot}}}{\partial P}\right)}_{T} \sim\)0.55 \(W{m}^{-1}{K}^{-1}{{GPa}}^{-1}\). It follows that the horizontal (lateral) variations of \(\varLambda\) within a subducting slab are significantly smaller than the vertical (radial) ones, as also reported by ref. 13. In previous models13,41,46,48, however, the contribution of \({\varLambda }_{{rad}}\) was not sufficient to compensate for the diminishing \({\varLambda }_{{lat}}\), and large vertical variations of thermal conductivity were mainly generated by \(\varLambda\) jumps across phase boundaries13. In our models, in contrast, \({\varLambda }_{{rad}}\) contribution attenuates the effects of temperature on slab’s total thermal conductivity due to the lattice component, thus making more apparent the effects of pressure given by \({\varLambda }_{{lat}}\) contribution.

Large-scale impact of \({{{\boldsymbol{\Lambda }}}}_{{{\boldsymbol{rad}}}}\) on slab dehydration and subduction dynamics

Previous models43,44 reported that \(P\)-dependent \({\varLambda }_{{tot}}\) significantly decelerate subduction rate by increasing slab’s thermal buoyancy. Potentially, including the contribution of strongly T-dependent \({\varLambda }_{{rad}}\) into \({\varLambda }_{{tot}}\) will dramatically alter the whole slab dynamics. In our models, the largest temperature difference between the two cases (\({\varLambda }_{{rad}}=0\) and \({\varLambda }_{{rad}} > 0\)) is located in the region extending from the former ocean floor to \({{\boldsymbol{ \sim }}}40 \ {km}\) towards the slab core (Fig. 3c, S15). At the trench, this region is the coldest part of the oceanic lithosphere, experiencing the largest temperature gradient once subduction begins. Moreover, the uppermost part of the slab is characterized by the greatest diversity of lithologies57,58 and the highest viscosity59. Most importantly, this region contains the hydrous minerals formed by seafloor alteration11.

Hydrous minerals are particularly temperature-sensitive and only survive in the cold regions of the slab8. Dehydration is thus triggered at shallower depth in slabs that heat up more rapidly. To quantify slab dehydration, we extracted the final \(P\), \(T\) profiles from our models, and tracked the sequence of phase reactions60 that take place in a simplified hydrous harzburgite mineral assemblage: olivine, enstatite, and antigorite (SM). Importantly, in the \({\varLambda }_{{lat}}\)-only models, the \(P\), \(T\) profiles taken at \(10 \ {km}\) from the slab surface (roughly corresponding to the Moho in the oceanic lithosphere) cross only partial dehydration reactions (i.e., the solid reaction products are not all anhydrous), allowing the slabs to deliver \(13\)–\(25\%\) of their original water content at the base of MTZ (\(660 \ {km}\)) for all tested subduction velocities. In contrast, models that include radiative heat transport are ~100–150 K warmer (Fig. S14). This has dramatic consequences: the most superficial part (\(0\)–10 km) of slow and moderately fast slabs (\({v}_{{sink}} < 10{cm}\) \({{yr}}^{-1}\)) dehydrate completely before reaching the MTZ (Fig. 4a). Nonetheless, the hydrous minerals can still reach the MTZ even in \({\varLambda }_{{rad}}\)-including models when they are hosted in the coldest region of the slab, \(\sim 30 \ {km}\) below its surface (Fig. 4b). Similar trends are seen in the models of \(60\)–\(80 \ {Myrs}\) old slabs, whereas 20–40 \({Myrs}\) old slabs – which are thinner and initially hotter – lose their entire water budget outside of the MTZ (Fig. 4, S18–S20; Extended Data Table S8). Compared to the reference model (Fig. 3b), the \(50 \ K\) hotter slabs produced with +30% \({\varLambda }_{{rad}}\) uncertainty lead to \(16 \ {km}\) shallower dehydration (Fig. S13). Similar dehydration depths are obtained when using \({\varLambda }_{{rad}}\) formulation from ref. 35. On the other hand, the \(50 \ K\) colder slab produced with the −30% \({\varLambda }_{{rad}}\) uncertainty leads to a \(25 \ {km}\) deeper dehydration.

Fig. 4: Depth of dehydration reactions as a function of sinking velocity vsink.
figure 4

Slab age is \({t}_{{slab}}=80\) \({Myrs}\), slab thickness \({H}_{{slab}}=120\) \({km}\). Different markers indicate the type of dehydration reaction (star: complete dehydration; triangle: partial dehydration), whereas the colors indicate the type of model (\({\varLambda }_{{rad}}=0\): blue; \({\varLambda }_{{rad}} > 0\): red). The light blue area indicates the MTZ (\(410\)–\(660\) \({km}\)). The two subplots indicate different positions of extracted P, T profiles inside the slab: (a) \(10\) \({km}\) from slab surface (Moho); (b) \(30\) \({km}\) from slab surface (coldest region within the slab corresponding to maximum hydration depth)11. Partial dehydration reactions determine the loss of 50–70% of the initial slab water content at <350 \({km}\) of depth. The remaining 30% is lost between 640 and 660 \({km}\) of depth (Fig. S13 in SM).

Full size image

Heating enhanced by the radiative mechanism can have important consequences for the dynamics of the slab, for example, an increase of +100 \(K\) due to radiative heating, reduces temperature-dependent viscosity59 by 2 orders of magnitude, thus favoring plastic deformation in the slab (Figs. S21–S22). Similarly, models that include radiative heat transport are characterized by \(15\)–\(35 \ {kg}\) \({m}^{-3}\) lower density compared to \({\varLambda }_{{lat}}\)-only models (Figs. S23–S24), which corresponds to a \(0.3\)–\(1.0\%\) decrease in total slab density (\({{\boldsymbol{ \sim }}}3500\) \({kg}\) \({m}^{3}\)). A \(+25\) \({kg}\) \({m}^{-3}\) density contrast is sufficient to maintain subduction61, therefore a \(-35 < \Delta \rho < \)-\(15\) \({kg}\) \({m}^{-3}\) contrast should considerably reduce slab descending velocity, as also reported by refs. 43,44. The effects of \({\varLambda }_{{rad}}\) on the dynamics of subducting lithosphere warrant detailed future thermo-mechanical modelling.

The experimental breakthrough presented in our work informed a new model of thermal conductivity that properly accounts for the lattice \({\varLambda }_{{lat}}\) and radiative \({\varLambda }_{{lat}}\) contributions in the upper mantle. This study demonstrates that olivine is infrared transparent even at the high \(P\), \(T\) conditions of the upper mantle. It follows that the radiative mechanism significantly contributes to heat diffusion in olivine-rich upper mantle, representing up to ~40% of its total thermal conductivity. For this reason, \({\varLambda }_{{rad}}\) plays an important role in slab heating, especially as temperature rises during subduction. Our thermo-kinematic models support this claim and provide solid evidence that \({\varLambda }_{{rad}}\) can have far-reaching effects on slab subduction. In particular, we show that the rapid heating enhanced by radiation reduces the breakdown depth of hydrous phases, potentially triggering intermediate-depth seismicity in the slab. In addition, the faster heating might also promote slab bending and stagnation62 in the MTZ by reducing rock density and viscosity. We therefore suggest that future numerical models of slab subduction should account for radiative heat transport.

Methods

Optical measurements in LH-DAC

We probed the temperature-dependence of the optical absorption coefficient of olivine in two types of experiments (SM). In the first, we measured \({\alpha }^{{Ol}}\) of non-oriented \(350\)-\(\mu m\)-thick double polished samples of natural peridotitic olivine (Fa9.1) at \(1\) \({atm}\) and up to \(1273\) \(K\) using a high-temperature LINKAM T1500 stage. In the second, we measured \({\alpha }^{{Ol}}\) of synthetic non-oriented ~\(20\)-\(\mu m\)-thick single crystals of fayalite (Fa100) in a laser-heated diamond anvil cell (LH-DAC) with an unpolarized broadband (supercontinuum) laser probe (\(4000\)–\(24000\) \({{cm}}^{-1}\))37,38,39. LH-DAC experiments were carried out at \(P \sim \!1.1,\sim \!5,\sim \!10\) \({GPa}\) and \(T\) up to \(1210 \ K\). The use of fayalite was necessary for reliable detection of the absorption bands in thin (\(\sim \!\!15\)–\(20\) \(\mu m\)) samples. The absorption coefficient of fayalite, scaled for iron content, shows that fayalite is a suitable model for mantle olivine63 (Fig. S1). Essential details of Fa100 synthesis and its suitability as a model for mantle olivine are described in the SM. Radiative thermal conductivity of a mantle (Fe-poor) olivine \({\varLambda }_{{rad}}^{{Ol}}\) was then computed using optical data from both types of experiments (SM).

Numerical model

To quantify the effects of olivine \({\varLambda }_{{rad}}^{{Ol}}\) on the thermal evolution of subducting slabs13, we performed a series of numerical models using a MATLAB (R2022a)64 code (SM), that solves the 2D heat diffusion equation with the finite difference method65. This code is based on a previously developed algorithm49. In our models, we defined a rectangular \({L}_{{slab}}\times {H}_{{slab}}\) slab, which was discretized into a mesh of square cells, with a constant grid spacing of \(500 \ {m}\). Slab length was set to \({L}_{{slab}}=600 \ {km}\), whereas slab thickness ranged between \(60\le {H}_{{slab}}\le 120 \ {km}\) as it was calculated from the slab age: \(20\le {t}_{{slab}}\le 80 \ {Myrs}\) (eq. S25). The initial temperature profile of the slab upon entering the mantle was computed with the half-space cooling equation66 (eq. S26). In the model, the slab subducted vertically (dip angle \({\theta }_{{dip}}\)= 90°) at a prescribed velocity (\({v}_{{sink}}\)), while heat flowed from the hot ambient mantle into the cold slab interior as a function of \(\varLambda (P,T)\). Subduction occurred from \(60\)–\(120 \ {km}\) depth (depending on \({t}_{{slab}}\)) down to \(660 \ {km}\), i.e., the base of the MTZ. In the model, we ignored the interactions between the subducting slab and the overriding plate. During subduction, \(P\), \(T\) were increased following the PREM53 and the adiabatic temperature profile of the mantle reported by ref. 54. The \(T\)-dependent heat diffusion equation is non-linear and requires the re-computation of the thermodynamic properties of the slab using the \(P\), \(T\) conditions of the current time-step. The equations used to compute \(P\), \(T\)-dependent thermal conductivity Λtot [Wm−1 K−1], density ρ [kg m−3], specific heat capacity Cp [J kg−1 K−1], and heat diffusivity \(\kappa \ [{m}^{2}{s}^{-1}]\) are reported in the SM (eqs. S6–S9; S13–S24; Extended Data Tables 1, 4). We restricted the mineralogical composition of the slab to olivine and its polymorphic modifications wadsleyite (\({Wd}\)) and ringwoodite (\({Rw}\)), which is adequate because the volumetric fraction of these minerals is \(\sim 60\%\) in the upper mantle12 and \(\sim 80\%\) in the oceanic lithosphere67. Therefore \({\varLambda }^{{slab}}=\) \({\varLambda }_{{lat}}^{i}(P,T)\) \(+\) \({\varLambda }_{{rad}}^{i}(T)\) where \(i\) indicates \({Ol}\) this study55,68, at shallow depth, \({Wd}\)36 between \(410\) and \(520\) \({km}\)12, and \({Rw}\)36,45 between \(520\) and \(660 \ {km}\)12. The presence of other upper mantle minerals, such as pyroxene or garnets, does not have a strong effect on the radiative thermal conductivity of the whole rock, because the optical properties of these phases are expected to be similar to that of olivine at high temperature35. Therefore, the inclusion of olivine and its polymorphs is adequate for estimating the aggregate radiative conductivity of slab lithologies. We also note that, in the computation of \({\varLambda }_{{tot}}\), we ignored the effect of grain size (i.e. light and phonon scattering) as they are unimportant to first order (SM). Overall, we ran a total of 40 models divided into two sets: 1) \({\varLambda }_{{rad}}=0\) and 2) \({\varLambda }_{{rad}} > 0\). In each set, we varied slab age \({t}_{{slab}}:\) [20; 40; 60; 80] \({Myrs}\) and \({v}_{{sink}}:\) [2.5; 5.0; 7.5; 10; 12.5] cm yr−1. The full description of the physical model and its numerical discretization, the governing equations, the limitations, and the benchmark of the numerical solution are reported in the SM.