1 Introduction

The Atlantic region exhibits pronounced climate variability on a wide range of timescales, from seasonal to centennial and beyond (Marshall et al. 2001a). A lot of research during the last years focused on the so-called Atlantic Multidecadal Oscillation AMO (Folland et al. 1986; Knight et al. 2005, 2006), which can be identified in proxy data and instrumental observations, and in climate models (Delworth and Mann 2000). The AMO is defined as the leading statistical (Empirical Orthogonal Function, EOF) mode of long-term North Atlantic (NA)-SST variability, with widespread climatic impacts, for example on the surface climate over the land areas adjacent to the NA. Here, we use the more general term Atlantic Multidecadal Variability (AMV) to describe the multidecadal variability (~ 30–80 years) in the Atlantic, as the timescale and spatial pattern vary considerably among climate models and because it is not clear, whether the multidecadal variability is truly oscillatory or not. Booth et al. (2012) argue that a major part of the multidecadal variance in detrended 1860–2005 NA sea surface temperatures can be explained by aerosol forcing, which would speak against an oscillation. The term AMV used here would encompass the AMO, if the latter were an oscillation.

Other modes of low-frequency NA-climate variability have also been documented in observations and climate models such as the tropical Atlantic Meridional Mode (~ 10–15 years, Chiang and Vimont 2004), which influences, for example, the rainfall over Northeast Brazil and Northwest Africa. The topic of this study is the Atlantic decadal-to-bidecadal variability (~ 10–20 years, ADV hereafter). ADV has been described from proxy data (Sicre et al. 2008; Chylek et al. 2011; Cronin et al. 2014) as well as instrumental observations (Deser and Blackmon 1993; Sutton and Allen 1997; Divine and Dick 2006; Deser et al. 2010; Vianna and Menezes 2013; Årthun et al. 2021). Here we describe the ADV that is simulated in a climate model, located in the extratropical NA and clearly distinguishable from the AMV. Significant ADV has also been found in a variety of other climate models (Grötzner et al. 1998; Álvarez-García et al. 2008; Frankcombe et al. 2010; Escudier et al. 2013; Ortega et al. 2015; Årthun et al. 2021).

The two major ocean-circulation systems in the NA are the subpolar gyre (SPG) and the Meridional Overturning Circulation (MOC). In a number of climate models, the SPG and the MOC play important roles in the AMV. The SPG and MOC also influence each other. The SPG affects the MOC by influencing the formation of North Atlantic Deep Water (Delworth et al. 1993; Lohmann et al. 2009). The MOC drives SPG variability through changing the ocean-heat transport into the subpolar NA (Zhang 2008). In comparison to the AMV, the role of the ocean circulation in the ADV that is the topic of this study is less studied, because the ADV appears to be overwhelmed by the more pronounced AMV (Nigam et al. 2018).

Understanding the ADV is important among others to enhance multiyear predictions of the SST over the extratropical NA and related surface-climate variability over the surrounding land areas (Latif et al. 2006; Yeager and Robson 2017; Reintges et al. 2020). The roles of local atmospheric forcing and ocean dynamics in driving the extratropical NA-SST variability are still poorly understood. Consistent with the stochastic climate-model theory (Hasselmann 1976), observed short-term seasonal-to-interannual SST variations in the midlatitude NA can be largely attributed to local surface-heat flux forcing, which is integrated by the oceanic mixed layer (Frankignoul 1985; Cayan 1992). On decadal and longer timescales, the ocean dynamics are thought to play an important role in the generation of the extratropical NA-SST anomalies (Bjerknes 1964), with the surface-heat fluxes acting as a damping on the SST anomalies (Gulev et al. 2013). On these timescales, the extratropical SST variability is significantly driven through variations in heat-transport convergence, as inferred from climate models and ocean reanalysis (e.g., Eden and Willebrand 2001; Williams et al. 2014; Ba et al. 2014; Delworth et al. 2017). Finally, there is an ongoing debate about the sensitivity of the atmosphere to extratropical SST anomalies and the importance of dynamical ocean–atmosphere coupling.

