<figure><img src="image_1.png" /></figure>

HAL
open science

# Turbulent isopycnal mixing dominates thermohaline transformations of intermediate ocean waters

Florian Sévellec, Nicolas Kolodziejczyk, Esther Portela

► To cite this version:

Florian Sévellec, Nicolas Kolodziejczyk, Esther Portela. Turbulent isopycnal mixing dominates thermohaline transformations of intermediate ocean waters. Nature Communications, 2025, 16 (1), pp.9838. 10.1038/s41467-025-64806-z. hal-05406508

HAL Id: hal-05406508

https://hal.univ-brest.fr/hal-05406508v1

Submitted on 9 Dec 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Turbulent isopycnal mixing dominates thermohaline transformations of intermediate ocean waters

Received: 18 July 2024

Accepted: 28 September 2025

Published online: 07 November 2025

<figure><img src="image_4.png" /></figure>

Check for updates

Florian Sévellec, Nicolas Kolodziejczyk & Esther Portela

The ocean, because of its ability to store heat and as the biggest mediator of water exchanges on Earth, modulates the climate from years to centuries. By combining Water Mass transformation and the Temporal-Residual-Mean frameworks, together with a direct estimation of turbulent diffusivity (computed through an artificial intelligence method), we determine temperature and salinity transformations along density surfaces. Because isopycnal transformations barely alter the density field, they remain hidden from the ocean dynamics. Here we unveil these hidden, isopycnal transformations from in situ observations at 1000-m depth. We show that these isopycnal transformations are 5 times as large as their diapycnal counterpart and are occurring for ~ 67% in the Southern Ocean. Our estimations stress the role of the spatially-varying and anisotropic characteristics of isopycnal diffusion, often not considered in ocean models. This has major consequences for our understanding of the role of the ocean in response to anthropogenic climate changes.

Ocean circulation transports a vast amount of water over the globe (up to 150 Sv, where $1 \text{ Sv} = 10^6 \text{ m}^3\text{s}^{-1}$, in the Antarctic Circumpolar Current, for instance)$^1$. This large-scale ocean circulation is also associated with transport of heat (up to 1.5 PW northward at $20^\circ\text{N})^2$ and freshwater (around 0.5 Sv southward at $35^\circ\text{N}$ in the North Atlantic)$^3$. The former is a fundamental component of the climate equilibrium and adjustment by transporting the excess of heat from the tropics to the poles. Whereas the latter indirectly reflects the geographical pattern of evaporation, precipitation, river discharge, and ice melting, all taking part in the water cycle. In particular, the ocean is the most important region of precipitation and evaporation (i.e., 75% and 85%, respectively), making it particularly relevant for the Earth water cycle$^4$.

Beyond this equilibrated view, the Earth is experiencing changes since the pre-industrial era, with a significant amount of anthropogenically-induced greenhouse gases released into the atmosphere$^5$. As a consequence, the surface atmospheric temperature is increasing$^5$, despite a significant mitigation by the ocean, with an ocean absorption of about 90% of the Earth's heat excess$^{5-7}$. Together with this warming, the water cycle is also affected$^8$. This is reflected by changes in the evaporation and precipitation patterns and amplitudes$^8$. These changes, either directly (when over the ocean) or

indirectly through runoffs (when over land), are imprinted in the ocean freshwater transport$^{9,10}$ through, for instance, modification of the ocean salinity distribution$^{11}$. Hence, the ocean heat and freshwater transports are central actors of the Earth equilibrium and climate changes.

Mesoscale eddies are ubiquitous worldwide in the ocean$^{12}$ and, unlike the mean transport, they contribute to the heat and freshwater transport in two-fold. Firstly, equivalently to the laminar component, they move ocean water properties (i.e., its temperature and salinity) in a conserving manner within the ocean. Secondly, eddies also induce a transformation of these properties through turbulent diffusion$^{13}$ (e.g., mixing cold and warm waters or salty and fresh ones). Transformations of seawater properties can occur either through interactions with the atmosphere at the ocean surface or through turbulent diffusion in the ocean interior. Whereas the former is easier to estimate$^{8,14}$, the latter has always been more difficult to assess and is often deduced using inverse methodology$^{15-18}$. Assessing the role of turbulent eddy diffusion for heat and freshwater transports and for the transformation of seawater properties is a timely question of ocean and climate sciences, as demonstrated by numerous numerical model studies$^{19-21}$.
There has been a large range of studies on Water Mass Transformations<sup>22-28</sup>. These studies often assess the transformations in terms of temperature and salinity changes, without specifically separating the respective role of isopycnal and diapycnal transformations. In parallel, the role of diapycnal diffusion on temperature and salinity transformations in the ocean interior has been well established<sup>16,29</sup>. More recently, studies have suggested the key role of isopycnal processes on heat and freshwater uptakes<sup>30,31</sup>, opening the possibility of an important role of the isopycnal diffusion. This has been further confirmed as potentially important in the North Atlantic for the maintenance of the mean temperature and salinity structure<sup>32</sup>. However, these studies rely on either numerical models or inverse methods, both intrinsically depending on a chosen closure scheme (and sometimes imposed parameters). Hence, it is key to use direct observations to assess whether isopycnal or diapycnal diffusion is dominant in the ocean interior. Indeed, in ocean science, where reproducible controlled experiments are not possible, the use of direct observations is the only way to validate or invalidate hypotheses derived from these numerical or inverse modeling approaches.

In this work, we use direct observations to assess the changes and transformations of temperature and salinity induced by mesoscale eddy turbulence, which acts as an equivalent isopycnal turbulent diffusion. We show that, at 1000 m depth, these temperature and salinity changes mainly occurred in the Southern Ocean and along the Gulf Stream. These isopycnal changes induce temperature and salinity transformations significantly more important than their diapycnal equivalents. We also conclude that the spatially variable and anisotropic nature of the isopycnal turbulent diffusion is keys for an accurate estimation of the temperature and salinity transformations.

## Results

The major difficulty of obtaining robust, direct estimations of heat and freshwater transformations by eddies comes from the need of knowing turbulent diffusive coefficients independently of temperature and salinity fields (to avoid the dependence of the results to a closure scheme)<sup>33,34</sup>. Here, to circumvent this difficulty, we have used the recently estimated diffusivity coefficients<sup>35</sup> derived from Argo float deep displacements<sup>36</sup>. This dataset has been obtained using a class of data-driven Artificial Intelligence methods: the Analog Method<sup>37,38</sup>. There, the spatially varying and anisotropic diffusivity coefficients were estimated through local tracer dispersion induced by a turbulent velocity field reconstructed from deep Argo displacements. This dataset provides the full diffusivity tensor, explicitly acknowledging the anisotropy, rather than a single diffusivity coefficient, as most commonly done<sup>33,34,39</sup>. These resulting diffusivity coefficients are strong where the flow is strong (Fig. 1d, e), but with a slight spatial offset with the core of the currents<sup>40,41</sup>, hence with the largest density gradients. This could be rationalized by the existence of Kinetic Energy transfers<sup>42</sup>. These coefficients (Fig. 1d, e) together with temperature and salinity fields from gridded ISAS17 product<sup>43</sup> (Optimal Interpolation<sup>44</sup> from Argo profiles, Fig. 1a, b) allow us to have a direct estimation of isopycnal heat and freshwater transformations and uncover an unobserved facet of deep ocean circulation.

## Effect of isopycnal turbulent diffusion

In the present study, we assess the global effect of the turbulent diffusion on temperature and salinity fields. However, because of observational constraints (i.e., the typical parking depth of Argo float deep displacements), our study focuses at 1000 m depth. This depth is particularly relevant for estimating heat and freshwater transports and seawater property transformations for two reasons. (1) It is close to the depth of the maximum of the Atlantic Meridional Overturning Circulation (depth of no-motion for this circulation, as estimated at 26.5°N in the North Atlantic by RAPID)<sup>45</sup>. (2) It is the depth with the largest

contrast in salinity<sup>46,47</sup>, but also includes signatures of the coldest waters<sup>47</sup> (Fig. 1).

The transformation of seawater properties works hand in hand with the concept of Water Masses<sup>48</sup>. Water Masses are volumes of water with consistent seawater properties (Fig. 1f) that are set through ocean-atmosphere exchange at the ocean surface, where water masses originate. They then spread within the ocean through advection and mixing and are transformed by the latter. At 1000 m depth, we will consider three intermediate water masses: Red Sea-Persian Gulf Intermediate Water (RSPGIW) formed in the Red Sea and Persian Gulf<sup>19</sup>, the Mediterranean Water<sup>50</sup> (MW) reaching the Atlantic through the Strait of Gibraltar, and the AntArctic Intermediate Water<sup>51</sup> (AAIW) intersecting the 1000-m depth along the Polar Front<sup>52</sup>. We will also consider more dense Water Masses, namely the North Atlantic Deep Water (NADW), its transformed southern extension, the Circumpolar Deep Water<sup>53</sup> (CDW), and the Antarctic Bottom Water<sup>54</sup> (AABW). Despite not being their depth of maximum intrusion, these Water Masses are present at 1000 m depth in the high latitudes of the North Atlantic and of the Southern Ocean (Fig. 1f).

