Heat transport into the interior ocean induced by water-mass subduction

The subduction of oceanic water masses provides a crucial pathway for anthropogenic heat to enter the subsurface ocean, thereby shaping deep-reaching warming signatures. Analyzing data from eight ocean and atmosphere reanalysis datasets, we show that the average annual subduction rate of the global ocean (excluding 10° S–10° N) is 312.4 ± 27.9 Sv, resulting in a mean heat transport of 20.2 ± 2.1 PW towards the subsurface ocean. This subduction-driven heat transport has exhibited an increase of 0.09 ± 0.08 PW/decade since 1970. The increase predominantly stems from the overall enhancement of subduction within the latitudes of 30° S–50° S, dictated by intensified westerly winds that lead to the deepening of the local mixed layer depth. Our findings underscore the essence of wind-driven changes in the Southern Ocean subduction, which wield considerable influence over the global climate by regulating the vertical transport of heat and carbon from the sea surface to the deep waters.


Introduction
More than 90% of the excess energy in the climate system has been absorbed by oceans since the mid 20th century, leading to prevailing warming of seawater (e.g.Levitus et al 2012, Trenberth et al 2014, Roemmich et al 2015, Cheng et al 2019, Fox-Kemper et al 2021).While oceanic warming has accelerated over time (e.g.Bindoff et al 2019, Cheng et al 2019, Johnson and Lyman 2020, Loeb et al 2021), the global mean surface temperature exhibited a slowdown in warming during the 2000s (e.g.Meehl et al 2011, Trenberth andFasullo 2013).This phenomenon has been attributed to the enhanced heat absorption by the subsurface ocean (e.g.Chen and Tung 2014, 2018, Liu et al 2016, Xie 2016, Zanna et al 2019).Notably, the Southern Ocean, which stored 75% of the excess heat and 43% of anthropogenic carbon in the global oceans between 1861and 2005(Frölicher et al 2015)), has experienced substantial warming (e.g.Roemmich et al 2015, Sallée 2018, Swart et al 2018, Lyu et al 2020, Li et al 2022), and its mid-depth warming is nearly double the global average (Gille 2002, Cai et al 2010, Durack et al 2014).Climate models suggest that anthropogenic changes in seawater properties can emerge in the interior ocean and will continue to intensify in the global oceans under greenhouse warming (Silvy et al 2020).Therefore, exploring the exact channel serving for heat and anthropogenic carbon sequestration in the ocean is crucial for improving the projection of the global warming rate.
Previous studies have documented a connection between the increase in ocean heat content and the widespread distribution of mode water (Häkkinen et al 2015, Gao et al 2018).For example, Gao et al (2018) found that more than 60% of the excess heat in the extratropical Southern Hemisphere Ocean can be ascribed to the Subantarctic Mode Water.Recently, Li et al (2023) showed that the mode and intermediate water layers can account for 89% of the heat uptake over the global oceans during 2005-2020.Subduction, describing the irreversible transfer of water mass into the deep ocean, is the sole oceanic process that links surface forcing to subsurface water masses.It acts as a conduit for the transfer of heat, freshwater, carbon and nutrients from the sea surface to the ocean interior (Bates et al 2002, Sabine et al 2004, Palter et al 2005, Hartin et al 2011, Ludicone et al 2011, Sukigara et al 2011, Sallée et al 2012, DeVries et al 2017).Studies have indicated that the volume of Subantarctic Mode Water and Antarctic Intermediate Waters formation plays an important role in determining the efficiency of anthropogenic carbon and excess heat sequestration (Mignone et al 2006, Bourgeois et al 2022).Thus, it is inferred that water-mass subduction, primarily occurring in the mid-and high-latitudes, may be the key process for the excess heat into the interior ocean under global warming.
However, previous studies have mainly focused on subduction rates, water mass characteristics, and their variability (e.g.Liu and Huang 2012, Liu and Wu 2012, Xu et al 2016, Gao et al 2018, Qu et al 2020, Li et al 2021, Qiu et al 2021).Yet, the subduction-induced heat transport into the interior ocean and its potential impact on subsurface heat uptake remain poorly understood, limiting our comprehensive understanding of ocean heat redistribution under global warming.In this study, we use five ocean state estimates and three atmospheric reanalysis datasets to pursue it.The rest of the paper is organized as follows: section 2 outlines the datasets and methods utilized in our analysis.Section 3 presents the characteristics, long-term variability, and mechanisms of subduction-induced heat transport during 1970-2016.Section 4 provides a summary of our main findings.

