A stretched polar vortex increased mid-latitude climate variability during the Last Glacial Maximum

The Arctic stratospheric polar vortex (SPV) is a key driver of winter weather and has been found to modify winter climate variability and its predictability in Eurasia and North America on inter-annual and decadal timescales. However, to what extent this relationship plays a role in driving climate variability on glacial–interglacial timescales is still unknown. Here, by systematically analysing SPV changes in four sets of PMIP4 (Paleoclimate Modelling Intercomparison Project – Phase 4) simulations for the Last Glacial Maximum (LGM) and the pre-industrial (PI) period, we explore how the SPV changed during the glacial climate and how it influenced climate variability. Our results show that under LGM conditions, the SPV stretched towards the Laurentide ice sheet; this was accompanied by anomalous upward wave propagation and enhanced SPV variability, which increased the likelihood of cold-air outbreaks into the mid-latitudes. During the LGM, this stretched SPV pushed cold Arctic air further equatorwards, increasing winter climate variability over the more southern mid-latitudes. In particular, in winter, this strengthened cooling over the mid-latitudes extended beyond the coverage of the Laurentide ice sheet (unlike in summer). SPV-induced temperature variability also explains the inter-model spread as removing SPV variation from the model results reduces the inter-model spread by up to 5°C over mid-latitude Eurasia. These results highlight the critical role of the SPV in connecting the polar region and mid-latitudes on glacial–interglacial timescales. These connections are reminiscent of intra-seasonal stratosphere–troposphere coupling.

1 Introduction

The stratospheric polar vortex (SPV) is an area of high-speed, cyclonically rotating winds present in both polar regions. Influenced by atmospheric waves propagating upwards from the troposphere, the formation and decay of the SPV occur on a seasonal timescale (Baldwin et al., 2003; Kolstad et al., 2010; Cohen et al., 2014). In the Arctic, it has been found that the SPV forms in autumn, when Arctic temperatures cool rapidly (Kolstad et al., 2010). The increased temperature difference between the polar region and the tropics causes strong winds to develop, and the Coriolis effect causes the vortex to spin up (Baldwin et al., 2003; Baldwin and Thompson, 2009). Interacting with the troposphere, the SPV strengthens when the latitudinal temperature gradient increases and reaches its maximum in winter (Cohen et al., 2021; Plumb, 1985; Takaya and Nakamura, 2001). The stratospheric polar vortex breaks down again in spring as the polar region warms up and the latitudinal temperature gradient decreases. These changes in stratospheric SPV strength can feed back, affecting the weather and climate closer to the Earth's surface by influencing Arctic air intrusions into the mid-latitudes (Baldwin and Dunkerton, 2001; Baldwin et al., 2003; Cohen et al., 2014). When the stratospheric vortex over the Arctic is strong, there is a single vortex with a jet stream that is well constrained near the polar front, and the Arctic air is well contained. When this northern stratospheric vortex weakens, it either breaks into two or more smaller vortices or is displaced away from the North Pole (Cohen et al., 2021; Kretschmer et al., 2018; Cohen and Jones, 2011). The flow of Arctic air then becomes more disorganized, and masses of Arctic air can push equatorwards (Kretschmer et al., 2018; Cohen et al., 2014, 2021). These seasonal dynamic changes make the SPV a major driver of winter weather over the mid-latitudes (Kolstad et al., 2010).

Apart from its seasonal changes, the Arctic SPV is also characterized by considerable inter-seasonal and inter-annual variability (Manzini et al., 2012; Reichler et al., 2012). Influenced by variations in atmospheric waves from the troposphere below, the stratosphere organizes chaotic wave forcings and creates long-lived changes (ranging from a week to several months) in hemispheric-scale circulation (McIntyre and Palmer, 1983; Cohen et al., 2014). This can occasionally trigger the breaking of stratospheric waves, analogous to ocean waves breaking on a beach, causing stratospheric airflow to become more disorganized and allowing masses of cold Arctic air to push further south, bringing with them a sharp temperature drop (Baldwin et al., 2003; McIntyre and Palmer, 1983). These sporadic occurrences of weak Arctic SPV events have significant impacts on surface weather and climate variability on inter-seasonal and inter-annual timescales (Cai et al., 2024; Kolstad et al., 2010; Manzini et al., 2012). In the troposphere, the negative phase of the North Atlantic Oscillation (NAO), defined as the surface sea-level-pressure difference between the subtropical high and the subpolar low, is observed following the weakening and warming of the stratospheric polar vortex (Yang et al., 2016). In this sense, Arctic SPV variation has been considered an important driver of weather predictability and climate variability over Eurasia and North America on inter-annual and decadal timescales (Kim et al., 2014, 2022; Zhang et al., 2022). However, whether it also plays a role in climate variability on longer timescales, such as on the scale of the glacial–interglacial cycle with dynamic continental ice sheets, has not been explored yet.

