Introduction

The Northern Hemisphere (NH) is characterized by several multidecadal climate trends that have been attributed to anthropogenic climate change1. These trends include the large-scale Atlantic warming2, strengthening of the positive North-Atlantic-Oscillation (NAO) or Northern Annular Mode (NAM)3,4,5, Artic stratospheric cooling4,6,7, multidecadal modulation in the Arctic sea-ice8, and changes in Eurasian and Mediterranean temperature and precipitation4,9,10. The recent weakening and reversing of several NH-climate trends, despite the ongoing anthropogenic forcing (Figs. 1 and 7)5,11,12, has increased the interests in understanding and predicting the decadal-to-multidecadal climate variability13,14,15. Several efforts have been made5,6,11,16,17,18,19,20,21, but a general physical framework linking these trends is still lacking.

The North-Atlantic is characterized by multidecadal fluctuations of sea surface temperatures (SSTs) (Methods) known as Atlantic-Multidecadal-Oscillation (AMO) or Variability (AMV, Fig. 1a–d)5,14,22,23. A positive (negative) AMV index (AMVI) describes basin-scale surface warming (cooling) that is largest in the subpolar region. The basin-scale AMV-temperature change is in general characterized by a large degree of spatial coherence across the entire North Atlantic basin24. The AMV has been linked to major climate shifts worldwide. Although stochastic atmospheric dynamics can cause low-frequency SST variability through thermodynamical processes25,26, long-term unperturbed coupled climate model simulations indicate that the AMV may be an internal mode of climate variability driven by fluctuations in the Atlantic meridional overturning circulation (AMOC)27,28. A comparison between slab ocean and fully coupled atmosphere-ocean models shows that the spatial coherence in the SST pattern associated with the AMV can be simulated only, if the AMOC-related ocean dynamics is included24. In addition, external forcing agents, such as changes in total-solar-irradiance as well as aerosol, can also play a role29,30.

Fig. 1: Multidecadal NAO, AMV and Arctic sea-ice variability in the observations and reanalysis.
figure 1

a Multidecadal fluctuations of the wintertime (JFM) standard AMVI (SAMVI, dashed blue, Methods), local AMVI (LAMVI, blue, Methods), NAOI (red), upward oceanic turbulent heat-fluxes over Labrador Sea areas (green, 45–55 W and 45–60 N) and scaled time-integrated NAOI (dashed red). b Comparison between the wintertime (JFM) NAOI (red), AMVI (blue, Atlantic SST between 0–70 N) and the averaged sea-ice concentration over the Arctic region 60W-50E and 50–90 N. c, d Wintertime local multidecadal SST anomalies associated with strong multidecadal North-Atlantic Warming (years 1886–1889, 1949–1952 and 2010–1013) and cooling (years 1916–1919 and 1980–1983) respectively (in °C). e, f Wintertime 1000 hPa multidecadal geopotential height anomalies associated with strong positive (1921–1924 and 1989–1992) and negative (1887–1890 and 1959–1962) multidecadal NAO-phase (in m). All indices in a and b have been standardized by dividing the resulting deviations from the mean by the standard deviation, in such way that the values express the portion of the standard deviation explained by the multidecadal changes.

The North-Atlantic region is also known for significant multidecadal fluctuations in the NAO5,31, which is strongly related to the NAM32,33 (Fig. 1a, b, e, f). Positive NAO- and NAM-phase are associated with strengthening of the high-latitude westerlies and Stratospheric-Polar-Vortex (SPV) as well as poleward shift of the eddy-driven jet, storm-track and weather systems over the North-Atlantic and European regions, with opposite conditions for the negative phase. A well-known manifestation of multidecadal NAO variability is the transition from negative to positive phase between the 1960s and 1990s12. It was associated with significant stratospheric cooling7, increase in precipitation, warming, and snow melting over the Northern Eurasian continent, as well as cooling and precipitation decrease over southern Europe and the Mediterranean4,33. The observed large-scale AMV-warming was accompanied by a multidecadal NAO-shift into negative phase. Vice versa for the large-scale AMV-cooling (Figs. 1 and 2). This multidecadal transition into positive (negative) NAO-phase was associated with a downward propagating positive (negative) NAM-signal from the stratosphere into the troposphere and by thermal wind balance with cooling (warming) in the lower polar stratosphere (Fig. 2). Studies using stand-alone atmospheric models driven by observed and coupled model AMV warming patterns show that such changes in the stratosphere/troposphere-coupled circulation can be forced by the North-Atlantic-SST5,31,34. The AMV-induced changes in the tropospheric baroclinicity, changes in the tropical SST and associated upward planetary-wave propagation can all play a role5,34,35. In particular, the changes in the SST-front and associated low-level baroclinicity along the Golf stream can also lead to wave-driven NAO-like circulation response34,35,36.

Fig. 2: Seasonal evolution of the stratospheric changes associated with the AMV in the NCEP-reanalysis.
figure 2

a, b Seasonal evolution of the multidecadal NAM indices between 1000 and 10 hPa for strong multidecadal positive (1989–1992) and negative (1959–1962) NAO-phases, respectively. c, d The evolution of the lower stratospheric temperature (100–70 hPa) associated with a and b, respectively (in °C).

Positive NAO-phases are also associated with enhanced upward turbulent heat flux over Labrador Sea and delayed large-scale Atlantic warming. Vice versa for the negative NAO-phases (Fig. 1a). The large-scale North-Atlantic SST changes can be explained by the delayed response of the AMOC to the NAO and the related poleward heat transport, as demonstrated by stand-alone ocean models20,37,38 and coupled atmosphere-ocean models18,19. The oceanic response to the NAO typically involves modifications in water density over the Labrador Sea due to changes in the turbulent heat-fluxes39. Consistent with the Hasselmann model, the AMV/AMOC response to the NAO can be seen as an accumulation of the time-integrated effect of the NAO. The NAO acts as forcing that drives the time derivative of the AMV/AMOC25,40,41. The time lag between the NAO and the AMOC is explained by the large oceanic inertia. The time-integrated NAO can therefore be seen as proxy for the AMOC40 and its simulated change support the role of the NAO in driving the AMOC and AMV (Fig. 3a).

Fig. 3: Stratosphere/troposphere/ocean coupled damped multidecadal oscillation in model simulation.
figure 3

