A delayed oscillator model for the quasi-periodic multidecadal variability of the NAO

Wavelet analysis of the annual North Atlantic Oscillation (NAO) index back to 1659 reveals a significant frequency band at about 60 years. Recent NAO decadal variations, including the increasing trend during 1960–1990 and decreasing trend since the mid-1990s, can be well explained by the approximate 60-year cycle. This quasi 60-year oscillation of the NAO is realistically reproduced in a long-term control simulation with version 4 of the Community Climate System Model, and the possible mechanisms are further investigated. The positive NAO forces the strengthening of the Atlantic meridional overturning circulation (AMOC) and induces a basin-wide uniform sea surface temperature (SST) warming that corresponds to the Atlantic multidecadal oscillation (AMO). The SST field exhibits a delayed response to the preceding enhanced AMOC, and shows a pattern similar to the North Atlantic tripole (NAT), with SST warming in the northern North Atlantic and cooling in the southern part. This SST pattern (negative NAT phase) may lead to an atmospheric response that resembles the negative NAO phase, and subsequently the oscillation proceeds, but in the opposite sense. Based on these mechanisms, a simple delayed oscillator model is established to explain the quasi-periodic multidecadal variability of the NAO. The magnitude of the NAO forcing of the AMOC/AMO and the time delay of the AMOC/AMO feedback are two key parameters of the delayed oscillator. For a given set of parameters, the quasi 60-year cycle of the NAO can be well predicted. This delayed oscillator model is useful for understanding of the oscillatory mechanism of the NAO, which has significant potential for decadal predictions as well as the interpretation of proxy data records.


Introduction
The North Atlantic Oscillation (NAO) is one of the dominant modes of atmospheric circulation variability in the Northern Hemisphere (NH), particularly in the North Atlantic region (Hurrell 1995;Hurrell et al. 2003, and references therein). It has been well documented in the literature that over interdecadal timescales the NAO has important impacts on regional and hemispheric climates in the NH (Hurrell 1995;Thompson and Wallace 2000;Hurrell et al. 2003;Li et al. 2013, and references therein). The NAO shows a remarkable upward trend over the second half of the twentieth century (Hurrell 1995;Thompson and Wallace 2000). This upward trend explains much of the observed warming trend over Eurasia and North America (Hurrell 1996), and has been linked to the interdecadal variations of Asian winter monsoon (Gong et al. 2001; Abstract Wavelet analysis of the annual North Atlantic Oscillation (NAO) index back to 1659 reveals a significant frequency band at about 60 years. Recent NAO decadal variations, including the increasing trend during  and decreasing trend since the mid-1990s, can be well explained by the approximate 60-year cycle. This quasi 60-year oscillation of the NAO is realistically reproduced in a long-term control simulation with version 4 of the Community Climate System Model, and the possible mechanisms are further investigated. The positive NAO forces the strengthening of the Atlantic meridional overturning circulation (AMOC) and induces a basin-wide uniform sea surface temperature (SST) warming that corresponds to the Atlantic multidecadal oscillation (AMO). The SST field exhibits a delayed response to the preceding enhanced AMOC, and shows a pattern similar to the North Atlantic tripole (NAT), with SST warming in the northern North Atlantic and cooling in the southern part. This SST pattern (negative NAT phase) may lead to an atmospheric response 1 3 et al. 2005). However, since the 1990s, the NAO has shown a significant decreasing trend (Cohen and Barlow 2005;Cohen et al. 2012). In particular, the winters of 2009/2010 and 2010/2011 were characterized by strong negative NAO phases (Taws et al. 2011), and the occurrence of negative NAO is suggested to be one of the main factors causing extremely cold temperatures in the NH . The weakening of the NAO in the past two decades also made a major contribution to the decrease in Arctic sea ice concentration (Deser and Teng 2008). Recently, Li et al. (2013) demonstrated that the NAO is connected to the NH mean surface temperature over multidecadal timescales, with the NAO leading by 15-20 years, and thus can act as a useful predictor. Therefore, a better understanding of the NAO decadal variability is helpful for the interpretation and future projections of decadal-scale climate variations.
Many studies have been conducted on the cause of the NAO decadal variability, but there is still little consensus. Several factors, such as greenhouse gas emissions and warming in tropical oceans, have been suggested to account for the NAO interdecadal variations. A number of modeling studies have revealed a change in the NAO toward a positive polarity in response to increasing greenhouse gas concentrations (Shindell et al. 1999;Fyfe et al. 1999), arguing that greenhouse gases were a contributing factor to the NAO upward trend during the 1960s to 1990s. Hoerling (2001) and Bader and Latif (2003) presented modeling evidence that the NAO increasing trend during the 1960s to 1990s was forced by a progressive warming of tropical sea surface temperatures (SST). Recently, Cohen and Barlow (2005) re-examined the various factors that have been related to the decadal variations in the NAO, and found that neither the warming in tropical oceans nor the increase of greenhouse gases could explain the NAO downward trend during the two most recent decades. Thus, the underlying dynamical mechanism of NAO decadal to multidecadal fluctuations remains elusive.
The North Atlantic Ocean features remarkable multidecadal variability and has received a considerable amount of attention since it shows the largest potential for decadal predictability (Latif et al. 2006a;Corti et al. 2012). The principle expression of the variability is the oscillation of North Atlantic SST between alternating warm and cold phases, which is commonly referred to as the Atlantic multidecadal oscillation, or AMO (Enfield et al. 2001). It has been identified in both observations and reconstructions that the AMO has a statistically significant spectral peak in the 50-70 year band (Schlesinger and Ramankutty 1994;Delworth and Mann 2000;Gray 2004;Tung and Zhou 2013). There is some consensus that the multidecadal variability in the North Atlantic SST is closely linked to the fluctuations in the Atlantic meridional overturning circulation (AMOC) (Delworth and Mann 2000;Latif et al. 2004;Knight et al. 2005;). Vellinga and Wu (2004) used coupled model simulations and fully diagnosed the mechanisms responsible for the AMOC fluctuations.
Previous studies have also suggested that the NAO significantly influences the oceanic variability in the North Atlantic, such as the strength of the ocean circulation and North Atlantic SSTs (Visbeck et al. 2003 and references therein). In particular, recent studies have shown that the NAO decadal variation can exert a persistent forcing on the AMOC, and thus leads the surface AMO (Latif et al. 2006a, b;Latif and Keenlyside 2011). Li et al. (2013) identified the lead time between the NAO and AMO (NAO leading by 15-20 years) using observations. In the present paper, we further show that there is a two-way interaction between the NAO and AMO; that is, the NAO has a positive forcing effect on the AMO, consistent with the previous findings, while the AMO in turn provides a delayed negative feedback on the NAO by inducing a meridional SST gradient pattern. These positive and negative feedback mechanisms act successively, leading to the quasi-periodic multidecadal variability of the NAO.
The remainder of the paper is organized as follows: The datasets and methodology are described in Sect. 2. Section 3 outlines the main spectral feature of the NAO multidecadal variability, its relation to the AMO, and the fundamental feedback mechanisms. In Sect. 4 we construct a theoretical model to explain the findings. Finally, conclusions and a discussion are presented in Sect. 5.

