Recent Water Mass Changes Reveal Mechanisms of Ocean Warming

: Over90%ofthebuildupofadditionalheatintheEarthsystemoverrecentdecadesiscontainedintheocean. Since 2006, new observational programs have revealed heterogeneous patterns of ocean heat content change. It is unclear how much of this heterogeneity is due to heat being added to and mixed within the ocean leading to material changes in water mass properties or is due to changes in circulation that redistribute existing water masses. Here we present a novel diagnosisof the‘‘material’’and ‘‘redistributed’’ contributions to regionalheatcontentchangebetween2006and 2017 that is based on a new ‘‘minimum transformation method’’ informed by both water mass transformation and optimal transportation theory. We show that material warming has large spatial coherence. The material change tends to be smaller than the redistributed change at any geographical location; however, it sums globally to the net warming of the ocean, whereas the redistributed component sums, by design, to zero. Material warming is robust over the time period of this analysis, whereas the redistributed signal only emerges from the variability in a few regions. In the North Atlantic Ocean, water mass changes indicate substantial material warmingwhile redistribution cools the subpolar region as a result of a slowdownin the meridional overturning circulation. Warming in the Southern Ocean is explained by material warming and by anomalous southward heat transport of 118 6 50 TW through redistribution. Our results suggest that near-term projections of ocean heat content change and therefore sea level change will hinge on understanding and predicting changes in ocean redistribution.


Introduction
Over the past 50 years, as atmospheric greenhouse gas concentrations have increased, the ocean has absorbed more than 10 times as much heat as all other components of the climate system combined (Rhein et al. 2013).This warming showed substantial spatial variability between 1993 and 2005, being up to 10 times as much in some regions as the global average (Zhang and Church 2012).It is unclear whether this variability is due to geographical variation in the interior propagation of surface warming versus redistribution of existing heat within the ocean.
Ocean warming is an important issue because ocean thermal expansion is the largest projected contribution to global mean sea level rise in the twenty-first century (Church et al. 2013).Numerical climate models disagree on the pattern and amplitude of ocean heat content (OHC) change and hence on sea level rise under anthropogenic greenhouse warming (Gregory et al. 2016).Understanding how heat has been taken up and redistributed by the ocean is essential for predicting future changes in sea level.
Numerical ocean models forced with historical atmospheric conditions have proved to be useful tools in quantifying how variability in atmospheric forcing can set variability in OHC (Drijfhout et al. 2014) and sea level (Penduff et al. 2011) at interannual to decadal time scales.However, such models can be unrealistic for simulating multidecadal climate change because of model drift and inaccuracies in long-term changes in atmospheric forcing, particularly global mean heat fluxes (Griffies et al. 2009).On the other hand, coupled ocean atmosphere climate models are routinely used to capture the effect of long-term climate forcing.But such models only accurately simulate past unforced variability in regional OHC when, by chance, their internal variability is in phase with the observed system.
An advance in terms of numerical ocean climate modeling has come from the separation of OHC change into an ''added'' and a ''redistributed'' component in climate model simulations, where the former is due to change in the surface heat flux, and the latter due to rearrangement of existing OHC because of altered ocean heat transports (Banks and Gregory 2006).This decomposition is analogous to the ''anthropogenic'' and ''natural'' decomposition that has revolutionized our understanding of oceanic carbon records (Khatiwala et al. 2013).Here we will present a novel method to diagnose the ''material'' component of OHC change, which we will show is closely related to the ''added'' component introduced by Banks and Gregory (2006).
Recent work has aimed to reconstruct the drivers of OHC change based on observationally derived air-sea boundary conditions.Zanna et al. (2019) for example used surface temperature anomalies combined with a tracer-based approach to reconstruct the role of anomalous surface heat fluxes in centennial heat content change.Roberts et al. (2017) estimated the contribution of air-sea heat flux changes in setting mixed layer and full-depth-integrated OHC budgets over recent decades and inferred the role of ocean circulation as a residual.Here we aim to circumvent reliance on such boundary conditions and infer the mechanisms of ocean heat content change directly based on water mass changes.
Water mass-based methods have been used to decompose local temperature and salinity changes into a dynamic ''heave'' component and an apparently material component at constant density based on a one-dimensional view of the water column (Bindoff and McDougall 1994).However, their analysis did not distinguish between material processes and horizontal advection, insofar as they affect the water mass properties of an individual water column.
Here we introduce a new method based on water mass theory, called the minimum transformation method, which we use to estimate recent drivers of three-dimensional OHC change.In section 2 we will review water mass theory and establish the relationship between changes in water masses as defined by their temperature and salinity and material changes in seawater temperature.We will describe in section 3 how this theory is translated into a practical method to estimate material changes in water masses and map these into geographical space.We present an application of this minimum transformation method to recent data over the Argo period in section 4 and give results in section 5. We discuss the results and compare them with existing work in section 6, and we give conclusions in section 7.