a The reconstructed NAOI (red), AMVI (blue), AMOCI (yellow) and scaled time-integrated NAOI (dashed red) from the coupled multidecadal stratosphere/troposphere/ocean variability using MSSA. b Running explained variance of the AMVI (blue), NAOI (red), NAMI at 20 hPa (dashed red) and AMOCI (defined as Atlantic Ocean Meridional Overturning Mass Stream function averaged, around the climatological maximum, 0–2000 m and 30to40N, yellow). c The averaged periodic evolution over 16 phases of multidecadal changes (using phase-composite, Methods) of the reconstructed multidecadal NAOI (red), AMVI (bleu), AMOCI (yellow) surface salinity (dashed green), and the averaged upward turbulent oceanic heat-fluxes (green) over Labrador Sea. d The same as c but for the reconstructed Gyre index (defined as difference in the barotropic mass stream function between subtropical and subpolar gyre (Methods, yellow)) and the reconstructed average over Labrador Sea of the Evaporation (Ev, blue), surface salinity (green) and difference between evaporation and precipitation (Ev–Pr, red). All indices in ad have been standardized by dividing the resulting deviations from the mean by the standard deviation, in such way that the values express the portion of the standard deviation explained by the multidecadal changes.

There is thus modeling evidence for two-way coupling between the stratosphere/troposphere-coupled system and the AMV (e.g., refs. 5,18,20,34). Whether this two-way coupling can lead to an oscillation in the climate system is one of our major questions in this study. Delayed oscillations can exist without external forcing40 arising from positive feedback that increases the fluctuation amplitude and subsequent delayed negative feedback that reverses the oscillation to an opposite phase. Such feedback chain has been used to construct conceptual oscillation model of El Niño-Southern Oscillation ENSO42, multidecadal NAO variability40 and the observed periodic fluctuations in the atmospheric baroclinic annular mode43. Observational and model data were previously used40 to derive a delayed oscillation model of the quasi-periodic NAO and North Atlantic SST-Tripole pattern (NAT). In this oscillation the positive feedback in the growing phase of the AMV is maintained by the impact of the time-integrated NAO-forcing on the AMV/AMOC following the Hasselmann model25 and the delayed negative feedback results from the delayed response of the NAT to AMOC-induced heat transport. The delayed NAT-response acts to reverse the NAO through the well-known NAT feedback on the NAO44. It was also demonstrated that both the magnitude and timing of the delayed multidecadal NAO-forcing are relevant for the period of the oscillation40. These results are consistent with previous model studies showing the importance of the NAO for the observed variability of the AMO/AMOC18,20.

There is observational (Fig. 1b45) and modeling46 evidence for multidecadal fluctuations in the Arctic sea-ice. A recent model study suggested that the observed NAO multidecadal fluctuations and the associated meridional overturning circulation response can explain a considerable part of the observed multidecadal Arctic sea-ice variability19. Finally, there is observational evidence that the anomalies associated with the AMV-signal including NAO can propagate throughout the Northern Hemisphere via a sequence of atmospheric and lagged oceanic teleconnections. This propagation has been termed “stadium wave” and can spread over the whole northern hemisphere leading to hemispheric impact of the AMV47,48,49,50. This hemispheric impact of the NAO/AMV coupled variability has been used to construct an NAO-based linear prediction model for the NH surface temperature41.

The purpose here is (1) to assess the extent to which a coupled delayed oscillator framework in transient forced climate can help to understand the observed coherent multidecadal variability in the NH stratosphere/troposphere/ocean-coupled system including sea-ice; (2) to identify the mechanisms connecting the different climate components involved and the climatic impacts of this oscillation; and (3) to demonstrate that the coupled atmosphere/ocean variability can enhance the skill of statistical decadal-to-multidecadal climate projection. Since the damped delayed quasi-oscillatory behavior in the coupled NAO/AMOC system can occur without time varying external forcing40,51, the focus will be on the physical processes and feedbacks maintaining the coupled stratosphere/troposphere/ocean variability and not on the role of external radiative forcing.

Results

Multidecadal stratosphere/troposphere/ocean-coupled oscillation

Because of the short observational records, long-term simulations with ocean/atmosphere-coupled models are highly needed. The multidecadal atmosphere/ocean-coupled variability and oscillations in the North-Atlantic has been considered mainly by excluding the impact of time varying external forcing and neglecting stratospheric processes (e.g., refs. 40,52,53). Such coupled oscillatory mode of variability has been identified between tropospheric NAO and AMOC through the Atlantic SST-tripole40. However, external forcing, like volcanic aerosol, can also be important. The AMV impact on the NAO through stratosphere/troposphere coupling, as detected in reanalysis (Fig. 2) and simulated with stratosphere-resolving Atmospheric-General-Circulation model (AGCM)-experiments5, and the possible role of sea-ice has not been considered within the framework of atmosphere/ocean-coupled oscillations. Therefore, we extracted the multidecadal stratosphere/troposphere/ocean-coupled oscillation from a millennial-scale transient climate simulation using the stratosphere-resolving atmosphere/ocean-coupled Max Planck Institute Earth System Model for Paleo applications (MPI-ESM, Methods). The simulation includes the effects of time varying historical natural and anthropogenic external forcing, following the protocol of the Paleoclimate Model Intercomparison Project phase 3 (PMIP3). To isolate the coupled stratosphere/troposphere/ocean-oscillation, we performed Multichannel Singular Spectrum Analysis (MSSA, Methods) on the combined coupled fields from stratospheric (20 hPa, NDJ) and tropospheric (1000 hPa, JFM) winter geopotential height, winter Atlantic SSTs (JFM) and annual mean Atlantic overturning mass stream function. These variables constitute the most relevant components of our coupled analysis. The reanalysis shows that the stratosphere/troposphere coupling is part of the multidecadal NAO variability (Fig. 2a, b) and the atmosphere-only model experiments5 show that the NAO can respond to the AMV only if the stratosphere is included. The inclusion of the stratosphere is thus important in capturing the observed dynamics of the NAO response to the AMV, which includes the wave-driven stratosphere/troposphere coupling. The inclusion of the AMOC helps in capturing the dynamical response of the SST that results from the accumulated or time-integrated impact of the NAO-forcing on the AMOC following the Hasselmann model25,40. All the other variables are used just for understanding the mechanisms and the implications of the detected multidecadal variability. Since the NAO (AMV) can respond to the AMV (NAO) without the involvement of the sea-ice forcing5,20, we include only the NAO, AMV and the components that are important for their dynamical feedback (stratospheric circulation and AMOC mass stream function) in the basic MSSA analysis. The MSSA is a multichannel extension of the singular spectrum analysis, which is used to extract non-linear oscillations. The MSSA is based on the time-space Empirical Orthogonal Function (EOF) analysis and extracts both temporal and spatial patterns of oscillatory modes embedded in the data.

