Interactions between the stratospheric polar vortex and Atlantic circulation on seasonal to multi-decadal timescales

Variations in the strength of the Northern Hemisphere winter polar stratospheric vortex can influence surface variability in the Atlantic sector. Disruptions of the vortex, known as sudden stratospheric warmings (SSWs), are associated with an equatorward shift and deceleration of the North Atlantic jet stream, negative phases of the North Atlantic Oscillation, and cold snaps over Eurasia and North America. Despite clear influences at the surface on sub-seasonal timescales, how stratospheric vortex variability interacts with ocean circulation on decadal to multi-decadal timescales is less well understood. In this study, we use a 1000 year preindustrial control simulation of the UK Earth System Model to study such interactions, using a wavelet analysis technique to examine non-stationary periodic signals in the vortex and ocean. We find that intervals which exhibit persistent anomalous vortex behaviour lead to oscillatory responses in the Atlantic Meridional Overturning Circulation (AMOC). The origin of these responses appears to be highly non-stationary, with spectral power in vortex variability at periods of 30 and 50 years. In contrast, AMOC variations on longer timescales (near 90-year periods) are found to lead to a vortex response through a pathway involving the equatorial Pacific and quasi-biennial oscillation. Using the relationship between persistent vortex behaviour and the AMOC response established in the model, we use regression analysis to estimate the potential contribution of the 8-year SSW hiatus interval in the 1990s to the recent negative trend in AMOC observations. The result suggests that approximately 30 % of the trend may have been caused by the SSW hiatus.

1 Introduction

Variations in the strength of the Northern Hemisphere (NH) stratospheric polar vortex associated with sudden stratospheric warming (SSW) events is the single largest source of interannual variability in the NH winter stratosphere. Extreme disruptions of the vortex during SSWs also represent a key mechanism for stratosphere–troposphere coupling and are widely acknowledged to lead to anomalies in NH midlatitude surface climate, particularly in the North Atlantic sector (Baldwin and Dunkerton, 2001). Equally, the absence of an SSW, when the vortex is relatively undisturbed by waves propagating from the troposphere so that the vortex is unusually strong and cold, also has significant impacts on surface weather (Shaw and Perlwitz, 2013; Lawrence et al., 2020). An understanding of SSW dynamics and vortex variability is, therefore, important for seasonal to sub-seasonal forecasting of surface weather (Domeisen et al., 2020a, b), as they provide a significant source of predictive skill. The majority of studies that have examined the associations between the stratospheric polar vortex and surface anomalies have considered their in-season impact. However, the NH winter stratosphere also exhibits variability on decadal to multi-decadal timescales (Dimdore-Miles et al., 2021). The interaction between this low-frequency variability in vortex strength and variations in surface climate on similar timescales is not fully understood.

The impact of stratospheric polar vortex (hereafter referred to as the vortex) variations on surface climate variability was highlighted in the seminal work of Baldwin and Dunkerton (2001), which demonstrates a steady downward propagation of Northern Annular Mode (NAM) anomalies from the middle stratosphere to the surface, following strong and weak vortex events in the ERA40 reanalysis dataset. The anomalies associated with a weak (strong) vortex were shown to persist at the surface intermittently for approximately 60 d after the event and projected significantly onto the negative (positive) pattern of the North Atlantic Oscillation (NAO). Subsequent modelling and observational studies have corroborated these results. Domeisen (2019) show, using ERA40 and ERA-Interim datasets, that approximately two-thirds of SSW events are followed either by a switch from positive to negative NAO or a persistent negative NAO pattern. Charlton-Perez et al. (2018) consider this coupling from the perspective of tropospheric weather regimes and find, in both ERA-Interim and the ECMWF Integrated Forecasting System, a 40 %–60 % increase in probability of transition to a negative phase of the NAO given a 1 standard deviation reduction in the polar vortex strength. This SSW influence on the NAO is relatively well represented in general circulation models (GCMs; Baldwin et al., 2021), and idealised modelling studies are also able to show a direct downward influence of events on the NAO (White et al., 2020; Gerber et al., 2009), although some studies have noted that simple models tend to overrepresent the persistence of surface anomalies (Gerber et al., 2008a, b),

The influence of SSW events on negative NAO phase probability has subsequently been shown to influence other features of NH midlatitude climate. Thompson et al. (2002) show, in the National Centers for Environmental Prediction (NCEP) and National Center for Atmospheric Research (NCAR) reanalysis, that SSWs are followed by a 60 d interval of anomalously low surface temperatures in eastern North America, northern Europe, and eastern Asia. Subsequent observational studies find similar response patterns to SSWs (Kolstad et al., 2010; King et al., 2019; Lehtonen and Karpechko, 2016), and the effect can also be seen in GCM simulations (Tomassini et al., 2012; Lehtonen and Karpechko, 2016). Further impacts of SSWs on tropospheric circulation include an equatorward shift and deceleration of the North Atlantic eddy-driven jet stream (Hitchcock and Simpson, 2014; Maycock et al., 2020). Links with persistent blocking events are also shown in observations and modelling studies (Davini et al., 2014; Vial et al., 2013), although Taguchi (2008) finds no significant link between the phenomena.

While the in-season influence of SSWs on tropospheric circulation and surface variability is discussed extensively in previous work, their coupling with modes on longer timescales, such as ocean variability, is less well understood. One of the primary features of Atlantic Ocean variability is the Atlantic Meridional Overturning Circulation (AMOC), which consists of a northward transfer of warm, saline water that occurs in the top 2 km of the Atlantic Ocean (the upper cell), accompanied by a corresponding return flow of southward transport at lower depths (the lower cell; Kuhlbrodt et al., 2007; Xu et al., 2014; Buckley and Marshall, 2016). The strength of the AMOC varies significantly on decadal–centennial timescales (Delworth et al., 1993; Biastoch et al., 2008; Tulloch and Marshall, 2012; Menary et al., 2012) and is thought to be a key driver of North American and European surface variability via modulation of the Atlantic sea surface temperatures (SSTs) and heat transport (Knight et al., 2005; Delworth and Mann, 2000; Frierson et al., 2013; Frankignoul et al., 2013). It has been shown to influence a diverse range of features, such as European summertime temperatures (Sutton and Hodson, 2005), biogeochemical conditions in the northwest Atlantic (Lavoie et al., 2019), and abrupt climate shifts in palaeoclimate records (Alley, 2007; Cheng et al., 2009).

Multiple potential drivers of AMOC variability at different timescales have been studied extensively in both observations and modelling studies. On intra-annual to inter-annual timescales, variability in the AMOC has been closely associated with wind variations over the North Atlantic region through Ekman transport anomalies or wind stress curl forcing (Wang et al., 2019; McCarthy et al., 2012; Mielke et al., 2013; Yang, 2015). On inter-annual to decadal timescales, AMOC variability has been associated with buoyancy anomalies in the subpolar region, particularly in the Labrador Sea (Delworth et al., 1993; Medhaug et al., 2012). This mechanism is linked to variability in mixed layer depth and the occurrence of deep convection over the same region, particularly in NH winter (Böning et al., 2006; Biastoch et al., 2008; Robson et al., 2012; Wang et al., 2015). Mixed-layer anomalies in the Labrador Sea are an indication of the strength of deep convection in this region, which has been shown to be associated with AMOC variations in modelling studies (Eden and Willebrand, 2001; Eden and Jung, 2001) and observations (Latif and Keenlyside, 2011).

An association between stratospheric polar vortex variability and the AMOC on decadal timescales has been previously investigated (Reichler et al., 2012; Schimanke et al., 2011), but the mechanism of its influence remains unclear. For example, Reichler et al. (2012) examine the response of the AMOC to strong and weak polar vortex events and show a lagged, oscillatory response in the AMOC. They propose a pathway involving alterations of wind stress and ocean–atmosphere heat flux anomalies in the western Atlantic due to the changed NAO patterns following the vortex events. The effect is prominent in a preindustrial (PI) control of a single model (Geophysical Fluid Dynamics Laboratory, GFDL, CM2.1) and, to some extent, in a suite of CMIP5 models. An impact of long-term changes in the NAO on the strength of the AMOC is supported by a number of studies (Visbeck et al., 1998; Delworth and Dixon, 2000; Delworth and Greatbatch, 2000; Eden and Willebrand, 2001; Lohmann et al., 2009; Robson et al., 2012). Most recently, Delworth and Zeng (2016) used a set of idealised GCM experiments in which they impose a perpetual ocean–atmosphere heat flux pattern associated with different NAO phases. They find significantly different AMOC mean states, depending on the imposed pattern (a stronger AMOC under positive NAO flux conditions than a control simulation). Haase et al. (2018) also analysed the in-season influence of SSW events on the NAO and ocean–atmosphere heat fluxes that then impacts the strength of deep convection in the North Atlantic, using the Community Earth System Model (CESM1) Whole Atmosphere Community Climate Model (WACCM). The study notes the presence of an anomalously shallow mixed layer depth in the Labrador Sea following an SSW event.

A key result from Reichler et al. (2012) is the decadal modulation of the SSW–AMOC co-variability. However, decadal to multi-decadal variability in vortex strength is not well understood. Some studies have focused on potential polar vortex impacts at the surface but suffer from low statistical significance due to the short observational record. Garfinkel et al. (2017, 2015) and Cohen et al. (2009) link decadal fluctuations in vortex strength with modulation of the global warming signal in Eurasian surface temperature in both reanalyses and CMIP3 models. Schimanke et al. (2011) demonstrate multi-decadal signals in SSW occurrence in a multi-century GCM simulation and propose an influence of these signals on similar period variability in Eurasian snow cover and Atlantic SSTs. However, results from this study are difficult to interpret, as the GCM used (EGMAM – ECHO‐G with Middle Atmosphere Model) exhibits significant bias in its vortex representation, with a mean SSW rate of only two events per decade compared to six events per decade in most reanalyses (Ayarzagüena et al., 2019). They note that repeating their study with a more advanced model is required to corroborate their findings. Manzini et al. (2012) examined decadal fluctuations in SSW events in a 260-year prescribed SST simulation of a GCM and analysed their impacts at the surface. They show that decadal vortex variability excites similar timescale variations in surface temperature and sea ice coverage between Greenland and Norway over the Atlantic sector. They propose this connection to be indicative of a delayed response of the AMOC to stratospheric forcing via the NAO, which subsequently influences northward Atlantic heat transfer and sea ice melt rates, as well as surface temperature anomalies.