Data
The observed NAO index (hereafter referred to as NAOI O ) used in this study is defined as the difference in the normalized SLP zonally averaged over the North Atlantic sector from 80°W to 30°E between 35°N and 65°N (Li and Wang 2003;Li et al. 2013), derived from the Hadley SLP dataset (Allan and Ansell 2006) for the period 1850-2012. Because of the shortness of available instrumental records, we also use the reconstructed NAO index (hereafter referred to as NAOI R ) derived by Luterbacher et al. (1999)

3
Climate models provide a useful approach to analyze multidecadal variability. A long-term pre-industrial control simulation of version 4 of the Community Climate System Model (CCSM4) is used. CCSM4 is a fully coupled model of the earth's physical climate system, and a general description is given in Gent et al. (2011). Several recent studies have used the CCSM4 pre-industrial control simulations to investigate the Atlantic multidecadal variability and found that the model performs reasonably well in capturing the variability (Danabasoglu et al. 2012;Zanchettin et al. 2014). The control simulation was integrated over 1,300 model years, and a 300-year segment between model years 800 and 1,100 is used in the present analysis. The reason for this choice is to avoid the initial adjustment period because of the model spin-up process (Danabasoglu et al. 2012). The other reason is that the multidecadal oscillation of the NAO in the 300-year segment is regular and very similar to the reconstruction (Fig. 2), ensuring our analysis is relevant to the real world. The CCSM4 simulation data used for the analysis include variables at the airsea interface (e.g. SLP and SST) and variables in the deep ocean (e.g. sea water temperature and meridional overturning circulation streamfunction).

Statistical method
Principal Oscillation Pattern (POP) analysis (von Storch et al. 1988) was applied to investigate the dominant SST patterns of the oscillatory multidecadal variability. In contrast to Empirical Orthogonal Function (EOF) analysis, which considers only the spatial co-variability in a dataset, the POP method considers the full spatiotemporal structure. POPs are in general complex, with a real part (P R ) and imaginary part (P I ). The corresponding complex coefficient time series satisfy the standard damped harmonic oscillator equation, so that the evolution of the system in the POP space spanned by the real and imaginary part can be interpreted as the following oscillatory sequence of the spatial pattern: · · · P R → −P I → −P R → P I → P R · · ·. Prior to the POP analysis, band-pass filtering (Butterworth band-pass filter) was applied to annual mean SST anomalies in order to highlight the multidecadal variability. An EOF analysis was then performed on the filtered data, and only the leading 10 EOF modes were retained in the POP analysis. Note that the SST data were area-weighted prior to the analysis and no normalization was applied. Similar approaches to POP analysis have been adopted in other studies (e.g. Park and Latif 2010).
We apply a local wavelet spectrum analysis with a Morlet wavelet base (Torrence and Compo 1998) to test whether the time series show significant periodic oscillations. The statistical significance of the linear regression coefficient and correlation between two auto-correlated time series is accessed via a two-tailed Student's t test using the effective number of degrees of freedom (Pyper and Peterman 1998;Li et al. 2013). The effective number of degrees of freedom N eff is given by the following approximation: where N is the sample size and ρ XX (j) and ρ YY (j) are the autocorrelations of two sampled time series X and Y at time lag j.