To estimate the impact of the turbulent diffusion on the temperature and salinity of these water masses, we follow the Temporal-Residual-Mean framework (see Methods – Isopycnal Turbulent Transformation Framework)<sup>13</sup>. This framework disentangles advective and mixing effects when eddy turbulence is averaged. This allows for an accurate attribution of eddy-induced transport and transformation. Also, this framework acknowledges that transformations act within isopycnal layers as epineutral flux (Fig. 2). To this purpose averaging are done along-isopycnal, so that fictitious diapycnal transport, due to incorrect estimates of the isopycnal slopes, can be limited<sup>55</sup>. This transformation-induced process has almost no signature on the density field (except for the weak cabbeling and thermobaricity components, although they remain crucial for the density transformations<sup>56</sup>), nor pressure, nor circulation. As a consequence, there is no impact on Potential Energy, and it can act passively, as long as Eddy Kinetic Energy remains to maintain the turbulent diffusion. In this sense, for the ocean dynamics, it remains hidden along the isopycnals (Fig. 1c). However, it still contributes to heat and freshwater transports from cold to warm and from salty to fresh regions, respectively. This process can be described in the form of temperature and salinity rates of change or of water-mass transformation rates (see Methods – Water Mass Transformation Rates). The rates of change of a material property (temperature or salinity) measure the amount of warming/cooling or salinification/freshening induced by the isopycnal turbulent diffusion; whereas the transformation accounts for the amount of seawater that have their properties transformed by this process.

## Rates of temperature and salinity change

As mentioned above, the isopycnal diffusion of temperature and salinity can be assessed as the induced rates of change in terms of temperature and salinity [see Methods – (10)]. High magnitudes of the rates of change are due to high diffusivity coefficients, sharp isopycnal temperature or salinity gradients, or a combination of both. At 1000 m depth, we find that they are of the order of kelvins per year or tenth of (g kg⁻¹) per year [or a few 10⁻⁸ K s⁻¹ and 10⁻⁹ g kg⁻¹ s⁻¹] for temperature and salinity, respectively. These rates of change are mainly located in the Southern Ocean, along the Antarctic Circumpolar Current, and, to a lesser extent, in the North West of the North Atlantic (Figs. 3a and 4a). Indeed, the Southern Ocean accounts for 67% of the global rates of change for temperature and salinity [as measured by the R² of the transformations described in (11)], despite representing only 31% of the global ocean area. When including the North Atlantic it reaches 95% for temperature and salinity rates of change. A local maximum of these rates of change can be seen in the Zapiola gyre, which corresponds to a local maximum of diffusivity coefficients (Fig. 1d, e). There is also significant rates of change along the Agulhas return current and its
<figure><img src="image_5.png" /></figure>

**Fig. 1 | Climatological tracer fields and horizontal diffusivity at 1000 m depth.**
a Conservative Temperature, (b) Absolute Salinity, and (c) neutral density at 1000 m depth from the climatological-annual mean of the gridded ISAS17 product⁴³. (d) Primary and (e) secondary diffusivity, as well as their respective directions, estimated from deep Argo displacements at 1000 m depth from the ANDRO dataset³⁶. (f) Six Water Masses are identified at 1000 m depth¹⁹ and denoted through color patches. These Water Masses are the Mediterranean Water (MW), the

Red Sea–Persian Gulf Intermediate Water (RSPGIW), the AntArctic Intermediate Water (AAIW), the North Atlantic Deep Water (NADW), the Circumpolar Deep Water (CDW), and the AntArctic Bottom Water (AABW). Four key regions are also highlighted and used in the rest of the study: North Atlantic (NA), Atlantic sector of the Southern Ocean (SO-A), Indian sector of the Southern Ocean (SO-I), and Pacific sector of the Southern Ocean (SO-P).

merging with the Antarctic Circumpolar Current, which also corresponds to a local maximum of diffusivity coefficients. Indeed, between South Africa and the Kerguelen plateau, there is a warming and salinification of the relatively cold and salty AAIW (Figs. 3a and 4a), a cooling and freshening of the relatively warm and fresh AAIW, and a warming and salinification of the warm and salty residual of RSGPIW (trapped in the Agulhas system⁵⁷, Fig. 1b).

Unlike the geographical map that reveals where the rates of change are located, the TS-diagram (i.e., representation with temperature-salinity axes) reveals which water masses are subjected to the isopycnal diffusion (i.e., change along isopycnal surfaces). This can be represented globally (insets in Figs. 3a and 4a) or displayed for each relevant basin (following the division of Fig. 1f), which allows a simpler relation with local Water Masses (Fig. 3b-e and 4b-e). In the North Atlantic, despite the weak isopycnal turbulent diffusion coefficients in the subtropical gyre interior (Fig. 1d, e), temperature and salinity rates of change within the MW (along the 27.6-27.8 kg m⁻³ density class) are visible with a slight cooling and a freshening of MW (Fig. 3b). A more

important cooling and freshening is also visible at the edge between MW and NADW. In the Southern Ocean, the geographical separation is quite revealing. It suggests different water masses at play for each of the three sectors (Atlantic, Indian, and Pacific). In the Atlantic sector of the Southern Ocean, we observe a cooling and freshening of CDW and a warming and salinification of AAIW and AABW (Figs. 3c and 4c). In the Indian and Pacific sectors, the rates of change correspond to a warming and salinification of the AABW, a cooling and freshening of the CDW, a warming and salinification of the coldest and freshest portion of the AAIW, and a cooling and freshening of the warmest and saltiest portion of the AAIW (Figs. 3d, e and 4d, e). In the Indian sector, there is also a warming and salinification of the warmer portion of the AAIW, which derives from RSPGIW⁵⁷.

## Water mass transformations

To further quantify these effects, the isopycnal turbulent diffusion can also be expressed as Water Masses Transformation (see Methods – Water Mass Transformation Rates). As discussed earlier, since the
<figure><img src="image_6.png" /></figure>

**Fig. 2 | Distribution of temperature and salinity along 25°W and potential transformation induced by the oceanic isopycnal turbulence.** a Conservative Temperature and (b) Absolute Salinity along 25°W from the climatological-annual mean of the gridded ISAS17 product¹³. The light gray lines indicate neutral density contours of 26.5, 27.25, 27.85, and 28.02 kg m⁻³. Because the existence of Temperature and Salinity (TS) gradients along isopycnal, the oceanic isopycnal turbulent diffusion (gray arrow) can transform TS-properties (black spiral). This is especially true at 1000 m depth (denoted by the horizontal black line), where

isopycnal temperature and salinity gradients are important given the range of different Water Masses. At 25 °W we have: AntArctic Bottom Water (AABW), Ant-Arctic Intermediate Water (AAIW), Mediterranean Water (MW), North Atlantic Deep Water (NADW), and Circumpolar Deep Water (CDW). As suggested by the schematic (a and b), it is possible to compute the rates of change, following (10), and the corresponding water mass transformation, following (11), associated to a tracer (e.g., temperature or salinity).

turbulent eddy-induced diffusion acts along isopycnals¹³ it cannot modify density (except for weak nonlinear effects, which affect density and hence the circulation)⁵⁶. Hence, we will focus on transformations induced by temperature and salinity changes. We use the traditional Water Mass Transformation framework⁴⁸,⁵⁸, often defined across density layers. Here, we compute transformation rates induced by temperature and salinity changes along the isopycnals [see Methods – (II)]. They are computed both globally and for each region defined earlier (Fig. 1f). These two transformations are evidence of heat and freshwater transport by isopycnal turbulent diffusion.

We find density transformations of the order of ±2Sv (i.e., density-class standard deviation) for density changes induced by either temperature and salinity, respectively (colored – red and blue – bars in Fig. 5a, c, respectively). This is important as these transformations correspond to a 70 m layer only (i.e., the typical density layer thickness due to the stratification at our studied depth, see Methods – Water Mass Transformation Rates and Supplementary Fig. 1). For comparison, the buoyancy transformations are of the order of a few tens of Sverdrup when computed over the full ocean depth⁵⁶. We rationalize the importance of the temperature and salinity transformations over the density ones because of their strong thermohaline compensating effect on density (i.e., the so called spiciness). The Water Masses located at our studied depth exhibit particularly intense spiciness contrasts¹⁸, although it is not necessarily the case at other depths. More broadly, at the ocean surface, warming regions are often regions of net evaporation, whereas cooling regions are often regions of net precipitation. This means that there is a non-negligible part of the ocean circulation that transports heat and freshwater from cold and fresh regions to warm and salty regions without transporting density.