The Last Glacial Maximum (LGM; ∼21 000 years ago) is the most recent global cold extreme and has been widely documented in various proxy records (Cleator et al., 2020b). The LGM world was very different from the present world, with ice sheets covering northern North America and Fennoscandia, as well as the Greenland and Antarctic ice sheets that are still present today (Clark and Mix, 2002; Peltier et al., 2015). These extensive ice sheets represented large changes in topography and modified the spatial pattern of surface temperatures (Harrison et al., 2015; Kageyama et al., 2021). This alteration is expected to have generated planetary waves that could propagate into the stratosphere and affect SPV variation through the active interaction between the troposphere and stratosphere (Baldwin et al., 2003; Baldwin and Thompson, 2009).

As for the climate of the LGM, many studies – from proxy record compilations to the Paleoclimate Modelling Intercomparison Project (PMIP) and data assimilation – have been carried out to investigate the spatial patterns of cold climate features and their driving mechanisms (Harrison et al., 2015; Kageyama et al., 2021; Cleator et al., 2020b; Tierney et al., 2020; Annan et al., 2022). The enhanced cooling at high latitudes and the increased land–sea temperature contrast have been identified as key features of the glacial climate during the LGM, resulting from polar amplification (Kageyama et al., 2021). Another common feature emerging from proxy data and PMIP model results is the enhanced winter cooling over the mid-latitudes. Proxy-based reconstructions show a temperature reduction that is 5–8°C greater in winter than in summer (Cleator et al., 2020b). Model results from PMIP4 (the fourth phase of the Paleoclimate Modelling Intercomparison Project) and data assimilation results also show enhanced cooling in winter relative to summer (Annan et al., 2022; Kageyama et al., 2021; Tierney et al., 2020). The detailed structure of this enhanced winter cooling at the mid-latitudes, however, exhibits large spatial variation and varies widely among models (Annan et al., 2022; Kageyama et al., 2021; Tierney et al., 2020). For instance, two groups of data assimilation results show different levels of LGM cooling, with zonal-profile differences of 3°C over the middle and high latitudes (Annan et al., 2022; Tierney et al., 2020). Therefore, the research questions to be answered in this paper are as follows: (1) how did the stratospheric polar vortex (SPV) vary during the LGM glacial climate? And, more importantly, (2) how did these SPV changes contribute to climate variability during the LGM?

2 Methods

2.1 PMIP4 LGM simulations

Given its representativeness of full glacial conditions, the LGM has been the focus of the Paleoclimate Modelling Intercomparison Project (PMIP) since its inception (Braconnot et al., 2012; Braconnot and Kageyama, 2015; Harrison et al., 2015; Kageyama et al., 2017). Key climate drivers considered in PMIP simulations include extensive continental ice sheets, lower atmospheric greenhouse gas concentrations, and changes in orbital parameters (Kageyama et al., 2017). Compared with its predecessors, the new PMIP4 experimental protocol also includes some changes, including newly added forcings. For instance, PMIP4 allows vegetation and atmospheric dust loadings to change accordingly. Therefore, we focus here on PMIP4 results as these represent the most up-to-date simulations of the LGM climate.

We searched the Earth System Grid Federation (ESGF) database for all available PMIP4 LGM simulations and found five suitable simulations in total (Table S1 in the Supplement). Given our interest in exploring the relationship between mid-latitude winter climate and polar vortices, geopotential height and air temperature up to the stratospheric level are the two most important variables. Extra variables, such as sea-ice extent and sea surface pressure, can help us pinpoint the reasons behind these changes. With these target model output variables, five models provide downloadable data from the ESGF database for the LGM and pre-industrial (PI) period. Notably, since CESM2-FV2 and CESM2-WACCM-FV2 are from the same model family, we include only CESM2-FV2 for two reasons: (1) to maintain comparability with other PMIP models, which do not include a chemical component, and (2) because a previous study found no significant climatic differences between the two versions with respect to the LGM (Zhu et al., 2022). In the end, we selected four simulations for our analysis: MIROC-ES2L (hereafter MIROC), AWI-ESM-1-1-LR (AWI-ESM), MPI-ESM1.2-LR (MPI-ESM), and CESM2-FV2 (CESM-FV).