More recently, Dimdore-Miles et al. (2021, henceforth referred to as DM21) examined long-term variability in the strength of the winter polar stratospheric vortex in a 1000 year preindustrial (PI) control simulation of the UK Earth System Model. They identified sequences of up to 11 consecutive years in which at least one SSW occurred every year (similar to that observed during the period 1998–2004) and also sequences of up to 12 consecutive years with a strong undisturbed vortex, as was observed during the 1990s (Manney et al., 2005; Pawson and Naujokat, 1999). They identified multi-decadal signals of 90-year periodicity in the latter that persisted for approximately 450 years of the 1000 years. Using wavelet and cross-spectral analysis, they associated this with a similar signal in the amplitude modulation of the quasi-biennial oscillation (QBO) and proposed that the vortex variability was driven by the QBO through modulation of the Holton–Tan relationship (Lu et al., 2008, 2014). However, the focus of that paper was primarily in the stratosphere, identifying and understanding the relationship between the long-term QBO and vortex variability. In this paper, we use the same 1000-year π control simulation to extend that study to examine the links between long-term vortex variability and the North Atlantic, including oceanic modes. We first examine the near-surface in-season response to extreme vortex events to demonstrate that the model is able to reproduce corresponding anomalies in mean sea level pressure (MSLP), ocean–atmosphere heat flux, and SSTs consistent with previous studies. We then explore the response on multi-decadal timescales, using a wavelet analysis method to examine non-stationary signals in the vortex and the AMOC variability. We find that multi-year intervals which exhibit the same type of persistent, anomalous vortex behaviour also exhibit co-variability with the AMOC across multiple timescales. The interactions and feedback mechanisms on these different timescales are explored in more detail, using a combination of lag/lead composite analysis and a wavelet spectral decomposition method. The amplitude and lag of the AMOC response to a multi-year period of strong polar vortex in the model is determined. This is then used to estimate the potential contribution of the observed 8 consecutive years with a very stable, undisturbed stratospheric vortex in the 1990s (Pawson and Naujokat, 1999) to the recent observed negative trend in the strength of the AMOC. The paper is structured as follows: Sect. 2 describes the GCM used in the investigation, the spectral analysis method (wavelet analysis), and the relevant climate indices. Section 3 presents results from the analysis. Section 4 provides a summary and discussion of the results.

2 Data and methods

2.1 Model configuration

This work utilises the same model as that analysed in DM21: The first version of the UK Earth System Model (henceforth referred to as UKESM). UKESM is a stratosphere resolving coupled ocean–atmosphere–land–sea ice model. It contains 85 vertical levels in the atmospheric domain simulated by the Global Atmosphere 7.1 component (GA7.1), 35 of which lie above 35 km in altitude (Walters et al., 2019; Williams et al., 2018). GA7.1 runs at N96 horizontal resolution (∼ 135 km at the Equator). Ocean circulation is simulated by GO6.0 (Storkey et al., 2018), which contains 75 vertical levels and runs at 1^∘ horizontal resolution. Additional interactive components simulating land surface, sea ice, and atmospheric chemistry processes are added via coupling with JULES (GL7.0), CICE (GSI8.1), and UK Chemistry and Aerosols (UKCA) models, respectively (Walters et al., 2019; Ridley et al., 2018; Mulcahy et al., 2018). The model's representation of the North Atlantic Climate system has been examined previously by Robson et al. (2020); this study shows that UKESM is able to realistically simulate key features, such as the AMOC. However, this study highlights the model's underrepresentation of the variability in the AMOC, an issue which is reported in other physical models (Roberts et al., 2014).

We utilise the same 1000-year PπI control simulation as examined in DM21. This simulation is spun-up to achieve model equilibrium before initialising, following the method outlined in Yool et al. (2020). The run is forced using CMIP6 preindustrial values for concentrations of major GHGs (greenhouse gases; global mean 284.317 ppm (parts per million) CO₂; 808.25 ppb (parts per billion) CH₄; 273.02 ppb N₂O). There are no volcanic eruptions in the simulation, but a background stratospheric volcanic aerosol level is imposed using climatological values between 1850 and 2014 estimated from satellite products and other model simulations (Menary et al., 2018). The simulation does not exhibit a solar cycle. We choose to analyse a π control simulation due to the length of integration performed (1000 years) compared to the timescales of stratospheric variations shown in DM21 (near 90-year variability). The length of this simulation provides a greater number of possible cycles of such variability available for analysis, compared to the use of a historical simulation. Furthermore, a π control allows us to analyse the internal variability in the stratosphere.

To estimate the contribution of stratospheric variations to recent observed AMOC trends, we also make use of observation-based datasets of the atmosphere and oceans. First, we utilise the reanalysis data from the European Centre for Medium-Range Weather Forecasts (ECMWF), with ERA5 (Hersbach et al., 2020) for assimilated observations for geopotential height (GPH) and MSLP fields downloaded from https://climate.copernicus.eu/climate-reanalysis (last access: 10 August 2021). Second, the RAPID array dataset, which provides time–depth profiles for the meridional overturning mass streamfunction in the Atlantic region at 26^∘ N (Moat et al., 2020), is used. These data are measured through a combination of ocean mooring and ship-based, satellite, and submarine telephone cable observations to estimate the strength of primary contributions to the meridional overturning circulation, i.e. Ekman transport (through wind stress), transport through the Florida Straits, and transport driven by the east–west density gradients between the American and African continents (McCarthy et al., 2015).

2.2 Wavelet analysis

We employ the same wavelet analysis method as DM21, based on Torrence and Compo (1998), to examine potential non-stationary spectral characteristics of time series data over a range of periods. This method is outlined in full in Sect. 2.3 of DM21 and is reproduced briefly here.

The wavelet transform of a 1-dimensional time series, x, of length N and uniform time step δt, is given by the convolution between the series and a wavelet function ψ, which has been scaled by the quantity s, as follows:

$$\begin{array}{}\text{(1)}& W_n\left(s\right)=\sum _{n^{\prime }=\mathrm{0}}^{N-\mathrm{1}}{x}_{n^{\prime }}{\mathit{\psi }}^{\ast }\left[(n^{\prime }-n){\displaystyle \frac{\mathit{\delta }t}{s}}\right],\end{array}$$

where s is the scale of the wavelet which indicates its period, and n is the time index. Varying s and translating along the time axis builds up a power spectrum for x in the time period domain given by $\left|W_n\right(s){|}^{\mathrm{2}}$. Following Torrence and Compo (1998), we vary the scale parameter in increasing powers of 2, such that s_j=s₀2^jδj and $j=\mathrm{0},\mathrm{1},\mathrm{\dots },J$, where j is the index for the wavelet scale, s₀ is the smallest resolvable scale for the series and J is the index corresponding to the largest scale given by the following:

$$\begin{array}{}\text{(2)}& J=\mathit{\delta }{j}^{-\mathrm{1}}lo{g}_{\mathrm{2}}\left({\displaystyle \frac{N\mathit{\delta }t}{s_{\mathrm{0}}}}\right).\end{array}$$

The translated and scaled wavelet evaluated at a given scale, s, has the following form:

$$\begin{array}{}\text{(3)}& {\mathit{\psi }}^{\ast }\left[(n^{\prime }-n){\displaystyle \frac{\mathit{\delta }t}{s}}\right]={\left({\displaystyle \frac{\mathit{\delta }t}{s}}\right)}^{\mathrm{1}/\mathrm{2}}{\mathit{\psi }}_{\mathrm{0}}\left[(n^{\prime }-n){\displaystyle \frac{\mathit{\delta }t}{s}}\right],\end{array}$$

and we select the form of the wavelet basis function ψ₀ following the recommendation of Torrence and Compo (1998) as a Morlet wavelet, an oscillatory function enveloped by a Gaussian process, which is expressed as follows:

$$\begin{array}{}\text{(4)}& {\mathit{\psi }}_{\mathrm{0}}\left(p\right)={\mathit{\pi }}^{-\mathrm{1}/\mathrm{4}}{e}^{i{\mathit{\omega }}_{\mathrm{0}}p}{e}^{\frac{p^{\mathrm{2}}}{\mathrm{2}}}.\end{array}$$

The advantages of using a Morlet wavelet for analysing signals in a climate time series are primarily due to its ability to resemble many of the features commonly observed in climate time series (e.g. changes in dominant period and amplitudes). The full argument for this choice is presented in Lau and Weng (1995) and DM21. To directly compare spectra of different indices, we normalise all time series by subtracting the mean and dividing by its standard deviation before performing the wavelet transform. In order to effectively compare spectral power across a range of frequencies, we additionally scale the power spectrum by dividing by the scale parameter (s_j, as defined above) associated with each frequency. This is done following the methodology of Liu et al. (2007), which shows that unscaled spectra exhibit a bias towards overestimated powers at longer periods and that an effective comparison across timescales is possible with such scaling. We also define a confidence interval for wavelet power observed for the series x by comparing the observed power to that produced by a time series modelled as a first-order autoregressive (AR1; red noise) process, r, given by the following:

$$\begin{array}{}\text{(5)}& r_n=\mathit{\alpha }{r}_{n-\mathrm{1}}+z_n,\end{array}$$

where α is the lag 1 autocorrelation of x, and z_n is a Gaussian white noise term. Torrence and Compo (1998) show that the power spectrum of this is χ² distributed and, therefore, can be used to define a 95 % confidence interval for any observed power.

2.3 Cross-wavelet spectra

We also define a measure of coincident spectral power between two time series, i.e. the cross-wavelet spectrum. This metric indicates whether two series exhibit power at the same time points and frequencies. The cross-wavelet spectrum of two time series x and y, with associated wavelet spectra $W_n^x$ and $W_n^y$, is calculated by projecting one spectrum onto the other as follows:

$$\begin{array}{}\text{(6)}& \left|W_n^{xy}\right(s\left)\right|=\left|W_n^{x\ast }\right(s\left)W_n^y\right(s\left)\right|,\end{array}$$