Water mass theory
Water mass analysis has long been used in physical oceanography to trace the origin of waters (Montgomery 1958).In the latter half of the twentieth century a quantitative framework emerged to describe the relationship between water masses, airsea fluxes, and mixing [Walin 1982; see the review by Groeskamp et al. (2019)].Recent work has seen this framework advanced in two ways specifically relevant to our work here: to multiple tracer dimensions to understand the thermodynamics of ocean circulation (Nycander et al. 2007;Zika et al. 2012;Döös et al. 2012;Groeskamp et al. 2014;Hieronymus et al. 2014) and to unsteady problems to understand the ocean's role in transient climate change (Palmer and Haines 2009;Evans et al. 2014;Zika et al. 2015a,b;Evans et al. 2017Evans et al. , 2018)).
An example of the utility of the water mass transformation framework in understanding transient change is provided by Zika et al. (2015a).They demonstrate that the distribution of water in salinity coordinates is influenced by the water cycle and turbulent mixing, the latter only being able to collapse the range of salinities the ocean covers.This means that changes in the width of the salinity distribution indicate an enhancement of the water cycle and/or a reduction in that rate at which salt is mixed.In this project we extend this concept to consider how changes in the temperature-salinity distribution relate to material changes in water masses.
Material changes in Conservative Temperature Q (hereinafter simply ''temperature'' or T) following the motion of an incompressible fluid are related to Eulerian changes and advection by where u is the 3D velocity vector and DT/Dt is the material derivative, which is related to sources and sinks of heat and irreversible mixing.Conservative Temperature is used here since it is a more accurate ''heat'' variable than potential temperature (McDougall 2003), although the later is still routinely used in ocean models, including the one analyzed in section a of appendix A.
Even if a perfect record of ›T/›t were available at a fixed location, we would not know the relative roles of advection (u Á =T) and material processes (DT/Dt).To separate them, we consider the water mass perspective as an alternative to the Eulerian perspective.The following theory draws directly from Hieronymus et al. (2014).
We characterize water masses by their T and Absolute Salinity S A (IOC/SCOR/IAPSO 2010; hereinafter simply ''salinity'' or S).The volume y of water per unit temperature and salinity and at temperature T* and salinity S* is where the integral is over elements dV of ocean volume that are cooler than T* and fresher than S*.An estimate of y that is based on recent observational analysis is given in Fig. 1a.(These data are described in detail in section 4.) Considering all of the water in the ocean and retaining the incompressibility assumption, the only way y can change is via transformation-that is, by making water parcels warmer, colder, saltier, or fresher as described by the following continuity equation [derived formally in Hieronymus et al. (2014)]: where _ T is the average material derivative of T within a water mass.That is, and likewise _ S is the average material derivative of S.An estimate of recent changes in y is given in Fig. 1b.
In Eq. (3) the terms y _ T and y _ S are the transformation rates in the temperature direction (Sv g 21 kg 21 ; 1 Sv [ 10 6 m 3 s 21 ) and salinity direction (Sv C8 21 ) respectively.Equation (3) states that the amount of water between two closely spaced isotherms (T and T 1 ›T) and isohalines (S and S 1 ›T) will go up if more water is made warmer at T than at T 1 ›T and/or more water is made saltier at S than at S 1 ›S.
When the system is in a statistically steady state the water mass distribution y remains constant such that where the overbar represents a sufficiently long time average.In this steady case, the vector field described by y _ T and y _ S can be characterized by a thermohaline streamfunction (Zika et al. 2012;Groeskamp et al. 2014).
Here, we will not attempt to estimate this steady-state component of water mass transformation [e.g., as Groeskamp et al. (2017) have done].Rather we will attempt to quantify only the component required to explain changes in y.That is, we aim to quantify the anomaly in the transformation rate (y _ T) 0 such that y _ T 5 y _ T 1 (y _ T) 0 , and likewise for (y _ S) 0 , with Note that a steady-state component like Eq. ( 5) can always be added to (y _ T) 0 and (y _ S) 0 such that Eq. ( 6) is still satisfied.
However, we seek only the net change in water mass transformation required to explain changes in y and therefore seek the smallest (in a root-mean-square sense) values of _ T 0 and _ S 0 that satisfy Eq. ( 6).That is, we seek the smallest change in air-sea heat and freshwater fluxes and mixing-in a net sense-that can explain changes in water masses.We call this the minimum transformation.
Here we will use changes in y to infer the minimum transformation and therefore estimate y _ T 0 .This will allow us to estimate the material processes influencing ocean temperature change.

The minimum transformation method
We now apply water mass theory to understand changes in a discrete set of water masses describing the ocean over two time periods.We will then describe the application of a minimum transformation method that exploits an ''earth mover's distance'' (EMD) algorithm to estimate the amount of material warming required to affect changes in those water masses.

a. Discrete water masses
Consider the set of N discrete water masses with the ith water mass defined by the limits [T min i , S min i , ]. Essentially, our water masses are hypercubes in T-S-x-y-z space (more arbitrary space-and time-dependent regions can be defined without affecting the method described below).To indicate whether water is within the ith water mass we define a boxcar function P i such that The volume of water in the ith water mass at time t is then ÐÐÐ P i (x, t) dV.
We consider two time periods: an early period (t 0 2 Dt # t , t 0 ) and a late period (t 0 # t , t 0 1 Dt).The average volume of the ith water mass over the early period is V1 i and the average volume of the jth water mass over the late period is V2 j such that V1 To change the set of volumes V1 i into the set of volumes V2 j requires a transformation of water in T-S space.When water transforms, it changes its T and S and can also move geographically.
To understand how water is transformed from the physical location and physical properties of one water mass to another we use the shorthand x(t 1 Dtjx, t) for the position of a water parcel at time t 1 Dt conditional on it previously being at position x at time t.That is, where, as previously, u is the 3D velocity vector.We describe the transformation rate between the early and late water masses with the matrix g.The ith column and jth row of this matrix g ij correspond to the average rate of transformation of water from early water mass i to late water mass j such that In Eq. ( 12) the term P i (x, t)P j [x(t 1 Dtjx, t), t] isolates water that was in the ith water mass at time t and was subsequently in the jth water mass at some time Dt later.The quantity g ij is therefore the average rate (m 3 s 21 ) at which water in the ith early water mass is transformed into the jth late water mass.
Since the total volume of water is conserved between the early and late periods all the water from the early water masses (V1 i ) must be transformed into late water masses.Likewise, all water masses from the late period (V2 j ) are made from water masses of the early period.That is, The average temperature change of water that transforms from V1 i to V2 j is then We can also write Eq. ( 14) as where T 2 ji is the volume-weighted average temperature of the water in the jth late water mass that was previously in the ith early water mass and T 1 ij is the volume-weighted average temperature of the water in the ith early water mass that is later in the jth late water mass.The transformation g ij involves a range of water parcels with a range of temperatures T(x, t), whose mean is T 1 ij , in the early period moving to a range of temperatures T[x(t 1 Dtjx, t), t], whose mean is T 2 ji , in the late period.To simplify this problem, we assume that in both periods the water masses are well mixed.This means that we expect that the mean temperature of any sample of water parcels from water mass i in the early period will equal the mean temperature of the water mass as a whole, and in particular this is true for the sample of parcels that ends up in water mass j in the late period.Thus T 1 ij 5 T1 i with this assumption.By a similar argument, T 2 ji 5 T2 j , and hence the average T and S change of water transforming from the ith early to the jth late water mass as the difference of the average T and S of the two water masses.That is, DT ij 5 T2 j 2 T1 j and DS ij 5 S2 j 2 S1 j .
This above approximation preserves the following equality relating the change in global volume-weighted temperature to the transformation matrix: and likewise for the volume-weighted salinity.
We have effectively discretized the continuum of trajectories from early to late water masses into a finite set of discrete trajectories.This discretization clearly leads to some information loss; however, such losses are unavoidable in any computationally feasible inverse method.
Note that, even if the ith water mass for the early period has the same temperature and salinity bounds as the ith water mass of the late period, the distribution of properties within the water mass can change.That is, in general T1 j 6 ¼ T2 i and S1 j 6 ¼ S2 i , so g ij is always a transformation, even with i 5 j.For example, assume the ith water mass has temperature bounds 18 and 28C and that the water between those bounds is on average at 1.98C in the early period and 1.18C in the late period.Groeskamp et al. (2014) called this a ''local effect'' and included it as an separate term in their formulation.Here, we find it convenient to consider the transformation from the ith early water mass at 1.98C to the ith late water mass at 1.18C to be yet another transformation-no different than between any other pair of water masses.
We relate the transformation rate to the average material temperature tendency required to warm the ith early water mass to form the range of destination water masses it arrives at in the late period.That is, We use _ T to define a 3D material temperature change field DT Material such that Note here that we are relating _ T i only to the anomaly of the Lagrangian tendency (i.e., DT 0 /Dt rather than DT/Dt) as it appears in Eq. ( 19).This is because our _ T i describes only the changes in the transformation rate required to explain changes in the water mass distribution [as in Eq. ( 6)].There can be (and indeed is) an additional ''mean'' transformation rate that leads to cycles of water in T-S space but does not lead to any changes in water mass inventories with time (Groeskamp et al. 2014).Implicit in Eq. ( 19) is the assumption that the anomalous warming of a particular water mass occurred evenly (in a volume-and time-weighted sense) over the regions and times during which that water mass existed in the early period.
We will contrast the inferred material warming at one location x against the total warming DT(x) 5 with the residual of the two being a redistribution component such that By construction, DT redistribution accounts for the advective redistribution of temperature (u Á =T), which does not affect the underlying water masses and therefore is not accounted for in DT material .

b. Finding the minimum transformation using an EMD algorithm
Our goal now is to estimate the transformation matrix g.Out of the infinite number of choices that could satisfy Eq. ( 13), we will look for the smallest (in a least squares sense) possible transformation required to change the distribution.We call this the minimum transformation.
Previous studies have diagnosed transformation rates from time-dependent changes in water mass distributions by searching for a minimum least squares solution on a regular T-S (Evans et al. 2014) or density-spiciness grid (Portela et al. 2020).Because of the dramatic variations in volume per unit temperature and salinity of the World Ocean (Fig. 1b) we choose to describe the distribution in an unstructured way.Furthermore, we exploit recent advances in the area of ''optimal transportation theory''-in particular, the EMD algorithm that is mentioned at the beginning of section 3 (Pele andWerman 2008, 2009).
The EMD solves the hypothetical problem of moving earth from a set of mounds, each with varying amounts of earth, into a set of holes with varying amounts of empty space to be filled, where the total volume of the mounds is equal to that of the holes.In our case the ''mounds'' are the early water masses and the ''holes'' are the late water masses.The optimization problem is to find the set of transfers (from a mound to a hole, or the early to late water masses) that gives the smallest possible total of mass-weighted distance (the product of the mass and the distance of a transfer) that needs to be traveled in order to empty the mounds and fill the holes.For the EMD algorithm, we require a distance metric d, which is a matrix whose ith column and jth row d ij is the cost of moving water from the ith early water mass to the jth late water mass.The EMD algorithm then estimates g such that Eq. ( 13) is satisfied and the following total mass-weighted ''distance'' is minimized: We use the following distance metric: where temperature and salinity differences are squared so that the distance is positive definite and long trajectories in T-S space are penalized more than short ones and a is a constant that scales the salinity change relative to the temperature change and whose choice is described in the next section.The intent of d ij is to permit movement between water masses that are adjacent geographically without additional penalty but at the same time to stop direct exchange between geographically disconnected water masses, for example between water masses in the Southern Ocean and the Arctic.To achieve this we set d ij 5 0 where the ith and jth water masses are in the same or adjacent geographical regions and ] 2 } otherwise (in practice we use d ij 5 10 6 in the latter case).Regions that share a meridional or zonal boundary are considered to be adjacent.The Arctic and North Pacific Oceans are not considered to be adjacent, whereas the Indian Ocean and equatorial Pacific regions are considered to be adjacent.Our motivation for using EMD is simply to find the smallest amount of transformation (in a least squares sense) required to explain observed water mass change.If T-S changes in the ocean could be explained purely by adiabatic redistribution of existing water masses, then our method would prioritize this solution.Our initial guess is therefore this adiabatic solution (i.e., where g ij 5 0 for all i and j).The EMD algorithm finds the smallest deviation possible from this adiabatic case.We cannot rule out larger compensating transformations having taken place.In principle, solutions given different initial guesses (e.g., an initial guess for g that is based on a numerical simulation) could be explored.We leave this to future work.
Figure 2 summarizes the minimum transformation method schematically.In the schematic just four early and four late Unauthenticated | Downloaded 10/30/23 10:56 AM UTC water masses are defined with two in one geographical area and two in another.The minimum transformation moves water from the ith early to the ith late water masses in all four cases (i.e., g ii 6 ¼ 0 for all i).In addition, a substantial amount of water is moved from the second early water mass to the first late water mass (g 21 ) and from the third early water mass to the fourth late water mass (g 34 ).The observed change in temperature is therefore explained by a material warming of 28 and 18C of the two warmer shallower water masses and of 0.58C for the cooler deeper water masses.The remainder of the Eulerian pattern of temperature change is explained by redistribution.This schematic representation is vastly simplified as compared to our actual implementation of the minimum transformation method, which is described in the next section.