These four selected simulations were run with the corresponding models, each belonging to a fully coupled Earth system model. For MIROC, the atmospheric module is CCSR AGCM, which has 128×64 grid points in latitude and longitude and includes 40 vertical layers, with the top layer reaching 3 hPa. The oceanic component is COCO4.9 with nested sea ice, utilizing tripolar coordinates with 360×256 grid points in latitude and longitude, as well as 63 vertical levels (Hajima et al., 2020). MPI-ESM represents the atmosphere with ECHAM6, which has 192×96 grid points in latitude and longitude and 47 vertical layers, with the top layer reaching 0.01 hPa. The marine module is MPIOM1.63, utilizing 256×220 grid points in latitude and longitude, as well as 40 vertical layers (Mauritsen et al., 2019). AWI-ESM comprises the atmospheric component ECHAM6 (the same as MPI-ESM), the ocean–sea-ice component FESOM, and the terrestrial carbon model JSBACH. The atmospheric and oceanic components have average resolutions of around 250 and 100 km, respectively (Shi et al., 2020). For CESM-FV, the atmospheric component is CAM6, which operates with 144×96 grid points in latitude and longitude and includes 32 vertical layers, with the top layer reaching 2.25 hPa. The oceanic component is POP2 with the sea-ice module CICE5.1; it has 320×384 grid points in latitude and longitude and 60 vertical layers (Danabasoglu et al., 2020).

All four simulations were performed with prescribed atmospheric greenhouse gases, following the protocol given in Kageyama et al. (2017). According to ice core records, the greenhouse gas (GHG) concentrations during the LGM were generally lower than present-day levels (Bereiter et al., 2015). The considered GHGs include CO₂, CH₄, and N₂O, which were prescribed at concentrations of 190 ppm, 375 ppb, and 200 ppb, respectively. As for the orbital parameters, all the simulations were run with a prescribed eccentricity of 0.018994, an inclination of 22.949°, and a perihelion of 114.42° (Kageyama et al., 2017). These orbital configurations cause a slight decrease in summer solar radiation at high latitudes in the Northern Hemisphere and an increase in winter insolation, with a total magnitude of 10 W m⁻² (Kageyama et al., 2017). For more detailed information, refer to Kageyama et al. (2017, 2021).

2.2 Polar-vortex analysis

To investigate the effect of the stratospheric polar vortex on mid-latitudinal climate, we constructed composites of the climate with weak and strong vortex activities and compared them with the average climate state. The composite procedure is based on the vortex strength index (VSI) defined by Kolstad et al. (2010), in which $\mathrm{VSI}=-\sum \left(\right(Z-\stackrel{\mathrm{‾}}{Z}\left)\cos \mathit{\phi }\right)/\sum \cos \mathit{\phi }$, where Z is the geopotential height, $\stackrel{\mathrm{‾}}{Z}$ is the climatological mean of the geopotential height, φ is the latitude, and the sum is performed over all grid points north of 65° N. The reason for the minus sign is that the vortex is weak when the pressure is high and vice versa (Kolstad et al., 2010). This VSI is a conventional quantity for measuring vortex variability and has been validated as a reliable indicator of vortex variation on both seasonal and inter-annual timescales (Zhang et al., 2022).

The yearly varying VSI time series data were calculated using winter season (December–February (DJF)) geopotential height data at 20 hPa since the Arctic stratospheric polar vortex shows strong seasonal variations. As indicated by the monthly changes in the VSI, shown in Fig. S1 in the Supplement, the VSI calculated from the PI simulations is much weaker during the summer, with a negative value that remains at a stable level of less than −1000 gpm (geopotential metres). By contrast, the VSI is strengthened during the winter, exhibiting a positive value, and shows larger inter-annual variability (Fig. S1). This seasonal variability in the PI simulation is similar to that observed in the results of the ERA5 reanalysis data, suggesting that models are able to capture these variations. Accordingly, we selected years with strong and weak SPVs based on the corresponding winter (DJF) VSIs, using 1 standard deviation (σ). The strong SPV years are represented by those with a VSI larger than σ, while the weak SPV years are indicated by those with a VSI below σ (shown as red and green dots in Fig. S2). It is worth noting that MIROC shows much less inter-annual variability than the other models. For the PI simulation, MIROC has a σ value of 70 gpm, which is only one-third that of the other models. Further comparison with ERA5 reanalysis data suggests that MIROC seems to underestimate its inter-annual variability (Fig. S2). Similar to the previous analysis, we composite all the weak and strong SPV years to denote the climate states of weak and strong vortices, respectively (Kolstad et al., 2010; Zhang et al., 2022).