The MSSA of the MPI-ESM-simulation isolates a coupled stratosphere/troposphere/ocean oscillation in the second and third MSSA modes (Fig. 3). This has a characteristic period of ca. 50-years. It has the properties of a damped oscillation, which is characterized by gradual changes in the amplitude (Fig. 3a, b). In particular, the running explained variance over about one and a half period of the oscillation (defined as ratio between the running variance of the extracted multidecadal time series and the running variance of the raw non-filtered time series, Methods) varies substantially through time (Fig. 3b). In this oscillation, the standardized AMOC-index features the highest maximal running explained variances, exceeding 43%. This reflects the inertia of the internal ocean dynamics. The standardized NAM-index (or NAO index) shows the lowest running explained variance, reaching only 10%. This reflects the energetic atmospheric internal variability and the high-frequency character of extratropical atmospheric dynamics. The running explained variance of the oscillatory component of the AMVI can account for up to about 30% of the total variance, which is in-between the values for the NAO (NAM) and the AMOC. This reflects the combined impact of both atmosphere and ocean on North-Atlantic SSTs. Overall, the strength of the simulated multidecadal oscillation is still comparable to the (short) observational record. The observed multidecadal fluctuation of the standard AMVI, which includes the impact of global SSTs, can explain around 48% of the variance over the last one and a half AMV fluctuations and the maximal standardized AMVI-values can reach over 0.5 (Fig. 1a). The maximal variance explained by the reconstructed AMVI and the associated standardized values in the model remain smaller. However, the maximal changes in the observed standardized Local AMVI (LAMVI, Methods), which exclude the impact of global SSTs, remain comparable to the MSSA-based AMV index extracted from the model. For the NAO index (NAOI) the running variance explained by the reconstructed multidecadal fluctuations (can reach up to 10% in the model) and associated maximal standardized values (can reach about 0.17) remain smaller than the corresponding reanalysis NAOI-values (14% and about 0.23, respectively). The AMOC and the time-integrated NAO vary in phase, which support the general view that the NAO is forcing the AMOC and that AMOC can be seen as an accumulation of the NAO-forcing following the Hasselmann model.

Life cycle of the oscillation

By decomposing the 50-year period of the oscillation into 16 phases and performing a composite analysis over each individual phase (each corresponding to about 3 years, Methods), one can disentangle the main mechanisms and feedbacks maintaining this oscillation (Fig. 3c, and d, Fig. 4 and Supplementary Fig. 8). Consistent with the observations and reanalysis, the phase composite indicates that positive NAO (phase 1–2) is associated with downward propagating positive NAM-anomalies from stratosphere into troposphere (Supplementary Fig. 1a). In agreement with reanalysis, the wintertime strengthening of the SPV during multidecadal NAM-intensification is associated with stratospheric cooling (Supplementary Fig. 1b). Because of the damped character of the oscillation, the averaged changes from the phase-composite (Fig. 4) are smaller than the changes seen in short observational record. However, the individual reconstructed stratosphere/troposphere coupled multidecadal changes can be comparable with the observations (Supplementary Fig. 2, Figs. 12). But the simulated standardized NAO index remains smaller compared to the observations. The positive NAO-phase results in a delayed AMOC-strengthening, which is caused by a thermal- and haline-driven increase of upper-ocean water density and mixed-layer depth, through changes in both NAO-induced upward turbulent heat fluxes and freshwater (Figs. 3c, d and 4b–d). The increase of upward surface turbulent heat fluxes and the associated deepening of the mixed-layer in the Labrador Sea can be explained by the strong NAO-induced high-latitude westerlies and associated transport of continental North-Canadian cold air. The strong high-latitude westerlies can also contribute to the mixed-layer deepening through wind-induced dynamical mixing. The enhanced high-latitude salinity following the positive NAO can result from the delayed strengthening of the wind-driven subtropical and subpolar Atlantic-gyres (Fig. 3d and 4e), which transport the saltier subtropical water-masses northwards. The enhanced salinity can also result from local NAO-induced evaporation increase (Fig. 3c, d), which acts to reduce the freshwater supply in the Labrador Sea. The delayed AMOC-strengthening reaches the maximum between phases 4 and 6 and is associated with an intensification of the North Atlantic subpolar gyre (Figs. 3c, d and 4e, f). Therefore, both the gyre and the overturning circulations transport more heat poleward (Fig. 4g), leading to North-Atlantic warming peaking between phases 6 and 7 (Fig. 3c and 4h). The NAO-induced North-Atlantic warming can also result from the high-latitude deepening of the mixed-layer (Fig. 4d), which can dampen the effect of upward turbulent heat flux and lead to relative warming by increasing the oceanic heat-capacity51.

Fig. 4: Dynamics of the damped coupled stratosphere/troposphere/ocean oscillation.
figure 4

a, b The MSSA-reconstructed 1000 hPa geopotential height (in m) and upward oceanic turbulent heat-fluxes (in w m−2) respectively at phase 2. c, d Surface salinity (in psu) and mixed-layer depth (in m) during phase 3, respectively. e Ocean barotropic mass stream function during phase 4 (in 10 8 Kg s−1). f Latitude-depth section of the ocean meridional overturning mass stream function in the Atlantic Ocean (north of 30° S) during phase 5 (in 108 Kg s−1). g The zonally averaged poleward Atlantic heat transport in W (during the phase 5) by the total oceanic circulation (red), the AMOC (yellow) and the oceanic gyre (blue). h The reconstructed SST anomalies at phase 7 (in °C), i the reconstructed 1000 hPa geopotential heights at phase 8 (in m).