Data and application of the minimum transformation method
Observational estimates of T and S come from the objective analysis provided by the Enact Ensemble (V4.0, hereinafter EN4; Good et al. 2013).EN4 has a 18 3 18 horizontal resolution with 42 vertical levels.We analyze each month between 2006 and 2017 inclusive.We split these data into two time periods: an early period between 2006 and 2011 inclusive and a late period between 2012 and 2017 inclusive (i.e., t 0 5 0000 1 January 2007 and Dt 5 6 years).
We then define a discrete set of water masses for each time period by splitting the ocean into nine geographical regions and within each region by splitting up the ocean according to T-S bins.Our nine geographical regions are the Southern Ocean south of 358S, the subtropical Pacific and Atlantic Oceans between 358 and 108S, the Indian Ocean north of 358S, the tropical Pacific and Atlantic Oceans between 108S and 108N, the North Pacific north of 108N, the Atlantic Ocean between 108 and 408N, and the Atlantic and Arctic Ocean north of 408N.To avoid discontinuities in our resulting analysis we transition linearly from one region to another over a 108 band (Fig. 5).
We define T and S bin boundaries ([T min , T max ] and [S min , S max ] respectively) using a quadtree.The quadtree starts with a single (obviously oversized) bin with T boundaries [26.48, 968C] and S boundaries [25.2, 46 g kg 21 ] in which the entirety of the ocean's seawater resides.The single bin is then split into four equally sized FIG. 2. Schematic describing a simplified hypothetical implementation of the minimum transformation method.(left) Between a late and an early period, surface waters warm, especially to the south, where the ocean is fresher and the upper ocean layer becomes thicker.(center) The ocean is split into a southern region containing water masses 1 and 3 and a northern region containing water masses 2 and 4. Between the early and late periods, water masses 1 and 4 increase in volume and 2 and 3 reduce in volume.Taking into account the changing temperatures, salinities, and volumes of the early and late water masses, the ''minimum transformations'' g ij are found using the EMD algorithm.These suggest modest warming of each water mass with some of early water mass 2 transforming to become late water mass 1 (g 21 ) and some of early water mass 3 transforming to become late water mass 4 (g 34 ).(right) The total temperature change is heterogeneous.A warming of 28C explains changes in water mass 1, a warming of 18C explains changes in water mass 2, and a warming of 0.58C explains changes in water masses 3 and 4.This warming is projected onto the location of those water masses in the early period to show the ''material change.''The residual of the total and material changes is then explained by a ''redistribution'' that involves intense subsurface warming in the southern region and intense subsurface cooling in the northern region.
bins with the same aspect ratio as the original bin.The same process of splitting into four is repeated for any bin whose volume change is greater than a threshold of 62 3 10 12 m 3 (equivalent to the volume of a 58 longitude by 58 latitude region at the equator with a depth of 200 m) or until the bin size is 0.48C by 0.2 g kg 21 .Average volumes for each water mass are shown in Fig. 3.In the supplementary text we show that changing the size of these bins by a factor of 2 does not substantially change our results.The quadtree is applied within each region and for the change between the late and early periods.This results in bin edges defining N 5 1447 water masses.These bins are then used to define both the early water masses and the late water masses.
We choose the constant a to be the ratio of a typical haline contraction coefficient to a typical thermal expansion coefficient (a 5 b 0 /a 0 5 4.28).This does not mean that transformations along density surfaces are necessarily preferred; rather, the squares in Eq. ( 22) mean that density-compensated changes in T and S are penalized as much as changes of the same magnitude where one of the signs is reversed.The inferred DT material for each water mass is shown in Fig. 4. We have tested the sensitivity of our method to varying a by a factor of 2 and found only negligible changes in inferred warming (see section b of appendix A).
In section a of appendix A, we compare the results of our method applied to synthetic data from a climate model simulation with an added-heat variable explicitly simulated by the model.We find good agreement between added heat and our inferred DT material and between simulated redistributed heat and our inferred DT redistribution when ocean temperature and salinity are fed in as ''data'' to the method.Section b of appendix A also explores sensitivity of our results to parameter choices.The uncertainties we place on OHC change are 62 standard deviations of a bootstrap ensemble, also described in section c of appendix A.
To produce maps of the total, material and redistributed contributions to the heat content we multiply the density and heat capacity of seawater by the respective temperature change and vertically integrate these through the entire water column.Our method also produces a material salinity change.We leave discussion of those data to future work.