where $W_n^{x\ast }\left(s\right)$ is the complex conjugate of the wavelet power spectrum of x (Grinsted et al., 2004). The complex argument of $W_n^{xy}\left(s\right)$ gives the local phase difference between signals in x and y in frequency–time space. The phase relationship between the two time series can be represented by a vector that subtends an angle representing the phase difference in radians. On all plots of the cross-spectra, the arrows to the right (left) denote signals in the two series which are in phase and correlated (anti-correlated). Vertical arrows indicate a phase relationship of $\frac{\mathit{\pi }}{\mathrm{2}}$ between the time series, so that the evolution of one is correlated with the time rate of change of the other. As for individual power spectra, we define a confidence interval for which the cross-power of a larger amplitude is deemed significant (>95 % confidence interval) by comparing the power exhibited by an actual series with a theoretical red noise process. The cross-power of two such AR1 processes is theoretically distributed, such that the probability of obtaining cross-power greater than a set of red-noise processes is as follows:

$$\begin{array}{}\text{(7)}& D\left({\displaystyle \frac{\left|W_n^{xy}\right(s\left)\right|}{{\mathit{\sigma }}_x{\mathit{\sigma }}_y}}<p\right)={\displaystyle \frac{Z_{\mathit{\nu }}\left(p\right)}{\mathit{\nu }}}\sqrt{P_k^x{P}_k^y},\end{array}$$

where σ denotes the standard deviation of the time series, Z is the confidence interval defined by p (Z=3.999 for 95 % confidence), ν is the degrees of freedom for a real wavelet spectrum (ν=2), and $P_k^x$ is the theoretical Fourier spectrum of the AR1 process for a given wavenumber k and is given by the following:

$$\begin{array}{}\text{(8)}& P_k={\displaystyle \frac{\mathrm{1}-{\mathit{\alpha }}^{\mathrm{2}}}{|\mathrm{1}-\mathit{\alpha }{e}^{\mathrm{2}i\mathit{\pi }k}{|}^{\mathrm{2}}}}.\end{array}$$

2.4 Model diagnostics

We utilise the Northern Annular Mode (NAM) as a metric for the strength of the vortex as used by Baldwin and Dunkerton (2001) and numerous subsequent studies. The NAM is defined as the first principal component (PC) of the zonal mean deseasonalised geopotential height (GPH) field evaluated at latitudes north of 20^∘ N over the NH winter season (Dec–Mar) on a given pressure level. While GPH fields used to calculate the NAM are generally defined on 2D grids (in the latitude–longitude plane), we utilise a metric based on the zonal mean GPH following the methodology of Baldwin and Thompson (2009, this work shows that an effective calculation of the NAM can be carried out with such zonal mean quantities). To measure the vortex strength, we evaluate the NAM at 10 hPa, which is the pressure level used to identify major SSWs, and the resulting index is henceforth known as NAM₁₀. An individual vortex event (either strong or weak) is recorded when the daily NAM₁₀ crosses a threshold of +1.5 (strong) or −2 (weak). The day on which these thresholds are crossed is referred to as the central date. After this date, the NAM₁₀ must return to values between the thresholds for at least 10 consecutive days (which is the approximate radiative timescale of the mid-stratosphere) before another event can be recorded. The strong threshold value for events is chosen in accordance with the methodology of Baldwin and Dunkerton (2001), and the weak threshold selected such that it results in approximately the same rate of weak events (SSWs) as is reported in DM21, using the same simulation of UKESM but with a zonal wind definition of SSWs (0.54 events in winter).

We also use the NAM₁₀ to derive an index for the appearance of intervals of consecutive winters which show persistent vortex behaviour. The persistent NAM₁₀ interval index is defined as follows. First, the NAM₁₀ is averaged over each NH winter season (December–March); this gives a measure of mean vortex strength for each winter. Second, this index is smoothed using a Gaussian filter which is carried out through a convolution of the time series with a 1D Gaussian kernel in the time domain given by the following:

$$\begin{array}{}\text{(9)}& f(t,\mathit{\sigma })={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }{\mathit{\sigma }}^{\mathrm{2}}}}}{e}^{-\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{t}{\mathit{\sigma }}\right)}^{\mathrm{2}}},\end{array}$$

where σ is the standard deviation of the distribution defined by the kernel. We choose σ=2 years, following the method of Reichler et al. (2012) and as a method analogous to the 5-year smoothing applied to an SSW time series in DM21. The selection of σ=2 years allows contributions to the smoothed value from values approximately 7 years either side of the central year, as the value of the Gaussian window decays to near zero approximately 3.5σ from its mean. However, the largest contributions come from 3–4 years either side of the central year. This allows the smoothing to capture instances of ∼ 6–8 consecutive years with persistent vortex behaviour, which is a similar length to intervals observed in reanalysis (e.g. the 1990s; Pawson and Naujokat, 1999).

We subsequently define persistent NAM₁₀ intervals, when the vortex exhibits the same type of behaviour for a number of consecutive years, using extreme values of the smoothed NAM₁₀ index. A persistent NAM₁₀ interval is recorded when the smoothed NAM₁₀ index value falls within the top five percentile values. Once such an interval occurs, another cannot be recorded for 15 years afterwards to avoid choosing multiple central years within the same interval. Using five percentile values gives approximately the same rate of persistent vortex intervals as is reported in Reichler et al. (2012), so we proceed with this threshold throughout for a direct comparison with this study. Tests were also carried out to assess the sensitivity of our results to this threshold and are reported in Sect. 3.

We define an AMOC index, following the procedure in Reichler et al. (2012). The AMOC is defined using the overturning streamfunction field in the Atlantic sector. At each time point, the AMOC index is the maximum stream function value at any depth at a chosen latitude. We evaluate the index at 30, 45, and 50^∘ N and measure the response and co-variability with the NAM₁₀ time series and other climate indices defined below. We derive the observed AMOC index from the RAPID array data as the maximum meridional overturning circulation (MOC) at each time point at 26^∘ N. We also utilise a definition of the North Atlantic Oscillation from Hurrell et al. (2003). The NAO index is defined as the first PC of the Dec–Mar MSLP in the region 20–80^∘ N and 90^∘ W–40^∘ E. The PC is calculated by taking the first empirical orthogonal function (EOF) of deseasonalised MSLP anomalies and projecting this EOF onto the anomaly field. We additionally derive an ocean–atmosphere heat flux field defined as the sum of latent and sensible heat fluxes between the ocean surface and the atmosphere (i.e. positive values indicate the exchange of heat from the ocean to the atmosphere). We derive an index for the occurrence of deep convection anomalies in the equatorial eastern Pacific region. This index is defined by the top-of-atmosphere outgoing longwave radiation (OLR) averaged over the box 10^∘ S–10^∘ N and 240–290^∘ E. The OLR field is utilised as it acts as a proxy for the occurrence of convection anomalies. When deep convection is enhanced, the cloud top height is increased and, therefore, OLR is reduced. The eastern Pacific box is selected following a sensitivity analysis to establish the region which exhibits 90-year timescale variations in OLR. It is also a similar region to studies that consider east Pacific El Niño–Southern Oscillation (ENSO) patterns which are identified as a separate mode of variability to the traditionally used central pacific ENSO region (Johnson, 2013).

We also utilise the same QBO metric as in DM21. This is defined as the zonal mean zonal wind (ZMZW) averaged over the latitude range 5^∘ S–5^∘ N and the pressure level range 15–30 hPa. This metric captures the degree of vertical coherence in the QBO, an attribute shown to be important for QBO teleconnections with the NH midlatitudes (Andrews et al., 2019). Following the method of DM21, we also define the instantaneous amplitude of the deep QBO, which exhibits significant modulation in the westerly phase (see DM21 Figs. 10 and 11), using the Hilbert transform of the QBO (Hil[QBO(t)]), defined as follows:

$$\begin{array}{}\text{(10)}& \mathrm{Hil}\left[\mathrm{QBO}\right(t\left)\right]={\displaystyle \frac{\mathrm{1}}{\mathit{\pi }t}}\ast \mathrm{QBO}\left(t\right),\end{array}$$

where ∗ signifies a convolution, and t is discretised time. The time-varying amplitude, A(t), can then be expressed in terms of QBO(t) and Hil[QBO(t)] such that, in the following:

$$\begin{array}{}\text{(11)}& \mathrm{QBO}\left(t\right)+\mathrm{Hil}\left[\mathrm{QBO}\right(t\left)\right]i=A\left(t\right)e^{i\mathit{\theta }},\end{array}$$

where θ(t) is the instantaneous phase angle, which is a measure of signal progression through a cycle at time t.

Finally, we analyse the relationship between the magnitude of smoothed stratospheric NAM₁₀ extremes and 17-year lagged AMOC anomalies using linear regression with a single predictor (the lagged AMOC). We analyse the strength of the linear relationship using a correlation coefficient, r, and estimate a significance level for this value using a bootstrapping which assesses the probability that such a value results if the phases in signals in the NAM₁₀ and AMOC are randomly assigned but the overall autocorrelation structure is retained. We do this by comparing the r value calculated with real data with those produced from a set of synthetic NAM₁₀ series. These synthetic data are generated by applying a Fourier transform to the smoothed NAM₁₀ index, randomly shuffling the Fourier phases and subsequently inverse Fourier transforming to generate a surrogate time series with the same Fourier power spectrum as the real data. Repeating this data generation and calculating the correlation between the magnitude of positive extremes in the surrogate NAM₁₀s and the 17-year lagged AMOC builds up a probability distribution function (PDF) for the r value which can be used to estimate the significance level for a real r value.

3 Results

3.1 In-season surface responses to anomalous polar vortex events

We begin by diagnosing the in-season response to anomalous vortex events exhibited by surface variability in the model to assess its suitability for studying interactions on longer timescales. Figure 1 shows the mean sea level pressure (MSLP) composite differences between strong and weak polar vortex years (Fig. 1; top row). The composites have been determined by selecting MSLP values associated with events in which the daily NAM₁₀ values cross the +1.5 (strong) or −2 (weak) threshold (see Sect. 2.4). The composite differences demonstrate a significant lagged MSLP response, with strong (weak) vortex years corresponding to a positive (negative) NAO pattern, in agreement with previous modelling and observational studies (see Sect. 1). Additionally, a strong significant positive anomaly over the North Pacific (the Aleutian Low – AL) and corresponding negative anomaly over the Siberian region is evident in the month leading up to the vortex anomaly (−1–0-month lags). The former signal has been widely studied (Rao et al., 2019) and links the intensity of the AL to the strength of vertically propagating planetary waves that subsequently interact with the stratospheric vortex and influence its strength. DM21 examined this coupling using the same π control simulation as presented here and found a similar statistically significant relationship between the AL and the frequency of SSWs, but the regression coefficients were small in comparison with the QBO influence. Here the association between the AL and the vortex strength appears marginally stronger (r=0.39 with the NAM), which may be due to the NAM's ability to capture both types of vortex anomalies (strong and weak). In the same study, DM21 found that the AL exhibited minimal decadal to multi-decadal variability that was coherent with the decadal to multi-decadal variability in the vortex, and for this reason, the role of the AL is not considered in detail in this study.