For reference, we compare these isopycnal transformations to the diapycnal ones, which induce a net density transformations⁵⁶. These diapycnal transformations are computed using reconstructed vertical diffusion coefficients⁵⁹,⁶⁰ (see Methods –

Alternative Diffusivity Coefficients.), which show the expected typical control by topography (Supplementary Fig. 2). To compare the isopycnal and diapycnal transformations, we measure their amplitude through their standard deviation along the density coordinate. We find that the temperature and salinity isopycnal transformations are 3 and 5 times as large as the diapycnal ones (colored-bars vs cyan line in Fig. 5a, c), respectively. Note that this result contrasts with density transformations for which, as expected, diapycnal transformations is the dominant factor over isopycnal ones, except for a few specific Water Masses (i.e., Subantarctic Mode Water, AAIW, and CDW), for which cabbeling and thermobaricity are important factors⁵⁶.

These dominant isopycnal transformations reveal a mean density gain and loss of +1.5 Sv and -1.9 Sv, respectively. The density transformations correspond to three regimes in density coordinates (Fig. 5): a warming and salinification of the lightest water, a cooling and freshening of density classes between 27.35 and 28 kg m⁻³, and a warming and salinification of the densest water. As expected within an isopycnal framework, transformations with respect to neutral density classes show that density gain and loss by cooling and warming are associated to density loss and gain by freshening and salinification, respectively, preserving the density field (Fig. 5e). Looking at water classes and regional separation, we find warming and salinification of the lightest class of AAIW and cooling of the densest class of AAIW. This mainly occurs in the Indian sector of the Southern Ocean. On the other hand, transformations associated with CDW and AABW exhibit a warming and salinification, which mainly occurred in the Atlantic sector of the Southern Ocean, and to a second degree in the Pacific sector for CDW. A transformation peak is also visible in the North Atlantic sector. It corresponds to a cooling compensated by a freshening at the edge between MW and NADW. Similarly, in the same sector, the full range of MW is also transformed through cooling and freshening, whereas a slight warming and salinification also occur for the NADW.
<figure><img src="image_7.png" /></figure>

**Fig. 3** | **Temperature rates of change at 1000 m depth.** Rates of Change in (a) Conservative Temperature at 1000 m depth induced by the isopycnal turbulence. The Inset shows the same rates of change but as a Temperature-Salinity diagram, with neutral density in contours (in kg m⁻³ with intervals of 0.2 kg m⁻³). **b-e** same as the inset but for regional estimations of the rates of change for Conservative Temperature. The four selected regions are North Atlantic (NA), Atlantic sector of Southern Ocean (SO-A), Indian sector of Southern Ocean (SO-I), and Pacific sector

of Southern Ocean (SO-P), as separated by the lines and represented in Fig. 1f. In the scatter plots, the values are ordered in absolute rates of change (potentially hiding weaker rates of change of the same Temperature-Salinity class). The color patches show properties of selected Water Masses³⁰: Mediterranean Water (MW), Red Sea-Persian Gulf Intermediate Water (RSPGIW), AntArctic Intermediate Water (AAIW), North Atlantic Deep Water (NADW), Circumpolar Deep Water (CDW), and AntArctic Bottom Water (AABW).

To better understand the role and to assess the importance of spatially varying anisotropic diffusivity coefficients³⁶, we have tested a range of less complex diffusivity coefficients (see Methods – Alternative Diffusivity Coefficients). To quantitatively test the misestimations of these simplest diffusivity coefficients, we used their relative error in comparison to the most complete version [see Methods – (12)]. The simplest one was a spatially-uniform and isotropic isopycnal diffusivity coefficient (light gray bars in Fig. 5a, b). It mainly preserves the signs of the transformations, but underestimates them by up to 53% for both temperature and salinity. Nonetheless, this approximation warm and salinified AAIW from 27.4 kg m⁻³ to 27.5 kg m⁻³, whereas it was cool and freshened with the most complete version of the diffusivity. Using spatially-varying and isotropic isopycnal diffusivity coefficients (dark gray bars in Fig. 5a, b) better

conserves the shape of the transformations but overestimates by up to 48% both of these temperature and salinity transformations.

## Robustness of the isopycnal diffusivity dataset

To test the robustness of our results, we have recomputed the temperature and salinity transformations using two alternative diffusivity coefficient datasets (see Methods – Alternative Diffusivity Coefficients). The first dataset, that we named Groeskamp, estimates diffusivity using salinity observations and assuming the validity of the mixing length theory³⁴. The second dataset, that we named Roach, uses the same deep Argo float displacements than us, but a different methodology: a pairwise dispersion³⁹.

Groeskamp and Roach diffusivities (Fig. 6a, b, respectively) have the same overall pattern as our reference diffusivities (Fig. 1d, e):
<figure><img src="image_8.png" /></figure>

**Fig. 4 | Salinity rates of change at 1000 m depth.** Rates of Change in (a) Absolute Salinity at 1000 m depth induced by the isopycnal turbulence. The inset show the same rates of change but as a Temperature-Salinity diagram, with neutral density in contours (in kg m⁻³ with intervals of 0.2 kg m⁻³). b-e same as the inset but for regional estimations of the rates of change for Absolute Salinity. The four selected regions are North Atlantic (NA), Atlantic sector of Southern Ocean (SO-A), Indian sector of Southern Ocean (SO-I), and Pacific sector of Southern Ocean (SO-P), as

separated by the lines and represented in Fig. 1f. In the scatter plots, the values are ordered in absolute rates of change (potentially hiding weaker rates of change of the same Temperature-Salinity class). The color patches show properties of selected Water Masses¹⁹: Mediterranean Water (MW), Red Sea-Persian Gulf Intermediate Water (RSPGIW), AntArctic Intermediate Water (AAIW), North Atlantic Deep Water (NADW), Circumpolar Deep Water (CDW), and AntArctic Bottom Water (AABW).

intensification in western boundary currents and within the Antarctic Circumpolar Current. Because Groeskamp and Roach diffusivity coefficients correspond solely to the cross-flow diffusion, they are more comparable to our secondary diffusivity coefficients (Fig. 1e), which are by construction, weaker than the primary ones. Groeskamp and Roach diffusivity coefficients range between 2000 m²s⁻¹ and 5000 m²s⁻¹ in active regions (as our secondary diffusivity coefficients). In comparison with the two other datasets, Roach diffusivity coefficients are less spatially smooth.

Because the Groeskamp and Roach datasets provide only a single local diffusivity coefficient, we restrict the comparison to the spatially varying and isotropic framework. Computing the temperature and salinity transformations with the different diffusivity coefficient datasets shows the good consistency of the results (Fig. 6c and d). When

using our spatially varying and isotropic diffusivity coefficients, we find stronger transformations than those obtained with Groeskamp and Roach diffusivity coefficients. However, the transformations are qualitatively well represented. As expected, the agreement is better when applying the spatially varying and isotropic diffusivity coefficients solely based on the secondary diffusivity coefficients. Here, the agreement is qualitative and quantitative. This stresses the robustness of our described temperature and salinity transformations with respect to the methodology used to infer the diffusivity coefficients.

## Discussion

The ocean, through its absorption of 91% of the heat imbalance⁵ and as the largest region of evaporation and precipitation on Earth³, is a key component of our changing climate. At first glance, the ocean is filled
<figure><img src="image_9.png" /><figcaption>**Fig. 5** | **Temperature and Salinity water mass transformations at 1000 m depth for density classes.** Water mass transformations in terms of density changes for (a and b) Conservative Temperature and for (b and d) Absolute Salinity, and (c and d) their sum for density classes. Positive or negative values show warming or salinification and cooling or freshening, respectively. **a**, **c**, and **e** Red and blue bars show positive and negative transformations for the anisotropic and spatially-varying diffusivity, respectively; dark and light gray bars show transformations for the isotropic and spatially-varying and for the isotropic and globally-uniform diffusivity (set as the local-mean and global-mean of the anisotropic diffusivity), respectively. The cyan lines are the transformations induced by diapycnal diffusion. **e**, **d**, and **f** The transformations are split in geographical sectors: global (black lines),</figcaption></figure>

North Atlantic (NA, orange line), Atlantic sector of Southern Ocean (SO-A, blue line), Indian sector of Southern Ocean (SO-I, yellow line), and Pacific sector of Southern Ocean (SO-P, green line), as represented in Fig. 1f. The curves have been smoothed using a moving average window of 0.06 kg m⁻³ and re-interpolated on regular 0.012 kg m⁻³ grid with a spline interpolation. The horizontal color patches show neutral density classes of selected Water Masses⁴⁹: Mediterranean Water (MW), Red Sea–Persian Gulf Intermediate Water (RSPGIW), AntArctic Intermediate Water (AAIW), North Atlantic Deep Water (NADW), Circumpolar Deep Water (CDW), and AntArctic Bottom Water (AABW). **e** and **f** Non-zero values for the sum of isopycnal temperature and salinity transformations correspond to the computational error.

by fast, intense mesoscale eddy features, ubiquitous worldwide¹². Hence, the role of these eddies and their associated turbulence in balancing the ocean circulation and its heat and freshwater transports remains a vivid question of current climate and ocean research.