Results
Patterns of total OHC change between early and late periods are heterogeneous (Fig. 5a).There are basin-scale patches of decreasing heat content in the western equatorial and tropical Pacific, in the Pacific sector of the Southern Ocean, in the subtropical south Indian Ocean, and in the subpolar North Atlantic.Warming is seen most strongly in the tropical eastern Pacific, South Atlantic Ocean, and subtropical North Atlantic.These changes are highly sensitive to the specific observation years chosen and the length of the epochs reflecting the regional time scale of variability associated with the redistributed component.Uncertainty is far larger than the signal in the majority of regions (stippling in Fig. 5a) and coincident with previously identified regions of large sea level anomaly variability (Penduff et al. 2011).
However, there are a few regions (e.g., patches of the Southern Ocean and North Atlantic) where the regional redistributed signal is robust and emerges from the uncertainty (Fig. 5b).The patterns of redistributed heat observed in the Pacific are consistent with interdecadal Pacific oscillation (IPO)-driven thermosteric sea level variability (Lyu et al. 2017).The IPO was typically positive in the late period and negative in the early period (see https://psl.noaa.gov/gcos_wgsp/Timeseries/ for these data).
Material heat content change shows a smaller amplitude but more coherent signal than redistributed heat (Figs.5b,c).Material warming is seen across almost the entirety of the globe, with maxima in the Southern Hemisphere and Atlantic subtropical convergence zones (Maximenko et al. 2009), consistent with model simulations of passive ocean heat uptake due to anthropogenic greenhouse warming (Gregory et al. 2016).In such model simulations, anomalous heat fluxes into the ocean predominate at mid-to high latitudes and this heat is distributed throughout the ocean largely passively via subduction (downwelling) in the North Atlantic and the Southern Ocean (Marshall et al. 2015).
Strikingly, the uncertainty in material heat content change is far smaller than that of total OHC change (stippling in Fig. 5c).This suggests that heat was added to and distributed within the ocean persistently over the Argo period and that this warming is not an artifact of a particularly warm year or years.
Zonally integrating the net OHC change reveals a signal of roughly the same magnitude as its uncertainty at all latitudes (Fig. 6a).Zonally integrated redistributed heat likewise has a small signal to uncertainty ratio except in the Southern Ocean (Fig. 6a).Accumulating the redistributed heat contribution from north to south gives the meridional heat transport due to redistribution.Broadly, heat is redistributed from north to south with a southward cross-equatorial transport of 73 6 60 TW between the two epochs (Fig. 6c).
Material heat content change (Fig. 6a) is larger than its uncertainty at most latitudes and shows a peak at 358S and at 158 and 358N.The material heat content change peaks at 358S and 358N are collocated with climatological wind stress curl minima, where material warming due to anomalous surface heat fluxes may be accumulating due to convergence of surface Ekman transport.
Table 1 shows material, redistributed, and total heat content changes by ocean basin.Material heat content change is distributed among the Indian, South Pacific, and South Atlantic basins approximately according to their area.However, the tropical and subtropical North Atlantic stores close to 20% of the global ocean's material heat content change despite representing less than 10% of its area (Table 1).An outsized role for the North Atlantic in storing material heat content change in the climate system has also been foreseen in numerical modeling studies (Lee et al. 2011).
We identify robust redistributed warming signals in the subtropical North Atlantic and Southern Ocean.Warming in the subtropical North Atlantic is compensated by cooling in the subpolar North Atlantic consistent with a 40 6 13 TW southward transport of heat across 448N (Fig. 6c).Southward heat redistribution across 328S brings 118 6 50 TW into the Southern Ocean.

Discussion
Recent anomalous southward heat transport in the North Atlantic has been well documented and has been attributed to a downturn in the Atlantic meridional overturning circulation The large apparent meridional heat transport we have identified in the Southern Ocean was previously identified by Roberts et al. (2017) based on the residual of observed OHC change and estimates of air-sea heat fluxes.Their approach captures additional heat in the system where it is fluxed into the ocean while our approach estimates how that heat is distributed.Nonetheless, the correspondence between our results and theirs is reassuring and perhaps not surprising if the redistribution signal is large as both approaches indicate.
The approach of Zanna et al. ( 2019) is more directly comparable to ours.They reconstruct the passive contribution to ocean warming since 1850 by propagating SST anomalies into the ocean interior using Green's functions.They report changes for a much longer time frame (1955-2017 as opposed to our 2006-17), and therefore magnitudes of warming estimates are not comparable, but a comparison of patterns of change is relevant.In terms of our zonally averaged material warming and their ''passive warming'' the two datasets share peaks at approximately 358S and 358N potentially attributable to surface Ekman convergence (see their Fig. 3).Zanna et al. (2019) report relatively small amounts of passive warming at low-latitude regions while we report a peak in material warming there.This may suggest that the material warming we estimate at low latitudes is in fact related to interannual to decadal variability.An explanation of this may be that the lower low-latitude SST corresponds to a predominance of a negative IPO (Lyu et al. 2017), leading to anomalous ocean heat uptake over our study period.This is a commonly cited explanation for the so-called global warming hiatus discussed in the 2010s (Whitmarsh et al. 2015).Zanna et al. (2019) compare their inferred passive warming between 1955 and 2017 to the warming observed in situ.Based on this they find evidence of a southward redistribution of heat in the Northern Hemisphere but no substantial southward redistribution in the Southern Hemisphere.This suggests that the southward redistribution of heat inferred by both Roberts et al. (2017) and this study in the Southern Hemisphere may be a more recent occurrence.Indeed, two recent studies have shown that the Southern Hemisphere dominance of ocean heat content change during the twenty-first century is not consistently represented in historical climate simulations and is likely linked to internal variability (Bronselaer and Zanna 2020;Rathore et al. 2020).
Here we have exclusively analyzed the Hadley Centre's EN4 dataset.Sensitivity to observational coverage is mitigated in part by our consideration of data during the Argo observing period .We consider uncertainties to have been reasonably estimated based on our bootstrapping approach, which subsamples those years (see section c of appendix A).Because of EN4's mapping approach, however, regions where minimal observations were made (e.g., the marginal ice zones in the Southern Hemisphere and below 2000 m) will likely have muted trend estimates.This issue will require special attention when our method is applied to the pre-Argo period and in particular with regard to salinity observations, which are less numerous than temperature observations (Clément et al. 2020).

Conclusions
In summary we have shown the following: tropical North Atlantic Ocean, consistent with simulations of the addition of heat into the ocean due to greenhouse forcing.
d The majority of the variance in ocean heat content change at scales of 18 3 18 over that period can be explained by a redistribution of existing water masses within the ocean.
d The inferred redistribution indicates a downturn in northward meridional heat transport into the subpolar North Atlantic of 40 6 13 TW and an anomalous southward heat transport into the Southern Ocean of 118 6 50 TW.
The material warming signal that we have inferred is generally weaker than redistribution, but the signal is far less sensitive to changes in the years over which the analysis was carried out.This suggests that material warming may be giving a robust indication of slow thermodynamic changes in the ocean.This could be a result of anthropogenic forcing, although that would be remarkable since the midpoints of the early and late periods are only 6 years apart.
We expect the strength of the material warming signal to increase into the future as the ocean warms.However, since the redistribution signal is so large, circulation changes and variability must be understood if near-term ocean temperature variability and regional sea level change are to be projected accurately.Accuracy of the analysis we have presented in this paper relies on the following assumptions: 1) the mapping from transformations in T-S space for each region to local changes in geographical space is accurate, 2) the ''minimum transformation'' inferred using the EMD algorithm, including our choice of distance metric, accurately estimates the net thermodynamic transformation, 3) the resolution of our T-S grid is sufficiently fine to capture relevant water masses, and 4) the density of observations and the procedure used to map them onto a regular grid are sufficiently accurate for us to quantify changes in water mass volumes.
We investigate the impact of each of these assumptions in appendix A. We investigate assumptions 1 and 2 using synthetic data from a climate model in which ''added heat'' is explicitly simulated (section a of appendix A), and we investigate assumptions 3 and 4 using sensitivity tests (sections b and c of appendix A).A bootstrap approach is taken in the latter case to derive uncertainty estimates.