The mechanisms driving ADV are under debate. By using climate models and observations, Årthun et al. (2021) propose that the decadal SST variability in the NA is driven by the atmosphere through air–sea heat fluxes and delayed ocean-circulation changes, which is consistent with some previous studies (Nigam et al. 2018; Häkkinen 2000). Some studies linked the interdecadal MOC variability (~ 20–30 years) to westward propagating baroclinic Rossby waves (e.g., Frankcombe et al. 2010; Sévellec and Fedorov 2013; Muir and Fedorov 2017). Moreover, variability linked to the Arctic sea-ice and the East Greenland Current may contribute to the ADV (Escudier et al. 2013; Ortega et al. 2015). Finally, external forcing as well could be a driver of the ADV (Swingedouw et al. 2015; Drews et al. 2022).

The NAO is the leading mode of atmospheric winter variability over the NA region (Hurrell 1995; Hurrell et al. 2013). On the one hand, the NAO directly forces the NA SST north of the equator and the surface climate over the adjacent land areas. On the other hand, the NAO, by its impacts on surface-wind stress and air-sea heat exchange, influences the NA-ocean circulation (Eden and Jung 2001; Eden and Greatbatch 2003; Dong and Sutton 2005; Delworth and Zeng 2016; Årthun et al. 2021) and in turn with a time delay the SST and surface climate in the NA sector. The NAO is largely stochastic temporally with a power spectrum exhibiting near-white noise character. It is not clear whether the interaction between the NAO and the NA is “one-way” or “two-way”. To be specific, the ocean may respond passively as an integrator to the atmospheric forcing linked to the NAO (one-way; e.g., Frankignoul et al. 2000). While on the other hand, the oceanic changes may feed back onto the atmosphere (two-way; Marshall et al. 2001b; Bellucci et al. 2008; Sun et al. 2021, S21 hereafter).

Here we investigate the internally generated ADV that is observed in a preindustrial control integration of a version of the Kiel Climate Model (KCM). This simulation has been analyzed in S21, addressing the multidecadal variability over the NA that is the KCM’s version of the AMV. It is shown that the AMV and the ADV coexist in the model, with well-defined spectral characteristics but similar spatial patterns. We are concerned with the following questions regarding the ADV in the KCM: first, what is the space–time structure of the ADV? Second, what is the role of the ocean circulation in the ADV, in particular the roles of the gyre and overturning circulations? Third, what is the role of dynamical ocean–atmosphere interactions in the ADV? Finally, fourth, to which extent do the AMV and the ADV interact with each other? The paper is organized as follows. Section 2 gives an overview of the materials and methods used in this study. The model’s ADV is described in terms of selected variables in Sect. 3, and the major conclusions are presented in Sect. 4.

2 Materials and methods

2.1 Model data

Data from a preindustrial control run of 3000 years length with the KCM (Park et al. 2009) is used in this study. The atmospheric component of a version of the KCM is ECHAM5 (Roeckner et al. 2003) with a horizontal resolution of T42 (2.8°×2.8°) and 19 vertical levels. NEMO (Madec 2008) is used as the ocean-sea ice component, employing a horizontal resolution of 2° (ORCA2) with an increased meridional resolution of 0.5° near the equator, and with 31 vertical levels. The atmospheric component is coupled to the ocean-sea ice component by OASIS (Valcke et al. 2006). Annual-mean anomalies relative to the long-term mean are used.

We correct the surface-freshwater flux over the NA (10° N–80° N), which largely eliminates the upper-ocean salinity biases and partly the upper-ocean temperature biases over the region (Park et al. 2016, S21). Fig. S1 shows exemplarily the 2 m-temperature and SST biases (model-observations) in the uncorrected and corrected experiment.

2.2 Statistical methods and definition of climate indices

