Turbulent isopycnal mixing dominates thermohaline transformations of intermediate ocean waters

Domain: earth_science
Published: 2025-11-07
Editorial outcome: major_revision

Abstract
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.

Introduction
Ocean circulation transports a vast amount of water over the globe (up to 150 Sv, where 1 Sv ≡ 10^6 m^3 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°N)2 and freshwater (around 0.5 Sv southward at 35°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 cycle4.

Beyond this equilibrated view, the Earth is experiencing changes since the pre-industrial era, with a significant amount of anthropogenically-induced greenhouse gazes released into the atmosphere5. As a consequence, the surface atmospheric temperature is increasing5, despite a significant mitigation by the ocean, with an ocean absorption of about 90% of the Earth’s heat excess5–7. Together with this warming, the water cycle is also affected8. These changes, either directly (when over the ocean) or indirectly through runoffs (when overland), are imprinted in the ocean freshwater transport9,10 through, for instance, modification of the ocean salinity distribution11. Hence, the ocean heat and freshwater transports are central actors of the Earth equilibrium and climate changes.

Mesoscale eddies are ubiquitous worldwide in the ocean12 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 diffusion13 (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 estimate8,14, the latter has always been more difficult to assess and is often deduced using inverse methodology15–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 studies19–21.

There has been a large range of studies on Water Mass Transformations22–28. 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 established16,29. More recently, studies have suggested the key role of isopycnal processes on heat and freshwater uptakes30,31, 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 structure32. 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)33,34. Here, to circumvent this difficulty, we have used the recently estimated diffusivity coefficients35 derived from Argo float deep displacements36. This dataset has been obtained using a class of data-driven Artificial Intelligence methods: the Analog Method37,38. 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 done33,34,39. 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 currents40,41, hence with the largest density gradients. This could be rationalized by the existence of Kinetic Energy transfers42. These coefficients (Fig. 1d, e) together with temperature and salinity fields from gridded ISAS17 product43 (Optimal Interpolation44 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)45. (2) It is the depth with the largest contrast in salinity46,47, but also includes signatures of the coldest waters47 (Fig. 1).

The transformation of seawater properties works hand in hand with the concept of Water Masses48. 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 Gulf49, the Mediterranean Water50 (MW) reaching the Atlantic through the Strait of Gibraltar, and the AntArctic Intermediate Water51 (AAIW) intersecting the 1000-m depth along the Polar Front52. We will also consider more dense Water Masses, namely the North Atlantic Deep Water (NADW), its transformed southern extension, the Circumpolar Deep Water53 (CDW), and the Antarctic Bottom Water54 (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)13. 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 limited55. 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 transformations56), 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^-1) per year [or a few 10^-8 K s^-1 and 10^-9 g kg^-1 s^-1] 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^2 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 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 system57, 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 isopycnals). 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^-3 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 RSPGIW57.

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 isopycnal temperature and salinity gradients are important given the range of different Water Masses. At 25°W we have: AntArctic Bottom Water (AABW), AntArctic 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).

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 factors56.

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^-3, 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.

To better understand the role and to assess the importance of spatially varying anisotropic diffusivity coefficients36, 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^-3 to 27.5 kg m^-3, 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 theory34. The second dataset, that we named Roach, uses the same deep Argo float displacements than us, but a different methodology: a pairwise dispersion39.

Groeskamp and Roach diffusivities (Fig. 6a, b, respectively) have the same overall pattern as our reference diffusivities (Fig. 1d, e): 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^2 s^-1 and 5000 m^2 s^-1 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.

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, (1)