a. Assessment of the minimum transformation method using synthetic data
We use synthetic data from the Hadley Centre Climate Model, version HadCM3 (Gordon et al. 2000), to assess the minimum transformation method.Specifically, we exploit the configuration used for the Flux Anomaly Forced Model Intercomparison Project (FAFMIP; Gregory et al. 2016).We will consider two specific model experiments used by FAFMIP: piControl, which is a reference experiment with no external forcing, and FAFheat, in which the ocean is warmed by an imposed surface heat flux.

TRACERS
In HadCM3, the Lagrangian derivative of seawater potential temperature, T (note here that we use potential temperature rather than Conservative Temperature because the HadCM3 conserves potential temperature), is set by sources and sinks of heat Q, predominantly at the air-sea interface, and the divergence of parameterized diffusive temperature fluxes F such that As we discussed in section 3, the minimum transformation method is used to estimate the anomaly in DT/Dt with respect to a statistically steady time average.This anomaly can be related to the anomaly in heat sources and sinks Q 0 and diffusive temperature fluxes F 0 such that In the HadCM3's FAFMIP simulations an ''added temperature'' tracer T added is simulated; T added is simulated as a passive tracer initialized at zero and forced at the ocean boundary by the imposed heat flux anomaly Q* and with time-evolving diffusive flux F added such that An additional ''redistributed temperature'' tracer T redist is furthermore defined such that T 5 In practice Q 0 6 ¼ Q* in the FAFMIP experiments discussed here.This is because the net surface flux responds to changes in T redist at the sea surface.This has a large influence in the North Atlantic where anomalous ocean warming leads to a slowdown in the AMOC and therefore to a reduction in T redist at subpolar latitudes (Gregory et al. 2016).Indeed, unlike the redistributed heat inferred using our method, T redist , as defined in FAFMIP, can be a net nonzero contributor to ocean heat content.Also, F 0 redist 6 ¼ 0 since changes in circulation lead to changes in the diffusive flux with time.Furthermore, we are not able to average DT added /Dt along the pathways connecting early and late water masses as would be required for a perfect comparison between model ''truth'' and the inferences of the minimum transformation method.Despite the above caveats, we consider it worthwhile to assess our method by comparing the average change in T added over water masses to our inferred DT material .

