# Global heat balance and heat uptake in potential temperature coordinates

Antoine Hochet, Rémi Tailleux, Till Kuhlbrodt, David Ferreira

► To cite this version:

Antoine Hochet, Rémi Tailleux, Till Kuhlbrodt, David Ferreira. Global heat balance and heat uptake in potential temperature coordinates. Climate Dynamics, 2021, 57, pp.2021-2035. 10.1007/s00382-021-05832-7. insu-03683254

HAL Id: insu-03683254

https://insu.hal.science/insu-03683254

Submitted on 31 May 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Global heat balance and heat uptake in potential temperature coordinates

Antoine Hochet<fnref n="1" /> · Rémi Tailleux<fnref n="2" /> · Till Kuhlbrodt<fnref n="2" /> · David Ferreira<fnref n="2" />

Received: 21 November 2020 / Accepted: 27 April 2021 / Published online: 13 June 2021
© The Author(s) 2021

## Abstract

The representation of ocean heat uptake in Simple Climate Models used for policy advice on climate change mitigation strategies is often based on variants of the one-dimensional Vertical Advection/Diffusion equation (VAD) for some averaged form of potential temperature. In such models, the effective advection and turbulent diffusion are usually tuned to emulate the behaviour of a given target climate model. However, because the statistical nature of such a “behavioural” calibration usually obscures the exact dependence of the effective diffusion and advection on the actual physical processes responsible for ocean heat uptake, it is difficult to understand its limitations and how to go about improving VADs. This paper proposes a physical calibration of the VAD that aims to provide explicit traceability of effective diffusion and advection to the processes responsible for ocean heat uptake. This construction relies on the coarse-graining of the full three-dimensional advection diffusion for potential temperature using potential temperature coordinates. The main advantage of this formulation is that the temporal evolution of the reference temperature profile is entirely due to the competition between effective diffusivity that is always positive definite, and the water mass transformation taking place at the surface, as in classical water mass analyses literature. These quantities are evaluated in numerical simulations of present day climate and global warming experiments. In this framework, the heat uptake in the global warming experiment is attributed to the increase of surface heat flux at low latitudes, its decrease at high latitudes and to the redistribution of heat toward cold temperatures made by diffusive flux.

**Keywords** Heat uptake · Simple Climate Model

## 1 Introduction

Ocean heat uptake is of great importance in climate change predictions: 90 % of the anthropogenic increase in heat stored in the climate system ends up in the oceans (Levitus et al. 2012), thus contributing to sea level rise via thermal expansion. The main effects controlling the heat balance include the upwelling of deep water driven by the Southern Ocean winds, cooling by deep water formation, as well as isopycnal and diapycnal mixing, most of which require

to be parameterized in current AOGCMs (see for instance Marshall and Zanna (2014) and references therein). The ocean heat uptake efficiency, defined as the ratio of net heat flux into the climate system over the global mean surface air temperature change (Gregory and Mitchell 1997; Raper et al. 2002), has been found to vary by a factor of 2 across CMIP5 models (Kuhlbrodt and Gregory 2012) outlining its high sensitivity to the parametrization choices made by the various modelling groups. A better understanding of the heat balance processes will thus also help constrain mixing parameters in these models.

One of the most common method to rationalize the heat balance in the ocean consist in studying its vertical structure from its horizontally-averaged properties. This method is justified by its simplicity but also by the interest in the vertical structure of the temperature which is linked with the idea of ocean heat storage.

The vertical heat transport described by the horizontally-averaged heat balance is often compared with the theory of early models of the deep circulation such as Wyrtki (1961)

<footnote>Rémi Tailleux
r.g.j.tailleux@reading.ac.uk

Antoine Hochet
antoine.hochet@univ-brest.fr</footnote>

<footnote n="1">¹ CNRS, Ifremer, IRD, Laboratoire d'Océanographie Physique et Spatiale (LOPS, UMR 6523), Univ of Brest, IUEM, Brest, France</footnote>

<footnote n="2">² Department of Meteorology, University of Reading, Reading, UK</footnote>
where dense water downwells at high latitude due to convection in very localised regions and upwells uniformly at mid and low latitudes. This leads to the classical view where the upwelling of cold/dense waters is balanced by downward diffusion of heat. This model is frequently referred to as the one-dimensional Vertical Advection/Diffusion (VAD hereafter) model in the literature. So far, however, it has proved difficult to reconcile the classical view of heat balance offered by the VAD model with that resulting from numerous studies of the horizontally-averaged heat balance such as Gregory (2000), Wolfe et al. (2008), Kuhlbrodt et al. (2015) and Exarchou et al. (2015). Indeed, in such studies the horizontally averaged advective heat fluxes are often found to be downward and the horizontally averaged diffusive or eddy-resolved heat fluxes (thus the average of a combination of iso and diapycnal diffusive fluxes) upward, which is seemingly the opposite of what the standard VAD model predicts (Gregory 2000).

Yet, the VAD model appears nevertheless successful at emulating the temperature variations of complex AOGCM (Raper et al. 2001). As a result, the VAD model has formed the basis for the one-dimensional representation of ocean heat uptake in Simple Climate Models (SCMs) such as MAGICC (Meinshausen et al. 2011). SCMs are used for instance to evaluate the amount of CO₂ that can be released in the atmosphere before reaching the 2°C limit (Meinshausen et al. 2009) and play an important role in policy making decisions about global warming mitigation strategies.

To reconcile these two approaches, Huber et al. (2015) proposed to calibrate the VAD equation (i.e. the set-up of vertical velocity $w$ and diffusive coefficient $K$) using a physical approach rather than the behavioral approach used in previous studies such as Raper et al. (2001). The two approaches differ in that the behavioural approach calibrates the VAD model parameters to mimic the temperature variations of complex AOGCMs using statistical techniques, whereas the physical approach seeks to calibrate such parameters by linking them to the processes that control them.

However, when horizontal averaging is used as the underlying basis for the physical calibration, the diffusion coefficient can occasionally be negative owing to the heat diffusion being occasionally upward in parts of the ocean. Moreover, the time variation of $K$ and $w$ were found crucial in emulating correctly the temperature of AOGCM thus complicating the practical implementation of the method. We have thus identified the two following points: (1) the possibility to justify the VAD model from horizontally-averaging the three-dimensional advection/diffusion equation for heat is far from obvious; (2) the occasional up-gradient nature of the horizontally-averaged heat flux complicates the construction of a one-dimensional VAD model because it does not act to reduce the

vertical temperature gradient as is expected physically. To circumvent this difficulty, we adopt a different approach: instead of averaging on constant depth surfaces we average on constant potential temperature ($\theta$ hereafter) surfaces, following an approach similar to that recently developed by Hochet et al. (2019). The averaged diapycnal diffusion is then automatically downgradient and we will further show that the advection through $\theta$ surfaces is in theory zero, leading to a much simpler equation than that obtained with constant depth surfaces. Ferrari and Ferreira (2011) has used a similar approach to study the ocean heat transport in order to filter out any recirculation of waters at constant temperature. Holmes et al. (2018) also used a similar approach to study the diathermal heat transport in a global ocean sea ice model.