where C is the tracer concentration and D_t is the material derivative (=∂_t + u∂_x + v∂_y + w∂_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-timescale by splitting each term in an Eulerian mean (¯) and its associated anomaly ('). This reads:

D_t ¯C = -∂_x ¯(u'C') - ∂_y ¯(v'C') - ∂_z ¯(w'C'), (2)

where D_t is the mean material derivative (=∂_t + ¯u∂_x + ¯v∂_y + ¯w∂_z). This shows that several mean (¯C, ¯u, ¯v, and ¯w) and turbulent (¯(u'C'), ¯(v'C'), and ¯(w'C')) terms contribute to the maintenance of the mean state.

In the Temporal-Residual-Mean framework, we can redistribute the terms13,64 so that the evolution of the mean tracer is:

D_b_t ˆC = d ˆC, (3)

where ˆC is the modified tracer, D_b_t is the modified material derivative [=∂_t + ˆu∂_x + ˆv∂_y + ˆw∂_z, and ˆu=¯u+˜u, ˆv=¯v+˜v, and ˆw=¯w+˜w, where ˜u, ˜v, and ˜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 velocity65 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 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 processes66, 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 flawed55. 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 closure13,67 (i.e., Fick’s law), the density-weighted isopycnal transformation reads:

d_i ˆC_so = ∂_z γ ∇_n · K · ∇_n (∂_z γ)^-1, (5)

where ∇_n refers to the gradient along the neutral density (γ) and K is the spatially-varying and anisotropic isopycnal diffusivity tensor. Here, the gradient is computed using the projected nonorthogonal approach68,69. The diffusivity tensor is expressed as:

K = [κ_x~x~ κ_x~y~; κ_x~y~ κ_y~y~], (6)

where x~ and y~ are the zonal and meridional direction along neutral density, respectively, and κ_x~x~, κ_y~y~, and κ_x~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_i ˆC_so = ∂_z γ [∂_x~ (κ_x~x~ ∂_x~ ˆC) + ∂_x~ (κ_x~y~ ∂_y~ ˆC) + ∂_y~ (κ_x~y~ ∂_x~ ˆC) + ∂_y~ (κ_y~y~ ∂_y~ ˆC)], (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 density69 (or epineutral gradient) as expressed in (7): computing the ∇_n (or ∂_x~ and ∂_y~). As an example, but it would be equivalent for the meridional direction, for the zonal direction, we have:

∂_x~ σ ≈ ∂_x σ + s_x ∂_z σ, (8)

where σ is any scalar quantity (e.g., Temperature or Salinity), s_x = δz/Δx is the zonal component of the isopycnal slope, Δσ|γ,x,z is the difference in σ at constant γ, x, and z, respectively, Δx~ is the distance at constant γ along the zonal direction, and Δx and δz are the zonal and vertical distances, respectively, covered by the neutral density surface. In this calculation, a few approximations have been done: (i) Δx~ ≈ Δx^2 + δz^2, assuming that the epineutral displacement can be geometrically approximated by the zonal and vertical displacement; (ii) δz ≪ Δx, assuming that the vertical displacement of the epineutral is weak compared to the zonal one (i.e., small slope approximation65), allowing a Taylor expansion of the denominator; and (iii) ∂_x σ and ∂_z σ are constant over Δx and δz, respectively. Equivalently, we have ∂_y~ σ ≈ ∂_y σ + s_y ∂_z σ, 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 important70,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 (Θ), by setting ∂_x~ γ = 0 and ∂_y~ γ = 0, such as:

s = - (α ∇_h Θ - β ∇_h S_A) / (∂_z γ), (9)

where ∇_h is the horizontal gradient, and α and β 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 ˆC as a function of Cartesian coordinates. This reads:

d_i ˆC_so ≈ ∂_z γ [∂_x (κ_x~x~ (∂_x ˆC + s_x ∂_z ˆC)) + ∂_y (κ_x~y~ (∂_x ˆC + s_x ∂_z ˆC)) + ∂_x (κ_x~y~ (∂_y ˆC + s_y ∂_z ˆC)) + ∂_y (κ_y~y~ (∂_y ˆC + s_y ∂_z ˆC)) + s_x ∂_z (κ_x~x~ (∂_x ˆC + s_x ∂_z ˆC)) + s_x ∂_z (κ_x~y~ (∂_y ˆC + s_y ∂_z ˆC)) + s_y ∂_z (κ_x~y~ (∂_x ˆC + s_x ∂_z ˆC)) + s_y ∂_z (κ_y~y~ (∂_y ˆC + s_y ∂_z ˆC))], (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 rates56,72 as:

Ψ_Θ = ∂_γ ∫Σ(γ0≤γ) α d_i Θ_so dΣ Δz, (11a)

Ψ_S_A = ∂_γ ∫Σ(γ0≤γ) β d_i S_A_so dΣ Δz, (11b)

where Σ(γ0≤γ) is the surface of a tracer class below the threshold value γ at 1000 m depth. Here Ψ_Θ and Ψ_S_A represent neutral density transformations through temperature or salinity changes over a layer of thickness Δz, respectively. These transformations are expressed in Sv, where 1 Sv ≡ 10^6 m^2 s^-1. In practice, the transformations are estimated using neutral density increment (we chose: Δγ = 0.03 kg m^-3). For consistency, the layer thickness can be estimated though the scaling: Δz ≈ Δγ (∂_z γ)^-1 (which derives from dγ = ∂_x γ dx + ∂_y γ dy + ∂_z γ 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 level40,41 and variations with Eddy Kinetic Energy36. Hence, we have chosen to restrict the layer thickness amplitude. We have used a uniform layer thickness based on the median value: Δz = 70 m (Supplementary Fig. 1b), regardless of the local stratification (∂_z γ). 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 (Ψ_Θ and Ψ_S_A) can be interpreted as the component of the total transformation induced by isopycnal turbulence. 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: Ψ_Θ + Ψ_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. ˆC and ∂_z ˆC (as well as ∂_z γ) at all locations [∀ (x, y)] along the depth, but primarily above the maximum Argo depth sampling of ~2000 m for robustness. In this study, C can be Conservative Temperature or Absolute Salinity (Fig. 1a, b);
2. K at all locations [∀ (x, y)] at the constant depth, which is determined by Argo deep displacements at their parking depth36 (Fig. 1d–f).

Hence, the main constraint occurs through 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 displacements36, which are available35 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 observations43 and available at https://doi.org/10.17882/52367. The diffusivity coefficients were obtained by reconstructing the 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 done33,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°×3° cells (i.e., the native cell of the diffusivity coefficients) on 1°×1° 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 study49. 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 κ_Uni&Iso_x~x~ = κ_Uni&Iso_y~y~ = (κ_x~x~ + κ_y~y~)/2 and κ_Uni&Iso_x~y~ = 0. These coefficients are applied to (5).
2. Spatially-varying and isotropic isopycnal diffusivity coefficients: the diffusivity coefficients are set to κ_Var&Iso_x~x~ = κ_Var&Iso_y~y~ = (κ_x~x~ + κ_y~y~)/2 and κ_Var&Iso_x~y~ = 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 = sign(log(Ψ_comp / Ψ_ref)) * |log(Ψ_comp / Ψ_ref)|, (12)

where Ψ_ref is the reference transformation (i.e., most complete diffusivity), Ψ_comp is the transformation used for comparison (i.e., less complete diffusivity), and the log term allow 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 coefficients34. This dataset is available at: https://figshare.com/articles/dataset/Groeskamp_et_al_2020_-_mixing_diffusivities/12554555. The Groeskamp diffusivity coefficients (κ_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: κ_Groeskamp_x~x~ = κ_Groeskamp_y~y~ = κ_Groeskamp and κ_Groeskamp_x~y~ = 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 coefficients39. This dataset is available at: https://github.com/croachutas/Isopycnal_Diffusivity. The Roach coefficients (κ_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: κ_Roach_x~x~ = κ_Roach_y~y~ = κ_Roach and κ_Roach_x~y~ = 0. This is applied to (5).

c. Since these two set of diffusion coefficients were computed as cross-flow diffusion coefficients, consistently, we have used our secondary diffusivity coefficient (κ_s) to have a fairer comparison with our diffusivity dataset39. 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: κ_Var&Isosec_x~x~ = κ_Var&Isosec_y~y~ = κ_s and κ_Var&Isosec_x~y~ = 0. These coefficients are applied to (5).

Finally, for comparison purposes, we have computed the diapycnal transformation as a vertical diffusion: d_i ˆC_ia = ∂_z (κ_z ∂_z ˆC), where κ_z is the vertical diffusivity coefficient. We have used spatially varying diffusivity coefficients59,60 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.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_3": "Summary\nThe authors apply a new estimate of epineutral mesoscale eddy diffusivities at 1000 dbar (based on statistics of Argo\nparking drift; published separately) to a T&S climatology (previously published) to estimate isoneutral eddy-driven water\nmass transformations in temperature and salinity coordinates. The key results are that (1) transformation rates are poorly\nestimated if one uses a globally-constant isoneutral eddy diffusivity (as most climate models presently do) and (2) most of\nthis transformation occurs in Southern Ocean water masses, where both diffusivities and mean isoneutral property gradients\nare large.\nMajor comments\nI have some major concerns about (1) the methodology, (2) the contextualization in terms of the broader literature, and (3)\nthe clarity of communication of the results. The results are interesting but I think not as ground-breaking as the abstract\nmakes it sound. I would recommend a decision of Reject and Resubmit.\nThe main result is the calculation of water mass transformation rates (in temperature and salinity coordinates) due to\nisoneutral eddy mixing. As the authors note several times, they are by observational constraints of the Argo parking depth\nconstrained to the 1000 dbar surface. What they do not sufficiently explain to the reader is that their choice of scaling these\npoint-wise transformation rates by a 100m layer is entirely arbitrary. In fact, it is not even really clear to me how the water\nmass transformation calculations are made. Water masses transformations are inherently three dimensional processes that\nthe authors are approximating by a two dimensional slice and then scaling up arbitrarily. The description in lines 410-414\n(plus equation 11) is confusing even to someone with a lot of experience thinking about water mass transformations, so I\nsuspect the general reader of Nature Comms will be lost. There are also some things I don’t understand about their results,\nsuch as: on L. 164-166, it doesn’t look like these actually line up that well to me? Can you just convert these to both be in Sv\nor Sv m-1 and add them together so we can see what the net transformation in density coordinates is? What does it mean\nabout your calculation if this is very different from zero?\nThe authors cite most (but not all) of the key references in the literature, but they never contextualize any of their results in\nterms of these references. For instance, they mention early on that isopycnal transformations have previously mostly been\ninferred indirectly using inverse models (I am thinking mostly of the Groeskamp papers), but they never circle back to discuss\nhow their results compare to these estimates! Are they very similar? Very different? Should we trust yours more, or theirs?\nDo the differences tell us anything? Also, they claim to have “unveiled” these “hidden” isopycnal transformations, but they\nare hardly the first to make estimates of these transformations–even if they may have a case for their estimate being the most\nsophisticated since it uses independent datasets with a spatially-variable eddy diffusivity. I recommend referencing papers\nthat discuss the impact of isoneutral mixing in climate models (two recent examples:\nhttps://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2019GL085208 and\nhttps://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2021MS002914).\nThe roles of eddy-bolus velocities (or the skew-symmetric part of the diffusivity tensor) and the diapycnal velocity in the\nocean & climate system are unjustifiably downplayed to make their isoneutral mixing sound more important than it is. The\nauthors in some places forget the qualifier that they are only discussing the isopycnal mixing role of eddies and make\nstatements like “since eddies induce a flow along isopycnals, they cannot modify density… hence the circulation”, which are\nnot true if one always considers the leading-order quasi-advective role that eddies play in the general circulation. Also, the\ncomparison with diapycnal mixing transformations (19x smaller!) feels disingenuous; the reason people usually talk about\ndiapycnal mixing is not because how it transforms the temperature or salinity of water masses, but how it does so\ndiapycnally in a way that no other process in the ocean interior can.\nSome of the key points are poorly articulated and difficult to decipher. In most cases I could eventually decipher what you are\ntrying to say, but not always (see minor comments).\nFinally, the braggadocious tone of the writing is off-putting and I encourage the authors to tone it down. The entire\nmanuscript is written as though the authors are the first to realize that eddies are rapidly mixing temperature and salinity\nalong isopycnals, or at least the first to make an estimate of it. Some other examples: 1) Water mass transformation is a\nframework that is more than 50 years old, so calling it “state-of-the-art” is a tough sell. 2) The other “state-of-the-art”\nframework is the authors’ novel AI-based method for estimating turbulent anisotropic diffusivities. It is indeed very interesting\nand–to my knowledge–a useful novel approach, but the authors don’t actually go any further in describing it in this\nmanuscript than stating its name (the “Analog method”) and citing their prior work! 3) Given how impactful this data set is in\nthis manuscript, it deserves at least a paragraph in the method section. (There is plenty of fluff in the introduction and\nconclusions that you could cut out to make space).\nMinor comments:\nEq. (4) – I understand where the authors are coming from, but it doesn’t really make sense to decompose the transformations\nin this way because isopycnal + diapycnal already accounts for all directions. Either choose a two-way decomposition\n(isopycnal vs. diapycnal) or explain that you are breaking these down by process: mixed-layer turbulence (strong in all\ndirections), interior mesoscale turbulence (isopycnal), and small-scale turbulence (in all three directions, but the isopycnal\ncomponent is swamped out by the mesoscale isopycnal so can be approximated by diapycnal).\nFigure 3. Colormaps are not very interpretable and I think very bad for people with colorblindness.\nFigure S1. Why not convert both to Sv so they can be directly compared? “For the isotropic diffusivity” I think should read\n“For the anisotropic diffusivity”. Caption should clarify these are “Water mass transformations in density space” or something.\nL. 374-375: I don’t understand this.\nWhy epineutral rather than isoneutral? Is epineutral more specifically a geometric descriptor?", "reviewer_4": "Response to the reviews of the manuscript\nTurbulent Isopycnal Mixing dominates Thermohaline Transformations of Intermediate\nOcean Waters\nby F. Sévellec, N. Kolodziejczyk, and E. Portela\nWe thank the two reviewers and the editor for their constructive comments, which have\nhelped us improving and clarifying the manuscript. Our detailed response is below.\nFont legend:\n1. Reviewer’s comments are in bold,\n2. Our responses are in normal font, and\n3. Text added in the manuscript is in italic.\nReviewer #3 (Remarks to the Author):\n2nd Review of Turbulent Isopycnal Mixing dominates Thermohaline Transformations of\nIntermediate Ocean Waters by Florian Sevellec, Nicolas Kolodziejczyk, and Esther Portela\nI appreciate the improvements the authors have made. Unfortunately, I still don’t feel there\nis enough clarity of the methods, embedding of the subject in literature nor clarity of the\npresentation that the manuscript is ready for publication. I recommend major revisions.\nThank you for your supportive comment. We hope that the modifications make this new\nversion now ready for publication.\nNovelty\nI’m now more enthusiastic about the work than I was before. I think the idea of comparing\nthe isopycnal transformation of salt and temperature (which have received relatively little\nattention) to their diapycnal counterparts, a worthwhile exercise. However, the reason it\nhas received less attention may also be because it is less interesting. The authors are trying\nto make the case that this is in fact very interesting, but I’m not yet convinced based on\nthese results. I’m still open to be proven otherwise. Below I will continue to make myself\nclear.\nAgain, thank you for this overall supportive comment. For initial clarification (below you will\nfind our point-by-point specific response), we would like to say that we do not want to\ndemonstrate that thermal- and/or haline- transformations are more important than density-\ntransformations. We are simply concentrating on these transformations because they are\nalso a key aspect of the ocean circulation. This is important in the context of heat and\nfreshwater transport (though it has little effect on overturning). Indeed they are key aspects\nof the ocean circulation that are currently under stressed within climate changes. We agree\nwith the reviewer that they have somehow received less attention.\n1 /20\nThe main result as presented\nBecause the ocean does move along and across neutral surfaces, so does T&S. To get local\nchanges in T or S we thus have to account for local air-sea-ice interaction and mixing\nprocesses (and maybe geothermal heating). These mixing processes are either mesoscale\nhorizontal (in mixed layer), mesoscale isoneutral (in deep ocean) or vertical. The results can\nbe used to calculate WMT in T, S or density coordinates. In the S&T framework this leads to\nlarger transformation (Sv) from one temperature to another (or salinity) due to isoneutral\nmixing than due to vertical mixing. That is the main finding of this paper. This result,\nhowever, is not put in context of other studies and is wrongly compared to density\noverturning.\nWe agree with the general summary of our work. Regarding the last part, we have improved\nthese aspects in the manuscript. Please see the detailed response below.\nNote that comparing density-transformation with thermal- and haline- transformation is\npossible. Higher value of the latter over the former does demonstrate that the ocean is also\nsignificantly “spicy”. Hence, solely concentrating on density transformation “hides” important\nheat and freshwater transports. This is what we try to demonstrate quantitatively.\nPutting this results in perspective\nFirst, other literature on this subject is provided by the vast amount of old and new studies\nthat have looked at transformation in T&S coordinates. This is essentially what has been\ndone here as well as this is exactly diahaline and diathermal transformation by these\ndifferent mixing processes. Many diagrams have been shown and never have I seen a factor\n20 difference. Maybe I missed something, but that should then be made clear by a proper\nembedding into this literature. Below is a list of references for this subject.\nThank you for this list of references. They are now included in the manuscript. We did not see\nany inconsistency with previous work. We have revised our estimate by acknowledging the\nspatial variation of diapycnal/vertical diffusion coefficient. We now have a factor 12. The fact\nthat our computation is solely located at a single depth can also explain such a large value\n(i.e., other depth may have weaker value).Hence, we still find that at 1000-m depth the\ntemperature and salinity transformations are widely dominated by isopycnal (but diathermal\nand diahaline) processes over diapycnal processes.\nWe have now acknowledged the literature in the introduction and used it to develop our\nscientific question. Also we used it to provide further context to our conclusion.\nSec.[The Role…] p.3 l.54-61: “There has been a large range of studies on Water Mass\nTransformations23-28. Current estimations disagree on the dominance of isopycnal diffusion29\nor of diapycnal diffusion67,30 for the temperature and salinity transformations in the ocean\ninterior. However these studies rely on either numerical models or inverse method, both\nintrinsically depending on a chosen closure scheme (and sometimes imposed parameters).\nHence it is key to use direct observations to assess whether isopycnal or diapycnal diffusion is\ndominant in the ocean interior. Indeed, the use of direct observations is the only way to\nvalidate or invalidate hypotheses derived from these numerical or inverse modelling\napproaches. This is the goal of the present study.”\n2 /20\nSec.[Conclusions] p.14 l.272-277: “In the literature there exists disagreement regarding the\nkey role of isopycnal diffusion in comparison to diapycnal diffusion for temperature and\nsalinity transformations6,297,30. Our direct estimates suggest dominance of isopycnal diffusion\nat 1,000-m depth. We hypothesize that this is because we fully acknowledged the anisotropic\nnature of the isopycnal diffusion within our framework, unlike in inverse or numerical\nmodelling studies. This dominance is particularly important since isopycnal transformations\nare key for heat uptake in the context of Global Warmin29.”\nThen there is the comparison to density-driven circulation. It is often repeated that the\ndensity-compensated T&S transformation is very important. I’m not 100% convinced why,\nbut here is one reason I can think of:\n- There is only so much mixing and air-sea-ice interaction that can change T&S\n(transformation potential). We see that a lot of this work goes into density-\ncompensated changes. As a result, this reduces the impact of these fluxes in the ocean\ncirculation, through density transformations.\nThis is a key point. The ocean structure is significantly density compensated. This come from\nthe fact that, to first order, warming regions are regions of net evaporation and cooling\nregions are regions of net precipitation. Hence there is a significant part of the ocean structure\nthat is “spicy” (warm and salty or cold and fresh). This is well established.\nBeyond this general comment, at 1000-m depth we have waters that are showing strong\nspiciness contrasts (e.g., there AAIW, MW, and RSPPGIW, as could be seen in Fig. 2). This\nmakes the studied depth particularly relevant to describe the “hidden” circulation.\nFollowing your comment, we have added a paragraph on this point in the text:\nSec.[Water Masses…] p.8 l.155-165: “We find transformations of the order of ±0.22 Sv m-1\nand ±0.13 Sv m-1 (i.e., latitudinal standard deviation) for temperature and salinity,\nrespectively (Fig. 5a and b, respectively). This is important as, for instance, this transformation\nwithin only a 100-m layer is non-negligible in comparison with the +100 Sv of total external\nbuoyancy forcing53. We rationalize the importance of the temperature and salinity\ntransformations over the density ones because of the strong thermohaline compensating\neffect on density (i.e., the so called spiciness). The water masses located at our studied depth\nexhibit particularly intense spiciness contrasts, although it is not necessarily the case at other\ndepths. More broadly, at the ocean surface, warming regions are often regions of net\nevaporation, whereas cooling regions are often regions of net precipitation. This means that\nthere is a non-negligible part of the ocean circulation that transports heat and freshwater\nfrom cold and fresh regions to warm and salty regions without transporting density.”\nYou can see this, because despite vertical mixing having caused little diahaline/diathermal\ntransformation, these small numbers still have a huge impact on the density overturning\n(for example when I looked at reference 39, black line, figure 2). This is of the same order\nas the number of Sverdrups of the diahline/diathermal transformations that are density\ncompensated. On top of that, the nonlinear part (Cabbeling/thermobaricity) causes up to\n3 /20\n20Sv of diapycnal transformations, especially in the Southern Ocean (Fig 3 of reference 39).\nThat is about half of the “hidden” isopycnal diahline/diathermal transformations. So these\nnonlinear parts that are often referred to as “small” in the manuscript actually have\nsubstantial numbers associated with them. Especially for the density overturning.\nA. The comparison of the results throughout the manuscript give the wrong impression\nof the importance of the results. They have comparable sizes (within a factor 3 up\nto the same magnitude) and this should be rectified.\nB. It is not put in perspective that the diapycnal component of diahaline/diathermal\ntransformation may be small, but it has a huge effect on density overturning. This\nneeds to be acknowledged and put in context for us all to understand why we care\nabout these results. So within diathermal/diahaline transformation they are small,\nbut their impact is huge for circulation. The impact of the isopycnal transformation\nis in fact smaller for circulation (half or so).\nC. This discrepancy can be explained by the “transformation potential”. Because you\nhave basically showed that a lot of the T&S transformation is density compensated\n(a sink for affecting overturning) and the part that affects overturning is small, but\nimportant for world ocean.\nThank you for this clear and detailed comment.\nWe think the way we have diagnosed the transformation by assuming a consistent 100 m\nlayer as led to some confusion. Our 20 Sv of T and S transformations was only for a layer of\n100 m. In this sense it is significantly larger than the transformation by cabbeling and\nthermobaricity that act on a larger layer.\nHowever we are not arguing that density transformation is less important than temperature\nand salinity transformations. The comparison with density transformations was only there to\ngive a sense of scale to the reader with the most usual framework. We have rephrased the\ntext to be clearer on the goal of the comparison and avoid misunderstanding.\nWe have modified this throughout the manuscript consistently with your comments and\nclarified that cabbeling and thermobaricity are important for density transformations, and\nthus overturning.\nSec. p.5 l.100-102: “This transformation-induced process has almost no signature on the\ndensity field (except for the weak cabbeling and thermobaricity components, although they\nremain crucial for the density transformations53)”\nImproving Methods\n• Fig. 4, Section on WMT and appendix on diffusivity data: I appreciate the use of\ndifferent isopycnal diffusivities. But I think you should not interchange isopycnal\ngradient with horizontal gradients. We know this is wrong and we can easily do it\nright. Instead, it is more insightful to change the diffusivities you use. Because these\nare values, we don’t know for sure. There are quite a few product out there that\nprovide diffusivities at 1000m depth and I think you could instead use these. So:\ntesting different diffusivities in an isopycnal framework. I advise using the work of\n4 /20\nyour references 18, 19 (Groeskamp et al and Cole et al) and in addition that of Roach\net al (see below).\no I now just read in your response to previous review comments that you are\nreluctant to do this, because you don’t like the way these diffusivities are\nobtained. I disagree that this is a reason not to use them. They are independent\nestimates compared to Sevellec et al 2022. None of them (neither Sevellec 2022)\nare perfect. They all have their own flaws and we are not yet at the stage to\ndetermine which one is “best”. Until we are, it is worth comparing them.\nThank you for this comment, we agree. Following your suggestion, we have now computed\nisopycnal diffusivity with two other isopycnal diffusion coefficient estimations. We have\nadded a new figure on this point (i.e., Fig. 6), reproduced below (Fig. R1). Despite some\ndifferences we find a consistent result within their intrinsic restrictions (i.e., cross-diffusion).\nThanks to your comment this further demonstrates the robustness of our results and the key\nrole of anisotropic diffusion.\nFigure R1: Other estimations of isopycnal diffusivity and induced Temperature and Salinity water mass\ntransformations at 1,000-m depth. (a and b) Isopycnal diffusivity coefficients at 1,000-m estimated with\nFigualrteernSat3iv:e Ometthhoedrs3e2,58 (idmenaottieodn ass Gorfoeisskoapmypc annad lRodaicffh,u rseisvpeitcytivealyn)d. Wiantderu mcaesds trTanesmfoprmeartaiotnusr feor a(cn) d\nSaliCnointsyervwaativtee rTemmpaesrastutrre aannds ffoorr m(d)a Atbisoonlustea Stalin1i,t0y.0 B0lu-me lindese sphtohw. tra(nasfaonrmdatbio)nIs sfoopr ythcen iasoltdroipffiuc sainvdit y\ncoeffispcaietinatllsy-avatr1yi,n0g0 0d-iffmuseivsittiym. Saotleidd awndit hdaashlteedr nblautei vlienems eshthowod rseYs,uZlts( dfoern tohtee davaesraGger obeestwkaemenp prainmdarRy oanadch ,\nrespseecctoivnedlayr)y. dWiffuastievirtiems aasnsdt sreacnosnfdoarrmy adtififounsisviftoiers,( rce)spCeoctnivseelryv. aOtriavnegTe eamndp reerda tliunrees sahnodw ftorran(sdfo)rAmabtsioolnus te\nfor the Groeskamp and Roach diffusivities, respectively. The curves have been smoothed using a moving\nSalinity. Blue lines show transformations for the isotropic and spatially-varying diffusivity. Orange\naverage of 0.06 Kg m-3 re-interpolated on regular 0.012 Kg m-3 grid with a spline interpolation. The horizontal\nand red lines show transformations for the Groeskamp and Roach diffusivities, respectively. The\ncolor patches show neutral density classes of selected Water Masses45.\ncurves have been smoothed using a moving average of 0.06 Kg m−3 re-interpolated on regular\n0.012 Kg m−3 grid with a spline interpolation. The horizontal color patches show neutral density\nclasses of selected Water Masses32.\nWe have added a figure (Fig. 6) and a text regarding this point:\nSec.[Robustenss…] p.11-12 l.210-242: “To test the robustness of our results, we have\nrecomputed the temperature and salinity transformations using two alternative diffusivity\ncoefficient datasets (see Methods – Alternative Diffusivity Coefficients). The first dataset, that\nwe named “Groeskamp” diffusivity coefficients, used a different set of observations to infer\n5 /20\n3\ndiffusivity coefficients than our estimates32. Here the diffusivity coefficients are estimated\nthrough salinity observations assuming the validity of the mixing length theory. The second\ndataset, that we named “Roach” diffusivity coefficients, used the same deep observations\nthan our estimates, but a different methodology: a pairwise dispersion of deep Argo float\ndisplacements59.\nWhen comparing our diffusivity estimation (Fig. 1d and e) from Groeskamp and Roach ones\n(Fig. 6a and b, respectively), we find the same overall pattern, with intensification in western\nboundary currents and within the Antarctic Circumpolar Current. Because Groeskamp and\nRoach diffusivity coefficients correspond solely to the cross flow diffusion, they are more\ncomparable to our secondary diffusivity coefficients (Fig. 1e), which is by construction weaker\nthan the primary one. This secondary diffusivity coefficients are between 2,000 m2 s-1 and\n6,000 m2 s-1 in active regions. This is comparable to Groeskamp and Roach diffusivity\ncoefficients. In comparison with the two other datasets, Roach diffusivity coefficients are less\nspatially-smooth.\nBecause Groeskamp and Roach datasets provide only a single local diffusivity coefficient, we\nrestrict the comparison to the spatially-variable and isotropic framework. Computing the\ntemperature and salinity transformations with the different diffusivity coefficient datasets\nshow the good consistency of the results (Fig. 6c and d). When using the spatially-variable and\nisotropic diffusivity coefficients, we find temperature and salinity transformations that are\nstronger than those obtained with Groeskamp and Roach diffusivity coefficients. However the\ntransformations are qualitatively well represented. As expected, the agreement is better when\napplying the spatially-variable and isotropic diffusivity coefficients solely based on the\nsecondary diffusivity coefficients. Here the agreement is qualitative and quantitatively good.\nThis stresses the robustness of our described temperature and salinity transformations with\nrespect to the methodology used to infer the diffusivity coefficients.\nThe good agreement of the three datasets regarding the salinity transformations (Fig. 6d) is\nparticularly interesting. Indeed Groeskamp diffusivity coefficients are based on an inversion\nmethod using salinity observations (i.e., somehow the diffusivity coefficients carry the\nknowledge of the salinity transformations). Hence the agreement of our analysis regarding\nsalinity transformations regardless of the diffusivity coefficient datasets confirms the\nrobustness of the overall methodology.”\n• Using a constant diapycnal diffusivity is as misleading as using a constant mesoscale\ndiffusivity. While the latter is something you test, the former is not. I advise using\nthe one by Lavergne et al. from reference 53. As with isopycnal mixing, also the\nvertical mixing is mostly impacted by peaks in diffusivities, and this can alter the\nrelative importance of their diahaline transports. Also, please mention the used\ndiffusivity in the text and not buried away in the last line of the appendix. I think\nyou have only used Lavergne as a “guide” to come up with the constant value. If you\nare not going to use this variable diffusivity, then it is important to discuss this more\nthoroughly in the paper.\n6 /20\nThank you for this comment. We agree. It is better and easier to use the spatially-variable\nvalues of Lavergne. We have updated the manuscript. This does not qualitatively change\nthe results, but only slightly change the quantitative values. The ratio of\nisopycnal/diapycnal transformations for T and S at 1000-m depth is now 13 and 8 (it was\nand 18 and 19 with spatially-constant vertical diffusivity coefficient), respectively. Please\nsee the new results in Fig. R2 and R3.\nThe manuscript has been updated accordingly, with the modification of Fig. 5 and addition\nof one supplementary figure (Fig. S1).\nFigure R2: Vertical diffusivity at 1,000-m depth. Vertical diffusivity coefficients55,56 at 1,000-m used to\nassess the temperature and salinity diapycnal transformations.\n(a) CONS TEMP WATER MASS TRANSFORMATION\nCDW MW\n0.4 AABW NADW RSPGIW\nFigure S1: Diapycnal diAffAIuWsivity at 1,000-m depth. Diapycnal diffusivity coefficientsX at 0.2\n1,000-m. 0\n-0.2\n-0.4\n-0.6\n-0.8\n-1\n-1 0 1 2 3 4 5 6 7 8 9\nCONSERVATIVE TEMPERATURE (°C)\nFigure R3: Temperature and Salinity water mass transformations at 1,000-m depth. Water mass\ntransformations for (a and c) Conservative Temperature and for (b and d) Absolute Salinity. Positive or\nnegative values show warming or salinification and cooling or freshening, respectively. (a and b) Red and\nblue bars show positive and negative transformations for the isotropic diffusivity, respectively; dark and\nlight grey bars show transformations for the isotropic spatially-varying and globally-uniform diffusivity\n1\n7 /20\n)1-m\nvS(\nNOITAMROFSNART\n(b) ABS SAL WATER MASS TRANSFORMATION\n0.5 AAIW MW\n0.4 AABW RSPGIW\nCDW NADW 0.3\n0.2 0.1\n0\n-0.1\n-0.2\n-0.3\nVar & Aniso - pos Var & Iso Horizontal -0.4\nVar & Aniso - neg Uni & Iso Diapycnal -0.5\n-0.6\n34.3 34.4 34.5 34.6 34.7 34.8 34.9 35 35.1 35.2 35.3\nABSOLUTE SALINITY (g kg-1)\n)1-m\nvS(\nNOITAMROFSNART\n0.4\n0.2\n0\n-0.2\n-0.4\n-0.6\n-0.8\n-1 0 1 2 3 4 5 6 7 8 9\nCONSERVATIVE TEMPERATURE (°C)\n)1-m\nvS(\nNOITAMROFSNART\n(c) REGIONAL DECOMPOSITION\nCDW MW 0.4\nAABW NADW RSPGIW\nAAIW 0.3\n0.2\n0.1\n0\n-0.1\nGLOBAL SO-A SO-P -0.2\nNA SO-I DIAPYCNAL\n-0.3\n34.3 34.4 34.5 34.6 34.7 34.8 34.9 35 35.1 35.2 35.3\nABSOLUTE SALINITY (g kg-1)\n)1-m\nvS(\nNOITAMROFSNART\n(d) REGIONAL DECOMPOSITION\nAAIW MW\nAABW RSPGIW\nCDW NADW\nFigure 6: Temperature and Salinity water mass transformations at 1,000-m depth.\nWater mass transformations for (a and c) Conservative Temperature and for (b and d) Absolute\nSalinity. Positive or negative values show warming or salinification and cooling or freshening, re-\nspectively. (aandb)Redandbluebarsshowpositiveandnegativetransformationsfortheisotropic\ndiffusivity, respectively; dark and light grey bars show transformations for the isotropic spatially-\nvarying and globally-uniform diffusivity (set as the local-mean and global-mean of the anisotropic\ndiffusivity), respectively. The purple lines are the transformations estimated in a horizontal dif-\nfusivity framework. The cyan lines are the transformations induced by diapycnal diffusivity. (c\nand d) The transformations are split in geographical sectors: global (black lines) North Atlantic\n(NA, orange line), Atlantic sector of Southern Ocean (SO-A, blue line), Indian sector of Southern\nOcean (SO-I, yellow line), and Pacific sector of Southern Ocean (SO-P, green line), as represented\nin Fig. 1f. Thecurves have been smoothed using a moving average of 0.4◦C and 0.044 g Kg−1 and\nre-interpolated on regular 0.01◦C and 0.0011 g Kg−1 grid with a spline interpolation, respectively.\nThe horizontal color patches show properties of selected Water Masses32.\n8\n(set as the local-mean and global-mean of the anisotropic diffusivity), respectively. The purple lines are\nthe transformations estimated in a horizontal diffusivity framework. The cyan lines are the\ntransformations induced by diapycnal diffusivity. (c and d) The transformations are split in geographical\nsectors: global (black lines) North Atlantic (NA, orange line), Atlantic sector of Southern Ocean (SO-A, blue\nline), Indian sector of Southern Ocean (SO-I, yellow line), and Pacific sector of Southern Ocean (SO-P, green\nline), as represented in Fig. 1f. The curves have been smoothed using a moving average of 0.4oC and 0.044\ng Kg-1 and re-interpolated on regular 0.01oC and 0.0011 g Kg-1 grid with a spline interpolation,\nrespectively. The horizontal color patches show properties of selected Water Masses45.\n• Calculating isopycnal gtradients is done in the “Redi” way. This method is known to\nhave problems in regions of weak stratification. This has a strong impact on\nisopycnal gradients and related WMT rates. Especially in the Southern Ocean en for\ncabbeling. I advise a better way of calculating these gradients using either\nisopycnally averaged values or a searching algorithm.\nThank you for this comment. Over the large-scale of the computation and depth considered,\nwe do not have issue with weak stratification (slope remains below 1%). Nonetheless, we\nhave added the references.\nMethods.[Computation…] p.27 l.506-508: “Note that there exists more accurate computation\nmethods of the neutral slopes and gradients, which are key when the slopes become\nimportant67,68.”\n• L144 – Layer thickness is important. Can you check sensitivity to this choice? The same\nfor the isopycnal binning you use in combination with this layer thickness.\nThe layer thickness was chosen to simplify the interpretation (volume transformation in Sv)\nand make the result easier to grasp. However we do not have any robust reason to set it to a\nspecific value. Given its questionable value and confusion it led to, we have now chosen to\nremove it. We feel that it is better to stay as close as possible to the actual computation and\nits limit (i.e., single depth leading to “surface” transformation in Sv m-1).\nOther comments\nIn the abstract it says: “This has major consequences for our understanding of climate\ndynamics and of the role of the ocean in response to anthropogenic climate changes.”\nThis is easy to say, but I’m not yet convinced. It requires a clearer explanation why.\nWe feel that observing (or more specifically determining through direct observations) the\nisopycnal temperature and salinity transformations has direct consequences on heat and\nfreshwater transport, hence on the role of the ocean in the climate.\nL91 – in the SO and some other places, cabbeling is not weak at all.\nWe have modified the sentence.\nL93 – add: through eddy transport and mixing.\nThese sentence have been modified.\nL95 – The former, you mean eddy transport? Please specify process.\nThis has been clarified\nL102 – of of\nThank you.\nL111 – Do you have an explanation for the dipole patterns in the Agulhas?\nIt corresponds to two water class: AAIW and residual of RSPIGW. The latter is quite present\nclose to the Agulhas at 1000-m depth (see the local maximum of T and S in Fig1a and 1b). This\ncan be inferred in the transformation within the TS diagram.\nSec.[Rate of Changes…] p.6-7 l.122-127: “There is also significant rate of changes along the\nAgulhas return current and its merging with the Antarctic Circumpolar Current, which also\ncorresponds to a local maximum of diffusivity coefficients. Indeed, between South Africa and\nthe Kerguelen plateau, there is a slight warming and salinification of the relatively cold and\nfresh AAIW (Fig. 3a and 4a) and a significant cooling and freshening of the relatively warm\nand salty residual of RSGPIW (trapped in the Agulhas system53, Fig. 1b).”\nL113 – You refer to S&T figures, but only introduce these in the sentence below…This is\ninconsistent.\nThank you. This is corrected.\nL116-L117 – You can remove this sentence.\nDone\nSection WMT – see main comments.\nThank you. We have modified the text according to your (main) comments.\nL151 – why is warming and freshening expressed as Sv? His is confusing. Do you mean, this\namount of cooling, requires a transformation of this amount of Sv with cooler water”\nThe unit is the outcome of the thermal transformations in Temperature space. This means\nthat a warming of X Sv of YoC water class corresponds to a X Sv rate of volume of water\n“created” (i.e., transformation) at this particular YoC temperature. (Since nothing is “created”\nthis is compensated by the same amount of water lost at another temperature class.)\nNote however that now our unit is in Sv m-1. Now, we are not dealing with volume change\n100 m around 1000-mn depth, but area change at 1000-m depth\nL166 / Fig S1 – what are the units? What exactly is done here?\nThis is strictly identical to the computation of temperature/salinity transformations, but now\nit is expressed in transformation per density class (vs T and S class). In this framework, we are\ncloser to “classical” density transformation, but have spitted the role of transformation in\nterm of temperature and salinity changes.\nMathematically this is applying density in the integral and derivative of (11), but keeping\ntemperature and salinity in the expression of the transformation term.\n9 /20\nL171 – Not sure if Veronis the right reference here. That is about getting the direction of the\nisopycnal surface correct, not about the variability of the coefficient (if I remember well). If\nyou mean to say that you have used horizontal gradients than that’s a different story. This\nwould be the correct reference, but it is not written clearly.\nWe agree with the reviewer. Here the key word was “horizontal”. We have rephrased to avoid\nmisunderstanding.\nFig 1 and Fig 2: I feel like these can be combined as one does not provide much more\ninformation than the other. If you do keep them, you should put the density labels into fig\n2.\nWe kept the two figures, because we have not reach the manuscript length limit. As figure 2\nis intended to be an schematic, density values for the isopycnals (gray lines) are given in the\ncaption to keep the figure simple .\nFig3 – The insets really are too small to read. You need separate (sub)figure for this.\nWhy do you “order them”. This way you focus on the highest value only, which is\nmisleading. Is the average over a WM not more important?\nThank you for this comment. Following your advice, we have split the figure in two. One for\ntemperature and one for salinity. This allows for larger subfigures and insets. We chose to\nkeep the geographical and TS representation in the same figure, since it is important to see\nboth at the same time to grasp the information of location and Water Masses together.\nThe order is to avoid feeling space with values with no transformations (there is plenty of\nthem). Also dot have roughly the same “values” in term of transformation (i.e., they represent\nthe same spatial surface). Finally this allows to easily comparing (top) geographical-space and\n(bottom) TS-space, all with the same colorbar. However overall this gives the good sense of\nthe transformation and high transformations are easily related to a density class.\nCOMMENT PREVIOUS REVIEW (we kept only your last comment/response):\nPage 3 first paragraph. This section acknowledges mixing obtained from inversions. This is\nregardless of your opinion about which method is best. So it is normal to first acknowledge\nthese effort. You can then argue nicely about their shortcomings and why you choose to\nuse a different method.\nSorry for misinterpreting your initial comment. We have now included your suggested\nreferences. They are indeed relevant in the context of the paragraph.\nSec.[The role of…] p.3 l.47-50: “Transformations of seawater properties can occur either by\nocean-atmospheric interactions or by turbulent diffusion in the ocean interior. Whereas the\nformer is somehow more easily accessible8,14, the latter has always been more difficult to\nassess and is often deduced using inverse methodology15–18”\n10 /20\nP18 – eq5: is this correct? In fact, the neutral tracer gradient is 3D, but sometimes\napproximated as 2D. See your reference 50 for a good explenation and the difference\nbetween them. Notation wise this means you need to use “small n” and explain that it is\nthe projected nonortogonal approach. If you do, it should be correct within that\napproximation.\nThank you for this comment. We have modified the equation and text accordingly to it.\nIndeed, as you suggested, in our framework, the gradient is already along (or “project” onto)\nthe density surface.\nMethods.[Isopycnal Turbulent…] p.25 l.482: “Here the gradient is computed using the\nprojected nonorthogonal approach65,66.”\nP18, last section into P19. You do have the Veronis effect. The fact that you have spurious\ntracer flux in a non-neutral direction can be translated as false mixing. You get this when\nyour slopes used to calculate neutral gradients are not exactly neutral, so they cause this\nspurious transport and lead to mixing that is not actually there.\nWe agree with the reviewer. We have added a comment on the fact that our assumption (but\nin reality any assumption) on the computation of the neutral density slope would lead to a\nVeronis effect.\nSec.[Effect…] p.5 l.99-100: “This implies that the transformations act within isopycnal layers\nas epineutral flux (Fig. 2), limiting the Veronis Effect51 in our estimation.”\n11 /20\nReviewer #4 (Remarks to the Author):\nSummary\nThe authors apply a new estimate of epineutral mesoscale eddy diffusivities at 1000 dbar\n(based on statistics of Argo parking drift; published separately) to a T&S climatology\n(previously published) to estimate isoneutral eddy-driven water mass transformations in\ntemperature and salinity coordinates. The key results are that (1) transformation rates are\npoorly estimated if one uses a globally-constant isoneutral eddy diffusivity (as most climate\nmodels presently do) and (2) most of this transformation occurs in Southern Ocean water\nmasses, where both diffusivities and mean isoneutral property gradients are large.\nThank you for this concise and accurate summary of our work.\nMajor comments\nI have some major concerns about (1) the methodology, (2) the contextualization in terms\nof the broader literature, and (3) the clarity of communication of the results. The results are\ninteresting but I think not as ground-breaking as the abstract makes it sound. I would\nrecommend a decision of Reject and Resubmit.\nThank you for your careful reading and assessment of our manuscript. Please see below our\npoint-by-point response to your comments, as well as the specific modifications they led to.\nMethodology\nThe main result is the calculation of water mass transformation rates (in temperature and\nsalinity coordinates) due to isoneutral eddy mixing. As the authors note several times, they\nare by observational constraints of the Argo parking depth constrained to the 1000 dbar\nsurface. What they do not sufficiently explain to the reader is that their choice of scaling\nthese point-wise transformation rates by a 100m layer is entirely arbitrary. In fact, it is not\neven really clear to me how the water mass transformation calculations are made. Water\nmasses transformations are inherently three dimensional processes that the authors are\napproximating by a two dimensional slice and then scaling up arbitrarily. The description in\nlines 410-414 (plus equation 11) is confusing even to someone with a lot of experience\nthinking about water mass transformations, so I suspect the general reader of Nature\nComms will be lost.\nThank you for this comment. This has significantly helped us clarify the manuscript and\navoid/remove unneeded arbitrary choices.\nFirst of all we have removed the arbitrary scaling by a layer-thickness of 100 m. So that now\nTS-transformation are in Sv m-1. All the figures have been updated (i.e., Figs. 5, 6 and S2.) We\nstill discuss this 100-m scaling in the text, to give a point of comparison with previous results,\nbut it is not key anymore.\nRegarding the computation of the transformation. The text in l.410-414 as well as the\nequation (11) are standard. We have provided references for them. Mathematically it\n12 /20\ncorresponds to the probability density function of cooling and warming (or freshening and\ngaining salt) per temperature classes (salinity classes). This shows how much volume of water\nat a given temperature is cooled/warmed. We have added this description in the text:\nSec.[Water Mass…] p.28 l.519-522: “Temperature or salinity transformations mathematically\nexpressed a probability density function of cooling/warming or loosing/gaining salt per\ntemperature or salinity classes, respectively. They physically correspond to the surface of\nwater at a given temperature that is cooling/warming or at a given salinity that is\nlosing/gaining salt, respectively.”\nThere are also some things I don’t understand about their results, such as: on L. 164-166, it\ndoesn’t look like these actually line up that well to me? Can you just convert these to both\nbe in Sv or Sv m-1 and add them together so we can see what the net transformation in\ndensity coordinates is? What does it mean about your calculation if this is very different\nfrom zero?\nAs explained above, they cannot be directly compared, even when plotted per density classes,\nsince they are temperature-transformations or salinity-transformations. Even if plotted with\nthe same unit (i.e., Sv), they correspond to two radically different physical processes: changes\nin temperature and changes in salinity.\nTwo alternative methods exist that have been acknowledged/cited in the manuscript:\n(1) Conversion to density changes: If one convert them in term of density-changes in density\nclasses, then they are comparable. They will be by construction close to zero (small change of\ndensity only due to thermobaricity and cabbeling). These has been done numerous times\n(these studies are cited in the manuscript), and we will not be able to disentangle the impact\nof turbulent isopycnal mixing on temperature and salinity transformation, which is the goal\nof the present study.\n(2) Computation of volume changes: It is possible to compute volume change in T-S space.\nThis a different approach to temperature and salinity (or density) transformations. Then the\ntransformation is often inferred through inverse method. This is not the chosen strategy here,\nsince we want to explicitly test (i.e., direct method) the role of isopycnal turbulence term (i.e.,\nisopycnal diffusion) on temperature and salinity.\nWe hope this clarifies our approach.\nContextualization\nThe authors cite most (but not all) of the key references in the literature, but they never\ncontextualize any of their results in terms of these references. For instance, they mention\nearly on that isopycnal transformations have previously mostly been inferred indirectly\nusing inverse models (I am thinking mostly of the Groeskamp papers), but they never circle\nback to discuss how their results compare to these estimates! Are they very similar? Very\ndifferent? Should we trust yours more, or theirs? Do the differences tell us anything?\n13 /20\nWe agree that the lack of comparison with previous studies made the manuscript weak and\nthe significance of our results hard to estimate.\nTo overcome this issue we have used two strategies:\n(1) We have both contextualized our scientific question and results in the introduction and\nconclusion, respectively;\n(2) We have re-computed the analysis using alternative diffusivity coefficient datasets given\nin the literature.\nPlease see below the details of (1) and (2) rationale and induced manuscript changes.\n(1) We have now acknowledged the literature in the introduction to develop our scientific\nquestion. Also we used it in the conclusion to give a wider context to our conclusion.\nSec.[The Role…] p.3. l.54-61: “There has been a large range of studies on Water Mass\nTransformations22-28. Current estimations disagree on the dominance of isopycnal diffusion29\nor of diapycnal diffusion16,30 for the temperature and salinity transformations in the ocean\ninterior. However these studies rely on either numerical models or inverse method, both\nintrinsically depending on a chosen closure scheme (and sometimes imposed parameters).\nHence it is key to use direct observations to assess whether isopycnal or diapycnal diffusion is\ndominant in the ocean interior. Indeed, direct observations is the only way to validate or\ninvalidate hypotheses derived from these numerical or inverse modelling approaches. This is\nthe goal of the present study.”\nSec.[Conclusions] p.14. l.272-277: “In the literature there exist disagreement regarding the key\nrole of isopycnal diffusion in comparison to diapycnal diffusion for temperature and salinity\ntransformations17,29,30. Our direct estimates suggest that the dominance of isopycnal diffusion\nat 1,000-m depth. We hypothesize that this is because we fully acknowledged the isotropic\nnature of the isopycnal diffusion within our framework, unlike in inverse or numerical\nmodelling studies. This dominance is particularly important since isopycnal transformations\nare key for heat uptake in the context of Global Warming29.”\n(2) Most of the studies are concentrated on density transformations (not temperature and\nsalinity ones) or Water Masses volume changes. This makes a direct comparison difficult. To\novercome this, we have redone the analysis but using isopycnal diffusivity coefficients\nestimated from previous studies. These coefficients were computed using significantly\ndifferent methods, hypotheses, and observational datasets. For fair comparison on the\ntemperature and salinity transformation, we apply them to the same temperature and\nsalinity field, though.\nWe find that the three tested diffusivity coefficient led to the same qualitative and\nquantitative transformations (Fig. R4). Thanks to your suggestion this confirms the robustness\nof our analysis.\n14 /20\nFigure R4: Other estimations of isopycnal diffusivity and induced Temperature and Salinity water mass\ntransformations at 1,000-m depth. (a and b) Isopycnal diffusivity coefficients at 1,000-m estimated with\nFigualrteernSat3iv:e Ometthhoedrs3e5,s58t (idmenaottieodn ass Gorfoeisskoapmypc annad lRodaicffh,u rseisvpeitcytivealyn)d. Wiantderu mcaesds trTanesmfoprmeartaiotnusr feor a(cn) d\nSaliCnointsyervwaativtee rTemmpaesrastutrre aannds ffoorr m(d)a Atbisoonlustea Stalin1i,t0y.0 B0lu-me lindese sphtohw. tra(nasfaonrmdatbio)nIs sfoopr ythcen iasoltdroipffiuc sainvdit y\ncoeffispcaietinatllsy-avatr1yi,n0g0 0d-iffmuseivsittiym. Saotleidd awndit hdaashlteedr nblautei vlienems eshthowod rseYs,uZlts( dfoern tohtee davaesraGger obeestwkaemenp prainmdarRy oanadch ,\nrespseecctoivnedlayr)y. dWiffuastievirtiems aasnsdt sreacnosnfdoarrmy adtififounsisviftoiers,( rce)spCeoctnivseelryv. aOtriavnegTe eamndp reerda tliunrees sahnodw ftorran(sdfo)rAmabtsioolnus te\nfor the Groeskamp and Roach diffusivities, respectively. The curves have been smoothed using a moving\nSalinity. Blue lines show transformations for the isotropic and spatially-varying diffusivity. Orange\naverage of 0.06 Kg m-3 re-interpolated on regular 0.012 Kg m-3 grid with a spline interpolation. The horizontal\nand red lines show transformations for the Groeskamp and Roach diffusivities, respectively. The\ncolor patches show neutral density classes of selected Water Masses45.\ncurves have been smoothed using a moving average of 0.06 Kg m−3 re-interpolated on regular\n0.012 Kg m−3 grid with a spline interpolation. The horizontal color patches show neutral density\nclasses of selected Water Masses32.\nWe have added a figure (Fig. 6) and a text regarding this point.\nSec.[Robustenss…] p.11-12 l.210-242: “To test the robustness of our results, we have\nrecomputed the temperature and salinity transformations using two alternative diffusivity\ncoefficient datasets (see Methods – Alternative Diffusivity Coefficients). The first dataset, that\nwe named “Groeskamp” diffusivity coefficients, used a different set of observations to infer\ndiffusivity coefficients than our estimates32. Here the diffusivity coefficients are estimated\nthrough salinity observations assuming the validity of the mixing length theory. The second\ndataset, that we named “Roach” diffusivity coefficients, used the same deep observations\nthan our estimates, but a different methodology: a pairwise dispersion of deep Argo float\ndisplacements59.\nWhen comparing our diffusivity estimation (Fig. 1d and e) from Groeskamp and Roach ones\n(Fig. 6a and b, respectively), we find the same overall pattern, with intensification in western\nboundary currents and within the Antarctic Circumpolar Current. Because Groeskamp and\nRoach diffusivity coefficients correspond solely to the cross flow diffusion, they are more\ncomparable to our secondary diffusivity coefficients (Fig. 1e), which is by construction weaker\nthan the primary one. This secondary diffusivity coefficients are between 2,000 m2 s-1 and\n6,000 m2 s-1 in active regions. This is comparable to Groeskamp and Roach diffusivity\n3\n15 /20\ncoefficients. In comparison with the two other datasets, Roach diffusivity coefficients are less\nspatially-smooth.\nBecause Groeskamp and Roach datasets provide only a single local diffusivity coefficient, we\nrestrict the comparison to the spatially-variable and isotropic framework. Computing the\ntemperature and salinity transformations with the different diffusivity coefficient datasets\nshow the good consistency of the results (Fig. 6c and d). When using the spatially-variable and\nisotropic diffusivity coefficients, we find temperature and salinity transformations that are\nstronger than those obtained with Groeskamp and Roach diffusivity coefficients. However the\ntransformations are qualitatively well represented. As expected, the agreement is better when\napplying the spatially-variable and isotropic diffusivity coefficients solely based on the\nsecondary diffusivity coefficients. Here the agreement is qualitative and quantitatively good.\nThis stresses the robustness of our described temperature and salinity transformations with\nrespect to the methodology used to infer the diffusivity coefficients.\nThe good agreement of the three datasets regarding the salinity transformations (Fig. 6d) is\nparticularly interesting. Indeed Groeskamp diffusivity coefficients are based on an inversion\nmethod using salinity observations (i.e., somehow the diffusivity coefficients carry the\nknowledge of the salinity transformations). Hence the agreement of our analysis regarding\nsalinity transformations regardless of the diffusivity coefficient datasets confirms the\nrobustness of the overall methodology.”\nAlso, they claim to have “unveiled” these “hidden” isopycnal transformations, but they are\nhardly the first to make estimates of these transformations–even if they may have a case\nfor their estimate being the most sophisticated since it uses independent datasets with a\nspatially-variable eddy diffusivity. I recommend referencing papers that discuss the impact\nof isoneutral mixing in climate models (two recent examples:\nhttps://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2019GL085208 and\nhttps://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2021MS002914).\nThank you for these references, they are now cited in the manuscript. This give a wider\ncontext on the key role of isopycnal diffusion, beyond transformations, as tested in numerical\nmodels. We have toned-down our statements. We did not intend to imply being the first to\nhypothesize that there will be temperature and salinity transformations along isopycnal\nmixing. Instead, we pointed to the fact that our study is the first observational study\nquantifying the role of T and S isopycnal transformations. It is key to test numerical hypothesis\nthrough independent observations to validate or reject them. This is what we did.\nSec.[The Role of…] p.3 l.51-53: “Assessing the role of turbulent eddy diffusion for heat and\nfreshwater transports and for the transformation of seawater properties is a timely question\nof ocean and climate sciences, as demonstrated by numerous numerical model studies19-21.”\nThe novelty of our study resides in the use of two independent set of global observations to\ncompute the transformations: (1) the diffusivity coefficients and (2) temperature and salinity\nfields. Then we used these two independent estimate to deduce the temperature and salinity\n16 /20\ntransformations. This is in this sense that we observe for the first time and globally this\ncirculation.\nOn the other hand, the word “hidden” referred to a property of the transformations (i.e., that\nit does not change much the density, and so the pressure, hence the circulation):\ntransformations occur without (much) change in the circulation, making this transformations\n“hidden”. In simpler terms this circulation is isopycnal leading to no change in density (except\nthe “small” nonlinear effect). This does not mean that it was “hidden” from our knowledge.\nWe have clarified this point in the text to avoid misinterpretation.\nSec.[Effect of…] p.5-6 l.100-105: “This transformation-induced process has almost no\nsignature on the density field (except for the weak cabbeling and thermobaricity components,\nalthough they remain crucial for the density transformations52), nor pressure, nor circulation.\nIn this sense it remains “hidden” along the isopycnals (Fig. 1c). However it still contributes to\nheat and freshwater transport from cold to warm and from salty to fresh regions,\nrespectively.”\nThe roles of eddy-bolus velocities (or the skew-symmetric part of the diffusivity tensor) and\nthe diapycnal velocity in the ocean & climate system are unjustifiably downplayed to make\ntheir isoneutral mixing sound more important than it is. The authors in some places forget\nthe qualifier that they are only discussing the isopycnal mixing role of eddies and make\nstatements like “since eddies induce a flow along isopycnals, they cannot modify density…\nhence the circulation”, which are not true if one always considers the leading-order quasi-\nadvective role that eddies play in the general circulation.\nThank you, we have clarified this key point. Note that in the Temporal-Mean-Framework, at\nequilibrium and away from forcing, the flow is indeed along isopycnals (i.e., 𝐷# 𝜌& = 0 ⇒ 𝐮+ ∙\n$\n𝛁𝜌& = 0). However, in the specific sentence commented by the reviewer, we were referencing\nto eddy mixing. We have modified the sentence to avoid confusion:\nSec.[Conclusions] p.13 l.249-251: “since eddy turbulence corresponds to a flow along\nisopycnals, its induced mixing cannot modify density […], hence the circulation.”\nWe have also further acknowledged the role of diapycnal mixing potentially at other depths\nand for density transformations at our studied depth diapycnal\nSec.[Conclusions] p.14 l.282-285: “However, given the intensity of the circulation and\ntransformation identified in our study, it is important to acknowledge it at other depths. In\nparticular, diapycnal mixing could become more prevalent near the seabed, mitigating the key\nrole of isopycnal mixing for temperature and salinity transformations.”\nRegarding the eddy-bolus velocities, since they do not represent mixing but advection, they\ndo not induce transformations, thus, we have disregarded them in the present study. We\nagree that they are key for others aspects of ocean and climate dynamics, but this is out of\nthe scope of the paper that focuses on transformations.\n17 /20\nAlso, the comparison with diapycnal mixing transformations (19x smaller!) feels\ndisingenuous; the reason people usually talk about diapycnal mixing is not because how it\ntransforms the temperature or salinity of water masses, but how it does so diapycnally in a\nway that no other process in the ocean interior can.\nThis is correct. We have further clarified our statement anout isopycnal mixing being 10 times\nas efficient as diapycnal mixing for temperature and salinity transformations at our studied\ndepth (~1000 dbar). This is not true for density transformations, which are not studied here.\n(Indeed, diapycnal mixing is the dominant over isopycnal mixing for density transformation,\nexcept for some specific Water Masses, as shown by Groeskamp et al, Geophys. Res. Lett.,\n2016, for instance). This has now explicitly discussed in the manuscript.\nSec.[Water Masses…] p.9 l.169-174: “We find that the temperature and salinity isopycnal\ntransformations are 13 and 8 times as large as the diapycnal ones (black vs cyan line in Fig. 5c\nand d), respectively. Note that this result contrasts with density transformations for which, as\nexpected, diapycnal transformations is the dominant factor over isopycnal ones, except for a\nfew specific Water Masses (i.e., Subantarctic Mode Water and AAIW, and CDW), for which\ncabbeling and thermobaricity are important factors52.”\nSome of the key points are poorly articulated and difficult to decipher. In most cases I could\neventually decipher what you are trying to say, but not always (see minor comments).\nFollowing your suggestions and main point below, we have significantly revised the text. Also,\nplease see below for modifications of the text following your minor comments. We hope that\nthis new version of the manuscript is easier to follow.\nClarity of communication\nFinally, the braggadocious tone of the writing is off-putting and I encourage the authors to\ntone it down. The entire manuscript is written as though the authors are the first to realize\nthat eddies are rapidly mixing temperature and salinity along isopycnals, or at least the first\nto make an estimate of it.\nWe understand your comment. We have done our best to tone-down the writing. As\nexplained above, and to the best of our knowledge, we are the first to make a global (even if\nrestricted to 1,000-m depth) independent observational estimate of these transformations.\nHowever we never claim that we are the first to hypothesize their existence.\nSome other examples: 1) Water mass transformation is a framework that is more than 50\nyears old, so calling it “state-of-the-art” is a tough sell.\n18 /20\nWe have clarified the text. The “state-of-the-art” refer to the Water Mass Transformation\ntogether with Temporal-Residual-Mean or Thickness-Weighted Average that, even if 10 year\nold, are often disregarded when making such key diagnostic.\nSec.[Abstract] p.1 l.10-13: “Here we use a robust framework, combining the Water Mass\nTransformation together with the Temporal-Residual-Mean, and a novel direct estimation of\nturbulent anisotropic diffusivity (computed through an innovative artificial intelligence\nmethod), to determine temperature and salinity transformations along density surfaces.”\n2) The other “state-of-the-art” framework is the authors’ novel AI-based method for\nestimating turbulent anisotropic diffusivities. It is indeed very interesting and–to my\nknowledge–a useful novel approach, but the authors don’t actually go any further in\ndescribing it in this manuscript than stating its name (the “Analog method”) and citing their\nprior work!\nThank you for this comment. We have clarified this point in the new version of the\nmanuscript. In particular, now that we are testing three different estimates, it is even more\ncrucial to describe them. This is now done in the manuscript.\nSec.[The role of…] p.4 l.65-69: “[…] we have used the recently estimated diffusivity coefficients\nderived from Argo float deep displacements34. This dataset has been obtained using a class of\ndata driven Artificial Intelligence method: the Analog Method35,36. There, the spatially-\nvariable and anisoptropic diffusivity coefficients was estimated through local tracer dispersion\ninduced by a turbulent velocity field reconstructed from deep Argo displacements.”\nMethods.[Tracer Fields,…] p.30 l.549-552: “The spatially-variable and anisotropic diffusivity\ncoefficients were obtained by reconstructed an analogue turbulent field from deep Argo\ndisplacement. This reconstructed turbulent field was used to determine the local dispersion of\na locally-injected spike of passive tracer, the rate of dispersion estimating the local diffusivity\ncoefficients along the zonal, meridional, and cross direction.”\n3) Given how impactful this data set is in this manuscript, it deserves at least a paragraph\nin the method section. (There is plenty of fluff in the introduction and conclusions that you\ncould cut out to make space).\nThank you. As stated above it is now clarified and thanks to your comments we have tested\ntwo other estimates (which lead to consistent results).\nMinor comments:\nEq. (4) – I understand where the authors are coming from, but it doesn’t really make sense\nto decompose the transformations in this way because isopycnal + diapycnal already\naccounts for all directions. Either choose a two-way decomposition (isopycnal vs. diapycnal)\nor explain that you are breaking these down by process: mixed-layer turbulence (strong in\nall directions), interior mesoscale turbulence (isopycnal), and small-scale turbulence (in all\n19 /20\nthree directions, but the isopycnal component is swamped out by the mesoscale isopycnal\nso can be approximated by diapycnal).\nThank you for this comment, we agree that this could be confusing. We were following the\nsuggested split of previous studies. Following your advice we have revised the split and only\nkept isopycnal + diapycnal (removing horizontal), since our study is always below the ocean\nmixed layer. Please see new equation (4):\n𝑑 = 𝑑123 +𝑑516, (4)\n/0 /0 /0\nFigure 3. Colormaps are not very interpretable and I think very bad for people with\ncolorblindness.\nThe colormap were designed for good interpretability for people with the most common\ntypes of colorblindness. If this is not the case, we will be happy to change them.\nFigure S1. Why not convert both to Sv so they can be directly compared? “For the isotropic\ndiffusivity” I think should read “For the anisotropic diffusivity”. Caption should clarify these\nare “Water mass transformations in density space” or something.\nWe agree that this would be ideal. Unfortunately, the units cannot be converted in Sv. They\nare the units from (11), when the T and S transformations are computed for density classes,\nwhich results in these rather complicated units.\nThank you we have modified the caption accordingly to your comment.\nL. 374-375: I don’t understand this.\nWe are simply restating that isopycnal transformation are by definition along constant\ndensity, and so by construction avoid density gradient. We have rephrase this sentence to be\nmore explicit:\nMethods.[Isopycnal Turbulent…] p.25 l.475-476: “Note that isopycnal heat and freshwater\ntransformations act along constant density surface, thus, by construction they avoid density\ngradients.”\nWhy epineutral rather than isoneutral? Is epineutral more specifically a geometric\ndescriptor?\nThis is the accepted terminology in the literature. To the best of our understanding this\nencompasses the idea that it acts along constant value of neutral-density. Note that “iso-\nneutral” is somehow odd because it suggests the juxtaposition “same” (iso) and “no change”\n(neutral). Hence epineutral is a bit more specific that isopycnal since it encompasses the\ndefinition of neutral density (i.e., isopycnal could be defined for potential density, for\ninstance).\n20 /20"}, "author_response": "Response to the reviews of the manuscript\nTurbulent Isopycnal Mixing dominates Thermohaline Transformations of Intermediate\nOcean Waters\nby F. Sévellec, N. Kolodziejczyk, and E. Portela\nWe thank the two reviewers and the editor for their constructive comments, which have\nhelped us improving and clarifying the manuscript. Our detailed response is below.\nFont legend:\n1. Reviewer’s comments are in bold,\n2. Our responses are in normal font, and\n3. Text added in the manuscript is in italic.\nReviewer #3 (Remarks to the Author):\n2nd Review of Turbulent Isopycnal Mixing dominates Thermohaline Transformations of\nIntermediate Ocean Waters by Florian Sevellec, Nicolas Kolodziejczyk, and Esther Portela\nI appreciate the improvements the authors have made. Unfortunately, I still don’t feel there\nis enough clarity of the methods, embedding of the subject in literature nor clarity of the\npresentation that the manuscript is ready for publication. I recommend major revisions.\nThank you for your supportive comment. We hope that the modifications make this new\nversion now ready for publication.\nNovelty\nI’m now more enthusiastic about the work than I was before. I think the idea of comparing\nthe isopycnal transformation of salt and temperature (which have received relatively little\nattention) to their diapycnal counterparts, a worthwhile exercise. However, the reason it\nhas received less attention may also be because it is less interesting. The authors are trying\nto make the case that this is in fact very interesting, but I’m not yet convinced based on\nthese results. I’m still open to be proven otherwise. Below I will continue to make myself\nclear.\nAgain, thank you for this overall supportive comment. For initial clarification (below you will\nfind our point-by-point specific response), we would like to say that we do not want to\ndemonstrate that thermal- and/or haline- transformations are more important than density-\ntransformations. We are simply concentrating on these transformations because they are\nalso a key aspect of the ocean circulation. This is important in the context of heat and\nfreshwater transport (though it has little effect on overturning). Indeed they are key aspects\nof the ocean circulation that are currently under stressed within climate changes. We agree\nwith the reviewer that they have somehow received less attention.\n1 /20\nThe main result as presented\nBecause the ocean does move along and across neutral surfaces, so does T&S. To get local\nchanges in T or S we thus have to account for local air-sea-ice interaction and mixing\nprocesses (and maybe geothermal heating). These mixing processes are either mesoscale\nhorizontal (in mixed layer), mesoscale isoneutral (in deep ocean) or vertical. The results can\nbe used to calculate WMT in T, S or density coordinates. In the S&T framework this leads to\nlarger transformation (Sv) from one temperature to another (or salinity) due to isoneutral\nmixing than due to vertical mixing. That is the main finding of this paper. This result,\nhowever, is not put in context of other studies and is wrongly compared to density\noverturning. \nWe agree with the general summary of our work. Regarding the last part, we have improved\nthese aspects in the manuscript. Please see the detailed response below.\nNote that comparing density-transformation with thermal- and haline- transformation is\npossible. Higher value of the latter over the former does demonstrate that the ocean is also\nsignificantly “spicy”. Hence, solely concentrating on density transformation “hides” important\nheat and freshwater transports. This is what we try to demonstrate quantitatively.\nPutting this results in perspective\nFirst, other literature on this subject is provided by the vast amount of old and new studies\nthat have looked at transformation in T&S coordinates. This is essentially what has been\ndone here as well as this is exactly diahaline and diathermal transformation by these\ndifferent mixing processes. Many diagrams have been shown and never have I seen a factor\n20 difference. Maybe I missed something, but that should then be made clear by a proper\nembedding into this literature. Below is a list of references for this subject.\nThank you for this list of references. They are now included in the manuscript. We did not see\nany inconsistency with previous work. We have revised our estimate by acknowledging the\nspatial variation of diapycnal/vertical diffusion coefficient. We now have a factor 12. The fact\nthat our computation is solely located at a single depth can also explain such a large value\n(i.e., other depth may have weaker value).Hence, we still find that at 1000-m depth the\ntemperature and salinity transformations are widely dominated by isopycnal (but diathermal\nand diahaline) processes over diapycnal processes.\nWe have now acknowledged the literature in the introduction and used it to develop our\nscientific question. Also we used it to provide further context to our conclusion.\nSec.[The Role…] p.3 l.54-61: “There has been a large range of studies on Water Mass\nTransformations23-28. Current estimations disagree on the dominance of isopycnal diffusion29\nor of diapycnal diffusion67,30 for the temperature and salinity transformations in the ocean\ninterior. However these studies rely on either numerical models or inverse method, both\nintrinsically depending on a chosen closure scheme (and sometimes imposed parameters).\nHence it is key to use direct observations to assess whether isopycnal or diapycnal diffusion is\ndominant in the ocean interior. Indeed, the use of direct observations is the only way to\nvalidate or invalidate hypotheses derived from these numerical or inverse modelling\napproaches. This is the goal of the present study.”\n2 /20\nSec.[Conclusions"}]

Decision letter