The Principal Oscillation Pattern (POP) method (Hasselmann 1988; von Storch et al. 1988, 1995, Appendix) is applied to the model data, which is a multivariate statistical method. The POPs may be regarded as the normal modes of a linearized system whose system matrix is estimated from data. Here, the POPs are calculated using a joint dataset comprised of three different variables over the Atlantic sector (0° N–70° N, 80° W–0° E): the barotropic streamfunction (PSI), representing the gyre circulation, the overturning streamfunction (MOC), representing the thermohaline circulation, and the sea level pressure (SLP), representing the low-level atmospheric circulation. Similar results are obtained when either only using SST and PSI or PSI and MOC in the POP analysis. The two leading POPs are complex. The most energetic POP mode (POP1), which is the KCM’s version of the AMV, is described in S21. In this study, we discuss the second most energetic POP mode (POP2) from the same POP analysis that describes the KCM’s ADV. We note that the patterns associated with POP2 share some similarities with those associated with POP1, suggesting that the two POP modes are closely linked to each other. We will return to this point in Sect. 3e.

The real-part coefficient time series of POP2 (PC2r) is used as an index of the ADV. The spatial patterns of oceanic and atmospheric variables linked to POP2 are derived by local linear regression, where the PC2r is normalized by its standard deviation. We also use correlation analysis, power-spectral and cross-spectral analysis. Details about these methods and the significance test are given in the Appendix. Different from S21, the background noise in the power-spectral analysis is reconstructed with autoregressive (AR) models following Mecking et al. (2014), because higher order AR process can better represent the MOC and SPG variability.

Besides PC2r, several other indices are used in this study. An NAO index is defined following Hurrell et al. (2013) as the PC of the leading EOF of the annual SLP anomalies over the extratropical NA region (20° N–80° N, 90° W–40° E), which accounts for 44.8% of the total variance. A MOC index is defined as the maximum of the overturning streamfunction at 40° N (Zhang 2008). An SPG index is defined as the inverted area average of the barotropic streamfunction over the subpolar NA (50° N–58° N, 42° W–26° W), similar to Lohmann et al. (2009).

3 Results

Figure 1 depicts the power spectra of the SPG index (Fig. 1a) and the MOC index (Fig. 1b). The background noise of the SPG index and MOC index is estimated by an AR(1) and AR(3) process, respectively. Both spectra exhibit statistically significant peaks above the 95% confidence level at multidecadal timescales, indicating that both SPG and MOC are essential to the AMV (S21). The power spectrum of the SPG index additionally exhibits statistically significant enhanced variability at decadal-to-bidecadal timescales.

Fig. 1
figure 1

a Power spectrum of the SPG index. b Power spectrum of the MOC index. The spectra are shown by blue solid curves, the 95% confidence level by orange dashed curves with reference to a fitted AR process (black dashed curves). The time series have been normalized by their respective standard deviations prior to computing the spectra

Full size image

POP2 describing the ADV and accounting for 15% of the total joint variance exhibits a rotation period of 18.9 years and a decay time of 6 years. The last 300 years of the real-part coefficient (PC2r, blue line) and imaginary-part coefficient (PC2i, orange line) time series are shown in Fig. 2a. The power spectra of PC2r (Fig. 2b) and PC2i (Fig. 2c), calculated from all 3000 years, both exhibit a broad and statistically significant peak in the decadal-to-bidecadal range relative to the fitted background AR(2) process. At these decadal-to-bidecadal timescales, the squared-coherence spectrum between PC2r and PC2i exhibits highly significant squared coherences (Fig. 2d, blue line) and a phase of about π/2 (90°, corresponding to 4.7 years), with PC2r lagging PC2i (Fig. 2d, orange stars). The statistically significant decadal-to-bidecadal peaks in the power spectra of the POP2 coefficient time series and in the squared-coherence spectrum confirm that the ADV in the KCM is clearly distinguishable from the AMV described by POP1.

Fig. 2
figure 2

a Real-part (PC2r) and imaginary-part (PC2i) time series of the second most energetic POP mode (POP2). Shown are the last 300 years. In b–d, all 3000 years of PC2r and PC2i are used. b, c Power spectrum of the normalized POP2 time series (blue solid curves) and the 95% confidence level (orange dashed curves) with reference to a fitted AR(2) process (black dashed curves). b PC2r, c PC2i. d Squared coherence (blue solid curve) and phase angle (orange stars) between PC2r and PC2i. A phase angle of zero indicates that the two time series vary in phase (orange dotted line), while a positive (negative) phase lag indicates that PC2r leads (lags). The phase lag is only shown when the squared coherence exceeds the 99% confidence level (black dotted line)

