The Signature of Oceanic Processes in Decadal Extratropical SST Anomalies

The relationship between decadal sea surface temperature (SST) and turbulent heat fluxes is assessed and used to identify where oceanic processes play an important role in extratropical decadal SST variability. In observational data sets and coupled climate model simulations from the Coupled Model Intercomparison Project Phase 5 archive, positive correlations between upward turbulent heat flux and SSTs indicate an active role of oceanic processes over regions in the North Atlantic, Northwest Pacific, Southern Pacific, and Southern Atlantic. The contrasting nature of oceanic influence on decadal SST anomalies in the Northwest Pacific and North Atlantic is identified. Over the Northwest Pacific, SST anomalies are consistent with changes in the horizontal wind‐driven gyre circulation on timescales of between 3 and 7 years, in both the observations and models. Over the North Atlantic, SST anomalies are also preceded by atmospheric circulation anomalies, though the response is stronger at longer timescales—peaking at around 20 years in the observations and at around 10 years in the models.


Introduction
The variability of sea surface temperatures (SSTs) has been linked with regional climate anomalies on decadal timescales. A prominent example is the Atlantic Multidecadal Oscillation (AMO), or Atlantic Multidecadal Variability, in the North Atlantic (Enfield et al., 2001), which exhibits a substantial influence on the climate over Europe (O'Reilly et al., 2017;Sutton & Dong, 2012;Sutton & Hodson, 2005), the United States (Nigam et al., 2011;McCabe et al., 2004;Zhang & Delworth, 2006), and the Sahel (Folland et al., 1986;Knight et al., 2006). Variability in the ocean circulation has been cited as an important factor controlling the AMO. Many modeling studies have suggested that variability in the Atlantic meridional overturning circulation plays an important role in governing the multidecadal SST variability in the North Atlantic (Delworth et al., 1993;Delworth & Mann, 2000;Drews & Greatbatch, 2017;Latif et al., 2004;Knight et al., 2005), while other studies have indicated that much of the observed SST variability in the subpolar North Atlantic is associated with changes in the horizontal gyre circulation (Häkkinen & Rhines, 2009;Häkkinen et al., 2011;Piecuch et al., 2017).
The study by Clement et al. (2015), however, argued that ocean circulation is not important for generating the AMO. Clement et al. showed that climate models coupled to a slab ocean (i.e., without varying ocean circulation) were able to generate multidecadal SST variability that resembled the observed AMO. However, on 10.1029/2018GL079077 decadal timescales there are positive (upward) turbulent heat flux anomalies in the warm phase of the AMO (and vice versa in the cool phase of the AMO) in both observations (Gulev et al., 2013) and fully coupled climate models with varying ocean circulation (O'Reilly et al., 2016), whereas turbulent heat fluxes do not significantly covary with North Atlantic SSTs in the climate models with a slab ocean (O'Reilly et al., 2016;Zhang et al., 2016). This indicates that the processes, which participate in the coupled evolution of SSTs and turbulent heat flux, rely on variability of ocean processes, such as variability in horizontal/overturning ocean circulation or horizontal/vertical mixing processes. However, the nature and timescales of these oceanic processes remain unclear.
Recently, Cane et al. (2017) used a simple stochastic model to demonstrate that uncorrelated white noise forcing from ocean processes-in addition to atmospheric white noise forcing-is sufficient to generate the coupled evolution of decadal SST and turbulent heat flux anomalies seen in observations and fully coupled models. This result is important as it indicates that although ocean processes are necessary for SST and turbulent heat flux anomalies to be correlated on decadal timescales, it is not necessary for the ocean processes themselves to exhibit coherent decadal forcing. Cane et al. showed that the low-pass correlation (r) of SST and turbulent heat flux is directly related to the fraction of the total forcing from oceanic processes (b 2 ), through the relation r 2 = b 2 . However, Zhang (2017) argued, through comparison of SST and salinity anomalies, that the ocean's role is somewhat underestimated in the study of Cane et al. (2017), due to the lack of oceanic damping processes. Nonetheless, the simple relation of Cane et al. (i.e., r 2 = b 2 ) shows that oceanic processes play a role in determining the nature of decadal SST anomalies where they covary significantly with turbulent heat fluxes (Zhang, 2017). Moreover, the magnitude of the correlation gives an approximate measure of the strength of the oceanic contribution in this framework. In addition to ocean processes, external forcing has been shown to explain some of the AMO variance over the observational period (Bellomo et al., 2018;Booth et al., 2012;Mann & Emanuel, 2006;Murphy et al., 2017;Ting et al., 2014).
While much work on the ocean's role in governing decadal SST anomalies in the extratropics has focused on the North Atlantic, decadal SST variability has also been extensively documented in North Pacific (Mantua et al., 1997;Minobe, 1997;Newman et al., 2016Newman et al., , 2003. Decadal SST variability has also been documented in the Southern Ocean (Gille, 2008;G. Wang & Dommenget, 2016), although to a lesser extent owing to the sparsity of observations over the twentieth century. In this paper we analyze the relationship between decadal SSTs and turbulent heat fluxes across the global oceans. We use the decadal SST and turbulent heat flux correlation, based on Cane et al. (2017), as a diagnostic to identify where oceanic processes are important for decadal SST variability. The variability of SST in these regions is analyzed and compared with the North Atlantic, using observational data sets and coupled climate models. Our results highlight the contrasting nature of oceanic control on maintaining decadal SST anomalies in the North Pacific and North Atlantic-with a particular emphasis on the role of large-scale atmospheric forcing and the characteristic timescales of the oceanic response.

Observation-Based Gridded Data Sets
We analyze observational data sets between 1880 and 2014. We use monthly SSTs from the National Oceanic and Atmospheric Administration Extended Reconstructed SST version 4 data set (Huang et al., 2015). We use monthly turbulent heat flux (i.e., sensible plus latent heat flux) and sea level pressure (SLP) data from the 20th Century Reanalysis (20CR) Project v2c (Compo et al., 2011). The 20CR is a reanalysis product produced by assimilating only surface pressure observations with a prescribed SST boundary condition, such that the turbulent heat fluxes are calculated by the model and should be interpreted cautiously. However, Gulev et al. (2013) showed that the relationship between decadal turbulent heat fluxes and SSTs in the North Atlantic are qualitatively similar to the turbulent heat fluxes calculated only from ship observations (although the latter are reconstructed from sparse observations in some periods). Therefore, the 20CR provides one plausible estimate of the turbulent heat fluxes over the observational period. All products were linearly detrended and regridded to a common 5 ∘ × 5 ∘ grid prior to the analysis to focus on relatively large-scale SST anomalies. We tested the sensitivity to the linear detrending method by quadratically detrending and also removing a warming trend based on historical CO 2 concentrations (Meinshausen et al., 2017), neither of which qualitatively changed the results. More quantitatively, the correlation values and anomalies in the key regions identified below all differ by less than 10% with different detrending methods.

Preindustrial Control and Historical Simulations
In addition to the observational products, we analyze output from the Coupled Model Intercomparison Project Phase 5 (CMIP5) archive (Taylor et al., 2012). The first set of experiments consist of the first member of each of the preindustrial control (piControl) simulations in the CMIP5 archive. Monthly turbulent heat flux, SST, and SLP data were used for the 38 models in the archive that had at least 200 years of preindustrial simulation output. The preindustrial control simulation has no changes in external forcing (i.e., greenhouse gas emissions, anthropogenic aerosols, and volcanic eruptions); in these simulations we are therefore analyzing the internal variability of the coupled models. The second set of experiments consist of the first member of each of the historical simulations in the CMIP5 archive. Monthly turbulent heat flux, SST, and SLP data were used for the 35 models (a subset of the 38 preindustrial models) in the archive that have historical simulation output that covers the period 1861-2005. The historical simulations include estimates of historical greenhouse gas concentrations, anthropogenic aerosols, and volcanic eruptions, and studies have argued that-through radiative effects-these external forcings may be responsible for significant decadal SST variability (e.g., Bellomo et al., 2018;Booth et al., 2012). As with the observational data sets, the data from the historical simulations were linearly detrended prior to the analysis (though the results are not sensitive to this particular detrending method). All model variables were regridded to a common 5 ∘ × 5 ∘ grid prior to the analysis to simplify the analysis across the models and to aid comparison with the observational data sets.
Results from the preindustrial control and the detrended historical simulations are found to be qualitatively very similar. The main difference between the two ensembles is that the agreement between the models is larger in the preindustrial control models. This is likely due to the greater sample size in each of the preindustrial control models, which have an average length of 512 years, compared to the data sets from the historical simulations, which are all 145 years. Both CMIP5 ensembles are discussed throughout but the preindustrial simulations are primarily presented in some figures, with the corresponding plots for the CMIP5 historical simulation included in the supporting information.

Significance Tests
The significance of the correlation and regressions calculated from the observational data sets is estimated using a Monte Carlo procedure. For each calculation, 1,000 surrogate time series were generated by taking the Fourier transform and randomizing the phase of each component before taking the inverse transform, such that the resulting time series replicate the spectral properties of the original (Ebisuzaki, 1997;Kaplan & Glass, 2012). The resulting distribution of correlation or regression coefficients was then used to estimate the probability that the observed magnitude could occur by chance.

Signature of Oceanic Processes in Regional Decadal SST Anomalies
To assess where oceanic processes appear important in influencing decadal SST anomalies, we first analyze the correlation between the decadal SST and turbulent heat flux anomalies (defined as positive out of the ocean).
Here and in the analysis that follows we use a simple 10-year moving average to define the decadal anomalies. Figure 1 shows maps of these grid point correlations for the observational data sets and the ensemble means of the CMIP5 preindustrial control simulations and the CMIP5 historical simulations. The correlation between the observational SST and turbulent heat fluxes is large (r > 0.6) over the North Atlantic. Despite the use of a different heat flux data set, this result is similar to that of Gulev et al. (2013), albeit peaking slightly further east. The CMIP5 models, which are very similar for the preindustrial and the historical simulations, also show positive correlations between decadal SST and turbulent heat flux anomalies over the North Atlantic, although the region of positive correlations extends further poleward than in the observational data sets. Therefore, in both the observations and models, the nature of decadal SST and turbulent heat flux variability indicates a dependence on active ocean processes, as previously highlighted by Cane et al. (2017), O'Reilly et al. (2016, and Zhang et al. (2016). Many studies have argued that in observations, the SSTs are also likely influenced by external forcing (e.g., Bellomo et al., 2018;Booth et al., 2012;Mann & Emanuel, 2006;Swingedouw et al., 2017;Ting et al., 2014). However, the relationship between SSTs and turbulent heat flux anomalies in the historical simulations are very similar to the pre industrial control simulations, suggesting that external forcing does not play an important role in determining this relationship, though it could impact the amplitude of the variability (e.g., Murphy et al., 2017;Swingedouw et al., 2017).
The most striking aspect of Figure 1, however, is that there are numerous regions across the global oceans where decadal SSTs and turbulent heat fluxes are positively correlated. While we emphasize that correlation is not causation, we will highlight possible mechanisms for the regional differences in the remaining parts of the paper. Relying on either the observational data sets and CMIP5 models alone might be unwise, due to observational sparsity (e.g., Deser et al., 2010) and various model biases (e.g., C. Wang et al., 2014). Therefore, we will only focus on regions where decadal SSTs and turbulent heat fluxes are positively correlated in both the observational data sets and CMIP5 model ensembles. One such region is the tropical Pacific, where the positive correlations indicate that ocean processes are important in the evolution of decadal SST and turbulent heat flux anomalies. This is likely associated with the decadal variability of the El Niño/Southern Oscillation, with El Niño (La Niña) periods characterized by warm (cold) SST anomalies and large turbulent heat flux anomalies out of (in) the tropical Pacific (e.g., C. Wang, Deser, et al., 2017), induced by the thermocline adjustment via Kelvin waves in the tropical Pacific (e.g., McPhaden, 1999).
Perhaps more interesting though are the regions in the extratropics where decadal SSTs and turbulent heat fluxes are positively correlated. Along with the North Atlantic, there are also broadly consistent signatures of active ocean processes in the Northwest Pacific, the South Pacific, and the South Atlantic, highlighted by the black boxes in Figure 1. The region in the Northwest Pacific is within the Kuroshio-Oyashio Extension region, where the confluent western boundary currents detach from Japan (Qiu, 2001). Similar covariability between decadal SST and turbulent heat flux anomalies in this region using observation-based data sets have been identified (e.g., Tanimoto et al., 2003) with the turbulent heat flux out of the ocean associated with ocean heat flux convergence (Qiu et al., 2017). Finally, the regions of positive correlation in the South Pacific and South Atlantic are both located on the equatorward side of the Antarctic Circumpolar Current and the surface westerlies, where there is wind-driven downwelling associated with the Southern Ocean overturning circulation (Marshall & Speer, 2012).

The Role of Decadal Large-Scale Extratropical Atmospheric Circulation Anomalies
The decadal SST anomalies in the extratropical regions highlighted in Figure 1 are clearly not directly forced by turbulent heat fluxes associated with anomalous large-scale atmospheric circulation. However, they may be linked to ocean processes (advection or mixing) driven by wind and/or buoyancy forcing that are associated with large-scale atmospheric circulation anomalies. To initially gauge the role of large-scale atmospheric circulation anomalies we defined decadal-averaged SST indices computed over the regions highlighted in Figure 1 (the SST anomalies associated with these indices are shown in Figure S1). We then regressed the annual SLP and turbulent heat flux anomalies onto the decadal SST anomalies using both the observational data sets and the CMIP5 models, which are shown in Figure 2 (n.b. equivalent results are found if the SLP and turbulent heat flux and also decadally averaged prior to the regression calculations). Due to the relatively small sample size, the observational maps are much noisier than the CMIP5 model ensembles, therefore the color scales and contour intervals for the panels related to observations in Figure 2 are chosen to be double those for the CMIP5 models. Nonetheless, the SLP and turbulent heat flux anomalies associated with the decadal SST indices in the observations are broadly consistent with the corresponding anomalies in the CMIP5 models.
Warm decadal SST anomalies in the Northwest Pacific are associated with an anticyclonic SLP anomaly over the extratropical North Pacific and a narrow band of strong positive turbulent heat flux anomalies in both the observations (Figure 2a) and CMIP5 models (Figures 2e and S2e). The anticylonic SLP anomaly corresponds to a poleward shift in the climatological westerlies. The Northwest Pacific SST response to large-scale atmospheric circulation anomalies has been noted in observations (e.g., Schneider & Miller, 2001). The latitudes of the Kuroshio and Oyashio Extensions both exhibit a northward shift (not necessarily coherent) in response to the poleward shift of the westerlies (and the zero wind stress curl line) a few years prior (Nonaka et al., 2006;Sasaki et al., 2013;Schneider & Miller, 2001;Seager et al., 2001;Qiu et al., 2017). These shifts are due to wind-driven changes in the horizontal gyre circulation in the extratropical North Pacific. In the CMIP5 models, the Kuroshio and Oyashio Extensions are not well resolved and there is a single boundary current that detaches from Japan at 40-45 ∘ N. This potentially explains the a band of negative turbulent heat flux anomalies to the south of the SST index region in the CMIP5 models, which likely reflects a shift in the region of maximum turbulent heat flux along the Kuroshio-Oyashio Extension larger than seen in observations, which is only 2-3 ∘ in latitude (e.g., O'Reilly & Czaja, 2015;Qiu & Chen, 2005). The SST and turbulent heat flux anomalies are therefore consistent with a wind-driven northward shift of the Kuroshio-Oyashio Extension-associated with a poleward shift of the westerlies in the CMIP5 models-similar to that seen in previous individual model studies (e.g., Kwon & Deser, 2007;Schneider et al., 2002;Seager et al., 2001). Therefore, Figures 2a, 2e, and S2e indicate that the Northwest Pacific decadal SST anomalies in the observations and CMIP5 models are consistent with being generated by changes in the horizontal wind-driven ocean circulation (though here this is only inferred and not demonstrated explicitly).
Unlike over the Northwest Pacific, the decadal SST anomalies over the North Atlantic do not seem to be associated with significant large-scale atmospheric circulation anomalies in the observations (Figure 2b) or in CMIP5 models (Figures 2f and S2f ). This indicates that the decadal SST signals are not generated by large-scale atmospheric circulation anomalies on subdecadal timescales, otherwise we would expect substantial SLP anomalies in the regression maps, as in the Northwest Pacific. It is also worth noting that in these regression maps the positive turbulent heat fluxes that are significant seem to be located in the mid-Atlantic in the observations but in the subpolar gyre region in the CMIP5 models. These are only contemporaneous regression maps on decadal timescales and do not offer clear insight into whether the atmosphere or ocean are leading. We will analyze the response timescales in the next subsection.
Over the South Pacific region, decadal SST anomalies are associated with a deepening of the Amundsen Sea Low (e.g., Raphael et al., 2016) and a poleward shift of the climatological westerlies-in both the observations and the CMIP5 models. Directly over the index region, there is a weakening of westerlies (i.e., of the climatological SLP gradient). Over the South Atlantic region, there is also a consistent weakening of the local westerlies; however, the hemispheric SLP anomalies have the opposite sign over the polar region in observations and the CMIP5 models. This inconsistency is perhaps not surprising, given the paucity of Southern Hemisphere observations and the significant inadequacies in the simulations of Southern Ocean climate in the CMIP5 models (e.g., Sallée et al., 2013). Nonetheless, the weakening of the westerlies over both the South Pacific and South Atlantic regions indicates that the evolution of the decadal SST, and covarying turbulent heat flux anomalies are linked to local wind forcing.
A dominant term in the Southern Ocean mixed-layer heat budget is the horizontal Ekman heat flux convergence term, that is, dT dt ≈ −u Ek ⋅ ∇ H T, where T is mixed-layer temperature, ∇ H is the horizontal gradient operator, and u Ek = (U Ek , V Ek ) is the Ekman transport (S. Dong et al., 2007;Tamsitt et al., 2016). Specifically, this is dominated by the meridional advection, which is proportional to the zonal wind stress, x (i.e., V Ek = x f , where is density and f the Coriolis frequency). In the climatological mean, the meridional Ekman advection transports cold water equatorward and acts to cool the SSTs in the southern midlatitudes. In both the South Pacific and South Atlantic SST regions, the weakening of the local westerlies reduces the meridional equatorward Ekman transport and thereby reduces the cooling of the midlatitude SSTs, which results in warm decadal SST anomalies in these regions along with positive turbulent heat flux anomalies.

Dominant Timescales of the Atmospheric Forcing and SST Anomalies in Regions of Active Ocean Processes
The contemporaneous regression maps in Figure 2 indicate that the SST anomalies over the North Atlantic are less clearly associated with large-scale atmospheric forcing on decadal timescales than in the Northwest Pacific or Southern Ocean regions. In an attempt to gain further insight into how atmospheric anomalies force the SST anomalies in the regions where ocean processes are important, we analyze lagged regressions. Cane et al. (2017) advocated that low-pass filtering can induce false lead-lag relationships; we therefore opt for an alternative method. Specifically, atmospheric SLP anomalies averaged over the interval t = [t 0 − , t 0 ] are regressed onto the normalized SST index averaged over the interval t = [t 0 + 1, t 0 + 1 + ]. This gives the SLP and turbulent heat flux anomalies that are leading the (normalized) SST indices. Here is the length of the averaging period (in years), used to characterize the timescale of the atmospheric forcing and SST response. The lagged regression maps for two example timescales, = 3 years and = 10 years, are shown for the Northwest Pacific and North Atlantic SST regions in Figure 3.
On the = 3-year timescale, in both the observations and CMIP5 models, there are anticyclonic anomalies over the extratropical North Pacific leading the SST index (Figures 3a, 3e, and S3e). This shows that the change in the wind-driven circulation by poleward shift in the westerlies leads to warmer SSTs in the Northwest Pacific region, consistent with increased poleward advection of warm water through the subtropical western boundary current. There are also positive turbulent heat flux anomalies in the Northwest Pacific region in both the observations and CMIP5 models, which is likely caused by the westerlies shifting poleward such that there is more upward turbulent heat flux over the Northwest Pacific region. In the North Atlantic, the SLP anomalies leading the North Atlantic SST index are relatively weak and inconsistent between the observations and CMIP5 models (Figures 3c, 3g, and S3g)-though there are some warm turbulent heat flux anomalies leading the SSTs in both, albeit only significant in the subpolar gyre region in the CMIP5 models. where the SLP anomalies are significant at the 5% level, stippling in (e)-(h) indicates where more than 75% of the models exhibit regression coefficients of the same sign. The turbulent heat flux anomalies in (a)-(d) are only shaded where they are significant at the 10% level and in (e)-(h) are only shaded where more than 80% of the models exhibit regression coefficients of the same sign. The black boxes indicate the region used to calculate the SST index in each panel. The lime green boxes indicate the regions used to define the Aleutian Low and Icelandic Low SLP indices. The bottom row shows the correlation between anomalous SLP indices averaged over the interval t = [t 0 − , t 0 ] regressed onto the SST index averaged over the interval t = [t 0 + 1, t 0 + 1 + ], for (i) the Northwest Pacific and (j) the North Atlantic. The correlation for each timescale-in number of years, -is shown for the reanalysis (blue); each CMIP5 preindustrial control simulation (gray points) and the median of the preindustrial control ensemble (black line); each CMIP5 historical simulation (pink points) and the median of the preindustrial control ensemble (red line). The blue circles indicate where the correlation in the reanalysis is significant at the 5% level. SST = sea surface temperature; SLP = sea level pressure; CMIP5 = Coupled Model Intercomparison Project Phase 5.

10.1029/2018GL079077
On the 10-year timescale, the anticyclonic anomalies over the North Pacific still lead the warm SST anomalies in the Northwest Pacific region (Figures 3b, 3f, and S3f ), although these circulation anomalies are not statistically significant in the observations and occur in a smaller proportion of the models. Over the North Atlantic, however, there is a striking difference between the 10-year and the 3-year timescales. In both the observations and CMIP5 models, strong negative SLP anomalies over the northern part of the North Atlantic basin-a deepening of the Icelandic Low-leads warm SST anomalies (Figures 3d, 3h, and S3h). One notable difference, however, is the large-scale hemispheric wave train pattern seen in the observations, whereas in the models the SLP pattern looks similar to the Northern Annular Mode (NAM). This difference could simply be due to aliasing of decadal signals in the relatively short observations. However, it is also possible that the models are missing a teleconnection between the Pacific and Atlantic basins on decadal timescales associated, for example, with the Interdecadal Pacific Oscillation (B. Dong & Dai, 2015;Henley et al., 2015;Mantua et al., 1997). Indeed, the recent study by Henley et al. (2017) showed that the CMIP5 models underestimate the amount of decadal variance in Pacific SSTs. A common feature of both the observations and models is the strong turbulent heat fluxes out of the ocean over the North Atlantic, associated with the anomalous westerly and northwesterly winds during the 10 years preceding the warm SSTs in the North Atlantic region, meaning that oceanic processes are likely acting to provide additional heat in the mixed layer in the following 10 years. This demonstrates that the dominant timescales governing the development of the North Atlantic SSTs are seemingly significantly longer than in the Northwest Pacific, indicating that different oceanic processes are responsible.
We performed similar analysis over the South Pacific and South Atlantic SST indices ( Figure S4). On the 3-year timescale, the atmospheric forcing closely resembles the decadal regression maps (i.e., Figures 2 and S2) in both the observations and CMIP5 models. These results show that the SST anomalies on these timescales follow shifts in the local westerlies and are most likely due to the accompanying changes in the meridional Ekman transport as described previously, though the turbulent heat flux response is quite inconsistent between the CMIP5 models (in both the preindustrial control and historical ensembles). On the 10-year timescale there is little coherence between the observations and the CMIP5 models. For example, the observations show very large anomalies over the South Pacific, but there are negligible anomalies in the CMIP5 models. These are likely due to spurious decadal trends in the observations owing to the very sparse observations in this region and the performance of model in the region, owing to these inconsistencies we will not discuss these results in further detail.
To further investigate the dominant timescale of the atmospheric forcing of multiyear SST anomalies in the Northwest Pacific and North Atlantic SST regions, we performed lagged correlation analysis for a range of . We calculated area-averaged SLP indices over the Aleutian Low region (for the Northwest Pacific SST index) and the Icelandic Low region (for the North Atlantic SST index), as shown by the green boxes in Figure 3. The SLP indices were then averaged over the interval t = [t 0 − , t 0 ] and correlated with the SST indices averaged over the interval t = [t 0 + 1, t 0 + 1 + ], such that the SLP indices are leading the SST indices. These lagged correlations are shown in Figures 3i and 3j for timescales ( ) from 1 to 40 years.
In the observations, the SLP anomalies in the Aleutian Low region-where positive anomalies represent a weakening of the Aleutian Low-are positively correlated with the Northwest Pacific SSTs that follow. The positive correlation peaks for timescales between = 2 and 6 years, but the correlation drops off substantially beyond = 10 years (Figure 3i). Similar behavior is seen in the ensemble mean of both the preindustrial control and historical CMIP5 models, albeit peaking at slightly shorter timescales of = 2-3 years. The individual models themselves span a wide range of correlation values in the preindustrial control simulations (gray dots) and even more so in the historical simulations (red dots); however, the observed correlations that are significant for = 6-7 years are close to the upper end of the model range. The SLP anomalies in the Icelandic Low region are negatively correlated with the North Atlantic SSTs, with the magnitude of the observed correlation being largest (r ≈ −0.6) at a timescale of = 20 years. The correlation of the CMIP5 models are strongest at a timescale of about = 10 years, somewhat less than in the observations. However, negative correlations are found across a wide range of timescales, which are clearly longer than the timescale in the North Pacific. Interestingly, in the North Atlantic none of the preindustrial control simulations and only one of the historical simulations exhibit a stronger correlation than the observations at = 20 years, despite the wide range of correlations. This suggests that while the models are qualitatively capturing the correct sign of covariability in the North Atlantic, which is encouraging, the processes that determine the timescale of the North Atlantic response to anomalous forcing by the Icelandic Low are not correctly represented in the CMIP5 models.
For completeness, we also performed this analysis on the South Pacific region. The meridional SLP gradient across the South Pacific region clearly leads SST anomalies on timescales between = 1 and 3 years in both the models and observations and weakens on longer timescales ( Figure S3). Again, however, none of the preindustrial control simulations and only one of the historical simulations exhibits a relationship as strong as in the observations for timescales of = 2 and 3 years.

Summary and Discussion
In this study, we have analyzed the relationship between decadal SST and turbulent heat fluxes to identify where oceanic processes may play a nonnegligible role in extratropical decadal SST variability. Positive correlations between decadal turbulent heat flux and SSTs indicate an active role of oceanic processes over regions in the North Atlantic, Northwest Pacific, Southern Pacific, and Southern Atlantic (Figure 1). Of particular interest is the contrasting nature of the oceanic influence on decadal SST anomalies in the Northwest Pacific and North Atlantic (i.e., Figure 2). Over the Northwest Pacific, SST anomalies respond to atmospheric forcing via changes in the horizontal wind-driven gyre circulation on timescales of 3-7 years, in both observations and CMIP5 models (Figure 3i). Over the North Atlantic, SST anomalies are also influenced by atmospheric circulation anomalies, though the response is stronger at longer timescales than in the Northwest Pacific-peaking at about 20 years in the observations, compared to only 10 years in the CMIP5 models ( Figure 3j).
The longer timescales in the North Atlantic are not consistent with solely a wind-driven gyre response, as in the Northwest Pacific. Instead, upward turbulent heat fluxes over the subpolar gyre (particularly the Labrador Sea region) in the preceding decades seem particularly important prior to decades with warm North Atlantic SST (i.e., Figure 3). Robson et al. (2012) and Williams et al. (2015) found that anomalous northward ocean heat transport in recent decades was largely a response to anomalous buoyancy forcing in the subpolar gyre, associated with turbulent heat flux out of the ocean. In addition, recent modeling simulations indicate that most of the decadal SST variability over the recent observational period can be reproduced in models that are forced with the large-scale observed buoyancy flux forcing (Delworth et al., 2017). While traditionally the North Atlantic subpolar gyre SST response to buoyancy forcing has often been interpreted as due to changes in the meridional overturning circulation (e.g., Delworth & Zeng, 2016), some of the observed variability has been attributed to changes in the horizontal gyre circulation (Piecuch et al., 2017;Williams et al., 2014), though their relative contributions are unclear. Finally, intrinsic variability from entrainment, overflow, and mixing affecting the deep mixed layer in the subpolar North Atlantic might be influencing the region on decadal timescales contributing to the upward heat fluxes.
Despite the qualitative consistencies between the CMIP5 models and the observations that have been highlighted here, there are some potentially important differences. The correlation between the decadal turbulent heat flux and SST over the North Atlantic in the CMIP5 models is weaker than in observations. More generally, the atmospheric circulation anomalies that lead the SST anomalies exhibit stronger relationships in the observations than the vast majority of models-in both the North Atlantic and Northwest Pacific-and also exhibit longer characteristic timescales of variability (i.e., Figure 3). The disparity between observations and the CMIP5 models could potentially be due to weaker atmospheric circulation anomalies on decadal timescales, as seen in most CMIP5 models (X. Wang, Li, et al., 2017) and also in other coupled models (Kim et al., 2018;Parker et al., 2007). While the lack of low-frequency atmospheric variability may explain why the dominant timescales of variability are too short, it does not account for the relatively weak response in the models. We hypothesize that the weak model response is related to deficiencies in the oceanic model components and to overly strong oceanic damping processes (including air-sea feedbacks and mixing)-particularly in the North Atlantic midlatitude and subpolar gyre region, where there are significant model biases (e.g., Heuzé, 2017). Further investigation of the mixed-layer dynamics in climate models-and the fidelity with which the relevant processes compare to observations-is necessary to isolate the source of the relatively weak decadal SST anomalies in response to large-scale atmospheric circulation anomalies.