Consistent with the reanalysis (Figs. 1 and 2) and stand-alone AGCM-experiments5,34, the atmosphere adjusts to the Atlantic warming with changes towards the negative NAO and with a downward propagating negative NAM from stratosphere into troposphere (Figs. 3c, 4i and 5a). The weakening of the SPV is associated, by thermal wind balance, with warming in the stratosphere (Fig. 5b). The momentum budget in the Transformed Eulerian mean (TEM) formalism (Methods, Fig. 5c–f) indicates that in the free atmosphere the seasonal weakening of the westerlies towards the negative NAM and the associated downward NAM-migration following the warm AMV phases are maintained mainly by the resolved atmospheric wave forcing against the non-resolved and residual Coriolis forcing in the upper troposphere and stratosphere. Near the surface, the negative NAM is maintained mainly by the residual Coriolis forcing. These conditions are reversed during the positive NAM phases (Supplementary Fig. 1). The negative NAO is associated with low-latitude westerly wind anomalies that weaken the trade winds and upward turbulent heat-fluxes (Fig. 4i, Supplementary Fig. 3a, b). This can initiate a wind-evaporation-SST-feedback54 warming the tropical SSTs further. Such tropical SST-warming has been shown to be crucial for further intensification of the negative NAO in response to the AMV5. The wind-evaporation-SST-feedback can also enhance the upward surface heat flux and therefore strengthen the lower-latitude cooling during positive NAO (Fig. 4a, b). The enhanced negative NAO in phases 7–10 weakens the upward turbulent heat-fluxes and salinity over the Labrador Sea that yields shallower mixed-layer and delayed AMOC-weakening (Fig. 3c, d, Supplementary Fig. 3 and Supplementary Fig. 8). This delayed negative feedback reverses the oscillation towards its cooling phase (phases 13–15, Supplementary Fig. 3h) and ultimately leads to trends towards positive NAO (Supplementary Fig. 3i).

Fig. 5: The dynamics of the seasonal evolution of the stratosphere/troposphere-coupled system in the coupled multidecadal stratosphere/troposphere/ocean oscillation.
figure 5

a The seasonal evolution of reconstructed (Methods) multidecadal NAM-index during the phase 9 at different levels between 1000 and 0.1 hPa. b The zonally averaged lower stratospheric temperature (100–70 hPa) during the phase 9 (in °C). c As in a but for the combined contribution of resolved wave and non-resolved parameterized forcing (i.e., total forcing, in m s−1/day). df As in c but for d resolved wave forcing, e non-resolved forcing and f the Coriolis torque of the residual circulation respectively (in m s−1/day).

Response of sea-ice and surface climate

The Arctic sea-ice strongly responds to the fluctuations of the coupled atmosphere/ocean system in the North-Atlantic19,45,46. The NAO can impact the Arctic sea-ice directly through the poleward atmospheric heat transport and wind-stress55, and indirectly through changes in the ocean circulation19. The change of sea-ice in phase with the positive NAO shows a dipole structure (Fig. 6a, b) with negative (positive) sea-ice concentration anomalies over Greenland and Barents Seas (Labrador-Sea), which is consistent with the direct NAO impact on sea-ice. However, the most prominent multidecadal Arctic sea-ice melting results mainly from the ocean adjustment to positive NAO between phases 3 and 4, where the poleward oceanic heat transport towards the Arctic is maximal (Fig. 6a, c). These changes are characterized by an overall sea-ice melting between the Greenland/Barents Sea and Labrador Sea along the East Greenland coast and result from the enhanced poleward heat transport by the AMOC and subpolar gyre (Fig. 6a). The persistent sea-ice melting over the North-Atlantic coincides with pronounced tendency towards negative NAO that accelerates the multidecadal transition of the NAO from the positive into negative phase. Since AGCMs tend to simulate negative NAO or NAM in response to the Arctic sea-ice melting56, the sea-ice melting may contribute to the acceleration of the multidecadal phase transition of NAO and NAM after phase 4.

Fig. 6: Impact of the coupled multidecadal stratosphere/troposphere/ocean oscillation on the sea-ice, surface temperature, and precipitation.
figure 6

a The standardized averaged evolution (using phase-composite) of the reconstructed sea-ice concentration (averaged over 60 W–50 E and 50–90 N, in yellow), reconstructed NAOI (red), poleward oceanic heat transport in the Atlantic basin (averaged between 55 N and 75 N, in solid blue), associated contribution from the overturning circulation (blue dotted) and gyre (blue dot dashed). b, c The sea-ice concentration (in %) in b at phase 1 (with positive NAO) and (c) phase 4 (with enhanced sea-ice melting), respectively. d, e Surface temperature changes (in °C) in d at phase 1 (associated with positive NAO) and e at phase 4 (associated with enhanced Arctic sea-ice melting). f, g as in d and e, but for the precipitation flux (in mg m2 s−1).

Consistent with the well-known NAO impacts on continental climates32,33, the positive NAO in phase 1–2 is associated with warming and enhanced precipitation over North Eurasia and cooling over the Northeast America, Northwest Atlantic, North Africa, and Middle East, whereas Southern Europe and the Mediterranean region feature reduced precipitation (Fig. 6d, f). The regional climate anomalies are reversed during the negative NAO phase (Supplementary Fig. 4). The highest amplification of the high-latitude near-surface atmospheric warming occurs between phases 4 and 6, at the peak of sea-ice melting due to the poleward oceanic heat transport (Fig. 6e). This Arctic warming is like the Arctic amplification and is maximal over land, where temperature-albedo-feedback associated with continental snow-cover may also play an important role. A near-surface atmospheric cooling of similar amplitude occurs in phases 12-13 (Supplementary Fig. 4), when sea-ice cover and Atlantic cooling through the reduced poleward oceanic heat transport are maximal. The ocean-atmosphere-land-ice including sea-ice feedback can thus be important in amplifying the AMV-induced temperature changes supporting the result shown in19. The Arctic warming amplification is associated with enhanced (reduced) precipitation in the North Eurasia (Mediterranean) (Fig. 6). These conditions are reversed for the cooling amplification (Supplementary Fig. 4).

Discussion