This question is even more intriguing when acknowledging that, since eddy turbulence corresponds to a flow along isopycnals, its induced mixing cannot modify density (except through weak – but potentially locally important⁵⁶ – nonlinear effects), hence the circulation. However, they interact with temperature and salinity fields, transforming seawater properties along their way¹³, without a trace on the circulation, but still acting on the heat and freshwater transports. From an energetic point of view, by transporting temperature and salinity properties, without disturbing the density (because of the compensating effect), these transformations do not impact the potential energy, remaining hidden for ocean dynamics and acting passively, as long as Eddy Kinetic Energy exists for maintaining the eddy turbulence.

The first result of the present study is the observational estimation of temperature and salinity transformations at 1000 m depth. By using observations of the ocean temperature and salinity together with a direct observational estimate of anisotropic diffusivity coefficients induced by mesoscale eddy turbulence³⁶, we have unveiled the existence of this intense, hidden isopycnal transformations fundamental for the transport of heat and freshwater in the ocean and, therefore, for setting the global temperature and salinity seawater properties. These isopycnal transformations are 3 and 5 times as large as the more classical diapycnal transformations in terms of temperature and salinity transformations, respectively. In particular, in the North Atlantic we find an important dominance of isopycnal transformations over diapycnal transformations at our studied, 1000 m depth, consistently with an independent recent observational study⁶¹. We have shown that these dominant isopycnal transformations occur for ~67% in the Southern Ocean (despite representing only 31% of the global ocean area) and is particularly relevant for AAIW, CDW, and AABW.
<figure><img src="image_10.png" /><figcaption>Fig. 6 | Other estimations of isopycnal diffusivity and induced Temperature and Salinity water mass transformations at 1000 m depth for density class. (a and b) Isopycnal diffusivity coefficients at 1000 m estimated with alternative methods³⁴,³⁹ (denoted as Groeskamp and Roach, respectively). Water mass transformations in terms of density changes for (c) Conservative Temperature and for (d) Absolute Salinity for density classes. Blue lines show transformations for the isotropic and spatially varying diffusivity. Solid and dashed blue lines show results for the average between primary and secondary diffusivities and secondary</figcaption></figure>

diffusivities, respectively. The former is strictly identical to the dark gray bars of Fig. 5. Orange and red lines show transformations for the Groeskamp and Roach diffusivities, respectively. The curves have been smoothed using a moving average window of 0.06 Kg m⁻³ re-interpolated on a regular 0.012 Kg m⁻³ grid with a spline interpolation. The horizontal color patches show neutral density classes of selected Water Masses⁴⁰: Mediterranean Water (MW), Red Sea–Persian Gulf Intermediate Water (RSPGIW), AntArctic Intermediate Water (AAIW), North Atlantic Deep Water (NADW), Circumpolar Deep Water (CDW), and AntArctic Bottom Water (AABW).

The second result is the demonstration of the importance of acknowledging the full complexity of the isopycnal diffusivity tensor. We have shown that considering the spatially varying and anisotropic diffusivity coefficients is key for the accurate transformation estimations. Simpler diffusivity coefficients (as often used in parameterization of laminar/eddy-less ocean component used in some current climate models)⁶² lead to misestimation of the seawater properties transformation (by up to 53%), thus impacting oceanic heat and freshwater transports. Although because our study is limited to a single-depth layer, we could not quantitatively estimate these impacts. Hence, this should be tested in state-of-the-art climate models to further identify the relevance and significance of our finding. In particular, the role of spatially varying and anisotropic diffusivity coefficients in modulating and redistributing ocean heat storage in climate simulations should be evaluated. Given the projected changes of ocean mesoscale turbulence over the 21st century⁶³, hence of the isopycnal diffusivity coefficients³⁶, the future response of our identified intense hidden isopycnal transformations is in question.

Whereas, the role of diathermal and diahaline water-mass transformations has been long established¹⁶,²⁹, there exists in the recent literature evidence of the importance of the isopycnal direction over the diapycnal one for these temperature and salinity transformations¹⁸,³²,⁶¹. Our direct estimates suggest the dominance of isopycnal diffusion over diapycnal diffusion at 1000 m depth. We hypothesize that this is because we fully acknowledged the anisotropic nature of the isopycnal diffusion within our framework, unlike in inverse or numerical modeling studies. This dominance is particularly important since isopycnal transformations are key for heat uptake in the context of Global Warming³⁰.

Because of observational constraints³⁶, our analysis is restricted at 1000 m depth. Hence, our results are only representative of a horizontal slice of the ocean, potentially masking geographical variations linked to the local layer thickness (relatively more important in the Southern Ocean and in the subpolar North Atlantic). However, this depth is particularly interesting because it is close to the depth of the maximum meridional overturning transport in the Atlantic (i.e., the depth of no-motion for this circulation), and it is where some of the most saline waters reside and where the coldest waters also exist. This makes this depth particularly relevant for heat and freshwater transports. However, given the intensity of the circulation and transformation identified in our study, it is important to acknowledge it at other depths. In particular, diapycnal mixing could become more prevalent near the seabed, mitigating the key role of isopycnal mixing for temperature and salinity transformations. Hence, building a robust 3-dimensional picture, fully constrained by observations, of the isopycnal turbulent transformations of the global ocean will be a direction for future work.

## Methods

### Isopycnal turbulent transformation framework

Starting from the advection equation for a tracer (which could be temperature or salinity) away from forcing (mainly occurring at the ocean surface) and for high Péclet number (Pe ≫ 1, measuring the advection to diffusion ratio), we have:

$$D_t C = 0, \qquad (1)$$

where C is the tracer concentration and $D_t$ is the material derivative (= $\partial_t + u\partial_x + v\partial_y + w\partial_z$, where t is time, x, y, and z are the zonal,
meridional, and vertical coordinates, and u, v, and w are the zonal, meridional, and vertical velocities).