Full size image

The real-part and imaginary-part spatial patterns of POP2 (POP2r and POP2i, respectively) are shown in Fig. 3. The POP2r-MOC pattern (Fig. 3a) is characterized by a well-developed overturning-streamfunction anomaly in the extratropical region (40° N–60° N) with a maximum in the depth range 1000–3000 m. There, explained variances are relatively large with values exceeding 40%. The PSI pattern (Fig. 3b) related to POP2r exhibits a dipole structure, with negative PSI anomalies in the subpolar NA and positive PSI anomalies in the midlatitude NA, indicating a northward shifted North Atlantic Current (NAC) and stronger SPG (see also Fig. 4h). Explained variances in the POP2r–PSI pattern are relatively large in the negative and positive center, with values of up to 40% in the negative and 60% in the positive center. POP2 involves statistically significant ocean-circulation changes only in the extratropical NA, which is one of the main differences to AMV (POP1) that is associated with basin-wide overturning changes. We note that the AMV derived from instrumental observations exhibits significant tropical SST anomalies.

The SLP component of POP2r (Fig. 3c) exhibits an NAO-like pattern, with anomalously low-pressure north of about 60° N and anomalously high-pressure to the south. The explained variances in the SLP are generally small (we note the contour interval of 0.1), as expected given the chaotic nature of the extratropical atmosphere. Nevertheless, the explained variances amount to up to 20% in the low-pressure anomaly center. In the positive SLP-anomaly center to the south, the explained variances are typically 10%.

Fig. 3
figure 3

Spatial patterns of POP2. a–c Real-part patterns and d–f imaginary-part patterns. a, d MOC (meridional overturning circulation streamfunction, Sv), b, e PSI (barotropic streamfunction, Sv), and c, f SLP (sea level pressure, Pa). Color indicates loadings and contours explained variances. The explained variance for MOC and PSI is shown with a contour interval of 0.2, and for SLP with a contour interval of 0.1

Full size image

The MOC, PSI and SLP imaginary-part patterns of POP2 (POP2i) are shown in Fig. 3d–f, respectively. Theoretically, the imaginary part is in quadrature with the real part and leads by π/2, which is confirmed by the cross-spectral analysis of the two POP-coefficient time series (Fig. 2d). The imaginary-part patterns are considerably weaker compared with their real-part counterparts. The POP2i–MOC pattern exhibits weak negative overturning-streamfunction anomalies centered near 45° N and in the depth range 1000–3000 m, where the explained variances are less than 20% (Fig. 3d). The POP2i–PSI pattern (Fig. 3e) only exhibits significant anomalies in the subpolar region. The positive anomalies are observed in the Labrador Sea and along the southern flank of the SPG (see the grey lines in Fig. 4f–j for the climatological zero-PSI contour), where the explained variances are up to 30% (Fig. 3e). We note the small region of negative loadings in POP2i–PSI southeast of the southern tip of Greenland, with explained variances up to 20% in the center. The SLP anomalies linked to POP2i (Fig. 3f) are small, with explained variances of less than 10%. There are weak positive anomalies over the subpolar NA that partly reach into the Arctic, and there is a small area of small negative anomalies over Southeast Greenland.

POP2r can be considered as the mature phase of the ADV, exhibiting large loadings in all three variables, MOC, PSI and SLP, which suggests strong ocean-atmosphere interactions taking place at this phase of the POP cycle. In POP2i, only the PSI pattern exhibits noticeable anomalies that account for a relatively large fraction of the variability, further supporting an important role of the gyre circulation in the ADV. In the subsequent correlation, regression and cross-spectral analyses, the real-part coefficient time series, PC2r, is used as an index of the ADV.

3.1 Ocean-circulation changes