Figure 1 Surface patterns associated with anomalous winter stratospheric NAM₁₀ events. The top row shows the monthly mean sea level pressure anomaly (hPa), the middle row shows the ocean–atmosphere heat flux, defined as the sum of latent and sensible heat fluxes (Wm⁻²), and the bottom row shows the SSTs (K). Coloured shading shows where the composite differences between strong and weak NAM₁₀ events are statistically significant at the 95 % level under a two-tailed Student's t test. The heading of each subfigure indicates the month range relative to the central date of each NAM₁₀ anomaly. Signals at negative times indicate that the surface anomaly leads the stratospheric NAM₁₀ anomaly. Signals at positive times indicate that the stratospheric NAM₁₀ anomaly leads the surface response.

The model also exhibits significant responses in ocean–atmosphere heat flux (Fig. 1; middle row). The largest flux anomalies are seen within 30 d (lag of 0–1 months), and their spatial pattern resembles that of a North Atlantic tripole with positive anomalies over the subpolar North Atlantic between approximately 50–65^∘ N, negative anomalies off the East Coast of the USA, and a second positive anomaly off the coast of northeast Africa. This pattern is consistent with the model response found by Reichler et al. (2012) to anomalous stratospheric NAM₁₀ events and the pattern associated with positive NAO phases in Delworth and Zeng (2016). As with the MSLP composites, there are visible anomalies in the 30 d leading up to the identified events (lag −1–0 months) both over the North Atlantic and Pacific regions. The Atlantic pattern may correspond to early responses to a disrupted or strengthened vortex and possible precursors to events. The Pacific anomalies preceding events are considerably smaller than the Atlantic anomalies and are concentrated over the Aleutian Low region.

The SST response to anomalous stratospheric NAM₁₀ events (Fig. 1; bottom row) over the North Atlantic lags behind the heat flux anomalies by around 2 months, with the largest amplitude anomalies at around 2–4-month lags. The anomaly pattern resembles that of the heat flux anomalies (with a change of sign), consistent with a mechanism in which the SSTs respond to the anomalous heat fluxes (as discussed in Hausmann et al., 2017). A prominent negative tropical east Pacific anomaly is obvious in the months leading up to anomalous vortex events, together with anomalies that resemble the Pacific Decadal Oscillation (PDO) in the region of the Aleutian Low (Mantua et al., 1997), and these features persist for several months. Variability in this region is dominated by El Niño–Southern Oscillation (ENSO) variations and a significant body of work (e.g. Domeisen et al., 2019) has proposed teleconnections between ENSO and vortex strength, consistent with the type of association exhibited here, i.e. negative (positive) SSTs or La Niña (El Niño) conditions associated with an anomalously strong (weak) vortex.

3.2 Surface impacts of persistent vortex anomalies

The in-season anomaly patterns associated with anomalous stratospheric NAM₁₀ events shown in Fig. 1 confirm that the model can reproduce the observed influence of vortex anomalies at the surface, particularly over the Atlantic region. We now extend the analysis to examine decadal-scale variability. Following the approach of Reichler et al. (2012), we smooth the NAM₁₀ index (Fig. 2) and then select the upper and lower five percentiles of this index to identify intervals with a persistent consecutive strong or weak polar vortex (see Sect. 2.4 for more details). The red and blue dots in Fig. 2 indicate the central year of intervals identified with persistent consecutive vortex anomalies (each dot represents the centre of intervals of approximately 8 years; see Sect. 2.4). Characteristic surface responses associated with these intervals are then analysed by compiling composites surrounding the central year of each positive and negative interval at lags of −40 (before the intervals) and 40 (after intervals) years. Calculating the (positive minus negative) composite difference can then be used to assess the potential surface impacts to observed intervals of persistent consecutive vortex anomalies, such as the consecutive strong anomalies throughout most of the 1990s and the consecutive weak anomalies in the early 2000s.

Figure 2 Time series of the December–March mean NAM₁₀ index level (green) and smoothed NAM₁₀ (purple) evaluated at the 10 hPa. The smoothed series is calculated by applying a Gaussian filter (σ=2 years) to the green series. Red and blue dots indicate the occurrence of persistent strong and weak vortex intervals, respectively, defined as extreme values (top and bottom five percentiles) in the filtered NAM₁₀ index. Interval central years are selected such that at least 10 years lie between consecutive intervals.

A lead–lag analysis of the composite differences in the AMOC strength at three different latitudes is shown in Fig. 3. The figure can be directly compared with Fig. 4c of Reichler et al. (2012), who suggest that decadal-scale variability in vortex strength acts to amplify a similar timescale of variability in the AMOC through resonance between the two signals. Similar to that work, an oscillatory AMOC response to the stratospheric anomalies is evident here, with significant positive anomalies in the AMOC at 45 and 50^∘ N at lags of approximately 3–5 years after persistent NAM₁₀ intervals, followed by negative anomalies at lags between 10 and 23 years. This response pattern is more clearly evident by taking the low-pass filtered versions of the AMOC responses (Fig. 3b, d and f). Even after the high-frequency signals have been removed, there are significant composite differences at ∼ 10–23-year lags, with the maximum responses of up to 1.5 Sv appearing between lags of 15–20 years. (We note, however, that this low-pass filtering of the AMOC time series reduces the overall variance so that the threshold for a composite difference to pass the significance test is lower, and this may increase the responses that are deemed significant, as shown in Fig. 3.)

Oscillatory response behaviour is not exhibited clearly by the AMOC at 30^∘ N, although extended negative anomalies at lags of ∼ 10–28 years after intervals are visible. Response patterns at this latitude are also significantly smaller than those at 45 and 50^∘ N in the filtered composites (Fig. 3a and b). One possible explanation of this is that the coupling mechanism between the NAM₁₀ and AMOC may involve an AMOC response that originates at higher latitudes and then propagates equatorward, which leads to less forcing of the AMOC evaluated further south. Zhang (2010) also note latitudinal differences in the AMOC response, so these differences between 30, 45, and 50^∘ N are not unexpected. The AMOC signals at 45 and 50^∘ N also exhibit significant positive anomalies preceding the persistent vortex intervals, at a lead of approximately 20 years, and this is also present, albeit smaller in magnitude, in the low-pass filtered indices. This precursor to persistent NAM₁₀ intervals is not found in corresponding results from Reichler et al. (2012), and the role of this feature is considered in more detail in Sect. 3.5.

Figure 3 Lagged response of the AMOC index to persistent NAM₁₀ intervals. Blue time series shows AMOC composite difference values between positive and negative NAM₁₀ intervals defined in Sect. 2.4. The x axis denotes the lead (negative values) or lag (positive values) relative to the interval's central year. Solid blue dots denote composite differences significant at the 95 % level under a two-tailed Student's t test. Panels (a), (c), and (e) show monthly AMOC composites, while panels (b), (d), and (f) show smoothed AMOC composites, using a Gaussian filter (σ=2 years).

We now examine this vortex–AMOC teleconnection in closer detail to explore possible physical pathways responsible for an AMOC response to persistent NAM₁₀ intervals. Figure 4b shows that there is also an oscillatory response in the NAO. This consists of a positive zero-lag response difference (consistent with Fig. 1), which is a significant negative NAO anomaly between lags of 10–18 years, followed by a positive NAO response at lags of around 28 years but with a smaller amplitude than the zero-lag response.

Also evident from Fig. 4b is that the oscillatory responses of the NAO and AMOC are similar, with the NAO leading the AMOC response by 2–3 years. Both responses vary with periods of 28–30 years, but, interestingly, the negative responses in both signals are larger and longer lasting at 10–20-year lags than those at 0 and 28–30-year lags. This difference in the magnitude and persistence suggests that the NAO and AMOC response patterns cannot be explained as a straightforward oscillatory response to the NAM₁₀ forcing at zero lag. If this were the case, then the response amplitude would be expected to decay with time, and the negative response at 10–20-year lag would be smaller than the initial positive response. Instead, a form of feedback mechanism is required to explain the amplified 10–20-year lagged responses or a resonant mechanism, as proposed in Reichler et al. (2012). If such a feedback mechanism were present, then one might also expect to see some oscillatory behaviour in the smoothed NAM₁₀ time series in response to the feedback from the surface. To investigate this, Fig. 4a shows the corresponding lead–lag difference analysis for the NAM₁₀ index (i.e. the smoothed NAM₁₀ composited around its own extreme values). This supports the presence of a feedback mechanism, since it also exhibits oscillatory behaviour with the same period of around 30 years. However, this oscillation is largely evident from the two positive peaks at 0- and 30-year lag, with the latter being substantially damped. There is no significant response at lags of 10–20 years, where the NAO and AMOC responses were largest. This suggests that the negative NAO and AMOC responses at these lags are unlikely to be due to resonance with the NAM₁₀ signal.

