Surface flux and ocean heat transport convergence contributions to seasonal and interannual variations of ocean heat content

,


Introduction
The ocean plays a crucial role in the global energy budget as it is Earth's primary heat reservoir on climaterelevant time scales [e.g., Palmer and McNeall, 2014;Von Schuckmann et al., 2016]. As such, variations in ocean heat content (H) play a fundamental role in global and regional climate variability [e.g., Jin, 1997;Meehl et al., 2011;Roberts et al., 2015] and the transient response to climate change [Kuhlbrodt and Gregory, 2012;Geoffroy et al., 2013]. In addition, ocean heat uptake is one of the major contributors to past, and predicted future, changes in global sea level due to the thermal expansion of seawater [Church et al., 2013;Kuhlbrodt and Gregory, 2012].
Regional variations in H are also important for the initialization of seasonal-to-decadal forecasts [Dunstone et al., 2011;Robson et al., 2012], with evidence that skillful predictions of societally relevant variables may be possible at lead times of several years or more in some locations [e.g., Griffies and Bryan, 1997;. Understanding the drivers of regional H variability is thus an important area of research to promote advances in our ability to make skillful climate predictions. For example, in regions where ocean variability is dominated by the response to stochastic forcing from the atmosphere, predictability beyond atmospheric decorrelation time scales is limited to the persistence of anomalies provided by the heat capacity of the ocean mixed layer. In contrast, heat content variations associated with ocean dynamical processes (e.g., advection of temperature anomalies by the mean circulation, baroclinic wave propagation, geostrophic circulation anomalies) offer the potential for skillful forecasts on seasonal-to-decadal time scales through accurate initialization of the ocean state [e.g., Dunstone and Smith, 2010;Robson et al., 2012;Roberts et al., 2016].
[2014], we consider heat content in the mixed layer (H mld ), which varies geographically and represents the part of the water column that communicates directly with the atmosphere above. Variations in H mld are (by definition) highly correlated with sea surface temperature (SST) changes and therefore particularly relevant for seasonal-to-interannual climate variability and prediction [e.g., Alexander and Deser, 1995;Cassou et al., 2007;Taws et al., 2011]. Second, we consider total column heat content (H tot ), the volume of which is generally dominated by the ocean interior that is away from direct atmospheric influences. H tot is particularly relevant for interannual changes in steric sea level in the low latitudes and midlatitudes that are dominated by the thermosteric response to wind-stress variations [K€ ohl, 2014;Forget and Ponte, 2015;Stammer et al., 2013;Roberts et al., 2016].
The importance of a dynamically active ocean is clear for seasonal-to-interannual variations of H mld associated with tropical modes of climate variability, such as the El Niño-Southern Oscillation (ENSO) [Jin, 1997;McPhaden et al., 2006;Mayer et al., 2014] and gives rise to useful predictive skill in a number of societally relevant variables, including sea level [Miles et al., 2014;Roberts et al., 2016]. Historically, a number of studies have suggested that the midlatitude oceans are essentially passive on seasonal-to-interannual time scales, acting primarily as a static heat reservoir that integrates stochastic atmospheric forcings [e.g., Frankignoul and Hasselmann, 1977;Frankignoul, 1985]. In this paradigm, there is limited scope for predictive skill originating from ocean dynamics in the midlatitudes on interannual time scales.
Previous studies have evaluated the nature of air-sea interactions in both observations and models by calculating simultaneous correlations between local SST anomalies and either surface flux or wind speed anomalies [e.g., Gulev et al., 2013;Kirtman et al., 2012;Chelton and Xie, 2010;Bryan et al., 2010]. These studies have emphasized the ocean influence on the atmosphere in regions of strong SST gradients and the sensitivity of this influence to model resolution [Kirtman et al., 2012;Bryan et al., 2010]. In the flux-SST framework, negative correlations indicate a negative feedback such that positive SST anomalies are associated with heat loss from the ocean, which is interpreted as oceanic forcing of the atmosphere. Note that throughout this manuscript we use the convention that air-sea heat fluxes are positive into the ocean. In contrast, zero (or positive) correlations are interpreted as evidence for atmospheric forcing of ocean temperatures. Gulev et al. [2013] performed such a correlation analysis for multidecadal reconstructions of latent and sensible heat fluxes for a limited region of the North Atlantic (40-608N 25-558W) and found that interannual variations of SST in this region reflected forcing from the atmosphere whereas decadal changes in SST must be controlled by the ocean.
Other studies [e.g., Buckley et al., 2014;Cunningham et al., 2013] have taken a budget-based approach to evaluating the drivers of ocean heat content variability and found a more significant role for ocean dynamics in the midlatitudes, particularly in regions with strong currents and substantial eddy activity. For example, several studies have used ocean model simulations or model-based ocean state estimates to examine North Atlantic heat budgets and found an important role for heat transport convergence anomalies in the generation of interannual variability [Grist et al., 2010;Piecuch and Ponte, 2012;Buckley et al., 2014Buckley et al., , 2015. Observation-based heat budget analyses have also concluded that interannual variability of ocean heat transports play a dominant role in the heat budget of the subtropical North Atlantic [e.g., Cunningham et al., 2013;Bryden et al., 2014;Kelly et al., 2014], particularly in the Gulf Stream where positive anomalies in upper ocean heat content have been shown to drive heat loss to the atmosphere [Dong and Kelly, 2004;Dong et al., 2007;Buckley et al., 2014]. However, many of the studies that have found a dominant role for ocean dynamics [e.g., Kelly et al., 2014;Grist et al., 2010;Cunningham et al., 2013] have focused on full-depth ocean heat budgets. In these cases, it is not clear that the diagnosed drivers of depth-integrated ocean variability are related to the processes that set variability in the mixed layer and therefore it is possible that budgets that include the deeper ocean will provide an inaccurate picture of the relevance of ocean dynamics for driving air-sea interaction.
Here we present an observation-based seasonal and interannual analysis of heat budget terms and evaluate the dominant drivers of both H mld and H tot for the global oceans covering the period 1985-2012. Our analysis of the mixed layer is similar to recent studies of the North Atlantic [Buckley et al., 2014[Buckley et al., , 2015 but extended to the global oceans and based on observational data sets rather than an ocean state estimate. We combine four surface heat flux products and two ocean heat content products using a novel Kalman smoother-based method and infer the regional contributions from ocean heat transport convergences as a residual. In addition, we quantify the locally forced contribution to ocean heat transport convergence Journal of Geophysical Research: Oceans 10.1002/2016JC012278 variability due to Ekman transports and calculate correlations between H mld and surface flux contributions to identify regions where the ocean is driving the atmosphere on interannual time scales. Section 2 describes the data and methods used in this study. Section 3 presents the results of our heat budget analyses and identifies the dominant drivers of interannual variability of ocean heat content. Section 4 discusses the sensitivity of our analysis to the evolving ocean observing system. Section 5 summarizes our main results and conclusions.

Heat Conservation
We calculate the heat content (H) of ocean layers between the surface g and a depth of D using

Tdz
(1) where T is ocean temperature, q 0 is a reference density, c p is the heat capacity of seawater, and D is either the maximum depth of the seasonally varying mixed layer derived from the Kara et al. [2000] monthly climatology (D mld ; Figure 1a) or the ocean bottom (D b ). The Kara et al. [2000] mixed layer is defined using a density criteria that gives the optimal estimate of turbulent mixing penetration. Therefore, when using using D mld , we are considering the heat content in the fraction of the upper ocean that is available to interact with the atmosphere on seasonal-to-interannual time scales (H mld ). This is the same definition of the upper ocean used for heat budget analyses by Buckley et al. [2014] and is chosen such that mixing of heat across D mld can be considered negligible in most regions [Buckley et al., 2015]. When using D b , we are considering the column-integrated heat budget that is important for variations in thermosteric sea level. Heat budgets for both layers satisfy conservation equations of the form (b-f) Amplitude of the climatological seasonal cycle for (b) net surface heat fluxes Q net , (c) heat storage rates in the mixed layer @ @t H mld , (d) heat storage rates over the full depth of the ocean @ @t H tot , (e) heat transport convergence integrated over the mixed layer C mld , and (f) heat transport convergence integrated over the full depth of the ocean C tot . Seasonal cycles are calculated using the ensemble mean of data products for the period 2006-2012.

Journal of Geophysical Research: Oceans
where C is the vertically integrated ocean heat transport convergence within the layer of interest and Q net is the combined influence of radiative and turbulent fluxes of heat at the air-sea interface. Our estimates of heat transport convergence implicitly include the impacts of advection at all scales, including variations that occur in response to wind forcing from the atmosphere, as well as contributions from mixing and diffusive processes. Note that heat transport convergences integrated over the ocean mixed layer (C mld ) include a contribution from vertical advection at the lower boundary whereas all vertical advection terms cancel when heat transport convergences are integrated to the ocean floor (C tot ). A complete derivation of the relevant conservation equations is provided as supporting information.

Surface Heat Fluxes
Monthly air-sea heat fluxes for the period 1985-2012 are estimated using two complementary ''direct'' and ''residual'' approaches based on constraints from data-assimilating atmospheric reanalyses. For the direct approach, estimates of net air-sea heat flux are derived by combining monthly mean values of latent heat flux, sensible heat flux, net shortwave radiation and net longwave from the standard output of the NCEP/ NCAR [Kalnay et al., 1996] and ERA-interim [Dee et al., 2011] atmospheric reanalyses. These products provide globally complete fields that are constrained by the model physics and assimilation of available observations from the atmosphere and air-sea interface. However, direct flux estimates from reanalyses can exhibit significant biases due to poorly constrained near-surface air temperature or humidity and/or deficiencies in the representation of low-level clouds [Josey et al., 2014[Josey et al., , 2013.
To reduce the impact of systematic biases in the direct flux estimates, we also estimate fluxes using a residual method. In this approach, air-sea heat fluxes for the period 1985-2012 are estimated by combining satellite-derived radiative fluxes at the top of atmosphere [Allan et al., 2014] with estimates of atmospheric energy tendencies and divergences from the ERA-Interim [Dee et al., 2011] and MERRA [Rienecker et al., 2011] atmospheric reanalyses. Atmospheric energy divergences are mass corrected [Mayer and Haimberger, 2012] and the land surface fluxes are constrained to ensure realistic energy balance . Adjustments to land surface energy fluxes are redistributed over the oceans to produce realistic meridional transports of total energy . Net air-sea fluxes estimated using the residual method are less sensitive to biases in the near-surface meteorology and cloud properties, but are more affected by deficiencies in the atmospheric circulation [Josey et al., 2013;Liu et al., 2015;Trenberth et al., 2001]. Also note that flux estimates from the residual method give the net flux out of the atmospheric column, which in areas of sea-ice cover is not necessarily the same as the net flux into the ocean. For this reason, we do not include regions with climatological sea-ice cover >10% in our analysis.
Monthly mean anomalies from all four products, and their correlation with the ERA-interim direct flux estimate, are shown for an illustrative location in the east equatorial Pacific in Figure 2. This comparison demonstrates that all four products are able to resolve interannual variability of air-sea fluxes associated with ENSO, and also illustrates that the differences between residual and direct flux estimates from the same reanalysis are comparable to the differences between direct estimates from two different reanalysis systems.

Ocean Heat Content
Monthly ocean heat content anomalies for 1985-2012 are estimated using two gridded analyses of in situ temperature profiles from version 4.1.1 of the EN4 database of quality controlled ocean observations [Good et al., 2013]. The gridded analysis produced as part of the EN4 quality control procedure is based on a local optimal interpolation algorithm combined with a background damped persistence forecast using anomalies from the previous month [Good et al., 2013]. Time-varying expendable bathythermograph (XBT) biases are corrected using the scheme described by Gouretski and Reseghetti [2010]. The Met Office statistical ocean reanalysis (MOSORA)  uses an optimal interpolation based on global covariances informed by climate model simulations to interpolate sparse observations. In the MOSORA analysis, XBT biases are corrected using an average of the adjustments proposed by Wijffels et al.

Estimation of Heat Transport Convergences
Global ocean heat transport convergences (C) are estimated as a residual by combining observations of Q net , H, and their associated uncertainties using a forward Kalman filter and backward Rauch-Tung-Striebel smoother (RTS) [Rauch et al., 1965;Wunsch, 1996] based on local heat conservation. This method allows meaningful propagation of uncertainties and is similar to that described in Kelly et al. [2014Kelly et al. [ , 2016. Estimates of H, Q net , and C are constrained via a weighted least-squares optimization that minimizes the misfit between observations and model predictions with weights depending on their associated uncertainties. A full description of the methodology is provided as supplementary material. Henceforth, estimates derived using the Kalman RTS smoother are identified as b H; d Q net , and b C.
Prior to the application of the Kalman smoother, observational data sets are regridded to a common 18 3 18 grid using bilinear interpolation. Grid boxes with more than 10% climatological ice-cover and/or bathymetry shallower than 200 m are not included in our analysis. In order to isolate variability on interannual time scales, observation-based data sets are deseasonalized using the mean seasonal cycle for the period 2006-2012 and smoothed using a 12 month low-pass Butterworth filter. Note that only the seasonal variations are removed from the deseasonalized data such that their time-mean values are preserved.
In the final step before application of the Kalman filter, multiproduct means are calculated from the smoothed and deseasonalized estimates of Q and H and spatially varying observational uncertainties are estimated using where n p is the number of observational data sets (n H p 5 2 and n Q p 5 4), n t is the length of the time series, x p t is the estimated value for a specified month and data product, and x t is the multiproduct mean. Heat budgets are then performed independently for each grid box using a Kalman RTS smoother based on a centered-difference approximation to equation (2) (see supporting information).
To illustrate our methodology, we present here a brief comparison with a well-documented event in the North Atlantic subtropical gyre (NASTG). Observational studies have shown that the upper 2000 m of the NASTG between 268N and 418N cooled throughout 2010 and remained anomalously cool until the end of 2011 [Cunningham et al., 2013;Bryden et al., 2014]. This large-scale cooling of the ocean was attributed to a heat transport deficit at the southern boundary associated with interannual variability of the Atlantic meridional overturning circulation (  H tot during 2010 that can only be explained by negative heat transport convergence anomalies b C tot (Figure 3c). Based on this comparison, we conclude that our method is capable of resolving large-scale interannual variations in the ocean heat budget. Our estimates of b C tot during this period are positively but weakly correlated (r 5 0.46) with observed variations in meridional ocean heat transport from the RAPID-MOCHA array at 268N [Johns et al., 2011] (Figure 3c). However, we do not expect a perfect relationship as the RAPID transports do not account for changes at the northern boundary of the domain and there are uncertainties in our estimated heat transport convergences.

Ekman Forcing
In order to interpret the drivers of heat transport convergence variability, we follow a similar approach to previous studies [e.g., Gill and Niller, 1973;Marshall et al., 2001;Buckley et al., 2014] and calculate the local forcing from horizontal Ekman transports and vertical Ekman pumping. Note that this forcing accounts for variations in ocean heat content that occur in the absence of further adjustments in the ocean. For this reason, we do not expect these terms to be a good explanation for heat content variations in regions where baroclinic waves can adjust on time scales of less than a year such as along equatorial and coastal wave guides. Ekman heat transport convergences are separated into two contributions: (i) an upper ocean contribution due to Ekman layer transports across horizontal temperature gradients (C mld ek ) that we assume is constrained to the mixed layer and (ii) an ocean interior contribution due to isopycnal heave associated with Ekman pumping across vertical temperature gradients (C int ek ). These two contributions are calculated as follows: where T ek is the temperature of the Ekman layer, T is the depth-averaged temperature of the ocean interior, f the Coriolis parameter, s is wind stress, U ek 5 s3b z q o f is the depth-integrated Ekman layer transport, and w ek 5 1 q 0 curl s f À Á is the vertical Ekman pumping velocity. We then combine these terms into a single total Ekman forcing term C tot ek 5C mld ek 1C mld ek . Note that since we are only interested in the component of C ek that can be attributed to local atmospheric variability, and not the changes due to interannual variations in temperature that we cannot directly attribute to the atmosphere, we calculate C ek using a seasonally varying temperature climatology. This approach is justified because Buckley et al. [2015] found that monthly variations of C ek in most regions of the North Atlantic are dominated by local wind variability rather than changes to the temperature field. A full derivation of equations (4) and (5) is provided as supporting information.

Diagnosing Drivers of Ocean Variability
We evaluate drivers of ocean variability using Pearson correlation coefficients between the different terms of the ocean heat budget. A process is diagnosed as a driver of ocean variability if it has a significant positive correlation with the ocean heat content tendency. In addition, we use a fraction of variance (FOV) skill score that accounts for the magnitude of explained variability that is defined as follows:

FOV512
VarðX2YÞ VarðXÞ (6) where X is the signal we are trying to explain and Y is the variability attributable to a specific process. A perfect score (FOV 5 1) occurs when X and Y are perfectly correlated and have the same variance. Since the magnitude of the signal is also taken into account, it is possible for a process to have a FOV skill score < 0 even it is otherwise well correlated with the target variable. This scenario is indicative of compensation between processes such that the magnitude of the residual variability is small compared to the magnitude of the compensating processes.

Seasonal Cycle
To provide context for our analysis of interannual variability, we first evaluate the relative importance of the different terms in the seasonal heat budget. Seasonal variations in C mld and C tot are estimated as a residual Journal of Geophysical Research: Oceans 10.1002/2016JC012278 from equation (2) following estimation of seasonal variations in heat storage rates. Seasonal cycle amplitudes and correlations between each term of the seasonal heat budget are shown in Figures 1 and 4 and the FOV in seasonal heat storage rates explained by Q net and C is shown in Figure 5.
An initial inspection of seasonal cycle amplitudes indicates that seasonal changes in @ @t H mld (Figure 1c) are largely a response to heat input at the ocean surface. Seasonal correlations and FOV skill scores (Figures 4a, 4c and 5a, 5c) confirm this general result, but also highlight the contribution from seasonal variations in C mld along the equator and in regions with large meridional sea-surface temperature (SST) gradients (e.g., in the Antarctic circumpolar current, ACC). This conclusion is consistent with an analysis of the Atlantic Ocean from Piecuch and Ponte [2012] who concluded that seasonal variations generally represent surface inputs, except near the equator, where advective changes are important.
In contrast, seasonal variations in @ @t H tot cannot be explained without substantial variations in C tot , particularly within 6208 of the equator and in the Southern Ocean (Figures 1d, 1f, 4b, 4d, and 5b, 5d). We interpret this result as evidence for the importance of seasonal changes in heat content below the mixed layer associated with the adiabatic heave of isotherms caused by seasonal variations in upwelling and the wind-driven circulation.

Mixed Layer Ocean Heat Budgets
To evaluate the magnitude of interannual variability (i.e., the signal, S) in each component of our mixed layer heat budget, we calculate variances for our smoothly varying monthly estimates of b Q net ; @ @t b H mld , and b C mld (Figures 6a, 6b, and 6d). In the North Atlantic region, our variance estimates can be compared with those from the ECCO ocean state estimate used by Buckley et al. [2014]. In general, our variances are lower since we are looking at interannual time scales rather than monthly means. However, there are some significant differences in the patterns of variability that cannot be attributed to temporal filtering. In particular, Buckley et al. [2014] find large magnitude variations in the northernmost subpolar gyre and coastal Labrador Sea that is not present in the raw monthly data used for the present study (not shown). Possible reasons for . Correlations between climatological seasonal cycles of (a) surface heat flux (Q net ) and mixed layer depth heat content tendency @ @t H mld À Á , (b) Q net and total ocean heat content tendency @ @t H tot À Á , (c) mixed layer heat transport convergence (C mld ) and @ @t H mld , and (d) total ocean heat transport convergence (C tot ) and @ @t H tot . Stippled areas indicate regions with jrj ! 0.576, the critical value corresponding to p 5 0.05 for a two-sided test assuming N22 degrees of freedom, where N 5 12.
Journal of Geophysical Research: Oceans 10.1002/2016JC012278 such a discrepancy include: (i) a systematic bias in the way that reanalysis products represent air-sea exchanges over the boundary currents in the subpolar gyre or (ii) imperfect adjustments to the surface boundary conditions of the ECCO product as part of the optimization procedure (see Figure 2 in Buckley et al. [2015]).
We also evaluate signal/noise (S/N) ratios for b Q net and b C mld , where N5 1 nt P nt t51 r t and r 2 t is the time-varying Kalman smoother prediction error variance. S/N ratios are not reported for @ @t b H mld or @ @t b H tot as our Kalman smoother methodology returns uncertainty estimates for heat content rather than heat content tendencies. This calculation shows that S/N > 1 in almost all locations (stippled regions in Figures 6a and 6d) indicating that our estimates of b Q net and b C mld are robust given our estimated uncertainties.
From Figure 6, it is evident that the magnitude and spatial patterns of b Q net and b C mld variability are very similar, with the largest magnitude signals in equatorial and frontal regions. The exception is the Southern Ocean, where signals in b Q net are generally smaller than those in b C mld indicating an important role for ocean dynamics in the mixed layer of the eddy-rich ACC. Correlations between mixed layer heat budget terms (Figures 7a and 7c) and FOV skill scores (Figures 8a and 8c) reveal that the relative importance of b Q net and b C mld varies strongly by region, as found by previous studies of the North Atlantic [Buckley et al., 2014].

Variations in b
Q net are an important contributor to b H mld variability in the interior of subtropical gyres, in the North Atlantic and North Pacific subpolar gyres, and at the northern limit of Antarctic sea ice extent. The pattern of b H mld variability explained by b Q net in the North Atlantic (Figure 8a) is very similar to that found by Buckley et al. [2014] in an analysis of monthly mean anomalies, although the magnitude of FOV skill is substantially lower in the present study. This could indicate an increased role for ocean dynamics on interannual time scales, or it could be a consequence of the inherent uncertainty associated with observational heat budgets.

Variations in b
C mld are particularly dominant in the central equatorial Pacific and the Antarctic circumpolar current (Figure 8c). However, there are some regions where we expect ocean dynamics to be important  Figure 5a mixed layer depth heat content tendency @ @t H mld À Á explained by surface heat flux (Q net ), (b) total ocean heat content tendency @ @t H tot À Á explained by Q net , (c) @ @t H mld explained by mixed layer heat transport convergence (C mld ), and (d) @ @t H tot explained by total ocean heat transport convergence (C tot ). Locations with FOV < 0 are masked white.

10.1002/2016JC012278
that have FOV skill scores below zero, such as the western boundary currents in the North Atlantic and North Pacific and the eastern tropical Pacific. In these locations, there is strong compensation between b Q net and b C mld (Figure 9) due to damping of mixed layer temperature anomalies by air-sea heat fluxes. Thus, the residual signal in @ @t b H mld is not well explained by b C mld . There is also some evidence of the opposite process in the subtropical Pacific, where variations in mixed layer temperature driven by b Q net are damped by b C mld . Finally, in the eastern North Atlantic and eastern North Pacific, b Q net and b C mld both play an important role in driving heat content changes in the mixed layer.
In Figure 10a, we have schematically summarized the drivers of b H mld by identifying the regions where correlations are significant in Figure 7. Consistent with previous studies [Jin, 1997;Dong et al., 2007;Dong and Kelly, 2004;Buckley et al., 2014], we find that ocean heat transport processes dominate the upper ocean heat budget in the equatorial oceans and in regions of strong currents and substantial eddy activity. In addition, there are large regions of the Atlantic and Pacific oceans where ocean dynamics are an important contributor, along with local air-sea fluxes, to the generation of mixed layer temperature anomalies. Based on this synthesis, we conclude that it is inaccurate to consider interannual variations of b H mld to be a passive  response to local atmospheric forcings [e.g., Frankignoul and Hasselmann, 1977;Frankignoul, 1985], except in selected parts of the high latitude and subtropical oceans where Q net is dominant (Figure 10a). This conclusion suggests there that there is potential for interannual predictability of SSTs in the midlatitudes originating from active ocean dynamics and initialization of the ocean state that goes beyond the previously described one-dimensional persistence mechanisms [e.g., Alexander and Deser, 1995;Cassou et al., 2007;Taws et al., 2011].

Surface Ocean Feedbacks on the Atmosphere
To evaluate the role of ocean feedbacks on the atmosphere on interannual time scales, we calculate correlations between b H mld and b Q net (Figure 11a). We also calculate correlations between b H mld and the turbulent flux contributions to b Q net using latent and sensible heat fluxes from the ERA-interim and NCEP/NCAR reanalyses (Figures 11b and 11c). This comparison demonstrates that the response of Q net to changes in b H mld is dominated by changes in the turbulent fluxes. Figures 10a and 11a reveals that there is a good correspondence between the zonally coherent regions of negative correlation between b Q net and b H mld and those regions where ocean dynamics is identified as the dominant contributor to the mixed layer heat budget such as the equatorial oceans, western boundary currents, and in the ACC (Figure 10a). Thus, in these locations, our analysis indicates that interannual variations in b C mld are responsible for the generation of b H mld anomalies, which then generate turbulent heat flux anomalies that drive the overlying atmosphere. Importantly, surface flux variations in these  Similarly, there is also a good correspondence between the approximately zonal regions of positive correlation between b Q net and b H mld (Figure 11a) and locations where local air-sea fluxes are determined to be important in the mixed layer heat budget, such as the subtropical Pacific and subpolar oceans (Figure 10a). In these locations, interannual b H mld anomalies are generated by turbulent flux anomalies forced by near- surface meteorological conditions such that negative surface temperature anomalies occur when heat is fluxed from the ocean into the atmosphere.

Full-Depth Ocean Heat Budgets
Following the same procedure as for the mixed layer, we now evaluate the dominant terms in the interannual heat budget for the full-depth ocean. Similar to our results for the seasonal heat budget, we find that variances of @ @t b H tot and b C tot are substantially larger than those in b Q net (Figures 6a, 6c, and 6e). In addition, variability in @ @t b H tot is substantially larger than in @ @t b H mld , which emphasizes the importance of heat content changes below the mixed layer on interannual time scales. These observations are consistent with the expectation that, due to its vertically integrated nature, contributions from heat transport convergences will become more important when heat budgets are performed for ocean layers of increasing depth. Correlations (Figures 7b and 7d) and FOV skill scores (Figures 8b and 8d) confirm that variations in b C tot dominate the full-depth heat budget on interannual time scales, except in selected locations in the high latitudes such as the Labrador Sea where b Q net explains more than half of the variance in @ @t b H tot . The near-global dominance of ocean dynamical processes for the generation of interannual variability in b H tot is also clear in our schematic summary plot (Figure 10b).
Our results emphasize that, away from regions of deep convection and water mass transformation, interannual variations in full-depth ocean heat content (and thus thermosteric sea level) are controlled by the rearrangement of heat within the ocean rather than the loss or addition of heat at the surface. This result is in agreement with previous studies of full-depth regional heat budgets [e.g., Grist et al., 2010;Cunningham et al., 2013;Bryden et al., 2014;Kelly et al., 2014] and consistent with work that has found interannual variability of thermosteric sea level to be governed by the adiabatic heave of isotherms in response to surface wind-forcing [K€ ohl, 2014;Forget and Ponte, 2015;Stammer et al., 2013;Roberts et al., 2016]. Importantly, this analysis demonstrates that the term-balance in a full-depth ocean heat budget is not always a useful predictor of the important balances in the mixed layer. In addition, our analysis suggests that variations in thermosteric sea level may be substantially more predictable than near-surface temperature anomalies due to the greater importance of ocean dynamics and therefore increased potential for skill arising from the initialization of the ocean state.

Impact of Local Ekman Forcing
Having identified the regions where heat transport convergences are important for interannual variations in b H mld and b H tot , we now evaluate the extent to which these transports represent an Ekman response to local wind-stress curl forcing. Using the expressions defined in section 2.5, we calculate estimates of C mld ek and C tot ek using monthly wind stress data from the ERA-interim reanalysis [Dee et al., 2011]. These data are then smoothed with a 12 month Butterworth filter for comparison with b C mld and b C tot , respectively. Note that we do not calculate Ekman forcing terms within 658 of the equator as f ! 0.
Correlation analysis (Figures 12a and 12b) indicates that there are some regions where Ekman forcings are an important contributor to interannual variability of b C mld and b C tot , particularly in the subtropical oceans and in the eastern North Pacific. However, the FOV explained by locally forced Ekman variations (Figures  12c and 12d) is typically much less than 50%, especially in areas where b C mld dominates the mixed layer heat budget (Figure 10a). This result implies that interannual variations in heat transport convergences must be dominated by non-Ekman dynamics such as baroclinic wave adjustments to remote atmospheric forcings, intrinsic ocean variability associated with eddies, and the advection of temperature anomalies by the mean flow. Importantly, these are processes that offer the potential for interannual predictability through the initialization of the ocean state, whereas locally forced Ekman transports require predictability of atmospheric wind stress curl anomalies.
This result contradicts with the analysis of Buckley et al. [2014], which found a more significant role for Ekman heat transport convergences in the generation of heat content anomalies in the interior of the North Atlantic ocean. The origin of this discrepancy appears to be a significant difference in the magnitude of Ekman heat transport convergence variability derived from ERA-interim winds compared to estimates from the ECCO ocean state estimate, particularly over the North Atlantic subpolar gyre. Our estimates of Ekman driven heat transport convergence are consistent with previous studies of the North Atlantic [Marshall et al., 2001] and are not particularly sensitive to our choice of atmospheric wind data set. It is possible that these Journal of Geophysical Research: Oceans ek are a result of the adjustments applied to wind stress boundary conditions as part of the ECCO optimization procedure.

Sensitivity to an Evolving Observing System
During the period 1985-2012 there were important changes in the observing systems that constrain estimates of ocean heat content and surface heat fluxes. In the early 2000s, the Argo network of autonomous profiling floats [Riser et al., 2016] replaced XBTs as the dominant source of temperature observations resulting in a corresponding increase in the depth of regular temperature observations from less than 700 m to about 2000 m [Lyman and Johnson, 2014]. In addition, from 2000 onwards, estimates of top-of-atmosphere radiation fluxes have been constrained by the Clouds and the Earth's Radiant Energy System (CERES) satellite observations [Loeb et al., 2012]. To evaluate the potential impact of changes to the global observing system, we rerun our Kalman filter using only data from the Argo period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) and repeat the correlation analysis presented in Figure 10. Note that repeating our Kalman RTS analysis on this reduced data set is necessary to prevent the propagation of information from the earlier period. The updated synthesis figure for the Argo period ( Figure 13) shows the same features that are present in our earlier analysis,  although the boundaries between regions are less well-defined due to reduced length of the time series used to calculate correlations. Based on this comparison, we believe our main conclusions are unaffected by systematic changes in the observing system since 1985.
Despite the near global coverage of the upper 2000 m of the ocean provided by Argo, temperature observations from the abyssal ocean are still limited to those recovered from geographically sparse hydrographic sections that are repeated every 1-10 years [Purkey and Johnson, 2010]. This dearth of observations from depths >2000 m means that our estimates of interannual variability in H tot may be systematically biased. However, until the planned extension of the Argo array to the abyssal ocean takes place [Johnson et al., 2015] we cannot accurately quantify the magnitude of this missing signal from observations alone. In addition, data-assimilating model syntheses show a wide range of behaviors in the deep ocean where they are unconstrained by observations . However, although the infrequent observations from hydrographic sections are insufficient to characterize interannual variability in the deep ocean, they can be used to estimate the magnitude of decadal trends in ocean basins Johnson, 2010, Desbruyeres et al., 2016]. Figure 14 shows ocean heat content trends from the EN4 analysis for three different depth layers compared with an analysis of repeat hydrographic sections [Desbruyeres et al., 2016] for the period 2000-2015. From this comparison is evident that the EN4 analysis fails to capture both the pattern and magnitude of abyssal heat content trends. However, the regional trends in the abyssal ocean are generally an order of magnitude smaller than those in the upper 2000 m. Based on this comparison, we expect interannual variability, like decadal variability, to be intensified within the upper 2000 m of the ocean. Therefore, although there are deficiencies in our ability to observe interannual variability in the deep ocean, our assessment is that these missing signals are unlikely to have a substantive impact on our results ( Figures  10 and 13).

Summary and Conclusions
In this study, we have presented an observation-based analysis of seasonal and interannual ocean heat budgets and evaluated the dominant drivers of H mld and H tot for the period 1985-2012. We combined four surface heat flux products and two ocean heat content products using a novel Kalman smoother-based method and inferred the regional contributions from ocean heat transport convergences as a residual. In addition, we quantified the locally forced contribution to ocean heat transport convergence variability due to Ekman transports and identified regions where the ocean is driving the atmosphere on interannual time scales.
For the climatological annual cycle, variations in H mld are dominated by seasonal variations in surface heat input, except along the equator and in frontal regions where heat transport convergences are important. In contrast, seasonal variations in H tot (and thus thermosteric sea level) cannot be explained without substantial heat transport contributions, particularly within 6208 of the equator and in the Southern Ocean. This result emphasizes the importance of seasonal changes in heat content below the mixed layer associated with the adiabatic heave of isotherms caused by seasonal variations in upwelling and the wind-driven circulation.
On interannual time scales, non-Ekman ocean heat transport processes are found to dominate H mld variations in the equatorial oceans and regions of strong ocean currents and substantial eddy activity. In these locations, surface temperature anomalies generated by ocean dynamics result in turbulent heat flux anomalies that drive the overlying atmosphere. In addition, both heat transport convergences and local air-sea fluxes are important for the generation of mixed layer temperature anomalies in large regions of the Atlantic and Pacific oceans. Based on these results, we concluded that ocean dynamics play an active role in the generation of H mld anomalies on interannual time scales, except in selected regions of the high latitudes and subtropical oceans where forcing from local air-sea fluxes dominates. Away from regions of deep convection and water mass transformation, interannual variations in full-depth ocean heat content (and thus thermosteric sea level) are controlled by the rearrangement of heat within the ocean rather than the loss or addition of heat at the surface. This result demonstrates that the term-balance in a full-depth ocean heat budget is not always a useful predictor of the important balances in the mixed layer.
These results have important implications for the predictability of H mld and H tot on interannual time scales. The importance of active ocean dynamics for the generation of mixed layer temperature anomalies across large regions of the extratropical oceans suggests that there is potential for interannual predictability of SSTs in the midlatitudes originating from initialization of the ocean state that goes beyond the previously described one-dimensional persistence mechanisms [e.g., Alexander and Deser, 1995;Cassou et al., 2007;Taws et al., 2011]. In addition, our analysis suggests that variations in thermosteric sea level may be substantially more predictable than near-surface temperature anomalies due to the greater importance of ocean dynamics and therefore increased potential for skill arising from the initialization of the ocean state. Finally, given the importance of the mixed layer budget for surface temperature anomalies, the budget-perspective we have presented could also be instructive for understanding and attributing changes in heat content associated with recently described global warming ''hiatus' ' and ''surge'' events [e.g., Roberts et al., 2015;Meehl et al., 2011;England et al., 2014;Kosaka and Xie, 2013].