This equation can be modified to show the evolution of the large-spatial scale or long-time scale by splitting each term in an Eulerian mean (.) and its associated anomaly (.'). This reads:

$$ \overline{D_t \hat{C}} = - \partial_x \overline{u' \hat{C}'} - \partial_y \overline{v' \hat{C}'} - \partial_z \overline{w' \hat{C}'} \quad (2) $$

where $\overline{D_t}$ is the mean material derivative ($=\partial_t + \bar{u}\partial_x + \bar{v}\partial_y + \bar{w}\partial_z$).

This shows that several mean ($C$, $\bar{u}$, $\bar{v}$, and $\bar{w}$) and turbulent ($\overline{u' \hat{C}'}$, $\overline{v' \hat{C}'}$, and $\overline{w' \hat{C}'}$) terms contribute to the maintenance of the mean state. In this context, we can redistribute the terms$^{13,64}$ so that the evolution of the tracer becomes:

$$ \hat{D_t \hat{C}} = d_C \quad (3) $$

where $\hat{C}$ is the modified tracer, $\hat{D_t}$ is the modified material derivative [$=\partial_t + \hat{u}\partial_x + \hat{v}\partial_y + \hat{w}\partial_z$, and $\hat{u}=\bar{u}+\bar{u}$, $\hat{v}=\bar{v}+\bar{v}$, and $\hat{w}=\bar{w}+\bar{w}$, where $\bar{u}$, $\bar{v}$, and $\bar{w}$ are the zonal, meridional, and vertical eddy-induced velocities, respectively], and $d_C$ is the transformation term. This expression suggests that the modified tracer is advected by both the mean and the eddy-induced velocity$^{65}$ and transformed. In this Temporal-Residual-Mean framework, the eddy-induced advection and the transformation are consequences of the skew-symmetric and symmetric components of the turbulent flux, respectively. This framework is particularly useful to accurately separate and attribute the advection and mixing effects of averaged eddy turbulence. The advection effect only displaces tracer properties, whereas the mixing effect changes tracer properties. Hence, this framework avoids the misinterpretation of eddy-induced transport by horizontal turbulence as a transformation.

In the ocean interior (i.e., below the ocean mixed layer), the transformation term can take the form$^{56}$ of

$$ d_C = d_C^{\text{iso}} + d_C^{\text{dia}} \quad (4) $$

where $d_C^{\text{iso}}$ and $d_C^{\text{dia}}$ are the transformation along the isopycnal and the diapycnal direction. Here, by being in the ocean interior (i.e., below mixed layer depth), we have neglected the potential effect of horizontal transformation (important within the mixed layer). In the rest of the study, we will focus on the isopycnal component of these fluxes and compare it to the diapycnal one. Note that isopycnal heat and freshwater transformations act along a constant density surface, thus, by construction, they avoid density gradients. Although the diapycnal component, induced by small-scale mixing processes$^{66}$, will be used as a reference. It is worth noting that the use horizontal diffusion instead of isopycnal diffusion outside of the mixed layer is fundamentally flawed$^{55}$. With this approximation, we would find transformations that are completely missing the correct warming/cooling and freshening/salinification. Also, as expected for such flawed approximation, the compensation between temperature and salinity changes are missed. This would result in unrealistic transformations.

Following a downgradient flux closure$^{13,67}$ (i.e., Fick's law), the density-weighted isopycnal transformation reads:

$$ d_C^{\text{iso}} = \partial_z \gamma \nabla_n \cdot \left( \mathbf{K} \cdot \frac{\nabla_n \hat{C}}{\partial_z \gamma} \right) \quad (5) $$

where $\nabla_n$ refers to the gradient along the neutral density ($\gamma$) and $\mathbf{K}$ is the, possibly spatially-varying, anisotropic diffusivity tensor acting along the neutral density (or anisotropic epineutral diffusivity tensor). Here, the gradient is computed using the projected

nonorthogonal approach$^{68,69}$. The diffusivity tensor is expressed as:

$$ \mathbf{K} = \begin{pmatrix} \kappa_{\bar{x}\bar{x}} & \kappa_{\bar{x}\bar{y}} \\ \kappa_{\bar{x}\bar{y}} & \kappa_{\bar{y}\bar{y}} \end{pmatrix}, \quad (6) $$

where $\bar{x}$ and $\bar{y}$ are the zonal and meridional direction along neutral density, respectively, and $\kappa_{\bar{x}\bar{x}}$, $\kappa_{\bar{x}\bar{y}}$, and $\kappa_{\bar{y}\bar{y}}$ are the zonal, meridional, and cross diffusivities along the neutral density, respectively. Hence, the isopycnal transformations for a spatially varying and anisotropic diffusivity tensor read:

$$ d_C^{\text{iso}} = \partial_z \gamma \partial_{\bar{x}} \left( \frac{\kappa_{\bar{x}\bar{x}}}{\partial_z \gamma} \partial_{\bar{x}} \hat{C} \right) + \partial_z \gamma \partial_{\bar{x}} \left( \frac{\kappa_{\bar{x}\bar{y}}}{\partial_z \gamma} \partial_{\bar{y}} \hat{C} \right) + \partial_z \gamma \partial_{\bar{y}} \left( \frac{\kappa_{\bar{x}\bar{y}}}{\partial_z \gamma} \partial_{\bar{x}} \hat{C} \right) + \partial_z \gamma \partial_{\bar{y}} \left( \frac{\kappa_{\bar{y}\bar{y}}}{\partial_z \gamma} \partial_{\bar{y}} \hat{C} \right). \quad (7) $$

This expression shows that the diffusivity is dependent on 3 scalars that measure both the intensity of the diffusivity and its orientation. Indeed, since the off-diagonal terms are strictly identical in (6), it simply induces a rotation of diffusivity, which will act independently along its two orthogonal natural directions (directions after rotation). Here, it becomes apparent that the use of a single coefficient cannot fully represent the complexity of these isopycnal transformations, both in terms of intensity and directionality.

## Computation of epineutral flux

One last difficulty is to estimate the gradient along neutral density$^{69}$ (or epineutral gradient) as expressed in (7): computing the $\nabla_n$ (or $\partial_{\bar{x}}$ and $\partial_{\bar{y}}$). As an example, but it would be equivalent for the meridional direction, for the zonal direction, we have:

$$ \begin{align} \partial_{\bar{x}} \sigma &\approx \frac{\Delta \sigma |_{\gamma}}{\Delta \bar{x}} \approx \frac{\Delta \sigma |_{z} + \Delta \sigma |_{x}}{\sqrt{\Delta x^2 + \Delta z^2}} \approx \frac{\partial_x \sigma \Delta x + \partial_z \sigma \Delta z}{\Delta x \sqrt{1 + \frac{\Delta z^2}{\Delta x^2}}}, \tag{8} \\ &\approx \partial_x \sigma + s_x \partial_z \sigma - \frac{1}{2} s_x^2 \partial_x \sigma - \frac{1}{2} s_x^3 \partial_z \sigma + \dots, \nonumber \\ &\approx \partial_x \sigma + s_x \partial_z \sigma, \nonumber \end{align} $$

where $\sigma$ is any scalar quantity (e.g., Temperature or Salinity), $s_x = \delta z / \Delta x$ is the zonal component of the isopycnal slope, $\Delta \sigma |_{\gamma, x, z}$ is the difference in $\sigma$ at constant $\gamma$, $x$, and $z$, respectively, $\Delta \bar{x}$ is the distance at constant $\gamma$ along the zonal direction, and $\Delta x$ and $\delta z$ are the zonal and vertical distances, respectively, covered by the neutral density surface.

In the last calculation, a few approximations have been done: (i) $\Delta \bar{x} \approx \sqrt{\Delta x^2 + \delta z^2}$, assuming that the epineutral displacement can be geometrically approximated by the zonal and vertical displacement; (ii) $\delta z \ll \Delta x$, assuming that the vertical displacement of the epineutral is weak compared to the zonal one (i.e., small slope approximation$^{65}$), allowing a Taylor expansion of the denominator; and (iii) $\partial_x \sigma$ and $\partial_z \sigma$ are constant over $\Delta x$ and $\delta z$, respectively. Equivalently, we have $\partial_y \sigma \approx \partial_y \sigma + s_y \partial_z \sigma$, where $s_y$ is the meridional component of the isopycnal slope. Note that there exist more complex but more accurate computation methods of the neutral slopes and gradients, which are key when the slopes become important$^{70,71}$. However, in our application, the slopes remain weak, with a maximum slope magnitude of 0.08 and a median value of 0.0001.

The slope is defined through the horizontal gradients of Absolute Salinity ($S_A$) and Conservative Temperature ($\Theta$), by setting $\partial_{\bar{x}} \gamma = 0$ and $\partial_{\bar{y}} \gamma = 0$, such as:

$$ \left( \frac{s_x}{s_y} \right) = - \frac{\alpha \nabla \Theta - \beta \nabla S_A}{\alpha \partial_z \Theta - \beta \partial_z S_A}, \quad (9) $$
where $\nabla$ is the horizontal gradient, and $\alpha$ and $\beta$ are the thermal expansion and haline contraction coefficients, respectively. $S^2$ is the magnitude of the slope measured through the components inner product (= $s_x^2 + s_y^2$),

This epineutral formulation for the gradients allows us to compute an approximation of the isopycnal diffusion, described in (5) and (7), for any tracer $\mathcal{C}$ as a function of Cartesian coordinates. This reads:

$$d_{\mathcal{C}}^{\text{iso}} \approx \partial_z \gamma \left\{ \begin{aligned} & \partial_x \left[ \frac{\kappa_{xx}}{\partial_z \gamma} (\partial_x \hat{C} + s_x \partial_z \hat{C}) + \frac{\kappa_{yy}}{\partial_z \gamma} (\partial_y \hat{C} + s_y \partial_z \hat{C}) \right] \\ & + s_x \partial_z \left[ \frac{\kappa_{xx}}{\partial_z \gamma} (\partial_x \hat{C} + s_x \partial_z \hat{C}) + \frac{\kappa_{yy}}{\partial_z \gamma} (\partial_y \hat{C} + s_y \partial_z \hat{C}) \right] \\ & + \partial_y \left[ \frac{\kappa_{yy}}{\partial_z \gamma} (\partial_x \hat{C} + s_x \partial_z \hat{C}) + \frac{\kappa_{yy}}{\partial_z \gamma} (\partial_y \hat{C} + s_y \partial_z \hat{C}) \right] \\ & + s_y \partial_z \left[ \frac{\kappa_{yy}}{\partial_z \gamma} (\partial_x \hat{C} + s_x \partial_z \hat{C}) + \frac{\kappa_{yy}}{\partial_z \gamma} (\partial_y \hat{C} + s_y \partial_z \hat{C}) \right] \end{aligned} \right\} \quad (10)$$

This is the formulation, which includes anisotropic, isopycnal diffusivities, that is used in our study.

## Water mass transformation rates

Building on (10), the action of the isopycnal turbulent term could be further estimated in terms of transformation rates$^{56,72}$ as:

$$\Psi_{\Theta}^{\#} = \partial_{\gamma} \iint_{\Sigma(\gamma' \le \gamma)} -\alpha d_{\Theta}^{\text{iso}} d\Sigma \Delta z, \quad (11a)$$

$$\Psi_{S_A}^{\#} = \partial_{\gamma} \iint_{\Sigma(\gamma' \le \gamma)} \beta d_{S_A}^{\text{iso}} d\Sigma \Delta z, \quad (11b)$$

where $\Sigma(\gamma' \le \gamma)$ is the surface of a tracer class below the threshold value $\gamma$ at 1000 m depth. Here $\Psi_{\Theta}^{\#}$ and $\Psi_{S_A}^{\#}$ represent neutral density transformations through temperature or salinity changes over a layer of thickness $\Delta z$, respectively. These transformations are expressed in Sv, where $1 \text{ Sv} = 10^6 \text{ m}^2 \text{ s}^{-1}$. In practice, the transformations are estimated using neutral density increment (we chose: $\Delta y = 0.03 \text{ kg m}^{-3}$). For consistency, the layer thickness can be estimated though the scaling: $\Delta z \approx \Delta y(\partial_z \gamma)^{-1}$ (which derives from $dy = \partial_x \gamma dx + \partial_y \gamma dy + \partial_z \gamma dz$). Computation over the globe at 1000 m depth show values from 10 m to more than 3000 m (Supplementary Fig. 1a). In particular, the thickness is significantly more important in the Southern Ocean (e.g., the Weddell Sea) and in the subpolar North Atlantic (e.g., Labrador Sea) than in the rest of the globe. However, extrapolating our computation over more than 3000 m is inconsistent with our knowledge of vertical variations of diffusivity coefficient: existence of a steering level$^{40,41}$ and variations with Eddy Kinetic Energy$^{36}$. Hence, we have chosen to restrict the layer thickness amplitude. We have used a uniform layer thickness based on the median value: $\Delta z = 70 \text{ m}$ (Supplementary Fig. 1b), regardless of the local stratification $(\partial_z \gamma)$. Hence, our results, despite being in Sverdrup, only reflect a slice of the ocean at 1000 m depth.

These transformations mathematically expressed a probability density function of density changes induced by cooling/warming or freshening/salinification per neutral density classes, respectively. They physically correspond to the volume of water at a given density that is gaining/loosing density by cooling/warming or by salinification/ freshening, respectively.

These transformation rates ($\Psi_{\Theta}^{\#}$ and $\Psi_{S_A}^{\#}$) can be interpreted as the component of the total transformation induced by isopycnal turbulence. Other components exist$^{23}$, such as the one induced by the vertical turbulence$^{73}$ [$w\mathcal{C}'$, visible in (2)] and the small-scale diffusion [neglected in (1), by assuming $\text{Pe} \gg 1$].

To summarize, (11) shows the existence of transformations when both isopycnal oceanic turbulence and isopycnal tracer gradients are present (Fig. 2). In this framework, one can compute transformations associated to temperature and salinity (Fig. 5). The sum of the two

transformations with respect to neutral density classes can also be defined as: $\Psi_{\Theta}^{\#} + \Psi_{S_A}^{\#}$. This sum is expected to be zero and, therefore, reflects the computation errors (Fig. 5e and f).

## Tracer fields, diffusivity data, and computation method

One can compute a tracer transformation and transformation rates by knowing two properties:

1. $\hat{C}$ and $\partial_z \hat{C}$ (as well as $\partial_z \gamma$) at all locations $[\forall (x, y)]$ along the surface (defined by $z = z_0$, where $z_0$ is a constant depth). This can be obtained after assuming $\hat{C} \approx \hat{C}$ from relevant reanalyzes at any depth, but primarily above the maximum Argo depth sampling of ~2000 m for robustness. In this study, $\mathcal{C}$ can be Conservative Temperature or Absolute Salinity (Fig. 1a, b);

2. $\mathbf{K}$ at all locations $[\forall (x, y)]$ along the same surface (defined by $z = z_0$, where $z_0$ is a constant depth), which can be estimated by Argo deep displacements at their parking depth$^{36}$ (Fig. 1d-f).

Hence, the main constraint occurs through $\mathbf{K}$ which is only determined at Argo float parking depth (mainly located at ~1000 m) and determines the focused depth of our study.

We have used the two-dimensional gridded direct estimate of the turbulent anisotropic, isopycnal diffusivity coefficients derived from Argo deep displacements$^{36}$, which are available$^{35}$ at https://doi.org/10.17882/91335, and the global three-dimensional ISAS17 gridded temperature and salinity data derived from objective analysis of in situ observations$^{43}$ and available at https://doi.org/10.17882/52367.

The spatially varying and anisotropic diffusivity coefficients were obtained by reconstructing an analog turbulent field from deep Argo displacement. This reconstructed turbulent field was used to determine the local dispersion of a locally-injected spike of passive tracer, the rates of dispersion estimating the local diffusivity coefficients along the zonal, meridional, and cross directions. This allows for the full determination of the diffusivity tensor, rather than a single coefficient as often done$^{33,34,39}$.

In practice, after conversion of the data to Conservative Temperature and Absolute Salinity, as well as the computation of neutral density, all the computations are done using $3^\circ \times 3^\circ$ cells (i.e., the native cell of the diffusivity coefficients) on $1^\circ \times 1^\circ$ grid (i.e., the native grid of the diffusivity coefficients). Data are vertically averaged over 40 m bins, and vertical levels (needed for the computation of vertical gradients) are separated by 20 m. All conversions to Conservative Temperature and Absolute Salinity, as well as computation of neutral density, are done using the Gibbs SeaWater Oceanographic Toolbox of TEOS-10 code.

For the water masses framework, we have identified six Water Masses located at the depth of the study$^{49}$. They are defined through their temperature and salinity range. They are the Mediterranean Water (MW), the Red Sea-Persian Gulf Intermediate Water (RSPGIW), the Antarctic Intermediate Water (AAIW), the North Atlantic Water (NADW), the Circumpolar Deep Water (CDW), and the Antarctic Bottom Water (AABW).

## Alternative Diffusivity Coefficients

To test the importance of using spatially-varying and anisotropic isopycnal diffusivity coefficients, we have tested a range of less complex diffusivity coefficients:

1. Spatially-uniform and isotropic isopycnal diffusivity coefficients: the diffusivity coefficients are set to $\kappa_{xx}^{\text{Uni \& Iso}} = \kappa_{yy}^{\text{Uni \& Iso}} = (\langle \kappa_{xx} \rangle + \langle \kappa_{yy} \rangle)/2$ and $\kappa_{xy}^{\text{Uni \& Iso}} = 0$, where $\langle \rangle$ reflects the global horizontal average. These coefficients are applied to (5).

2. Spatially-varying and isotropic isopycnal diffusivity coefficients: the diffusivity coefficients are set to $\kappa_{xx}^{\text{Var \& Iso}} = \kappa_{yy}^{\text{Var \& Iso}} = (\langle \kappa_{xx} \rangle + \langle \kappa_{yy} \rangle)/2$ and $\kappa_{xy}^{\text{Var \& Iso}} = 0$. These coefficients are applied to (5).
To compare transformations using the different versions of the diffusivity, we used a relative error (RE) using the most complete version as the reference. This reads:

$$RE = \text{sign} \left[ \log \left( \frac{\Psi_{\text{comp}}^{\#}}{\Psi_{\text{ref}}^{\#}} \right) \right] \frac{\text{var} \left( \Psi_{\text{comp}}^{\#} - \Psi_{\text{ref}}^{\#} \right)}{\text{var} \left( \Psi_{\text{ref}}^{\#} \right)}, \quad (12)$$

where $\Psi_{\text{ref}}^{\#}$ is the reference transformation (i.e., most complete diffusivity), $\Psi_{\text{comp}}^{\#}$ is the transformation used for comparison (i.e., less complete diffusivity), and the log term allows to estimate the sign of the mis-estimation (positive and negative for over- and under-estimation, respectively).

In addition to ensuring the robustness of our results, we have tested two other diffusivity coefficient datasets. They have been computed using different methods and/or different observations. Hence, we will also test a modified version of our dataset, which is more aligned with the constraints of these two datasets and allows a fairer comparison between the three datasets.

a. The first dataset will be named Groeskamp and corresponds to spatially-varying coefficients³⁴. This dataset is available at: https://figshare.com/articles/dataset/Groeskamp et al_2020--mixing-diffusivities/12554555. The Groeskamp diffusivity coefficients ($\kappa_{\text{Groeskamp}}$) are obtained using a mixing length closure theory to extract cross-flow diffusivity coefficients from salinity ocean observations. These coefficients cover the full depth of the ocean but we extracted their values at 1000 m (Fig. 6a). Then, to be applied in our framework, we set this single spatially-varying coefficient as spatially-varying but isotropic coefficients: $\kappa_{\text{Groeskamp}}^{\text{Var}} = \kappa_{\text{Groeskamp}}^{\text{Iso}} = \kappa_{\text{Groeskamp}}$ and $\kappa_{\text{Groeskamp}}^{\text{Var}} = \kappa_{\text{Groeskamp}}^{\text{Iso}} = 0$. Then we apply this set of coefficients to (5).

b. The second dataset will be named Roach and also correspond to spatially varying diffusivity coefficients³⁹. This dataset is available at: https://github.com/croachutas/Isopycnal_Diffusivity. The Roach coefficients ($\kappa_{\text{Roach}}$) use a pairwise distance from deep Argo float displacements to infer the cross-flow diffusivity coefficients. By construction, this dataset is restricted to our 1000 m studied depth (Fig. 6b). As for the Groeskamp diffusivity coefficients, we set this single spatially varying coefficient as isotropic spatially varying coefficients: $\kappa_{\text{Roach}}^{\text{Var}} = \kappa_{\text{Roach}}^{\text{Iso}} = \kappa_{\text{Roach}}$ and $\kappa_{\text{Roach}}^{\text{Var}} = \kappa_{\text{Roach}}^{\text{Iso}} = 0$. This is applied to (5).

c. Since these two sets of diffusion coefficients were computed as cross-flow diffusion coefficients, consistently, we have used our secondary diffusivity coefficient ($\kappa_s$) to have a fairer comparison with our diffusivity dataset³⁹. The secondary diffusivity coefficient is the second eigenvalue of the diffusivity tensor (K) described in (6). It is, by construction, weaker than the primary eigenvalue and is more consistent with the cross-flow diffusivity coefficient (i.e., avoiding subgrid stirring and shearing). In principle, this secondary diffusivity is associated with a particular isopycnal direction (i.e., the second eigenvector of the diffusivity tensor). However, here for fair comparison with Groeskamp and Roach diffusivity coefficients described above, we have also applied it as a spatially varying and isotropic isopycnal diffusivity coefficient. Hence, it is set as: $\kappa_{\text{Var}}^{\text{Var & Iso}} = \kappa_{\text{Var}}^{\text{Var & Iso}} = \kappa_s$ and $\kappa_{\text{Iso}}^{\text{Var & Iso}} = 0$. These coefficients are applied to (5).

Finally, for comparison purposes, we have computed the isopycnal transformation as a vertical diffusion: $d_{\hat{C}}^{\text{dia}} = \partial_z(\kappa_s \partial_z \hat{C})$, where $\kappa_s$ is the vertical diffusivity coefficient. We have used spatially varying diffusivity coefficients⁵⁹,⁶⁰ available at https://doi.org/10.17882/73082, and have extracted the corresponding values for 1000 m depth (Supplementary Fig. 2). This shows the expected topographic control of the vertical diffusion coefficients.

## Data availability

The global gridded temperature and salinity product is ISASI7 (release key: 86441) and available at: https://doi.org/10.17882/52367. Turbulent anisotropic, isopycnal diffusivity coefficients derived from ANDRO Argo deep displacements are available at https://doi.org/10.17882/91335. ISAS and ANDRO Argo floats displacements products are made available by SNO Argo France at LOPS Laboratory (supported by UBO/CNRS/Ifremer/IRD) and IUEM Observatory (OSU IUEM/CNRS/INSU), and was funded by Ifremer, Coriolis, SOERE-CTDO2, SNO Argo France, and CPER Bretagne Obsocean. Groeskamp isopycnal diffusivity coefficients are available at https://figshare.com/articles/dataset/Groeskamp et al_2020--mixing_diffusivities/12554555. Roach isopycnal diffusivity coefficients are available at https://github.com/croachutas/Isopycnal_Diffusivity. Vertical diffusivity coefficients are available at: https://doi.org/10.17882/73082.

## Code availability

The Gibbs SeaWater Oceanographic Toolbox of TEOS-10 package used for the computation of Neutral Density, Absolute Salinity, and Conservative Temperature is available at http://www.teos-10.org/pubs/gsw/html/gsw_contents.html. The M_Map package for MATLAB, used for global mapping of data, is available online at https://www-old.eoas.ubc.ca/-rich/map.html. The code developed for the analyses and for producing the figures are available from the corresponding author on request.

## References

1. Xu, X., Chassignet, E. P., Firing, Y. L. & Donohue, K. Antarctic circumpolar current transport through drake passage: What can we learn from comparing high-resolution model results to observations? J. Geophys. Res. Oceans 125, e2020JC016365 (2020).

2. Forget, G. & Ferreira, D. Global ocean heat transport dominated by heat export from the tropical pacific. Nat. Geosci. 12, 351–354 (2019).

3. Mignac, D., Ferreira, D. & Haines, K. Decoupled freshwater transport and meridional overturning in the south atlantic. Geophys. Res. Lett. 56, 2178–2186 (2019).

4. Stephens, G. L. et al. Earth’s water reservoirs in a changing climate. Proc. Math Phys. Eng. Sci. 476, 20190458 (2020).

5. IPCC. Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change (Cambridge University Press, 2021).

6. von Schuckmann, K. et al. Heat stored in the Earth system 1960–2020: where does the energy go? Earth System Science Data 15, 1675–1709 (2023).

7. Li, Z., England, M. H. & Groeskamp, S. Recent acceleration in global ocean heat accumulation by Mode and Intermediate Waters. Nat. Commun. 14, 6888 (2023).

8. Hegerl, G. C. et al. Challenges in quantifying changes in the global water cycle. Bull. Amer. Meteor. Soc. 96, 1097–1115 (2015).

9. Skiris, N. et al. Salinity changes in the world ocean since 1950 in relation to changing surface freshwater fluxes. Clim. Dyn. 43, 709–736 (2014).

10. Zika, J. D. et al. Maintenance and broadening of the ocean’s salinity distribution by the water cycle. J. Clim. 28, 9550–9560 (2015).

11. Durack, P. J. & Wijffels, S. E. Fifty-year trends in global ocean salinities and their relationship to broad-scale warming. J. Clim. 23, 4342–4362 (2010).

12. Chelton, D. B., Schlax, M. G., Samelson, R. M. & de Szoeke, R. A. Global observations of large oceanic eddies. Geophys. Res. Lett. 34, L15606 (2007).

13. McDougall, T. J. & McIntosh, P. C. The temporal-residual-mean velocity. part II: Isopycnal interpretation and the tracer and momentum equations. J. Phys. Oceanogr. 31, 1222–1246 (2001).
14. Tomita, H., Kunio, K., Masahisa, K. & Tsutomu, H. Advances in the estimation of global surface net heat flux based on satellite observation: J-ofuro3 v1.1. *Front. Mar. Sci.* **8**, 612361 (2021).

15. Zika, J. D., McDougall, T. J. & Sloyan, B. M. A tracer-contour inverse method for estimating ocean circulation and mixing. *J. Phys. Oceanogr.* **40**, 26–47 (2010).

16. Groeskamp, S., Sloyan, B. M., Zika, J. D. & McDougall, T. J. Mixing inferred from an ocean climatology and surface fluxes. *J. Phys. Oceanogr.* **47**, 667–687 (2017).

17. Hautala, S. L. The abyssal and deep circulation of the northeast pacific basin. *Prog. Oceanogr.* **160**, 68–82 (2018).

18. Portela, E., Kolodziejczyk, N., Maes, C. & Thierry, V. Interior water-mass variability in the southern hemisphere oceans during the last decade. *J. Phys. Oceanogr.* **50**, 361–381 (2020).

19. Sévellec, F. & Fedorov, A. V. Stability of the Atlantic meridional overturning circulation in a zonally-averaged ocean model: the effects of freshwater flux, wind stress, and diapycnal diffusion. *Deep Sea Res.* **58**, 1927–1943 (2011).

20. Jones, C. S. & Abernathey, R. P. Isopycnal mixing controls deep ocean ventilation. *Geophys. Res. Lett.* **46**, 13144–13151 (2019).

21. Holmes, R. M., Groeskamp, S., Stewart, K. D. & McDougall, T. J. Sensitivity of a coarse-resolution global ocean model to a spatially variable neutral diffusivity. *J. Adv. Model. Earth Syst.* **14**, e2021MS002914 (2022).

22. Speer, K. & Tziperman, E. Rates of water mass formation in the North Atlantic Ocean. *J. Phys. Oceanogr.* **22**, 93–104 (1992).

23. Zika, J. D., England, M. H. & Sijp, W. P. The ocean circulation in thermohaline coordinates. *J. Phys. Oceanogr.* **42**, 708–724 (2012).

24. Döös, K., Nilsson, J., Nycander, J., Brodeau, L. & Ballarotta, M. The world thermohaline circulation. *J. Phys. Oceanogr.* **42**, 1445–1460 (2012).

25. Groeskamp, S., Zika, J. D., McDougall, T. J., Sloyan, B. M. & Laliberté, F. The representation of ocean circulation and variability in thermodynamic coordinates. *J. Phys. Oceanogr.* **44**, 1735–1750 (2014).

26. Pemberton, P., Nilsson, J., Hieronymus, M. & Markus Meier, H. E. Arctic ocean water mass transformation in S–T coordinate. *J. Phys. Oceanogr.* **45**, 1025–1050 (2015).

27. Evans, D. G., Zika, J. D., Naveira Garabato, A. C. & Nurser, A. J. G. The cold transit of Southern Ocean upwelling. *Geophys. Res. Lett.* **45**, 3386–3395 (2018).

28. Zika, J. D., Gregory, J. M., McDonagh, E. L., Marzocchi, A. & Clément, L. Recent water mass changes reveal mechanism of ocean warming. *J. Phys. Oceanogr.* **34**, 3461–3479 (2021).

29. Hieronymus, M., Nilsson, J. & Nycander, J. Water mass transformation in salinity-temperature space. *J. Phys. Oceanogr.* **44**, 2547–2568 (2014).

30. Clement, L. et al. Mechanisms of ocean heat uptake along and across isopycnals. *J. Phys. Oceanogr.* **35**, 4885–4904 (2022).

31. Zou, S., Lozier, M. S., Li, F., Abernathey, R. & Jackson, L. Density-compensated overturning in the Labrador Sea. *Nat. Geosci.* **13**, 121–126 (2020).

32. Evans, D. G., Holliday, N. P., Bacon, S. & Le Bras, I. Mixing and air-sea buoyancy fluxes set the time-mean overturning circulation in the subpolar north atlantic and nordic seas. *Ocean Sci.* **19**, 745–768 (2023).

33. Cole, S. T., Wortham, C., Kunze, E. & Brechner Owens, W. Eddy stirring and horizontal diffusivity from Argo float observations: Geographic and depth variability. *Geophys. Res. Lett.* **42**, 3989–3997 (2015).

34. Groeskamp, S., LaCasce, J. H., McDougall, T. J. & Rogé, M. Full-depth global estimates of ocean mesoscale eddy mixing from observations and theory. *Geophys. Res. Lett.* **47**, e2020GL089425 (2020).

35. Sévellec, F., Colin de Verdière, A. & Kolodziejczyk, N. Deep horizontal turbulent diffusivity. *SEANOE https://doi.org/10.17882/91335* (2022).

36. Sévellec, F., Colin de Verdière, A. & Kolodziejczyk, N. Estimation of horizontal turbulent diffusivity from deep Argo float displacements. *J. Phys. Oceanogr.* **52**, 1509–1529 (2022).

37. Lorenz, E. N. Atmospheric predictability as revealed by naturally occurring analogues. *J. Atmos. Sci.* **26**, 636–646 (1965).

38. Ayet, A. & Tandeo, P. Nowcasting solar irradiance using an analog method and geostationary satellite images. *Solar Energy* **164**, 301–315 (2018).

39. Roach, C. J., Bawada, D. & Speer, K. Global observations of horizontal mixing from Argo float and surface drifter trajectories. *J. Geophys. Res. Oceans* **123**, 4560–4575 (2018).

40. Smith, K. S. & Marshall, J. Evidence for enhanced eddy mixing at middepth in the Southern Ocean. *J. Phys. Oceanogr.* **39**, 50–69 (2009).

41. Abernathey, R., Marshall, J. & Mazloff, M. Enhancement of mesoscale eddy stirring at steering levels in the Southern Ocean. *J. Phys. Oceanogr.* **40**, 170–184 (2010).

42. Sévellec, F., Colin de Verdière, A. & Kolodziejczyk, N. Global observations of deep ocean kinetic energy transfers. *J. Phys. Oceanogr.* **54**, 1633–1654 (2024).

43. Kolodziejczyk, N., Prigent-Mazella, A. & Gaillard, F. ISAS temperature, salinity, dissolved oxygen gridded fields. *SEANOE https://doi.org/10.17882/52367* (2023).

44. Bretherton, F. P., Davis, R. & Fandry, C. Technique for objective analysis and design of oceanographic experiment applied to MODE-73. *Deep Sea Res. Oceanogr. Abstr.* **23**, 559–582 (1976).

45. Stepanov, V. N., Iovino, D., Masina, S., Storto, A. & Cipollone, A. Methods of calculation of the atlantic meridional heat and volume transports from ocean models at 26.5°n. *J. Geophys. Res. Oceans* **121**, 1459–1475 (2016).

46. Talley, L. D. Freshwater transport estimates and the global overturning circulation: Shallow, deep and throughflow components. *Prog. Oceanogr.* **78**, 257–303 (2008).

47. Talley, L. D., Pickard, G. L., Emery, W. J. & Swift, J. H. Descriptive Physical Oceanography: An Introduction (sixth edition). (2011).

48. Walin, G. On the relation between sea-surface heat flow and thermal circulation in the ocean. *Tellus* **34**, 187–195 (1982).

49. Emery, W. J. *Encyclopedia of Atmospheric Sciences* 1556–1567 (Academic Press, Oxford, 2018).

50. Iorga, M. C. & Lozier, M. S. Signatures of Mediterranean outflow from North Atlantic climatology – 1. Salinity and density fields. *J. Geophys. Res.* **104**, 25985–26009 (1999).

51. Carter, L., McCave, I. N. & Williams, M. J. Circulation and water masses of the southern ocean: A review in anatractic climate evolution. *Dev. Earth Environ. Sci.* **8**, 85–114 (2009).

52. Orsi, A. H., Whitworth, T. & Nowlin, W. D. On the meridional extent and fronts of the Antarctic Circumpolar Current. *Deep Sea Res.* **42**, 641–673 (1995).

53. Defant, A. *Physical Oceanography* (Pergamon Press, 1961).

54. Emery, W. J. & Meincke, J. Global water masses: summary and review. *Oceanol. Acta* **9**, 383–391 (1986).

55. Veronis, G. The role of models in tracer studies, *Numerical models of the Ocean circulation*. *Natl. Acad. Sci.* **133**–146 (1975).

56. Groeskamp, S., Abernathey, R. P. & Klocker, A. Water mass transformation by cabbeling and thermobaricity. *Geophys. Res. Lett.* **43**, 10,835–10,845 (2016).

57. Roman, R. E. & Lutjeharms, J. R. E. Red Sea Intermediate Water at the Agulhas current termination. *Deep Sea Res.* **54**, 1329–1340 (2007).

58. Groeskamp, S. et al. The water mass transformation framework for ocean physics and biogeochemistry. *Annu. Rev. Mar. Sci.* **11**, 271–305 (2019).

59. de Lavergne, C. et al. A parameterization of local and remote tidal mixing. *J. Adv. Model. Earth Syst.* **12**, e2020MS002065 (2020).

60. de Lavergne, C. Global tidal mixing maps. *SEANOE https://doi.org/10.17882/73082* (2020).
61. Fernández Castro, B. et al. Isopycnal eddy stirring dominates thermohaline mixing in the upper subpolar north atlantic. *J. Geophys. Res.* **129**, e2023JC020817 (2024).

62. Eyring, V. et al. Overview of the coupled model intercomparison project phase 6 (CMIP6) experimental design and organization. *Geosci. Model Dev.* **9**, 1937–1958 (2016).

63. Beech, N. et al. Long-term evolution of ocean eddy activity in a warming world. *Nature Clim. Change* **12**, 910–917 (2022).

64. Young, W. R. An exact Thickness-Weighted Average formulation of the boussinesq equations. *J. Phys. Oceanogr.* **42**, 692–707 (2012).

65. Gent, P. R. & McWilliams, J. C. Isopycnal mixing in ocean circulation model. *J. Phys. Oceanogr.* **20**, 150–155 (1990).

66. de Lavergne, C., Madec, G., Le Sommer, J., Nurser, A. J. G. & Naveira Garabato, A. C. The impact of a variable mixing efficiency on the abyssal overturning. *J. Phys. Oceanogr.* **46**, 663–681 (2016).

67. Redi, M. H. Oceanic isopycnal mixing by coordinate rotation. *J. Phys. Oceanogr.* **12**, 1154–1158 (1982).

68. Starr, V. P. A quasi-lagrangian system of hydrodynamical equations. *J. Meteor.* **2**, 227–237 (1945).

69. McDougall, T. J., Groeskamp, S. & Griffies, S. M. On geometrical aspects of interior ocean mixing. *J. Phys. Oceanogr.* **44**, 2164–2175 (2014).

70. Groeskamp, S., Barker, P. M., McDougall, T. J., Abernathey, R. P. & Griffies, S. M. Venm: An algorithm to accurately calculate neutral slopes and gradients. *J. Adv. Model. Earth Syst.* **11**, 1917–1939 (2019).

71. Shao, A. E., Adcroft, A., Hallberg, R. & Griffies, S. M. A general-coordinate, nonlocal neutral diffusion operato. *J. Adv. Model. Earth Syst.* **12**, e2019MS00199 (2020).

72. Maze, G., Forget, G., Buckley, M. & Marshall, J. Using transformation and formation maps to study the role of air-sea heat fluxes in north atlantic eighteen degree water formation. *J. Phys. Oceanogr.* **39**, 1818–1835 (2009).

73. Holmes, R. M., Zika, J. D. & England, M. H. Diathermal heat transport in a global ocean model. *J. Phys. Oceanogr.* **49**, 141–161 (2019).

## Acknowledgements

This work was supported by ARVOR and OASIS projects from LEFE IMAGO and GMMC program, respectively, by the OceaniX project funded through the French ANR program, and by the ISblue project, Interdisciplinary graduate school for the blue planet (ANR-17-EURE-0015) co-funded by a grant from the French government under the program "Investissements d'Avenir". These data were collected and made freely available by the International Argo Program and the national programs that contribute to it (https://argo.ucsd.edu and https://www.ocean-ops.org). The Argo Program is part of the Global Ocean Observing System. F.S. is supported by the EERIE project (Grant Agreement 101081383) funded by the European Union. Views and opinions expressed are, however, those of the authors only and do not necessarily

reflect those of the European Union or the European Climate Infrastructure and Environment Executive Agency (CINEA). Neither the European Union nor the granting authority can be held responsible for them.

## Author contributions

All three authors (F.S., N.K., and E.P.) contributed equally to the experimental design and the writing of the manuscript. F.S. initiated the analysis and developed the data analysis.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-025-64806-z.

**Correspondence** and requests for materials should be addressed to Florian Sévellec.

**Peer review information** *Nature Communications* thanks Henri Drake and the other anonymous reviewer(s) for their contribution to the peer review of this work. A peer review file is available.

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

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

**Open Access** This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

© The Author(s) 2025