Alternatively, a possible physical pathway could involve an amplifying feedback mechanism between the NAO and AMOC responses. In this scenario, the positive zero-lag NAO response would drive a positive ocean–atmosphere heat flux anomaly over the subpolar North Atlantic, as seen in the response patterns to individual vortex events in Fig. 1. This heat flux anomaly would then lead to persistent negative anomalies in the near-surface ocean temperatures as heat is removed from the ocean via variations in wind stress and evaporation. The black line in Fig. 4b shows the lead–lag difference for ocean–atmosphere heat flux in the region encompassing the subpolar North Atlantic (45–65^∘ N, 15–60^∘ W). This region was selected to encompass the region with the largest heat flux response to individual vortex events in Fig. 1 (see middle row). Figure 4c shows the corresponding depth profile of the ocean temperature response from the same region. A positive heat flux response, and upper ocean (0–200 m depth) cooling, is evident at 0–1-year lags. This heat flux perturbation would, in turn, drive a positive AMOC anomaly at 2–3-year lags via changes in the mixed layer depth and deep convection in the subpolar North Atlantic, an effect discussed in Delworth et al. (1993) and Medhaug et al. (2012). This increase in AMOC strength would subsequently increase the Labrador Sea temperature via poleward transport of heat. This is confirmed by the positive, deep (down to 2000 m) ocean temperature anomaly at a lag of 10–20 years in Fig. 4c. In turn, the reversal of the Labrador Sea temperatures can feed back onto the NAO (see, e.g., Frankignoul et al., 2013), inducing a negative NAO phase at 10-year lags as the increased Labrador Sea heat content alters the ocean–atmosphere heat fluxes in the same region. Finally, this switch in the NAO phase would lead to a subsequent negative AMOC anomaly via the same heat flux mechanism outlined above for an opposite NAO phase. This sequence of feedbacks would thus act to enhance the persistence and magnitude of the secondary extreme in the NAO and AMOC. Reichler et al. (2012) also briefly suggest a similar mechanism to account for the AMOC response in their simulations but by involving a negative feedback of the AMOC onto itself and a role for the NAO.

Figure 4 (a) Composite differences in the smoothed NAM₁₀ index around extreme NAM₁₀ intervals (positive minus negative intervals). (b) As in Fig. 3f but for the AMOC at 50^∘ N (blue), the Dec–Mar mean NAO index (green) and the ocean–atmosphere heat flux (sum of latent and sensible fluxes; black) are averaged over an Atlantic box defined by 45–65^∘ N, 15–60^∘ W. All indices are smoothed with a Gaussian filter (σ=2 years). (c) Lagged responses of the ocean temperature anomaly depth profiles to persistent NAM₁₀ intervals. Composite differences between strong and weak intervals are shown for the same Atlantic box as the heat flux index. Hatching indicates composite differences that are not significant at the 95 % level under a two-tailed Student's t test.

3.3 Non-linear responses to strong and weak vortex intervals

So far, we have considered composite difference responses to persistent strong and weak vortex intervals. However, we know that the vortex evolution during strong and weak vortex years is very different, and this is likely to lead to differing interactions with surface and ocean variability. The surface responses to the two extremes are therefore unlikely to be equal and opposite. For example, weak vortex winters are mostly associated with SSWs whose impact at the surface is observed on average 0–60 d after their central date (Baldwin and Dunkerton, 2001). Furthermore, the vortex often exhibits a preconditioned state (see, e.g., Charlton and Polvani, 2007, and Bancalá et al., 2012) in which it becomes anomalously strong in the weeks running up to an SSW. So, the timing of SSW events within a given season will dictate both the overall strength of the NAM₁₀ measured over the winter season (which we use to construct the persistent NAM₁₀ index) and the overall strength of the subsequent surface response. In contrast, winters that exhibit an anomalously strong vortex will, on average, exhibit such behaviour throughout the whole season so that the impact on the surface will be present for a larger fraction of the winter season.

To assess the influence from each type of vortex extreme separately, Fig. 5 shows the lead–lag composite analysis of the NAO, AMOC, NAM₁₀, and heat flux signals for the persistent positive (strong vortex) and persistent negative (weak vortex) NAM₁₀ intervals separately. The AMOC patterns associated with each NAM₁₀ type are slightly different. The persistently strong vortex composite shows clear oscillatory behaviour with a period of approximately 28 years. These patterns resemble many of the features observed in the AMOC composite differences (Fig. 4d). Specifically, a positive AMOC anomaly is present at lags of 2–3 years, a negative AMOC anomaly at 10–23 years (with maximum amplitude at ∼ 15–20 years), and a second positive anomaly at approximately 30 years. On the other hand, the persistently weak vortex composites exhibit a more complicated response, with double-peaked minima at lags of −20 and −11 years and double-peaked maxima at 14 and 25 years. The weak vortex composites also exhibit no significant AMOC response at 2–3-year lags, unlike the strong vortex composites. Both event types are associated with NAO anomalies at zero lag, but the response to strong vortex intervals is larger in magnitude, which is consistent with the larger 2–3-year lagged AMOC response to this event type. The zero-lag NAO response is followed by an extreme of the opposite sign at approximately 16- and 14-year lags for strong and weak intervals, respectively. As with the AMOC composites, the NAO response to persistently strong vortex intervals exhibits a pronounced oscillatory behaviour of periods around 28 years. The corresponding NAM₁₀ analysis also shows oscillatory behaviour with periods of around 28 years. (We note that the NAM₁₀ results show statistical significance at both lead and lag times. However, the lead/lag interpretation is less meaningful in this case since the NAM₁₀ is used both as the signal and in the selection of the composites. The significance at both lead and lag times simply confirms that there is oscillatory behaviour.) The double-peaked behaviour of the AMOC associated with weak intervals is also reflected somewhat in the subpolar North Atlantic heat flux response (Fig. 5c), with positive response peaks at approximately 11- and 20-year lags. There is also a zero-lag heat flux anomaly associated with weak vortex intervals corresponding to the negative NAO response but the corresponding response to strong intervals is not significant.

To address the lack of heat flux signal, we also examine the responses of a more sophisticated metric of ocean–atmosphere heat flux (Fig. 5e) which captures the time variation of a more complex spatial structure of heat flux anomalies associated with vortex variations than a simple box average, i.e. the familiar tripole structure evident in the responses to individual vortex events shown in Fig. 1 with centres over the subpolar North Atlantic, off the East Coast of the USA, and the coast of northeastern Africa. This index is calculated by projecting the wintertime (Dec–Mar) annual mean flux field onto the pattern defined as the average 0–60 d composite response to individual vortex events (i.e. the mean of the patterns shown in Fig. 1, middle row ,at a lag of 0–1 and 1–2 months). This pattern is provided in Fig. 6. This index's responses recover some of the features of the AMOC composites, such as the ∼ 0 lag response to strong events and the ∼ 28-year period variation in the strong composite response. Double-peaked behaviour in the weak responses is also present in this flux metric, but, as with the box averaged metric, it is not as clear as the feature in the AMOC responses. More study is required into the precise patterns of heat flux anomalies associated with different vortex interval types to account for the patterns evident from these composites.

The asymmetry between the AMOC and NAO responses to persistently strong and persistently weak vortex intervals and the complexity of the separate responses show that the interactions between the NAM₁₀ and these surface modes are complex, with some suggestion of oscillatory behaviour on different timescales. In the following sections, we address these complexities in more detail by analysing the frequency spectra of the time series and show that some of these complexities can be explained in terms of the non-stationarity of the signals.

Figure 5 (a) Composites of (a) Gaussian smoothed NAM₁₀, (b) NAO, (c) subpolar NA heat flux, (d) AMOC anomalies at 50^∘ N, and (e) a NA heat flux index derived by projecting the Dec–Mar heat flux onto the response pattern in Fig. 6 that are associated with persistent vortex intervals of different types. On each subfigure, the red (blue) plots show the lead/lag responses to composites of strong (weak) persistent NAM₁₀ intervals. Solid dots denote composite anomalies are significant to the 95 % level under a two-tailed Student's t test.

Figure 6 Surface ocean–atmosphere heat flux pattern associated with anomalous winter stratospheric NAM₁₀ events averaged between 0–60 d, following vortex events. Coloured shading shows the composite differences between strong and weak NAM₁₀ events.

3.4 Non-stationary variability

In an analysis of this same UKESM π control simulation, DM21 showed that variability in the stratospheric polar vortex occurs on a range of timescales and is highly non-stationary. Although the composite analysis presented above shows oscillatory behaviour with periods of approximately 30 years, the results are likely complicated by the presence of non-stationary variability at other periodicities. We therefore analyse the frequency characteristics of the filtered NAM₁₀ index. Figure 7 shows the wavelet power spectrum of this index. It reveals variability on a range of timescales. As expected from the composite analyses, the spectrum exhibits intermittent power throughout the whole simulation at periods between 15 and 40 years. There is also significant power corresponding to a period of approximately 90–100 years that persists for ∼ 300 years of the simulation (year numbers 520–820; approximately three cycles) and power at the 50-year timescale that persists for 120 years (year numbers 500–620; approximately two cycles).

We note that much of the ∼ 30-year periodicity seen in the strong composite analysis is likely associated with significant wavelet power in the interval between years 300 and 410 because this interval displays the largest number of positive extremes in the smoothed NAM₁₀ time series (compare the 30 year wavelet power in Fig. 6 with the red dots in Fig. 2). The large number of elevated NAM₁₀ extremes selected in this interval is also likely affected by the presence of extremely long timescale variability. Qualitative inspection of Fig. 2 shows that the underlying NAM₁₀ amplitude increases from around year 200, reaches a peak at ∼ year 380, and thereafter declines. This means that years within this interval are more likely to reach the top five percentiles and qualify as an anomalously strong vortex. The origin of this multi-centennial variability is unclear, and robust analysis of such a low-frequency signal is difficult as the wavelet power of this multi-century variability is located mostly outside the so-called cone of influence that marks the boundary at which edge effects become significant.

Figure 7 (a) Dec–Mar NAM₁₀ values smoothed with a Gaussian filter (σ=2 years). (c) Wavelet power spectrum of time series in panel (b). Hatching represents areas outside the cone of influence in which edge effects are significant and power should not be considered. Blue contours represent the 95 % confidence level, assuming mean background AR1 red noise. (b) Morlet wavelet used for the wavelet transform in the time domain. (d) Global power spectrum, the wavelet power averaged over the whole simulation (blue line), and global 95 % confidence spectrum (red dashed line).

We also analyse the frequency characteristics of AMOC variability. Figure 8 shows the corresponding wavelet spectra of the AMOC at 50^∘ N and also the cross-wavelet spectra between the filtered NAM₁₀ and the AMOC, which gives a time-varying measure of co-variability between the two indices at different periods. The wavelet spectrum for the AMOC at 50^∘ N exhibits a peak in spectral power, corresponding to approximately 130 years that persists for nearly 400 years of the simulation (and also at longer periods up to 250 years, but boundary effects are an issue at these multi-centennial timescales, as discussed above). There are also portions of significant power at approximately 30- and 50-year periods, both of which persist for approximately two cycles (∼ 60 and ∼ 100 years, respectively) which are also apparent in the global power spectrum. We also note that the main 30-year power comes from the 300–400-year interval, coinciding with the interval of most activity in both the NAM₁₀ analysis (Fig. 7) and the selected high NAM₁₀ percentiles (Fig. 2).