In order to obtain more insight into the dynamics of the ADV, we first investigate the changes in ocean-circulation related to the ADV. Maps of local regression coefficients of MOC and PSI upon PC2r are shown at different time lags (Fig. 4a–j), ranging from lag − 4 years to lag + 4 years, representing approximately one half of the POP2 cycle. The MOC evolution (Fig. 4a–e) describes the transition from a normal (Fig. 4a) to an anomalously strong overturning (Fig. 4e). As expected, the pattern at lag zero (Fig. 4c) closely resembles the real-part MOC pattern of POP2 (Fig. 3a), with positive anomalies in the extratropical NA in the region 40° N–60° N. The MOC pattern at lag − 4 years (Fig. 4a) resembles the POP2i–MOC pattern (Fig. 3d), which is consistent with the phase lag between PC2r and PC2i (Fig. 2d).

Fig. 4
figure 4

Lead–lag regression maps upon PC2r from lag − 4 years to lag  4 years (at 2-year intervals) of (a–e) meridional overturning streamfunction (MOC), f–j barotropic streamfunction (PSI) with top 50–700 m depth-integrated ocean currents superimposed (arrows), k–o SLP. The grey lines in h indicate PSI climatology with a contour interval of 10 Sv. The grey lines in f–j indicate the zero contour of the PSI climatology. Dots indicate regressions that are not significant at the 95% confidence level. Explained variances in a–e and k–o are indicated by black contours. The contour interval in a–e is 0.2, in k–o 0.01

Full size image

The PSI regression maps (Fig. 4f–j) are superimposed with the depth-integrated upper-ocean (50–700 m) currents. At lag − 4 years, negative PSI anomalies develop in the Irminger Sea (Fig. 4f). Thereafter, the negative anomalies strengthen and expand until an anomalously strong SPG has fully developed at lag zero (Fig. 4h). Further to the south, positive PSI anomalies develop concurrently. Figure 4h is overlapped with the complete PSI climatology, from which the northward shift of the NAC can be inferred. Around lag zero, both SPG and NAC are enhanced. The dipole pattern, with negative anomalies in the subpolar NA and positive anomalies to the south, weakens after lag zero (Fig. 4i–j).

3.2 SLP changes

The real-part SLP-anomaly pattern of POP2 is the positive phase of the NAO (NAO+, Fig. 3c). Figure 4k–o depicts the lag-regression maps of SLP upon PC2r. From lag − 4 years to lag zero, the regression maps demonstrate gradually strengthening NAO + patterns (Fig. 4k–m), where the explained variances also become larger. After lag zero (Fig. 4n, o), the SLP-anomaly patterns fundamentally change exhibiting extremely small anomalies. Small low-pressure anomalies are located over the extratropical NA with a center over Southern Greenland, which only account for a minor fraction of the variance but are still statistically significant.

We calculate the cross-correlation function between PC2r and the NAO index (Fig. 5a). The cross-correlation function is highly asymmetric, exhibiting larger correlation coefficients when the NAO index leads and a maximum of 0.45 at lag zero. The correlation coefficients become very small and mostly insignificant when the NAO index lags PC2r. The switch in the structure of the SLP anomalies after lag zero (Fig. 4n, o) is consistent with the sharp drop in the cross-correlation function after lag zero. Overall, the regression and correlation analyses suggest a leading role of the NAO in the ADV. When using winter-mean SLP data, the relationship between PC2r and the NAO index is similar but more pronounced (not shown).

Fig. 5
figure 5

Cross-correlation functions of PC2r with a the NAO index, b SPG index, c MOC index, and d PC1r, the real-part coefficient time series of POP1 describing the AMV in S21. The time lag is in years. Filled dots indicate that the correlations are statistically significant at the 95% confidence level

Full size image

The development of the wind-stress anomalies related to PC2r features strengthening westerly wind stress between the two poles of the NAO + pattern up to lag zero (not shown). This intensification goes along with strengthening positive wind-stress curl over the subpolar and strengthening negative wind-stress curl over the midlatitude NA. After lag zero, the wind stress and wind-stress curl anomalies become rather weak, which is in line with the very weak SLP anomalies at these time lags.

