A key question in climate modeling is to what extent sea surface temperature and upper-ocean heat content are driven passively by air–sea heat fluxes, as opposed to forcing by ocean dynamics. This paper investigates the question using a climate model at different resolutions, and observations, for monthly variability. At the grid scale in a high-resolution climate model with resolved mesoscale ocean eddies, ocean dynamics (i.e., ocean heat flux convergence) dominates upper 50 m heat content variability over most of the globe. For deeper depths of integration to 400 m, the heat content variability at the grid scale is almost totally controlled by ocean heat flux convergence. However, a strong dependence on spatial scale is found—for the upper 50 m of ocean, after smoothing the data to around 7°, air–sea heat fluxes, augmented by Ekman heat transports, dominate. For deeper depths of integration to 400 m, the transition scale becomes larger and is above 10° in western boundary currents. Comparison of climate model results with observations show that the small-scale influence of ocean intrinsic variability is well captured by the high-resolution model but is missing from a comparable model with parameterized ocean-eddy effects. In the deep tropics, ocean dynamics dominates in all cases and all scales. In the subtropical gyres at large scales, air–sea heat fluxes play the biggest role. In the midlatitudes, at large scales >10°, atmosphere-driven air–sea heat fluxes and Ekman heat transport variability are the dominant processes except in the western boundary currents for the 400 m heat content.
Classically, the ocean has been considered as an integrator of atmospheric weather noise, leading to ocean variability on time scales longer than those of the atmosphere (Frankignoul and Hasselmann 1977; Barsugli and Battisti 1998, hereafter BB98). A recently presented challenge to this paradigm is that ocean-driven processes can dominate over atmosphere noise in governing ocean mixed layer temperature variability in certain regions and scales of time and space. These ocean-driven processes include intrinsic variability (i.e., mesoscale and submesoscale eddies and other frontal variability; Sérazin et al. 2015, 2017; Nonaka et al. 2016; Smirnov et al. 2014; Pak et al. 2016) and ocean dynamical processes forced by the atmosphere due to local Ekman advection (Larson et al. 2018), basin-scale variability (Kwon and Deser 2007), and frontal variability forced remotely (Qiu and Chen 2005).
Many of these processes are encapsulated in extensions to the stochastic BB98 coupled system by Wu et al. (2006), Zhang (2017), and Bishop et al. (2017). As explained by Zhang (2017), it is important to include both ocean damping terms and ocean noise in the stochastic framework, and Bishop et al. (2017) show that this gives rise to the observed positive covariability of ocean temperature and air–sea heat flux (out of ocean) in ocean-eddying regions, such as the western boundary currents and Antarctic Circumpolar Current (ACC). These are regions where the ocean is prominent in driving SST variability.
The idealized stochastic models are useful for interpreting observations and model results but have certain limitations. For example, they typically do not include Ekman advection in the ocean, a source of SST variability that is (partially) independent of air–sea heat fluxes but is essentially driven by atmosphere motions. An exception to this is Lee et al. (2008), who introduced an Ekman term to explain the fact that AMIP simulations with prescribed SST gave more air temperature variability in some midlatitude regions compared to coupled runs, which goes against the BB98 expectation that coupling introduces “reduced thermal damping.”
The role of Ekman advection has also been explored in coupled model experiments where the variability of surface wind stress has been decoupled from the ocean model component (and replaced by climatological wind stress) by Larson et al. (2018). They referred to the simulation as being “mechanically decoupled” and found that removing the wind-driven part led to reduced SST variability in midlatitudes, and almost zero variability in the deep tropics (due to removal of ENSO winds), but actually enhanced SST variability in the subtropics.
In a pair of papers closely related to this current study, Buckley et al. (2014, 2015) examined the processes driving ocean heat content variability in the North Atlantic in a data assimilating ocean model. In particular, the importance of surface heat flux forcing and Ekman and geostrophic advection was determined using an upper-ocean heat budget. They found that in the subtropical gyres, the former two processes dominated heat content variability, while in the Gulf Stream/North Atlantic Current region geostrophic variability played a major role, especially on time scales longer than a month. Those results, obtained from the 1° ECCO model, will be expanded upon in this paper, and may also be considered as a benchmark due to its use of observational data in the state estimate. Our expansion of their work includes consideration of spatial scale, as well as discussion of a climate model that resolves mesoscale eddies in most of the world’s oceans, and a standard-resolution climate model that parameterizes eddies.
The role of intrinsic variability in the ocean has previously been addressed using forced ocean model simulations by Sérazin et al. (2015, 2017), Pak et al. (2016), and Nonaka et al. (2016), and in observed data by Roberts et al. (2017). The latter paper in particular noted a dominant role of intrinsic variability [i.e., variability of ocean heat flux convergence (OHFC)] in full-depth heat content variability as well as a significant role in mixed layer temperature variability. Their paper computed the OHFC as a residual in the heat budget due to lack of detailed observations of the OHFC: in contrast, our approach in the current paper is to use full heat budgets from a model to augment findings from analysis of observations.
In this paper we examine the causes of ocean temperature variability in coupled climate models of different resolutions, and compare against observations where available. The paper is a companion to Bishop et al. (2017) and Small et al. (2019) and follows on from work by Kirtman et al. (2012) and Putrasahan et al. (2017). In those papers it was found that low-resolution climate models mainly show a response of SST to atmosphere-driven air–sea heat flux variability, while high-resolution models show a weak response. Instead the air–sea heat fluxes in high-resolution models damp the SST, which is also seen in observations in ocean-eddying regions. This paper expands on this body of work to investigate causes of upper-ocean heat content and relationship to surface fluxes and ocean dynamics, by combining inferences from observational data with full heat budgets of a climate model.
The main questions addressed here are the following:
How important is OHFC versus air–sea heat fluxes in driving upper-ocean heat content variability, and what is the sensitivity to spatial scale?
How well do climate models of different resolutions compare to observations of heat content and surface heat flux variability, where available?
A related question is as follows:
3) Why does the SST variability have such large magnitude in high-resolution models in western boundary currents (WBCs), despite the strong damping effect of air–sea fluxes?
The low-resolution simulation presented here used version 4 of the Community Climate System Model (CCSM4; Gent et al. 2011). CCSM4 comprises the Community Atmosphere Model, version 4 (CAM4; Neale et al. 2013), Community Land Model, version 3.5 (CLM3.5; Oleson et al. 2008), Community Ice Code, version 4 (Hunke and Lipscomb 2008), and the Parallel Ocean Program, version 2 (Smith et al. 2010; Danabasoglu et al. 2012). CAM4 uses the finite-volume dynamical core and is discretized on a regular longitude–latitude grid of approximately 0.5° spacing and has 26 hybrid pressure levels in the vertical. CLM3.5 runs on the same horizontal grid as CAM. POP2 is on a nominal “1°” grid (actually 1.11° in zonal direction and meridionally between 0.27° (at equator) and 0.54° at higher latitudes). It has 60 vertical levels, with 10 m grid spacing in upper 100 m. Air–sea fluxes are computed in the coupler with the Large and Yeager (2009) scheme. The simulation used here was a short run of 10 years initialized from year 863 of a long preindustrial control run of CCSM4, which saved additional variables for ocean heat budgets. The model simulation is simply referred to as low resolution (LR).
The high-resolution climate simulation uses a successor to CCSM4, namely, version 1.1 of the Community Earth System Model (referred to here as CESM1; Hurrell et al. 2013). The main differences with CCSM4 in our case is the use of the Community Atmosphere Model 5.2, using a spectral element dynamical core (Park et al. 2014), and the Community Land Model, version 4 (Lawrence et al. 2011). The grid spacing is approximately 0.25° in the atmosphere. The ocean has a tripolar grid, with spacing of [0.1°, 0.1° cos(latitude)] for the longitude and latitude, respectively (a Mercator projection with isotropic spacing), except for grid distortion near the poles. Small et al. (2014) described the main control simulation, run for 100 years under “present day” (year 2000) greenhouse gas conditions. In this present work we use an 8-yr extension, branched from year 90 of the control run, which had extra heat budget variables saved at a monthly frequency. We refer to the short extension as high resolution (HR).
Note that, in contrast to Small et al. (2019) and Kirtman et al. (2012), the low-resolution model used here differs not just in its ocean resolution but also in atmosphere resolution and in climate model version. The reason for this is that the above simulations were unusual in saving extra data for the full heat budget described below. A number of studies have shown that the nature of the kind of differences that we identify in this paper, between HR and LR, are almost solely due to ocean resolution (Bryan et al. 2010; Kirtman et al. 2012; Sérazin et al. 2015). In particular, Small et al. (2019) compared CESM1 with identical 0.25° spacing in the atmosphere and 1° versus 0.1° in the ocean. A key metric from that paper, namely, the correlation between latent heat flux and SST on monthly time scales, was repeated for our LR and HR simulations and gave very similar results to those in Small et al. (2019).
b. Observational products of surface heat flux, SST, and SSH
The Japanese Ocean Flux Datasets with Use of Remote Sensing Observations, version 3 (JOFURO-3), is the evolution of the original J-OFURO dataset (Kubota et al. 2002; Tomita et al. 2019). The new version is available for 1988–2013 with daily and monthly mean temporal resolution at 0.25° spatial resolution. The focus of this study is on the monthly mean output. The dataset is derived solely from satellite data except for 2 m air temperature taken from NCEP–DOE reanalysis (Kanamitsu et al. 2002). Daily averaged SST is an ensemble median of multiple satellite data sources and of Reynolds et al. (2007) SST. The fluxes are computed using the COARE 3.0 bulk flux algorithm (Fairall et al. 2003). We use the latter years 2002–13 when more satellite data were available. Full details of the dataset are given in Tomita et al. (2019). The J-OFURO3 dataset was found to give consistent results with other datasets in Small et al. (2019) regarding the metric of SST–surface heat flux relationship, but it is acknowledged that considerable uncertainty exists between products for other metrics such as heat flux variance and trend (Zhang et al. 2018).
Sea surface height (SSH) anomalies are obtained from the CNES AVISO products, which merge data from altimeters on board TOPEX/Poseidon, ERS-1, Envisat, Jason-1, Jason-2, and CryoSat-2 to give a continuous record from 1992 (Ducet et al. 2000). The product used here was a monthly gridded 0.25° dataset (see acknowledgments).
c. Upper-ocean heat budget analysis
The heat budget at a particular level is written as
where T is potential temperature; ρ0 is a reference ocean density; cp is heat capacity of the ocean; ∇ is the 3D gradient operator; u is the 3D velocity; κ is the vertical diffusivity, and κΓ is the KPP countergradient flux of temperature; and P includes convergence of parameterized transport, such as Gent–McWilliams (Gent and McWilliams 1990), Redi (1982), and submesoscale (Fox-Kemper et al. 2008) parameterizations; and HMIX is horizontal mixing. Vertically integrating to a depth H gives
where Q is the surface heat flux, with positive values denoting heat gain; Qp is the penetrative heat flux (accounting for the absorption profile with depth of shortwave radiation, now separated from the rest of the mixing term), denoted Qp,−H at depth H; term vi is the interior mixing at depth H. For LR, HMIX is zero and the above terms are known exactly. For HR, P is zero, and all the terms except HMIX [horizontal (biharmonic) diffusion] and penetrative shortwave radiation were saved, and the latter two lumped together as a small residual, while for LR the residual is negligible. Note that the tendency term is exact (i.e., difference between instantaneous SSTs at beginning and end of month); also the budget terms, where known, are exact (e.g., for advection, covariance terms are accumulated each time step and averaged at the end of each month). Term vi will be referred to as interior vertical mixing (VMIX) and terms v and vi will be lumped together as VDIFF below. Further, terms ii and iii will be combined and referred to as OHFC.
The heat budget was initially computed for individual grid cells, and then the budget terms are linearly regressed onto the tendency term following Doney et al. (2007). Monthly anomalies are defined for all variables as departures from the average seasonal cycle of the simulation.
For the spatial-scale analysis of section 4, the budget terms on the 0.1° grid are smoothed with boxcar averages. (The total budget term is smoothed, not the state variables that go into the budget term.) Although the boxcar filter is very simple (and does have some problems associated with aliasing into higher wavenumbers; Emery and Thomson 2001), it retains the exact budget balance when taken over an integer number of grid cells. For more sophisticated techniques, the reader is referred to, for example, Arbic et al. (2014), Torres et al. (2018), and Laurindo et al. (2019) for distinguishing motions based on spectra in wavenumber and/or frequency space.
For simplicity, a smoothing of 1° full width refers to boxcar smoothing over a region of 10 × 10 grid points. As the model grid is on a Mercator projection, , where re is the radius of Earth and ϕ is latitude, the actual area of 10 × 10 points goes like cos(lat)2, and the area thus substantially diminishes toward high latitudes. In Fig. S1a in the online supplemental material, comparison is made of a 10° full-width smoothing done using the above definition [Mercator, i.e., 10° cos(lat)2 smoothing] with that for (10°)2. As expected, in higher latitudes, the latter gives a much bigger change than the former when compared against the original unsmoothed data, because the area of smoothing is greater, but for latitudes lower than 50° the differences are quite small. In the text we use the Mercator definition. Note also that near the model poles, the grid is deformed and so the smoothing area is not consistent, but we neglect that here. Finally, for the larger-scale smoothing (e.g., 10°), the results are plotted only every 5th or 10th point of the original grid to reduce computation time. This will not significantly affect the results as the budget results are generally coherent over scales larger than 5 or 10 ocean grid points.
3. Upper-ocean heat budget
It has been shown in Small et al. (2019) that in most regions of the globe, on monthly time scales, air–sea latent heat flux variability is driven by wind and air humidity variability, with the exception of the WBC regions. Further, in the Larson et al. (2018) mechanically decoupled low-resolution coupled model, in the midlatitudes and on monthly time scales, buoyancy fluxes alone (mainly due to heat fluxes) are responsible for 60%–80% of the SST variance away from frontal zones (Ekman heat transport variability contributes much of the remaining variability—see also section 3b).
Do these results imply that air–sea fluxes are the most important drivers of SST and heat content variability in the real world? In this section, the question is examined using a high-resolution climate model and compared to the standard-resolution model and to available observations. Specifically, we analyze upper 50 m heat content variability, which is typically dominated by mixed layer processes and is closely related to SST, and then a fuller depth heat content to 400 m, which, with a few exceptions (such as northern Atlantic high latitudes and sub-Antarctic frontal zone in winter), is much deeper than the mixed layer.
a. Heat budget to 50 m depth
A classic example of surface heat flux–driven upper-ocean temperature variability is shown in Fig. 1a, from the LR simulation, at a location in the subtropical North Atlantic. The tendency term is dominated by the VDIFF term (sum of surface heat flux and interior mixing). In contrast, the OHFC term is quite weak and not in phase with tendency.
Two further example time series from the LR run show different balances. In the Gulf Stream region (Fig. 1b), both OHFC and VDIFF contribute substantially. This is also the case in the Southern Ocean (Fig. 1c), and here it is notable that OHFC and VDIFF often act in concert—that is, they appear to be correlated, and both are correlated with the tendency.
The corresponding time series from HR show, in general, a larger role for OHFC (Fig. 2). The heat content tendency in the subtropical North Atlantic is now due to a combination of VDIFF and OHFC (Fig. 2a), while OHFC is dominant in the Gulf Stream (and note the much larger range of heat content tendency in Fig. 2b than Fig. 1b). In the Southern Ocean the HR results (Fig. 2c) are more similar to LR (Fig. 1c) except that the three terms (tendency, OHFC, and VDIFF) are less clearly correlated than in LR. This is because the OHFC is now a combination of Ekman and intrinsic processes—see below.
These time series results are now generalized to extend to the whole globe, by computing regressions between the budget terms and the heat content tendency at each grid point. Here, regression values near one indicate dominance of that term in the heat budget, while negative values indicate that the term must be counteracted by another term (Doney et al. 2007).
For a shallow depth integral (50 m here) the LR heat content tendency is mainly driven by the VDIFF term over large regions of the globe, especially the subtropical gyres (Fig. 3a). However, OHFC dominates in the eastern deep tropical Pacific and Atlantic (Fig. 3b; we loosely refer to deep tropical as within 10° of latitude of equator). Meanwhile, in the Southern Ocean and midlatitude North Pacific and North Atlantic, both OHFC and surface heat flux can be important in the vicinity of SST fronts. The VDIFF–heat content tendency regression is dominated by the surface heat flux contribution to VDIFF as seen in Fig. S2 (the regression between heat content tendency and surface heat flux).
In contrast, the HR heat content tendency is dominated by OHFC (regression >0.6 over most of the globe; Fig. 3d). The regions where the regression drops to between 0.4 and 0.6 are where the VDIFF term becomes important and includes the subtropical North and South Atlantic, northeast and southeast Pacific, and parts of the tropical warm pool (Fig. 3c). Note that the VDIFF regression becomes negative in parts of the Southern Ocean—likely where the mixed layer depth (MLD) often exceeds the 50 m and so the VMIX term at 50 m becomes important. This is also confirmed in Fig. S2, which shows that the regression between heat content tendency and surface heat flux is positive everywhere.
How realistic are the HR and LR results for 50 m heat content?
Using observations, it is difficult to estimate the OHFC on monthly time scales and with sufficient resolution to compare with the model. However, comparisons can be made for terms involving the surface heat flux term, which dominates VDIFF and has been estimated from satellite observations, in situ data, and reanalysis. Further, we use observed SST as a proxy for 50 m temperature (which is found to be a reasonable approximation in the model; Fig. S3 compares the metric discussed below using SST with that based on the 50 m heat content—they give similar results.) The correlation between SST tendency and surface heat flux for the observational analysis J-OFUROv3 was shown in Small et al. (2019),1 as was the same quantity from HR and CESM1 with low-resolution ocean. The latter exhibited systematically overstrong correlation over the whole globe: for example, reaching magnitudes of 0.6 to 0.8 in the subtropical gyres (compared to 0.4 to 0.6 in J-OFUROv3 and HR) and between 0.5 and 0.6 in the Southern Ocean (compared to typical values between 0 and 0.4 in J-OFUROv3 and HR). J-OFUROv3 was found to be similar to other products such as OAFlux and SeaFlux. We can conclude that HR has realistic drivers of ocean heat content tendency in the upper 50 m, on monthly time scales, but the comparable model with low-resolution ocean is too strongly heat flux driven.
b. Role of Ekman heat transport
In HR, OHFC is an important term in the 50 m heat balance, and it is expected that much of the OHFC variability is due to intrinsic oceanic variability, that is, variability that is not directly and locally forced by the atmosphere. The OHFC is a less important term in LR, but in the deep tropics it often dominates, and in some regions of the extratropics it contributes up to a half of the heat content tendency (Fig. 3b). Here we focus on the extratropics and ask what drives the OHFC variability in CESM.
A first clue is that some of the largest OHFC regressions in LR occur close to ocean fronts such as the Gulf Stream, Kuroshio Extension, and Agulhas Return Current. This is also true of HR—but the dynamical reasons are different, as we show below.
It is expected that Ekman heat transport will play a role in regions of variable winds and ocean temperature gradients. Anomalies of Ekman heat transport (EHT) can be written as
where it is assumed that the temperature in the Ekman layer is equal to the SST. The anomalies of EHT are due to anomalous winds acting on climatological SST gradients [first term of (2)], or climatological winds acting on anomalous SST gradients [second term of (2)]. This separation of these terms is explored below, but first the full EHT anomaly is diagnosed in the model simulations.
Analysis of LR reveals that much of the OHFC variability is due to EHT, as seen by regressing (at each grid point) the two terms: for example, under the North Pacific storm track and a large part of the Southern Ocean, there are large regions where EHT contributes 80% of the total OHFC (Fig. 4a). As a consequence, the regression of ocean heat content tendency (to 50 m depth) and EHT (Fig. 4b) looks very similar to that between ocean heat content tendency and OHFC, but with somewhat weaker values (cf. with Fig. 3d). In contrast, the equatorial band of the eastern Pacific and Atlantic exhibits a different character of OHFC variability, with only a small part driven by Ekman transport. This is expected in this region where a nonlocal response to wind anomalies through wave propagation is an important and well-understood component of the variability in upper-ocean heat content (cf. Figs. 4b and 3b).
Moving on to HR, the EHT contribution to OHFC is, in general, much weaker than seen in LR (cf. Figs. 4a,c), but there are some small localized regions where 80% of the total OHFC is contributed by EHT, such as in the Southern Ocean and south of the equator in the eastern basins. In contrast, the Ekman terms plays a notably small role in the western basins such as the Gulf Stream, Kuroshio Extension, as well as the Agulhas Return Current (Fig. 4c). This last result confirms the importance of intrinsic variability, such as quasigeostrophic eddies in HR, but also some larger-scale motions as discussed in section 4.
The regression of ocean heat content tendency and EHT in HR (Fig. 4d) is more complex and less coherent in sign than in LR (Fig. 4b), and globally the contribution of EHT to heat content tendency in HR is much weaker than that of full OHFC (Fig. 3d) confirming the importance of motions arising from intrinsic oceanic variability.
The decomposition of the Ekman term reveals interesting differences between the simulations, which is diagnosed by regressing OHFC separately onto the Ekman part due to wind stress anomalies (first term in brackets on RHS of 2), and onto the part due to SST gradient anomalies (second term in brackets on RHS of 2). For LR, the wind stress anomaly part dominates over SST gradient part in all regions except the eastern tropical Pacific and Atlantic Oceans (Figs. 5a,b). In contrast, the SST gradient anomalies dominate the EHT in HR (inferred from Figs. 5c and 5d), including in the same eastern tropical Pacific and Atlantic Oceans, but also in the high-latitude Southern Ocean, due to the large SST gradients and their variability that are expected in a high-resolution model. Interestingly, the EHT variability due to SST gradients in the WBCs and ARC regions does not contribute substantially to the OHFC variability (see Fig. 5d) presumably because it is overwhelmed by the intrinsic (quasi)geostrophic variability in those regions.
In the midlatitude storm-track regions, strong westerly wind anomalies can induce EHT by creating meridional flow across SST gradients of zonal currents, and more generally, if the wind variability is aligned parallel to an ocean front, an EHT anomaly is formed. As the temperature gradient here is from warm subtropical water to cold subpolar water, a negative meridional gradient, the sign is such that it creates relative warming under weaker westerly winds and relative cooling under stronger westerly winds. It is thus in phase with the surface heat flux effect where stronger winds create higher surface heat loss, as shown in Fig. 1c. In contrast, in subtropical easterly trade wind regions stronger easterly winds cause relative warming through EHT and thus act against surface heat flux effects (Larson et al. 2018)—but this is often hard to see as the Ekman term is small there.
A traditional approach is to separate ocean flow into two dominating parts: the Ekman, wind-forced part (giving rise to EHT), and the geostrophic, ocean dynamic part. The above results show that much of the OHFC variability in LR is due to EHT. In contrast, in HR the Ekman part only plays a significant role in some localized regions of the Southern Ocean and the southeastern basins. It follows that geostrophic flow dominates most of the OHFC in HR. Indeed, the regression maps between OHFC and heat content tendency in HR (Fig. 3d) have a very different spatial pattern to that seen in LR, with instead a strong resemblance to maps of EKE as seen, for example, in McClean et al. (2011). Does this all mean the EHT is absent or weak in HR? This will be addressed in section 4 below, where its importance is found to depend on spatial scale.
c. Heat budget to 400 m depth
For a deeper depth integral, such as to 400 m (actually 426 m in our case due to the vertical grid spacing of the model, but referred to as 400 m below for simplicity), the surface heat flux term generally plays a weaker role in driving heat content variability, as expected since this depth is much deeper than most mixed layers in the world’s oceans, excepting high latitudes in winter. Having said that, there are still evident differences between LR and HR in the heat content tendency as discussed next.
For the subtropical North Atlantic location, OHFC now plays a bigger role in LR than seen in the upper 50 m (cf. Fig. 6a with Fig. 1a), but VDIFF is still an important term. However, in the Gulf Stream and Southern Ocean locations in LR, VDIFF is typically smaller than OHFC (Figs. 6b,c). The latter dominates the tendency, but it can be noted again that VDIFF often acts in concert with the OHFC, with an influence of EHT, but much weaker than in the 50 m integral (not shown).
For the high-resolution model, time series of heat content tendency to 400 m at the same three locations show a complete dominance by OHFC, and VDIFF is a very small term (Fig. 7). Further, the OHFC in HR is an order of magnitude bigger than in LR in the subtropical North Atlantic and Gulf Stream locations (Figs. 6 and 7). These results are generalized to the globe in Fig. 8. Regressions between heat content tendency and surface heat flux (i.e., VDIFF) in HR are near zero (Fig. 8c) and instead the variability globally is almost entirely driven by OHFC (Fig. 8d). In LR, this is also true for a wide swath of the tropics and subtropics (Fig. 8b), but in midlatitudes, both the surface flux (Fig. 8a) and OHFC are important to heat content variability.
d. How realistic are the HR and LR results for 400 m heat content?
Observational data of upper heat content tendency to ~400 m of sufficient resolution in time and space to compare with the model results are lacking. However, SSH represents the deep structure of the ocean, and observations of SSH variability from satellite have existed since 1992. The correlation between SSH and ocean heat content is reasonably high in observations (Willis et al. 2004) and very high in CESM-HR (Fig. 9), especially in regions of strong eddy variability.
The instantaneous correlation between SSH anomalies and surface heat flux anomalies is similar to that between SST and surface heat flux shown in Small et al. (2019): high magnitude of correlations in ocean frontal zones in HR and in observations, and much weaker correlations for LR (Fig. 10). Further, in LR there is a widespread reversal of sign to positive in the Southern Ocean (Fig. 10c), something that is seen only in isolated regions of observations and HR, notably in the sub-Antarctic zone of deep mixed layers (Figs. 10a,b).
Further, Fig. 10 uses net surface heat flux for the observation (Fig. 10a) but latent heat flux for the model (Figs. 10b,c). In Fig. S4 it is shown that there is very little difference between using latent heat flux and net heat flux. Further, as expected from Fig. 9, very similar correlation patterns with surface heat flux were found in the model when considering ocean heat content integrated to 400 m depth, instead of SSH (Fig. S5).
Most important for the heat budget results of section 3c is the relationship between SSH tendency and surface heat flux, an allowable proxy for the heat content tendency–surface heat flux relationship since SSH and heat content are so highly correlated. It can be seen from Figs. 11a and 11b that in both observations and HR there is a very small correlation between these variables in most part of the globe, whereas LR has quite high correlations, for example, in the eastern Atlantic off Europe, consistent with the heat content regressions of Figs. 8a and 8c. The results confirm that at the grid scale, the tendency of 400 m heat content and SSH is not driven by air–sea heat fluxes in HR and observations, while in LR there is still an influence of surface heat fluxes.
4. Scale dependence of upper-ocean heat budget
All the above results were presented for gridpoint relationships, which do not distinguish between processes acting at different spatial scales. The dependence on spatial scale of processes driving SST variability is now explored by successively smoothing the high-resolution ocean heat budget terms with boxcar filters of various widths (see section 2c for details). The HR and smoothed results are shown in Figs. 12–16 and compared with LR results. The smoothing width varied from 1° to 10°—wider smoothing than 10° was not attempted due to problems near the boundary that particularly affect the results in WBCs, eastern boundary currents, etc. Further, the 1° smoothing of budget terms led to results quite similar to the unsmoothed data and are thus not shown. (Note that this implies that gridscale noise is not a dominant process in the unsmoothed regression maps of section 3—otherwise, 1° smoothing would have a substantial effect.)
For ocean heat content tendency and VDIFF averaged to 50 m depth, smoothing of HR to 5° converts the regression pattern to something similar to LR over most of the globe (Fig. 12). The implication is this: for example, for the subtropical gyres, LR heat content tendency is mostly driven by air–sea heat flux (atmosphere driving) and so is HR on large spatial scales. However, in WBC and ACC regions, which have near-zero regressions in HR, 7° to 10° smoothing is required to match HR values to LR (Fig. 12).
A similar story applies to ocean heat content tendency and OHFC to 50 m (Fig. 13)—5° smoothing of HR recovers the LR result except in the WBC regions where smoothing to at least 10° is required (e.g., in ACC). As discussed in section 3b, the OHFC term in LR under the midlatitude storm tracks is mostly the Ekman transport: the overall similarity of LR with 10° smoothed HR then suggests that it is also driven by the Ekman heat transport on large scales.
This latter point was confirmed explicitly by applying the spatial smoothing to the EHT anomaly term in HR, as described in section 2c for other budget terms. [Note the smoothing was applied to the first and second terms of RHS of (2), or their combination, rather than being applied to wind stress anomalies and SST gradient anomalies separately.] The regression of OHFC and EHT (both 10° smoothed from HR; Fig. 14a) looks similar to the LR result (Fig. 4a), as does ocean heat content tendency and EHT (cf. Fig. 14b with Fig. 4b). Further, when the 10° smoothed EHT was decomposed into the parts due to wind stress anomaly and part due to SST gradient anomaly (Figs. 14c,d), the results look very similar to those for LR (Figs. 5a,b) with the wind anomaly part dominating. The only exception to the general similarity of the Ekman term in 10° smoothed HR to that in LR is in the western boundary currents and ACC, where the Ekman term in HR contributes negligibly to the OHFC even after 10° smoothing.
For averages to deeper depth of 400 m, and for ocean heat content tendency and VDIFF, smoothing of HR to 7°–10° (Fig. 15) reduces the regression pattern to something similar to LR in most regions, (i.e., OHFC-driven in tropics: a mixture of OHFC and surface heat flux driving in extratropics), but in WBCs the regression values even after 10° smoothing are weaker than in LR. The dominance of ocean dynamics in these regions in HR at scales of 10° is seen more clearly in the OHFC results (Fig. 16)—whereas smoothing to 7°–10° recovers the low-resolution result in most regions (Fig. 16, bottom right), in the WBCs and ACC the regression of heat content tendency with OHFC is still higher than 0.8 compared to typical values of 0.4 to 0.7 in LR. Further, the spatial pattern of regressions in WBCs is quite different in HR smoothed to 10° and in LR, suggesting that the large-scale processes acting in HR are not local Ekman processes but more likely coherent shifts of the WBC systems due to internal dynamics (Sérazin et al. 2015, 2017; Nonaka et al. 2016) or due to remote forcing from the atmosphere (e.g., Qiu and Chen 2005).
It would be optimal to ground truth these smoothed heat budget results with an observed proxy, as was done in sections 3a and 3d for the gridpoint budgets. In section 3a, SST tendency and LHF from J-OFUROv3 was used to benchmark (Small et al. 2019). However, when smoothing the same J-OFUROv3 data (Fig. S6), it is found that the correlations between SST tendency and latent heat flux are generally weaker in magnitude than in smoothed CESM-HR, and in unsmoothed CESM-LR (e.g., subtropical gyres have correlation magnitudes of 0.6 to 0.8 in the model, 0.4 to 0.7 in smoothed J-OFURO3). This may imply that the climate model is more atmosphere driven than J-OFURO at large scales—but it could also be impacted by any errors in J-OFURO due to sampling limitations, algorithms, etc. that will necessarily reduce correlations compared to the regularly sampled model data.
In section 3d SSH was used as a proxy for ocean heat content for gridpoint analysis, which was acceptable as correlations between SSH and 400 m heat content were frequently above 0.8 and 0.9. However, when the heat content of the model is smoothed spatially, we found that the correlation with model SSH (similarly smoothed) decreased quite rapidly (Fig. S7). This may relate to the fact that for small scales, the baroclinic ocean variability (such as eddies) dominates both the heat content and SSH variability, whereas on larger scales, barotropic variability can dominate, which would be reflected in SSH but not in heat content. The latter possibility was hinted at by animations of SSH (not shown). Consequently, it is not easy to ground truth the smoothed heat content results—having said that, the dominance of OHFC in WBCs at large scales is consistent with results from the relatively coarse grid analyses of Roberts et al. (2017) and Buckley et al. (2014), and especially with Sérazin et al. (2015), who, in a high-resolution ocean model forced by observed state, found that intrinsic ocean variability contributed a large fraction (from 30% to 80%) of large-scale sea level variance in WBC regions.
a. Relationship of intrinsic versus forced variability
The above findings can be compared to a recent analysis by Roberts et al. (2017), which surveyed the ocean heat budget based on observations. They found that although the air–sea heat flux was globally the primary driver of seasonal temperature variability in the mixed layer, and also dominated the full-depth-average temperature variability in midlatitudes, for interannual variability a very different picture emerged. The interannual variability of mixed layer temperature was more strongly correlated with OHFC than with air–sea heat flux, and the full-depth-averaged temperature variability was almost fully explained by OHFC. Further, the variance of OHFC in the mixed layer was similar in magnitude to the air–sea heat flux variance, while the full-depth variance of OHFC was much larger, reaching values of 10 000 W2 m−4 (equivalent to a standard deviation of 100 W m−2) compared to values of <400 W2 m−4 (equivalent to a standard deviation of 20 W m−2) for surface heat flux variance. These results are consistent with Fig. 6, although we use 400 m instead of full-depth averages.
The results are also consistent with recent analysis of the contribution of ocean internal variability to total ocean heat content variability using high-resolution ocean models. For example, Sérazin et al. (2015) compared the intrinsic ocean variability in a climatologically forced ocean simulation against the variability arising in a simulation using interannually and monthly varying forcing. They found that intrinsic ocean variability of sea level completely dominated the low-frequency (>18 months) and small-scale (length scale L < 6°) variability everywhere in the extratropics. For high frequency (<18 month) and small scales there was total dominance by intrinsic variability. For low frequency and large scales (L > 12°), significant fractions (>0.5) of the sea level variance in the Southern Ocean were due to intrinsic variability, while in the Northern Hemisphere this ratio reached 0.3 to 0.5 in many parts of the subtropics and midlatitudes. Sérazin et al. (2017) showed that intrinsic variability of ocean heat content (OHC) was as large as forced variability on interannual to decadal time scales over much of the Southern Ocean and parts of the Northern Hemisphere oceans—with the intrinsic contribution growing larger as the depth of integration for OHC was increased (Sérazin et al. 2017, their Fig. 1).
Note much of the previous two paragraphs refers to interannual and lower-frequency motions, whereas our paper is for monthly variability. Previous results (Buckley et al. 2014, 2015; Bishop et al. 2017; Small et al. 2019) indicate that ocean dynamics plays a larger role with increasing time scale in many regions of the globe including but not limited to WBCs.
The drivers of heat content variability in the North Atlantic based on the ECCO state estimate (Wunsch and Heimbach 2007) were examined by Buckley et al. [2014, 2015; see also Fig. 7 of Buckley and Marshall (2016) for a summary]. They explored the frequency dependence of the heat budget (as opposed to a space-scale dependence of the current paper) and integrated in depth to the maximum climatological MLD in the seasonal cycle. Consistent with our results, they found a dominance of air–sea heat flux and Ekman (combined in their “local” term) in driving subtropical variability, whereas a large role for ocean geostrophic motions was found in the Gulf Stream regions (equivalent to the intrinsic ocean dynamics in our results). They also found the damping of variability in the Gulf Stream region by surface heat fluxes, as found in Bishop et al. (2017) and Small et al. (2019).
In contrast to our results, Buckley et al. find that the geostrophic motions at WBCs only become important at time scales of 1 year or more, eventually contributing about 50% of the OHC variability for 10-yr time scales. Our current paper suggests that the intrinsic ocean dynamics already dominates for monthly time scales in WBCs (Fig. 3) and this is reasonably consistent with the observations shown in Small et al. (2019). A possible reason for this difference between HR and the ECCO results is the model grid, 0.1° and 1°, respectively. Indeed, the largest initial misfits of ECCO model to data occur in WBC regions [Fig. 2 of Wunsch and Heimbach (2007)], where low-resolution models have large biases. On the flip side, for large scales that are resolved by ECCO, the fact that ECCO assimilates data, and can thus be regarded as an estimate of the observed state, gives us confidence in our finding that heat content variability in WBCs is generally driven by ocean geostrophic flows.
b. Role of Ekman heat transport
Our results are also consistent with previous findings regarding the role of EHT. In both LR and HR, the local Ekman process combines with the local air–sea heat fluxes when the atmosphere plays the driving role: in midlatitude storm-track regions they act in concert (leading to an increase in SST variance when Ekman is included; Larson et al. 2018), whereas in the subtropics, EHT counters the surface heat flux effect (and hence reduces SST variance). Lee et al. (2008) explored this further in an extension of the BB98 stochastic model and found that coupling led to regional-dependent differences of air temperature variability compared with atmosphere-only simulations, and it was not the case that coupling led to “reduced thermal damping” everywhere. Of more direct relevance to our model results, the superposition of air–sea heat flux and Ekman can lead to misleading correlations between air–sea heat flux and heat content tendency—too high in the midlatitudes and too low in the subtropics—compared to more “truthful” results from regressions of the budget terms. Further, the fact that spatially smoothed HR gives similar OHFC regressions to LR—and hence is representing the EHT at large scales—is an encouraging result,2 and one that is not easily seen at the grid scale because of the masking by the small-scale geostrophic motions.
Another finding of this paper is that ocean frontal zones are regions of both enhanced Ekman heat transport variability and of intrinsic ocean variability, as the strong SST gradients promote both of these processes. For LR, the Ekman heat transport variability dominates, while for HR, the intrinsic variability is most important at small scales and Ekman variability dominates at larger scales.
c. Air–sea feedback and SST damping
The HR heat budget results of section 3a partly answers a major question: why does the SST variability have such large magnitude in HR in WBCs, despite the damping effect of air–sea fluxes? Figure 2b reveals that this damping term (VDIFF, dominated by the surface heat flux) is much smaller than the ocean heat flux convergence term that maintains the SST anomalies in the Gulf Stream. Thus, even though surface heat flux variability is much stronger in CESM-HR than CESM-LR over WBCs (Small et al. 2019), it is still not enough to strongly damp the SST anomalies. However, more analysis is needed to fully answer question 3, namely, quantifying the air–sea feedback term (Frankignoul et al. 1998) and how it relates to SST variability.
d. Limitations of standard-resolution models
The results of section 4 showed that over much of the ocean, it is sufficient to smooth the HR data to around 5°–7° to recover the same type of ocean heat content budget as LR. This suggests that the eddy effects are mostly removed when smoothed over that scale in this model and resolution, and it also implies that the EHT variability that is seen on the grid scale in LR is an underlying feature of HR (but is masked at small scales by eddy effects). However, over WBCs and ACC it was found that even a 10° smoothing of HR could not recover the LR results, especially for the 400 m heat content. Clearly 10° scales should be resolvable by the 1° LR ocean model, but the above results show that LR does not represent these scales sufficiently. This suggests that the small-scale eddy motions (which are resolved by HR but not by LR) significantly impact motions on scales of 10° and hence the latter are poorly represented by LR and other models with 1° ocean resolution or coarser.
The parameterization of mesoscale and submesoscale, that is included in the P term in LR, is a deterministic function of the background state and does not add any extra variability: rather it is intended just to give the net effect of an ensemble of eddies on the model state. Indeed, the Gent–McWilliams (Gent and McWilliams 1990) and Fox-Kemper et al. (2008) parameterizations do not take into account inverse cascades of energy back up to the large-scale. More recent parameterizations (e.g., Jansen et al. 2015; Grooms et al. 2015; Bachman 2019) are aimed at introducing stochasticity into the parameterizations and/or transferring energy back to large scales.
e. Limitations of high-resolution, ~0.1°, models
A limitation of the current study is that scales below the ocean mesoscale are not considered (except within the context of the submesoscale parameterization of LR, which, as with the mesoscale parameterization in LR, plays a negligible role in our metrics of interest). The smallest scales resolved in HR are ~20–30 km at high latitudes (assuming around 6 points needed to resolve a wavelength: the smallest grid cells are 3–5 km) and ~50–60 km in the tropics. Recent literature has found an important role for explicitly resolving some parts of the submesoscale wavenumber spectrum, such that vertical processes at these scales, and eddy-driven restratification, can be as important as other terms in the budget such as air–sea heat flux (Su et al. 2018; Yu et al. 2019; Klein et al. 2019). These findings were obtained from ocean models and limited observations; it remains to be seen how these processes affect a fully coupled system.
On monthly time scales OHFC plays a significant role in driving ocean heat content variability averaged to 50 m depth in the high-resolution climate model, and results are quite consistent with observations (Bishop et al. 2017; Small et al. 2019). In contrast, for standard-resolution climate models, air–sea heat fluxes, driven in turn by the atmosphere variability, dominate the 50 m heat content variability, except in some atmosphere storm-track and SST gradient regions where Ekman advection also plays a sizeable role, and in the deep tropics where ocean dynamics is mostly resolved and drives heat content variability.
For ocean heat content to a significantly deeper depth, 400 m, variability in the high-resolution climate model is almost purely driven by OHFC. The fact that the relationship between SSH (a proxy for deep heat content) and turbulent heat fluxes in observations is very similar to that seen in the high-resolution model supports this. For the standard-resolution model, OHFC is also the dominant player in the tropics and subtropics, while in the midlatitudes the surface heat fluxes play an equally important role.
All the above results refer to pointwise analysis of time series at the model horizontal grid scale. Even in regions of low EKE, such as subtropical gyres, the OHFC can dominate the local heat budget in observations and HR. However, exploring the spatial scales of the high-resolution CESM ocean heat budget reveals that for ocean heat content to 50 m, averaging over boxes of 3°–5° square brings the budget close to the low-resolution model and a mainly atmosphere-driven state except in WBCs where 10° smoothing is required. For 400 m ocean heat content, averaging over 5°–7° is required to gain the low-resolution model state, excepting that in some regions of WBCs and ACC the heat content is considerably more driven by OHFC even when averaged over 10° squares, compared to the low-resolution model. The implication is that LR is missing intrinsic ocean variability on scales of 10° or more that naively should be resolvable by the 1° model.
Although the quantitative scales discussed here are specific to CCSM/CESM, one should expect qualitatively similar results for other climate models run at eddy-resolving versus eddy-parameterized resolutions. However, further investigation would be needed of eddy-permitting models (not examined here, but see Sérazin et al. 2015, 2017; Roberts et al. 2016) to see to what extent the intrinsic motions dominate.
The broader implications of our results for the Coupled Model Intercomparison Project (CMIP) class (~1° here) model experiments are as follows: for monthly time scales, the CMIP models can capture the upper-ocean (400 m or less) heat content variability at 10° scales or more, except for deep heat content to 400 m in the WBCs that are key regions of intrinsic ocean variability. For smaller scales than 5°–10°, the CMIP models do not accurately represent the nature of the variability. However, the increasing use of higher-resolution climate models for select experiments (e.g., HighResMIP) allows for the possibility of proper representation of the nature and statistics of small-scale heat content variability. Further, for longer time scales than monthly, it is likely that ocean intrinsic variability becomes even more important, especially in the WBCs (Buckley et al. 2014).
This work was funded by NASA Physical Oceanography Awards NNX16AH60G and 80NSSC18K0769. Young-Oh Kwon and Hyodae Seo (WHOI) are thanked for discussing the results. John Truesdale (NCAR) is thanked for extending the high-resolution simulation to gather extra data. Sea surface height data were gathered from https://www.aviso.altimetry.fr/en/data/products/sea-surface-height-products/global/gridded-sea-level-anomalies-mean-and-climatology.html. The altimeter products were produced by Ssalto/Duacs and distributed by Aviso+, with support from Cnes (https://www.aviso.altimetry.fr). Computing and data storage resources, including the Cheyenne supercomputer (https://doi.org/10.5065/D6RX99HX), were provided by the Computational and Information Systems Laboratory (CISL) at NCAR. The control HR simulation was previously performed at NCAR–University of Wyoming Supercomputer Center and supported by DOE Scientific Discovery for Advanced Computing (SCIDAC) Award SC0006743. The extension was performed as part of the Blue Waters sustained-petascale computing project, supported by the National Science Foundation (Awards OCI-0725070 and ACI-1238993) and the state of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana–Champaign and its National Center for Supercomputing Applications. This work is also part of the “High Resolution Earth System Modeling Using Blue Waters Capabilities” PRAC allocation supported by the National Science Foundation (Award ACI-1516624). The CESM project is supported primarily by the National Science Foundation (NSF). This material is based upon work supported by the National Center for Atmospheric Research, which is a major facility sponsored by the NSF under Cooperative Agreement 1852977.
This article is included in the Climate Implications of Frontal Scale Air–Sea Interaction Special Collection.
Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JCLI-D-19-0295.s1.
We should expect Ekman heat transport to play a role at large scales regardless of model resolution.