Data and methods
The heat transport into the interior ocean induced by subduction, Q sub , is estimated as where C p and ρ represent the heat capability and density of seawater, respectively.T sub denotes the temperature of subducted water masses.Based on the 'Stommel demon' , the density and depth of mixed layer reach their local maxima in late winter, and a demon works that can select the late winter water mass properties to transfer into the interior ocean (Stommel 1979).Here T sub is identified as the average temperature of the late-winter mixed layer, typically March for the Northern Hemisphere and September for the Southern Hemisphere, respectively.The annual subduction rate, S ann , is defined as the volume of water irreversibly moving from the mixed layer and entering the pycnocline within one year (e.g.Woods 1985, Williams et al 1995).The widely used Lagrangian algorithm in previous studies was adopted to diagnose the annual subduction rate.In this method, the first subducted trajectory of water particles, released at the bottom of mixed layer in late winter, is tracked for one year (blue dashed arrow in figure S1).The difference in depth between this trajectory and the local mixed layer base is then assigned as the annual subduction rate (figure S1; Qiu and Huang 1995): ) . (2) Here, W ek = 1 ρ curl ( τ f ) is the Ekman pumping velocity, where τ denotes the surface wind stress.β is the meridional derivative of the Coriolis parameter f.The first term on the right-hand side indicates the averaged vertical pumping along the 1 year Lagrangian trajectory, while the second term represents the contribution from lateral induction.T is one year, h m,0 is the late-winter (March for Northern Hemisphere and September for Southern Hemisphere) mixed layer depth, and h m,1 indicates the mixed layer depth at the ending location of the 1 year Lagrangian trajectory.The mixed layer depth is estimated as the depth where the potential density is larger than that at the 10 m depth by 0.03 kg m −3 (de Boyer Montégut et al 2004).
The observational ocean temperature product from the Institute of Atmospheric Physics (IAP; Cheng et al 2017) and version 4 of the Met Office Hadley Centre 'ENSEMBLES' series of datasets (EN4; Good et al 2013) are used to identify the increase in ocean temperature during the past several decades.The monthly gridded temperature, salinity, and horizontal currents from the ocean state estimates were adopted to evaluate the subductioninduced heat transport into the interior ocean, and wind datasets were utilized to calculate the Ekman pumping velocity.Specifically, we analyzed five ocean state estimates (refer to table S1 for detailed information) and three wind datasets: the Japanese 55 year Reanalysis (JRA-55; Kobayashi et al 2015), the National Centre for Environmental Prediction (NCEP)/National Centre for Atmospheric Research Reanalysis 1 (Kalnay et al 1996), and European Centre for Medium-Range Weather Forecasts Reanalysis v5 (Hersbach et al 2019).Ensembles from 15 results were employed to examine the subductioninduced heat transport into the interior ocean.It is noted that we excluded the equatorial region within 10 • S-10 • N due to the absence of Ekman pumping velocity.