The NAO + phase not only enhances the wind-stress curl but also oceanic heat loss of the subpolar NA, which will drive stronger deep convection and eventually stronger overturning. The cross-correlation functions between PC2r and the SPG index (Fig. 5b) and PC2r and the MOC index (Fig. 5c) elucidate the relationship between NAO+ (as given by POP2r, Fig. 3c) and ocean circulation. The correlation coefficient of PC2r with the SPG index is 0.52 at lag zero and 0.54 when PC2r leads by 1 year. The correlation coefficients decrease to zero at lag ± 5 years and exhibit a minimum of − 0.35 at lag + 9 years. Thus, SPG strength is almost in phase with the NAO. The cross-correlation function of PC2r with the MOC index is quite different to that of PC2r with the SPG index, exhibiting a maximum of 0.56 when PC2r leads by 5 years, while the correlation coefficients are rather small around lag zero. Thus, the MOC exhibits a delayed response to the low-frequency part of NAO+, with a time lag of several years. In summary, the relationships between SLP, SPG and MOC support the notion that in the KCM, persistent changes in the NAO drive the ocean-circulation changes related to the ADV, both gyre and overturning circulation, consistent with Eden and Jung (2001).

3.3 SST changes and the origins

Generally, SST anomalies linked to the KCM’s ADV (Fig. 6a–e) go along with heat-flux anomalies of the opposite sign (Fig. S2), indicating that the SST anomalies in most regions are driven by the ocean. Further, major SST and heat-flux anomalies are restricted to the middle and high-latitude NA. At lag zero, the SST-anomaly pattern exhibits a tripolar structure in the meridional direction with two regions of positive SST anomalies, one in the mid-latitude and the other in the subpolar NA, and a region of negative SST anomalies in between (Fig. 6c). At this time, the largest SST anomalies are observed in the mid-latitude NA near 45° N, which extend from the coast of North America to the central NA. The explained variances are the largest close to the North American coast with values up to 50%. At lag + 4 years, the SST anomaly pattern has changed to a monopole located in the subpolar NA (Fig. 6e).

The upper-ocean (0–700 m) heat content in the subpolar (western half of the midlatitude) NA is decreasing (increasing) until lag zero (Fig. 6f–h). We note that no significant SST and heat-content changes are observed in the low latitudes, which is a major difference to the KCM’s AMV. At the same time, the mixed layer depth around southern Greenland is increasing until lag zero (Fig. 6k–m), indicating enhanced deep convection eventually spins up the overturning (Fig. 4c–e).

Fig. 6
figure 6

Lead–lag regression maps upon PC2r from lag − 4 years to lag + 4 years (at 2-year intervals) of a–e sea surface temperature, f–j top 700 m ocean heat content, k–o mixed layer depth. Dots indicate regressions that are not significant at the 95% confidence level. Explained variances are indicated by contours. The contour interval is 0.1

Full size image

The total northward heat transport (Fig. 7, black curves) is decomposed into an overturning component (MOC, blue curves) and a gyre component (gyre, orange curves) by considering the zonal mean and deviations from the zonal mean of the meridional velocity and temperature field, respectively (Bryden and Imawaki 2001). The heat-transport convergence (divergence) determining the temperature tendency is the negative (positive) meridional derivative of the heat transport. From lag − 4 years to lag zero, the total northward heat transport is dominated by the gyre component. The development of positive SST anomalies up to lag zero in the subpolar and western midlatitude NA (Fig. 6a–c) is consistent with the gyre heat-transport convergence in these regions (Fig. 7a–e). The concurrent development of cold SST anomalies between the positive SST anomalies is consistent with the gyre heat-transport divergence around 50° N. After lag zero, convergence of both gyre and overturning heat transports (Fig. 7f–i) contribute to the development of positive SST anomalies in the subpolar NA (Fig. 6d, e). This is another major difference to the AMV, in which the overturning heat transport plays a dominant role after lag zero (Fig. S3).

