Past and future ocean warming

Changes in ocean heat content (OHC) provide a measure of ocean warming, with impacts on the Earth system. This Review synthesizes estimates of past and future OHC changes using observations and models. The top 2,000 m of the global ocean has significantly warmed since the 1950s, gaining 351 ± 59.8 ZJ (1 ZJ = 1021 J) from 1958 to 2019. The rate of warming increased from <5 to ~10 ZJ yr−1 from the 1960s to the 2010s. Observed area-averaged warming is largest in the Atlantic Ocean and southern oceans at 1.42 ± 0.09 and 1.40 ± 0.09 × 109 J m−2, respectively, for the upper 2,000 m over 1958–2019. These observed patterns of heat gains are dominated by heat redistribution. Observationally constrained projections suggest that historic ocean warming is irreversible this century, with net warming dependent on the emission scenario. By 2100, projected warming in the top 2,000 m is 2–6 times that observed so far, ranging from 1,030 [839–1,228] ZJ for a low-emission scenario to 1,874 [1,637–2,109] ZJ for a high-emission scenario. The Pacific is projected to be the largest heat reservoir owing to its size, but area-averaged warming remains strongest in the Atlantic and southern oceans. Ocean warming has extensive impacts that pose risks to marine ecosystems and society. The projected changes necessitate a continuation and improvement of observations and models, along with better uncertainty estimation. Greenhouse gas forcing has increased ocean heat content, with large impacts on the Earth system. This Review outlines observed and projected global and regional changes, revealing an observed 0–2,000 m global increase of 351.4 ± 59.8 ZJ from 1958 to 2019, and a projected increase of 1,874 [1,637–2,109] ZJ by 2100 under SSP5-8.5.

Owing to the large thermal inertia of the ocean, subsurface warming represents the slow response to external influences, in particular to GHG forcing (Box 1). In response to past and current carbon emissions, future ocean warming is therefore committed for many centuries 5,9 , and is related to the current acceleration of ocean warming 10 (Fig. 1a). For instance, under representative concentration pathway (RCP) 8.5, it is projected that the total upper-2,000-m OHC increase from 2017-2100 will be ~5-7 times that observed from 1970 to 2017 (rEF. 6 ). The irreversibility of this ocean warming on centennial timescales creates additional requirements for climate policy, particularly considering the widespread impacts 10,11 .
Indeed, ocean warming has a multifaceted role in the Earth system via its links to the energy, water and carbon cycles, and resulting feedbacks (Fig. 1). For instance, the vast majority of radiative energy trapped by GHGs is ultimately absorbed at the ocean surface. Part of the energy is stored in the ocean, but this warming leads to increased longwave radiation and surface evaporation, cooling the ocean 12,13 (Fig. 1a, F2). Increased evaporation, in turn, leads to moistening of the atmosphere, released as latent heat during precipitation 14 , invigorating the hydrological cycle 15 (Fig. 1a, F3). In response, the sea surface salinity pattern is amplified, with salty regions getting saltier and fresh regions getting fresher 16 . These temperature and salinity changes also alter ocean density and circulation, triggering global ocean heat uptake (OHU) 6,17 . Ocean warming also affects the Earth's surface albedo [18][19][20] , and sea-ice and ice-sheet 21,22 melt processes, increasing OHU 11,23 (Fig. 1a, F4 and F5).
ocean heat redistribution processes, including circulation, mixing and convection, eventually spread this additional heat to the ocean interior [24][25][26] (Fig. 1a, F6). With this heat, anthropogenic CO 2 is also absorbed and transported, leading to a synchronous increase of OHC and ocean acidification 27,28 . However, through increased stratification and reduced uptake capacity and efficiency, warming decreases the efficiency of the oceanic carbon sink, further enhancing atmospheric CO 2 and ocean warming, a positive feedback 28,29 (Fig. 1a, F7).

Earth energy imbalance
EEi. The net downwelling radiation at the top of the atmosphere, represented as the balance of absorbed solar radiation (allowing for reflection and scattering) and outgoing longwave radiation.
Ocean heat content oHC, or ocean heat storage (oHs). A change or anomaly of the thermal energy of the ocean assumed to have a fixed volume (V x,y,z , in units of J), and vertical integration (V z , in units of J m −2 ). Calculated as oHC (x,y,z) = C TdV x y z ( , , ) p V x y z ( , , ) ∭ ρ following TEos-10 standards, where c p is a constant of ~3,991.9 J (kg K) −1 , ρ is potential density in kg m −3 and T is conservative temperature measured in degrees Celsius.

0123456789();:
In addition, ocean ecosystems and biological processes are also affected by ocean warming, and these then alter the carbon uptake and storage 6,30 . Thus, ocean warming threatens the fundamental conditions that make the Earth conducive to sustaining life (Fig. 1), necessitating better understanding of its changes.
In this Review, we outline past and future ocean warming, its drivers and consequences. We begin by outlining the current ocean observing system for monitoring ocean warming, followed by discussion of contemporary global and regional OHC changes. To support adaptation and mitigation, model projections of ocean warming are then provided. The far-reaching consequences of ocean warming on physical, human and biological systems of the Earth system are subsequently outlined, before ending with a discussion of the remaining challenges and outlook for monitoring and understanding ocean warming.
From 1999, however, the international Argo system of autonomous profiling floats revolutionized oceanographic observations, measuring the upper 2,000 m at unprecedented resolution, and dominating ocean observations since about 2005 (rEFs. 34,35 ). Currently, ~3,900 profiling floats constitute the array, providing a system to make near-global measurements of the Earth's open ocean 36 , with a target resolution of one profile every 3° × 3° every 10 days. However, the Argo data are much more limited in polar areas, shallow and coastal regions, in some major current systems, including the Indonesian throughflow (ITF), and in the water column below 2,000 m. Because the Argo network has incomplete global coverage, the global upper-2,000-m OHC rates based on Argo-only products can underestimate the rate of change by 10-20% [37][38][39] . Deep Argo has been developed to better capture measurements at 4,000-6,000 dbar (rEF. 40 ) with a resolution target of one profile every 5° × 5° every 10 days 41 . So far, ~170 Deep Argo floats have been regionally implemented out of the necessary 1,200. At a minimum, Argo floats measure temperature and salinity, but are increasingly also measuring biological and chemical properties 42 .
These Argo data are complemented by other observation platforms, including XBTs, MBTs, Nansen bottles, ship-based CTDs (conductivity-temperature-depth sensors), gliders and moored arrays, offering opportunities for OHC comparisons and data acquisition in regions where Argo does not operate. For instance, the global XBT network still provides continuous and high-resolution (eddy-resolving) temperature profile data along repeated transects and critical channels 43 . The Marine Mammals Exploring the Ocean Pole to Pole (MEOP) programme -which equips marine animals, primarily southern elephant seals or other pinnipeds, with measurement devices -has also provided 600,000 vertical profiles in high-latitude coastal regions since about 2004 (rEF. 44 ). Ice-tethered profilers offer further data in high latitudes, obtaining temperature and salinity profiles with excellent vertical resolution.
Synthesizing all these in situ measurements into estimates of OHC requires consideration of quality control, bias correction, gap-filling (mapping) and uncertainty 33,45 . For example, the first estimate of OHC indicated long-term ocean warming in the upper 3,000 m over 1948-1998, and a particularly warm period from the 1970s to the early 1980s 46,47 (Supplementary  Table 1). However, this pronounced warm period was found to be a spurious consequence of systematic errors in XBT data 48 , as similarly identified in MBT 49 and Argo data 50 . These data quality issues and error corrections must be taken into account (Supplementary Table 1). The mapping method -how a global map is created from incomplete and inhomogeneous observationsalso introduces uncertainty 36,51,52 . Many mapping strategies based on an objective analysis approach result in a conservative bias toward low-magnitude changes by assuming no change in less-sampled regions [53][54][55][56][57] (Supplementary Table 1). In doing so, some OHC estimates exhibit unduly weak long-term trends 5,53,58 . Several methods have been proposed to resolve these mapping challenges arising from data scarcity, including: the use of model simulations to guide the reconstruction 59 ; incorporating SST 60 and/or sea level information 61 ; combining five years of data for a pentadal estimate 62 ; and using machine learning 63 (Supplementary Table 1).
In addition to these direct methods to observe and estimate OHC, several alternative approaches are also available, each with their own limitations (Supplementary Table 1). Indirect methods can be used to infer OHC change from physical relationships between OHC and other variables 33,64,65 . For instance, altimetry and space gravimetry have been used since the early 2000s to estimate OHC based on a sea level budget constraint 64,[66][67][68] , and SST used based on an assumption that temperature Ocean heat uptake oHU. The accumulated contribution of heat added into the ocean (heat gain) or removed from the ocean (heat loss) through heat fluxes at the air-sea, ice-sea and land-sea interfaces (in units of W m −2 ). globally, it is synonymous with 'oHC change, trend, rate, tendency' . regionally, oHU and redistribution contribute to local oHC change.