2) ASSESSMENT BASED ON SYNTHETIC DATA
There are two aspects of the minimum transformation method that we aim to assess using these data: the uncertainty introduced by 1) projecting an inferred warming signal from temperature and salinity classes (water masses) to the geographical location of those water masses and 2) using the EMD algorithm.
The FAFMIP protocol does not describe historical climate change but rather an idealized increase in ocean heat content as would be expected from a doubling in atmospheric CO 2 .Our observational record is centered on the beginning of 2012 when the global atmospheric CO 2 concentration reached 392 ppm (Conway et al. 1994), which is approximately 40% above preindustrial levels of approximately 280 ppm.Although no comparison can be perfect, we consider this reasonable motivation to choose years 35-46 of the FAFMIP experiments to test our method.

PROJECTION
Figure A1a shows the column integral of the added-heat tracer for years 41-46 for the HadCM3 FAFheat experiment (the tracer is represented in kelvins but is here converted to the more familiar unit of watts per meter squared by multiplying by the heat capacity and density and dividing by 43 years).As was done to the EN4 data, we selected water mass bins using a quadtree approach.Figure A1b shows column-integrated added-heat change between years 41 and 46, but in this case the added-heat tracer is first averaged within each water mass within each of the nine geographical regions and then is projected back into the location of those water masses.What this projection amounts to is simply homogenizing the added-heat tracer within each water mass in each region.If added-heat change varies substantially within a water mass this method will smooth out those variations.
The information loss in the reprojection is difficult to discern between Figs.A1a and A1b, particularly in the Southern Ocean and Indian and Pacific Ocean basins.In the North Atlantic Ocean, simulated added heat is concentrated farther north than in the homogenized fields.In the zonal mean (Fig. A1c) the reprojected added heat has an RMS error of 0.5 TW (8lat) 21 with differences of up to 2 TW (8lat) 21 in the subtropical Northern Hemisphere.The mismatch in the North Atlantic is possibly due to water masses with the same T-S properties being distributed between the subpolar and subtropical regions, and it may be fruitful to distinguish between water masses in alternative ways in future.