The coherent multidecadal changes in different NH-climate components can be understood within the framework of a damped coupled stratosphere/troposphere/ocean-oscillation including Arctic sea-ice (Supplementary Fig. 8). The existence and the understanding of such coupled quasi-periodic oscillations can enhance the long-term prediction of climate-system components involved and help improving Atmosphere-Ocean General Circulation Models (AOGCM) and developing statistical prediction methods. A lagged-multi-linear regression model (Methods) shows that the performance of the statistical model in simulating the observed detrended climate variability in projection mode can be improved, if the coupled co-variability of the ocean, atmosphere, and sea-ice is considered (Supplementary Fig. 5 and 6). While the observed multidecadal variability of wintertime sea-ice and surface temperature depends on both long-term climate change and coupled multidecadal variability, the changes of the NAO and associated extremes depend mainly on the coupled AMV. During multidecadal positive (negative) phases, the occurrence of the extreme positive (negative) events (based on 95% (0.5%) percentile) is enhanced (Fig. 7). In particular, the multidecadal NAO variability has led to persistent extreme positive NAO-values, stratospheric cooling (Fig. 2d) (manifested in extreme low Sudden-Stratospheric-Warming (SSW)-frequency57), persistent and enhanced north Eurasian warming3 (Fig. 7) and extreme precipitation decrease in the Mediterranean area9,10 in the 1990s. These have affected the economy of the North-African and Southern-European countries58. The end of the last century was also marked by AMOC-strengthening59 and multidecadal melting of the wintertime Arctic sea-ice over the North-Atlantic region that lasted until the beginning of the 21st-century and accelerated the long-term Arctic sea-ice melting and global warming (Fig. 1a, b, Fig. 7). All these coherent multidecadal changes are physically consistent within the damped coupled stratosphere-troposphere-ocean-oscillation framework (Supplementary Fig. 8).

Fig. 7: Implication of the coupled atmosphere/ocean variability for the near-future climate.
figure 7

The non-normalized wintertime (JFM) indices of historical (black line) and ensemble projection (Methods) computed as: a averaged North Atlantic SST (over the same region as used for the AMV index in °C), b NAO index (standardized), c averaged Arctic sea-ice concentration (60° W–50° E and 50–90° N, in %), d averaged Eurasian surface temperature (10° W–130° E and 45–85° N, in °C), e global mean surface temperature (in °C). All projections were performed based on the observed coupled atmosphere/ocean variability using lagged multi-regression method (Methods). The future projections were first performed on detrended data and the trends (that were extrapolated into the future) were added afterward. This allowed us to understand the interaction between the multidecadal fluctuations and the long-term climate trends. Yellow lines represent predictions with sea-ice concentration index that was detrended using a linear trend, and the gray lines represent the predictions from the sea-ice that was detrended using a quadratic trend (Methods). The ensemble mean (red line) is computed from the individual members using both linear and parabolic approximations of sea-ice trend. The smooth black and red lines represent 50–70 year smoothing of the raw data. Horizontal green lines represent the 5 and 95% percentiles.

The multidecadal AMV-component reached recently the maximal warming (Fig. 1a, Fig. 7a) like the 1940s-to-1950s. This was followed by a negative NAO, North Atlantic cooling, Eurasian cooling and sea-ice extension, which resulted in long-term hiatus in wintertime global warming (1950s–1970s). In line with the life cycle of the coupled Atmosphere/Ocean oscillation (Fig. 3c, d, Fig. 4, Supplementary Fig. 8), the multidecadal NAO-component is currently moving into the negative phase (Figs. 1 and 7b), the SSWs are becoming more frequent57, the precipitation and associated extremes in the western Mediterranean are increasing12, and there is an indication of an ongoing weakening of the AMOC59,60. The introduced coupled oscillation and associated statistical climate projection suggests that the ongoing negative NAO-trend will continue (Fig. 7), which can lead to further increase (decrease) of the precipitation in the Southern (Northern) Europe, a weakening of the AMOC and subpolar Atlantic-gyre and therefore to multidecadal sea-ice extension and Atlantic cooling trend, like the period 1950s–1970s (Figs. 1b and 7). This can offset the long-term Arctic sea-ice melting, north of the Atlantic basin, leading to sea-ice melting hiatus (Fig. 1b and 7). The net result is NAO- and sea-ice-induced increase of the precipitation over the Mediterranean region and Eurasian cooling (Supplementary Fig. 4c–f) that can offset the Northern Hemispheric and global warming together with the North Atlantic cooling and lead to a warming hiatus like the one seen in the 1950s–1970s (Fig. 7). The impact of the NAO on the NH surface temperature via delayed AMV-response (according to the Hasselmann model) was previously used to construct an NAO-based linear prediction model for the NH surface temperature41. This model predicts decadal cooling and is therefore consistent with our prediction. The impact of the AMV on the NH-hemisphere and global temperature can be seen as broad impact of the AMV on climate system61. Such broad AMV impact can connect different mode of climate variability and different climate components of the climate system together through for example the stadium-wave climate regime47,48,49,50. Despite the consistent coherence among the different climate components in the model projection, the statistical model projects an enhanced decadal excursion. This decadal excursion can be seen as numerical artifact (Methods) that can result from the detrending method used before learning process, low accuracy of the observational data in the first half of last century, short learning record and possible limitation of statistical prediction. Unlike the NAO, whose multidecadal trends depend mainly on the multidecadal fluctuations, the surface temperature has shown a strong long-term anthropogenic warming and sea-ice has shown a strong long-term melting trend (Fig. 7). Because of these trends, the multidecadal warming and sea-ice melting have become much more extreme and the multidecadal cooling and sea-ice extension have become much less extreme (Fig. 7a–e).

The periodicity and the coupled character of this oscillation can thus enhance the decadal-to-multidecadal climate predictability, benefiting risk and impact assessments on the same time scales. State-of-the-art AOGCMs agree on the high decadal predictive skill of the North-Atlantic SSTs. However, they still underestimate the decadal predictive skill of the NAO/NAM, which can be partially reached through very large ensemble simulations62. This shortcoming can be caused by the weak representation of the predictable signal in the model or/and enhanced background noise due to model bias and internal variability, which is known as signal-to-noise paradox63. The instantaneous standardized NAO index in the introduced coupled oscillation remains in general smaller compared to observations, which can be seen as indication of a reduced signal-to-noise ratio in the MPI-ESM model. This can result from the still existing model biases such as the North Atlantic SST bias. An increase of the decadal prediction skill of the NAO in large ensemble prediction using CESM-DPLE-model has been attributed to the AMV64 indicating the ability of this model to simulate the coupled NAO/AMV variability in prediction mode. However, the reduced signal-to-noise ratio can strongly mask the multidecadal quasi-periodic variability making its detection and prediction in model-simulations difficult. Our work shows that the lack of the multidecadal coupling among atmosphere, Atlantic SST and Arctic sea-ice can reduce the skill of the predictable components of the climate system (Supplementary Figs. 5 and 6). More specifically, the AOGCM models that have difficulties in reproducing the observed NAO-AMV-link27 and the associated possible impact on the Arctic sea-ice can also show reduced long-term predictive skill for the large-scale atmospheric circulation. Several reasons for such shortcomings have been proposed, e.g. lack of stratospheric resolution and the Atlantic cold SST-bias14.