The cross-power spectrum between the filtered NAM₁₀ and the AMOC shows three distinct features corresponding broadly to the three timescales prominent in the individual spectra of both indices. Significant cross-power is evident at 90–100-year periods for approximately 350 years (between 450 and 800 years; three–four cycles). The phase relationship between the signals (indicated by arrows in Fig. 8) within this portion of the cross-spectrum show a mixture of left-pointing arrows in the earlier portion, which indicate an anti-correlated relationship (π out of phase) and downward-pointing arrows in the later portion that indicate a $\frac{\mathit{\pi }}{\mathrm{2}}$-phase relationship. This later phase relationship can be interpreted in a number ways, with maxima in the AMOC leading to maxima in the NAM₁₀, minima in the NAM₁₀ leading to maxima in the AMOC, or maxima in the NAM₁₀ index coinciding with maxima in the rate of change of the AMOC (see Sect. 2.2 for more details of the cross-spectra arrows and how they are derived). There is also significant cross-spectral power centred around 30 years (between 300–400 years; three cycles) and 50 years (between 500–600 years; two cycles). In contrast to the 90-year periodicity, the phase arrows point to the right and slightly upwards, indicating that maxima in the NAM₁₀ index lead maxima in the AMOC by a small fraction of the cycle. This phase relationship is consistent with the composite analysis presented in Figs. 3 and 5 that indicated that a positive (negative) NAM₁₀ leads to a positive (negative) AMOC response approximately 2–3 years later. The wavelet spectra for the AMOC at 30 and 45^∘ N are provided in Figs. A1a and A2a. They show broadly the same features as that of the AMOC at 50^∘ N; however, it is notable that the AMOC at 30^∘ N does not exhibit significant variability on the 50- and 30-year timescales. This is reflected in the cross-spectrum with the smoothed NAM₁₀ index which shows minimal cross-power on these timescales (Fig. A1b).

To understand these non-stationary signals in the context of the proposed mechanism for vortex-AMOC interactions involving the NAO, we also analyse the power spectrum for the Dec–Mar NAO (Fig. 9). The wavelet power spectrum for the NAO (Fig. 9a) exhibits a portion of significant power at periods of 30 years, between ∼ 300–400 years, and a feature at 50 years, between ∼ 500–600 years, similar to the NAM₁₀ and the AMOC. Furthermore, the cross-power spectrum between the NAO and the NAM₁₀ indicate that signals in the two indices on the ∼ 30- and ∼ 50-year periods are also coincident in time for the 70–100 years they persist for. The phase relationship between these signals is small (arrows pointing to the right), indicating an in-phase relationship between the indices, which is consistent with the zero-lag relationship between the NAO and filtered NAM₁₀ extremes presented in the composite analysis (Figs. 4 and 5). The NAO wavelet analysis shows no significant power on the 90–100-year timescale. This suggests that the co-variability between the NAM₁₀ and AMOC on these longer timescales does not involve the NAO and is likely to arise through different mechanisms. We return to examine this feature in more detail in Sect. 3.5.

Results from this wavelet analysis may also explain the different behaviour of the AMOC around persistent positive and negative NAM₁₀ intervals (Fig. 5). First, we note that the contribution to the composites from the prime interval exhibiting ∼ 30-year oscillatory behaviour (∼ 300–400 years) comes solely from a collection of persistent strong NAM₁₀ intervals (see the red dots in Fig. 2). In contrast, a high proportion of the weak NAM₁₀ contributions to the composite analysis (8 out of 13) occur within the interval between ∼ 500–600 years, which exhibits variability at both 50- and 90-year periodicity. The complicated double-peaked behaviour of the AMOC response following persistent weak vortex intervals can now be better understood. The double minima in AMOC response at 10- and 20-year leads and at 15- and 25-year lags can now be explained as manifestations of the 50- and 90-year AMOC responses, e.g. a half-cycle between the minimum at 10-year lead and maximum at 15-year lag (a half-cycle of 25 years) gives a periodicity of 50 years, while a half-cycle between the minimum at 20-year lead and 25-year lag (thus a half-cycle of 45 years) gives a periodicity of 90 years. Additional support for this interpretation comes from the fact that the NAO composite analysis (Fig. 5b) shows a response that corresponds broadly with the AMOC minimum at a 10-year lead and maximum at a 15-year lag, suggesting a mechanism that involves the NAO on the 50-year timescale, but there is no corresponding response in the NAO at the 90-year periodicity, which is in agreement with the wavelet spectra and cross-spectra in Fig. 9, which indicates that the NAO does not exhibit variability on timescales of ∼ 90 years.

Figure 8 (a) AMOC time series at 50^∘ N (top), wavelet power spectrum (shaded contours represent wavelet power and dark blue contours the 95 % significance level compared to an AR1 process; bottom left), and global wavelet power spectrum (blue) and 95 % confidence level (dashed red; bottom right). (b) Cross-spectra between filtered Dec–Mar NAM₁₀ series and the AMOC index, NAM₁₀, and AMOC time series (top) and the cross-power spectrum (bottom). Shading indicates cross-power, dark blue contours the 95 % confidence interval, and arrows the relative phase angle between signals in the time series (to the right – in phase; vertically upwards – $\frac{\mathit{\pi }}{\mathrm{2}}$ out of phase, with positive peaks in the NAM₁₀ leading those in the AMOC; to the left – π out of phase; vertically downwards – $\frac{\mathit{\pi }}{\mathrm{2}}$ out of phase, with positive peaks in the AMOC leading those in the NAM).

Figure 9 Same as Fig. 8 but for the Dec–Mar NAO index. Panel (a) shows the wavelet power spectrum of the NAO, and panel (b) shows the cross-power spectrum between the NAO and the NAM₁₀ index.

3.5 Surface forcing of the stratosphere

The absence of an NAO signal that corresponds to the 90-year periodicities seen in both the NAM₁₀ and AMOC wavelet spectra and in the AMOC–NAM cross-spectrum (Fig. 8) indicates that the mechanism for the stratosphere–AMOC teleconnection on this longer timescale may be distinct in nature to those observed at the 30- and 50-year periodicities. As noted earlier, the phase relationship associated with this feature is also different to the shorter timescales ($\frac{\mathit{\pi }}{\mathrm{2}}$ out of phase), suggesting that the direction of causality is also switched, i.e. that the AMOC leads the stratosphere response on these timescales rather than vice versa. To study this more closely, we investigate possible pathways involving variability on this timescale.

In an earlier analysis of this UKESM π control simulation, DM21 highlighted the ∼ 90-year variability in the frequency of SSWs and demonstrated that this was closely associated with similar timescale variations in the amplitude of the QBO (and, in particular, the westerly QBO phase). Figure 10a shows the wavelet spectrum of the same QBO index employed by DM21 after smoothing with the same Gaussian filter utilised throughout this work (see DM21 Fig. 12 and also Sect. 2 for a description of how the QBO index was derived). A portion of significant power at ∼ 90-year periods persists for approximately 200 years of the simulation between years 600–800, which coincides with the significant response in the same interval of the NAM₁₀ and AMOC spectra at 90-year periodicity. Cross-spectra of the QBO index with the smoothed NAM₁₀ index (Fig. 10b) also corroborates the findings of DM21, with coincident signals at the 90-year timescales and right-pointing arrows that indicate an in-phase relationship, so that an interval of persistently strong positive (westerly) QBO anomalies coincides with an interval of persistent positive (strong) polar vortex anomaly. The sign of this teleconnection is consistent with the well-known Holton–Tan teleconnection (Lu et al., 2008, 2014), but it is present on much longer timescales. DM21 showed that this long-term variability originates primarily from long-term variations in the strength of the westerly QBO phase.

While the long timescale–QBO vortex teleconnection was demonstrated by DM21 (and confirmed here), the cause of the long-term westerly QBO variability was not established. In a preliminary investigation, DM21 performed a wavelet analysis of equatorial SSTs in various regions, including the equatorial east Pacific, to explore whether the ∼ 90-year QBO variability could be explained by SST triggering of convective activity that generates the gravity and other equatorial waves that contribute to the QBO. However, no 90-year periodicity in equatorial SSTs was found. DM21 also performed a corresponding analysis of an index representing the strength of the Aleutian Low, to explore whether that could explain the 90-year signals in both the QBO and the polar vortex via a modulation of the strength of large-scale planetary wave forcing, but no periodicity at 90 years was found in the AL index. Nevertheless, the extremely long period of the QBO and vortex variations suggests a driving mechanism that is likely linked to the ocean because of the characteristic long oceanic timescales. We therefore extend the investigation of DM21 to pin down the mechanism that links the observed variability in 90-year timescales in the AMOC, the QBO, and the polar vortex.

To extend this investigation, we examine the variability in the east Pacific top-of-atmosphere outgoing longwave radiation (OLR) as a proxy for deep convection (instead of using the east Pacific SSTs, as in DM21). When deep convection is enhanced, cloud-top height is increased, and therefore, OLR is reduced. Figure 10c shows the wavelet analysis of the Sep–Nov OLR in the east Pacific. It exhibits 90-year periodicity, significant cross-power with the smoothed QBO amplitude, and the NAM₁₀ index (Fig. 10d, e). The signals in the OLR and both the QBO and NAM₁₀ are anti-correlated (left-pointing arrows indicating a π phase difference). This is consistent with reduced OLR (increased deep convection) leading to greater QBO amplitude through increased wave forcing. The corresponding cross-spectra between the AMOC and the OLR metric (Fig. 10f) also indicate a significant portion of cross-power in the interval ∼ 600–800 years co-located with the feature seen in the NAM₁₀ spectrum (see the dashed contours in Fig. 10 which indicate the region with significant power in the NAM₁₀ spectrum). The phase relationship, in this case, is mostly $\frac{\mathit{\pi }}{\mathrm{2}}$ (the majority of arrows pointing upwards), indicating that one of the quantities depends on the time rate of change of the other. This result is similar to the study of Timmermann et al. (2005), who found a sensitivity of the equatorial Pacific region to periodic forcing of the AMOC. Their study showed a dependence of the Pacific thermocline on the rate of change of the AMOC. (We note, however, that the ∼ 7.5 Sv AMOC perturbation imposed in their study was considerably larger than the AMOC variations in our simulation.) Similarly, a lagged, cross-basin connection between the NA overturning circulation and Pacific sea surface height (SSH) was proposed by Cessi et al. (2004), who interpreted it in terms of the propagation of oceanic Kelvin and Rossby waves with anomalies communicated between Atlantic and Pacific via the Indian Ocean, as well as through the Drake Passage. We therefore suggest this as a possible pathway for the influence of the AMOC on the polar vortex at 90-year timescales in this simulation via modulation of deep convection in the East Pacific that influences the amplitude of the QBO.