Figure 1 (a) Geopotential height at 20 hPa in the PI simulations (upper row) and LGM simulations (lower row), illustrating the shape and strength of the polar vortex. The dashed contour lines denote intervals of 200 gpm, ranging from 24 600 to 25 200 gpm. The red and blue lines (lower row) correspond to 25 000 gpm for the PI period and LGM, respectively. The black lines in the upper row correspond to the ERA5 reanalysis data for the period 1940–2024. (b) The upper row shows the December Plumb wave activity flux (WAF) anomalies (expressed in terms of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}$ m² s⁻²) at 100 hPa in the PI simulations, with positive values indicating an upward wave energy flux that vertically propagates from Rossby waves and negative values representing downward-reflected planetary waves. The lower row presents LGM WAF anomalies for the PI period (LGM-PI). The solid and dashed contours represent positive and negative WAFs, respectively, at intervals of $\mathrm{0.1}\times {\mathrm{10}}^{-\mathrm{3}}$ m² s⁻², based on the PI simulations for reference.

3 Results and discussions

3.1 The Arctic polar vortex stretched over ice sheets during the LGM

As demonstrated by troughs of low geopotential height in the stratosphere (e.g. at 20 hPa), the winter polar vortex was situated over the Arctic in the PI climate. The smallest geopotential height corresponds to around 24 500 gpm, with inter-model variations ranging from 24 400 gpm in MIROC to 24 600 gpm in AWI-ESM. Triggered by the asymmetric troposphere–stratosphere wave flux between Asia and North America, the SPV is not perfect circular but instead is stretched towards North America (Fig. 1a). An upward wave of energy over Asia, which is reflected downwards over North America, drags the centre of the SPV slightly towards the Atlantic side (Jones and Cohen, 2011; Kretschmer et al., 2018). This overall pattern fits the ERA5 reanalysis data, as shown by the similar shape of the 25 000 gpm contour (black line in Fig. 1a). A slightly more extensive SPV in the PI climate, compared to the SPV of the more recent climate used for the reanalysis, is expected, given that the PI climate was a bit colder than the 1940–2024 climate used for ERA5. This is consistent with the results of three out of the four models (with CESM-FV being the exception) (Fig. 1).

Compared to the PI period, one significant feature of the SPV during the LGM is that it was stretched towards the American continent. For instance, the 25 000 gpm contour line was stretched towards North America by 4–8° of latitude. Accordingly, the centre of the SPV shifted towards America. This ice-sheet-related stretching of the SPV seems to have increased temporal VSI variability, as indicated by the large standard deviations in three of the LGM simulations compared to the corresponding PI simulations. The exception is CESM-FV, which exhibits greater temporal variability during the PI period (much larger than that observed and seen in other models) than during the LGM (Fig. S2). When compared among different models, the magnitudes of the SPV responses during the LGM are consistent with the standard deviations of the inter-annual variabilities (Fig. S2). For instance, the relatively small LGM response of the SPV in MIROC aligns with its smaller inter-annual variability, which has a standard deviation of 50 gpm.

It has been illustrated that polar stratospheric variability is largely dominated by vertically propagating Rossby waves of tropospheric origin (Charney and Drazin, 1961; Cohen et al., 2007). During the LGM, the presence of the Laurentide ice sheet in the troposphere, up to 2–3 km height, may have generated planetary waves in the troposphere by modifying the topography (McIntyre and Palmer, 1983; Polvani and Waugh, 2004). Therefore, the ice-sheet-related planetary wave changes may have affected the SPV, which, in turn, descended and influenced the surface climate.

The stretched SPV during the LGM climate seems compatible with other studies on polar-vortex variations relevant to climate states. For instance, previous studies on decadal SPV variability have shown that the Arctic polar vortex shifted towards the Eurasian continent and away from North America in response to Arctic warming and sea-ice loss, particularly over the Barents–Kara seas in recent decades (J. Zhang et al., 2016). Our study reveals that the SPV shrank during the deglaciation from the LGM to the present, suggesting that the SPV evolved in the same direction as it has over the past few decades, despite the presence of different mechanisms (Manzini et al., 2012). Therefore, in line with previous relevant studies, our results suggest that a warmer climate favours SPV shrinkage and a shift towards Eurasia.