Subduction-induced heat transport and its change
The heat transport into the interior ocean induced by subduction, averaged over 1970-2016, is depicted in figure 1(a), with the maximum exceeding 300 W m −2 found in the tropics.This heat transport gradually decreases with the increasing latitude (figure 1(a)), different from the pattern of subduction rate (figure 1(b)) owing to the colder subducted waters at higher latitudes (figure 1(c)).For instance, although the Southern Ocean and northern North Atlantic are the key regions with high subduction rate, they are not the primary heat transport zones via subduction due to the cold surface water.
The climatological mean subduction rate, integrated over the global oceans (excluding the equatorial region), is estimated at 312.4 ± 27.9 Sv (1 Sv = 10 6 m 3 s −1 ; ensemble mean (EM) ± 1 standard deviation (STD) among datasets), resulting in a heat transport into the ocean interior of 20.2 ± 2.1 PW (1 PW = 10 15 W; EM ± 1 STD among datasets).The EM represents the average of 15 results calculated using five oceanic and three atmospheric datasets.The Southern Ocean (30 • S southward) is the largest contributor of subduction, accounting for 37.5% of the global subduction rate.However, it only caused a heat transport of 4.0 PW (19.8%).44.6% of the heat transport (9.0 PW) occurs in the Pacific Ocean (north of 30 • S), corresponding to a high subduction rate of 102.0 Sv.Moreover, the subduction rate in the Atlantic Ocean (north of 30 • S) is 71.3 Sv, leading to a heat transport of 5.1 PW.The Indian Ocean (north of 30 • S), contributing the least with 7% of the subduction rate, accounts for 10.4% of the heat transport, table 1.
Both the subduction rate and water temperature exhibit variability on multiple scales, including quasibiennial, interannual, decadal, and multi-decadal variations (e.g.  and then the EM is computed.To suppress interannual variability, a 5 year running mean filter was applied to the time series.The global subductioninduced heat transport shows an increasing trend of 0.09 ± 0.08 PW/decade (p < 0.05), depicted by the green line in figure 2(a).Overall, the subduction-induced heat transport into the interior ocean increased by 2.1 ± 1.9% over the period 1970-2016.It is worthwhile noting that ECDA v3.1 was based on the fully coupled climate model, while the other 4 oceanic datasets were based on different oceanic models with various forcing fields and resolutions (table S1).These differences lead to spread in the estimated heat transport.In addition, discrepancies in wind datasets are another error source.Adopting different wind datasets can also lead to different trends (figure S2).Generally, the trend in subduction-induced heat transport is accompanied by a considerable uncertainty across various ocean and wind datasets, with a STD of 0.25 PW/decade.Here, the average of various results derived from five ocean state estimates and three wind datasets, referred to as EM, is adapted to represent the change of subduction and its induced heat transport.
The changes in subduction-induced heat transport into the interior ocean, denoted as ∆Q sub , represent the difference between the decade-average Q sub during 2007-2016 and the average of 1970-1979 throughout the paper.It shows a large spatial variability, with an amplitude of approximately 100 W m −2 (figure 3(a)).Both increases and decreases in subduction-induced heat transport are observed globally, and this pattern is common in the results from the ocean state estimates despite differences in detail (figure S3).The most prominent increase in heat transport induced by subduction is located in the Southern Ocean (30-50 • S; figure 3(c)), with hotspots at higher latitude as they extend eastward in the Southern Indian and Pacific Oceans (figure 3 The change in subduction-induced heat transport into the interior ocean, ∆Q sub , can be expressed as: (3) Here, the first component represents the contribution of the annual subduction rate change ∆S ann , the second term is ∆Q sub caused by the subducted water temperature change ∆T sub , and the third term is the combined effect of ∆S ann and ∆T sub .As shown in figures 3(a)-(c), the contribution of ∆S ann shows a nearly identical spatial distribution to the subduction heat transport change.Particularly in the major zones of increased subduction heat transport (30-50 • S), the correlation coefficient (r) between subduction-induced heat transport and the subduction rate reaches 0.97.Under global warming, the contribution of water temperature change is positive, leading to an increase in heat transport into the interior ocean; however, its magnitude is much smaller (figure 3(d)).Similarly, the combined effects of ∆S ann and ∆T sub also make a minor contribution (figure 3(e)).Consequently, the change in subduction-induced heat transport into the interior ocean is dominated by the subduction rate change.Between 2007Between -2016Between and 1970Between -1979, the change in subduction rate is notably prominent in the Southern Ocean and North Atlantic, with the maximum exceeding 100 m yr −1 (figure 4(a)).Similarly, both increased and decreased subduction rates coexist in the global oceans, and this pattern is common across various ocean state estimates despite differences in detail (figure S4).However, there is no robust trend in the EM global subduction rate (figure 2(c)), which is different from the earlier estimate of 5 Sv/decade based solely on SODA2.2.4 (Liu and Huang 2012).Our estimate from SODA2.2.4 (9.7 ± 2.0 Sv/decade) also deviates somewhat from the estimate by Liu and Huang (2012), due to the exclusion of the equatorial region and usage of different vertical velocity.Nevertheless, the 5 year running means of subduction rate within 30-50 • S exhibited an upward trend of 1.66 ± 0.69 Sv/decade, resulting in an increase of 12.2 ± 5.1% from 1970 to 2016 (figure 2(d)).The increasing trend is consistent despite different rates in previous estimates (Liu and Huang 2012, Qu et al 2020), and it may be related to the differences in time periods, regions and datasets.
The following discussions will focus on the Southern Ocean within 30-50 • S due to its dominance in both the change of subduction rate and the induced heat transport into the interior ocean.
The variability of the subduction rate (black) is primarily controlled by the lateral induction (blue), showing a correlation of 0.96 for the global oceans and 0.99 for the Southern Ocean.In comparison, the variability of the vertical pumping is much smaller, and its correlation coefficient with the subduction  2016 and 1970-1979 also demonstrates that the lateral induction change is the major contributor to the subduction rate change, while the contribution from vertical pumping is minor (figure 4).The dominance of lateral induction can further be confirmed by the nearly identical spatial distribution of subduction rate change for a given ocean state estimate when utilizing various wind stress datasets, which are only used to calculate the Ekman pumping velocity (figure S4).As described in equation ( 2), the late-winter mixed layer depth is a crucial component of lateral induction, and previous studies have documented that it is the key factor inducing the changes in regionintegrated subduction rate (Liu and Huang 2012, Qu et al 2020).Indeed, there exists a linear relationship between the change in subduction rate and that in the late-winter mixed layer depth.The correlation coefficients are greater than 0.63 for 12 of the 15 results in the Southern Ocean (30-50 • S), all significant at a 99% confidence level.Taking the result from SODA2.2.4 and JRA-55 for example (figure 5(a)), r between the subduction rate changes and late-winter mixed layer depth changes reaches 0.81.Additionally, the regression coefficient, 0.6, means that if the latewinter mixed layer deepens by 1 m, the vertical water mass transfer into the interior ocean will increase by 0.6 m in one year.
The mixed layer depth change is associated with many physical processes, with the turbulent mixing of water mass caused by wind stress and air-sea buoyancy forcings being the major driver.As depicted in figure 5  increase of near-surface stratification in recent decades (Yamaguchi and Suga 2019).The buoyancy forcing is therefore not a driver of deepening mixed layer depth.Within 30 • -50 • S, the surface winds show a noticeable intensification of the westerly winds (figure 5(d)), which may be the major contributor to thicken the mixed layer (Liu and Huang 2012).Actually, the relationship between the mixed layer depth and its forcing for a specific ocean state estimate may better represent the role of wind stress in deepening the mixed layer.Take SODA2.2.4 for example, the 5 year running means of wind stress amplitude averaged over the Southern Ocean (30 • -50 • S) shows a strong correlation with the mixed layer depth, with a correlation coefficient of 0.77 (statistically significant at the 99% confidence level).

Conclusions and discussions
Here we present the subduction-induced heat transport into the interior ocean, which increases from the poles to the equator.Aggregated across five ocean state estimates and three atmospheric reanalysis datasets, the results show that subduction induces 20.2 ± 2.1 PW of heat transport into the global interior oceans (excluding the equatorial regions).This heat transport slightly increased from 1970 to 2016, with a growth rate of 0.09 ± 0.08 PW/decade.More importantly, the Southern Ocean within 30 • -50 • S exhibited a growth rate of 0.08 ± 0.03 PW/decade, accounting for approximately 89% of the global increase, due to the enhanced subduction rate driven by intensified westerly winds.
The excess downward heat transport induced by subduction can lead to the increase of subsurface ocean heat content.Since 1970, the accumulative excess heat transport into the subsurface ocean, ∆Q yr , can be simply estimated as where Q sub is the subduction heat transport and Q sub0 is Q sub at the year 1970.A growth rate of 0.09 ± 0.08 PW/decade can give rise to that the accumulative excess heat into the subsurface ocean driven by subduction reaches 306.8 ± 272.7 ZJ (1ZJ = 10 21 J) in 2016.Observation shows that the ocean heat content accumulation since 1960 reached 464 ZJ in 2023, with 60% stored in the interior ocean deeper than 300 m (Cheng et al 2024).Based on the same ocean warming rate, we can roughly deduce that there was an approximate 204 ZJ increase in the subsurface ocean heat content (>300 m) in 2016 compared to 1970.Thus, the enhancement of subduction-driven downward heat transport can fully account for the increase of subsurface ocean heat content.Moreover, the excess heat transport into the interior ocean induced by the Southern Ocean subduction exhibits a roughly similar geographic pattern as the ocean subsurface warming (figure S5), reflecting the potential role of subduction in altering the heat content in the interior ocean.However, they are not completely consistent because ocean circulation can redistribute heat away from the locations where it enters the ocean at the surface (Manabe et al 1991, Cai et al 2010, Kuhlbrodt and Gregory 2012, Bryan et al 2014, Exarchou et al 2015).
This work confirms that subduction serves as a primary channel for excess heat transport into the interior ocean under global warming.It operates as an active process whereby the oceanic heat storage varies significantly due to variations in ocean circulation by regulating vertical heat exchange between the surface and subsurface ocean.Our study is only the first step toward quantifying the excess heat transport into the interior ocean induced by subduction; there are several shortcomings.Firstly, the equatorial ocean (10 • S-10 • N) excluded due to the constraints related to Ekman pumping velocity.The result based on SODA2.2.4 has suggested that subduction within the equatorial ocean is non-negligible and has exhibited a decreased trend since 1970 (Liu and Huang 2012).Consequently, our result might overestimate the growth rate of subduction-induced heat transport due to the absence of equatorial ocean.Secondly, there exists a significant uncertainty stemming from the usage of different oceanic or atmospheric datasets.That is, the ensemble of various datasets may yield different growth rates for subduction and its induced heat transport.Specifically, the ocean state estimates used in our study are based on oceanic models with relatively low or eddy-permitting resolutions, thus incapable of capturing mesoscale eddies.Previous studies have documented that the mesoscale eddies exhibit water properties that differ from those of the surrounding water (e.g.Frenger et al 2015, Sun et al 2019, Wang et al 2021), and they can significantly enhance subduction by inducing high-frequency variability in mixed layer depth (Qu et al 2002, Uehara et al 2003, Rainville et al 2007, Nishikawa et al 2010, Kouketsu et al 2011, Xu et al 2014, 2016).These factors are likely to exert a significant impact on the subduction rate and the induced heat transport, including their magnitude, spatial distribution, and variability.Therefore, a more accurate estimate of subduction-induced heat transport and its variability, based on ensembles generated from a wide range of eddy-resolving datasets, is an urgent need to obtain a comprehensive understanding.

Figure 1 .
Figure 1.(a) The subduction-induced heat transport into the interior ocean (Q sub , unit: W m −2 ), (b) the annual subduction rate (Sann, unit: m yr −1 ), and (c) the subducted water temperature (T sub ) with the unit of • C averaged over 1970-2016.All the results are based on the average of 15 results calculated using five oceanic and three atmospheric datasets.The left panels show the meridional distribution of the zonal-mean values, with the standard deviation (STD) among datasets indicated by shadings.
Qu and Chen 2009, Liu and Huang 2012, Cheng et al 2019, Qu et al 2020, Auger et al 2021).Consequently, the anomaly in subductioninduced heat transport, Q ′ sub , displays considerable variability from year to year (figure 2(a)).This anomaly is calculated as follows: for each ocean state estimate and wind stress data, Q ′ sub is defined as Q sub minus the mean value during 1980-2000,

Table 1 .Figure 2 .
Figure 2. The EM subduction-induced heat transport anomaly (Q ′ sub ) of the global ocean (equatorial region excluded; (a)) and the Southern Ocean within 30-50 • S (b), with the shading denoting the ±1 STD among datasets.The green lines indicate the trend of the time series filtered by a 5 year running mean.The EM annual subduction rate anomaly (black) and its two components: lateral induction anomaly (blue) and vertical pumping anomaly (red) in the global oceans (c) and in the latitudinal band 30-50 • S (d), with the shading denoting the ±1 STD among datasets for the subduction rate anomaly.The climatology is defined as the average during 1980-2000.
Figure 3. (a) Decadal changes in the subduction-induced heat transport into the interior ocean between 2007-2016 and 1970-1979 (∆Q sub ).(b) ∆Q sub from the contribution of the subduction rate change ∆Sann.(d) ∆Q sub from the contribution of the subducted water temperature change ∆T sub , and (e) the combined effects of ∆Sann and ∆T sub .(c)/(f) panels: the meridional distribution of the zonal-mean values.Stipplings indicate significant changes at 90% confidence level of the Student's t-test.All results are from EM and in the units of W m −2 .

Figure 4 .
Figure 4.The EM decadal changes in (a) subduction rate and its two components: (b) lateral induction and (c) vertical pumping between 2007-2016 and 1970-1979.The left panels show the meridional distribution of zonal-mean values, with the shading denoting the ±1 STD among datasets.
(b), changes in late-winter mixed layer depth show a large spatial pattern, characterized by deepening in most regions of the Southern Ocean.Due to global warming (figure 5(c)) and surface freshening at high latitudes, changes in buoyancy forcings act to suppress turbulence, as suggested by the global

Figure 5 .
Figure 5.The relationship between the subduction rate changes versus late-winter (September) mixed layer depth changes derived from SODA2.2.4 and JRA-55 winds, computed with all data points in the latitudinal band 30-50 • S. The linear regression coefficient (K) and the correlation coefficient (r) are shown.(b) The EM September mixed layer depth change and its zonal-mean with the shading denoting the ±1 STD among datasets in the right panel.The EM changes in September sea surface temperature (∆SST; c) and wind stress ∆τ (arrows for wind stress vectors, color shading for wind stress magnitude (d), and the zonal mean values (black for ∆τ and blue for ∆SST), with the shading denoting the ±1 STD among datasets.

L
Liu et al Swart N C, Gille S T, Fyfe J C and Gillett N P 2018 Recent Southern Ocean warming and freshening driven by greenhouse gas emissions and ozone depletion Nat.Geosci.11 836-41 Trenberth K E and Fasullo J T 2013 An apparent hiatus in global warming?Earth's Future 1 19-32 Trenberth K E, Fasullo J T and Balmaseda M A 2014 Earth's energy imbalance J. Clim.27 3129-44 Uehara H, Suga T, Hanawa K and Shikama N 2003 A role of eddies in formation and transport of North Pacific subtropical mode water Geophys.Res.Lett.30 1705 Wang X, Zhang S, Lin X, Qiu B and Yu L 2021 Characteristics of 3-dimensional structure and heat budget of mesoscale eddies in the South Atlantic Ocean J. Geophys.Res.: Oceans 126 e2020JC016922 Williams R G, Marshall J C and Spall M A 1995 Does Stommel's mixed layer "demon" work?J. Phys.Oceanogr.25 3089-102 Woods J D 1985 The physics of thermocline ventilation Coupled Ocean-Atmosphere Models (Elsevier Oceanography Series) vol 40, ed J C J Nihoul (Elsevier) pp 543-90 Xie S P 2016 Oceanography: leading the hiatus research surge Nat.Clim.Change 6 345-6 Xu L X, Li P L, Xie S P, Liu Q Y, Liu C and Gao W D 2016 Observing mesoscale eddy effects on mode-water subduction and transport in the North Pacific Nat. Commun.7 10505 Xu L X, Xie S P, Mcclearn J L, Liu Q Y and Sasake H 2014 Mesoscale eddy effects on the subduction of North Pacific mode waters J. Geophys.Res.: Oceans 119 4867-86 Yamaguchi R and Suga T 2019 Trend and variability in global upper-ocean stratification since the 1960s J. Geophys.Res.: Oceans 124 8933-48 Zanna L, Khatiwala S, Gregory J M, Ison J and Heimbach P 2019 Global reconstruction of historical ocean heat storage and transport Proc.Natl Acad.Sci.USA 116 1126-31