The heat balance averaged in temperature coordinates can be expected to be quite different from the well studied horizontally averaged heat balance in depth coordinates. Indeed, because nearly all isotherms outcrop at the ocean surface, heat fluxes through the coldest temperature classes may either reflect processes at great depth or at high latitudes. In the standard VAD model heat fluxes through the coldest horizontally averaged temperature only pertain to processes at great depth. It might be useful to keep in mind the results from horizontal averages of AOGCM outputs in Control Run (CR hereafter) with constant present day CO₂ concentration and warming climates (see for instance Gregory 2000; Huang et al. 2003; Brierley et al. 2010; Kuhlbrodt et al. 2015). In CR, the strongest downward heat transport comes from the mean advection while the largest upward heat transport comes from eddy induced advection (resolved or parametrized). In warming climates, the heat uptake takes place mostly in the Southern ocean and is due to the reduction of along-isopycnal mixing and of deep convection. We analyse the outputs of the ocean component of the HiGEM1.2 coupled atmosphere ocean general circulation model (AOGCM), which include a detailed set of temperature tendency diagnostics. HiGEM1.2 is a CMIP5-type model and this study thus contributes to the understanding of heat uptake in this class of models. To analyse the processes controlling ocean heat uptake, we study the heat balance in temperature coordinates first in a control run of the HiGEM model that we then compare to a warming climate run where the pre-industrial CO₂ has been doubled.

The article is organized as follows: in Sect. 2, we derive an alternative one dimensional equation of heat uptake using potential temperature coordinates and show that it allows to remove the effect of advection and to obtain a downgradient diffusion. In Sect. 3, we apply this new method to the study heat uptake first in the CR of HiGEM, then on a simulation where CO₂ concentration is doubled. The last section concludes and discusses the results.

<figure><img src="image_2.png" /></figure>
## 2 Method

Because AOGCM outputs are generally averaged over a period of time (1 month here) all terms of the temperature budget are decomposed into time mean and anomalies:

$$X = \bar{X} + X' \qquad (1)$$

where X represents any term of the heat budget, $\bar{X}$ the monthly average and $\bar{X}'$ the deviation from this monthly average so that $\bar{X}' = 0$. The time mean potential temperature $\bar{\theta}$ conservation can be written as:

$$\frac{\partial \bar{\theta}}{\partial t} + \bar{v} \cdot \nabla \bar{\theta} + \bar{\nabla} \cdot \mathbf{v}' \theta' = \nabla \cdot (\mathbf{K} \nabla \bar{\theta}) + \bar{V} \mathbf{M} + \bar{Q}_{\text{net}} \qquad (2)$$

For clarity we will drop the overline notation in what follows and keep it only when it involves anomalies. $\mathbf{v}$ is the 3D velocity vector, $\bar{\nabla} \cdot \mathbf{v}' \theta'$ is a term representing the effect of sub-monthly advection, $\mathbf{K}$ a diffusion tensor representing the effect of unresolved advection and small-scale irreversible mixing, $\mathbf{K} \nabla \theta \cdot \mathbf{k}$ with $\mathbf{k}$ the upward unit vector is thus zero at the surface. $\mathbf{K}$ thus contains the parameterization of both the isopycnal and diapycnal mixing terms. VM is a term representing all parameterized non-diffusive terms like convection and $Q_{\text{net}}$ the net heat flux through the surface. For comparison, the equations of the physical calibration of the VAD using the horizontal average of Eq. (2) are derived in Appendix A. Building on Winters and D'Asaro (1996)'s work, we first define a reference level $z_r$ of the temperature $\theta$. The use of a reference level will be useful to obtain an 1D evolution equation for the temperature along surfaces of constant reference depth as will become clear below. $z_r$ is the depth of isotherm $\theta$ in the reference state which is obtained

after an adiabatic rearrangement of each fluid parcel so that isotherms are horizontal and in ascending order. Note that unlike the reference state described in Winters and D'Asaro (1996), this reference state is not a state of rest because the density is here also a (non-linear) function of salinity and pressure. Such a rearrangement being volume conserving, the reference depth $z_r$ can thus be computed using the fact that the volume of water with temperature larger than $\theta$ is the same after and before the adiabatic rearrangement i.e.:

$$\int_{z_r}^{0} A(z) dz = \int_{V(\theta,t)} dV, \qquad (3)$$

where $V(\theta, t)$ is the volume of ocean with temperature $\theta_l$ satisfying $\theta < \theta_l < \theta_{max}$ with $\theta_{max}$ the maximum temperature in the ocean and $A(z)$ is the ocean area at depth $z$. The definition (3) of the reference depth makes it possible to rewrite the temperature $\theta(x, y, z, t)$ as a function of $z_r$: $\theta_r(z_r, t) = \theta(x, y, z, t)$. $\theta_r$ can be inverted to yield $z_r = z_r(\theta, t)$ or $z_r = z_r(x, y, z, t)$. Note that Eq. (3) shows that the volume $V(\theta, t) = V(z_r)$ of water of temperatures greater than $\theta$ is a function of $z_r$ alone and hence that it can be treated as a constant independent of time at fixed $z_r$. An alternative definition of $z_r$, that can be found for instance in Winters and D'Asaro (1996), is:

$$z_r(x, y, z, t) = \int_V \mathcal{H}[\theta(\tilde{x}, \tilde{y}, \tilde{z}, t) - \theta(x, y, z, t)] d\tilde{V} \qquad (4)$$

where $\mathcal{H}$ is the Heaviside step function and $V$ represents the ocean volume. The schematic shown on Fig. 1 summarizes the calculation of the reference depth as explained above.

<figure><img src="image_3.png" /><figcaption>**Fig. 1** Schematic showing how the reference depth $z_r$ associated with temperature $\theta$ is obtained. On the left is the physical space, on the right the reference space which is obtained through an adiabatic rearrangement of all fluid parcel in the physical space. Isotherms in the reference space are horizontal and only depend on $z_r$, temperature in the reference space is described by the function $\theta_r(z_r, t)$. The volume</figcaption></figure>

for all water parcels warmer than $\theta(X, t) = \text{const.}$ is shown by black stripes in both physical and reference space. This volume is the same in both spaces, and this property is used in formula (3) with $A(z)$ the ocean area at depth $z$ to compute the reference depth associated with the temperature
We now seek an evolution equation for $\theta_r(z, t)$ by integrating (2) over the volume $V(z_r)$, which after some manipulation yields:

$$ \int_{V(z_r)} \frac{\partial \theta_r}{\partial t} dV + \theta_r(z_r, t) \int_{z_r=\text{const.}} \mathbf{v} \cdot \mathbf{n} dS + \int_{V(z_r)} \nabla \cdot \overline{\mathbf{v}'\theta'} dV \\ = \int_{z_r=\text{const.}} \mathbf{K} \nabla \theta \cdot \mathbf{n} dS + \int_{V(z_r)} \text{VM} dV + \int_{V(z_r)} Q_{\text{net}} dV, \quad (5) $$

where $\mathbf{n} = -\nabla\theta/|\nabla\theta| = -\nabla z_r/|\nabla z_r|$ is the outward unit normal vector to the isothermal surface $\theta = \text{constant}$, which at fixed time coincides with the surface $z_r(x, y, z, t) = \text{constant}$. The diffusion term can be written as:

$$ \int_{z_r=\text{const.}} \mathbf{K} \nabla \theta \cdot \mathbf{n} dS = -K_{\text{eff}}(z_r, t) A(z_r) \frac{\partial \theta_r}{\partial z_r}(z_r, t) \quad (6) $$

with $K_{\text{eff}} = \frac{1}{A(z)} \int_S K_{\text{eff}}^{\text{loc}} dS'$ and $K_{\text{eff}}^{\text{loc}}$ a positive quantity independent of $\partial_{z_r} \theta_r$, indeed:

$$ \mathbf{K} \nabla \theta \cdot \mathbf{n} = -\frac{\partial \theta_r}{\partial z_r} (K_i(\nabla z_r - (\nabla z_r \cdot \mathbf{d})\mathbf{d}) + K_d(\nabla z_r \cdot \mathbf{d})\mathbf{d}) \cdot \frac{\nabla z_r}{|\nabla z_r|} \\ = -\frac{\partial \theta_r}{\partial z_r} \underbrace{(K_i \sin^2(\nabla z_r, \mathbf{d}) + K_d \cos^2(\nabla z_r, \mathbf{d}))}_{K_{\text{eff}}^{\text{loc}}} |\nabla z_r| \quad (7) $$