Ocean heat redistribution
The transport of heat within the ocean without involving any net global ocean warming or cooling through advection, convection, eddy mixing and small-scale diffusion.

NATure revIeWS | EARTH & EnvIROnmEnT
anomalies behave as a passive tracer 69,70 . These indirect estimates rely on simplified assumptions, introduce additional uncertainty (Supplementary Table 1) and exhibit larger spread (Fig. 2), but can complement direct estimates.
Reanalyses can further be used to infer OHC change, wherein data assimilation of various observations constrains a model 71,72 . These ocean reanalyses produce 4D (latitude, longitude, depth and time) gridded estimates of all variables, making them very suitable for mechanistic research. However, most reanalyses do not conserve properties (for example heat and salt budget) during the assimilation step and so are not dynamically consistent 73,74 . Comparing these reanalysis datasets indicates substantial spread in OHC changes, implying that many are hampered by model and assimilation biases 75 . Yet reanalysis has the potential to be superior to objective analysis after ~2005 (with better consistency with top-of-atmosphere net radiation observations) after properly dealing with model and assimilation biases 76,77 , but with large and spurious variability before this time 75 (Supplementary Fig. 1). Given that long-term changes are the focus of this Review, ocean reanalyses are therefore not included. Instead, subsequent discussion and analysis of OHC will focus on direct and indirect datasets that have complete bias correction and a valid evaluation procedure (Supplementary Table 1 (TABlE 1), albeit with some variability in trends and uncertainty owing to differences in data selection, trend estimates and methodological approaches 6,38,82 (Supplementary Information). These warming rates have increased over time (Fig. 2b,d). For instance, OHC at 0-700 m increased from ~0-3.5 ZJ yr −1 in the 1960s to ~4.8-8.0 ZJ yr −1 by the 2010s. OHC at 0-2,000 m increased from ~0-4.8 ZJ yr −1 to ~8.0-11.5 ZJ yr −1 over the same time periods. This quantification of acceleration is valid for >15-year periods, but its detection can be hampered by short-term fluctuations in OHC. At interannual timescales, global OHC fluctuations are dominated by El Niño-Southern Oscillation (ENSO) events 79,80,83 . During and after El Niño events, OHC decreases as some heat makes its way into the atmosphere, contributing to a temporal global surface warming. Major volcanic eruptions have also caused interannual to decadal cooling phases (Agung in 1963, El Chichón in 1982 and Pinatubo in 1991) 84,85 . Nevertheless, a shorter time series of EEI also reported an increase in rates since the early 2000s 7,66,86 , and the interplay between long-term change and interannual natural variations still needs to be elucidated.
A robust long-term acceleration is also evident in CMIP5 and CMIP6 simulations 85,87,88 , highlighting the anthropogenic connection to this accelerated warming (Fig. 2b,d). Indeed, detection and attribution research indicates that observed multidecadal-scale changes in global OHC are primarily the result of anthropogenic GHG emissions (partly offset by anthropogenic aerosols and volcanic forcings) 85,[89][90][91] , and that they have already emerged from background natural internal climate variability 92,93 . The OHC rate increase since 2001 is also thought to be driven by a large decrease in reflected solar radiation with a small increase in emitted infrared radiation associated with anthropogenic forcing 94 .

Global deep and full-depth changes.
Ocean warming is also evident at depth (Fig. 2e,f and TABlE 1). For instance, from 1958 to 2019, global OHC gain below 2,000 m is estimated at 26.0 ± 16.6 ZJ, or 0.43 ± 0.27 ZJ yr −1 . This estimate assumes that anthropogenic influence is not detectable before 1991 owing to the long response time of the deep ocean 38,58,59 , but such deep layers might be affected by past climates that could offset a portion of the anthropogenic signal 95 . This post-1991 uptake assumption has an observational basis. For example, highquality hydrographic survey data indicate notable ocean warming of 0.06 ± 0.03 W m −2 (0.96 ± 0.48 ZJ yr −1 ) since 1992 (rEFs. 3,96 ). These warming signals are confirmed by other direct and indirect estimates 97 , albeit with variability in magnitude, including from machine-learning 63 and variational minimization 98 approaches (Fig. 2e,f). Before the 1990s, observational estimates of OHC changes are more inconsistent 63,69,70 (Fig. 2e). However, the CMIP6 multimodel mean suggests a continuous increase in the deep ocean warming rate since the 1950s (Fig. 2f), but with substantial uncertainty due to diffusivity parameters, feedbacks and aerosols 99 .
OHC estimates at 0-2,000 m and below-2,000 m can be combined to yield total OHC changes at full depth. In doing so, it is estimated that ocean warming was 378.4 ± 64.5 ZJ from 1958 to 2019 (6.20 ± 0.98 ZJ yr −1 )

Box 1 | Surface and subsurface ocean warming
The land, atmosphere and upper ocean (surface and mixed layer) respond relatively quickly (normally in years 206 ) to surface radiative forcing, whereas the deep ocean typically adjusts over centuries to millennia 11 . Surface warming is approximately determined by cumulative emissions up to a given point in time, owing to a nearcancellation between the positive climate commitment and negative carbon cycle commitment 9,270 . This relationship implies that zero Co 2 emissions (assuming constant non-Co 2 emissions) lead to near-constant global surface temperature 270 -a finding that has critical implications for climate-change mitigation. However, the subsurface ocean responds far more slowly to anthropogenic forcing, in part owing to the slow nature of large-scale ocean circulation. Future subsurface ocean warming therefore has a much longer commitment time 9,11 . even after net zero emissions are reached, the ocean subsurface will continue to warm as heat is transported downwards into deeper ocean waters, and a positive earth energy imbalance remains until ocean and land carbon uptake sufficiently reduces atmospheric Co 2 concentrations. ocean heat content (oHC) and sea level rise both integrate ocean temperature changes below the surface, and they continue to increase long after surface temperature stabilizes. This fundamental distinction between surface warming and oHC increases has important implications for both mitigation and adaptation strategies. Furthermore, as the measurements of sea surface temperatures (SSTs) rely on different instruments, the observations of SSTs are fundamentally different in their temporal and spatial coverage and temporal extent compared to subsurface observations and oHC estimates 271 Table 3); shading represents the model spread (±1σ). Preindustrial control runs are used to adjust for simulation drift (Supplementary Information). c | As in a but for 0-2,000 m (rEFs. 37,52,55,56,59,60,62,63,65,67,69,269 ). d | As in b but for 0-2,000 m, and also including rates of change from CMIP5 (dashed lines). e | As in a but for below 2,000 m (rEFs. 63,69,70,[96][97][98] ). f | As in b but for below 2,000 m. The ocean has been warming at an increasing rate since the 1950s, with warming irreversible this century.  Table 2). Thus, as a global average, OHC has increased dramatically at the surface, at depth, and over the entire water column. These changes are particularly evident after 1990 and are accelerating over time. However, these global estimates mask regional contrasts between ocean basins.
Observed regional patterns of change At the regional scale, there is distinct variability in surface and deeper OHC changes, with stronger area-averaged warming in the low-and middle-latitude Atlantic Ocean and southern oceans compared with other basins (Fig. 3).
Here, OHC changes are assessed in six regions: the Pacific Ocean (35° S to 60° N), Atlantic Ocean (35° S to 64° N), Indian Ocean (35° S to 30° N), southern oceans (78°-35° S; the broad region south of Africa), the Mediterranean Sea and the Arctic Ocean. Given that the Cheng et al. 59 dataset minimizes regional errors through both mapping and instrumental bias correction [100][101][102] , it is used to calculate these regional OHC changes for the upper 2,000 m (see Supplementary Fig. 1 for other datasets). Quantification of deep OHC below 2,000 m is limited by data availability and is thus not thoroughly discussed.
Pacific Ocean. Changes in area-averaged Pacific Ocean OHC tend to resemble those of the global mean. Total OHC increased by 71.6 ± 8.2 ZJ from 1958 to 2019 for the upper 2,000 m (~20% of the global upper-2,000-m OHC increase), with an area-averaged change of 0.49 ± 0.06 × 10 9 J m −2 (Fig. 3c). This area-averaged change is comparatively smaller than other basins ( Fig. 3a) owing to the large Pacific Ocean volume, and the less effective downward transfer of heat arising from an absence of deep water formation north of 35° S and limited mode water and intermediate water formation in the North Pacific 103 .
In addition to these long-term area-averaged changes, there is pronounced spatial and temporal variability. Prior to 1990, a statistically insignificant trend of 0.05 ± 0.21 ZJ yr −1 is observed, followed by a substantial OHC increase of 2.42 ± 0.29 ZJ yr −1 since 1991 (Fig. 3c). Moreover, the Pacific Ocean warming pattern is characterized by weak and statistically insignificant trends in the western basin, and weak but robust warming (~0.25-0.5 W m −2 ) in the east (Fig. 3a). Northeast Pacific regional warming is also significant and characterized by prolonged marine heat waves (MHWs) from 2014 to 2021 (rEFs. 104,105 ).
This spatial and temporal variability emerges due to the presence of strong modes of climate variability, particularly ENSO 79,106-108 , the Interdecadal Pacific Oscillation (IPO) or Pacific Decadal Oscillation [109][110][111] , and Atlantic Multidecadal Variability (AMV) 112,113 . For instance, in the tropical Pacific, OHC decreases by ~0.1-0.3 W m −2 K −1 during the transition from El Niño to La Niña. This change is linked to ocean heat release associated with anomalous warming of surface waters 79 and subsequent evaporative cooling, in turn causing warmer global-mean surface temperatures during El Niño years 114 . Besides, ENSO also causes regional OHC variability from month to month associated with an adiabatic redistribution of heat both laterally and vertically in the tropical ocean, resulting in a local monthly OHC change up to 100 W m −2 K −1 and an opposite OHC variation between 0-100-m and 100-300-m layers 79 (Fig. 3c).
Although interannual Pacific OHC fluctuations are primarily influenced by ENSO, decadal-scale changes are dominated by the IPO. During a negative phase, stronger trade winds and an enhanced shallow tropical  81,259 , and then averaged. Trends represent averages over the Earth's surface (5.1 × 10 14 m 2 ). b 90% confidence intervals are given, whereby uncertainty is quantified by estimating internal error and structural error separately, and assuming they are independent (see Supplementary  Information). If instead it is assumed they are not independent, simply summing the two errors gives roughly 1.4 times larger error bars, as an upper limit.