Atmospheric general circulation model (AGCM)
We perform idealized numerical experiments to investigate the atmospheric circulation response to a given SST pattern. All the experiments have been designed and implemented on the Simplified Parameterizations, primitivE-Equation DYnamics (SPEEDY) model developed at the Abdus Salam International Centre for Theoretical Physics, Italy (Kucharski et al. 2013). SPEEDY is an intermediate complexity AGCM based on a spectral primitive-equation dynamical core in sigma coordinates, with semi-implicit treatment of gravity waves. The resolution used here is triangular 30 (T30, about 3.75° × 3.75° horizontal resolution), with 8 vertical levels. The model is computationally inexpensive, due to the simplified parameterizations, as well as the relatively low horizontal and vertical resolutions. Nevertheless, the simulation performance of the SPEEDY model compares very well with that of more complex AGCMs with respect to extratropical large-scale variability (Kucharski et al. 2013 and references therein). The control run was integrated for 15 years and forced with the climatological SST, and the last 10 years of the integrations were used to provide the basic annual mean state. The sensitivity experiment was integrated for 12 years and forced with the given SST anomaly pattern imposed on the climatological SST. A 10-member ensemble mean was constructed from the last 10 years of the integrations to reduce the uncertainties arising from different initial conditions. The atmospheric response to the SST forcing is defined as the difference between the ensemble mean for the sensitivity and control runs. (1) Throughout the whole period, the NAO fluctuates between positive and negative phases with time scales on the order of decades, and no significant trend is found. The correlation between the NAOI O and NAOI R over the 1850-2000 overlap period is high (r = 0.86, significant at the 95 % confidence level based on the 21-year running mean data), indicating the NAO index reconstructed by Luterbacher et al. (1999) is reasonably reliable. This is also demonstrated by previous studies on the comparison among different NAO reconstructions (Luterbacher et al. 1999;Cook 2003). We further applied a Morlet wavelet analysis to the NAOI C . As can clearly be seen in Fig. 1b, the wavelet spectrum for the full period shows significant energy over a frequency band from 50 to 70 years, suggesting a strong quasi-periodicity of about 60 years in the annual NAO index. In particular, the quasi-periodic fluctuations with periods of around 60 years are most pronounced from the mid-1700s to the end of the record. A quasi-periodic cycle of about 50-70 years in the NAO has also been reported in previous studies based on other NAO reconstructions (Cook et al. 1998;Olsen et al. 2012;Dieppois et al. 2013). Thus, our study further identifies the quasiperiodic component in terms of the combination of the NAO reconstruction and observation. The 50-70 year band-pass filtered NAOI C is strongly correlated with the 21-year running mean during the whole period (r = 0.79, significant at the 95 % confidence level), suggesting that a significant proportion of the multidecadal variability can be interpreted by the approximate 60 year cycle. In particular, the recent NAO interdecadal variations, including the upward trend during 1960-1990 and the downward trend since the mid-1990s, are also well explained. Therefore, the quasi-periodic oscillation with time scales of 50-70 years seems to be a real and robust aspect of the NAO variability, and is important for understanding the NAO interdecadal variations. Exploring the possible factors that contribute to the NAO quasi-periodicity is hindered by relatively short global climate records and insufficient ocean data. Thus, long-term simulations from state-of-the-art climate models provide an alternative tool to investigate the mechanisms of the multidecadal variability. Here, we use a pre-industrial control simulation with the CCSM4 climate model. The simulated annual NAO index over the simulation period is shown in Fig. 2a. The NAO index in CCSM4 is defined similarly to that in the observation, and is highly correlated with the principal component (PC)-based index of the NAO, defined as the leading EOF of annual SLP anomalies over the North Atlantic sector (20°-80°N, 90°W-30°E), yielding a correlation The dotted area indicates the 95 % significance level against red noise. The dashed line in the right-hand panel is the significance for the global spectrum. c As in (b), but for the annual AMO index coefficient of 0.92 for the entire simulation period. Consistent with the NAOI C , the wavelet analysis of the simulated annual NAO reveals significant power at the frequency band around 60 years ( Fig. 2b), prevailing through most of the simulation. The wavelet spectrum also shows some other peaks at about 4-8 and 16 years, but they are not statistically significant. The band-pass (50-70 years) filtered version of the simulated annual NAO index is significantly correlated with the 21-year running mean (r = 0.74), further confirming the significant variability with periods of about 60 years. The 50-70 year quasiperiodicity of the NAO is considerably well simulated by the CCSM4 model, suggesting that we can use the simulations to further interpret the physical mechanisms that control the NAO multidecadal variability. Figure 2a also displays the 21-year running mean of the AMO index in the CCSM4 simulations. We define the simulated AMO index as the area-weighted average of the annual SST anomalies over the extratropical North Atlantic north of 20°N (20°-65°N, 75°-7.5°W). Different from the usual definitions of the AMO index, the tropical North Atlantic (0°-20°N) has been omitted here. Previous studies have suggested that the tropical North Atlantic SST is significantly influenced by the well-known ENSO teleconnections (van Oldenborgh et al. 2009) and shows relatively weak multidecadal variability (Trenberth et al. 2007;Delworth et al. 2007). Meanwhile, in both observational and simulation studies, the AMO signal is found to be concentrated in the high latitudes of the North Atlantic (Sutton and Hodson 2005;Knight et al. 2005). The effect of the different domain size in the AMO index is small; the SST area average north of the equator is highly correlated to the average from 20°N to 65°N (r = 0.91) for the 21-year running mean data. The AMO appears to fluctuate between warming and cooling phases on multidecadal timescales. Wavelet analysis of the AMO index shows a band of significant variability near periods of around 50-70 years persisting through most of the simulations (Fig. 2c), indicating that the model produces a number of cycles of a repeating mode of North Atlantic SST variability with a 50-70 year time scale. This time scale is consistent with that identified in observation and reconstruction records (Schlesinger and Ramankutty 1994;Delworth and Mann 2000;Gray 2004), indicating that the spectral feature of the AMO is realistically reproduced in the CCSM4 simulations. The quasi-periodic multidecadal variability in the AMO has also been reported by previous studies with CCSM4 (Danabasoglu et al. 2012) as well as other climate model simulations (Delworth and Mann 2000;Knight et al. 2005). Moreover, in this study, the time scale of simulated AMO is basically identical to that found in the NAO index, implying a potential linkage between the AMO and NAO multidecadal variability.

