The rate at which the ocean moves heat from the tropics toward the poles, and from the surface into the interior, depends on diabatic surface forcing and diffusive mixing. These diabatic processes can be isolated by analyzing heat transport in a temperature coordinate (the diathermal heat transport). This framework is applied to a global ocean sea ice model at two horizontal resolutions (1/4° and 1/10°) to evaluate the partioning of the diathermal heat transport between different mixing processes and their spatial and seasonal structure. The diathermal heat transport peaks around 22°C at 1.6 PW, similar to the peak meridional heat transport. Diffusive mixing transfers this heat from waters above 22°C, where surface forcing warms the tropical ocean, to temperatures below 22°C where midlatitude waters are cooled. In the control 1/4° simulation, half of the parameterized vertical mixing is achieved by background diffusion, to which sensitivity is explored. The remainder is associated with parameterizations for surface boundary layer, shear instability, and tidal mixing. Nearly half of the seasonal cycle in the peak vertical mixing heat flux is associated with shear instability in the tropical Pacific cold tongue, highlighting this region’s global importance. The framework presented also allows for quantification of numerical mixing associated with the model’s advection scheme. Numerical mixing has a substantial seasonal cycle and increases to compensate for reduced explicit vertical mixing. Finally, applied to Argo observations the diathermal framework reveals a heat content seasonal cycle consistent with the simulations. These results highlight the utility of the diathermal framework for understanding the role of diabatic processes in ocean circulation and climate.
The ocean plays a critical role in the climate system by transferring heat from the tropics to higher latitudes (Trenberth and Caron 2001) and from the surface to depth (Gregory 2000). How this heat transport may change in response to future changes in surface forcing, both due to anthropogenic influences (Trenberth et al. 2014) and due to natural variability acting on interannual (Roemmich and Gilson 2011) and decadal (England et al. 2014) time scales, remains a first-order question in climate science. Critical to understanding and predicting these changes is an understanding of the role of the various diabatic processes responsible for transferring heat into and around the ocean and an accurate representation of these processes in climate models (e.g., Bryan 1987; Ferrari and Ferriera 2011).
Heat enters the surface ocean predominantly in the tropics, particularly in the eastern equatorial upwelling regions at sea surface temperatures (SSTs) between 22° and 29°C (Large and Yeager 2009; Valdivieso et al. 2017; Fig. 1a). In contrast, surface cooling predominantly occurs poleward of 25° latitude over the poleward flowing western boundary currents at SSTs below ~22°C (Fig. 1a). Therefore, in order to connect these regions of surface heat gain and loss, the ocean must transport heat not only from low to high latitudes (the meridional ocean heat transport) but also across isotherms from warm water toward colder water (the diathermal ocean heat transport). Indeed, in the model simulations discussed here, the tropical regions where the annual-mean SST exceeds 21.5°C receive a net positive heat flux of 1.6 PW, which must thus be transported across the annual-mean 21.5°C isotherm to be lost at cooler temperatures. This diabatic heat transfer, which acts to oppose and balance the surface forcing, ultimately occurs via small-scale diffusive mixing processes within the ocean (Niiler and Stevenson 1982; Döös et al. 2012; Zika et al. 2012; Hieronymus et al. 2014).
Two factors could alter conclusions about across-isotherm fluxes made from maps of the annual-mean surface heat flux and surface temperature: 1) the penetration of shortwave radiation into the ocean that warms colder subsurface layers (Iudicone et al. 2008) and 2) the net contribution of correlations between SST and surface heat flux variability at seasonal and subseasonal time scales (e.g., see Figs. 1b,c,d). However, when taken into account in the model simulations discussed here these factors do not substantially alter the requirement stated above that diabatic mixing processes must transport ~1.6 PW of heat across the 21.5°C isotherm. Without such diabatic mixing processes, the SST distribution, surface heat flux, interior water mass structure, and heat transport (Speer and Tziperman 1992) of the ocean would all be different because of the requirement that, to maintain a steady state, the net surface warming of each temperature class must be zero.
The realization that significant amounts of heat must be transferred across isotherms motivates an analysis of the ocean heat transport in a temperature coordinate. Such an approach, pioneered by Walin (1982), isolates the diabatic processes associated with surface forcing and mixing that alter the temperature of individual fluid parcels. This approach contrasts with, for example, studies that focus on the vertical heat transport budget (e.g., Gregory 2000; Wolfe et al. 2008; Exarchou et al. 2015; Kuhlbrodt et al. 2015), which, due to the averaging coordinate used, is complicated by contributions from many processes, including adiabatic advection. Zika et al. (2013) and Zika et al. (2015a) did employ a water mass framework to help understand the processes setting vertical heat transport, but precise budgets were not achieved.
By applying the temperature coordinate framework to observations of surface temperature, density, and heat and freshwater fluxes, several studies have highlighted the necessity of mixing within the warm, low-density tropical waters to close the oceanic heat budget (Speer and Tziperman 1992; Zhang and Talley 1998; Toole et al. 2004; Groeskamp et al. 2017). Modeling studies support these observational results, showing that mixing by parameterized small-scale three-dimensional turbulence (as opposed to parameterized mixing arising from mesoscale along-isopycnal stirring, which is thought to be more important at higher densities) provides the diathermal heat transport out of the pool of warm tropical waters necessary to balance the net warming of those waters by surface fluxes (McWilliams et al. 1996; Nurser et al. 1999; Hieronymus et al. 2014). However, these studies used low-resolution models and did not examine in detail the role of different parameterized mixing fluxes. Also, most attention has focused on the annually averaged, steady-state heat transport, with less attention paid to seasonal variability.
In this study, we revisit the water mass framework to analyze the annual-mean and seasonally varying diathermal heat transport within a modern global ocean sea ice model (introduced in section 2) at two different fine resolutions (1/4° and 1/10°) and with precise tendency diagnostics. We introduce the concept of internal heat content, which obeys a simpler budget than the full heat content and is independent of the arbitrary reference temperature typically used to define heat content (see section 3). The annual mean diathermal heat transport budget (section 4) provides an intuitive and simple measure of the contribution of different mixing processes. The main contributing processes include explicit vertical mixing associated with parameterized background diffusion, boundary layer mixing, shear-instability mixing, and bottom-intensified internal tide mixing, as well as numerical mixing associated with truncation errors in the model’s advection scheme. Numerical mixing, which we precisely quantify without the need for targeted experiments, makes a significant contribution to the diathermal heat budget and increases to compensate for reduced (or zero) explicit background diffusion (section 5). Examination of the seasonal cycle of the diathermal heat transport (section 6) highlights the important role played by the eastern equatorial Pacific cold tongue, where strong, seasonally varying shear instability makes a globally significant contribution to transferring heat toward cooler waters. Our results, discussed and summarized in section 7, highlight the utility of the diathermal heat transport framework for understanding the role of various diabatic processes in ocean circulation and their impact on climate.
2. Model and data
a. Global ocean–sea ice model
We use two global ocean–sea ice configurations based on the ocean and sea ice components of the GFDL CM2.5 and CM2.6 climate models respectively (Griffies et al. 2015; Delworth et al. 2012). The two configurations, hereafter referred to as MOM025 and MOM01, differ only in their horizontal and vertical resolution. Specifically, MOM025 has a 1/4° Mercator horizontal resolution and 50 vertical levels while MOM01 has a 1/10° Mercator horizontal resolution and 75 vertical levels [optimally positioned following Stewart et al. (2017)]. Both models are forced by a repeating seasonally varying atmospheric state based on version 2 of the Coordinated Ocean-Ice Reference Experiment Normal Year Forcing (CORE-NYF; Large and Yeager 2004). Shortwave radiation is distributed into the ocean interior using the scheme of Manizza et al. (2005), modulated by the climatological chlorophyll distribution compiled by Sweeney et al. (2005). There is no representation of geothermal heating.
Both models are considered eddy-resolving in the tropics, and neither includes a mesoscale eddy parameterization or explicit isopycnal or lateral diffusion. The suppression of large lateral tracer gradients is thus achieved by the multidimensional piecewise parabolic tracer advection scheme [Colella and Woodward (1984), with a monotonicity-preserving flux limiter following Suresh and Huynh (1997)]. We will refer to diffusion associated with the advection scheme as numerical mixing (discussed in more detail in section 3e). Vertical diffusion of both tracers and momentum is parameterized using the K-profile parameterization (KPP; Large et al. 1994), which governs mixing within the surface boundary layer, as well as interior convection, Richardson number–based shear instability, and double diffusion. There is also a bottom-enhanced internal tide mixing scheme (Simmons et al. 2004). The contributions of the different components of the vertical mixing are discussed in more detail in section 4c. Horizontal friction is implemented with a biharmonic operator and a Smagorinsky scaling for the viscosity coefficient (Griffies and Hallberg 2000). Momentum advection is achieved via a second-order centered operator with third-order Adams–Bashforth time stepping. The two model configurations are used to examine the sensitivity of our results to resolution (section 5). More information on the MOM025 and MOM01 configurations can be found in Spence et al. (2017) and Stewart et al. (2017).
b. Sensitivity to background diffusivity
We will also assess the sensitivity of our results to the choice of background vertical diffusivity included in MOM025 (section 5). We consider four cases:
No explicit background diffusivity . This case has been spun up for 513 years, with analysis coming from the last 5 years.
A background diffusivity of m2 s−1. This case has been spun up for 81 years initialized from the case, with analysis coming from the last year.
A background diffusivity of m2 s−1. This case has been spun up for 31 years initialized from the case, with analysis coming from the last year.
A latitude-dependent background diffusivity set to m2 s−1 poleward of 15° latitude, reducing linearly to m2 s−1 by 5° latitude. This case has been spun up for 80 years with analysis coming from the last 5 years.
A value of 10−5 m2 s−1 is typically used in many climate models, being representative of the average levels of internal wave–driven mixing in the ocean interior (e.g., Waterhouse et al. 2014). The m2 s−1 case can be considered a weak mixing case, and as will be discussed below, shows similar results to the case. The latitudinally dependent case follows the implementation used in the GO5.0 configuration of NEMO (Megann et al. 2014) and represents the observed reduction of internal wave–driven background mixing close to the equator (Gregg et al. 2003; also see Zhu and Zhang 2018). We will focus on the latitudinally dependent case, referred to as MOM025 Control, as we consider it the most representative of the observations. In all cases a background vertical viscosity of 10−4 m2 s−1 is used. For MOM01, we consider only a simulation with zero background diffusivity , which has been spun up for 72 years, with analysis coming from the last year.
c. Argo observations
We supplement the model results with observations of water mass structure down to 2000 m from the global network of Argo floats (Argo 2000). The Argo observations span the period 2004–14 and can be used to estimate the volume within different temperature classes and how this changes with time, as discussed in the next section.
3. Diathermal framework
Here we describe the diathermal framework, following Walin (1982), that we will use to analyze the model’s heat budget.
a. Eulerian heat budget
The heat content (relative to 0°C) at each location within the ocean is governed by the Eulerian heat budget:
where is the three-dimensional velocity, is potential temperature, is a constant reference density, and is a constant so-called specific heat capacity.1 We assume a Boussinesq fluid, as in the simulations. The tendency T represents the change in heat content within a unit volume with time. This heat content changes due to the convergence of advective heat fluxes A, vertical mixing processes M including vertical diffusion (with vertical diffusivity κ) and any nonlocal heat fluxes associated with processes such as surface boundary layer convection, surface forcing F, and horizontal diffusion H (explicitly zero in our simulations). In the appendix we provide more details on the numerical implementation of this heat budget, including the numerical values for and and how the heat budget and its diathermal counterpart are diagnosed.
b. Diathermal volume budget
It is convenient to analyze the heat budget not in geographical coordinates, but in a temperature coordinate, as this allows us to isolate the diabatic processes associated with mixing and surface forcing. We consider the volume2 of all water warmer than a given temperature (Walin 1982; see Fig. 2),
where the integral is performed following the temporally and spatially varying temperature field of the global ocean . Note that we will use script characters to denote variables that have been integrated over temperature layers as in Eq. (2), and therefore are functions only of temperature and time. There are two processes that can change the volume (blue arrows in Fig. 2),
where is the volume flux across the isotherm, or the water-mass transformation, and is the surface volume flux,
where PER is the net volume flux per unit area (m s−1) into the ocean associated with precipitation, evaporation, river runoff, and ice melt.
c. Diathermal heat budget
The heat content contained within the volume can be expressed using a Cartesian integral or an integral in temperature space,
where represents the volume within a given temperature interval , with the negative sign present because is defined in terms of the fluid warmer, not colder, than . The budget for the heat content is
where the processes leading to changes in (illustrated by the red arrows in Fig. 2) are as follows:
Surface heat fluxes [corresponding to the integral over the volume of its Eulerian counterpart F in Eq. (1)], associated with radiative, latent, and sensible surface heat fluxes and any ice–ocean heat fluxes associated with, for example, the formation of frazil ice.
- The heat input associated with surface volume fluxes . In the model simulations all surface volume fluxes are considered to enter and exit the ocean at the SST, and so do not alter the temperature of the surface layer, only its volume. Therefore can be related to the surface volume flux [Eq. (4)] via
The heat flux across the isotherm due to explicit vertical mixing [corresponding to the integral over of M in Eq. (1)], associated with the parameterized vertical diffusion and convection. We will also diagnose the contributions to from different processes separately.
Numerical mixing , which is associated with the advection scheme in the model and will be calculated by residual (see section 3e).
The heat flux associated with across-isotherm volume fluxes or water-mass transformation .
Note again that in our model simulations there is no explicit horizontal or along-isopycnal diffusion [captured in H in Eq. (1)], and so we do not include such a term (this role is taken by the numerical mixing term).
d. Internal and external heat content
The heat content can be split into two components by integrating Eq. (6) by parts,
The first component , associated with the volume of the layer and the boundary temperature, will be referred to as the external component of the heat content. The second component , which will be referred to as the internal heat content, can be expressed in terms of how much the volume average temperature of the layer exceeds [equivalent to the area-integrated relative heat content of Palmer and Haines (2009)]:
In this way, the internal heat content is independent of the zero reference chosen for the temperature scale. A budget for the tendency of internal heat content,
The budget for the internal heat content, Eq. (12), is simpler than the full heat content budget [Eq. (7)] and no longer explicitly contains the heat flux associated with water-mass transformation . The heat flux associated with surface volume fluxes has also been replaced by its equivalent internal component . Only those surface volume fluxes that enter at temperatures greater than increase the average temperature of the layer and contribute to the internal heat content , quantified by .
We will focus on the internal heat content and its budget [Eq. (12)] and will only briefly discuss the external heat content and the total heat content tendency . We focus on Eq. (12) because it allows us to isolate the direct influence of the mixing and surface-driven diffusive heat fluxes from the water-mass transformation that they induce. The heat fluxes control the water-mass transformation through a temperature derivative of Eq. (12):
Therefore, the total heat content tendency budget [Eq. (7)] effectively contains two contributions from each diabatic process: one contribution from their heat flux and one contribution from the across-isotherm volume flux (and its associated heat transfer) induced by that heat flux. Because of the temperature derivatives in Eq. (14), this second contribution is more susceptible to small variations in the heat fluxes with temperature and time, and thus can suffer from noise (see the appendix). In contrast, the heat fluxes themselves are more robust to analyze in the model simulations. In any case, at steady state Eqs. (7) and (12) are equivalent, as the global flux of volume (and heat) across isotherms associated with must be zero in order to conserve volume within each temperature class (apart from the small influence of the steady-state volume fluxes ). Therefore, the net transfer of heat across isotherms required to balance the surface forcing can only be achieved through the diffusive heat fluxes and , as most cleanly expressed in the internal heat content budget [Eq. (12)].
The terms on the rhs of Eq. (12) are obtained by binning the Eulerian heat budget diagnostic terms [Eq. (1), and its numerical counterpart Eq. (A1) in the appendix] within 0.5°C temperature bins. This binning process is carried out “online” (i.e., at every model time step) and monthly averages of each term are then saved. The heat and volume content tendencies are calculated by applying the same online temperature binning process to snapshots of the temperature field at the beginning and ending of each month. Likewise, we also estimate the internal heat content tendency from Argo observations, by calculating the change in global volume within each temperature class () and using Eq. (11). The numerical details, including various sensitivity tests performed to check robustness to different methods and numerical choices, are discussed in the appendix.
e. Numerical mixing
Our model simulations do not include any explicit parameterization for lateral mixing, or mixing along density surfaces, and thus rely on the numerical advection scheme to smooth sharp lateral tracer gradients. Moreover, the background vertical diffusivity is specifically chosen to be small (or zero) in several of our configurations.3 Such “numerical mixing,” which is present to varying degrees even in simulations with explicit lateral mixing and larger vertical diffusivities, has traditionally been difficult to characterize and quantify in ocean models (e.g., see Griffies et al. 2000; Ilicak et al. 2012; Gibson et al. 2017). The diathermal framework discussed here provides a useful quantification of numerical mixing by isolating its effects on the heat and volume budgets of temperature layers. By quantifying all the other terms in the internal heat content budget [Eq. (12)], numerical mixing can be calculated by residual:
Since the advection operator included in the simulations is three-dimensional, the numerical mixing diagnosed here may include the artificial suppression of both lateral and vertical temperature gradients, processes that cannot be evaluated separately. Similarly, we cannot evaluate what proportion of the numerical mixing is associated with spurious diapycnal mixing versus that portion that represents the small-scale dissipation of along-isopycnal temperature gradients. This framework provides information on the amount of numerical mixing occurring globally, as a function of both time and temperature, without the need to turn off other sources of diabatic changes in temperature, such as surface forcing (as is commonly done in other studies of numerical mixing; e.g., Griffies et al. 2000; Gibson et al. 2017).
4. Annually averaged heat fluxes
a. Global diathermal heat budget
The annually averaged global diathermal heat budget shows a balance between surface forcing (black line in Fig. 3) and mixing (red and blue lines in Fig. 3). Surface forcing warms the warmest waters, adding 1.6 PW of heat to waters warmer than 21.5°C. Heat is removed by surface forcing at colder temperatures, with 0.7 PW lost in the polar regions from waters cooler than 2°C, and some lost at intermediate temperatures (0.4 PW between 13° and 21°C). Thus, surface forcing acts to broaden the temperature distribution of the ocean, making the warmest waters warmer and the coldest waters colder (Nurser et al. 1999). This broadening is opposed by mixing, which acts to compress the temperature distribution by homogenizing water masses (Hieronymus et al. 2014; Zika et al. 2015b). In this model, the downgradient heat fluxes that achieve this homogenization are explicit vertical mixing (red line in Fig. 3) and numerical mixing (blue line in Fig. 3). The vertical mixing flux peaks, like the surface forcing, at warmer temperatures driving 0.9 PW of heat toward colder temperatures at 21.5°C. The diathermal heat transfer achieved by this mixing is critical for the ocean’s meridional heat transport as it allows the heat that enters the ocean in the tropical regions at warmer temperatures (>21.5°C) to be transferred to the cooler temperatures at which heat is lost back to the atmosphere in the midlatitudes.
Numerical mixing accounts for the remaining downgradient heat flux required to balance the surface forcing, contributing 0.7 PW at 21.5°C, and being relatively uniform across temperatures (blue line in Fig. 3). The contribution of numerical mixing exceeds that of explicitly parameterized vertical mixing at cooler temperatures between 14° and 2°C. At these cooler temperatures, numerical mixing may be achieving the role of along-isopycnal eddy mixing, which in lower-resolution models is parameterized using an explicit Redi diffusion (e.g., see Fig. 4a of Hieronymus et al. 2014). However, as mentioned in the previous section it is not possible to separate the numerical mixing into along-isopycnal and diapycnal components, and so a significant or even dominant diapycnal contribution cannot be discounted. The redistribution of shortwave radiation from the surface into the interior also plays a significant role (also see Hieronymus and Nycander 2013), accounting for a 0.5-PW flux of heat across the 21.5°C isotherm (gray dashed lines in Fig. 3, included in ).
The tendencies of both the total heat content and its internal component (magenta dashed and solid lines respectively in Fig. 3) are relatively weak. As discussed above, the total heat content tendency (and its external component ) tends to vary much more from year to year and from temperature to temperature than the internal heat content tendency. While weak relative to the other terms in the budget, there is a net heat flux into the ocean of 0.2 PW, indicating that the oceanic heat content has not yet equilibrated. This heat flux drives warming of the deep ocean within the temperature range from 1° to 8°C (thick magenta line in Fig. 3). Most of this tendency in heat content is due to the short 80-yr equilibration to the change in background diffusivity. This will be discussed in more detail in section 5.
b. Spatial structure of the heat fluxes
The spatial structure of the heat fluxes, which underlie the annually averaged diathermal heat budget, is heterogeneous. The annually averaged net surface heat flux is largest in the eastern equatorial Pacific, where it can exceed 150 W m−2 (Fig. 1a). This region is a hot spot for heat input into the ocean due to the substantial solar heating associated with cool SSTs, reduced evaporation, and clear skies (Large and Yeager 2009). There is a total flux of 0.63 PW into the cold tongue region of the eastern equatorial Pacific (defined by the Niño-3 region, 5°S–5°N, 150°–90°W). This is a large fraction of the net surface heat flux of 1.93 PW into the equatorial region (10°S–10°N), equivalent to the meridional heat flux out of the same region. It is also large compared with the global maximum diathermal heat flux of 1.6 PW (Fig. 3). The eastern equatorial Pacific region therefore plays a globally significant role in determining the pathways of heat transport through the ocean.
Mixing is responsible for moving heat from one temperature class to another within the ocean interior, and thus the spatial structure of mixing also plays an important role in determining heat transport pathways (Ferrari and Ferreira 2011; Döös et al. 2012; Zika et al. 2012). The vertical mixing heat flux across the 21.5°C isotherm (chosen as it hosts the maximum diathermal flux) also shows a strongly heterogeneous structure (Fig. 4a), reflecting the spatial structures of the vertical temperature gradient and the different parameterized processes that contribute to vertical mixing. In particular, the strongest vertical mixing fluxes occur in the eastern equatorial Pacific cold tongue due to the strong near-surface thermocline and shear-driven mixing in that region, as discussed next.
c. Components of mixing
A number of vertical mixing processes contribute to the diathermal heat flux (Fig. 6). In MOM025 Control, which has a background diffusivity of 10−5 m2 s−1 poleward of ±15° latitude linearly scaled down to 10−6 m2 s−1 within 5° of the equator, background mixing makes the dominant contribution to the flux over temperatures between 7° and 26°C, peaking at 20°C (green solid line in Fig. 6). This temperature corresponds to the subtropical thermocline, where vertical temperature gradients are large over a broad area of the global ocean. At 21.5°C, the heat flux due to background mixing is distributed over a wide band of the subtropics between 10° and 40° latitude, but is strongest in the subtropical North Pacific and Atlantic (Fig. 4b), where vertical temperature gradients below the mixed layer are largest.
Shear instability–driven mixing, parameterized with a Richardson number–dependent diffusivity (Large et al. 1994), is focused at temperatures between 18° and 27°C (purple line in Fig. 6). Shear-driven mixing occurs almost exclusively in the eastern equatorial Pacific in the strongly sheared and stratified region above the Equatorial Undercurrent (Figs. 4c and 5, which shows the latitude–depth structure at 110°W) where turbulence levels are high and the gradient Richardson number remains low throughout much of the year (Gregg et al. 1985; Smyth and Moum 2013; Moum et al. 2013). This component of vertical mixing is responsible for the largest fluxes (with an annual average exceeding 150 W m−2) and makes a significant contribution to the diathermal heat budget, highlighting the global importance of the eastern equatorial Pacific cold tongue.
Mixing in the surface boundary layer makes a significant contribution across the full temperature range associated with mixing near the surface isotherm outcrops (blue line in Fig. 6). Boundary layer mixing is peaked near 26°C, associated with broad mixing across the tropical regions where entrainment into the surface layer from the highly stratified waters just below the mixed layer base is common. At 21.5°C, boundary mixing is distributed over the seasonally shifting outcrop bands in the midlatitudes, and in the eastern Pacific cold tongue (Fig. 4d). The boundary layer mixing also shows a significant peak at the coldest temperatures near −2°C associated with deep water formation in the deep mixed layers at high latitudes.
Other sources of mixing, notably double diffusion, make a substantial contribution only at colder temperatures below 4°C (brown line in Fig. 6). Finally, mixing driven by the bottom-intensified internal tide scheme is relatively uniform in temperature space (gray line in Fig. 6) despite the parameterization giving a highly heterogeneous spatial distribution of mixing (Fig. 4e).
5. Sensitivity to background diffusivity and resolution
a. Background diffusivity
The magnitude of the diathermal heat transport was found to be sensitive to the level of background diffusion included in the model, and also more weakly to the model resolution. When a strong, spatially uniform background diffusivity of 10−5 m2 s−1 (a typical value used in many ocean models) is used, background diffusion dominates over other sources of mixing (cf. dotted green line to other solid lines in Fig. 6), and the total diathermal heat flux due to vertical mixing across 21.5°C is increased from −0.93 to −1.27 PW (cf. dotted and solid red lines in Fig. 7). In contrast, when background mixing is reduced to 0, the vertical mixing flux across 21.5°C is reduced to −0.47 PW (cf. thick dashed and solid red lines in Fig. 7). The resulting changes in the ocean’s temperature structure drive only minor changes to the other components of vertical mixing (not shown). However, the reduction in temperature gradients when background mixing is increased reduces the response of the background mixing heat flux to the background diffusivity (e.g., compare dotted and dashed green lines in Fig. 6, where the order-of-magnitude change in results in only a factor of ~6 change in the heat flux).
Changes in vertical mixing also induce changes to the other terms in the diathermal heat budget. When vertical mixing is increased, numerical mixing throughout much of the warmer temperature range decreases (cf. blue lines with red lines in Fig. 7 between 15° and 28°C). This may be a consequence of reduced vertical temperature gradients under increased background mixing, which reduces the spurious downgradient numerical smoothing of vertical gradients as those gradients are better resolved. This also suggests that when background mixing is weak (or zero; thick dashed lines in Fig. 7), numerical mixing takes the role of background mixing, acting to limit the sharpness of vertical temperature gradients in the thermocline. However, it is also possible that changing the background diffusivity alters the horizontal temperature gradients and thus numerical mixing associated with lateral advection.
Changes in numerical mixing do not completely compensate for changes in vertical mixing. Instead, the surface forcing diathermal heat flux also changes, increasing when the background diffusivity is increased (cf. black lines in Fig. 7). Comparing the m2 s−1 and cases, which show the largest difference, the change in surface forcing arises because the additional background mixing results in weak but widespread surface cooling throughout the tropical regions, in particular in the central Pacific [contours in Fig. 8, consistent with Zhu and Zhang (2018)]. The net surface heat flux into the ocean is then increased in these regions due to its dependence on SST (color in Fig. 8), resulting in a subsequent increase in the total flux of heat into waters warmer than 21.5°C. Neither the total surface area of the ocean warmer than 21.5°C nor the influence of solar penetration shows significant changes (not shown). Note that these responses may be altered if the same change in background diffusion was performed in a coupled model (e.g., Richards et al. 2009), as feedbacks onto the atmospheric state are not permitted in our model setup.
There are also minor changes in the internal heat content tendency (magenta lines in Fig. 7), associated with the differing equilibration times. The zero background diffusivity case has been equilibrated for 513 years at which point it shows a net heat flux into the ocean of 0.07 PW. In contrast, MOM025 Control and the 10−5 m2 s−1 case (which both show a larger net heat flux into the ocean of ~0.18 PW) have been run for only 80 years since the change in background diffusivity, which is enough to largely equilibrate the ocean warmer than ~10°C in the sense that the tendency in heat content at these temperatures is weak (magenta lines in Fig. 7).
Finally, the global diathermal heat budget is relatively insensitive to resolution, as shown by comparing the 1/10°, 75 vertical levels MOM01 configuration with MOM025 (cf. thin and thick dashed lines in Fig. 7, both of which have zero background diffusivity). Interestingly, the overall heat content tendency is smaller in MOM01 (0.05 PW; thin magenta dashed line in Fig. 7 at −2°C) than in MOM025 (0.07 PW). The drift may be weaker in MOM01, despite the much shorter equilibration time (72 years compared to 513 years for the MOM025 case), because of the influence of mesoscale and submesoscale eddies. Eddies are better resolved in MOM01 and generally act to reduce vertical heat transport into the deep ocean, keeping the surface layers warmer and reducing the air–sea heat flux (Griffies et al. 2015). Alternatively, this could be because MOM01 forms and maintains more Antarctic Bottom Water than MOM025 due to its higher vertical and horizontal resolution around Antarctica (and thus reduced numerically induced entrainment; Winton et al. 1998), resulting in reduced deep warming.
MOM01 also shows reduced numerical mixing compared to MOM025. This difference could be due to 1) a more realistic cascade of tracer variance from mesoscale stirring through to dissipation by vertical mixing (Smith and Ferrari 2009) in MOM01 and/or 2) reduced spurious diapycnal mixing in MOM01. As discussed in the previous section we are not able to separate the spurious diapycnal and along-isopycnal components of the numerical mixing, and so we cannot further evaluate the causes of model differences.
6. Seasonal cycle
a. Surface forcing and tendency
We now examine how the diathermal heat fluxes vary over a seasonal cycle. Spatially, the surface heat flux undergoes a strong seasonal cycle associated with the north–south migration of shortwave insolation (Figs. 1b,c,d). The globally integrated surface heat flux also varies seasonally by ±5 PW (dotted black line in Fig. 9a). While this globally integrated surface heat flux (equivalent to Fig. 10a at the coldest temperature) generally has a single seasonal peak in DJF, the total flux of heat due to surface forcing into all but the coldest temperature classes generally has two seasonal peaks in February and September–October that are well separated by a decrease in heating in November [red peaks in Fig. 10a; also compare example time series at 21.5°C (black solid line in Fig. 9a) to the globally integrated flux (black dashed line in Fig. 9a)]. The seasonal cycle of is driven not only by the seasonal variations in the surface heat flux but also by the seasonal variations in the temperature structure of the sea surface and upper ocean. The seasonal cycle in surface forcing drives a seasonal cycle in the internal heat content tendency that is similar in most respects (cf. Figs. 10b and 10a). The difference between the surface forcing and heat content tendency is accounted for by the seasonal cycle in mixing, discussed next. It is worth highlighting that the annual-mean surface forcing, peaking at 1.6 PW near 21.5°C (Fig. 3), arises from the rectification of a large seasonal cycle of order ±5 PW (Fig. 10a; see also black solid line in Fig. 9a).
The seasonal cycle of the internal heat content tendency can also be estimated from Argo observations of water-mass structure (Fig. 10c). This seasonal cycle agrees reasonably well with the modeled seasonal cycle (cf. Figs. 10b and 10c; also compare solid and dashed magenta lines in Fig. 9a), particularly at warmer temperatures, giving us confidence in the model results. There are some differences, particularly at cold temperatures where the Argo seasonal cycle is stronger and acquires two warming peaks in February and November. These differences may be due to the relatively sparse Argo data in the high latitudes (particularly south of 60°S). Also, the lack of data below 2000 m means that any seasonal adiabatic heaving of isotherms at these depths manifests here as water mass changes.
It is theoretically possible to directly estimate the seasonal cycle of mixing from observations by quantifying from observational air–sea heat flux products and taking the residual with the heat content tendency estimated from the Argo ocean temperature observations. However, preliminary analysis has indicated that uncertainties in air–sea flux products are presently too large for this purpose (see also Toole et al. 2004; Valdivieso et al. 2017).
b. Vertical mixing
While smaller than the seasonal cycles in surface forcing and tendency, there are also seasonal variations in the heat fluxes due to explicit vertical mixing (Fig. 11a) and numerical mixing (Fig. 11b). For vertical mixing, these seasonal variations can be further decomposed into the contributions from different processes (Fig. 12). At cool temperatures, below ~15°C, the seasonal variations in vertical mixing are mostly associated with variations in boundary layer mixing (Fig. 12c), due to seasonal variations in surface winds and the near surface temperature gradients. There is a strong peak during Southern Hemisphere winter at the coldest temperatures, near −1.5°C, associated with deep mixed layer convection driven by strong surface cooling in this season. Variations in background mixing also contribute to the seasonal cycle, with two seasonal peaks in February–March and August–September throughout most of the temperature range (Fig. 12a). In contrast, the seasonal cycle in bottom-intensified internal tide mixing is weak, with a slight seasonal peak in August (Fig. 11d).
Over warm temperatures between 17° and 27°C the vertical mixing–driven heat flux peaks in June–October and is weakest in April and May (Fig. 11a; also see example time series at 21.5°C, red line in Fig. 9b). Both shear-driven mixing (Fig. 12b; see also purple line in Fig. 9b at 21.5°C), and boundary layer mixing (Fig. 12c; see also blue line in Fig. 9b at 21.5°C) show large seasonal cycles at these temperatures, with peaks in the months of June–September. In particular, shear-driven mixing varies from near zero in February–May to a maximum of −0.25 PW in August. These seasonal variations are associated with the seasonal cycle of the tropical trade winds over the eastern Pacific (Horel 1982), and the associated wind- and shear-driven mixing in the upper Equatorial Undercurrent. Recent observational studies have indicated that rapid surface cooling in the months of June–September in the eastern equatorial Pacific cold tongue is associated with intense, seasonally varying turbulent mixing (Moum et al. 2013, their Fig. 3a). The seasonal variations in shear-driven mixing in this study are consistent with these observations and drive 50% of the global seasonal cycle in the flux of heat across the 21.5°C isotherm due to vertical mixing (cf. purple and red lines in Fig. 9b).
c. Numerical mixing
The seasonal cycle of the numerical mixing heat flux gives us some indications as to its nature and location. At warmer temperatures between 17° and 27°C the seasonal cycle of numerical mixing is similar to that of the explicit vertical mixing (cf. Figs. 11b and 11a). This suggests that much of the numerical mixing may occur in the tropical and subtropical thermocline, potentially due to spurious mixing due to the presence of high-frequency velocity and temperature variability, and high vertical temperature gradients, in the later part of the year (e.g., due to tropical instability waves and upper Equatorial Undercurrent variability; Jing et al. 2014; Holmes and Thomas 2015; Moum et al. 2009). However, there are also some notable differences between the numerical mixing and vertical mixing seasonality, particularly at cooler temperatures. Here, other processes may be contributing, such as the dissipation of along-isopycnal temperature variance generated by eddies through numerical processes that in lower-resolution models would be parameterized using a Redi diffusivity. Another possibility is numerical mixing associated with strong horizontal velocity variance and temperature gradients (e.g., Gibson et al. 2017), which in particular may occur in the year-round eddying western boundary current areas. The seasonal variations in numerical mixing in MOM01 (not shown) and MOM025 are similar, but the overall amplitude of numerical mixing is reduced in MOM01 (cf. thin and thick blue dashed lines in Fig. 7).
7. Discussion and summary
Using a Walin (1982)-based water-mass transformation framework we have quantified the processes contributing to the global transport of heat from warm to cold waters (or the diathermal heat transport) within a global ocean sea ice model. We introduced the concept of internal heat content [ in Eq. (6)], which is independent of the arbitrary reference temperature typically used to define heat content, and which also obeys a simpler budget than the full heat content as the heat fluxes associated with across-isotherm volume fluxes are removed. In this way, the internal heat content budget allows a clean and robust quantification of the partitioning of the diathermal heat transport between different processes.
Over an annual average, the diathermal heat transport is dominated by a balance between surface forcing, which acts to warm the warmest waters and cool the coldest waters, and mixing, which acts to homogenize ocean temperatures. In MOM025 Control the peak diathermal heat transport reaches 1.6 PW at 21.5°C (Fig. 3), comparable to the peak meridional heat flux out of the equatorial regions of 1.9 PW.
The dominant diathermal heat transport mechanisms are summarized in Fig. 13. Heat enters the ocean predominantly in the equatorial oceans (particularly the eastern equatorial Pacific; Fig. 1) at SSTs above 21.5°C (orange surface heating arrows in Fig. 13). However, heat leaves the ocean at cooler temperatures, predominantly over the midlatitude western boundary currents at SSTs below 21.5°C (orange upward arrows in Fig. 13; also see Fig. 1). To maintain a steady state and achieve the ocean’s portion of the climate system’s meridional heat transport, heat must be transferred downgradient, out of the volume of waters warmer than 21.5°C toward colder temperatures. This diathermal heat transfer is achieved via diffusive mixing, as advection cannot drive a net global transfer of heat across isotherms in steady state (besides the small contribution of surface volume fluxes), and globally surface forcing transfers heat upgradient.
Diffusive mixing includes contributions from a number of different processes (Fig. 6). In the control simulation, background vertical mixing (green arrows in Fig. 13) is the largest term, peaking around 20°C where thermocline temperature gradients are strongest. Mixing within the surface boundary layer also contributes substantially (aqua arrows in Fig. 13; Fig. 6), peaking at approximately 26°C. Surface boundary layer mixing accounts for most of the seasonal cycle in vertical mixing at temperatures below 15°C (Fig. 12c). Vertical mixing associated with shear instability is almost exclusively focused in a small region above the Equatorial Undercurrent in the eastern equatorial Pacific (Figs. 4c and 5, and purple arrows in Fig. 13). This shear-driven mixing makes a globally significant contribution to the diathermal flux, accounting for half of the seasonal cycle in the vertical mixing heat flux across the peak temperature of 21.5°C (Figs. 9b and 12b; also see Moum et al. 2013). Bottom-intensified tidal mixing (gray arrows in Fig. 13) also makes a contribution that is relatively uniform across temperature classes and not focused only in the deep ocean at cooler temperatures (Fig. 6).
Finally, the method we have presented also provides a useful way to diagnose the impact of numerical mixing (resolved in temperature and time; Fig. 11b) in realistic model configurations without the need for targeted experiments. In our model simulations numerical mixing makes a significant contribution to the diathermal heat flux (blue line in Fig. 3), and was found to be sensitive to the background diffusivity included in the model (section 5a). The prominant role of numerical mixing is a consequence of configuration choices made to minimize explicit sources of mixing (e.g., our configurations do not include any explicit horizontal or along-isopycnal diffusion, and the background diffusivity was small in several cases considered). These choices allow some aspects of the ocean circulation (e.g., eddy variability) to be better represented and minimize the overall levels of mixing. Numerical mixing was reduced when both horizontal and vertical resolution were increased (Fig. 7), although the relative contributions of vertical versus horizontal grid refinement are not yet known. Previous work suggests that the horizontal resolution may be the most important factor (Ilicak et al. 2012; Gibson et al. 2017). Future work will aim to further establish the processes contributing to numerical mixing, their sensitivity to model resolution, and the consequences of mixing numerically versus explicitly.
Our study has focused on model simulations where mixing processes are parameterized, and there remains significant uncertainty around these parameterizations. The framework that we have introduced, based on the internal heat content budget, provides an intuitive and robust method in which to quantify the sensitivity of heat flow within ocean models to different parameterizations and parameter choices. Our work highlights key regions, such as the eastern equatorial Pacific cold tongue, where specific parameterized processes play an important role in global budgets and thus should be the focus of future parameterization and model–observation comparisons.
Our framework can also be used to understand ocean heat uptake and its variability. Adiabatic rearrangement of heat within the ocean by variability in ocean circulation is important for driving changes in regional sea level (Yin et al. 2010), global surface temperature (England et al. 2014), and the global air–sea heat flux (Roemmich and Gilson 2011; Johnson and Birnbaum 2017). However, these adiabatic changes are largely reversible given a reversal of ocean circulation changes. Our water mass framework filters out adiabatic rearrangement and specifically isolates diabatic redistribution of heat within the ocean. Diabatic redistribution cannot be easily reversed as it relies on slower diffusive processes within the ocean interior to return water masses to their original state. Diabatic processes are thus pivotal in setting the long-term rate of heat and carbon uptake by the ocean, and thus also in setting future projections and climate sensitivity of coupled models. Future work should focus on using frameworks such as this to analyze the sensitivity of modeled ocean heat uptake and its associated time scales to model parameters.
We thank A. Heerdegen, A. Hogg, S. Griffies, and P. Spence for assistance with the modeling and diagnostics. We thank S. Groeskamp, T. McDougall, and two anonymous reviewers for useful comments. This research was undertaken with the assistance of resources and services from the National Computational Infrastructure (NCI), which is supported by the Australian Government. This project was also supported by the Earth Science and Climate Change Hub of the Australian Government’s National Environmental Science Programme (NESP) and the Centre for Southern Hemisphere Oceans Research (CSHOR), a joint research centre between QNLM, CSIRO, UNSW, and UTAS. The model output and configuration files are available upon reasonable request to the authors. The Argo data were collected and made freely available by the International Argo Program and the national programs that contribute to it. (http://www.argo.ucsd.edu; http://argo.jcommops.org). The Argo Program is part of the Global Ocean Observing System.
Numerical Implementation of the Diathermal Heat Budget
In this appendix we discuss in more detail the numerical implementation of the heat budgets.
a. Model heat budget
In MOM025 and MOM01 heat conservation is approximated at each grid cell by way of conservation of potential temperature [see Griffies (2012) for more information]:
where J kg−1 K−1, kg m−3, is the (time-variable) volume of the grid cell, and is potential temperature (where by using a constant for our simulations assume that potential temperature is conserved, which is not strictly the case; McDougall 2003). This equation, equivalent to Eq. (1), represents the change in heat content within each grid cell with time in Watts. This heat content changes due to advection of heat by the resolved model velocity field and the eddy-induced circulation from the submesoscale parameterization of Fox-Kemper et al. (2008) [both of which are captured by A in Eq. (1)], explicit vertical diffusion and the nonlocal fluxes of heat within the mixed layer associated with the KPP parameterization [both included in M in Eq. (1)] and surface forcing [captured by F in Eq. (1)]. There is no lateral or along-isopycnal diffusion term as these processes are not active in either MOM025 or MOM01. The SW penetration term redistributes the solar radiation from the surface layer to depth. The SSH smooth term captures heat content changes due to a smoothing of the sea surface height (SSH) applied at each time step for numerical stability. The term P captures the effective heat flux associated with freshwater inputs from rivers, ice melt, precipitation, and evaporation, which alter the volume of a surface grid cell (through ) and hence its heat content [the integral of P over the volume gives the surface volume flux contribution to the diathermal volume budget in Eq. (7)]. The heat budget Eq. (A1) closes exactly over each grid cell.
b. Numerical implementation of temperature based binning and integration
The diathermal heat budget framework described in section 3 was implemented online, performing the temperature binning integral [e.g., Eq. (2)] of each of the budget terms in Eq. (A1) at every time step. However, the advection and tendency terms A and T suffered from noise when binned online, associated with the fast time scale of advective changes in the temperature at each grid cell (although their sum was still exactly equal to the sum of the other terms in the temperature-binned budget, ). To avoid these problems we instead used instantaneous snapshots of the temperature field at the beginning and end of each month (the period over which the online diagnostics are accumulated) to calculate the heat content tendency and volume across each month. The internal heat content tendency was then obtained by rearranging the time derivative of Eq. (9):
We tested several different alternative methods of calculating [and therefore , which is found by residual using Eq. (15)], all of which were very consistent (Fig. A1). In particular, using Eq. (A2) agrees well with the internal heat content tendency calculated offline using monthly averages of the Eulerian tendency diagnostic T and the temperature field, which ignores the rectification of submonthly processes (cf. solid and dashed magenta lines in Fig. A1). A third possible method for calculating is to use only the volume tendency from the temperature snapshots and then integrate it in temperature following Eq. (11). This method also agreed well, although suffered from significantly more noise and did not match the surface forcing when integrated over all temperatures, resulting in a numerical mixing heat flux that did not sum to zero over all temperatures (cf. dotted and solid blue and magenta lines in Fig. A1). Finally, we also tested the sensitivity of our results to the size of the temperature bins, which had little impact (cf. solid lines, using a 0.5°C bin size, and dashed lines with circles, using a 1°C bin size in Fig. A1; note that any difference is likely due to the difference in year simulated).
Also shown in Fig. A1 is the total heat content tendency for the same January month, which evidently shows a large amount of noise and variability from temperature to temperature. In contrast, the internal heat content tendency is much smoother and more robust, as the noise present individually in the heat content and volume tendencies largely cancel. This highlights the utility of examining the internal heat content budget [Eq. (12)] as opposed to the full heat content budget [Eq. (7)].
The present framework could be extended to a non-Boussinnesq fluid by considering the total mass, rather than volume, of fluid warmer than .
It is thought that these choices allow for the maximum preservation of tracer gradients by the model and thus in a sense allow the most from their resolution.