Fig. 7
figure 7

a–i Regressions upon POP2r from lag − 4 years to lag + 4 years (at 1-year intervals) of the total northward heat transport (black), MOC component (blue) and gyre component (orange) in the North Atlantic. Units in PW

Full size image

3.4 Phase reversal

What determines the timescale and cyclic nature of the KCM’s ADV? As the heat content of the subpolar NA has a strong influence on the SPG strength, which is also the case in the AMV (S21), we hypothesize that the heat storage in the upper subpolar NA is a key to the mechanism underlying the ADV. What differs between the ADV and the AMV is the respective role of SPG and MOC in changing the heat content: the SPG plays a much more important role in the ADV than in the AMV in which the MOC is dominant (Fig. S3). Because gyre and overturning heat transports are both important in the ADV (Fig. 7g–i), the time needed to change the heat content in the subpolar NA is much shorter relative to the AMV. This gives rise to the shorter period of the ADV. The adjustment of the ocean circulation to the atmospheric variability linked to the NAO also provides through heat transport a delayed negative feedback that enables oscillatory behavior. However, the lack of positive ocean-atmosphere feedback leaves the ADV strongly damped.

3.5 ADV and AMV

The picture that emerges is that the ADV, as described by POP2, constitutes a mixed oceanic gyre-overturning eigenmode that is excited by the atmospheric variability linked to the NAO. The oceanic patterns of POP2 exhibit some similarities to those associated with the AMV (POP1, S21), which by itself suggests that the two modes are connected in some way. The POP2–MOC patterns, however, exhibit anomalies only in the extratropical NA, which differs from POP1 that comprise basin-wide MOC changes. Besides, the atmospheric patterns differ between POP1 and POP2: POP1 is characterized by an SLP-anomaly pattern that is a monopole centered over southern Greenland, whereas it is the NAO pattern in POP2.

We next calculate the cross-correlation function between the real-part time series of POP1 and POP2, PC1r and PC2r, respectively. There is a maximum of 0.52 at lag zero (Fig. 5d), which is statistically significant at the 95% confidence level, indicating a close link between the decadal-to-bidecadal and the multidecadal variability. Cross-spectral analysis is utilized to understand the behavior of the SPG and MOC on the decadal-to-bidecadal and multidecadal timescales. There is some overlap of the peaks in the squared coherence spectra of PC1r and PC2r with the MOC index (Fig. 8a, c) and with the SPG index (Fig. 8b, d), indicating that the MOC and SPG are important in both the AMV (POP1) and ADV (POP2). The squared coherence between PC1r and the MOC index (Fig. 8a) is the highest on multidecadal timescales (30–50 years, orange rectangular window) and beyond, and the two indices vary almost in phase at POP1’s rotation period of 38 years. The squared coherence between PC1r and the SPG index (Fig. 8b) exhibits a relatively narrow but pronounced peak at a timescale considerably shorter than POP1’s rotation period. The comparably long timescale of POP1 is thus mainly due to the MOC (S21).

Fig. 8
figure 8

Squared coherence (blue solid curves) and phase spectrum (orange stars) between a PC1r and the MOC index, b PC1r and the SPG index. c PC2r and the MOC index, and d PC2r and the SPG index. The period range 20–30 years is highlighted by blue vertical lines and the range 30–50 years by red vertical lines

Full size image

PC2r and the MOC index (Fig. 8c) exhibit the highest squared coherence in the decadal-to-bidecadal timescale range (blue rectangular window) near PC2’s rotation period of about 19 years, with PC2r leading by about π/2 (about 5 years). This time lag is consistent with the cross-correlation function between PC2r and the MOC index (Fig. 5c). The squared coherence of PC2r and the SPG index (Fig. 8d) exhibits a broad peak at decadal-to-bidecadal timescales with little phase lag. Consistent with the heat transports (Fig. 7), the cross-spectral analysis suggests that the MOC and SPG are equally important in the ADV.

4 Conclusions