Figure 10 (a) Wavelet power spectrum for the Sep–Nov deep QBO amplitude (see Sect. 2.4). (b) Cross-power spectra between deep QBO amplitude (Sep–Nov) and smoothed NAM₁₀ (Dec–Mar). (c) Wavelet power spectrum for the Sep–Nov area-weighted average equatorial east Pacific OLR (see Sect. 2.4). (d–f) Cross-power spectra between combinations of east Pacific OLR, deep QBO amplitude (both Sep–Nov means), AMOC at 50^∘ N (annual mean; all months), and NAM₁₀ (Dec–Mar) indices, which have all been smoothed with the same Gaussian filter that is used throughout (σ=2 years). (g–h) Cross-spectra between a Pacific SST index derived using the loading pattern in Fig. 11 and the east Pacific OLR and AMOC at 50^∘ N. Indices involved are indicated by the subfigure titles. Shading indicates the cross-power, using the same colour scale as in Figs. 6–8. Solid contours indicate the 95 % confidence interval for the power spectrum, and dashed contours show the 95 % confidence interval for the NAM₁₀ spectrum, for ease of comparison. Arrows indicate the relative phase angle between the signals in the indices.

An ongoing issue with this proposed pathway is the absence of a 90-year signal in the spectrum of equatorial east Pacific SSTs. One would expect a corresponding 90-year signal in this quantity, as sea surface heating of the atmosphere can lead to modulation of deep convection (Tompkins, 2001). However, DM21 showed the absence of a 90-year signal when examining a set of indices based on a box-averaged SST time series. One possible explanation for this lack of SST signal is that the spatial pattern of variability associated with the AMOC and OLR variation is more complex than a simple box average. A correlation map of the AMOC at 50^∘ N and the 45-year lagged equatorial Pacific SSTs (Fig. 11) shows significant correlations across the central and eastern Pacific region with centres on both sides of the Equator. Using this correlation map as a loading pattern and projecting it back onto the SST field allows us to derive an SST index which captures the temporal variability in this spatial pattern. The cross-spectrum of this index against the AMOC (Fig. 10g) and the east Pacific OLR (Fig. 10h) shows co-variation on the 90–100-year timescale, which may indicate that the variability in the SST pattern shown in Fig. 11 is more important on these timescales than a simple box average (as used in DM21). The phase relationship between the SST metric and the OLR is unclear, however. A more detailed examination of the intermediate steps in this proposed physical pathway is required to confirm its mechanisms, but this is outside the scope of the study. For example, The SST pattern in Fig. 11 may be linked to variability in the Pacific branch of the Walker circulation (the system of lower (upper) tropospheric easterlies (westerlies), with rising motion over the western Pacific and downwelling over the eastern Pacific) through changes in deep convection driven by SST heating of the atmosphere (Tompkins, 2001).

Figure 11 Correlation map of the AMOC at 50^∘ N vs. annual equatorial Pacific SSTs lagged by 45 years. Coloured shading indicates correlations significant to the 95 % level. Both indices are smoothed with the same Gaussian filter (σ=2 years) as used throughout this analysis.

3.6 Contribution of the stratosphere to recent AMOC changes

Our analysis of the UKESM simulation has identified co-variability between modes of variability in stratospheric circulation and the AMOC. Intervals in which the winter stratospheric polar vortex is consistently strong are, on average, followed by an extended negative anomaly in AMOC strength which peaks in magnitude at a lag of approximately 15–20 years (Fig. 3f). Recent observations of the AMOC have shown a negative trend in circulation strength of approximately −2.7 Sv between 2004 and 2012 (Smeed et al., 2018) before a marginal recovery after 2012 (Smeed et al., 2019). Modelling studies have proposed a key role for anthropogenic forcing in the AMOC slowdown over the 20th century and into the future (Liu et al., 2017; Bakker et al., 2016; Liu et al., 2019). However, the drivers of observed AMOC trends in the 21st century are not well understood. The results shown in the previous sections suggest the possibility that stratospheric variability has also contributed to the observed AMOC changes in response to nearly a decade of strong vortex years in the 1990s followed by a sequence of years with a weak, disturbed vortex in the early 2000s.

Using the relationship seen in the model between the modelled NAM, NAO, and AMOC, we can now compare the observed NAM, AMOC, and NAO indices for the interval 1979–2020 from the ERA5 and RAPID array datasets and assess the potential contribution of stratospheric variability to the observed AMOC trend (Fig. 12). The NAM₁₀ index (red bars) is characterised by an interval of strong vortex winters between 1988 and 1997 in which all but two winters exhibited a positive NAM₁₀. This interval also contains no SSWs (Pawson and Naujokat, 1999). This is followed by a run of winters between 1998 and 2005 which exhibit anomalously weak NAM₁₀ values with SSWs almost every year (Manney et al., 2005). The filtered NAM₁₀ index (red dashed line) reflects the presence of these time intervals, with a peak in positive values centred around 1995 followed by a negative extreme centred around 2003. The smoothed NAO (green dashed line) from the same dataset reflects some of these variations in the NAM₁₀, with positive NAO extremes in the 1990s. The long-term envelope of the NAO does not remain positive for as long as the NAM, primarily due to the anomalous negative NAO in 1996 (Halpert and Bell, 1997). The presence of this anomalous negative NAO in 1996 and the absence of a clear NAO anomaly in the following year, despite the presence of strong positive anomalies in the stratospheric NAM₁₀, is indicative of the range of other factors that influence the NAO in addition to the vortex. The AMOC strength estimated from the RAPID array observations between 2005 and 2019 is also shown (blue curve) and shows a negative trend between 2005 and 2012, followed by a recovery from 2012 onwards (Smeed et al., 2018, 2019).

Figure 12 Time series of the Dec–Mar NAO index (green bars) and NAM₁₀ index (red bars) from the ERA5 dataset. Dashed lines correspond to indices shown by bars smoothed with a Gaussian filter (σ=2 years). Also included is the annual AMOC time series at 26^∘ N from the rapid array dataset (blue).

The interval of observed consecutive strong NAM winters in the 1990s is anomalous in the reanalysis period, although the available data record is too short to allow a robust assessment. The amplitude and longevity of the observed anomaly are also large when compared to the UKESM simulation. Only two intervals in the UKESM simulation exhibit at least as many consecutive winters with strong (high NAM₁₀) conditions. These two intervals occur in the 300–400-year interval (centred around years 349 and 376), as shown in Fig. 13a. They each exhibit a sequence of 10 consecutive Dec–Mar anomalously positive NAM₁₀ values. The second interval exhibits 14 strong or marginally weak consecutive winters, with an allowance for the small negative NAM₁₀ value at year number 379. The presence of these two intervals is reflected in the smoothed NAM₁₀ values (Fig. 13b), and they represent the two of the three largest values of the filtered NAM₁₀ index (0.66 and 0.67). The corresponding smoothed AMOC index during these two intervals (blue curve in Fig. 13b) shows a positive AMOC response at lags of 2–3 years, followed by a negative response at 17–20 years, which is in good agreement with Fig. 3f. The negative responses at 17–20 years, following these two intervals of strong NAM₁₀ years, are the first and third greatest in magnitude compared to all other responses to persistent strong intervals. This is confirmed by Fig. 13c which shows the lagged AMOC response following all of the identified intervals with persistent positive NAM₁₀ anomalies (the two identified around years 349 and 376 are shown in black).

To estimate the response amplitude of the AMOC to an interval of persistently strong vortex winters, Fig. 13d shows a scatterplot of the central NAM₁₀ index of these intervals against the AMOC anomaly at 50^∘ N lagged by 17 years. This reveals a strong linear relationship ($r=-\mathrm{0.908}$) between the size of the persistent vortex anomaly and the subsequent negative anomaly in the AMOC at 50^∘ N 17 years later. The PDF of surrogate correlations used to assess the statistical significance (following the method outlined in Sect. 2) is displayed in Fig. A3a along with the correlation generated by the observed NAM₁₀ data. It shows that the r value lies well outside the distribution of surrogate correlations, indicating the high level of significance of the linear relationship.

A linear regression analysis on these data yields an estimated relationship between the variables, which satisfies the following:

$$\begin{array}{}\text{(12)}& {\mathrm{AMOC}}_{+\mathrm{17}}^{\prime }=-\mathrm{6.54}\mathrm{NAM}\max +\mathrm{3.11},\end{array}$$

where AMOC${}_{+\mathrm{17}}^{\prime }$ is the 17-year lagged AMOC anomaly at 50^∘ N, and NAM max is the magnitude of the positive extreme in smoothed NAM₁₀ at the centre of each interval. We can then use this relationship to predict the AMOC response to the observed sequence of the strong vortex years in the 1990s. The maximum smoothed NAM₁₀ occurs in 1996, so, using this relationship, the maximum AMOC response associated with the stratosphere would be expected 17 years later (2013), with an amplitude of −0.89 Sv (Fig. 13d; purple point).

Figure 13 (a) Dec–Mar NAM₁₀ index (red bars) from the UKESM simulation between year numbers 300 and 400. (b) AMOC (blue) and NAM₁₀ (red) indices smoothed with the Gaussian filter (σ=2 years) between year numbers 300 and 400 of the UKESM simulation. Black vertical lines show the location of the largest smoothed NAM₁₀ (c) AMOC response to persistent strong NAM₁₀ intervals. Light blue lines denote lagged AMOC responses of the AMOC from the whole UKESM simulation, and black lines show AMOC responses to NAM₁₀ intervals marked in (b) by vertical black lines. (d) Scatterplot of filtered NAM₁₀ index values occurring at persistent strong NAM₁₀ intervals throughout the whole UKESM simulation (y axis) against AMOC anomalies at 50^∘ N lagged 17 years after the persistent intervals' central year (x axis). Blue points indicate persistent intervals, and black dots represent the two intervals displayed in panel (b). The dotted line represents the linear line of the best fit for the points in black and blue. Also included is the 17-year lag AMOC anomaly predicted by projecting the regression coefficients used to construct the linear fit onto the maximum smoothed NAM₁₀ index in the ERA5 dataset (purple point).