The simple delayed oscillator model described in previous work40 is based on the annual AMV and NAO indices. In this model, the negative feedback results mainly from the delayed response of the NAT to the AMOC and the associated NAO-response. In this study, we focused mainly on the wintertime atmospheric circulation where the NAO or NAM can be forced by the AMV. The wintertime NAO-response to the AMV can be seen as additional negative feedback that precedes the delayed negative feedback associated with the NAT-response40. The NAT varies in phase with the NAO and can be therefore partially seen as instantaneous thermodynamical response to the NAO44. This instantaneous response can enhance the persistence of the NAO-response to the AMV through the positive NAO/NAT feedback and the oceanic thermal inertia. The damped couple quasi-periodic multidecadal variability of the NAO and AMV/AMOC has been shown also in control model simulations without transient external forcing40,51. Such damped character can be controlled by the long-lasting damping timescale of the ocean40. The current study shows that a damped coupled atmosphere/ocean oscillation can also occur in transient forced model experiment, which mimics the real world.

Finally, in the assessment of ozone depletion, the possible impact of multidecadal ocean/atmosphere-coupled-variability on the ozone trends, which have been mainly attributed to anthropogenic forcing65, has been overlooked. Since both stratospheric temperature and circulation changes can strongly alter the ozone concentration65 and, as shown here, these changes can result from the coupled stratosphere/troposphere-ocean multidecadal variability, a reassessment of multidecadal ozone variability is needed to quantify the contribution of atmosphere/ocean-coupled processes.

Methods

Observations and reanalysis

In this work, we use the Hadley center SST and SLP66 from 1870 to 2019 and NOAA Twentieth Century Reanalysis (20CR), version V2c67, to identify some characteristic of the multidecadal variability and associated coherent changes in NH climate system. The 20CR-data product provides global atmospheric circulation data set spanning 1850–2014 and it assimilates only surface pressure from monthly SST and sea-ice distributions as boundary conditions with a ‘deterministic’ Ensemble Kalman Filter. The 20CR-data can thus not be used to evaluate the multidecadal variability of the stratosphere. To assess the stratospheric variability, we use the standard NCEP-reanalysis, which assimilates both surface and free atmospheric data including the stratosphere68, but the time record is shorter than in the 20CR and Hadley center data. The 20CR data69 uses sea-ice boundary conditions from the COBE-SST270, which are available until 2019. For surface temperature we used HadCRUT (1870–2019), for 1000 hPa geopotential height and precipitation we used NOAA-20CR (1870–2014). Statistical prediction based on the coupled variability is performed using the available observational data. The latest starting year for our prediction is 2020 for data available until 2019.

Model experiment

Observational data and stand-alone atmospheric and oceanic model experiments indicate that the AMV and stratosphere/troposphere-coupled NAM (or NAO) can maintain each other. However, the short observational records and the lack of oceanic data limits strongly the identification and the understanding of the oscillatory character in the coupled stratosphere/troposphere/ocean system. Thus, we used long-term forced atmosphere/ocean-coupled simulation with stratosphere-resolving model and applied MSSA-method71,72, which is designed to isolate spatial and temporal properties of the most relevant oscillatory mode of variability. In this work, we used the last millennium simulation, past1000-r1i1p1, with Max Planck Institute Earth System Model for Paleo applications (MPI-ESM-P)73,74. The MPI-ESM consists of the atmosphere general circulation model (AGCM) ECHAM675 coupled to the global ocean/sea-ice model (Max Planck Institute Ocean Model, MPI-OM)76,77. ECHAM6 is run by default in a stratosphere-resolving configuration, which extends to 0.01 hPa (~80 km). ECHAM6 is run at a horizontal resolution of spectral truncation T63 (1.875°) and 47 vertical levels. The MPI-OM model features a conformal mapping grid with nominal 1.5° resolution and 40 vertical levels (GR1.5L40). In MPI-ESM-P, land-cover change maps78 and tables of annual values for eccentricity, obliquity, and perihelion are used. The past1000 simulation was performed according to the protocol of the Paleoclimate Modeling Intercomparison Project, phase 3 (PMIP3)79. Volcanic forcing is implemented using volcanic aerosol optical properties reconstructed by Crowley and Unterman80. Changes in global land-cover and agricultural areas are reconstructed from Pongratz et al.78. The solar radiation is the combined reconstruction of Vieira et al.81 over the Holocene and the recommended solar forcing for the CMIP5 20th-century (1850–2005) simulations. The 11-year solar cycle of varying amplitude is imposed artificially over the pre-industrial period79. Reconstructed total-solar-irradiance TSI-data between 850 and1849 were adjusted to the spectral bands of ECHAM6’s radiation scheme. The monthly ozone concertation between 850 and 1849 were calculated using the monthly ozone climatology between 1850–1860 from AC&C/SPARC data base and from the ozone dependency on solar irradiance. This dependency is calculated from regression coefficients between historical ozone concentrations between 1850–2005 and the annual 180.5 nm solar flux. The boundary conditions follow the CMIP5 protocol but using land-cover changes from Pongratz et al. for the period 1850–200578.

Indices

The observed NH-multidecadal fluctuation were extracted by defining various indices and using a wavelet bandpass filter for time scales between 50 and 70 year82. For the model, the indices are computed from the reconstructed multidecadal atmosphere/ocean-coupled MSSA-variability. All indices are divided by the corresponding standard deviation of the row data so that the values reflect the associated explained variance. Such standardization enables also the comparison among the indices from different variables in both observation and model simulation. The climate projection in Fig. 7 has been performed with non-standardized predictand, but the predictors were all standardized.

The observed wintertime (JFM) large-scale North-Atlantic SST variability is indexed by the detrended SST anomalies averaged over the region 0–60° N and 75–7.5° W for winter months. This index is called the Atlantic Multidecadal Variability Index AMVI83. The standard AMVI includes (in general) the direct radiative impact of external forcing (i.e., solar radiation, volcanic eruptions, and anthropogenic aerosols). To enhance the signature of internal variability associated with ocean circulation in the AMV, we define the Local AMVI (LAMVI) by removing the global temperature anomalies from the Atlantic SST at each month84,85. The modeled AMVI is computed from the atmosphere/ocean coupled variability based on the joint MSSA71,72. The modeled AMVI contains only the part of the AMV, which is coupled to both ocean and atmosphere and is therefore more comparable to the LAMVI than the standard AMVI.