TRANSFORMATION
We will test the minimum transformation method in the following three scenarios: 1) added heat only: heat is added to the ocean and water masses are not redistributed, 2) redistribution only: no heat is added and water masses are redistributed, and 3) added and redistributed heat: heat is added and water masses are redistributed.
Table A1 details the way data from piControl and FAFheat are used for these scenarios.
(i) Scenario 1 In this scenario there is no explicit ''redistribution'' signal in the model data.The purpose of this validation is to see how much of the change is attributed to material heat content change using the minimum transformation method.In the zonal mean (Fig A2a) the difference between the simulated and inferred added heat (which is precisely the inferred redistributed heat) has an RMS of 1.8 TW (8lat) 21 .
(ii) Scenario 2 In this scenario there is no explicit added-heat signal in the model data.This is simply a climate control run with no variations in forcing (solar, aerosol, etc.).There are, however, some very small changes in ocean heat uptake due to natural variability in the fluxes of heat at the air-sea interface.The purpose of this validation is to see how much of the change is attributed to our redistributed heat using the minimum transformation method.In the zonal mean (Fig. A2b) the difference between the simulated heat content change and the inferred redistributed heat (which is precisely the inferred added heat) has an RMS of 0.4 TW (8lat) 21 .

(iii) Scenario 3
In this scenario there is both an explicit added-heat signal in the model data and the model redistributes heat in response to both natural variability and the imposed warming.Despite the inclusion of a nonzero global mean net surface heat flux in FAFMIP redistributed heat (as described above), it is instructive to see how well our material and redistributed heat estimates compare to the directly simulated added and redistributed heat variables.In the zonal mean (Fig. A2c) the difference between both the simulated FAFMIP added heat content and the inferred material heat content change and between the simulated FAFMIP redistributed heat and our water mass based redistributed heat, has an RMS of 2.4 TW (8lat) 21 .We emphasize that this difference should not necessarily be directly attributed to an inaccuracy in our method considering the differing meanings of redistributed heat between the model simulations and our method.In broad terms, we consider the stated differences between directly simulated and inferred changes to be acceptable.We made no attempt to tune method parameters to optimize correspondence with the simulated variables, but this could be pursued in the future.

b. Parameter sensitivity
Here we test the sensitivity of the results, in particular the zonally integrated added heat, to parameter choices within the minimum transformation method.The two choices were (i) the choice of relative penalty on temperature versus salinity changes (i.e., parameter a) and (ii) the number of water masses in T-S space used to represent the early and late ocean states.We discuss sensitivity to these choices here.
The reference case for a is the ratio of a constant haline contraction coefficient [b 0 5 7.55 3 10 24 kg (g kg 21 ) 21 m 23 ] to a constant thermal expansion coefficient [a 0 5 1.76 3 10 24 kg K 21 m 23 ; i.e., a 0 5 b 0 /a 0 5 4.3 K (g kg 21 ) 21 ].This choice implies that a transformation by 1 g kg 21 in Absolute Salinity is penalized equivalently to a transformation of 4.3 K in temperature.A larger a will cause the method to favor transformation along the S axis, and a smaller a will favor transformation along the T axis.We test the method in three cases-a 5 a 0 , a 0 /2, and 2a 0 (Fig. A3a)-and find RMS differences of 0.3 TW (8lat) 21 between the reference case and the doubling and halving cases.
In terms of T-S resolution, our reference case has a minimum bin size of 0.2 g kg 21 and 0.4 K. Using the quadtree, the grid is refined until either this resolution is achieved or the volume within a particular bin falls below 62 3 10 12 m 3 .We test the sensitivity of this choice by both refining and coarsening the resolution by a factor of 2 in both the salinity and temperature dimensions and reducing the FIG.A2.(a) Zonally integrated simulated added heat (solid red) and inferred material heat content change (dashed red) based on the minimum transformation method for years 41-46 of the FAFheat experiment, comparing the simulation with and without added heat.(b) Zonally integrated simulated heat content change (solid blue) and inferred redistributed heat (dashed blue) based on our minimum transformation method, comparing years 35-40 and 41-46 of the piControl experiment.(c) Zonally integrated simulated added heat (solid red) and redistributed heat (solid blue) in the FAFheat experiment and inferred material heat content change (dashed red) and redistributed heat (dashed blue) based on our minimum transformation method applied to the model data.volume threshold by a factor of 4 also.Decreasing the resolution induces an RMS change in estimated zonally averaged OHC of 0.5 TW (8lat) 21 , and increasing the resolution induces an RMS change of 0.4 TW (8lat) 21 (Fig. A3b).

c. Robustness of the twenty-first-century trend
To quantify the sensitivity of our trend results to the time period chosen and the specific observations made and mapped in that period, we carry out a bootstrap calculation.Our aim here is not to determine how accurate our trend is but rather to determine how representative it is of time period as a whole or whether specific years strongly influence the result.
We chose to subsample the data by including and excluding entire years from the analysis.Six years are used for the early  and late (2012-17) periods of our analysis of EN4.We therefore considered all possible permutations of the numbers 1-6 and reran our analysis of EN4, subsampling the years corresponding to those six numbers.For example, in the case [1,3,3,4,5,6] the early-period data were replaced with the years 2006, 2008 and 2008 repeated, 2009, 2010, and 2011 and the late-period data were replaced with 2012, 2014 and 2014 repeated, 2015, 2016, and 2017.There are 46 656 uniquely ordered permutations of the numbers 1-6 when repetition is permitted.Since the calculation is insensitive to the order of the six years for either the early or late period, in practice we only need to consider the 462 unique permutations (ignoring order) and weight each by its frequency in the larger set of ordered permutations.
Figure 5 shows the mean and Fig. A4 shows the standard deviation of the bootstrap ensemble; 62 standard deviations of the spread in estimates of zonally averaged heat content change are shown in Fig. 6.Since these error estimates are generally larger than our other parameter sensitivity tests, we use them as our formal uncertainties throughout the main text.