Mode water
Formed when winter mixed layers are convectively overturned owing to gravitational instability, and characterized by low potential vorticity and nearly vertically homogeneous temperature, salinity and density.
overturning increase OHU and subsurface heat storage 115 . This greater OHC increase is associated with a slowdown in global-mean surface air temperature increases 115,116 .
The IPO, along with ENSO, also influences the ITF and thus heat redistribution between the Pacific and Indian Oceans 107,117 . Again during a negative IPO phase (or a La Niña event), model and observational evidence indicates that a substantial amount of heat can be transported from the Pacific to the Indian Ocean 107,118-120 . In the northwest Pacific, however, AMV governs decadal variability in North Pacific subtropical mode water and thereby OHC change 112,121 . For instance, during a positive AMV, as observed since ~1980, warming anomalies in the North Atlantic Ocean led to a poleward expansion of the westerly wind belt 114 and the subtropical ocean circulation 122,123 , and a resulting large OHC anomaly in the northwest Pacific Ocean 112 . As a result of this decadal variability, OHC trends calculated over these timescales are influenced by the start and end points of the analysis, and are therefore somewhat uncertain.
Atlantic Ocean. OHC changes in the Atlantic are larger than those observed in the Pacific (Fig. 3a,d). Indeed, in the upper 2,000 m, OHC increased by 117.2 ± 7.5 ZJ over 1958-2019 (~33% of the global upper-2,000-m OHC increase), with an area-averaged warming of 1.42 ± 0.09 × 10 9 J m −2 , ~3 times that for the Pacific Ocean (Fig. 3c). Spatial and temporal variability is also evident. For example, the Atlantic OHC trend doubled from 1.29 ± 0.19 ZJ yr −1 in 1958-1990 to 2.66 ± 0.27 ZJ yr −1 in 1991-2019, indicating an earlier onset of warming in the Atlantic than the Pacific. The Atlantic is further host to some of the largest rates of regional OHC changes, >1 W m −2 in most regions between about 35° S and 40° N, but exceeding 2 W m −2 in the Gulf Stream of the North Atlantic (Fig. 3a). Yet, whereas warming dominates the Atlantic average and regional signatures, cooling and a resultant reduction in heat content is evident at subpolar latitudes of the North Atlantic (~45-70° N, ~70° W to 0°) down to >800 m depth 124 (Fig. 3a). The mean rate of OHC change is −0.6 W m −2 in this region since 1958. This 'cold blob' and the apparent cooling/warming dipole are a predicted consequence of circulation changes and heat redistribution in response to GHG and aerosol forcing 4,24,27,90,125,126 . Notably, these circulation responses include a slowdown of Atlantic Meridional Overturning Circulation (AMOC) 127,128 , the fingerprint of which closely matches that of the OHC pattern 129,130 . However, there is a lack of long-term observational evidence for this AMOC slowdown, with some direct and indirect estimates suggesting a subtle reduction in meridional heat transport at 26° N in the Atlantic 128,131-133 , in concert with fluctuations in the North Atlantic Oscillation 134 (NAO). Moreover, modelling results are inconsistent in terms of the spatial pattern and timing of the cold blob in response to AMOC changes 135,136 , partly because of the incoherence of AMOC at different latitudes 137 . As such, changes in the AMOC and connections to the cold blob are still difficult to reconcile.
Atmospheric and air-sea interactions are vital processes in determining these North Atlantic OHC changes 138 . There have been substantial changes in the atmospheric circulation over the subpolar North Atlantic, including a northward movement of the jet stream and increased storminess 139,140 . Moreover, in the subarctic region where the cold blob is located, there is substantial heat loss attributable to intensified local winds associated with an NAO-like atmospheric circulation pattern during those winters 141,142 . A positive NAO favours strong westerly winds across the Atlantic mid-latitudes, bringing cooler drier air from North America and leading to increased winter surface heat and moisture fluxes into the atmosphere 143 . This elevated surface flux cools the ocean, increases deep water formation in the Labrador Sea area, and leads to a stronger AMOC several years later 143,144 . North Atlantic subpolar cooling has also been traced to the remote impact of Indian Ocean via atmospheric teleconnections 145 .
Indian Ocean. In the Indian Ocean, OHC changes have been more subdued until the late 1990s. Here, total upper-2,000-m warming is 32.1 ± 4.0 ZJ from 1958 to 2019 (accounting for ~9% of the global upper-2,000-m OHC increase), and area-averaged warming 0.70 ± 0.08 × 10 9 J m −2 . Although the total Indian OHC increase is smaller than other basins because of its area, this area-averaged warming is slightly larger than the Pacific Ocean but much smaller than the Atlantic Ocean. Observational estimates suggest a statistically insignificant trend of −0.05 ± 0.10 ZJ yr −1 from 1958 to 1990 (Fig. 3e), consistent across different datasets 146 ( Supplementary Fig. 1). However, from the late 1990s, rapid warming is evident, reaching 1.21 ± 0.14 ZJ yr −1 from 1991 to 2019, about half of the Pacific warming during the same period. The weak warming trend before 2000 has been attributed to interdecadal changes in the Indo-Pacific Walker circulation and a corresponding upper-layer cooling of the tropical Indian Ocean, moderating GHG-induced warming 147 . By contrast, the unprecedented rapid warming since 2000 has been linked to increased wind-driven heat transport in the ITF during the negative phase of the IPO 120,148 . Local wind and heat flux forcing also had a substantial contribution 146 .
In addition to these decadal-scale fluctuations, there is also substantial year-to-year variability. For example, 12 ZJ of heat was lost from 2015 to 2018, with a similar magnitude of recovery from 2019 to 2020 (Fig. 3e). These pronounced interannual fluctuations in Indian Ocean heat content are partly controlled by ENSO via its atmospheric teleconnection to the Indian Ocean basin mode 149,150 and impacts on the ITF 106,117,151 . During the 2015/2016 super El Niño event, for example, an unprecedented reduction of ITF volume and heat transport was responsible for the sharp reduction of OHC in the Indian Ocean 107,119 . The Indian Ocean is also modulated by Southern Ocean variability and climate trends, including the observed poleward migration of westerly winds, a southward shift of the subtropical gyre and resulting multidecadal trends in regional OHC 152-154 over the subtropical southern Indian Ocean. 1958 to 2019 (accounting for ~36% of the global upper-2,000-m OHC increase), equating to 1.40 ± 0.09 × 10 9 J m −2 (Fig. 3f). The rate of southern oceans warming has almost doubled from 1.49 ± 0.23 ZJ yr −1 over 1958-1990 to 2.72 ± 0.29 ZJ yr −1 over 1991-2019, and is consistent across datasets (Supplementary Fig. 1). Hence the NATure revIeWS | EARTH & EnvIROnmEnT southern oceans act as the largest heat reservoir since the 1950s, holding twice the excess anthropogenic heat as the Pacific. This OHC trend is broadly consistent with that of the Atlantic, which collectively reveal more intensive local warming than other basins. In the southern oceans, these strongest regions of warming are concentrated on the northern flank of the main fronts of the Antarctic Circumpolar Current (ACC) (Fig. 3a) (note also distinct surface cooling north of the Ross Sea 158,159 ). These long-term OHC changes are primarily attributable to greenhouse gas forcing, with a secondary role of stratospheric ozone depletion 160,161 , via several processes. For example, wind-driven upwelling of unmodified deep water keeps surface ocean cold, and the cold water absorbs anthropogenic heat which is then exported to the northern flank of the ACC by the background overturning circulation. There, it is subducted in the deep mixed layers formed north of the ACC 162,163 . These deep mixed layers are associated with the formation of subantarctic mode water 164,165 , and, in turn, Antarctic Intermediate Water -it is these waters that are absorbing significant quantities of anthropogenic heat 166 . However, model limitations along with sparse measurements have hampered understanding of the southern oceans' changing hydrography [167][168][169] . For example, theoretical predictions and high-resolution ocean modelling suggest that the observed strengthening of the westerly winds would drive an increase in eddy activity without changing the mean ACC transport, because the energy cascades rapidly into small-scale eddy fields [170][171][172] . This effect should, in turn, impact ocean interior mixing 173 and ocean heat content trends, but modern-day climate models do not generally resolve mesoscale eddies over the southern oceans.
Other processes are also thought to be important in driving OHC changes to the north of the ACC. For instance, lateral currents driven by buoyancy and steered by topography are associated with an observed acceleration of the zonal geostrophic currents within 45° S to 60° S. The change of currents leads to anomalous heat transport that relates to the pattern of OHC 174 . Moreover, freshening, linked to an amplification of the global hydrological cycle and sea-ice and ice-shelf changes, has weakened entrainment of subsurface warm water and could regulate warming patterns in the southern oceans 175,176 . Furthermore, the change in Antarctic sea-ice extent (a decrease since late 2016 after decades of increase) is probably associated with atmospheric circulation variations (the southern annular mode), and then could regulate the surface heat flux and OHC 13,177,178 .
In addition, the Tasman Sea has also witnessed substantial OHC increases 179,180 , particularly near southeastern Australia (Fig. 3a). Changes here are potentially linked to a southward migration of the East Australia Current extension 119 and a reduction in the ITF, linked to heat redistribution by variability in ocean circulation 125 .
Abyssal ocean warming is consistent with a reduced volume of Antarctic Bottom Water entering the deepest ocean layers. This abyssal warming has been linked to a meltwater-induced slowdown of the lower cell of the global meridional overturning circulation 181 , which would be consistent with abyssal OHC changes being largely due to heat redistribution, and not via anthropogenic heat uptake within the lower overturning cell.

