Global heat balance and heat uptake in potential temperature coordinates

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.


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 horizontallyaveraged heat balance is often compared with the theory of early models of the deep circulation such as Wyrtki (1961) 1 3 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 2 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 behavorial 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 ( hereafter) surfaces, following an approach similar to that recently developed by . The averaged diapycnal diffusion is then automatically downgradient and we will further show that the advection through 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 2 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 2 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 2 concentration is doubled. The last section concludes and discusses the results.

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: where X represents any term of the heat budget, (.) the monthly average and (.) � the deviation from this monthly average so that X � = 0 . The time mean potential temperature conservation can be written as: For clarity we will drop the overline notation in what follows and keep it only when it involves anomalies. is the 3D velocity vector, ∇ ⋅ � � is a term representing the effect of sub-monthly advection, a diffusion tensor representing the effect of unresolved advection and small-scale irreversible mixing, ∇ ⋅ with the upward unit vector is thus zero at the surface. 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 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 . 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 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 is the same after and before the adiabatic rearrangement i.e.: where V( , t) is the volume of ocean with temperature l satisfying < l < max with 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 (x, y, z, t) as a function of z r : r (z r , t) = (x, y, z, t) . r can be inverted to yield z r = z r ( , t) or z r = z r (x, y, z, t) . Note that Eq. (3) shows that the volume V( , t) = V(z r ) of water of temperatures greater than 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: where 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. Schematic showing how the reference depth z r associated with temperature 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 r (z r , t) . The volume for all water parcels warmer than (X, t) = 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 r (z, t) by integrating (2) over the volume V(z r ) , which after some manipulation yields: where = −∇ ∕|∇ | = −∇z r ∕|∇z r | is the outward unit normal vector to the isothermal surface = constant , which at fixed time coincides with the surface z r (x, y, z, t) = constant . The diffusion term can be written as: 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: 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: This equation links the volume integral on V(z r ) of the time derivative of the temperature to the diffusive flux, the submonthly 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 r : where we have used: The possibility to obtain a 1D equation for 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 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× CO 2 run of a climate model.

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 (10) atmospheric CO 2 concentration is doubled ( 2× CO 2 ). The control run length used in this article is 50 years and the 2× CO 2 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, = const. surfaces are exactly the same as the z r = const. surfaces. However, the fact that 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 2 run. They are then averaged over a 50 years period for the CR and on the 70 years of the 2 ×CO 2 run.

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 the HiGEM grid and used in the reference depth calculation (see formula 3) is shown on the right panel of Fig. 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: 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: 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 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 straigthforward 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 where 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 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 emphasize the importance of these two ranges of reference depth for the ocean heat budget.

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 (15) forcing flux z r = z r ∫ V(z r ) forcing dV 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: where Ae represents the effect of the imperfect formulation of the free surface boundary condition in HiGEM , and D the diffusion. We integrate it over the volume, V(z r ) , calculated from the daily outputs and average the result:  Figure 6 shows that the time mean of the volume integral of a monthly mean term is very similar when calculated on daily outputs V(z r ) or monthly outputs V(z r ) , and that the volume integral of ′ t , VM ′ and F ′ are negligible, giving all equalities added in Eq. (17). Comparing Eqs. (9) and (17) we deduce that: which shows that the residual advection appearing when volume integrating with monthly means is approximately equal to the higher frequency diffusion of temperature.
To sum up, we showed in this section that part (the other part is associated with Ae) of the volume integrated advection is associated with the sub-monthly diffusion.