The NAOI is computed by projecting the monthly anomalies of geopotential height at 1000 hPa (Z1000 hPa) level on first Empirical Orthogonal Function (EOF) of the wintertime (JFM) Z1000 hPa anomalies over the North-Atlantic sector (20–80° N; 90° W–40° E)33. The NAM indices (NAMI) are computed for different atmospheric levels. The NAMI is defined at each level as the projection of the monthly anomalies of geopotential height on its first EOF for the NH (20–90° N)32. To identify the observed stratospheric coupling with the surface NAO or NAM on the multidecadal timescale, we applied the multidecadal wavelet filtering on the vertical profile of the seasonal NAM evolution. The seasonal downward NAM propagation in the model is reconstructed from the multidecadal component of the coupled MSSA analysis. In the seasonal NAM evolution, the individual months represent the center of running 3 month average.

To understand the mechanisms of the coupled stratosphere/troposphere/ocean oscillations, we defined several regionally averaged indices (see figure captions) describing the fluctuations of the ocean circulation, heat and moisture fluxes and Arctic sea-ice. The definitions of these additional indices are provided in the text/figure caption where relevant.

Multichannel Singular Spectrum Analysis (MSSA) and phase-composite

The MSSA is a multichannel (or multi-grid) extension of the singular spectrum analysis (SSA)71,72, which is a spectral method for univariate time series. SSA and MSSA are also used to enhance the signal-to-noise ratio and extract different signals, including oscillatory modes. Unlike classical spectral methods, which use prescribed functions (for example Fourier transform), SSA is a data adaptive technique based on eigenvector analysis using lag-covariance matrix of the time series. In this method, a set of basis vectors or modes is objectively identified to uniquely describe the data. Considering a time series \(X_t\left( {t = 1,..,N} \right)\), with N the number of time steps. To perform SSA-analysis of Xt, we first embed Xt in vector space of dimension M < N (seen as window) by setting: \({{{\mathbf{Y}}}} = (X_{t^\prime },X_{t^\prime + 1},..,X_{t^\prime + l - 1},..,X_{t^\prime + M - 1})^T\) for each t. Afterward we compute the lag-covariance matrix \({{{\mathbf{C}}}} = \frac{1}{{N - M + 1}}{{{\mathbf{Y}}}}^{{{\mathbf{T}}}}{{{\mathbf{Y}}}},{{{\mathrm{Finally}}}}\,{{{\mathrm{the}}}}\,{{{\mathrm{lag}}}}\,{{{\mathrm{convariance}}}}\,{{{\mathrm{matric}}}}\,{{{\mathrm{is}}}}\,{{{\mathrm{diagonalized}}}}\,{{{\mathrm{into}}}}\,{{{\mathbf{\Lambda }}}} = {{{\mathbf{E}}}}^T\,{{{\mathbf{C}}}}\,{{{\mathbf{E}}}}\), where Λ is M × M diagonal matrix having the eigenvalues in decreasing order in the diagonal and E is M × M matrix having the eigenvectors ek as columns (where k = 1,.., M are the modes). Compared to the standard EOF the eigenvectors are in time dimension and are called temporal EOF (T-EOF) and the principal component (PC) are called T-PC and are computed for each mode k as \(\alpha _{t\prime k} = \mathop {\sum}\nolimits_{i = 1}^{i = M} {X_{t\prime + i}e_{k,i}}\). The signal of interest can be reconstructed from corresponding mode of eigenvalue decomposition, and an oscillatory mode can be reconstructed only from two neighboring modes k and k+1 having T-PCs nearly in phase quadrature and same oscillation-period. Here we are interested in multidecadal time scales. The reconstructed signal for two neighboring modes (T-RCs) is calculated as \(r_{k,k + 1} = \frac{1}{{M_t}}\left( {\mathop {\sum}\nolimits_{i = 1}^{i = M} {\alpha _{t - i,k}e_{k,i}} + \mathop {\sum}\nolimits_{i = 1}^{i = M} {\alpha _{t - i,k + 1}e_{k + 1,i}} } \right)\), where Mt is a normalization factor. For the MSSA, the time series are given in each channel or grid-point c (c = 1,.., L) by Xt,c and the embedded vector is constructed for time series of each channel as \({{{\mathbf{Y}}}} = (X_{t\prime ,1},X_{t\prime + 1,1},..,X_{t\prime + M - 1,1},X_{t\prime ,2},X_{t\prime + 1,2},..,X_{t\prime + M - 1,2},..,X_{t\prime ,L},X_{t\prime + 1,L},..,X_{t\prime + M - 1,L})^T\). The eigenvalue decomposition of the lag-covariance matrix \({{{\mathbf{C}}}} = \frac{1}{{N - M + 1}} = {{{\mathbf{Y}}}}^T{{{\mathbf{Y}}}}\), and the reconstruction of the oscillatory signal is performed in similar way as in SSA. The eigenvectors of the lag-covariance matrix contain both spatial and temporal pattern and are called spatio-temporal EOFs (ST-EOFs). Similar to SSA, we can determine ST-PCs and ST-RCs. As proposed by Vautard et al. 1992, the SSA and MSSA are successful in analyzing the oscillations with periods between M/5 and M, where M is the window used. Our interest is on multidecadal time scales. We therefore chose the window of M = 70 year. Our results are not sensitive to the choice of M.