This study focuses on the analysis of the NA decadal-to-bidecadal variability, termed here ADV, which is simulated in a multimillennial preindustrial control integration of a version of the KCM. To derive the leading modes of long-term variability over the NA sector, a POP analysis was applied to the barotropic streamfunction, the overturning stream function, and the atmospheric sea level pressure. In a previous paper (S21), we discussed the multidecadal variability simulated in the same integration that is described by the leading POP mode (POP1) accounting for 23.7% of the total joint variance, while the second most energetic POP mode (POP2, 15%) corresponds to the ADV which is investigated in this study.

Both POP1 and POP2 are oscillatory with a relatively well-defined rotation period, as exemplified by the statistically significant peaks in the power spectra of the corresponding POP-coefficient time series and in the power spectra of selected indices. The two modes involve significant variability of the SPG and of the MOC, i.e., they can be considered as mixed gyre-overturning modes. In both modes, the SPG and MOC play an important role in driving the upper ocean-heat content in the subpolar NA, which is a key variable in the phase reversal of the modes.

Both types of variability can be regarded as damped eigenmodes of the ocean circulation in the NA that are excited by the atmospheric noise, where the surface-heat flux acts as damping on the long-term SST anomalies. The ADV involves long-term variability of the NAO (Fig. 3c), which is a major difference to the AMV where the prevailing SLP-anomaly pattern is a monopole centered over Southern Greenland. In both modes, the delayed negative feedback is provided by the ocean circulation and associated meridional heat transports into the subpolar NA (Fig. 7 and Fig. S3). In the ADV, both gyre and overturning circulation equally contribute to the heat transport. On the contrary, the MOC dominates the heat transport in the multidecadal mode. The shorter period of variability, the ADV, in comparison to the AMV can be explained by the faster heat built-up in the upper ocean in the subpolar NA region that slows the gyre circulation and eventually the MOC. Finally, the decadal-to-bidecadal variability and the multidecadal variability are significantly correlated, which would require long time series to distinguish the two modes and to study their interaction provided the KCM results carry over to the real world.

The key elements of the mechanism underlying the ADV simulated by the KCM are consistent with previous studies on decadal-to-bidecadal variability. The important role of the atmospheric noise in exciting the ADV is consistent with Saravanan and McWilliams (1997), who were using an idealized ocean-atmosphere model. The major role of the ocean dynamics in long-term variability of the NA SST is consistent with Wills et al. (2019) and Zhang et al. (2019). The importance of the gyre circulation in setting the decadal timescale has also been suggested. By using observational data, Nigam et al. (2018) find that the ADV timescale is set by the low-frequency portion of the NAO and cross-basin transit of the Gulf Stream’s detached eastern section. The meridional heat transport into the subpolar region connected with both the MOC and SPG is important in the KCM, where the roles of the MOC and SPG have also been discussed before (Menary et al. 2015; Oldenburg et al. 2021). In recent climate-model studies, MOC variability with a timescale of 20–30 years is dominated by a weakly damped oceanic mode associated with westward propagating baroclinic Rossby waves (Sévellec and Fedorov 2013; Muir and Fedorov 2017). The mechanism described in this study is different, as there are no westward propagating signals of the SST gradients observed in the extratropical NA (not shown). This may be due to the different atmospheric forcing in the KCM or to the diverse representations of upper-ocean temperature biases in other models.

Our results highlight the importance of the ocean-heat content in the subpolar NA and the roles of the gyre and overturning circulation in the phase reversal. The specific mechanism of the ADV in our study has not been described previously. However, the results may be model dependent, but analyses of the decadal-to-bidecadal timescale variability in other models is beyond the scope of this study. Longer instrumental data or proxy records are needed to determine whether the processes and mechanisms operating in the KCM also occur in the real world, especially regarding the coexistence of the ADV and the AMV.

Finally, the ADV accounts for a considerable amount of the SST variability over the extratropical NA in the KCM, especially in the NAC-system region and subpolar NA. These areas may thus be regions of relatively high potential SST predictability. It is still an open question, however, to which extent the extratropical NA SST influences the atmosphere and thus the surface climate over the adjacent land areas.