Previous studies have demonstrated that vertical wave activity flux (WAF) is a key determinant of the strength of the troposphere–stratosphere interaction (Baldwin et al., 2003; Jones and Cohen, 2011; Polvani and Waugh, 2004). To investigate the wave energy interface between the troposphere and stratosphere, we calculated Plumb WAF anomalies at 100 hPa for December by following Plumb (1985). The Plumb WAF results of the PI simulation show an overall positive vertical WAF over the North Atlantic and Eurasia, indicating upward propagation from Rossby waves. The downward WAF over the Pacific and North America denotes downward-reflected planetary waves in this region (Fig. 1b). This simulated Eurasia–America asymmetric pattern in the PI climate state is roughly consistent with ERA5 reanalysis data, showing a similar location for the zero contours (Fig. S3). The main difference is that the ERA5 reanalysis shows a larger magnitude for both positive and negative WAFs compared to the PI simulations, and this underestimation is a common issue in coupled climate models. This simulated Eurasia–America asymmetric spatial pattern aligns with our current understanding of troposphere–stratosphere coupling (Cohen et al., 2007; Jones and Cohen, 2011; Kretschmer et al., 2018). Compared to the PI simulation, the vertical WAF in the LGM simulation is enhanced over the Eurasian continent in all model results, implying that more wave energy propagates upwards to the stratosphere. Meanwhile, the vertical WAF during the LGM is reduced over the Pacific Ocean, indicating that less wave energy propagates in this region. This wave energy, which propagates into the stratosphere, is supposed to lead to the anomalous warming and weakening of the SPV (Polvani and Waugh, 2004). This seems unusual, given the overall cooler background climate, as previous studies have suggested that a warming climate may weaken SPV strength by enhancing turbulent heat flux (e.g. Kim et al., 2014; Kug et al., 2015). If a warming climate weakens the SPV, we would expect the SPV to strengthen in a cooler climate. However, we found that the anomalous upward WAF response during the LGM likely occurred due to the existence of ice sheets on the continents, which altered the topography and surface properties, thereby potentially outweighing the effect of the background climate. Previous studies have found that alterations in topography and snow cover may induce anomalous WAF and SPV responses (Cohen et al., 2007; Allen and Zender, 2010; Cohen and Jones, 2011; White et al., 2018; Pan and Duan, 2023).

Figure 2 Seasonal LGM temperature anomalies for the PI period (LGM-PI) (a) and their seasonal differences (b), defined as the DJF values minus the JJA (July–August) values, both taken from panel (a). The first four columns of plots show results from the simulations, while the last column presents the data assimilation from Cleator et al. (2020) (Cleator20). The magenta contour lines mark the extent of specific ice sheets during the LGM.

3.2 The impact of the polar vortex on mid-latitude climates

3.2.1 A stretched SPV enhanced winter cooling during the LGM

Compared to the PI period, the LGM climate was significantly cooled by the ice sheets, lower GHG levels, elevated atmospheric dust, and related feedback processes between different components of the climate system. The simulated temperature was reduced by more than 15°C over the ice sheets. However, a closer look at the temperature changes in the simulations and the proxy-based reconstructions shows stronger LGM cooling in winter than in summer (Fig. 2). The LGM summer cooling was strictly constrained over the ice sheets, highlighting the controlling effect of ice sheets on the summer climate. The primary mechanisms by which the ice sheets cool the climate include elevated altitudes, an enhanced ice–albedo feedback, and modified atmospheric circulation (Renssen et al., 2009; Y. Zhang et al., 2016). By contrast, this enhanced cooling extends further into the mid-latitudes of the continents during winter. For instance, the −10°C isotherm line of the LGM winter temperature anomalies extends south to nearly 30° N in North America, far beyond the coverage of the Laurentide ice sheet (Fig. 2). This implies that extra processes play a role in cooling the mid-latitudes during winter, in addition to the direct ice sheet cooling effects. From a climate-forcing perspective, orbital-scale insolation could potentially induce seasonal and latitudinal changes reaching up to this magnitude. Nevertheless, an LGM orbital setting very similar to that of the present day (as discussed in Sect. 2.1) excludes this possibility.

Figure 3 Winter surface air temperature (SAT) differences (colour-filled areas; expressed in degrees Celsius) between the weak polar-vortex composites (low VSI) and strong polar-vortex composites (high VSI) for both the PI period and LGM, with the significance denoted by dots. The sea-level pressure (SLP; contours) for DJF, with an interval of 120 hPa, is also shown.

The enhanced LGM cooling over the mid-latitudes is probably linked to the stretched polar vortex, which is an important driver of mid-latitude climate. First, planetary waves generated by the presence of continental ice sheets could extend the southern boundary of regions where cool air may reach further south and induce cooling at the mid-latitudes Kolstad et al., 2010). Furthermore, during propagation, the elliptical shape of the SPV, stretched by the ice sheet, can create irregular waves that are less stable than those of the round-shaped SPV of the present (Zhang et al., 2022). Both of these factors may lead to enhanced winter cooling at the mid-latitudes during the LGM. The large SPV-related variations in winter cooling among these four models contribute to the inter-model spread (to be discussed in Sect. 3.3.2).

3.2.2 Linkage between polar-vortex variation and winter climate

The composite analysis (i.e. the composited winter temperature anomalies) reveals that the effect of a weak SPV on climate is indicated by a dipole pattern between the middle and high latitudes of the Eurasian and North American continents, respectively (Fig. S3). This effect can be further illustrated by the surface air temperature (SAT) difference between the weak and strong SPV composites (Fig. 3). The weak SPV causes positive temperature differences of ∼1°C in the south and negative differences of up to 2–3°C in the north, with a boundary at 40–45° N over Eurasia. In North America, a weak SPV causes warm conditions at high latitudes and cold conditions at lower latitudes, with a boundary ranging from 45° N in MPI-ESM to 30° N in MIROC.

These patterns can be explained by the atmospheric circulation, as indicated by sea-level-pressure changes (Fig. 3). A weak SPV can result in a weakening of the subpolar low by inducing positive anomalies in sea-level pressure and an increase in the subtropical high by prompting negative anomalies over the North Atlantic. This increased polar high facilitates the flow of colder air from the Arctic into Eurasia, while a decreased subpolar low can reduce heat carried by air from the south. Together, these two processes cause dipole responses on the Eurasian continent. For the North American continent, the upstream region of the Pacific shows opposite responses, and the temperature shows a contrasting dipole pattern (Fig. 3). The results of the winter temperature and surface circulation responses to SPV variation align with observational analyses of inter-seasonal and inter-annual variations. The observational analyses show a strong connection between the strength of the stratospheric polar vortex and the dominant pattern of surface weather variability, such as the North Atlantic Oscillation (NAO) (Yang et al., 2016). In this sense, the slowly varying stratospheric signal may help predict changes in the North Atlantic Oscillation and the weather over the coming months (Baldwin and Dunkerton, 2001).

During the LGM, the overall pattern of composite temperature anomalies for weak and strong VSIs was similar to that of the present day (Fig. 3). However, the strength of the temperature anomalies is generally enhanced over Eurasia in most of the models, with the exception of CESM-FV. For instance, the cooling anomalies increase from 2 to 3°C in MIROC. Other visible differences include the spatial extent of warming and cooling. For instance, the negative temperature anomalies over Eurasia extend slightly further south during the LGM compared to the PI period. Therefore, the difference in surface winter temperatures between strong and weak SPV conditions during the LGM generally reaches up to 2–4°C, which is slightly stronger than in the PI period.

Figure 4 Winter month climate variability, indicated by the standard deviations (SDs) of winter temperatures for the PI period and LGM (a); their LGM anomalies (b), calculated using panels (a) and (b); and the SD reductions (“PI_redct” and “LGM_redct”) (c). The SD reductions are each defined as the difference between the total standard deviation and the SD when the strong and weak SPV years (used for the composites) are removed. The magenta and green lines in panels (a) and (b) indicate the sea-ice margins present during winter for the simulations, where available (MIROC and MPI-ESM in this case).

3.3 SPV variation increases climate variability and uncertainty

3.3.1 Large inter-annual climate variability at the mid-latitudes

The spatial distribution of the temperature standard deviations shows large climate variations over mid-latitudinal continents and near the margins of sea-ice extensions (Fig. 4a). Compared to the PI period, climate variation during the LGM is generally enhanced over both land and ocean at the mid-latitudes. The region of large climate variability moves southwards during the LGM. This leads to dipole patterns in the temperature variability differences between the LGM and the PI period, with significant enhancements in the south and a reduction in the north (Fig. 4b). The enhanced climate variation over the North Atlantic during the LGM is in line with the southward-extended sea-ice margin (Fig. 4a and b). This implies that dynamic sea ice plays a controlling role in climate variability.

On land, the enhancement of LGM climate variability at the mid-latitudes of Eurasia and North America is visible and probably related to the stretched SPV (Fig. 4a and b). Removing the strong and weak SPV composites can clearly reduce climate variability over the mid-latitude continents (Fig. 4c). A common feature of this reduction across the models is its positive anomalies, despite some variations in the detailed patterns. Compared to the PI period, the region with the large reduction is located further south during the LGM, which is consistent with a stretched SPV. This enhanced climate variability during the winter over the mid-latitudes also appears in previous studies on the seasonal climate of the LGM (Cleator et al., 2020; Kageyama et al., 2021). Winter temperature variability has been found to be enhanced over the middle and high latitudes, further contributing to a large portion of LGM climate variability (Annan et al., 2022; Cleator et al., 2020; Kageyama et al., 2021).

Figure 5 Inter-model spread of LGM temperature anomalies for winter (DJF). Panel (a) shows the inter-model spread -of LGM temperature anomalies, defined as the root mean square deviation of LGM anomalies from the ensemble mean. Panel (b) shows the difference between the total inter-model spread of LGM temperature (a) and the inter-model spread determined after removing the weak and strong VSI composites. Panel (c) shows the latitudinal profile from panel (b) over the continents, with the latitudinal positions of ice sheets during the LGM overlaid on the profile. Dev. reduction: deviation reduction.

3.3.2 Large inter-model spread of winter temperatures

Although the four models generally exhibit similar patterns of LGM temperature anomalies, there are still visible differences. The root mean square deviation of individual models' temperatures shows that winter exhibits a larger inter-model spread than summer (Fig. 5). The largest inter-model spread is found over the Nordic seas, which is due to significant climate uncertainty induced by dynamic sea ice. Over the continents, the inter-model spread is generally larger than over the oceans, and large values are found over the mid-latitudes. This wide spread over the mid-latitudes is the primary contributor to climate uncertainty.

The large spread of multi-model temperatures over the mid-latitudinal continents can be significantly reduced by removing the strong and weak SPV composites (Fig. 5b). For zonal mean temperatures, the winter inter-model spread decreases by almost 0.8°C for the latitude range of 30–50° N when removing the weak and strong SPV years. In particular, the reduction over North America reaches up to 5 °C and is mainly distributed over 30–45° N (Fig. 5b). For the Eurasian continent, a similar reduction is found over the latitude range of 40–65° N.

Previous studies have found a less consistent winter climate across different models compared to the summer climate (Harrison et al., 2015; Kageyama et al., 2021). Our analyses identify a previously unknown source of mid-latitudinal climate variability and inter-model spread, suggesting that SPV variations need to be taken into account for realistic simulations of climate variability during the glacial–interglacial cycle.

4 Conclusions

Our analysis of changes in the stratospheric polar vortex (SPV) during the LGM suggests that it played a key role in linking the polar region and mid-latitudes, even in glacial climates, and in adapting to the glacial–interglacial cycle. SPV weakening caused two dipole structures in continental winter temperature during the LGM that were similar to those observed in the PI period, suggesting that the mechanisms by which the SPV affected the tropospheric climate operated during both periods. A comparison of the SPV between the PI period and the LGM shows that during the LGM, the SPV was stretched towards the American continent, which pushed the regions affected by outbreaks of cold air further south. The difference in surface winter temperatures between strong and weak SPV conditions during the LGM reached up to 2–4°C, which is slightly greater than that observed during the PI period. During the LGM, the stretched SPV pushed cold Arctic air further equatorwards, increasing mid-latitudinal winter climate variability. Removing the SPV variations can reduce inter-annual variability in winter temperature by up to 5°C over mid-latitude Eurasia. SPV-induced temperature variability also explains the inter-model spread as removing SPV variation persistently reduces winter temperature variation (i.e. the root mean square deviation) across the mid-latitudes by 0.8°C. These results highlight the critical role of the SPV in linking the polar region and mid-latitudes, even on glacial–interglacial timescales.