Possible mechanisms involved
In order to help understand the physical mechanisms involved, we applied POP analysis to detect the dominant SST patterns in association with the North Atlantic quasi-periodic multidecadal variability in CCSM4, and a 50-70 year band-pass filtering was performed on the annual mean North Atlantic SST data. Note that using a broader filter (e.g., a 30-90 year band-pass filtering) gives qualitatively similar results. The leading POP mode (Fig. 3), accounting for about 60 % of the filtered variability, is oscillatory, and characterized by a period of 60 years. This indicates that the 50-70 year quasi-periodic variability is mainly explained by POP1. The real-and imaginary-part PCs (Fig. 3c) show, as expected from the definition of POPs, a phase shift of about 90° during most of the analyzed period. The real-part pattern of POP1 ( Fig. 3a) features positive SST anomalies covering most of the North Atlantic region, similar to the basin-wide uniform SST pattern associated with the AMO (Enfield et al. 2001;Knight et al. 2005). Moreover, the real-part PC1 is strongly correlated with the 21-year running mean AMO index (r = 0.82, significant at the 95 % confidence level). Therefore, the real-part POP1 in fact represents the AMO.
The imaginary part (Fig. 3b) exhibits a strong meridional gradient of SST across the North Atlantic, with cooling in the subpolar North Atlantic and warming in midlatitude North Atlantic regions. Moreover, the imaginary part is reminiscent of the well-known North Atlantic SST tripole (Visbeck et al. 2003, and references therein;Deser et al. 2010). A simple North Atlantic tripole (NAT) index is therefore defined as the difference between the sum of averaged SST in two negative-value boxes (55°-65°N, 60°-20°W and 15°-30°N, 45°-15°W) and averaged SST in the positive-value box (30°-50°N, 80°-50°W). This definition is similar to Wu et al. (2009) but the locations of the three boxes are slightly shifted northward in the simulations. As expected, the 21-year running mean NAT index is highly correlated with the POP1 imaginary PC (r = −0.81, significant at the 95 % confidence level), indicating the POP1 imaginary part corresponds closely to the NAT pattern over multidecadal timescales. Park and Latif (2010) conducted POP analyses on the whole NH SST (north of 20°N) from a simulation by the Kiel Climate Model. Their POP1 is pronounced over the Pacific basin and corresponds to Pacific decadal variability, and their POP2 patterns resemble our POP1 patterns (AMO and NAT) in the North Atlantic. In this study, we focus on the possible influence of the AMO and NAT patterns on NAO quasi-periodic multidecadal variability, which has not been well addressed in the literature. The cross-correlation of the real-and imaginary-part PCs of POP1 (AMO and NAT indices) shows a quadrature relation 1 3 ( Fig. 4), with the NAT leading the AMO by 15-20 years and negative NAT lagging the AMO by 15-20 years. This lead-lag correlation supports the oscillatory sequence of the real and imaginary parts of the POP, which is as follows: · · · NAT → AMO → −NAT → −AMO → NAT · · ·.
Since the SST can influence the atmospheric circulation through changes in heat fluxes across the surface, we investigate the correlation patterns (maps of local correlation coefficients) of the North Atlantic SLP with the AMO and NAT in the CCSM4 simulations. Figure 5a shows the correlations with the AMO, exhibiting a basin-wide homogeneous pattern dominated by negative values. The negative correlations are most pronounced over Western Europe and the adjacent ocean. The AGCM simulated SLP response to the AMO SST anomalies also shows a basin-wide uniform pattern dominated by negative SLP anomalies (Fig. 5b), largely due to the strong surface heating induced by SST warming. This is consistent with the SLP patterns shown in previous observational and modeling studies focusing on the climatic impacts of the AMO (Delworth and Mann 2000;Sutton and Hodson 2005;Li and Bates 2007). In contrast, as shown in Fig. 5c, the SLP pattern in association with the NAT variability features a meridional dipole with significantly negative correlations around Greenland and significantly positive correlations in the region of 30°-50°N. This pattern strongly projects onto the NAO. Moreover, the AGCM simulated response to the NAT SST anomalies is also evidently NAO-like (Fig. 5d).
We further examine the lead-lag correlations between NAO, AMO and NAT in the CCSM4 simulations. The correlation curves between NAO and AMO based on decadally  (Fig. 6b) peak at lag zero, and a significant positive correlation (0.61 for unfiltered data and 0.72 for 21-year running means) is found, indicating that the NAT is largely in phase with the NAO. Based on the lead-lag correlations and the simulated SLP responses to AMO and NAT SST anomalies, we hypothesize three possible mechanisms involved in the oscillatory sequence of NAT and AMO. First, the NAT has a direct effect on the NAO; second, the NAO exerts some forcing effects on the AMO; and third, the AMO in turn provides some negative feedback on the NAT. These positive and negative feedback mechanisms may act successively, leading to the quasi-periodic multidecadal variability. In the following analysis, we provide more details about the physical processes of these mechanisms using the CCSM4 simulations.
Previous observation and modeling studies have provided evidence that the NAT is in close relationship with the NAO and that there is a positive feedback between them (Wu and Gordon 2002;Czaja and Frankignoul 2002;Kushnir et al. 2002, and references therein;Wu and Rodwell 2004;Gastineau et al. 2013). The gradient SST anomalies associated with the NAT induce changes in the position and strength of the storm tracks through thermal forcing, and the transient eddy-forcing associated with the storm tracks tends to excite and maintain the NAO-like response. Figure 7a shows the upward surface heat flux anomalies in association with the NAT. The heat flux anomaly has a similar pattern to the NAT, reflecting that the atmosphere acts to damp the SST anomalies over multidecadal timescales (Gulev et al. 2013). As a result, the low-level air temperature is affected and shows anomalous cooling and warming in the subpolar and subtropical basins, respectively (Fig. 7b), leading to changes of the baroclinicity in the lower troposphere and the synoptic perturbation growth. As shown by the maximum Eady growth rate (defined by 0.31 f|∂u/∂z|N −1 , where f is the Coriolis parameter, ∂u/∂z is the zonal wind vertical shear, and N is the Brunt-Vӓisӓlӓ frequency) at 850-hPa (Fig. 7c) and the 500-hPa geopotential height standard deviation (Fig. 7d), the enhanced meridional SST gradient of NAT further contributes to the increase of the storm track intensity and shifts the storm track northward, leading to a positive NAO phase. Therefore, the physical processes responsible for the strong coherence between the NAO and NAT over multidecadal timescales are further identified in the CCSM4 simulations.
There is substantial modeling evidence that NAO-related wind anomalies can lead to multidecadal fluctuations of the AMOC, which in turn produce the SST signatures of the AMO (Visbeck et al. 1998;Eden and Jung 2001;Visbeck et al. 2003 and references therein). Figure 8a shows the long-term mean streamfunction of AMOC for the CCSM4 simulation. The simulated position of maximum overturning (about 15 Sv) lies at the depth of 1,000-1,200 m, and 40°-50°N, which is basically consistent with the hydrographic data (Zhang and Wang 2013). Figure 8b shows the lagged regression of the streamfunction of AMOC onto the normalized NAO index when the NAO leads by 15 years. The regression pattern shows significantly positive anomalies in the North Atlantic Ocean, and largely resembles the long-term mean structure of the AMOC. This indicates that the AMOC tends to strengthen (weaken) about 1-2 decades after the NAO phase turns positive (negative). The enhanced AMOC leads to stronger northward heat transport in the upper ocean, and as a result, causes anomalously warm SST anomalies in the North Atlantic (Knight et al. 2005;Zhang and Wang 2013). An AMOC strength index is defined as the maximum of the annual mean AMOC streamfunction at 40°N, similar to a previous study (Zhang 2008). The AMOC varies strongly in phase with the AMO over multidecadal timescales (Fig. 8c) and is significantly correlated with the 15-year-ahead NAO (r = 0.23 for the original data and r = 0.54 for the decadally smoothed). Therefore, in accordance with previous studies (Latif et al. 2006a, b;Park and Latif 2010;Latif and Keenlyside 2011), the forcing effect of the NAO on the AMO through AMOC changes may be a key mechanism of the North Atlantic multidecadal variability.
Previous modeling studies have suggested that the key processes associated with the AMOC multidecadal fluctuations mainly consist of the advection induced by the overturning and the slow adjustment of the ocean circulation (Huang and Chou 1994;Greatbatch and Zhang 1995;Drbohlav and Jin 1998;Lee and Wang 2010;Zhang and Wang 2013;Ba et al. 2013;Huang et al. 2014). As shown in Fig. 9a, the warm (cold) phase of the AMO is associated with a uniform surface warming (cooling) and a subsurface cooling (warming) in the North Atlantic Ocean. Since the warm phase of the AMO coincides with the strengthening of the AMOC (Fig. 8c), the warming upper ocean water continues to be advected northward to the downwelling region of the AMOC and then sinks, while the cooling subsurface ocean water continues to be moved southward. the OHT shows maximum values around 40°N. This may lead to a heat convergence (divergence) north (south) of these latitudes, inducing an ocean column warming (cooling) tendency over the entire upper ocean down to about 2,000 m. As a result, the enhanced AMOC is expected to be decelerated due to this developing positive meridional temperature gradient anomaly. Thus, the OHT associated with the enhanced AMOC acts as a negative feedback to prevent the AMOC from further strengthening by weakening the meridional temperature gradient after a delayed time of ocean adjustment.
The ocean temperature signal following the AMO is illustrated by a correlation of temperature with the normalized AMO index as a function of time lag (Fig. 9c-e). The positive correlations are at first located in the upper North Atlantic and then propagate into the subpolar region, expanding downward; the negative correlations are shifted southward. This is consistent with the clockwise heat transport pathway of the OHT and further indicates that the ocean temperature is progressively modulated by the OHT associated with AMOC intensification. Previous studies have also shown a similar delayed response of ocean temperature to the preceding enhanced AMOC in both ocean (Greatbatch and Zhang 1995;Drbohlav and Jin 1998) and air-sea coupled (Ba et al. 2013) models. Moreover, the lagged correlation pattern at lag 18 is strikingly similar to the ocean temperature pattern associated with the NAT variability (Fig. 9f), implying that the AMOC/AMO may provide a delayed effect on the NAT, and the NAT index is negatively correlated with the 18-year-ahead AMO To further identify this delayed negative feedback, we provide the physical explanation for the time delay. The time delay can be qualitatively interpreted as the time it takes for transport along the OHT pathway. Thus, the time delay corresponds to the ratio L V (L and V are the pathway distance and velocity, respectively). We further separate L and V into vertical (L z and V z ) and meridional (L y and V y ) components. Due to the orthogonality of the two components, we have L V ≈ L y V y ≈ L z V z . V y is the meridional velocity of the ocean current, and it can be approximated by Ψ W ×D according to the formula of the AMOC streamfunction (Ψ, D and W are the maximum of the overturning streamfunction, the depth of the maximum, and the east-west width of the North Atlantic Ocean, respectively). Therefore, the time delay can be estimated by the ratio L y ×W ×D Ψ . As indicated in Fig. 9c-e, the OHT shifts the warming center northward by about 15 degrees of latitude (from 45°N to 60°N), so we set L y = 2 × 10 6 m. The North Atlantic Ocean comprises 1/6 of the latitudinal circumference, corresponding to W = 5 × 10 6 m. A previous multimodel analysis suggested that the position of the maximum AMOC transport occurs at 500-1,500 m depth and the strength of the AMOC is about 20 Sv (Zhang and Wang 2013); thus, we set D = 10 3 m and Ψ = 20 Sv. Dividing L y × W × D by Ψ, the time delay is about 16 years. The scale of the time delay is consistent with the overturning time on the order of one decade or more, which represents the timescale of overturning advection (Winton and Sarachik 1993;Ou 2011). Moreover, the time delay estimated by the ratio L y ×W ×D Ψ is very close to the time lag of the NAT relative to the AMO (shown in Fig. 4), further suggesting that the advection plays a key role in the SST evolution from positive AMO to negative NAT. As mentioned above, the negative NAT corresponds to the negative NAO phase in the atmosphere, and consequently, the negative NAO is expected to force the AMOC perturbation to its negative phase and the oscillation proceeds as discussed above, but in the opposite sense. Therefore, the delayed negative feedback of the AMOC/AMO on the NAT/NAO through the ocean advection process may be another important mechanism for the quasi-periodic multidecadal variability of the NAO.