The Mediterranean Sea.
In the Mediterranean Sea, upper-2,000-m OHC increased by 3.3 ± 0.4 ZJ from 1958 to 2019, with area-averaged warming of 1.11 ± 0.12 × 10 9 J m −2 . This OHC increase exceeds the global, Pacific, Indian and Arctic Ocean changes (Fig. 3g). As in these other basins, warming trend is superimposed by multidecadal variability 182,183 , with OHC rate increasing from −0.01 ± 0.01 ZJ yr −1 for the period of 1958-1990 to 0.12 ± 0.01 ZJ yr −1 after 1991.
The Mediterranean Sea is linked with the adjacent North Atlantic Ocean via their overturning circulation systems 184,185 , which exert control over the basin's ocean heat budget in addition to the local surface heat flux 186,187 .
In the upper ocean, the Atlantic water (0-150 m) flows into the Mediterranean Sea and moves eastward, so the climatic conditions along this path are crucial in determining the heat content. A stronger OHC increase in the eastern basin (where warmer intermediate waters formed) than the western basin is linked to the anomalous atmospheric conditions and reduced river freshwater input 124,188 . These warmer (and also saltier) intermediate waters then spread towards the western basin at sea subsurface (150-450 m) on their way back to the North Atlantic 189,190 . The increasing transport of heat from the eastern to the western Mediterranean affects the deep water formation process in the Gulf of Lion, enhancing the tendency of this site to produce warmer and saltier deep waters 188,191 . When the heat of the intermediate and deep waters flows out of the Strait of Gibraltar (Mediterranean Outflow Water), it, in turn, influences the meridional overturning circulation in the Atlantic Ocean 188,192,193 .
Arctic Ocean. The Arctic Ocean has been steadily warming. Total upper-2,000-m OHC increased by 8.6 ± 1.3 ZJ from 1958 to 2019, equating to area-averaged changes of (0.59 ± 0.09) × 10 9 J m −2 . Accordingly, about 2.5% of excess heat in the global upper-2,000-m ocean is stored in the Arctic basin since 1958 (3.3% since 1991) 124,194 , much of it being in the Atlantic Water layer 20,194 . Most of this OHC increase has occurred since the 1990s 195 (Fig. 3h), with statistically significant warming rates of 0.29 ± 0.05 ZJ yr −1 since 1991 compared with Fig. 3 | Observed and projected regional OHC changes. a | Observed OHC trends from 1958 to 2019 (rEF. 59 ), with non-stippled regions depicting regions with statistically significant trends at the 90% confidence level. b | Total (left axis) and per-area (right axis) annual OHC changes from observations 59 (black), and CMIP6 ensemble means under shared socioeconomic pathways SSP1-2.6 (blue) and SSP5-8.5 (red), all relative to a 2005-2019 baseline. Constrained and unconstrained projections for 2081-2100 are provided in bars in b, and in all other places unconstrained results are presented. Shading depicts model spread of ±1σ. Red dashed lines illustrate the linear slope of the heating rate (for global, the Earth's surface is used; for each basin, its area is used). The lower panel depicts temperature anomalies as a function of depth and time from observations (1958-2019) and SSP1-2.6 (2020-2100), relative to a 2005-2019 baseline (see Supplementary  Fig. 8 for SSP5-8.5 equivalents). c | As in b, but for the Pacific Ocean. d | As in b, but for the Atlantic Ocean. e | As in b, but for the Indian Ocean. f | As in b, but for the southern oceans. g | As in b, but for the Mediterranean Sea. h | As in b, but for the Arctic Ocean. Although present in all basins, ocean warming is unevenly distributed, with future warming dependent on socioeconomic scenarios. ◀ www.nature.com/natrevearthenviron 0.01 ± 0.03 ZJ yr −1 from 1958 to 1990. While statistically significant, uncertainty of the Arctic OHC time series is much larger than for the other five basins (Fig. 3h) owing to reduced subsurface data availability in sea-ice conditions. Nevertheless, observational datasets and a reanalysis product 194 exhibit consistent long-term changes for Arctic OHC (Supplementary Fig. 1).
OHC changes in the Arctic Ocean are strongly affected by rapid environmental changes, ocean heat transport and local air-sea/ice-sea heat exchanges 20,194 . For instance, ocean warming can be directly affected by heat transport into the region associated with the rapid decline of sea ice and the concurrent expansion of Atlantic waters into the Arctic Ocean (' Atlantification') [196][197][198][199] . However, direct observations reveal no significant trend in heat transport from the subtropical North Atlantic to the Arctic from 2014 to 2018 119,200 , although others suggest an increase from 1993 to 2016 201 , highlighting disparities related to both data uncertainty and substantial interannual variability 20 . Other environmental changes also affect Arctic Ocean OHC, including amplified surface warming since the 1960s at a rate twice the global value, driven by positive surface heat flux into the ocean associated with reduced outgoing longwave radiation and a reduction in albedo from the sea-ice-albedo feedback 23,197 . However, the air-sea/ice-sea heat exchanges are currently poorly quantified because of the data insufficiency 20 .
Thus, over 1958-2019, the southern oceans and Atlantic Ocean are the main heat reservoir, accounting for 36% and 33% of the global OHS, respectively; the Pacific and Indian oceans account for the other 20% and 9%. However, the percentages change with time and ocean area, given strong decadal variability of OHC change in each basin. Area-averaged warming is larger in the southern oceans, Atlantic Ocean and Mediterranean Sea (~1.11-1.42 × 10 9 J m −2 ) compared with the other basins (~0.49-0.70 × 10 9 J m −2 ), indicating very intensive warming mainly associated with the ocean circulations.

Ocean warming projections
In addition to the observed OHC changes, past GHG emissions have also committed the global ocean to future warming, the magnitude of which is dependent on the future socioeconomic pathway 6,202 (Fig. 3). Using drift-corrected models from the CMIP6 archive (Supplementary Table 3; Supplementary Fig. 2), these projected OHC changes are now discussed for the globally averaged ocean and for individual basins.  (Fig. 3b). However, these projections are strongly influenced by individual models' equilibrium climate sensitivity (ECS). Specifically, some CMIP6 models exhibit ECS > 5 K (rEF. 92 ) (Supplementary Figs. 3-7), and those models with higher ECS tend to project stronger ocean warming by the end of the century owing to strong correlations between those metrics 203,204 . As such, there is strong motivation to use observational trends from 2005-2019 to better constrain future projections 205 and produce more realistic OHC estimates (Supplementary Information).
In doing so, global OHC trends are slightly enhanced and uncertainties reduced, with the weight of both very low and very high ECS models being reduced. Total observationally constrained upper-2,000-m warming by the end of the century is estimated at 1,030 [839-1,228] ZJ and 1,874 [1,637-2,109] ZJ for SSP1-2.6 and SSP5-8.5, respectively (Fig. 3b). These projections are broadly consistent with CMIP5 results 5,205 (Fig. 2d), and represent OHC changes ~2-4 and ~4-6 times the observed 1958-2019 value (  Table 4 and 5). Importantly, these observationally constrained projections result in a substantial reduction in the uncertainty range, especially the upper bound (Supplementary Table 4 and 5). Thus, the largest CMIP6 projections from both very low and very high ECS models are very unlikely as they overestimate ocean warming over 2005-2019.
In addition to differences in total OHC between scenarios, rates of warming also vary. For instance, the projections (here unconstrained) indicate that, in the upper 2,000 m, the ocean warming rate will probably peak between about 2030-2040 at ~0.8 W m −2 and then decrease to ~0.5 W m −2 at the end of this century for low-emission scenarios (SSP1-2.6 and RCP2.6; Fig. 2d). In contrast, for a high-emission future (SSP5-8.5 and RCP8.5), ocean warming continues to increase throughout the twenty-first century to ~2.0-2.2 W m −2 by the 2090s, ~4 times the current level of OHU (Fig. 2d). Such scenario differences are also evident in the upper 700 m; OHC rates peak at ~0.5-0.6 W m −2 around 2030 and decrease to ~0.2-0.3 W m −2 during the 2090 s for SSP1-2.6, and continue increase to ~1.5 W m −2 during the 2090s for SSP5-8.5 (Fig. 2b). At depth (<2,000 m), these scenario differences are much smaller (Fig. 2f).
The vertical structure of global ocean warming is further distinct under different scenarios. For SSP1-2.6, global sea surface warming is stabilized at ~0.6 °C (relative to a 2005-2019 baseline) after ~2050 (Fig. 3b), consistent with the Paris Agreement of limiting global surface temperature to 2 °C above the preindustrial level. However, subsurface (below ~100 m) warming continues, and anthropogenic heat penetrates the deep layers in all ocean basins (Fig. 3b). For SSP5-8.5, both global surface and subsurface warming continues to the end of this century and beyond, with a surface warming of >2.8 °C and upper-600-m warming of >1°C above the 2005-2019 average (Supplementary Fig. 8).
Transient climate response and the vertical structure of OHC. The considerable heat penetration into the ocean interior critically affects surface temperature and therefore the global climate. The rate of OHC increase and

Shared socioeconomic pathway
The ssPs are a set of plausible trajectories of societal development and radiative forcing. ssP1-2.6 is a relatively low-emission scenario, representing the pathways to limit the global surface warming below 2 °C, which requires immediate, rapid and large-scale reductions in greenhouse gas emissions. ssP5-8.5 is a higher emissions scenario with projected warming >3 °C by 2100.

NATure revIeWS | EARTH & EnvIROnmEnT
the efficiency with which heat is transported from the upper ocean into the deeper ocean are critical in setting the strength and timescale of the climate adjustment in response to anthropogenic influences 24,203,206 . One important climate-change metric is the transient climate response (TCR), the surface warming under a 1% per year CO 2 increase from preindustrial values at the time of CO 2 doubling (70 years). The rate and patterns of ocean heat content are linked to TCR via the coupling at the sea surface 203,207 .
Both low-and high-TCR models reveal strong warming from the Equator to 30° N, and warming to greater depths in higher-latitude ocean areas (Fig. 4a,b). However, in low-TCR models, there is a greater warming at depth compared with models with high TCR (Fig. 4c). These differences exist at all latitudes down to 1,000 m depths, although differences are weakest in regions of strong warming from the Equator to 30° N.
When integrated across the volumes of water for OHC (Fig. 4d), the Southern Hemisphere indicates stronger intermodel contrasts than the Northern Hemisphere. The strongest contrasts are found near the equatorial undercurrent, the origin of which remains an open question. The differences in the North Atlantic intermodel contrasts are likely to be due to the overturning circulation response to a warming climate.
A key role for changes in upper-ocean salinity driven by both the surface freshwater flux and winds has been identified, with consequences for both the patterns and magnitudes of warming in coming decades 208,209 . Given the role of salinity and the fact that key cryosphere processes are absent from the current generation of climate a | Zonal mean ocean warming with depth normalized by global-mean near-surface air temperature warming from years 1 to 70 of 1% CO 2 simulations from CMIP6 for five models with the lowest TCR. b | As in a, but the five CMIP6 models with the highest TCR. c | The difference between a and b. d | Associated differences in heat content change between the lowestand highest-TCR models. Higher efficiency of ocean heat transport from the upper into the deeper ocean has a negative effect on the surface temperature increase for transient climate response.
www.nature.com/natrevearthenviron models used for future projections 210 , important caveats exist regarding the ability to accurately estimate TCR and the latitudinal structure of future ocean warming. Nevertheless, the OHC differences between low-and high-TCR models reflect the higher OHU of low-TCR models (Fig. 4d).

Regional trend projections.
As discussed, the pattern of ocean warming has been non-uniform 4 (Fig. 3). Accordingly, although all six regional basins are expected to warm throughout the twenty-first century -exceeding 0.2 °C relative to a 2005-2019 baseline in the upper 1,200 m and at least 0.05 °C at 2,000 m -regional contrasts in rates of warming can be anticipated. These regional trends are now discussed, basing analyses on drift-corrected CMIP6 model projections without observational constraints, as the emergent constraint method only applies for global OHC change (Fig. 3c-h). Thus, regional biases in models have not been fully accounted for 211,212 .
In the Pacific, total upper-2,000-m OHC change by the end of this century is projected to be 344 ± 74 ZJ for SSP1-2.6 and 575 ± 95 ZJ for SSP5-8.5 (Fig. 3c), ~3-6 times and ~6-10 times the observed changes over 1958-2019. Accordingly, the Pacific is expected to become the biggest heat reservoir owing to its large volume. However, the historical condition of low area-averaged warming continues. By the end of the century, area-averaged 0-2,000-m OHC reaches 2.37 ± 0.51 × 10 9 J m −2 and 3.97 ± 0.66 × 10 9 J m −2 for SSP1-2.6 and SSP5-8.5, respectively (Fig. 3c). By 2100, the 0.2 °C warming isotherm, a rough indicator of the vertical extension of warming, is projected to reach depths of ~1,200 m under SSP1-2.6. This depth is shallower than all other basins, likely because the Pacific Ocean is dominated by shallow tropical and subtropical overturning cells, and ventilation of thermocline water (~150-200 m) via subduction in the subtropical gyres 213 . However, as current climate models are not skilful at capturing ENSO 214,215 and IPO 124 variability, future projections remain uncertain 211,214 .
Similar to its past changes, the Atlantic Ocean is projected to undergo substantial warming. Total upper-2,000-m OHC increases are estimated to be 299 ± 84 ZJ for SSP1-2.6 and 546 ± 122 ZJ for SSP5-8.5, ~1-4 and ~3-6 times the observed 1958-2019 change, respectively (Fig. 3d). This heat gain is broadly consistent with that of the Pacific, despite representing only 57% of its area. Accordingly, the Atlantic witnesses some of the largest area-averaged rates of warming at 0-2,000 m: 3.62 ± 1.02 × 10 9 J m −2 for SSP1-2.6 and 6.61 ± 1.48 × 10 9 J m −2 for SSP5-8.5 (Fig. 3d). Moreover, the 0.2 °C anomaly isotherm reaches ~2,000 m and maximum basin-mean subsurface warming exceeds ~1 °C at 200-600 m by the end of the century under SSP1-2.6 (~1.6-2.8°C under SSP5-8.5). These large changes have been attributed to a decline of aerosol impacts in future Atlantic changes, reinforcing the greenhouse effect 216,217 .
In the Indian Ocean, total 0-2,000-m warming is projected to be 108 ± 34 ZJ and 183 ± 43 ZJ by 2100 for SSP1-2.6 and SSP5-8.5, respectively (Fig. 3e), ~2-5 times and ~4-7 times that observed over 1958-2019. The area-averaged changes resemble those of the Pacific, and are projected to be 2.35 ± 0.73 × 10 9 J m −2 for the low-emission scenario, and 3.99 ± 0.95 × 10 9 J m −2 for the high-emission scenario (Fig. 3e). However, the uncertainty range is larger than projected for the Pacific, partly associated with substantial natural variability, and consistent with a large mismatch between observed and simulated Indian Ocean warming from 1958 to 2019. By 2100, the 0.2 °C warming isotherm reaches depths of ~1,500 m, deeper than the Pacific Ocean but shallower than the Atlantic Ocean (Fig. 3e).
As in the Atlantic, projected changes in the southern oceans are also large. Net upper-2,000-m warming is projected to reach 220 ± 60 ZJ for SSP1-2.6 and 505 ± 95 ZJ for SSP5-8.5 (Fig. 3f), ~1-3 times and ~3-5 times the 1958-2019 change, respectively. These net increases equate to area-averaged changes of 2.44 ± 0.67 × 10 9 J m −2 and 5.60 ± 1.05 × 10 9 J m −2 , stronger than the Pacific and Indian Oceans, but increasingly smaller than the Atlantic Ocean, contrasting with observed changes. This large ocean warming is linked to accelerating subpolar westerlies and significant overturning of surface heat due to the formation of subantarctic mode water 164,165 , Antarctic intermediate water 218 in the subpolar ocean, and Antarctic Bottom Water around the Antarctic continental margin 219 , enhancing OHU and accumulation of heat just north of the ACC 163,216,220 . These water mass formations and vertical heat transports are reflected in deep-reaching warming by 2100, with the 0.2 °C anomaly isotherm expected to extend below 1,400 m (Fig. 3f). The uptake and drawdown of heat to depth is expected to slow around the Antarctic margin as meltwater and warming stratifies the surface mixed layer in those regions 181 and slows the formation of dense shelf water, but warming is expected to continue over much of the southern oceans.
The Mediterranean OHC increase at 0-2,000 m is expected to reach 6 ± 2 ZJ and 12 ± 2 ZJ for SSP1-2.6 and SSP5-8.5, respectively 221 . These changes are ~1-3 times and ~2-5 times observed values over 1958-2019, and reflect area-averaged changes of 2.06 ± 0.60 × 10 9 J m −2 and 4.09 ± 0.72 × 10 9 J m −2 (Fig. 3g). The 0.2 °C warming anomaly is projected to extend to ~1,200 m, which is deeper than for the global ocean 182 (Fig. 3g), owing to the thermohaline circulation in the Mediterranean Sea and efficient heat propagation to deeper layers 222,223 . Finally, in the Arctic Ocean, upper-2,000-m warming by 2100 is projected to reach 36 ± 22 ZJ for SSP1-2.6 and 76 ± 35 ZJ for SSP5-8.5 (Fig. 3h). These changes correspond to surface warming ~1-7 times and ~4-13 times that observed, totalling 2.50 ± 1.54 × 10 9 J m −2 and 5.25 ± 2.45 × 10 9 J m −2 for the low and high scenarios, respectively. Moreover, the 0.2 °C additional heat isotherm extends to ~1,800 m by the end of the century, and the >1 °C additional warming to within 100-800 m, spanning a larger depth than the Atlantic (Fig. 3h). The subsurface warming in the Arctic Ocean is much stronger than the global mean 224 (Fig. 3b,h). However, both the observational record and future warming in the Arctic are less reliable, given poor data coverage at the subsurface and the biases in Arctic Ocean hydrography and ice simulation in CMIP6 models 224 .
Drivers of future OHC patterns. The patterns of ocean warming are driven by the heat uptake at the ocean surface and the transport by currents, mixing and stirring in the ocean interior. The patterns of ocean warming can be thought of as the sum of the added heat from the warming pattern set by surface heat fluxes and by the unperturbed ocean circulation (analogous to the uptake and transport of passive tracers), and the redistributed heat as the transport of the OHU by the perturbed circulation 27,127 . These added and redistributed contributions to the total regional heat content changes can be quantified using a subset of models with carbon cycle 27 (Fig. 5).
Over the historical period (here, 1958-2015), the total heat patterns are dominated by the redistribution, with a pattern correlation between total and redistributed heat of R = 0.84 on average (Fig. 5d,j and Supplementary  Fig. 9). Indeed, warming along the Kuroshio extension, the southern oceans and to some extent the high-latitude North Atlantic are associated with redistribution of the pre-1958 temperatures from the changing ocean circulation in these regions, in broad agreement with observational inferences during 2000-2015 (rEFs. 125,225 ). Note, however, that the redistribution pattern varies widely across CMIP models. In contrast, warming in the subtropical Atlantic basin, which has emerged over the   models. c | As in b, but centred on 2100, representing the 5-year mean from 2095 to 2100. d | As in a, but normalized by global ocean heat uptake. e | As in b, but normalized by global ocean heat uptake. f | As in c, but normalized by global ocean heat uptake. g | As in a, but for added OHC normalized by global ocean heat uptake. h | as in b, but for added OHC normalized by global ocean heat uptake. i | As in c, but for added OHC normalized by global ocean heat uptake. j | As in a, but for redistributed OHC normalized by global ocean heat uptake. k | As in b, but for redistributed OHC normalized by global ocean heat uptake. l | as in c, but for redistributed OHC normalized by global ocean heat uptake. Stippling indicates where the multimodel mean anomaly is less than the intermodel standard deviation. See Supplementary Information for list of models used. The patterns of ocean warming are mainly driven by the redistributed heat as the transport of the OHU by the perturbed circulation in the historical era, but are projected to be dominated by added heat in the twenty-first century, set by surface heat fluxes and the unperturbed ocean circulation.
www.nature.com/natrevearthenviron cooling effect of redistribution, is primarily associated with added heat in the CMIP6 ensemble (Fig. 5g).
Between 2015 and 2050, under SSP5-8.5 forcing, a warming signal is expected to emerge in most regions of the ocean (Fig. 5b,e and Supplementary Figs. 10 and 11). This overall warming pattern is primarily set by the passive uptake of added heat, with pattern correlations between total and added heat being 0.67 (Fig. 5e,h). The added heat pattern is likely to be dominated by warming within the subtropical gyres in all basins, on the equatorward flank of the ACC due to Ekman transport 156,163,226 , and in the North Atlantic between ~30° N and 55° N due to heat uptake at high latitudes and southward transport at depth by the overturning circulation 69,227 . Exceptions to this overall warming pattern are evident in the North Atlantic warming hole and parts of the southern oceans, where cooling prevails (Fig. 5e). Trends in both these regions continue to exhibit a large spread across CMIP6 models and are dominated by redistribution patterns (Fig. 5k). Moreover, some of the added warming is partially offset by redistribution within the subtropics (~10°-30° N/S), consistent with CMIP5 27 .
By 2100, all regions of the ocean are projected to exhibit warming under SSP5-8.5 (Fig. 5c,f and Supplementary Figs. 10,11). At this time, the southern oceans and the Atlantic basin, along with most other regions, are heavily influenced by added heat (Fig. 5i). However, the redistribution of heat and circulation changes reduce net warming in the lower-latitude subtropics between ~10° and 30° N/S, and enhance warming in sectors of the mid-latitudes (between ~40° and 50° N/S) and along the equator (between ~10° S and 10° N) (Fig. 5l). This redistribution at 2100 is broadly consistent with the regional effects projected to emerge by 2050, with correlation patterns of 2050 and 2100 redistribution being 0.82. Idealized forcing experiments (flux perturbations under the flux-anomaly-forced model intercomparison project protocol [228][229][230] ; 1% yr −1 in CMIP5 231 ), also reinforce this finding. Pattern correlations between 2050 and 2100 OHU (Fig. 5e,f), added heat (Fig. 5h,i) and redistribution (Fig. 5k,l), are also high, at >0.8, >0.85 and >0.75, respectively, indicating that future OHC patterns have already emerged by 2050.
Thus, ocean heat redistribution, although dominant in the past, becomes less important in the future (at least by 2050). Future patterns become more dominated by added heat associated surface heat flux and mean ocean circulation. As such, future ocean warming patterns might become more predictable in the coming decades 27 . However, further work is necessary to assess these results in high-resolution eddy-rich models.

Impacts and consequences of ocean warming
Observed and projected ocean warming has substantial impacts across major Earth system components and scales (Fig. 1). For example, ocean warming accounts for more than 1/3 of global-mean sea level rise through thermal expansion, and thus dominates regional sea level patterns 82,202 . Sea level rise, in turn, increases the risks for coastal infrastructures and coastal habitats from salt water intrusion, coastal erosion and flooding in low-lying regions 232,233 (Fig. 1b). Ocean warming also decreases ocean density and increases upper-ocean stratification by 5.3% since 1960 (rEF. 234 ) (Fig. 1b), affecting the vertical and lateral exchanges of heat, carbon, oxygen, nutrients and other substances. The stratification increase, solubility reduction and circulation changes drive deoxygenation in the ocean interior by ~0.5-3.3% since the 1960s (rEFs. 6,235 ). By changing sea water buoyancy, ocean warming impacts ocean currents, for example, accelerating the zonally averaged Southern Ocean zonal flow in the upper layer 174 . Warmer water also reduces the efficiency of oceanic carbon uptake and storage 28,236 (Fig. 1b). The compounded effect of each of these impacts, especially following extreme events [236][237][238] , poses more substantial stress on the environment than their individual effects 236 (Fig. 1b), driving, for example, changes in net primary and export production 239,240 with socioeconomic impacts on marine fisheries and aquaculture systems 241,242 . Indeed, it is projected that, driven by multiple stressors, the maximum catch potential of tropical fish stocks in some tropical exclusive economic zones will decline by up to 40% by the 2050s under the RCP8.5 242 .
MHWs offer a strong example, whereby the relentless increase in OHC has direct implications for the frequency, intensity and extent of MHWs and other 'hot spots' within the ocean (Fig. 1b). With human-induced global warming and higher ocean heat content, it is inevitable that MHWs become more abundant, extensive and longer-lasting 237 . The highly anomalous ocean waters, including SSTs and upper OHC, often persist for more than a month 243 , resulting in large impacts on ocean ecosystems and marine life. Effects from thermal stress causes mass mortality of benthic communities, including coral bleaching, changes in phytoplankton blooms, adverse effects on kelp forests and sea grasses, toxic algal blooms, shifts in species composition and geographical distribution, and decline in fish and fisheries and seabirds 105,238 . As such, MHWs modify ecosystem assemblages, biodiversity, population extinctions and redistribution of habitat 244,245 .
One example is the prolonged MHW known as 'the blob' in the northeast Pacific and Gulf of Alaska from 2014-2016. This event greatly affected the ocean food web, shrinking phytoplankton blooms, which, in turn, diminished copepods, zooplankton and krill, and small fish, leading to the demise of ~1 million birds (notably murres), ~100 million cod and hundreds of humpback whales 104,105 . A similar unprecedented MHW in the south Tasman Sea in 2015-2016, where SSTs were up to 2.9 °C above normal owing to an ENSO-related alteration of the ITF and East Australian Current 119 , also had substantial impacts 246 . Ecosystem impacts ranged from new disease outbreaks in farmed shellfish (oysters, abalone) and salmon, to mass mortality of wild mollusks, and many out-of-range observations of several fish species. Thus, with continued warming, MHW events and their impacts are expected to worsen: climate models project that the frequency of MHW might increase 50 times by 2081-2100 relative to 1850-1900 under RCP8.5 (rEF. 247 ).
Ocean warming also intensifies tropical cyclones 248 , and associated changing ocean surface currents can indirectly affect pathways of storms 139 (Fig. 1b). In August 2017, the Gulf of Mexico became the warmest on record NATure revIeWS | EARTH & EnvIROnmEnT to that point in the summertime. There was a strong link between upper OHC and record high rainfalls of over 60 inches (1,520 mm) over five days and extensive flooding in hurricane Harvey over parts of Texas 249 . Other processes at the air-sea interface affected by ocean warming include an increase of surface evaporation and rainfall 12,13 , and an increase in precipitation anomalies tied to ENSO, and associated extreme weather events 250,251 (Fig. 1b).
Implications of ocean warming are also widespread across the Earth's cryosphere 195 , and have in turn affected the ocean itself 252 (Fig. 1b). Examples include the thinning of floating ice shelves and marine terminating glaciers from basal ice melt 21,253 , and the retreat and speedup of ice-sheet outlet glaciers in Greenland and Antarctica 254 and of tidewater glaciers in South America and in the Arctic 255 .
Other particular concerns are the potential abrupt changes associated with warming, such as ocean circulation (for example, AMOC) 256 and ocean ecosystems 6,248 . Ocean warming is a key driving element for AMOC changes. The potential for abrupt AMOC collapse as a 'low-probability, high-impact' event cannot be ruled out in the future 6,247 . Each species of marine organisms has an optimal temperature window for functioning; most organisms are therefore vulnerable to warming 257 . It is projected that most tropical coral reefs will be threatened 258 .

Summary and future perspectives
In summary, OHC has changed substantially since the 1950s and is projected to continue to do so in the future. In the upper 2,000 m, net global increases of 351.4 ± 59.8 ZJ (0.36 ± 0.06 W m −2 ) have been observed from 1958 to 2019, with the rate of warming accelerating from ~0.0-0.3 W m −2 in the 1960s to ~0.5-0.7 W m −2 in the 2010s. The pattern of ocean warming has been non-uniform in this historical era, including strong warming in the Atlantic and southern oceans, and overall is dominated by the redistribution of ocean heat by currents. Relative to 2005-2019, future warming is projected to reach 1,030 [839-1,228] ZJ for SSP1-2.6 and 1,874 [1,637-2,109] ZJ for SSP5-8.5 at the end of this century, ~2-4 (SSP1-2.6) to ~4-6 times (SSP5-8.5) the observed 1958-2019 change. On these timescales, added heat has an important role for OHC projections. Moreover, low GHG emissions would be likely to lead to a detectable and lasting reduction in ocean warming rate, with noticeable reductions in climate-change impacts. Indeed, given that ocean warming has already led to pervasive impacts and consequences, monitoring, understanding, adapting to and mitigating ocean warming must continue to be a high priority. Nonetheless, several gaps remain in measuring, estimating and understanding ocean warming.
First, the current ocean observing system needs to be sustained and extended to monitor OHC change at various spatiotemporal scales. The critical question of 'how adequate is the ocean observing system for monitoring the OHC changes at various timescales?' is still not fully answered. The international community have cast their eyes toward the future to improve not only the Argo fleet but also other measurement methodologies [259][260][261] . The ongoing and planned efforts include the maintenance of the current GOOS, and shipboard in situ measurements for calibration, validation and quality control of the Argo array; a drive toward spatial completion, including polar sea-ice zones, marginal seas and complicated channels; increased resolution in critical areas such as boundary currents and coastal areas; incorporation of deeper measurements (for example below 2,000 m); and inclusion of biological and chemical signals, along with temperature and salinity (Deep Argo and BioGeoChemical Argo programmes). The scientific community and funding agencies will need to be mindful of ensuring a continuous and comprehensive measurement network for the world's ocean -a network that incorporates new technologies as they are developed and retires old technologies that have outlived their usefulness, but with a good understanding of the handoff between technologies.
Second, uncertainty for OHC estimate needs to be better understood and quantified. In addition to data coverage, uncertainty also stems from mapping methods associated with data sampling, instrumental bias correction, choice of climatology, quality-control and other data processing procedures. These error sources are not independent of each other and are likely to lead to biases in current analyses. Thus, uncertainty in OHC estimate is method-dependent and producer-dependent. The contribution of each error source to the total OHC uncertainty is not fully understood, and the error range given by different datasets results in roughly 10-fold differences 33 . New approaches can be used to better quantify the uncertainty: for example, synthetic data for understanding and evaluating the mapping method, and exploiting ensembles (applying different techniques and forming an ensemble of OHC estimates) 59,[262][263][264] .
Third, syntheses of multisource (direct and indirect) observations and models are recommended to improve OHC estimates and mechanistic understanding. Synthesis of in situ observations with satellite-based observations (sea level altimetry, ocean colour, surface wind stress) and full atmospheric analyses shows the most promise. The indirect datasets can either serve as a cross-evaluation tool or constrain direct estimates, such as closing energy, water and sea level budgets. Attempts show promising results. As a coupled system, the separate impacts of atmospheric, ocean and ice (and other components) dynamics cannot be easily isolated 265 , and identifying the coupling mechanism driving OHC patterns remains a high research priority. Capabilities for integrating different sources of Earth system observations and information (for example, model-based data assimilation and simulations) for a comprehensive quantification of the energy budgets should be built. For example, the integration of atmospheric and oceanic data leads to a quantification of MHT time series at all latitudes, capable of resolving interannual variability 133 .
Fourth, the current generation of climate and Earth system models still contain non-negligible uncertainties in representing past and future ocean warming trends 202,205,212 . For example, there is substantial uncertainty in CMIP6 future projections in the Antarctic margin region, due to biases in simulated stratification 158 , hydrography around the Antarctic shelf 266 , and missing processes, including eddies, tides, ice-shelf cavities, and ocean-ice-shelf interactions. In the tropics, it is known www.nature.com/natrevearthenviron that tropical cyclones mix the ocean through substantial depths. Moreover, they form in hot spots and cool the ocean, yet they are largely absent from global models. In addition, interannual to decadal variability such as ENSO and the IPO remains poorly represented in modern-day climate models, yet they strongly control the pattern and timing of OHC anomalies. As ocean and atmospheric circulation have a key role in shaping the climate response (including ocean warming pattern and OHU efficiency) 127,203 , evaluating the accuracy of the wind, atmosphere and ocean circulation projections should continue to be a research priority. A continuous process-based understanding of model performance is recommended including the understanding and identification of persistent biases in simulations, especially with respect to the diagnostics for Earth's major system cycles.
Fifth, understanding of extreme OHC events, their compound effects and past changes in OHC should be strengthened. For example, MHWs have been identified from surface conditions, but subsurface components are also important, and OHC as an indicator offers a way to integrate these aspects. The simultaneous occurrence of ocean warming extremes with other extremes (deoxygenation, acidification) requires special attention 236 . A more complete understanding of OHC changes in the deep past before the widespread availability of instrumental records is also critical to put the current changes in the context of a longer timescale 267 . Research on past climate change also helps us to understand how natural drivers and human influence have changed the Earth's climate system. The difficulty is a lack of full-depth temperature proxy data. Several methods have been developed to derive OHC change back to the past 20,000 years (Last Glacial Maximum), but the uncertainty is large 95,268 .

Data availability
The observation and model data used in this review are available at http://www.ocean.iap.ac.cn/. CMIP6 model data is available at https://esgf-node.llnl.gov/search/ cmip6/.
Published online 18 October 2022