Mechanisms of Decadal North Atlantic Climate Variability and Implications for the Recent Cold Anomaly

: Decadal sea surface temperature (SST) ﬂuctuations in the North Atlantic Ocean inﬂuence climate over adjacent land areas and are a major sourceof skill in climate predictions. However,the mechanismsunderlying decadal SST variability remain to be fully understood. This study isolates the mechanisms driving North Atlantic SST variability on decadal time scales using low-frequency component analysis, which identiﬁes the spatial and temporal structure of low-frequency variability. Based on observations, large ensemble historical simulations, and preindustrial control simulations, weidentifyadecadalmodeofatmosphere–oceanvariabilityintheNorthAtlanticwithadominanttimescaleof13–18years. Large-scale atmospheric circulation anomalies drive SST anomalies both through contemporaneous air–sea heat ﬂuxes and through delayed ocean circulation changes, the latter involving both the meridional overturning circulation and the horizontal gyre circulation. The decadal SST anomalies alter the atmospheric meridional temperature gradient, leading to a reversal of the initial atmospheric circulation anomaly.The time scale of variability is consistent with westward propagation of baroclinic Rossby waves across the subtropical North Atlantic. The temporal development and spatial pattern of ob- served decadal SST variability are consistent with the recent observed cooling in the subpolar North Atlantic. This suggests that the recent cold anomaly in the subpolar North Atlantic is, in part, a result of decadal SST variability. observed mode of decadal SST variability). We then solve for the linear combination of the leading ﬁve model-based LFPs that maximizes the pattern correlation with the observational LFP. Similar results are obtained if we simply use the individual LFP (and associated LFC) that has the highest pattern correlation with the observational LFP. Before performingthelow-frequencycomponentanalysis,allmodeloutput is interpolated onto the HadISST (1 8 For CESM-LE, the ensemblemeanSST also from each ensemble at each prior to in on internal (unforced)


Introduction
The North Atlantic Ocean displays pronounced decadal variability ( Fig. 1; Deser and Blackmon 1993;Czaja and Marshall 2001). Decadal variations in North Atlantic sea surface temperature (SST) influence climate over adjacent continents and are a major source of skill in climate predictions (Hermanson et al. 2014;Msadek et al. 2014;Årthun et al. 2017;Yeager and Robson 2017). Understanding the physical mechanisms responsible is thus important for attributing current anomalies and predicting future changes in the North Atlantic sector.
The proposed mechanisms differ in the relative roles of ocean circulation changes and atmospheric forcing in setting upperocean heat content and SST anomalies, the importance of ocean gyre and overturning circulation, and whether decadal ocean temperature variability is part of a coupled atmosphereocean mode of variability or an ocean-only mode. The specific time scale of variability also ranges from less than 10 years to more than 30 years. One potential source of this discrepancy between studies is the frequent use of low-pass or bandpass filters to isolate the time scale of interest, which imposes a priori assumptions about the temporal structure of the variability. The use of low-pass filtered data can also lead to spurious low-frequency signals and complicate the detection of lead-lag relationships (Cane et al. 2017).
Here, we identify the time scale, spatial pattern, and mechanisms underlying decadal variability in the North Atlantic Ocean by using low-frequency component analysis (Wills et al. 2018(Wills et al. , 2019. This is a statistical method that identifies spatial anomaly patterns with the largest ratio of low-frequency variance to total variance without making a priori assumptions about the spatial pattern or the specific time scale (described in section 2). We apply this method to North Atlantic SST Denotes content that is immediately available upon publication as open access.
Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JCLI-D-20-0464.s1. anomalies in observations (section 3), a large ensemble of historical simulations (section 4), and unforced preindustrial control simulations with coupled climate models (CMIP5; section 5), and we use it to isolate the associated mechanisms of decadal SST variability (summarized in section 6). Importantly, the results of the analysis are not low-pass filtered, allowing insight into interactions across different time scales. We furthermore show that the identified mode of decadal SST variability has likely played an important role in the recent North Atlantic cold anomaly (Josey et al. 2018) (section 7). Conclusions are presented in section 8.

a. Observational data
To assess observed SST variability in the North Atlantic we use the Hadley Centre Sea Ice and Sea Surface Temperature dataset, version 1.1 (HadISST; Rayner et al. 2003). The North Atlantic Ocean is defined as 08-708N, 808W-158E (Fig. 1). Monthly data are available on a 18 grid starting in 1870, but are understood to be less reliable in the data-sparse periods before 1947. We thus consider data between 1948 and 2017; however, our results are similar for the full time series. The analysis has also been performed for other observation-based SST products (ERSST and EN4) and found to be consistent. We also use sea level pressure (SLP), net surface heat fluxes (turbulent and radiative), zonal winds at 200 hPa, and atmospheric temperature at 850 hPa from the NOAA Twentieth Century Reanalysis (Compo et al. 2011). Surface heat fluxes from the Objectively Analyzed Air-Sea Fluxes for the Global Ocean (OAFlux; Woods Hole Oceanographic Institution) were also investigated and showed consistent results. Indices of the North Atlantic Oscillation (NAO) and east Atlantic pattern (EAP) are defined as the two leading principal component time series of winter (November-April) SLP anomalies over the Atlantic sector.

b. CESM-LE
Further insights into the mechanisms of decadal atmosphere-ocean variability over the North Atlantic Ocean are obtained from the Community Earth System Model large ensemble (CESM-LE; Kay et al. 2015). These simulations use the fully coupled CESM1 model, which consists of the Community Atmosphere Model version 5; the Parallel Ocean Program, version 2 (POP2); the Community Land Model, version 4; and the Community Ice Code, version 4 (CICE4) (Hurrell et al. 2013). The spatial resolution of the ocean model is nominally 18 longitude 3 18 latitude, whereas the atmospheric model is 0.98 3 1.258. The CESM-LE includes 40 ensemble members for the time period 1920 to 2100. Here we use 30 ensemble members (#2-31) from the historical time period , that is, analyzing a total of 2580 years. All the analysis of CESM-LE presented in this paper is done for each ensemble member separately and then the results are averaged across the ensemble members. In addition to SST and SLP, we analyze ocean heat content, net surface heat fluxes, the ocean barotropic streamfunction, and the Atlantic meridional overturning circulation (AMOC). The AMOC strength is defined as the maximum of the zonally integrated meridional overturning streamfunction in the Atlantic basin.

c. CMIP5 models
We use output from preindustrial control (piCTRL) simulations from six global climate models available from phase 5 of the Coupled Model Intercomparison Project (CMIP5; Taylor et al. 2012): ACCESS1.0, CanESM2, CCSM4, CNRM-CM5, MPI-ESM-P, and NorESM1-M. The specific models were chosen based on the availability of SST, barotropic streamfunction, and meridional overturning streamfunction as output variables, and a simulation length of $500 years. We use the last 500 years of each simulation. The piCTRL simulations have fixed preindustrial levels of greenhouse gases and other external forcings, and thus allow for the investigation of internal climate variability. The aim of this study is not to present a full intercomparison of all CMIP5 models, but the six models nevertheless represent a wide range of model groups, ocean and atmosphere components, and resolutions.

d. Low-frequency component analysis
To identify decadal variability in North Atlantic SST we use low-frequency component analysis (LFCA; Wills et al. 2018Wills et al. , 2019, see also Schneider and Held 2001). LFCA isolates the low-frequency variability in a dataset by finding low-frequency patterns (LFPs)-that is, linear combinations of the leading empirical orthogonal functions (EOFs) that maximize the ratio of low-frequency to total variance in their corresponding time series, which are called low-frequency components (LFCs). The LFCs are sorted by the ratio of low-frequency to total variance, such that the leading LFCs describe modes of lowfrequency (decadal) variability. The LFCs are required to be orthogonal (uncorrelated), but the LFPs are not. A detailed description of LFCA is provided in Wills et al. (2019).
There are two free parameters in LFCA; the low-pass cutoff and the number of EOFs included. Here we define lowfrequency variance based on a linear Lanczos filter with a 10yr low-pass cutoff and reflected boundary conditions, in order to focus on decadal variability. For observed SST, 20 EOFs were used, retaining 94% of the total variance. The analysis has been performed with different numbers of EOFs and with different low-pass cutoff frequencies, and the spatial patterns and time scales are robust [see also Wills et al. (2018) for a detailed discussion of the robustness of LFCA to the choice of parameters]. It is important to note that although a low-pass filter is used to define the linear combination of EOFs, the resulting LFCs are unfiltered and can thus be used to identify lead-lag relationships at annual time scales.
For simulated SST (CESM-LE and CMIP5 models) we use 50 EOFs, such that the fraction of variance explained is the same for observations and models (;95% of the total SST variance is retained in each model). To investigate decadal SST variability that resembles narrow-band variability found in observations (see section 3), we solve for LFPs/LFCs that maximize the ratio of decadal to total variance based on a 10-30-yr bandpass filter (noting that using such a bandpass filter for the observational LFCA does not influence the characteristics of the observed mode of decadal SST variability). We then solve for the linear combination of the leading five model-based LFPs that maximizes the pattern correlation with the observational LFP. Similar results are obtained if we simply use the individual LFP (and associated LFC) that has the highest pattern correlation with the observational LFP. Before performing the low-frequency component analysis, all model output is interpolated onto the HadISST grid (18). For CESM-LE, the ensemble mean SST is also subtracted from each ensemble member at each grid point prior to analysis in order to focus on internal (unforced) variability.
All the analysis is based on winter (November-April) averages. We focus on winter because SST captures upper-ocean heat content variability; the average correlation between winter SST and upper-ocean (0-800 m) heat content in the subpolar North Atlantic in CESM-LE is 0.90.
When computing lead-lag regressions between the LFCs and other variables, significance levels are computed based on a random phase test (Ebisuzaki 1997). For CESM-LE, phase randomization is applied to the concatenated multimember time series such that it also randomizes phase across different ensemble members.

Observed North Atlantic decadal variability
The two leading LFPs/LFCs of observed North Atlantic SST anomalies correspond to basin-wide multidecadal variability ( Fig. S1 in the online supplemental material), closely resembling those identified in Wills et al. (2019) using monthly observations (ERSST) between 408S and 758N. Wills et al. (2019) found that the leading mode represents the impact of global warming, whereas the second mode corresponds to the Atlantic multidecadal oscillation (AMO). For the observations presented here, the correlation between LFC-2 and the AMO (calculated as detrended SST between 08 and 608N) is 0.70. The third LFP/LFC (Fig. S1) is also associated with multidecadal variability (;30 yr), but has a pattern more centered near the boundary between the subtropical and subpolar gyres. The time scale and pattern of LFP-3 are similar to the interdecadal SST variability discussed in Eden and Jung (2001). The leading three modes of low-frequency variability explain 38%, 15%, and 8% of the total low-frequency (.10 years) variance, respectively, with local explained variance of .50% (Fig. S1).
The fourth mode of observed low-frequency SST variability in the North Atlantic (Fig. 2) is associated with a decadal (13-16 years) time scale, consistent with that found for variations in ocean temperature along the North Atlantic Current (Årthun et al. 2017). This mode explains 10% of the total low-frequency variance in North Atlantic SST, with local explained lowfrequency variance of ;40% where the LFP has its maximum amplitude, such as in the eastern subpolar North Atlantic (Fig. S1). We will refer to this LFP/LFC as Atlantic decadal variability (ADV) throughout the rest of the text.
A positive phase of the ADV is characterized by warm anomalies in the subpolar North Atlantic (SPNA) and in the eastern subtropical Atlantic, and a cold anomaly in the western subtropical gyre. The spatiotemporal SST development associated with the ADV index is shown in Figs. 3a-e (all lags between 66 years are displayed in Fig. S2). Following an ADV maximum (Fig. 3c), the SPNA cools whereas the western subtropical gyre warms, leading to a reversal of the pattern after approximately 6-8 years (Fig. 3e).
The atmospheric circulation associated with the ADV shows a dipole SLP pattern with anomalous high pressure over the SPNA and a low pressure anomaly over the subtropics at lag 0, and vice versa at 66 years (Figs. 3f-j). The evolution of the SLP pattern leading up to an ADV maximum projects onto a combination of the NAO and EAP, whereas the lagged response is mainly captured by the EAP (Fig. 4).
The air-sea heat flux anomalies associated with the ADV are shown in Figs. 3k-o. Over the eastern SPNA, the warming leading up to an ADV maximum is associated with anomalous heat fluxes from the ocean to the atmosphere (Figs. 3b,l), suggesting that ocean circulation changes and anomalous ocean heat transport are important drivers of decadal SST anomalies in this region. This is also the case south of the Gulf Stream extension region (308N), where the warming following an ADV maximum is associated with ocean heat loss (Figs. 3d,n). However, surface heat fluxes also play an active role in driving the SST anomalies, such as by acting to warm the western SPNA leading up to an ADV maximum (Fig. 3m).
The surface heat flux patterns also hint at a reorganization of the circulation at the boundary between the subpolar and subtropical gyre (the intergyre region) as a response to atmospheric circulation changes: The anomalous ocean heat gain to the north of the Gulf Stream extension region and anomalous heat loss to the south leading up to an ADV maximum (the approximate position of the Gulf Stream is indicated by the maximum climatological heat loss in Fig. 3m) is consistent with a southward shift of the North Atlantic Current as a result of atmospheric circulation anomalies ( Fig. 3h; e.g., Marshall et al. 2001;Nigam et al. 2018). This will be further explored in section 4.
The reversal of the SST and SLP patterns after 6 years suggests the possibility of a coupled mode of atmosphereocean variability over the North Atlantic Ocean, where the ocean modulates the low-frequency evolution of the SLP pattern. To assess this hypothesis and to identify the underlying mechanisms, we now investigate low-frequency SST variability in coupled climate models.

a. Drivers of decadal SST variability
The pattern of North Atlantic decadal SST variability between 1920 and 2005, as simulated by CESM-LE (see section 2), is shown in Fig. 5. The spatial correlation with the observed ADV pattern is high (0.86), and the simulated ADV index has a time scale (15-18 years) similar to observations. We note that the ADV amplitude in CESM is weaker than in observations, consistent with North Atlantic low-frequency variability in CESM being underestimated (Kim et al. 2018). The explained low-frequency variance is also lower than in observations, but with a similar spatial pattern (maximum values of .20% in the eastern SPNA; Fig. S3). The spatiotemporal development of simulated SST, SLP, and air-sea heat fluxes associated with the ADV index is shown in Figs. 6a-o. For SST, all lags between 68 years are displayed in Fig. S4. As in the observations, the model shows a dipole pattern with anomalous low pressure over the midlatitudes and a high pressure anomaly over the SPNA and Nordic seas during a positive phase of the ADV (Fig. 6c). The low-frequency evolution of the SLP pattern associated with the ADV is mainly captured by the EAP, and not by the NAO (Fig. 4).
We note that upper-ocean salinity associated with the ADV shows a similar pattern to that of SST ( Fig. S5), consistent with covarying low-frequency temperature and salinity in the North Atlantic (Hall and Manabe 1997;Zhang 2017). Density anomalies in the North Atlantic are dominated by temperature anomalies (Menary et al. 2015a;Gastineau et al. 2018;Danabasoglu et al. 2019), and ADV-related salinity variations are thus not further assessed here.
The simulated surface heat fluxes associated with the ADV (Figs. 6k-o) show many common features to the spatial patterns from observations (Fig. 3). Specifically, the surface heat flux anomalies during the mature phases of the ADV (lags 68 and 0 years) are characterized by opposite signs over the western SPNA and western subtropical gyre regions. In contrast to observations, anomalous fluxes in CESM-LE also extend into the eastern SPNA. As in observations, the relationship between SST and surface heat fluxes suggests that both the atmosphere and ocean contribute to the formation of SST anomalies. In support of ocean-driven SST variability, we note especially that warm (cold) SST anomalies in the Gulf Stream extension region are associated with anomalous oceanic heat loss (gain).
To quantify the roles of oceanic and atmospheric forcing, in  (Tulloch and Marshall 2012;Buckley and Marshall 2016). We calculate the ocean heat supply (heat transport convergence) in these two regions by subtracting the surface heat fluxes from the heat content tendency (a centered difference on seasonal means used to estimate the time derivative). During the years leading up to an ADV maximum (lag 0), heat is supplied by ocean heat transport to the SPNA (Figs. 6p,q and 7a). Note in particular that the ocean is warming the SPNA 3-7 years before the maximum ADV index, while the surface fluxes act to reduce the warming (Fig. 7a). The surface heat flux then changes sign and acts to amplify the upper-ocean warming three years before the ADV maximum. The subsequent cooling is mainly driven by the ocean. This result is consistent with a recent heat budget analysis with a reanalysis-forced ocean-sea ice CESM simulation that finds that the heat transport convergence is the dominant term in the upper-ocean heat budget in the SPNA on multiyear time scales (Yeager 2020).
The ocean is also active in driving ADV-related temperature changes in the Gulf Stream extension region (Fig. 7b). Leading up to an ADV maximum, the Gulf Stream extension region cools mainly as a result of anomalously low ocean heat supply.
Surface heat input then contributes to the subsequent warming (lags 0-3 years), after which ocean supply once again dominates the heat budget. In agreement with previous studies (e.g., Buckley et al. 2014), the Gulf Stream extension region is a more heavily damped system where ocean heat transport convergence and surface heat fluxes are strongly anticorrelated (20.75). These results demonstrate that the spatial SST pattern associated with the ADV is not just the time-integrated response to surface heat flux anomalies, and that ocean heat transport plays a key role in decadal SST variability.
The ocean heat supply term in the calculated heat budget includes all mechanisms of heat supply by the ocean (advective heat transport, heat transport by parameterized eddies, and diffusive fluxes). A further decomposition of the advective heat transport into components associated with the mean and timevarying velocity and temperatures is beyond the scope of the present study. However, Gervais et al. (2018) find that changes in ocean heat supply in the subpolar North Atlantic in CESM-LE are almost entirely due to changes in (resolved) ocean velocities. The importance of anomalous advective transport to ocean heat supply in the North Atlantic is also corroborated by other models (e.g., Eden and Jung 2001;Menary et al. 2015a;Foukal and Lozier 2018), although anomalous temperatures become important along the Gulf Stream (Dong and Kelly 2004;Doney et al. 2007).
Ocean circulation changes in the North Atlantic result from changes in both the Atlantic meridional overturning circulation (AMOC) and the wind-driven gyre circulation (e.g., Williams et al. 2014). During an ADV minimum (lag 28), the . Dots indicate where the correlation is significant at the 90% (observations) and 95% (CESM-LE) confidence level (Ebisuzaki 1997). The NAO and EAP index were calculated by an EOF analysis of sea level pressure. The inset figures show the simulated NAO and EAP patterns, respectively (SLP contours spaced at 25 Pa). depth-integrated circulation in the Gulf Stream extension region is anomalously strong (Figs. 8a and 9a) and the intergyre boundary is shifted northward, consistent with anomalous ocean heat transport by the subtropical gyre into the western SPNA (Figs. 6p and 7a). At the same time, the AMOC is in a weak state. In subsequent years the AMOC strengthens whereas the circulation in the Gulf Stream extension region weakens, although remaining anomalously strong for another four years (Fig. 9a). The ocean heat supply to the SPNA is thus highest 7 to 3 years before an ADV maximum (Fig. 7a) when both the gyre and overturning circulations are anomalously strong.
The increased northward heat transport in the Gulf Stream extension region as a result of the atmospheric circulation anomalies at lag 28 resembles the fast barotropic response to NAO-like winds (e.g., Eden and Willebrand 2001). The strengthening of the AMOC from lag 28 to lag 24 (Figs. 8f,g and 9b) can be understood as a response to the anomalous atmospheric forcing at lag 28 years, which drives a large surface heat loss in the SPNA (Fig. 6k). This anomalous heat loss influences the AMOC through deep convection in the Labrador Sea, in CESM (Maroon et al. 2018;Danabasoglu et al. 2019) and other models (e.g., Medhaug et al. 2012;Delworth and Zeng 2016), which is reflected in a deepening of the winter mixed layer (Fig. 9c). The influence of the AMOC at 408N on ADV is maximum at lag 22, lagging the surface heat loss and convection by 5-6 years. The latter time lag is consistent with the southward propagation of density anomalies from the Labrador Sea to 408N and consequent AMOC response (Tulloch and Marshall 2012;Zhang and Zhang 2015;Jackson et al. 2016).

b. Ocean-atmosphere interaction
The dipolar SLP pattern associated with the ADV reverses sign between a positive and negative phase of the ADV (Figs. 6f-j). Although the NAO plays a role in forcing the ADV similar to that of the EAP, the SLP response is more closely related to the EAP (Fig. 4). The lagged response of the EAP to ADV suggests a possible feedback between decadal SST variability and the atmospheric circulation.
To further assess the potential influence of SST on the atmosphere we show in Figs. 10a-c the meridional atmospheric temperature gradient at 850 hPa, which provides the energy source for baroclinic eddies and is a dominant factor in setting the strength of the upper-tropospheric zonal winds (e.g., Brayshaw et al. 2011). The climatology exhibits a negative gradient (dT/dy , 0), which reaches a maximum in the Gulf Stream extension region. During a positive ADV, the western subtropical gyre is anomalously cold and the SPNA anomalously warm, thus weakening the meridional temperature gradient (Fig. 10a) and zonal winds (Fig. 10c) in the central North Atlantic. In the following years the ADV pattern reverses sign and the atmospheric temperature gradient strengthens over the Gulf Stream extension region, consistent with increased westerly winds (Figs. 10b,d) and a shift toward positive NAO/EAP conditions. The change in temperature gradient and zonal winds over the Gulf Stream extension region associated with one standard deviation of the ADV index is 13% and 14%, respectively, of the interannual variability (standard deviation). However, this is most likely a lower bound as climate models are known to underestimate the strength of ocean-atmosphere interaction (Kim et al. 2018;Czaja et al. 2019). In line with this, observations show stronger ADV-related changes in the meridional atmospheric temperature gradient and zonal winds (Fig. 10); 34% and 25%, respectively, of the interannual variability over the Gulf Stream extension region. We note, however, that the observed and simulated spatial patterns of the atmospheric temperature gradient and zonal winds associated with the ADV are very similar. The sea level pressure response to SST anomalies is also stronger in observations (Fig. 3) than in CESM-LE (Fig. 6), reflected in weaker correlations between the ADV index and the NAO/EAP in CESM-LE (Fig. 4). In observations, the SLP response over the SPNA to changes in the ADV is 150 Pa for SST anomalies of approximately 0.38C (Fig. 3), which is in general agreement with Czaja and Frankignoul (2002), who estimated the strength of SST forcing on the NAO to be 200-300 Pa K 21 . change in ADV is associated with adjustments of both the gyre and overturning circulation. The ocean adjustment-taken as the time from a positive to a negative circulation anomaly (Figs. 9a,b)-is 7-8 years; that is, half the period of variability.
Several mechanisms have been proposed to be responsible for the time scale of ocean adjustment to variable wind and buoyancy forcing, one of which is the westward propagation of baroclinic Rossby waves across the North Atlantic Frankcombe et al. 2010;Muir and Fedorov 2017). To assess this mechanism, we apply a complex principal component analysis (Horel 1984) to the SST data (keeping in mind that winter SST reflects upper-ocean heat content). This analysis detects traveling waves in the input time series and orders the dataset into modes of phase propagation in space and time according to the variance explained. Note that a 10-30-yr bandpass filter was applied before the analysis in order to focus on SST anomalies associated with ADV.
In both observations and the CESM-LE, westward-propagating SST anomalies are captured by the second mode of variability, explaining 27% and 30% of the variance, respectively (the first mode is well separated and represents eastward propagation in the direction of the mean current). Temperature anomalies propagate from east to west, from 208 to 608W, across the subtropical North Atlantic (308-408N) in 6-7 years (Fig. 11), in agreement with the transit time of Rossby waves at this latitude (Sturges et al. 1998;Tulloch and Marshall 2012). In support of these westward-propagating SST anomalies playing a role in the evolution of the ADV, there is a strong correlation between the ADV index and the temporal development of the identified mode of propagation (Fig. 11d).
Baroclinic Rossby waves can be excited by the Ekman response to large-scale atmospheric circulation anomalies (Anderson and Gill 1975;Sturges et al. 1998) and by internal ocean dynamics (Colin de Verdière and Huck 1999;te Raa and Dijkstra 2002;Arzel et al. 2018). Evidence in support of ADV-related winds forcing westwardpropagating Rossby waves is presented in Fig. 12, which shows that the spatial wind stress curl patterns associated with changes in the ADV and westward-propagating SST anomalies are highly similar. Note that the patterns reverse sign if the wind stress curl leads/lags by 8 years. A change from negative to positive wind stress curl anomalies leading up to a positive ADV is consistent with a weakening of the depth-integrated circulation in the Gulf Stream extension region (Fig. 9a) and a cooling of the western subtropical North Atlantic. In agreement with, among others, Grötzner et al. (1998) and Wu and Liu (2005), our results are thus suggestive of wind-driven Rossby waves setting the time scale of ocean adjustment, and, hence, the time scale of the ADV. Dedicated model sensitivity experiments would be required to establish that this mechanism is the one that gives rise to the ADV, rather than internal ocean dynamic. However, we note that a one-quarter phase lag between upper and subsurface temperature anomalies in the region of Rossby wave formation (here 308-408N, 208-308W; Fig. 11)-a key fingerprint of internally generated Rossby waves by large-scale baroclinic instability )-is not detected.

Robustness of mechanism across CMIP5 models
The robustness of the mechanism driving decadal Atlantic SST variability in observations and CESM-LE is further assessed by analyzing five CMIP5 preindustrial control simulations. Although the simulated ADV patterns differ among the models, all models show a pattern with a warm SPNA and a cold western subtropical gyre (Figs. 13a-f). The strength of the  (Ebisuzaki 1997). spatial correlation with the observed ADV pattern ranges from 0.50 (NorESM1-M) to 0.67 (MPI-ESM-P). The time scale of variability associated with the ADV in CMIP5 models is 14-18 years. Note that the LFCA was also performed for ACCESS1.0, but the ADV pattern did not resemble that in observations (spatial correlation of 0.26). The underlying mechanism was therefore not assessed for this model. The local explained low-frequency variance is similar to that in CESM- LE with maximum values of .20% mainly in the eastern SPNA (Fig. S6).
As in observations and CESM-LE, a positive ADV is associated with a dipolar SLP pattern, with anomalous high (low) pressure over the subpolar (subtropical) region. Consistent with CESM, all models also show a weakened subtropical gyre when the ADV index is at its maximum (Figs. 13g-l), although the magnitudes of the anomalies vary (note also that the substantial differences in the climatological gyre and overturning circulations). The reverse patterns are found when the ADV is at a minimum (not shown). CMIP5 models also show the strongest AMOC anomalies around 408N 2-3 years before a maximum ADV. The exception is CNRM-CM5, which has the strongest AMOC anomalies 8 years before, and negative values at lag 22. The reason for this different behavior is not pursued, but we note that, consistent with a weak climatological AMOC (Fig. 13p), convection in the Labrador Sea is weak in CNRM-CM5 (Heuzé 2017), which could influence how decadal buoyancy anomalies associated with the ADV are manifested in the AMOC.
As in observations and CESM-LE, the time scale of half an ADV cycle (7-9 years) is consistent with the travel time of westward-propagating SST anomalies across the North Atlantic (Figs. 13s-x). Westward propagation speed, however, varies significantly with longitude in all of the models. This was also noted in the analysis by Muir and Fedorov (2017), and may be a result of strong eastward flow in some parts of the basin, which could mask the westward-propagating anomalies.
In general, the multimodel mean patterns of SST, SLP, and ocean circulation are similar to those found in CESM-LE. CMIP5 preindustrial control simulations thus support the FIG. 11. Westward propagation of (a) observed (HadISST) and (b) simulated (CESM-LE) ocean temperature anomalies (in standard deviations) across the North Atlantic (208-608W) at 308-408N inferred from a complex principal component (CPC) analysis (Horel 1984). (a),(b) The spatiotemporal SST variability associated with the second mode of variability for observations and one CESM ensemble member. The second mode of variability explains 27% and 30% of the total variance in observations and CESM-LE, respectively. (c) The corresponding spatial phase plots (averaged across all CESM ensemble members) illustrating the characteristic westward propagation for all SST anomalies. A 10-30-yr bandpass filter was applied before the CPC analysis. (d) Lead-lag correlation between the ADV index and westward-propagating SST anomalies at 208W in CESM-LE. The blue line shows the multimember mean and the shading indicates the interquartile spread. Dots indicate where the correlation is significant at the 95% confidence level (Ebisuzaki 1997). existence of a mode of decadal North Atlantic variability, which is characterized by contrasting SST changes in the SPNA and western subtropical gyre, and that involves changes in both gyre and overturning circulation.

Summary and discussion of the proposed mechanism
The mechanism underlying the mode of decadal North Atlantic variability identified here-both observed and modeled-is summarized in Fig. 14. During a negative phase of the ADV the SPNA and eastern subtropical gyre are anomalously cold, whereas the western subtropical North Atlantic is anomalously warm. The atmospheric response to these SST anomalies is a strengthened meridional temperature gradient and an SLP dipole with centers over the subpolar (low pressure) and subtropical North Atlantic (high pressure). The associated changes in air-sea fluxes first act to reinforce the SST anomalies, but the anomalous atmospheric circulation also leads to an expansion of the subtropical gyre and a strengthened circulation in the Gulf Stream extension region. The latter increases the poleward heat transport, acting to warm the SPNA. The reversal of the initial SST pattern is also driven by a strengthened meridional overturning circulation, as a result of anomalous surface heat loss and deep convection in the Labrador Sea. The reversed SST pattern (positive ADV) weakens the meridional temperature gradient and the westerly winds, consistent with a reversal of the anomalous SLP dipole. The oceanic adjustment to atmospheric circulation anomalies-and, hence, the time scale of the decadal oscillation-is consistent with westward propagation of Rossby waves across the North Atlantic. Decadal SST variability in the North Atlantic is thus partly forced by the atmosphere, but the ocean is responsible for setting the preferred time scale of variability.
There is, in general, good agreement between the observed and simulated response to ADV. This includes the spatial patterns of SST, SLP, surface heat fluxes, and changes in the atmospheric temperature gradient, although the simulated atmospheric response is weaker than observed. The evolution of the simulated patterns is, in general, also more symmetric (around maximum ADV) than in observations. One possible explanation is that in CESM-LE, internal variability has been separated from the forced response by removal of the ensemble mean, whereas observations will be a mix of both. The CESM-LE analysis is also based on a much longer (combined) time series. The consistency between observed and simulated decadal SST variability and its drivers nevertheless provides confidence in the results presented here.
The ingredients of the mechanism presented in Fig. 14 are largely consistent with previous findings on North Atlantic decadal (10-20-yr time scale) variability, including the key role of ocean gyre adjustment to variable wind forcing (Grötzner et al. 1998;Wu and Liu 2005;Nigam et al. 2018) and the importance of meridional overturning circulation anomalies as a response to air-sea heat fluxes in the SPNA (Häkkinen 2000;Martin et al. 2019). Consistent with an important role for the AMOC in driving decadal SST variability, the antiphase SST relationship between the SPNA and Gulf Stream extension region projects well on the observed decadal AMOC fingerprint (Zhang 2008). More generally, the mechanisms identified here support an important role of ocean dynamics in low-frequency North Atlantic SST variability (Zhang 2017;Wills et al. 2019;Zhang et al. 2019;Yeager 2020).
In agreement with previous studies, our results suggest that the time scale of the SST variability is set by the transit time of baroclinic Rossby waves across the basin (Grötzner et al. 1998;Häkkinen 2000;Marshall et al. 2001;Wu and Liu 2005). Other studies have related the time scale to propagation of temperature (heat content) anomalies along the Gulf Stream and North Atlantic Current Ruiz-Barradas et al. 2018) or around the subpolar gyre (Menary et al. 2015a), or to the accumulation time (memory) of heat anomalies in the SPNA (Martin et al. 2019), which we do not rule out here.
Our results highlight the importance of ocean-atmosphere interaction in maintaining decadal SST variability (Grötzner FIG. 12. Regression of wind stress curl in CESM-LE onto (a) the ADV index and (b) westward-propagating SST anomalies at 208W (as identified from the CEOF analysis in Fig. 11). Units are 10 28 N m 3 per standard deviation of the respective indices. The black contour line in (a) shows the climatological zero line for the barotropic streamfunction, indicating the boundary between the subtropical and subpolar gyres. Häkkinen 2000;Wu and Liu 2005;Nigam et al. 2018;Martin et al. 2019). Although the feedback of SST on the atmosphere has not been assessed in detail, the decadal changes in SST and SLP associated with ADV are consistent with the influence of North Atlantic SST on the NAO through changes in the meridional SST gradient (Czaja and Frankignoul 2002;Gastineau and Frankignoul 2015;Baker et al. 2019). The importance of the feedback of the atmosphere onto the SST anomalies can also not be quantified here. However, in section 4c above we show that decadal changes in winds related to the ADV drive an ocean circulation response, mediated by westward-propagating Rossby waves. The ADV thus has the potential to be a coupled mode of variability, but a comparison between coupled and forced experiments (e.g., Farneti and Vallis 2011;Gastineau et al. 2018;Larson et al. 2018) would be required to further strengthen this hypothesis.

Implications for recent observed changes in the subpolar North Atlantic
Large changes have been observed in the SPNA during recent years (Robson et al. 2018;Chafik et al. 2019), manifested in a reversal of upper-ocean and SST trends between the periods 1994-2004 (warming) and 2005-15 (cooling) (Robson et al. 2016;Piecuch et al. 2017;Ruiz-Barradas et al. 2018). Consistent with the recent observed cooling over the SPNA, the ADV index shows a large drop from a positive phase in the late 2000s to a negative phase in 2015 (Fig. 2). The spatial pattern of observed cooling also resembles the ADV pattern (Figs. 15a,b), suggesting that the recent observed cooling is, in part, related to decadal-scale SST variability as described here. We note that none of the three leading modes of low-frequency SST variability show a large cooling over this period (Fig. S1). The warming in the SPNA from the mid-1990s to the early 2000s (Robson et al. 2016) is also well captured by the ADV index.
Further evidence of Atlantic decadal variability playing a role in the recent cooling of the SPNA is a strengthening of the subpolar gyre circulation and a weakening of the AMOC (Figs. 15c,d) (Smeed et al. 2018;Chafik et al. 2019), which are signatures of the ADV transitioning to a negative phase (the development between lag 0 to 18 years in Figs. 8 and 9). The recent Atlantic cold anomaly is thus consistent with ocean circulation changes associated with the decadal mode of North Atlantic SST variability identified here.
The important role of ocean circulation and associated heat transport changes in driving the recent cooling is in agreement with Robson et al. (2016) and Bryden et al. (2020). However, in addition to this ocean-driven cooling, the shifts to a positive NAO index in 2014 and 2015 lead to anomalous surface heat loss (Josey et al. 2018), which could explain why the ADVdriven cooling is less than that observed (note the different color scales in Figs. 15a and 15b). Our results nevertheless show that significant changes in ocean circulation preceded the cold anomaly and that surface heat fluxes simply enhanced it.
As the ADV index has been negative in recent years, and, hence, the SPNA in a cold state, it follows that in the absence of all other forcing we might now expect a transition to a positive ADV that would contribute to a warming of the SPNA over the next few years. A shift to a positive ADV will also-according to the mechanism presented here-be accompanied by a strengthening of the AMOC. Observations at 458N indeed suggest that the AMOC may already be increasing (Desbruyères et al. 2019).

Conclusions
There has recently been a large focus on identifying the mechanisms responsible for Atlantic multidecadal variability (AMV; e.g., Clement et al. 2015;O'Reilly et al. 2016;Drews and Greatbatch 2017;Wills et al. 2019;Zhang et al. 2019). However, decadal-scale variability embedded within the AMV has received less attention, despite being a prominent feature of observed North Atlantic temperature (Fig. 1) (Deser and Blackmon 1993;Häkkinen 2000;Czaja and Marshall 2001;Chafik et al. 2016;Nigam et al. 2018) and important for the climate of adjacent continents (Årthun et al. 2018;Robson et al. 2018). These decadal fluctuations in the North Atlantic FIG. 14. Schematic evolution of decadal ocean and atmosphere variability over the North Atlantic Ocean. The time line shows the processes acting to shift the ADV from a negative to a positive state and back again (see text for details). The map shows the SST anomalies during ADV maximum and the relevant atmospheric and oceanic circulation features (not to scale). Squiggly arrows indicate atmosphere-ocean heat fluxes. Red (blue) contours and letter H (L) denote high (low) sea level pressure anomalies.
Ocean are also a key source of skill in decadal climate predictions (Yeager and Robson 2017;Årthun et al. 2017).
In this study we have used low-frequency component analysis (Wills et al. 2019) to identify a decadal mode of North Atlantic SST variability [referred to as Atlantic decadal variability (ADV)]. The spatial pattern of the observed ADV resembles what has been referred to as the North Atlantic horseshoe pattern (Czaja and Frankignoul 2002), with a positive phase characterized by warm anomalies in the SPNA and in the eastern subtropical Atlantic and a cold anomaly in the western subtropical gyre. The ADV pattern also shows some similarity to the AMV pattern (e.g., Ting et al. 2011;Kim et al. 2018); this may be because, as shown in Wills et al. (2019), the traditional AMV/AMO definition mixes different time scales and mechanisms, including that of the ADV.
To determine the mechanisms underlying decadal SST variability we have used observations, large ensemble historical simulations (CESM-LE) and preindustrial control simulations (CMIP5). We find that the ADV is driven by large-scale atmosphere-ocean interaction that leads to anomalous SSTs both through concomitant air-sea heat fluxes and through delayed ocean circulation changes, the latter setting the time scale of variability (summarized in Fig. 14).
Coupled climate models are known to have mean-state biases in North Atlantic temperature and salinity, which could affect their representation of decadal variability (Menary et al. 2015b).
By using low-frequency component analysis, similar patterns of decadal SST variability can be identified across models and observations; the associated patterns of atmospheric and oceanic gyre circulation anomalies are fairly consistent across the different models. This provides confidence in the results and the proposed mechanism.
North Atlantic SST variability is determined by a wide range of processes acting on different temporal and spatial scales. The footprint of variability across the full range of time scales therefore needs to be considered when interpreting observed SST changes and trends. For example, the recent observed decadal cooling in the SPNA is not fully captured by the AMV index, as the subpolar cold anomaly is accompanied by a warm anomaly in the subtropics (Frajka- Williams et al. 2017). Rather, the cooling in the SPNA is consistent with the temporal development and spatial pattern of the ADV mode described here. This strongly suggests that the recent cold anomaly in the subpolar North Atlantic is, in part, a result of decadal SST variability, and that we might expect it to become less pronounced over the next few years. 727852). R.C.J.W was supported by the U.S. National Science Foundation (Grant AGS-1929775). H.L.J was supported by a RCN visiting researcher grant (310392). L.C. acknowledges support from the Swedish National Space Agency through the FiNNESS project (Dnr: 133/17). We thank the CESM Large Ensemble Community Project for making their data publicly available. We also acknowledge the World Climate Research Programme's Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modeling groups for producing and making available their model output.
Data availability statement. All data used in this study are freely available. HadISST data are available at http://hadobs. metoffice.com/. Data from the 20th Century Reanalysis are available at https://psl.noaa.gov/data/gridded/data.20thC_ ReanV2.html. The CESM-LE output is available via the Earth System Grid at https://www.earthsystemgrid.org. The CMIP5 data are available from the Earth System Grid Federation (ESGF) at esgf-node.llnl.gov/search/cmip5.