A delayed oscillator model for the quasi-periodic multidecadal variability of the NAO
Based on the analyses in Sect. 3, we can summarize three key mechanisms involved in the quasi-periodic North Atlantic multidecadal variability (Fig. 10): the coherence between the NAT and NAO; the forcing of the AMOC/ AMO by the NAO; and the delayed negative feedback of the AMOC/AMO. These three mechanisms can be represented by the following three equations, respectively: where NAO, NAT and AMO are normalized dimensionless index series and are functions of time t (in years) only. Equation 3 is similar to the Hasselmann climate model (Hasselmann 1976;Thompson et al. 2009) used by Li et al.
(2) NAO(t) ≈ NAT(t), (2013) to model the NH multidecadal variability in response to NAO variations. C is the estimated effective heat capacity, in the range of 13 to 35 W year m −2 K −1 , and the linear damping coefficient β is estimated in the range from 0.6 to 1.2 K (W m −2 ) −1 (Knutti and Hegerl 2008). According to previous studies, the effective heat capacity actually represents the thermal inertia of the climate system, which is broadly defined to include oceanic processes of long time scales (e.g. ocean heat uptake and transport) ). Here, β and C are set to 0.9 K (W m −2 ) −1 and 25 W year m −2 K −1 , respectively. In fact, the results are not sensitive to the choice of β or C provided that the values fall within ranges that are physically reasonable.
The parameter α represents the magnitude of the forcing effect of the NAO on the AMO and its units are K −1 W m −2 . Using the CCSM4 simulation data, the parameter α can be determined empirically. The 50-70 year band-pass filtered versions of the CCSM4-simulated AMO and NAO indices are used and both indices are normalized to have unit standard deviation, such that the indices are dimensionless and the effects of other time scales of variability are minimized. The Hasselmann climate model (Eq. 3) was initialized with the AMO index value at simulation year 1 and integrated with the time series of the NAO forcing term, and the output for simulation years 850-1,100 was used in the analyses. Figure 11a shows the root-meansquare error (RMSE) ratios between the CCSM4-simulated AMO and the output from the Hasselmann climate model as a function of parameter α. The optimal parameter α is determined at 2.4 K −1 W m −2 when the RMSE ratio is minimized (below 0.5). As shown in Fig. 11b, the CCSM4simulated AMO is largely in phase with the output from (e) (f) Fig. 9 a Correlations of the North Atlantic zonal mean (75°-7.5°W) ocean temperature with the AMO index in the CCSM4 simulations as a function of depth and latitude. b Regression of the OHT (10 14 W, see text for the definition) onto the normalized AMO index. c-e As in (a) but for lagged correlations. The lag (in years) is positive when AMO leads. f Correlations between the zonal mean ocean tempera-ture and the imaginary-part PC of POP1 (NAT index). For convenience of comparison with (e), the sign of the NAT index is reversed. In (a-f), a 21-year running mean was applied to the data prior to the correlation/regression analysis, and dots indicate correlations/regressions significant at the 95 % confidence level the Hasselmann climate model (r = 0.91, significant at the 95 % confidence level), and the two time series also show comparable amplitudes. The parameter τ in Eq. 4 denotes the time lag in years between the AMO and NAT and is on the order of one decade or more, consistent with the overturning timescale.
A simple delayed oscillator model for the NAT and NAO based on Eqs. 2-4 can be further derived as follows: Note that a similar delayed oscillator model can also be derived for the AMO. Analogous to the normal mode analysis of the linear delayed oscillator model to help understand the quasi-periodic behavior of the ENSO phenomenon (Battisti and Hirst 1989;Power 2010), we examine the natural frequencies of Eq. 5. Substitution of NAO(t) = NAO(0)e σt into Eq. 5 gives σ = − 1 βC − α C e −σ τ , where NAO(0) represents the initial value and σ = σ R + iσ I is the complex natural frequency. Equating the real and imaginary parts respectively gives: Since β and C are set to 0.9 K (W m −2 ) −1 and 25 W year m −2 K −1 , the solutions of Eqs. 6-7 or natural frequencies are functions of parameters α and τ. Figure 12 shows the solutions of Eqs. 6-7 (σ R and σ I ) for α in the range of 0.1-4 K −1 W m −2 and τ ranging from 1 to 30 years. The growth rate and period of oscillation can be determined easily by 1/σ R and 2π/σ I , respectively. A zero value in Fig. 12b indicates that the low-frequency mode is non-oscillatory. As shown in Fig. 12, the structure of the growth rate and period of oscillation is quite rich, and includes non-oscillatory growth and decay, as well as growing, decaying and stable oscillatory behavior. In general, solutions for parameter combinations of large values (large α and τ) often appear to be low-frequency and quasi-stable, while solutions for combinations of small values (small α and τ) often show high decay rates. In the above analyses, we found a parameter combination (α = 2.4 K −1 W m −2 , τ = 18 years) from the CCSM4 simulations. Figure 12 indicates that this combination corresponds to a weakly The enhanced AMOC continues to affect the heat transport, and due to slow ocean adjustment, the North Atlantic Ocean shows a delayed response (after about 18 years) to the preceding enhanced AMOC with an SST pattern that resembles the negative NAT phase. The negative NAT phase coincides with the negative NAO phase in the atmosphere, and thus the cycle proceeds, but in the opposite sense. Blue (black) text indicates oceanic (atmospheric) phenomena (a) (b) Fig. 11 Empirical determination of the parameter α for the NAO term in the Hasselmann climate model (Eq. 3). a Ratios of rootmean-square errors (RMSE) between the CCSM4-simulated AMO (real-part PC of POP1) and the output from the Hasselmann climate model to the root-mean-square of CCSM4-simulated AMO as a function of coefficient α. The optimal coefficient α is determined at 2.4 when the RMSE ratio is minimized (below 0.5). b CCSM4-simulated AMO (blue) and the output from the Hasselmann climate model (AMO NAO , red) for CCSM4 simulation years 850-1100 1 3 damped oscillation with a period of 65 years and a decay time (e-folding scale) of 180 years. This period is quite close to the multidecadal period of the NAO, indicating that, given realistic parameter values, this simple delayed oscillator model is capable of reproducing the quasi-periodic multidecadal variability of the NAO. Figure 13 shows some examples of the numerical solution of Eq. 5. Numerical solutions for τ = 18 years, and α = 2.4, 2 and 3 K −1 W m −2 are presented in Fig. 13a. They correspond to quasi-stable oscillatory (α = 2.4 K −1 W m −2 ), oscillatory decay (α = 2 K −1 W m −2 ), and oscillatory growth (α = 3 K −1 W m −2 ) behaviors, respectively. Consistent with the normal mode analysis, α = 2.4 K −1 W m −2 gives the best-fit to the oscillatory behavior of the NAO. Because the parameter α influences the magnitude of the delayed feedback term, the solution appears to grow (decay) if the parameter α increases (decreases). This indicates that the oscillation can be sustained when the parameter α is sufficiently large. Figure 13b shows numerical solutions for α = 2.4 K −1 W m −2 , and τ = 18, 10 and 30 years. The solutions exhibit quasi-stable oscillatory (τ = 18 years), oscillatory decay (τ = 10 years), and oscillatory growth (τ = 30 years) behaviors. As shown in Fig. 13b, the oscillatory period becomes long (short) as the delay time increases (decreases). Meanwhile, the solution oscillates only when the delay time is long enough. A long delay time allows alternating actions of amplification and stabilization and thus gives rise to an oscillation. Therefore, the time delay in the negative feedback of the AMO on the NAO is a key element of the delayed oscillator.
In order to further understand the important role of the delayed negative feedback, we rewrite Eq. 2, adding a stochastic forcing term as where F denotes the strength of the feedback and W(t) is white noise with zero mean and unit variance. The time series of NAO can be obtained from the numerical solutions of Eqs. 3, 4 and 8. The equations are integrated for 200 years with white noise forcing and the parameters α and τ are set to 2.4 K −1 W m −2 and 18 years, respectively. Figure 14 shows the numerical solutions of NAO for various values of the parameter F (F = 1 and 0.5). When F = 1, the time series of NAO shows very strong multidecadal variability with the oscillatory period around 60 years, and the variance percentage explained by the multidecadal oscillation (50-70 year band-pass filtered) is about 60 %. When F = 0.5, the NAO multidecadal variability is also evident but becomes weak, and the variance percentage of the multidecadal oscillation is reduced to 20 %. We also consider the limit case F = 0, which corresponds to no feedback from the ocean to the atmosphere, and the NAO is simply a white noise. Previous theoretical studies have also suggested that the feedback of ocean onto the NAO is critical for the NAO spectrum peak at the decadal frequency band (Marshall et al. 2001;Czaja and Marshall 2000). Therefore, our results are consistent with the previous findings and further highlight the importance of ocean-atmosphere coupling in the NAO multidecadal variability.