where $K_i$ and $K_d$ are the isoneutral and dianeutral turbulent diffusivities respectively. Using the non-divergence of the velocity field and neglecting the contribution of the freshwater fluxes (whose expression is derived in Appendix B) so that $w = 0$ at the surface, we have:

$$ \int_S \mathbf{v} \cdot \mathbf{n} dS = \int_{V(z_r)} \nabla \cdot \mathbf{v} dV = 0. \quad (8) $$

This equation holds even under a non-steady state and means that a closed volume cannot increase or decrease due to advection by a non-divergent velocity through its boundaries. The more general case for which $w(z=0) = E - P + R$ with $E, P$ and $R$ respectively the evaporation precipitation and river run-off is discussed in details in Hochet and Tailleux (2019) and described briefly in Appendix B. Using Eqs. (6) and (8) in Eq. (5) gives:

$$ \int_{V(z_r)} \frac{\partial \theta_r}{\partial t} dV = -K_{\text{eff}}(z_r, t) A(z_r) \frac{\partial \theta_r}{\partial z_r}(z_r, t) \\ - \int_{V(z_r)} \nabla \cdot \overline{\mathbf{v}'\theta'} dV + \int_{V(z_r)} \text{VM} dV + \int_{V(z_r)} Q_{\text{net}} dV. \quad (9) $$

This equation links the volume integral on $V(z_r)$ of the time derivative of the temperature to the diffusive flux, the sub-monthly advection, the vertical mixing and the surface heat

flux. Calculating the derivative of Eq. (9) with respect to $z_r$ and dividing by $A(z_r)$ gives an evolution equation for $\theta_r$:

$$ \frac{\partial \theta_r}{\partial t} = \frac{1}{A(z_r)} \frac{\partial}{\partial z_r} \left( K_{\text{eff}}(z_r, t) A(z_r) \frac{\partial \theta_r}{\partial z_r} \right) \\ - \frac{1}{A(z_r)} \frac{\partial}{\partial z_r} \int_{V(z_r)} \nabla \cdot \overline{\mathbf{v}'\theta'} dV - \frac{1}{A(z_r)} \frac{\partial}{\partial z_r} \int_{V(z_r)} \text{VM} dV \\ - \frac{1}{A(z_r)} \frac{\partial}{\partial z_r} \int_{V(z_r)} Q_{\text{net}} dV, \quad (10) $$

where we have used:

$$ \frac{\partial}{\partial z_r} \int_{V(z_r)} \frac{\partial \theta_r}{\partial t} dV = \frac{\partial}{\partial z_r} \int_{z_r}^{0} A(z) \frac{\partial \theta_r(z, t)}{\partial t} dz \\ = -A(z_r) \frac{\partial \theta_r(z_r, t)}{\partial t} \quad (11) $$

The possibility to obtain a 1D equation for $\theta_r(z_r, t)$ as given by Eq. (10) is one of the main advantage of the use of a reference level. Note that Eq. (10) is similar to Eq. (17) in Winters et al. (1995) with vertical mixing, sub-monthly advection and heating terms added. In agreement with Hieronymus et al. (2014), Eqs. (9) and (10) establish that the time evolution of the reference potential temperature is only a function of the effective diffusion, of the sub-monthly advection, of the forcing and of the vertical mixing. The (resolved) monthly advection does not play any role in the evolution of $\theta_r$ and the diffusive part is only due to the divergence of the downgradient diffusive flux. In the remaining of this paper we use Eq. (9) and (10) to study heat uptake in the Control Run and $2 \times \text{CO}_2$ run of a climate model.

# 3 Results

## 3.1 Model

HiGEM1.2 is an AOGCM pertaining to the CMIP5-type models. It is based on the UK MetOffice coupled AOGCM HadGEM1, but has a higher spatial resolution, of 0.83° lat. × 1.25° lon. (N144) in the atmosphere and 1/3° × 1/3° with 40 levels in the ocean. An implicit linear free surface scheme based on Dukowicz and Smith (1994) with explicit fresh water fluxes is used. Lateral mixing of tracers uses the isopycnal formulation of Griffies (1998), and the Gent and McWilliams (1990) (GM) adiabatic mixing scheme is not used. A detailed description of this model can be found in Shaffrey et al. (2009). We use two different runs of HiGEM1.2: (1) a Control Run (CT hereafter) where present-day boundary conditions are used, in particular, the atmospheric CO$_2$ concentration is set to 345 ppm, reflecting conditions in the 1980s, and (2) a perturbed run where
atmospheric CO₂ concentration is doubled (2× CO₂). The control run length used in this article is 50 years and the 2× CO₂ perturbation run length is 70 years.

The HiGEM diagnostics used here consist in monthly means of the potential temperature tendencies i.e. all terms at each grid point contributing to local changes in potential temperature. These terms comprise potential temperature change due to advection, diffusion (separately in the x, y and z directions), convection, mixed layer physics, ice physics, penetrating solar radiation and other surface fluxes. Note that there is no GM parameterisation, the advection diagnostic thus contains both the mean and resolved eddy-induced advection. We regroup in what follows convection and mixed layer dynamics into a vertical mixing (VM) term and penetrative solar, surface fluxes, ice physics into a forcing term. We are thus left with four terms: diffusion, advection, vertical mixing and forcing.

As seen from Eq. (5) the integral is performed on volumes defined by surfaces of constant $z_r$. For each time t, $\theta$ = const. surfaces are exactly the same as the $z_r$ = const. surfaces. However, the fact that $\theta_r(z_r, t)$ is also a function of time implies that the temperature associated with a given reference level is time dependent. Practically it means that we need to calculate the reference level for every monthly mean outputs and then perform the volume integral of the tendencies. The method used to calculate the volume integral of the heat tendencies is described in Appendix C. The reference levels and volume integral of heat tendencies are calculated for monthly means for both the Control Run and the 2×CO₂ run. They are then averaged over a 50 years period for the CR and on the 70 years of the 2×CO₂ run.

## 3.2 Control Run

### 3.2.1 Reference level

The 50 years mean reference level is shown on the left panel of Fig. 2. As expected, it is a monotonic function of temperature, deepest (shallowest) $z_r$ correspond to coldest (warmest)

temperatures. Because most of the volume of the ocean has small temperatures below 5 °C, the range of temperature between -1000 and 0 m is much larger (~ 25 °C) than at deeper depth: ~ 7 °C between -5500 and -1000 m.

The reference temperature gradient will therefore be much larger at shallow reference depth than at deep reference depth. The ocean area as a function of depth A(z) calculated for the HiGEM grid and used in the reference depth calculation (see formula 3) is shown on the right panel of Fig. 2.

### 3.2.2 Time mean of the volume integral of the heat tendencies as a function of the reference depth

At each grid cell, heat tendencies are decomposed using the following equation:

$$ \frac{\partial \theta}{\partial t} = \text{advection} + \text{diffusion} + \text{VM} + \text{forcing} \quad (12) $$

where “advection”, “diffusion”, “VM”, “forcing” are respectively the three dimensional heat tendencies due to advection, diffusion, VM, and forcing described in the last section. Equation (12) is then integrated on volume $V(z_r)$ described in Sect. 2:

$$ \begin{aligned} & \int_{V(z_r)} \frac{\partial \theta}{\partial t} dV' \\ &= \int_{V(z_r)} \text{advection} dV' + \int_{V(z_r)} \text{diffusion} dV' \\ & \quad + \int_{V(z_r)} \text{VM} dV' + \int_{V(z_r)} \text{forcing} dV' \end{aligned} \quad (13) $$

Figure 3 shows the time mean volume integral of the heat tendencies as a function of the reference level for the CR.