We apply MSSA to data from a long-term stratosphere-resolving coupled atmosphere/ocean GCM simulation. This is to better understand the multidecadal relationship among the stratosphere, troposphere and AMV as shown previously in the observations, and stand-alone atmosphere and ocean models. We use a last millennium run performed with the MPI-ESM model. Since models tend in general to simulate different AMV periodicity than observed and since the reconstructed forcing is imperfect, we don’t expect to synchronize the coupled multidecadal atmosphere/ocean variability with observations. The focus is on understanding the coherent coupled stratosphere/troposphere/ocean variability in the framework of a damped delayed oscillation and understanding the dynamical processes controlling this oscillation. In order to extract the coupled variability, we first apply the MSSA on the combined coupled field from winter geopotential height in the stratosphere (20 hPa, NDJ) and troposphere (1000 hPa, JFM), Atlantic SST (JFM) and yearly Atlantic meridional overturning mass stream function. The spatial degrees of freedom and the noise are reduced by applying the MSSA on the first 20 PCs obtained from the standard EOF analysis. The data are normalized and concatenated before performing the EOF and MSSA analysis. In order to understand the mechanisms of the coupled atmosphere/ocean oscillation, we rerun the MSSA each time again after adding one additional variable like the oceanic and turbulent heat-fluxes, sea-ice, vertical structure of the NAM, different terms of the atmospheric momentum budgets and other variables. The results of the MSSA does not change by adding one variable to the main four initial variables. Because of the memory limit, the data are interpolated to a coarser resolution before performing the joined EOF and MSSA analysis.

To understand the mechanisms maintaining the coupled atmosphere/ocean oscillation, we apply the phase-composite72. For this purpose, we reconstructed the multidecadal PC = \(c_{k,k + 1}\) (from the second and third MSSA modes, k = 1) of the leading channel of the 20 standard EOFs (used as impute for MSSA). This PC is sinusoidal and has the same period of the total reconstructed signal \(r_{k,k + 1}\) (t) (from the 20 EOFs). It explains the highest amount of variance and can therefore be used for composite analysis72. The index \(c_{k,k + 1}(t)\) is first normalized into b(t) using standard deviation and complexified as \(z\left( t \right) = b\left( t \right) + b\prime \left( t \right)i = A(t)e^{i\theta (t)}\), where the imaginary part b'(t) is the time derivative of b(t), \(A\left( t \right) = \sqrt {b(t)^2 + b\prime (t)^2}\) is the amplitude and \(\theta \left( t \right) = arctg\left( {\frac{{b\prime (t)}}{{b(t)}}} \right)\) is the phase of z(t). From the phase time series θ(t), we selected 16 phases \({{{\mathrm{in}}}}\,{{{\mathrm{the}}}}\,{{{\mathrm{phase}}}}\,{{{\mathrm{intervals}}}}\left[ {0,2\pi } \right]\) according to (m1) π/8 < θ(t) < mπ/8, where m varies between 1 and 16. We call the average of the reconstruction of each variable over each of these phases the phase-composite. The evolution of each reconstructed variable over these phases will give the mean evolution over one period of the oscillation.

Dynamical analysis in the atmosphere

To understand the atmospheric circulation changes associated with NAM and NAO, we use the Transformed Eulerian Mean (TEM) formalism of the zonally averaged quasi-geostrophic momentum budget. In the TEM, the wind changes or tendency can result from the interplay between the resolved wave forcing, non-resolved subscale processes, and residual Coriolis forcing86,87. The resolved wave forcing is proportional to the divergence of the Eliassen-Palm (EP) flux vector. A wave propagation into (from) the jet is characterized by wave convergence (divergence) near the jet and thus jet deceleration (acceleration). The non-resolved forcing includes the parametrized subscale processes like friction near the bottom, orographic and non-orographic gravity wave in the free atmosphere. The residual Coriolis forcing, which is proportional to the meridional residual velocity, tends to balance the resolved and non-resolved forcing and it accelerates (decelerates) the westerly wind by poleward (equatorward) mass-transport. Its sign can also be used to describe the atmospheric (residual) meridional overturning circulation.

Statistical climate projection

The future projections of different indices in Fig. 7 were computed using multi-linear regression model (via scikit-learn88) in the following way. First, the raw data were detrended (for details see below). Second, we took time series of NAO, AMO, sea-ice, and the variable (if applicable) we are projecting and found the lag of maximal correlations between the variable and NAOI, AMVI, and sea-ice (index). We do this so that all variables that we input into the statistical model are covarying (i.e., maximum correlation between them is then shifted to lag 0). Not that for the projection of the NAOI, AMVI, and sea-ice, we used only these three variables as predictors. To project one year ahead (i.e., for lead time +1), we input ten past time steps in the multi-linear regression model (i.e., lags −9, −8, −7, −6, −5, −4, −3, −2, −1, 0 years instead of using only lag 0 years). For larger lead times, we lagged the data by the lead time (i.e., we additionally lagged the above-listed lags by 9 years if lead time was 10 years). Note that we construct a different multi-linear regression model for each lead time (for up to 35-year lead time). Once we computed the projection for all lead times considered, we computed correlation skill score to test the skill of the statistical model in reproducing the observed data. Note that the training and testing datasets were the same due to data constraints, hence the skill is likely higher than expected. We also used the lagged-ensemble method for generating projections beyond year 2020 by using projections starting between years 2010–2019 following the same method as described above. This gave us a 10-member ensemble prediction for up to year ~2045. To test the importance of the ocean/atmosphere/sea-ice coupling in the projection mode, we computed the projection correlation skill using the following predictors: (1) predicted index only, (2) NAOI, AMVI and predicted index without sea-ice, and (3) NAOI, AMVI, predicted index and sea-ice. To understand the superposition between the long-term trends and multidecadal fluctuations, the long-term trend was extrapolated into the future. Note that because of the lack of the sea-ice measurements in the Arctic, the sea-ice shows almost no trend in the early 20th century. It is therefore difficult to separate the long-term trend from the multidecadal fluctuations, since most of the trend started at the end of the last century. Removing linear trend in sea-ice leads to a very strong multidecadal fluctuation, and removing the quadratic trend leads to a very weak multidecadal fluctuation in both historical and predicted sea-ice index. This sensitivity to the detrending method results in a decadal excessive excursion of the future sea-ice concertation when the linear trend is removed. All predictions (of all variables) are performed with linearly and quadratically detrended sea-ice separately and the ensemble means are computed from both ensemble (ultimately leading to 20 ensemble members per variable prediction between 2020 and 2045). All other variables (i.e., except for the sea-ice) were detrended linearly. Since statistical predictions were developed to mainly make predictions with possible future scenarios that consider the coupled-mode dynamics, we prefer using the term projection instead of prediction. Since the stratospheric and tropospheric decadal-to-multidecadal variability are in phase (Figs. 3b and 4 and Supplementary Fig. 1), adding stratosphere as predictor in the multi-regression model can lead to overfitting. We used therefore only tropospheric NAO, since the long-term stratospheric measurements are also not available.