This prediction suggests that approximately 30 % of the observed reduction in AMOC strength between 2005 and 2013 (0.89 Sv compared with 2.9 Sv in total) could be due to the response of the ocean to persistent forcing from consecutive strong vortex winters that occurred during the 1990s. We note, however, that our derivation of the modelled vortex–AMOC relationship is based on the AMOC response at 50^∘ N, where the response amplitude is largest, whereas the RAPID array dataset provides a measurement of AMOC strength at 26^∘ N. Figure 14 shows the scatterplot of filtered NAM₁₀ extreme magnitudes and lagged AMOC responses from the model at 30^∘ N. At this latitude, the linear relationship is significantly weaker than at 50^∘ N (r=0.652 for 50^∘ N vs. r=0.908 for 50^∘ N), but the correlation coefficient remains significant at the 95 % level. The predicted contribution from the strong vortex interval in the 1990s to the AMOC strength at 30^∘ N is reduced to −0.49 Sv (Fig. 14; purple dot), suggesting that approximately 17 % of the negative trend in the RAPID AMOC data may be due to stratospheric forcing from the 1990s. This is consistent with the composite analysis in Fig. 3a and b, which indicate that the modulation of the AMOC by the smoothed NAM₁₀ is less pronounced at 30^∘ N than at higher latitudes.

A similar analysis of the relationship between the magnitude of smoothed negative (weak) NAM₁₀ extremes and the lagged AMOC response (not shown) yields a much weaker relationship ($r=-\mathrm{0.21}$) that is not statistically significant. This asymmetry in the vortex–AMOC relationship between extreme positive and negative NAMs is perhaps not surprising, given that the surface impact of SSWs (that give rise to the negative NAM events) depends on the timing of the SSW within the winter season, whereas strong positive NAM intervals exhibit strong vortex conditions throughout the winter.

Figure 14 Same as Fig. 13d but for AMOC responses at 30^∘ N. A scatterplot of filtered NAM₁₀ index values occurring at persistent strong NAM₁₀ intervals throughout the whole UKESM simulation (y axis) against AMOC anomalies at 30^∘ N lagged 17 years after the persistent intervals' central year (x axis). Blue points indicate persistent intervals, and black dots represent the two intervals displayed in panel (b). The dotted line represents the linear line of the best fit for the points in black and blue. Also included is the 17-year lag AMOC anomaly predicted by projecting the regression coefficients used to construct the linear fit onto the maximum smoothed NAM₁₀ index in the ERA5 dataset (purple point).

4 Summary and conclusion

In this study, we have analysed the influence of persistent polar vortex extremes on surface and ocean circulation in a 1000-year π control simulation of UKESM1. Persistent vortex anomalies are identified using a smoothed NAM₁₀ index which characterises intervals of approximately 6–8 years during which the NH winter vortex is anomalously strong (positive NAM₁₀ anomaly) or weak (negative NAM₁₀ anomaly). While the surface impacts of stratospheric extremes in individual winters has received much attention (Baldwin and Dunkerton, 2001; Domeisen, 2019; Charlton-Perez et al., 2018), the surface impacts of consecutive sets of persistently anomalous winters have been less well studied, and teleconnections between the stratospheric vortex and ocean variability are not well characterised or understood.

We examine the AMOC response to long-term variations in the stratospheric polar vortex using composite analysis of the AMOC strength, following persistent anomalous NAM₁₀ intervals. We find oscillatory responses in the NAO and the AMOC consistent with previous work (Reichler et al., 2012). A diagnosis of the model supports a mechanism in which a persistently strong vortex (positive NAM) perturbs a positive NAO anomaly, which subsequently induces a positive AMOC response at 2–3-year lags via an increase in the subpolar North Atlantic ocean–atmosphere heat flux. This, in turn, feeds back onto the NAO to drive a reversal in NAO phase (to negative), which leads to a subsequent negative AMOC anomaly which peaks in magnitude at 15–20-year lags. The integrated effect of long-term oscillatory signals in the NAM₁₀ and the associated NAO variations may thus act as a metronome for the AMOC, which is a natural mode of oscillation in ocean circulation that varies on similar timescales. Further diagnosis of the separate impacts from intervals of persistent positive NAM (strong vortex) and persistent negative NAM (weak vortex with repeated SSW occurrences) showed that persistent strong vortex intervals had a much larger impact on the AMOC, which is perhaps not surprising because the vortex anomaly is consistently present throughout the whole winter.

We additionally found prominent non-stationary variations across multiple timescales in the AMOC, NAO, and NAM₁₀. Wavelet analysis revealed extended intervals in which 30-, 50-, and 90-year periodicities were dominant, so the composite response patterns were complicated by the superposition of contributions from intervals exhibiting different timescale behaviour.

Interestingly, while all three indices (NAM₁₀, NAO, and AMOC) co-varied at the 30- and 50-year periodicities, only the NAM₁₀ and AMOC co-varied at the 90-year periodicity, and co-spectra analysis suggested that the AMOC leads the vortex signal. This suggests a feedback of the AMOC variability onto the vortex that does not involve the NAO. A recent study (DM21) using the same UKESM simulation found long-term (90-year) co-variability between the vortex and the QBO. This suggests the possibility that the 90-year AMOC–vortex relationship could act via an influence of the AMOC on equatorial wave forcing of the QBO, which then influences the vortex through the well-known Holton–Tan relationship. This was explored through wavelet and co-spectra analysis of variations in tropical east Pacific deep convection and QBO amplitude. These showed similar 90-year co-variability, indicating that this is a plausible mechanism for the source of long-term variability in the QBO and the vortex, but further analysis of individual steps in the process, such as amplitude modulation of the various tropical waves that give rise to the QBO, would be required to confirm this.

Finally, we have applied the model results that link a lagged AMOC response to the presence of persistent vortex anomalies to assess the possible contribution of the interval of persistently strong vortex in the 1990s to the recent observed changes in the AMOC. Our analysis suggests a maximum AMOC response at 50^∘ N, with a lag of approximately 17 years. The lagged vortex–AMOC relationship is statistically significant, using the AMOC response at both 50^∘ N ($r=-\mathrm{0.908}$) and also at 30^∘ N ($r=-\mathrm{0.652}$), with the latter being closer in latitude to the RAPID array observations. Using a regression technique, we estimate that −0.49 Sv of the observed RAPID array AMOC trend by 2012 can be associated with the interval of persistent strong vortex behaviour centred on 1995. This represents nearly 17 % of the total decrease in AMOC transport between 2005 and 2013 in the RAPID array data. The observed negative trend in the AMOC has been attributed to a range of factors, including the influence of anthropogenic forcings (Caesar et al., 2018, 2021), but, to our knowledge, the potential role of vortex variability has not previously been considered. The origin of the interval of persistently strong vortex in the 1990s is unknown, but it was most likely due to internal variability. There is currently no consensus amongst climate models on how the vortex will respond to anthropogenic climate change (Ayarzagüena et al., 2020), and the 1990s appear to have been an anomalous period with no clear long-term trend emerging. As a result, our findings may indicate a significant role for internally generated signals in the recent negative AMOC trend.

There are several caveats to these results. First, the results come from a single model, although there has also been limited analysis of vortex–AMOC interactions in CMIP5 models (Reichler et al., 2012). On the other hand, there is evidence that GCMs underrepresent the influence of SSW events on the midlatitude tropospheric jet, the NAO, and surface temperatures, which is part of the signal-to-noise problem identified by Scaife and Smith (2018). Nevertheless, the model showed a clear NAO signal for up to 3 months following vortex anomalies, indicating a reasonable representation of stratosphere–surface interaction. A previous study using the UKESM (DM21) also noted an underestimation of SSW frequency compared to the ERA-Interim dataset, which indicates a positive bias in the mean vortex strength and the NAM₁₀, leading to the possible overrepresentation of positive NAM₁₀ intervals. The variability in the AMOC may also be underrepresented in the simulation, as Roberts et al. (2014) suggests that the models' decadal variability is smaller than that in the RAPID dataset. Additional analysis using a suite of CMIP6 models would clearly be useful to assess the robustness of the results presented here.

Composite and wavelet analysis of time series data presented here is effective for identifying co-variability between the stratosphere and ocean. However, these techniques are less capable of demonstrating causality within modes of variability, and so additional targeted experiments would be useful to establish this. Nevertheless, we have suggested physical pathways for vortex–ocean interactions on multi-decadal timescales which rely on well-established teleconnections (e.g. the in-season vortex–NAO connection) and demonstrate a possible key role for these interactions in recent AMOC behaviour. Our results stress the role and importance of non-stationary signals for understanding long-term variability in the climate system. This complexity can be detrimental to analysis which relies on common stationary methods such as the Fourier analysis. As a result, improved understanding and diagnosis of non-stationary climate variations and their underlying mechanisms will be key to overcoming these difficulties.

Appendix A

Figure A1 Same as Fig. 8 but for the AMOC at 30^∘ N. Panel (a) shows the wavelet power spectrum of the NAO, and panel (b) shows the cross-power spectrum between the AMOC and the NAM₁₀ index.

Figure A2 Same as Fig. 8 but for the AMOC at 45^∘ N. Panel (a) shows the wavelet power spectrum of the NAO, and panel (b) shows the cross-power spectrum between the AMOC and the NAM₁₀ index.

Figure A3 Probability distribution functions (PDFs) for correlations between the magnitude of NAM₁₀ extreme positive (a) and negative (b) values from surrogate NAM₁₀ data and anomalies in the AMOC at 50^∘ N evaluated 17 years later. Each NAM₁₀ surrogate is generated by applying a Fourier transform to the smoothed NAM₁₀ index, randomly shuffling the Fourier phases, and inverse transforming. Each PDF is built with 10 000 surrogates, and the correlation between the NAM₁₀ extreme magnitude and 17-year lagged AMOC anomaly at 50^∘ N is shown by vertical black lines in both subfigures.