Effective diffusivity
In what follows, we return to the analysis of the monthly means. The above results motivates us to include the nonvanishing advection term as part of our definition of effective diffusivity. The effective diffusivities associated with sub-monthly diffusion via the advection and associated to the monthly mean diffusion are calculated using the two following formulas: with: where the double overline denotes here the time mean over the 50 years of the CR. K eff is shown on the left panel of Fig. 7 and is seen to increase with depth from values around 1 × 10 −6 m 2 s −1 for z r = 0 to 2 × 10 −3 m 2 s −1 at −3500 m . It then decreases to 5 × 10 −4 m 2 s −1 at approximately −4500 m and increases again to 2 × 10 −3 m 2 s −1 for the deepest z r . Note that these values of the diathermal diffusive coefficient are at least one order of magnitude larger than the values 0(10 −5 m 2 s −1 ) commonly observed in the thermocline. This is mainly because the temperature gradient is generally not parallel to the neutral direction so that part of the large isoneutral mixing occurs in the diathermal direction. Warm waters associated with very shallow reference depth ( > −55 m ) have an effective diffusivity smaller than 10 −5 m 2 s −1 down to 10 −6 m 2 s −1 . This is partly explained by the large temperature gradient (i.e. r z r ) found at these reference depths as can be seen on Fig. 2. In the following section, we study the heat balance under a warming climate using a HiGEM run where the atmospheric concentration of CO 2 is doubled. Fig. 6 Comparison of the monthly outputs (plain line) average vs daily outputs (dash line) average of the volume integral of the heat tendencies (in K year −1 ) respectively as a fonction of z r (calculated on monthly means ) and z r (calculated on daily mean of ). VM is in blue, forcing in orange, diffusion in green, advection in red, sum of all terms in purple and the sum of diffusion and advection is in brown. Also shown (dotted, indistinguishable from plain) are the monthly means of the heat tendencies integrated on the volumes defined by z r instead of z r , but both are very close and thus indistinguishable on the figure. Note that the difference with Fig. 3 for the monthly outputs is due to the time mean performed on only one year compared to the 50 years of Fig. 3. The time mean of the temperature as a function of z r is shown on the right

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 2 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).

Effective diffusivity
The time mean between years 50 and 60 of the 2 ×CO 2 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 diff eff and K adv eff have similar variation because their reference depth dependence is mainly controlled by the variation of reference temperature gradient ( r z r ) −1 (not shown). K diff eff and K adv eff have very similar magnitude except between −500 m and 0 m, where K adv eff is almost one order of magnitude smaller than K diff eff . The total effective diffusivity during the 2 ×CO 2 run remains nearly constant in the upper 1000 m, increases in the range 3500-4750 m, and decreases everywhere else. Figure 9 shows the difference between the 70 years time mean of the volume average of the temperature tendencies in the 2 ×CO 2 and in the CR as a function of reference depth. The temperature increase found at all reference depths, as and reference depth ( z r ) in the 2 ×CO 2 run. Right: temperature difference (in °C) between the end and the beginning of the 2 ×CO 2 run as a function of z r shown on Fig. (8), can mainly be attributed to the increase in forcing found in 2 ×CO 2 . The volume integral of the forcing in 2 ×CO 2 is indeed much larger than that of the CR, with a difference close to 4 × 10 −3 K year −1 .

Volume average of the tendencies and heat flux convergence
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 substract the CR heat flux convergence from the 2 ×CO 2 heat flux convergence to understand what term drives the increase in temperature as shown on  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 2 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 latiutdes is reduced. The diffusion intensity increases for low reference depths in the 2 ×CO 2 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 2 run while the gradient of theta ( r 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 × 10 −3 K year −1 ) where diffusion Fig. 9 Volume average of the difference between the 2 ×CO 2 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 2 run is shown on the right Unit is 10 −3 K year −1 . The corresponding time mean temperature (over the 2 ×CO 2 ) is shown for each reference depth range on the right ( 7.12 × 10 −3 K year −1 ) and forcing ( 17.34 × 10 −3 K year −1 ) act together to increase the temperature. This range represents 72% of the difference between 2 × CO 2 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 2 compared to the CR.
To sum up, the increase of temperature in the 2× CO 2 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 2 , associated with the increased temperature gradient, that results in a higher transfer of heat from high to low temperatures.

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 ′ , defined as departures from a horizontal mean (see Appendix A). Physically, we know that ′ 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 surfaces. Indeed, by definition of this average, deviation from isotherms are zero (i.e. � = 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 2 concentration is set to 345 ppmv reflecting conditions in the 1980s (the control run), and one where the CO 2 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 × 10 −3 m s −2 between −5500 and −2000 m and increases from 1 × 10 −6 m s −2 to 1 × 10 −6 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 2 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 (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 classes.
where P, E and R represents respectively precipitation, evaporation and river runoff, s is the temperature at the surface, S(z r ) is the outcropping surface corresponding to the surface z r = const. and where we have use the fact that ∇ ⋅ = 0 imposes at each time:

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 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 min and maximum max of are obtained. An array vec = min , min + Δ , min + 2Δ , … , max is then constructed with Δ = max − min N and N = 1000 . For each temperature i vec (with i ∈ 1, … , 1000 ) in this array we sum the volume times the integrated term of all parcels with satisfying > i vec : where ΔV j is the parcel's volume and N i the number of parcels with satisfying > i vec . Using continuous notation for clarity we now have: And we make use of the previously calculated z r ( , t) to obtain: