A traceable physical calibration of the vertical advection-diffusion equation for modelling ocean heat uptake

The classic vertical advection-diﬀusion (VAD) balance is a central concept in studying the ocean heat budget, in particular in simple climate models (SCMs). Here we present a new framework to calibrate the parameters of the VAD equation to the vertical ocean heat balance of two fully-coupled climate models that is traceable to the models’ circulation as well as to vertical mixing and diﬀusion processes. Based on temperature diagnostics, we derive an eﬀective vertical velocity w ∗ and turbulent diﬀusivity k ∗ 𝜈 for each individual physical process. In steady state, we ﬁnd that the residual vertical velocity and diﬀusivity change sign in middepth, highlighting the diﬀerent regional contributions of isopycnal and diapycnal diﬀusion in balancing the models’ residual advection and vertical mixing. We quantify the impacts of the time evolution of the eﬀective quantities under a transient 1% CO 2 simulation and make the link to the parameters of currently employed SCMs


Introduction
The vertical ocean heat balance plays a fundamental part in the Earth's energy balance, where the latter describes the partitioning of additional radiative forcing in the climate system into ocean heat uptake and a temperature response [e.g., Gregory and Forster, 2008;Knutti and Hegerl, 2008;Church et al., 2011;Huber and Knutti, 2011].For example, the change in the vertical ocean heat balance under transient climate change is crucial in setting the magnitude and temporal evolution of the global temperature response: the smaller the heat uptake by the ocean, the larger the global temperature response, and vice versa.
A traditional framework for discussing the vertical ocean heat balance and ocean heat uptake is based on the following vertical advection-diffusion equation with temperature Θ, upwelling velocity w, and diffusivity k  : (1) Equation ( 1) is perhaps most famous for serving as the basis for the abyssal recipes papers by Munk [1966] and Munk and Wunsch [1998], where the latter study also considered a depth-dependent k  .Physically, such a model regards the vertical heat balance as a competition between the cooling effect due to the upwelling of cold abyssal waters induced by high-latitude deep water formation on the one hand and the downward diffusion of heat on the other hand.Applying this model to a bounded region of the Pacific Ocean at middepth resulted in typical values of the two parameters of w ∼ 0.7 × 10 −7 m/s and k  ∼ 1.0 × 10 −4 m 2 ∕s [Munk and Wunsch, 1998].
The above upwelling-diffusion balance is often used as a basis for simple climate models (SCMs), for example, the MAGICC model [e.g., Meinshausen et al., 2011].Key goals of such SCMs are to emulate the ocean heat uptake under climate change of more complex climate models and to compute probabilistic estimates of future temperature changes under a range of emission scenarios that provide crucial information for policy makers [Raper et al., 2001;Meinshausen et al., 2009Meinshausen et al., , 2011;;Rogelj et al., 2012].Accordingly, projections of future surface-air and ocean temperatures change strongly depend on the two parameters w and k  of the upwelling-diffusion equation, highlighting the importance and need of an accurate physical and statistical calibration of the parameters.
Currently, most SCMs employ a behavioral calibration approach, which refers to the specification of the parameters such that the SCM reproduces the behavior of the more complex model that it aims to emulate.
HUBER ET AL.
Typically, the parameters w and k  , assumed to be constant in time and uniform with depth, are fitted to an idealized forcing simulation of the complex model, e.g., an abrupt step CO 2 forcing.In this method, it is difficult to assess how the estimated values of w and k  relate to the physics of the more complex model, for example, to the model's residual circulation or its isopycnal and diapycnal diffusion.In addition, a key assumption when using the SCMs is that the calibrated parameters can be employed for various other forcing scenarios.Considering the complexity of the ocean dynamics, one wonders whether the calibration should be scenario dependent.
Over the past decade, a number of modeling studies have seemingly challenged the classical upwelling-diffusion view of the processes affecting the oceanic heat balance [Gregory, 2000;Wolfe et al., 2008;Hieronymus and Nycander, 2013].These studies find that in a quasi steady state, advection by the mean flow in both eddy-parameterizing and eddy-resolving climate models is associated with a downward heat flux with eddy fluxes tending to be in the opposite direction and that isopycnal diffusion is a major term in balancing the advective heating.Note that these model results describe globally averaged values at a particular depth level, while the classical upwelling-diffusion balance by Munk [1966] initially focused on the Pacific Ocean in middepth.
Based on the temperature diagnostics of an eddy-parameterizing and an eddy-permitting climate model, we present a new framework for achieving a physical calibration of w and k  by linking the two parameters to the actual physical processes of the models, resulting in an effective vertical velocity w * and an effective turbulent diffusivity k *  for each advective, diffusive, and mixing process.Such temperature diagnostics have been analyzed previously in the context of idealized transient and abrupt climate change scenarios [e.g., Gregory, 2000;Exarchou et al., 2014].Here we focus on the implications of the models' vertical heat balance on the magnitude and the sign of the vertical velocity w and diffusivity k  of equation ( 1) that form the physical basis of many SCMs.We additionally discuss the role of eddy advection and isopycnal diffusion which are missing in SCMs.We further present the time evolution of the effective quantities under an idealized transient climate change simulation and demonstrate that these spatial and time variations are key to evaluating the transient ocean heat uptake.

Climate Model Data
We use annual mean output of a 70 year control simulation and a 70 year transient climate change simulation with an annual 1% increase in CO 2 from two coupled climate models: HadCM3 and HiGEM1.2.HadCM3 is a version of the Hadley Centre coupled climate model and includes an atmospheric model with horizontal resolution of 2.5 • × 3.75 • and 19 vertical levels [Gordon et al., 2000].The ocean model is a rigid lid, primitive equation general circulation model with a horizontal resolution equal to 1.25 • × 1.25 • with 20 vertical levels.The isopycnal scheme of Griffies et al. [1998] is used for isopycnal diffusion with a diffusion coefficient equal to 1000 m 2 s −1 .The eddy-induced tracer transport is parameterized following Gent et al. [1995] with an eddy-induced diffusion coefficient with values between 300 and 2000 m 2 s −1 that depends on the local stratification.The vertical mixing of tracers is based on the Richardson number-dependent formulation by Pacanowski and Philander [1981] with a background diffusivity equal to  b = 10 −5 m 2 s −1 that increases with depth.Mixed-layer processes are parameterized by the Kraus-Turner scheme, and convection is parameterized as complete mixing according to Rahmstorf [1993].
HiGEM1.2 uses a higher-resolution compared to HadCM3, with 0.83 • latitude × 1.25 • longitude in the atmosphere and 1/3 • × 1/3 • with 40 levels in the ocean and is considered eddy permitting.Similar to HadCM3, its lateral tracer mixing is based on isopycnal scheme of Griffies et al. [1998] with a constant isopycnal diffusivity of 500 m 2 /s.In contrast to HadCM3, present-day boundary conditions were chosen for the control run with an atmospheric CO 2 concentration of 345 ppmv, reflecting conditions in the 1980s.More information on the two models, in particular with regard to their simulated ocean heat balance, can be found in Exarchou et al. [2014].We use a MATLAB loess filter to smooth the annual mean values in time.

Temperature Tendencies
In general, the change of potential temperature Θ of a model grid cell is given by the sum of its temperature tendencies associated with different physical processes: HUBER ET AL.
denoting the rate of heat change due the divergence of heat fluxes (F) related to advection (ADV), diffusion in the horizontal and vertical direction (DIFF), vertical mixing processes (VM, including mixed-layer processes and convection), and surface processes (SURF, composed of fluxes due to solar penetration, the sea surface and sea ice physics).For HadCM3, the advection of heat is a sum of the resolved mean flow and parameterized eddy advection.The partitioning is similar to Exarchou et al. [2014].The vertical diffusion tendency is decomposed into the individual contributions by isoneutral and diapycnal fluxes using an Interactive Data Language (IDL) routine called Partial Ocean Temperature Tendency Emulator (POTTE), described in Exarchou et al. [2014].Based on the temperature and salinity fields of the model, POTTE derives the isoneutral and diapycnal diffusion by reconstructing potential density surfaces in the ocean.For a horizontally averaged ocean column (denoted by the overbar), which is the focus here, the contribution by horizontal diffusion cancels, viz., (3)

Effective Vertical Velocity and Diffusivity
We start by considering the horizontally averaged advection-diffusion balance in the following form: where Q denotes a source term.Comparing this expression with the temperature evolution in equation ( 3), one can relate the individual process-based temperature tendencies with the advective and diffusive part of the classic upwelling-diffusion balance; that is, z Thus, the effective velocity w * is associated with adiabatic processes that redistribute heat in the ocean such as advection by the resolved flow, whereas the effective diffusivity k *  includes the diabatic, tracer-variance-destroying process of vertical diffusion.We explicitly note that physically w * is not related to the horizontally averaged vertical velocity, as the latter must be zero by continuity.The effective velocity w * should be regarded as a parameter in the simple model representation of the vertical advection of the horizontally averaged temperature profile Θ(z), chosen to reproduce the vertical divergence of the advective temperature fluxes of the complex climate model in a physically traceable manner.The fluxes represented by w * include the eddy fluxes  z (w ′ Θ ′ ), which are not explicitly represented in the SCMs.We consider here vertical mixing as an advective process.However, the method allows us to regard each process either as advective or diffusive.For example, in a model where convection is parameterized with an enhanced diffusivity, the term  z F VM could be represented in equation ( 6) instead of equation ( 5).In addition, choosing vertical mixing as an effective advective quantity facilitates the link to the commonly used partitioning of the advection-diffusion balance in SCMs (see section 3).External forcings resulting from penetrating solar radiation and from sea surface fluxes are represented in the source term Q.
Since the expressions on the right-hand side of the equations are given model outputs both for the control and climate change simulation, equations ( 5) and ( 6) constitute a simple equation for w * and a differential equation for k *  , which we denote henceforth as effective vertical velocity and diffusivity since they are related to the actual circulation and diffusion and mixing processes in the climate models.The values for w * are discretized at the center of each depth level and k *  is defined at the boundaries and can be obtained by integrating equation ( 6) from the bottom of the ocean upward, similar to Wolfe et al. [2008] for the Tropical Ocean.The supporting information provides further details regarding the computation of the effective quantities as well as a discussion regarding the discretization and the choice of boundary conditions.

Results
We start by considering the zonally and horizontally integrated steady state heat tendencies associated with vertical diffusion shown in Figure 1.Using POTTE, we further show the individual isopycnal and diapycnal components.Both models show a warming tendency in the deep ocean dominated by diapycnal diffusion (∼3-5 km), whereas the total vertical diffusion in the upper ocean below the mixed layer exhibits a cooling tendency due to the strong isopycnal diffusive cooling in the high latitudes, especially in the Southern Ocean.The middepth region in the Tropical Ocean as considered by Munk and Wunsch [1998] is highlighted in Figure 1 and features a general warming due to diapycnal diffusion.

The associated effective diffusivities k *
are shown in Figures 1c and 1d.In both models, the total effective vertical diffusivity features strong vertical variations and even a sign change in middepth between ∼1000 and 1500 m.Additionally, k *  decreases by about an order of magnitude from the bottom toward the HUBER ET AL.
©2015.The Authors.surface.The effective diffusivity k * −DIA for diapycnal diffusion almost linearly decreases toward the surface, whereas the mostly negative isopycnal diffusivity k * −ISO is more confined to the upper and middepth oceans.To compare these values with the estimates of Munk and Wunsch [1998] and common values in SCMs, we volume average the individual diffusivities over the depth range between 1 and 4 km which are illustrated in the grey-shaded box at the bottom of Figures 1c and 1d.For the total vertical diffusivity, we derive averaged values of about 0.4 × 10 −4 m 2 s −1 (HadCM3) and 0.9 × 10 −4 m 2 s −1 (HiGEM1.2).Note that these values are based on the horizontally averaged potential temperature and temperature tendencies, which include the high-latitude oceans as well.The corresponding values for k * −DIA are 0.67 × 10 −4 m 2 s −1 (HadCM3) and 1.2 × 10 −4 m 2 s −1 (HiGEM1.2).In both models, the derived isopycnal value of k * −ISO is about −0.27 × 10 −4 m 2 s −1 .HUBER ET AL.
©2015.The Authors.For a steady state to prevail, the temperature tendencies associated with advection by the resolved and eddy-parameterized flow as well as with vertical mixing need to balance the warming and cooling tendencies of vertical diffusion.Figure 2 shows the corresponding zonally integrated warming and cooling rates for advection.For the eddy-parameterizing model HadCM3, values for the resolved Eulerian advection and eddy-parameterized advection are shown individually.Generally, both models show an advective warming below the mixed layer down to depths of around 3500 m and a cooling below.Similar to previous findings by Gregory [2000], Hieronymus and Nycander [2013], and Exarchou et al. [2014], the resolved advection in the eddy-parameterizing model exhibits a warming tendency over most of the water column which is offset by a general cooling associated with the parameterized eddy fluxes, which results in a change of sign at a depth of approximately 3500 m (Figure 2a).In terms of the classic advection-diffusion balance, the strong advective warming in the high-latitude ocean dominates the cooling in tropical latitudes, thus resulting in an overall warming in the top 3000 m.
The corresponding effective vertical velocities w * of the two models are shown in Figures 2c and 2d.By definition, a positive velocity is associated with upwelling and cooling of a particular depth level, whereas a negative velocity corresponds to downwelling and warming.In terms of total advection, both models show a similar behavior of downwelling in the upper 3000 m and an upwelling in the abyssal ocean with velocities mostly ranging between −1 and 1 × 10 −7 m/s.We also associate the effect of vertical mixing with an effective vertical velocity, resulting in an upwelling (cooling) of most depths levels in the two models.Due to the change of sign in the deep ocean, volume averaging the values for the total advection over 1-4 km results in a rather small downwelling (warming) effective velocity of −5.5 × 10 −9 m/s (HadCM3) and −3.7 × 10 −8 m/s (HiGEM1.2),respectively.Adding the values of vertical mixing results in a net upwelling (cooling) with an effective velocity of 0.13 × 10 −7 m/s for HadCM3 and a net downwelling (warming) of −0.25 × 10 −7 m/s for HiGEM1.2.Overall, Figures 1 and 2 show that each of the different physical processes considered here shows a distinct vertical profile of their effective quantities w * and k *  .The time evolution of the effective quantities under a transient 1% CO 2 simulation is illustrated in Figure 3  physical processes.The evolution of w * and k *  is computed for each year of the transient simulation using the vertical profile of potential temperature and the temperature tendencies of the two models to compute w * (t, z) and k *  (t, z) using equations ( 5) and ( 6).In terms of effective velocity, both models show a decrease of w * associated with advection and vertical mixing, representing an increased downwelling of the flow and a reduced upwelling associated with vertical mixing, which is likely dominated by changing vertical stability of the water column in the high latitudes and the associated changes in convection.Overall, the total effective velocity in middepth (comprising both advection and vertical mixing) is increased by about 76% (HadCM3) and 78% (HiGEM1.2) at the time of double CO 2 (year 70).The change in the upper ocean is dominated in both models by the decrease in the upwelling w * VM from vertical mixing.With regard to the effective diffusivity, Figures 3c and 3d show an increase of the total vertical diffusivity k *  during the transient simulation, which is more pronounced for HadCM3 where it even changes sign from slightly negative to positive values.For both models, the diapycnal diffusivity remains fairly constant, whereas the isopycnal diffusivity shows a reduction over time that dominates the change of the total vertical diffusivity which is likely due to the change in isotherm slopes relative to isopycnals in the high latitudes during the transient CO 2 forcing.For the upper ocean, the diffusivities do not show significant changes during the transient simulation.
To assess the implications of the effective quantities for modeling ocean heat uptake, we built a forward model of the advection-diffusion balance (equation ( 4)) using the effective velocities w * and diffusivities k *  , which we force with the transient temperatures of the top 30 m and the vertical source terms Q(t, z) of the transient 1% CO 2 simulation (see supporting information for further information).Note that the source terms include the increased energy input through the ocean surface due to increased CO 2 forcing.We focus here on HadCM3 as the control drift in HiGEM1.2likely affects the transient values of k *  due to the integration from the bottom to the top.
To establish the link to the currently employed SCMs which are based on a positive diffusivity aimed to represent diapycnal mixing, we exploit the above mentioned mathematical ambiguity in w * and k *  to redefine k *  such that it only represents diapycnal mixing and combine all other physical processes into the advective part w * , thereby avoiding numerical issues with vertical diffusivity that changes HUBER ET AL.

Geophysical Research Letters
10.1002/2015GL063383 sign (refer to Figure 1).The resulting setup is shown in Figures 4a and 4b and shows a similar behavior as the classic picture of downward diffusion and upwelling over the whole water column.However, the effective quantities still feature strong vertical variations.Using the time-dependent values of the transient simulation, the forward model manages to reproduce the temperature change at the time of double CO 2 for the top 3000 m (Figure 4c).The model shows a negative drift in the abyssal ocean which might be related to control drift of HadCM3.
To assess the importance of the time evolution of the effective quantities, we replace the transient values for individual processes, i.e., the effective diapycnal diffusivity, with the control values and run the model forward.Figure 4c shows that the forward model is most sensitive to the time evolution of the effective velocity associated with vertical mixing as the temperature increase in the top 1000 m at double CO 2 is strongly underestimated using the control values for w * VM .Other processes also affect the prediction of temperature, notably in different depth regions.For example, using control values for advection, the heat uptake in middepth between ∼750 and 2000 m is underestimated.The temperature prediction is more sensitive to the change in isopycnal diffusion (expressed via the corresponding change in its effective w * ) than to the change in diapycnal diffusivity, which is fairly constant under the transient simulation (refer to Figure 3).Overall, Figure 4 suggests that not accounting for the time evolution of individual effective quantities can result in errors of about 1 • C for the upper ocean and up to ∼ 0.1 • C for middepth levels at time of CO 2 doubling, corresponding to roughly ∼50% of the forced signal in the middepth ocean.

Discussion and Conclusion
Using temperature tendencies of individual physical processes that set the ocean heat balance of two coupled climate models, we introduce a new framework to relate the vertical heat budget to the classic upwelling-diffusion balance of Munk [1966], which is the basis of many simple climate models.The effective velocities w * and vertical diffusivities k *  constitute a traceable physical calibration of the horizontally averaged advective, diffusive, and mixing processes of the two models.The effective advection-diffusion balance considered here is based on the horizontally averaged profile of potential temperature: various regions with distinct physical mechanisms such as the Tropical Ocean with its classical upwelling-diffusion balance and the Southern Ocean with its strong isopycnal diffusive and eddy advective fluxes are being spatially averaged.We find that the effective quantities w * and k *  change sign in middepth.As such, the framework presented here takes into account that the horizontally averaged ocean heat balance is dominated by the extratropics [e.g., Exarchou et al., 2014], a region which is treated in SCMs mostly implicitly.
We show that the effective quantities w * and k *  change under an idealized transient climate change simulation, most notably the effective velocity related to vertical mixing processes (Figure 3) dominated by changing stability properties and the associated convective heat exchanges in the high latitudes under transient CO 2 forcing.With the use of a simple forward model we demonstrate that the time evolution of each physical process is crucial in simulating the vertical profile of global ocean heat uptake under a transient 1% CO 2 simulation (Figure 4).
Note that a time dependency of the two parameters is also accounted for in the current generation of SCMs, for example, by empirically relating the decrease in upwelling w to the reduction in deep water formation rate under greenhouse gas warming [Raper et al., 2001] or by defining a warming-dependent vertical diffusivity k  to account for the enhanced vertical stratification under climate change [Meinshausen et al., 2011].In that context, Figures 3 and 4 suggest that the time dependency of the advective effective quantities as well as the time evolution of the isopycnal effective diffusivity play a larger role under transient CO 2 -induced ocean heat uptake than the corresponding change of the diapycnal diffusivity over time.
We emphasize that the partitioning of the physical processes into an advective and diffusive part can be chosen differently, for example, one could represent all processes as purely diffusive, as previously in Raper et al. [2001] as an application of their upwelling-diffusion model.As the main motivation of this study is to link the various physical processes of more complex models with the parametric and numerical setup of simple upwelling-diffusion climate models, we argue that the partitioning of the different processes into advective and diffusive effective quantities as outlined in equations ( 5)-( 7) based on their adiabatic or diabatic nature as well as their implementation in the more complex models (i.e., isopycnal diffusion is associated with an effective diffusivity) provides a physically traceable framework to relate the mechanisms of ocean heat uptake across different models of varying complexity.Interestingly, considering isopycnal HUBER ET AL.

Geophysical Research Letters
10.1002/2015GL063383 diffusion as an effective advective quantity as shown in Figure 4, mimicking the setup of SCMs, minimizes the time dependence of the two parameters w * and k *  , as the diapycnal diffusivity is almost constant under a transient CO 2 simulation.
We conclude that the framework presented here provides a way to illustrate the control and transient vertical oceanic heat budget in a climate model using familiar concepts of vertical velocity and vertical diffusivity, thus providing also a means for model intercomparison regarding the magnitude and vertical distribution of processes governing ocean heat uptake under climate change.As the effective process-based quantities show a strong spatial and time dependency, future work is required to understand their steady state and transient behavior for them to be employed for a wider range of emission scenarios.

Figure 1 .
Figure 1.Zonally and horizontally integrated steady state temperature tendencies associated with vertical diffusion in (a) HadCM3 and (b) HiGEM1.2.Units are TW (10 12 J/s).The total vertical diffusion tendency (solid lines) is decomposed into the diapycnal and isopycnal contribution (dashed and dotted lines).(c, d) The corresponding effective diffusivities k *  .The grey-shaded boxes of Figures 1c and 1d show the volume-averaged values of k *  between 1 and 4 km depth that are compared to the values of Munk and Wunsch [1998] (green lines).

Figure 2 .
Figure 2. (a, b) Similar to Figure 1 but for the zonally and horizontally integrated temperature tendencies associated with advection.(c, d) The effective vertical velocities w * .For HadCM3, the heat fluxes of the resolved flow and eddy-parameterized flow are shown separately.

Figure 3 .
Figure 3.Time evolution of the effective vertical velocities w * and diffusivities k *  for the control simulation (black) under a transient 1% CO 2 simulation (colored lines) averaged over the middepth ocean (300-2000 m) for (a, c) HadCM3 and (b, d) HiGEM1.2.Different processes are indicated with different colors.
. Volume-averaged values over the middepth ocean (chosen to be 300-2000 m) are shown for different HUBER ET AL. ©2015.The Authors.

Figure 4 .
Figure 4. Partitioning of the physical processes into an (a) effective vertical velocity w * and an (b) effective diffusivity k *  used to simulate the time evolution of the control and transient ocean heat content of HadCM3 by forward stepping equation (4).(c) The prediction of the temperature change at time of CO 2 doubling at the end of the transient 1% CO 2 simulation is shown.The colored dashed lines denote the cases where the control steady state values of the effective quantities are employed during the transient forward run.