The time mean of the integral of the tendencies is negative at all the reference depth for the diffusion, advection and vertical mixing and always positive for the forcing. This shows that for all $z_r$, diffusion, advection and vertical

<figure><img src="image_4.png" /><figcaption>**Fig. 2** Left panel: time mean of the reference temperature (in °C) in the CR as a function of the reference depth in meters. Right panel: surface of the ocean A(z_r) normalized by its maximum value as a function of the reference depth</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 3 Volume integral of heat tendencies (K year⁻¹) on V(z_r) associated with vertical mixing (convection + mixed layer dynamics) in blue, forcing (surface fluxes + ice) in orange, advection in red, diffusion in green as a function of the reference depth z_r in meters. The time mean of the temperature as a function of z_r is shown on the right. The sum of all terms is shown in purple. The time mean of the temperature (over the CR) as a function of z_r is shown on the right</figcaption></figure>

mixing act together to reduce the temperature of the volume of water parcels with z'_r larger than z_r while the forcing acts to increase it. Figure 3 is similar to Fig. 3 of Holmes et al. (2018) where the budget for the internal heat content of a global ocean sea ice model is expressed in terms of surface forcing, vertical mixing and “numerical” mixing (which is calculated as a residual and thus contains the isopycnal mixing). Our forcing term looks similar to theirs, the detailed comparison for the two other terms is less straightforward because they do not represent the same processes as ours but overall the sum of our VM, diffusion and advection terms act as the sum of their “numerical” and vertical mixing i.e. in opposition to the forcing.

The effect of a given tendency term over the entire volume of the ocean is given by its value at the deepest reference depth i.e. -5500 m. At this depth, the diffusion and vertical mixing are both zero, while the forcing is positive and the effect of advection is negative. The volume integral

of the advection is negative because of the imperfect way the free surface boundary condition is formulated in the model as explained in Kuhlbrodt et al. (2015). As explained in Sect. 2, the advection made by the monthly mean velocity on the monthly mean temperature is zero when volume integrated on V(z_r) and is therefore not part of the advection term in Eq. (13). The volume integral of the forcing on the entire volume of the ocean is positive because of the small control run drift.

All of the four terms have a large slope change at very shallow reference depth, around -55 m. It is explained by the fact that low and mid latitudes have shallow reference depths because their surface temperature is mostly contained between approximately 10 and 30 °C whereas the deeper reference depths are confined to high latitudes regions (Fig. 4).

The heating thus only occurs for reference depths shallower than -55 m, while the cooling occurs on a much larger range of reference depths: between -5500 and -55 m.

The negative sign of the volume integrated tendency due to diffusive processes (see Fig. 3) is consistent with the downgradient nature of heat diffusion. Indeed, writing the diffusion term as the divergence of a downgradient heat flux $F_{diff}$ as $-\nabla \cdot F_{diff}$, with $F_{diff} \cdot \nabla \theta < 0$, shows that:

$$ \int_{V(z_r)} \text{diffusion} dV = - \int_{z_r=\text{const.}} \mathbf{F}_{\text{diff}} \cdot \mathbf{n} \, dS = \int_{z_r=\text{const.}} \frac{\mathbf{F}_{\text{diff}} \cdot \nabla \theta}{|\nabla \theta|} \, dS < 0 \quad (14) $$

where $F_{diff}$ is the diffusive flux.