Summary and discussion
This study investigated the quasi-periodicity of the NAO over multidecadal timescales. The annual NAO index derived from the reconstructions and observations (NAOI C ) for the past 300-400 years shows a significant frequency band with a period of close to 60 years. The recent NAO interdecadal variations, including the increasing trend during  and the decreasing trend since the mid-1990s, are well explained by the approximate 60-year cycle in the NAO. This quasi 60-year oscillation was realistically reproduced in a long-term control simulation with the CCSM4 model. The simulated annual NAO and AMO indices both show a significant spectrum peak for periods around 60 years.
The possible mechanisms involved were further studied using the CCSM4 simulations. POP analysis revealed two dominant SST patterns over the North Atlantic in association with the 50-70 year cycle, and this cycle is manifested by the alternating appearance of these two SST patterns. One is basin-wide and uniform, corresponding to the AMO, while the other shows a strong meridional gradient and resembles the NAT pattern. It was found that the NAT covaries strongly with the NAO over multidecadal timescales. Further analyses identified two possible mechanisms responsible for the quasi 60-year oscillation: forcing of the AMO by the NAO, and delayed negative feedback exerted by the AMO. A positive NAO forces enhancement of the AMOC, leading to a positive AMO, and the enhanced AMOC continues to affect ocean heat transport, inducing an SST pattern with a time delay of the overturning timescale that resembles the negative NAT phase. The negative NAT phase coincides with the negative NAO phase in the atmosphere, and thus the cycle proceeds, but in the opposite sense.
Based on the above mechanisms, a simple delayed oscillator model has been constructed to explain the quasi-periodic multidecadal variability of the NAO. The magnitude of the NAO forcing effect (α) and the time delay of the AMO feedback (τ) are two key parameters of the delayed oscillator. Normal mode analysis was applied to assess the importance of the parameter values in determining the behavior of the model solutions. We showed that solutions can be oscillatory or non-oscillatory, and that they can grow, remain stable or decay, depending on the parameter values chosen. Our best estimates of the parameters of the delayed oscillator model are τ = 18 years and α = 2.4 K −1 W m −2 . A solution with these values is presented in Fig. 13, where the NAO shows a strong oscillatory behavior (with a period of 65 years and damping time of 180 years). This behavior is consistent with multidecadal variability of the NAO. Thus, this delayed oscillator model can be a useful tool to help understand the underlying physics of the multidecadal variability of the NAO. It is also noted that one third of the multidecadal variance in the reconstructed NAO is left unexplained by the quasi-60 year cycle, indicating that the delayed oscillator is not the only mechanism responsible for the NAO multidecadal variability.
The mechanisms proposed in this study are based mainly on the results from the CCSM4 model. There is some consensus in observational studies that the NAO leads to the NAT nearly instantaneously through surface heat flux changes (Wu et al. 2009;Deser et al. 2010) and that the AMO follows the NAO with a 15-20 year delay through AMOC changes (Latif et al. 2006a, b;Latif and Keenlyside 2011;Li et al. 2013). These two mechanisms may provide an alternative interpretation for the oscillatory sequence and could partly explain the lead-lag relationships found 1 3 in the present study. However, these two mechanisms cannot explain the quasi-60 year period of the NAO without a feedback from the ocean onto the NAO. Previous studies have suggested that the oceanic feedback is difficult to detect in the observations and model dependent. For example, Scaife et al. (2009) have shown in a multi-model study that the late twentieth century NAO trend cannot be reproduced in models forced by the observed SSTs. This indicates either that the SSTs may be a weak factor in the NAO variability or that these models may show much too weak impact of SSTs on the NAO. The present analysis with the CCSM4 model provides some modeling evidence for the NAT feedback on the NAO, which is critical for reproducing the NAO multidecadal variability. Since the performance of a numerical climate model is always limited, a limitation of our study is that the results are obtained only by one single model. Thus, a multi-model comparison of the NAO multidecadal variability for a better understanding of how robust the mechanisms are among different climate models warrants further investigation.
Previous studies have suggested several potential sources of predictability for the NAO. For instance, one potential source of predictability comes from the downward control exerted by the stratospheric polar vortex on the tropospheric NAO (Baldwin et al. 2003). It has also been shown that the autumnal snow extent over Eurasia significantly influences the winter NAO, with an above (below) normal October snow extent leading to a negative (positive) phase of the winter NAO (Cohen 2003). These sources of predictability mainly affect the NAO variability at timescales ranging from intraseasonal to interannual. The present study has implications for predicting the NAO over interdecadal timescales. Assuming that the approximate 60-year cycle in the NAO continues, we should be entering a negative phase of the cycle, as suggested by Fig. 1a. This may correspond to more negative NAO events during the next decade. Considering that the negative NAO usually leads to cooling over most of the NH continent (Hurrell 1995(Hurrell , 1996, the anthropogenic warming effect may be partially offset. A limitation in this study is that we do not take into account the impact of "weather noise" in the delayed oscillator model. Previous studies have suggested that such weather noise is needed to excite and energize the North Atlantic multidecadal variability (Delworth and Greatbatch 2000). In fact, within a linear framework, the noise simply deflects the solutions, and in this sense, the solutions presented here can be regarded as the ensemble mean of a large ensemble driven by noise. The ensemble mean is the predictable part of the multidecadal variability, and thus it is of special interest here.