Comparison with Atlantic Meridional Heat Transport
Trend at 26°N We compare our estimate of the contribution of redistribution to anomalous meridional heat transport north of 268N in the Atlantic (MHT Redist ; Fig. 6c) with conventional meridional heat transport (MHT) data reported by Bryden et al. (2020) (Table B1).Our (b) Zonally integrated inferred material heat content change for cases in which the T 2 S bins are shrunk using quadtree until they either contain a volume of seawater less than 62 3 10 12 m 3 or have a bin size of 0.48C by 0.2 g kg 21 (black), cases in which the minimum volume is 15.5 3 10 12 m 3 and the minimum bin size is 0.28C by 0.1 g kg 21 (blue), and cases in which the minimum volume is 248 3 10 12 m 3 and the minimum bin size is 0.88C by 0.4 g kg 21 (red).We have considered the difference in OHC between two 6-yr periods (2006-11 vs 2012-17)

FIG. 1 .
FIG. 1. Portrait of changing ocean water masses: (a) inventory of ocean volume in Conservative Temperature vs Absolute Salinity coordinates (mean of 2006-17 inclusive) and (b) change in water mass volume between the early half and late half of the period divided by the six years (Sv).According to water mass theory, changes in air-sea heat and freshwater fluxes and/or changes in rates of diffusion are required for these changes to occur.

FIG. 3 .
FIG. 3. Gray lines show Conservative Temperature T and Absolute Salinity S bounds of each water mass (or ''bin'') generated using a quadtree for each geographical region.The average T and S of the water found within each bin are shown by the location of each marker, and the volume is represented by the color scale (log 10 m 3 ).Inventories and mean T and S values represent the entire period (2006-17 inclusive).Inset panels show masks associated with each geographical region.

d
FIG. 5. Heterogeneous pattern of total and redistributed heat content change contrast against robust material heat content change: (a) change in depth-integrated ocean heat content between 2006-11 and 2012-17 inclusive, (b) inferred redistributed heat, and (c) inferred material heat content change based on changing water masses for the same period.Regions where the magnitude of the signal is less significant (less than 2 standard deviations of a bootstrap ensemble) are stippled.

FIG. 6 .
FIG. 6. Material heat content change is accumulating in the tropics and subtropics, whereas existing heat is being redistributed southward.(a) Total heat content change (gray), redistribution contribution (blue), and material contribution (red).(b) Contributions to material heat content change from the Indian (green), Pacific (orange), and Atlantic (yellow) Oceans.(c) Meridional heat transport due to redistribution in the Southern Ocean (blue), Atlantic (cyan), and Indian plus Pacific Oceans (magenta).Shaded areas represent 62 standard deviations of a bootstrap ensemble.
FIG. A1.(a) Directly simulated added heat by the FAFheat experiment averaged over years 41-46 of the experiment.(b) Inferred added heat when the same FAFheat data are first homogenized in water masses (bins in temperature-salinity coordinates) and then remapped into the locations of those water masses over the same period.(c) Comparison of the zonal integration of the two quantities shown in (a) and (b).
FIG. A3. (a)Zonally integrated inferred material heat content change for cases in which the parameter a is set at a reference value of a 0 5a 0 /b 0 5 4.3 K (g kg 21 ) 21 (black) and then reduced (red) and increased (blue) by a factor of 2. (b) Zonally integrated inferred material heat content change for cases in which the T 2 S bins are shrunk using quadtree until they either contain a volume of seawater less than 62 3 10 12 m 3 or have a bin size of 0.48C by 0.2 g kg 21 (black), cases in which the minimum volume is 15.5 3 10 12 m 3 and the minimum bin size is 0.28C by 0.1 g kg 21 (blue), and cases in which the minimum volume is 248 3 10 12 m 3 and the minimum bin size is 0.88C by 0.4 g kg 21 (red).
FIG. A4.(a) One standard deviation of the heat content change inferred on the basis of subsampling early (2006-11) and late (2012-17) years of the EN4 dataset.Also shown are 1 standard deviation of the ensemble of inferred (b) material heat content change and (c) redistributed heat on the basis of our minimum transformation method applied to the same subsampled data as in (a).

TABLE 1 .
Material, redistribution, and total contributions to heat content change by ocean basin in terawatts and area as fraction of global ocean area.Heat content change estimates are based on differences between the periods 2006-11 and 2012-17 inclusive.Uncertainties are 62 standard deviations.The Southern Ocean is defined as the entire ocean south of 328S.The South Pacific, South Atlantic, and Indian Ocean estimates exclude the ocean south of 328S.The North Atlantic is split into a region south of and a region north of 448N.The latter includes the Arctic Ocean.

TABLE A1 .
Summary of data used for three validation scenarios: T ref and S ref are the temperatures and salinities from the piControl experiment, respectively; T added is the added-heat variable; T redist is the redistributed heat variable from the FAFheat experiment; and S heat is the salinity variable from the FAFheat experiment.The numbers in parentheses are the experiment years chosen [e.g., T ref(41)(42)(43)(44)(45)(46)is temperature from years 41 to 46 of the piControl experiment].
. Hence our OHC change and MHT Redist are related through " ð t 0 1Dt where t 0 is midnight 31 December 2012 and Dt is 6 In practice we have averages of MHT covering April-March (see Table B1); we approximate Eq. (B2) using 6-yr running means of MHT anomalies and then averaging these between 2009/10 and 2014/15.Our uncertainties are 62 times the standard deviation of the 6-yr running means.

TABLE B1 .
Atlantic meridional heat transport (MHT; PW) at 268N(Bryden et al. 2020), MHT anomaly relative to 2006-17, and 6-yr running mean of the MHT anomaly.The mean of 6-yr running means is relevant to the difference in OHC between2006-11 and  2012-17.