Finally the sum of the advective and diffusive terms almost completely balance the forcing term because the VM is small compared to the three other terms. This is in contrast with the horizontally-averaged heat balance, for which the mean diffusive flux may occasionally be upward and balanced by a mean downward advection (see for instance Kuhlbrodt et al. (2015). Here the main balance is between a downward (toward deeper z_r) diffusion (and advection) and

<figure><img src="image_6.png" /><figcaption>Fig. 4 Time average of the temperature (in °C) at the surface in the CR. -2000, -400, -55 and -10 m contours of the reference depth are shown with dashed black contours</figcaption></figure>
an upward (toward shallower $z_r$) forcing flux where forcing flux can be defined as follows:

$$ \frac{\partial \text{forcing flux}}{\partial z_r} = \frac{\partial}{\partial z_r} \left( \int_{V(z_r)} \text{forcing dV} \right) \quad (15) $$

Figure 5 shows $\frac{\partial}{\partial z_r} (\int_{V(z_r)} \text{term})$ with “term” replaced by either forcing, advection, diffusion or VM. Positive values act to increase the temperature while negative values decreases the temperature.

To facilitate the interpretation of these noisy terms we also show their integration over three ranges of reference depth: [-5500 m, -5000 m], [-5000 m, -55 m] and [-55 m, 0 m]. In [-55 m, 0 m], advection, diffusion and VM all act to decrease the temperature and are balanced by the forcing. In [-5000 m, -55 m] the forcing decreases the temperature and is almost entirely balanced by the diffusion. The sum of all terms in this range of $z_r$ is slightly positive because of the CR drift. In [-5500 m, -5000 m] the forcing is negative and balances the sum of the remaining terms. The magnitude of the diffusion and forcing values in the shallowest and deepest ranges are respectively about two times and one third of that found in the intermediate range although both correspond to a much smaller volume (50 m and 500 m of reference depth vs almost 5000 m). This emphasizes the importance of these two ranges of reference depth for the ocean heat budget.

### 3.2.3 Advective term

In this section we show that the non-zero advection appearing in the above budget (Eq. 12) is approximately balanced by sub-monthly diffusion. To understand what

term balance this sub-monthly advection term, we have run the control run of HiGEM on a year with daily means outputs and repeated the calculation that led to Fig. 3. The comparison between results from monthly means and daily means outputs for the same year is on Fig. 6.

We first notice that the differences between the two time resolution for the forcing and the VM are very small. Secondly, as expected, the advection term is closer to zero when daily means are used rather than monthly means. Recall that it cannot be zero because of the problem in HiGEM with the free surface boundary condition. The diffusive flux is larger with daily means outputs than with monthly means so that the sum of the diffusion and of the advection remains approximately constant between the two outputs frequency (Fig. 6). To understand this, we first write Eq. (2) using monthly mean and anomalies:

$$ \frac{\partial \bar{\theta}}{\partial t} + \frac{\partial \theta'}{\partial t} + \mathbf{v} \cdot \nabla \theta + Ae = \bar{D} + D' + \bar{VM} + VM' + \bar{F} + F', \quad (16) $$

where $Ae$ represents the effect of the imperfect formulation of the free surface boundary condition in HiGEM (Kuhlbrodt et al. 2015), and $D$ the diffusion. We integrate it over the volume, $V(z_r)$, calculated from the daily outputs and average the result:

<figure><img src="image_7.png" /><figcaption>Fig. 5 Left panel: divergence of the volume integrated tendency terms as a function of reference depth ($z_r$) (derivative with respect to $z_r$ of terms shown on Fig. 3). Right panel: vertical integral of the left panel quantities over three different ranges of $z_r$. The three ranges</figcaption></figure>

are: [-5500 m, -5000 m] (third row), [-5000 m, -55 m] (second row) and [-55 m, 0 m] (first row). Unit is $10^{-3}$ K year$^{-1}$. The corresponding time mean temperature (over the CR) is shown for each reference depth range on the right
<figure><img src="image_8.png" /><figcaption>**Fig. 7** Left panel: 10-years time average of the effective diffusivity centered at year 5 (orange) and 55 (blue) of the 2×CO₂ run. Right panel: comparison of the 10 years averaged $K^{\text{diff}}$ (green) and $K^{\text{adv}}$ (red) centered at year 55 of the 2×CO₂ run. $K^{\text{diff}}$ is associated with</figcaption></figure>

monthly mean diffusion and $K_{\text{eff}}^{\text{adv}}$ with sub-monthly diffusion because sub-monthly advection is approximately balanced by sub-monthly diffusion, as explained in Sect. 3.2.3

## 3.3 2×CO₂ run

### 3.3.1 Time evolution of the reference temperature

The left panel of Fig. 8 shows the time evolution of the isotherms’ reference depth in the 70 years of the 2×CO₂ run.

All isotherms are seen to progressively deepen with time as the ocean is getting warmer. The net warming at the end of the 70 year period, defined as the difference between the temperatures at the end and beginning of the period, is depicted in the right panel. This shows that the largest increase (≈ 2.5 °C) occurs at shallow reference depth i.e. at high temperature. The temperatures between reference depths of –4000 m and –2000 m remains almost constant while the temperature between between –5500 and –4000 m i.e. the coldest waters, increases significantly (≈ 0.4 °C).

### 3.3.2 Effective diffusivity

The time mean between years 50 and 60 of the 2×CO₂ run of the effective diffusivities associated with diffusion

(Eq. 20) and advection (Eq. 21) are shown on the right panel of Fig. 7. Despite the differences between the volume integral of the temperature tendencies due to advection and diffusion (see Fig. 3), $K_{\text{eff}}^{\text{diff}}$ and $K_{\text{eff}}^{\text{adv}}$ have similar variation because their reference depth dependence is mainly controlled by the variation of reference temperature gradient $(\frac{\partial \theta_r}{\partial z_r})^{-1}$ (not shown). $K_{\text{eff}}^{\text{diff}}$ and $K_{\text{eff}}^{\text{adv}}$ have very similar magnitude except between –500 m and 0 m, where $K_{\text{eff}}^{\text{adv}}$ is almost one order of magnitude smaller than $K_{\text{eff}}^{\text{diff}}$. The total effective diffusivity during the 2×CO₂ run remains nearly constant in the upper 1000 m, increases in the range 3500–4750 m, and decreases everywhere else.

### 3.3.3 Volume average of the tendencies and heat flux convergence

Figure 9 shows the difference between the 70 years time mean of the volume average of the temperature tendencies in the 2×CO₂ and in the CR as a function of reference depth. The temperature increase found at all reference depths, as

<figure><img src="image_9.png" /><figcaption>**Fig. 8** Left: contours of $\theta_r(z_r, t)$ as a function of time (years) and reference depth ($z_r$) in the 2×CO₂ run. Right: temperature difference (in °C) between the end and the beginning of the 2×CO₂ run as a function of $z_r$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 9** Volume average of the difference between the 2×CO₂ and Control Run of the temperature tendencies due to VM (blue), forcing (orange), diffusion (green) and advection (red) as a function of reference depth ($z_r$). The time mean of the temperature in the reference space over the 2×CO₂ run is shown on the right</figcaption></figure>

shown on Fig. (8), can mainly be attributed to the increase in forcing found in 2×CO₂. The volume integral of the forcing in 2×CO₂ is indeed much larger than that of the CR, with a difference close to $4 \times 10^{-3}$ K year⁻¹.

The heat flux convergences are studied below to understand the time evolution of the temperature at each reference depth. Following Eq. (10) the heat flux convergence at reference depth $z_r$ are obtained by calculating the $z_r$ derivative of the volume integral of the tendencies. Then, we subtract the CR heat flux convergence

from the 2×CO₂ heat flux convergence to understand what term drives the increase in temperature as shown on Fig. 8. As for Fig. 5, the convergence terms are noisy and we thus integrate the results on five different ranges of reference level to facilitate the interpretation. This ranges are: [-5500 m, -5000 m], [-5000 m, -4000 m], [-4000 m, -2000 m], [-2000 m, -55 m] and [-55 m, 0 m] and are chosen to represent the vertical variation of the convergence terms. Right panel of Fig. 10 shows that, as expected, the sum of the four processes (advection, diffusion, VM and forcing) is always positive.

The difference between the forcing of the 2×CO₂ and the CR is positive for all ranges of reference depth while the diffusion is negative everywhere except in [-2000 m, -55 m]. The heat flux from the atmosphere to the ocean thus increases at shallow reference depths and low latitudes (see Fig. 4) whereas the ocean loss of heat to the atmosphere that occurs at deeper reference depth and at mid and high latitudes is reduced. The diffusion intensity increases for low reference depths in the 2×CO₂ run: a larger amount of heat is diffused toward low temperatures than in the CR resulting in a cooling of the ocean for $z_r > -55$ m and in a warming for $z_r$ between -2000 m and -55 m. As shown on Fig. 7, $K_{eff}$ remains approximately constant in the 2×CO₂ run while the gradient of theta ($\frac{\partial\theta}{\partial z_r}$) increases at shallow reference depth (see Fig. 8). The increase in diffusive flux toward low reference depth is thus explained by the increase in the temperature gradient at low reference depth. The largest value of the sum of all terms is found at shallow reference depth in the [-2000, -55] range ($25.77 \times 10^{-3}$ K year⁻¹) where diffusion

<figure><img src="image_11.png" /><figcaption>**Fig. 10** Left panel: difference between the 2×CO₂ and Control runs of the temperature tendencies terms divergence (positive values mean that temperature in the 2×CO₂ increases compared to CR) as a function of reference depth ($z_r$). VM is in blue, forcing in orange, diffusion in green, advection in red and the sum of all terms in purple. Right panel: vertical integral of the left panel terms over five ranges</figcaption></figure>

of reference level. The five ranges are: [-5500 m, -5000 m] (fifth row), [-5000 m, -4000 m] (fourth row), [-4000 m, -2000 m] (third row), [-2000 m, -55 m] (second row) and [-55 m, 0 m] (first row). Unit is $10^{-3}$ K year⁻¹. The corresponding time mean temperature (over the 2×CO₂) is shown for each reference depth range on the right

<figure><img src="image_12.png" /></figure>
(7.12 × 10⁻³ K year⁻¹) and forcing (17.34 × 10⁻³ K year⁻¹) act together to increase the temperature. This range represents 72% of the difference between 2 × CO₂ and CR total heating rate. At the surface the [-2000 m, -55 m] range is approximately located between 60°S and 30°S in the Southern hemisphere, and between 30 and 60°N in the Northern hemisphere (see Fig. 4). In the deepest range ([-5500 m, -5000 m]) the warming effect of the forcing (due to a weaker heat transfer to the atmosphere) is almost entirely balanced by the reduced VM found in the 2× CO₂ compared to the CR.

To sum up, the increase of temperature in the 2×CO₂ is mainly attributed to the increase of heat flux from the atmosphere to the ocean at low reference depth and to the decrease of the heat flux from the ocean to the atmosphere at deeper reference depth, particularly between [-2000 m, -55 m]. Diffusion acts to decrease the temperature at shallow reference depth ([-55 m, -0 m]) and to increase the temperature in the range below ([-2000 m, -55 m]). This is explained by the intensification of the diffusive flux in the upper reference depths of the 2×CO₂, associated with the increased temperature gradient, that results in a higher transfer of heat from high to low temperatures.

# 4 Conclusions

Following up on Huber et al. (2015), this paper explores an alternative way to develop a physical calibration of the classical VAD for the purposes of representing the ocean heat balance and ocean heat uptake in SCMs. VADs based on an Eulerian horizontal average—which represent the majority of existing SCM VAD—are not well suited to the development of physical calibrations, because one of the key terms controlling their time evolution involves the correlation between $w'$ and $\theta'$, defined as departures from a horizontal mean (see Appendix A). Physically, we know that $\theta'$ and $w'$ must be controlled both by surface buoyancy fluxes, wind forcing and interior mixing processes, meaning that we should expect their correlation to be partly advective, partly diffusive. How to perform such a separation in practice is not understood. However, interpreting such a term as purely diffusive reveals that it generally tends to act anti-diffusively, which in Huber et al. (2015) was found to be responsible for occasionally making the effective diffusivity negative, thus explaining the behaviour seen in studies such as Kuhlbrodt et al. (2015). This effect is in theory suppressed when the average is performed along constant $\theta$ surfaces. Indeed, by definition of this average, deviation from isotherms are zero (i.e. $\theta' = 0$) and thus cannot influence the evolution equation. The temperature time evolution is then only due to diathermal diffusion (toward low temperature), to surface heat

fluxes and to parameterized convection/mixed layer dynamics, while the temperature advection plays no role.

Using this new framework, we studied the heat balance and heat uptake in two HiGEM runs, one where the CO₂ concentration is set to 345 ppmv reflecting conditions in the 1980s (the control run), and one where the CO₂ concentration is doubled. In the CR, the balance is mainly between the downward (i.e. toward colder temperatures) sum of advection and diffusion and the upward forcing. Heat flux convergences (i.e. the reference depth derivative of the total heat fluxes) show that above a reference depth of approximately -55 m, the diffusion, advection, VM cool the ocean while the forcing heats the ocean and compensate almost completely this cooling. Below this reference depth (i.e. for most of the ocean volume), the main equilibrium is between the sum of advection and diffusion that heats the ocean and the forcing that cools the ocean. We showed that the advection term is in theory zero in this framework but that in practice it is true only when the outputs frequency is large enough (smaller than a month here). However, we showed that the advection term that appears when the monthly means are used can conveniently be linked to the higher frequency diffusion and that the total diffusion that would have been obtained with high frequency outputs is very close to the sum of the diffusion and advection obtained from monthly average. Further work need to be done to understand if this result can be generalized to other models or to the ocean. In HiGEM, sub-monthly forcing is negligible compared to sub-monthly advection and sub-monthly diffusion, however the question whether this is true in models with a more realistic representation of sub-monthly forcing remains to be addressed.

The effective diffusivity coefficient ($K_{eff}$) of the diathermal diffusion has then be calculated using the sum of tendencies from advection and diffusion. $K_{eff}$ is around $1 \times 10^{-3} \text{ m s}^{-2}$ between -5500 and -2000 m and increases from $1 \times 10^{-6} \text{ m s}^{-2}$ to $1 \times 10^{-6} \text{ m s}^{-2}$ between 0 and -2000 m. The fact that these values are at least one order of magnitude larger than the prescribed vertical diffusion in HiGEM indicates that isopycnal mixing for reference depth below -2000 m plays an important role in the heat budget. In the 2×CO₂ run temperature increases at every reference depth, particularly at shallow reference depth. This temperature increase is attributed to the increase in forcing at all reference depths: the heat flux from the atmosphere to the ocean increases (low reference depths, high temperatures) while the heat flux from the ocean to the atmosphere (mid and deep reference depths, low temperatures) decreases. The diffusive flux increases for reference depths between -2000 and 0 m which results in a cooling above -55 m and a warming between -2000 and -55 m. It contrasts with the results obtained with the horizontal average as in Kuhlbrodt et al. (2015) where the
warming of the horizontally averaged temperature is attributed to the vertical mixing in the top 1000 m and from increased downwelling below -1000 m. Vertical mixing plays no significant role in heat uptake using θ-coordinates except for the deepest reference depths (between -5500 and -5000 m) where it almost balances the warming due to the reduction of heat transfer to the atmosphere. Similarly, downwelling (i.e. advection) is in theory zero as explained above and can in practise be linked to the diffusion so that it does not play any significant role in the θ-averaged model.

Even if the evolution equation for the θ-averaged model is simpler than for the horizontal average model, some work remains to be done before the θ averaged model can be used to calibrate SCM. The two main point that we identify are the time evolution of the effective diffusivity coefficient $K_{eff}$ and of the reference level at the surface (i.e. $z_r(x, y, z = 0, t)$) under a warming climate. Indeed, once the time evolution of the reference level $z_r$ at the surface is known, the surface temperature can be deduced from the $\theta(z_r, t)$. The heat exchanges between the atmosphere and the ocean could then be deduced from the knowledge of this temperature. Understanding this two points would help to predict the evolution of the diffusivity and of the forcing in different $\theta$ classes.

## Appendix A: Revisiting the horizontally-averaged interpretation of the VAD in the light of the θ-based framework

In this appendix, we propose an alternative construction of the horizontally-averaged temperature previously considered by Huber et al. (2015) aimed at making it more easily comparable to the θ-based framework considered in this paper.

To obtain an equation for the horizontal average of the temperature, $\theta$ is first decomposed into its horizontal average plus departure from it:

$$ \theta = \langle \theta \rangle (z, t) + \theta'(x, y, z, t). \quad (22) $$

where $\langle \cdot \rangle$ is the horizontal average and $(x, y, z, t)$ are respectively the longitude, latitude, depth and time coordinates. The temperature departure from the horizontal average $\theta'$ is obtained using $\theta' = \theta - \langle \theta \rangle$ and satisfies: $\langle \theta' \rangle = 0$. Substituting $\theta$ in Eqs. (22) and volume integrating between the ocean surface and some depth $z$ yields:

$$ \begin{aligned} & \int_{V(z)} \frac{\partial \theta'}{\partial t} dV + \frac{\partial}{\partial t} \int_{V(z)} \langle \theta \rangle dV \\ & \qquad + \langle \theta \rangle \int_{x,y} w \, dxdy + \int_{x,y} w \theta' \, dxdy \\ &= \int_{x,y} \mathbf{K} \nabla \langle \theta \rangle \cdot \mathbf{k} \, dxdy + \int_{x,y} \mathbf{K} \nabla \theta' \cdot \mathbf{k} \, dxdy \\ & \qquad + \int_{V(z)} \mathrm{VM} \, dV + \int_{V(z)} \mathrm{Q}_{\mathrm{net}} \, dV \end{aligned} \quad (23) $$

where $V(z)$ is the volume between the surface and depth $z$ and $\int_{x,y} dxdy$ the horizontal integral on the whole ocean surface. The non-divergence of the velocity field has been used to transform the volume integral of the advection into a surface integral at depth $z$. The term involving the surface heating is zero in a statistical steady-state but positive in global warming experiments.

By definition of $\theta'$, the first term on the left hand side is zero and the third term of the rhs is also zero because of volume conservation, Eq. (23) then becomes:

$$ \begin{aligned} & \frac{\partial}{\partial t} \int_{V(z)} \langle \theta \rangle dV + \int_{x,y} w \theta' dxdy \\ &= \int_{x,y} \mathbf{K} \nabla \langle \theta \rangle \cdot \mathbf{k} \, dxdy + \int_{x,y} \mathbf{K} \nabla \theta' \cdot \mathbf{k} \, dxdy \\ & \qquad + \int_{V(z)} \mathrm{VM} \, dV + \int_{V(z)} \mathrm{Q}_{\mathrm{net}} \, dV. \end{aligned} \quad (24) $$

To make the calculation concrete, we assume that $\mathbf{K} = K_i(\mathbf{I} - \mathbf{d}\mathbf{d}^T) + K_d\mathbf{d}\mathbf{d}^T$ is a rotated diffusion tensor (Redi 1982), with $\mathbf{I}$ the identity tensor, $\mathbf{d}$ the unit normal vector pointing in the dianeutral direction, $K_i$ and $K_d$ the turbulent isoneutral and dianeutral mixing coefficients respectively. As a result, the projection of the diffusive flux of the horizontally average $\theta$ i.e.: $\mathbf{K} \nabla \langle \theta \rangle \cdot \mathbf{k}$ may be written as:

$$ \begin{aligned} \mathbf{K} \nabla \langle \theta \rangle \cdot \mathbf{k} &= \frac{\partial \langle \theta \rangle}{\partial z} \mathbf{K} \mathbf{k} \cdot \mathbf{k} = \frac{\partial \langle \theta \rangle}{\partial z} (K_i(\mathbf{k} - (\mathbf{k} \cdot \mathbf{d})\mathbf{d}) \\ & \qquad + K_d(\mathbf{k} \cdot \mathbf{d})\mathbf{d}) \cdot \mathbf{k} \\ &= \frac{\partial \langle \theta \rangle}{\partial z} \underbrace{(K_i \sin^2(\mathbf{k}, \mathbf{d}) + K_d \cos^2(\mathbf{k}, \mathbf{d}))}_{=K_{\mathrm{eff}}^{\mathrm{loc}}} \end{aligned} \quad (25) $$

$K_{\mathrm{eff}}^{\mathrm{loc}}$ is the local effective mixing coefficient and is always a positive quantity. Let us now take the derivative of (24) with respect to $z$ and divide the result by $A(z)$, the depth-dependent ocean area at depth $z$, which yields the following equation for $\langle \theta \rangle$:
$$ 
\begin{aligned}
& \frac{\partial \langle \theta \rangle}{\partial t} + \frac{1}{A(z)} \frac{\partial}{\partial z} \int_{x,y} w \theta' dxdy \\
& = \frac{1}{A(z)} \frac{\partial}{\partial z} \left( K_{\text{eff}} A(z) \frac{\partial \langle \theta \rangle}{\partial z} \right) + \frac{1}{A(z)} \frac{\partial}{\partial z} \left( \int_{x,y} \mathbf{K} \nabla \theta' \cdot \mathbf{k} \, dxdy \right) \\
& \quad + \frac{1}{A(z)} \frac{\partial}{\partial z} \left( \int_{V(z)} \text{VM} \, dV \right) + \frac{1}{A(z)} \frac{\partial}{\partial z} \left( \int_{V(z)} \text{Q}_{\text{net}} \, dV \right)
\end{aligned}
\quad (26) $$

where

$$ K_{\text{eff}} = \frac{1}{A(z)} \int_{x,y} K_{\text{eff}}^{\text{loc}} dxdy \quad (27) $$

is the horizontally-averaged $K_{\text{eff}}^{\text{loc}}$, and where we have used:

$$ \frac{\partial}{\partial z} \frac{\partial}{\partial t} \int_{V(z)} \langle \theta \rangle \, dV = \frac{\partial}{\partial z} \int_{z}^{0} A(z^v) \frac{\partial \langle \theta \rangle}{\partial t} (z^v) dz^v = -A(z) \frac{\partial \langle \theta \rangle}{\partial t} \quad (28) $$

Note that the derivation of (26) could have been done using directly an horizontal average on (2) instead of a volume integral followed by a vertical derivation. This method is preferred here to emphasize the similarity with the average along $\theta$ surfaces as will become clear in the following section.

In Huber et al. (2015), the VAD equation is written in the form:

$$ \frac{\partial \langle \theta \rangle}{\partial t} + w^* \frac{\partial \langle \theta \rangle}{\partial z} = \frac{\partial}{\partial z} \left( k^* \frac{\partial \langle \theta \rangle}{\partial z} \right) + Q, \quad (29) $$

where $k^*$ is an effective vertical diffusive coefficient, $w^*$ an effective vertical velocity and $Q$ a source term. Comparing this equation with Eq. (26) and identifying like for like terms suggests the following associations:

$$ 
\begin{aligned}
& \frac{1}{A(z)} \frac{\partial}{\partial z} \left( K_{\text{eff}} A(z) \frac{\partial \langle \theta \rangle}{\partial z} \right) \\
& \quad + \frac{1}{A(z)} \frac{\partial}{\partial z} \left( \int_{x,y} \mathbf{K} \nabla \theta' \cdot \mathbf{k} \, dxdy \right) \\
& \quad \Leftrightarrow \quad \frac{\partial}{\partial z} \left( k^* \frac{\partial \langle \theta \rangle}{\partial z} \right)
\end{aligned}
\quad (30) $$

$$ \frac{1}{A(z)} \left( \frac{\partial}{\partial z} \left( \int_{V(z)} \text{VM} \, dV \right) + \frac{\partial}{\partial z} \left( \int_{x,y} w \theta' dxdy \right) \right) \Leftrightarrow w^* \frac{\partial \langle \theta \rangle}{\partial z} \quad (31) $$

$$ \frac{1}{A(z)} \frac{\partial}{\partial z} \left( \int_{V(z)} \text{Q}_{\text{net}} dV \right) \Leftrightarrow Q \quad (32) $$

The diffusive terms in Eq. (26) are identified with the diffusive part of Eq. (29), the advective and vertical mixing terms with the advective part of (29) and the forcing with the forcing term of (29). This identification can be used to obtain a physical calibration of the VAD equation $w^*$ and $k^*$

as in Huber et al. (2015). As noted by Huber et al. (2015), this choice is not unique and the VM term could also be attributed to the diffusive part of the VAD for instance. The main novelty here is that the divergence of the diffusive flux is explicitly split into a part involving the vertical gradient of $\langle \theta \rangle$ which is always downgradient and a part involving the departure of $\theta$ from the horizontal average which can a priori be negative or positive. The fact that the horizontal mean diffusive heat flux can occasionally transport heat upwards, as first showed by Gregory (2000), means that the latter term may occasionally counteract the effect of the former term.

Due to volume conservation, the horizontally-averaged vertical velocity must vanish at all depths (i.e. $\langle w \rangle = 0$), and therefore cannot contribute to the effective advection $w^*$ of $\langle \theta \rangle$. As shown by Eq. (31) $w^*$ is rather associated with the advection of $\theta'$ through horizontal surfaces and to VM. Negative values around $w^* \approx -0.5 \times 10^{-7} \text{ms}^{-1}$ are found in Huber et al. (2015) for the resolved and eddy parametrized advection: advection of $\theta'$ transports heat downward.

To sum up, averaging the heat budget along depth levels introduces two new non-negligible terms involving horizontal temperature anomalies $\theta'$, which complicates the description of the horizontal mean temperature $\langle \theta \rangle$. This motivates us to seek an approach that avoid the introduction of such anomalies.

From Eq. (26), the only way to remove the terms involving $\theta'$ in (31) and (30) is to perform the average along temperature surfaces instead of horizontally. Indeed, $\theta'$ is then zero by definition.

## Appendix B: General framework accounting for freshwater fluxes

In this appendix, we consider the more general case where the effect of the free surface and of freshwater fluxes are not neglected in Eq. (5). Note that the full derivation of this equation is given in more details in Hochet and Tailleux (2019), for the sake of conciseness, we only give the main result here.

At the top, the ocean is bounded by a free surface of equation $z = \eta(x, y, t)$. With the effect of freshwater fluxes and of the free surface included, Eq. (5) of Sect. 2 becomes:

$$ 
\begin{aligned}
& \int_{V(z_r)} \frac{\partial \theta_r}{\partial t} dV + \int_{S(z_r)} (\theta_s - \theta_r)(E - P + R) dS \\
& \quad + \int_{V(z_r)} \nabla \cdot \overline{\mathbf{v}' \theta'} dV \\
& = \int_{z_r=\text{const.}} \mathbf{K} \nabla \theta \cdot \mathbf{n} dS + \int_{V(z_r)} \text{VM} \, dV + \int_{V(z_r)} \text{Q}_{\text{net}} \, dV,
\end{aligned}
\quad (33) $$
where P, E and R represents respectively precipitation, evaporation and river runoff, $\theta_s$ is the temperature at the surface, $S(z_r)$ is the outcropping surface corresponding to the surface $z_r = \text{const.}$ and where we have use the fact that $\nabla \cdot \mathbf{v} = 0$ imposes at each time:

$$ \int_{z_r=\text{const.}} \mathbf{v} \cdot \mathbf{n} \, dS = \int_{S(z_r)} (E - P + R) \, dS. \quad (34) $$

## Appendix C: Volume integral calculation

In this appendix we describe how the volume integral of the different terms in the temperature tendencies budget (i.e. Eq. 13) is calculated.

For each monthly mean output, each grid cell is vertically divided into ten smaller volumes at the center of which $\theta$ is linearly interpolated. The value of the term we want to integrate is divided by 10 and attributed to each of the ten sub-volumes corresponding to each grid cell. This procedure allows one to have a better resolution and to conserve the volume integral of the term. We experimentally found that using no subdivision leads to a larger amount of noise when calculating the $z_r$ derivative of the integrated term and that a larger number of vertical subdivision (> 10) has no significant effect on the results. Then, for each time step $t$, the minimum $\theta_{min}$ and maximum $\theta_{max}$ of $\theta$ are obtained. An array $\theta_{vec} = [\theta_{min}, \theta_{min} + \Delta\theta, \theta_{min} + 2\Delta\theta, ..., \theta_{max}]$ is then constructed with $\Delta\theta = \frac{\theta_{max} - \theta_{min}}{N}$ and $N = 1000$. For each temperature $\theta_{vec}^i$ (with $i \in 1, ..., 1000$) in this array we sum the volume times the integrated term of all parcels with $\theta$ satisfying $\theta > \theta_{vec}^i$:

$$ \sum_{j=1}^{N_j} \text{term}_j \times \Delta V_j \quad (35) $$

where $\Delta V_j$ is the parcel's volume and $N_j$ the number of parcels with $\theta$ satisfying $\theta > \theta_{vec}^i$. Using continuous notation for clarity we now have:

$$ \int_{V(\theta, t)} \text{term} \, dV \quad (36) $$

And we make use of the previously calculated $z_r(\theta, t)$ to obtain:

$$ \int_{V(z_r)} \text{term} \, dV \quad (37) $$

**Acknowledgements** This work was supported by the OUTCROP NERC grant NE/R010536/1. A. Hochet is supported by EU Marie Curie IF grant number 749924. We thank the two anonymous reviewers who helped improve a first version of this manuscript.

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

## References

Brierley CM, Collins M, Thorpe AJ (2010) The impact of perturbations to ocean-model parameters on climate and climate change in a coupled model. Clim Dyn. 34(2-3):325-343

Dukowicz JK, Smith RD (1994) Implicit free-surface method for the bryan-cox-semtner ocean model. J Geophys Res Oceans 99(C4):7991-8014

Exarchou E, Kuhlbrodt T, Gregory JM, Smith RS (2015) Ocean heat uptake processes: a model intercomparison. J Clim 28(2):887-908

Ferrari R, Ferreira D (2011) What processes drive the ocean heat transport? Ocean Model 38(3):171-186

Gent PR, Mcwilliams JC (1990) Isopycnal mixing in ocean circulation models. J Phys Oceanogr 20(1):150-155

Gregory JM (2000) Vertical heat transports in the ocean and their effect on time-dependent climate change. Clim Dyn 16(7):501-515

Gregory JM, Mitchell JF (1997) The climate response to co2 of the hadley centre coupled aogcm with and without flux adjustment. Geophys Res Lett 24(15):1943-1946

Griffies SM (1998) The Gent-McWilliams skew flux. J Phys Oceanogr 28:11

Hieronymus M, Nilsson J, Nycander J (2014) Water mass transformation in salinity temperature space. J Phys Oceanogr 44(9):2547-2568. https://doi.org/10.1175/jpo-d-13-0257.1

Hochet A, Tailleux R (2019) Comments on “diathermal heat transport in a global ocean model”. J Phys Oceanogr 49(8):2189-2193

Hochet A, Tailleux R, Ferreira D, Kuhlbrodt T (2019) Isoneutral control of effective diapycnal mixing in numerical ocean models with neutral rotated diffusion tensors. Ocean Sci 15(1):21-32

Holmes RM, Zika JD, England MH (2018) Diathermal heat transport in a global ocean model. J Phys Oceanogr 49(1):141-161. https://doi.org/10.1175/JPO-D-18-0098.1

Huang B, Stone PH, Sokolov AP, Kamenkovich IV (2003) The deep-ocean heat uptake in transient climate change. J Clim 16(9):1352-1363

Huber M, Tailleux R, Ferreira D, Kuhlbrodt T, Gregory J (2015) A traceable physical calibration of the vertical advection-diffusion equation for modeling ocean heat uptake. Geophys Res Lett 42(7):2333-2341. https://doi.org/10.1002/2015GL063383

Kuhlbrodt T, Gregory J (2012) Ocean heat uptake and its consequences for the magnitude of sea level rise and climate change. Geophys Res Lett 39(18):18608

Kuhlbrodt T, Gregory J, Shaffrey L (2015) A process-based analysis of ocean heat uptake in an aogcm with an eddy-permitting ocean component. Clim Dyn 45(11-12):3205-3226

Levitus S, Antonov JI, Boyer TP, Baranova OK, Garcia HE, Locarnini RA, Mishonov AV, Reagan J, Seidov D, Yarosh ES et al (2012) World ocean heat content and thermosteric sea level change (0-2000 m), 1955-2010. Geophys Res Lett 39(10):1-5
Marshall DP, Zanna L (2014) A conceptual model of ocean heat uptake under climate change. J Clim 27(22):8444–8465. https://doi.org/ 10.1175/jcli-d-13-00344.1

Meinshausen M, Meinshausen N, Hare W, Raper SC, Frieler K, Knutti
R, Frame DJ, Allen MR (2009) Greenhouse-gas emission targets
for limiting global warming to 2 c. Nature 458(7242):1158–1162

Meinshausen M, Raper S, Wigley T (2011) Emulating coupled atmos-
phere-ocean and carbon cycle models with a simpler model,
magicc6-part 1: Model description and calibration. Atmos Chem
Phys 11(4):1417–1456

Raper S, Gregory JM, Osborn T (2001) Use of an upwelling-diffusion energy balance climate model to simulate and diagnose a/ogcm results. Clim Dyn 17(8):601–613

Raper SC, Gregory JM, Stouffer RJ (2002) The role of climate sen-
sitivity and ocean heat uptake on aogcm transient temperature
response. J Clim 15(1):124–130

Redi MH (1982) Oceanic isopycnal mixing by coordinate rotation. J Phys Oceanogr 12(10):1154–1158

Shaffrey LC, Stevens I, Norton WA, Roberts MJ, Vidale PL, Harle JD, Jrrar A, Stevens DP, Woodage MJ, Demory ME, Donners J, Clark DB, Clayton A, Cole JW, Wilson SS, Connolley WM, Davies

TM, Iwi AM, Johns TC, King JC, New AL, Slingo JM, Slingo
A, Steenman-Clark L, Martin GM (2009) UK HiGEM: the new
UK high-resolution global environment model: model description
and basic evaluation. J Clim 22(8):1861–1896. https://doi.org/10.
1175/2008JCLI2508.1

Winters KB, D'Asaro EA (1996) Diascalar flux and the rate of fluid mixing. J Fluid Mech 317:179-193

Winters KB, Lombard PN, Riley JJ, D'Asaro EA (1995) Available potential energy and mixing in density-stratified fluids. J Fluid Mech 289:115–128. https://doi.org/10.1017/S002211209500125X

Wolfe CL, Cessi P, McClean JL, Maltrud ME (2008) Vertical heat transport in eddying ocean models. Geophys Res Lett 35(23):23605. https://doi.org/10.1029/2008GL036138

Wyrtki K (1961) The thermohaline circulation in relation to the general circulation in the oceans. Deep Sea Res (1953) 8(1):39-64

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