The Last Interglacial (LIG; 130–115 ka) is an important test bed for climate science as an instance of significantly warmer than preindustrial global temperatures. However, LIG climate patterns remain poorly resolved, especially for winter, affected by a suite of strong feedbacks such as changes in sea-ice cover in the high latitudes. We present a synthesis of winter temperature and precipitation proxy data from the Atlantic seaboard of Europe, spanning from southern Iberia to the Arctic. Our data reveal distinct, opposite latitudinal climate trends, including warming winters seen in the European Arctic while cooling and drying occurred in southwest Europe over the LIG. Climate model simulations for 130 and 120 ka suggest these contrasting climate patterns were affected by a shift toward an atmospheric circulation regime with an enhanced meridional pressure gradient and strengthened midlatitude westerlies, leading to a strong reduction in precipitation across southern Europe.

The Last Interglacial (LIG; ca. 130–115 ka; roughly coincident with marine isotope stage 5e) is an important case study for the response of the climate system to changes in climate forcing, with negligible human impact and only slight changes in total greenhouse gas forcing, but some of the largest anomalies in orbitally driven insolation forcing during the past million years (Berger and Loutre, 1991). This combination of factors led to the last known instance of warmer than preindustrial global temperatures in geological history, with average global annual sea-surface temperature estimated at +0.5 ± 0.3 °C above preindustrial (1870–1889 CE) levels (Hoffman et al., 2017) and summer surface warming reconstructed at up to 4–5 °C over Arctic land areas (CAPE–Last Interglacial Project Members, 2006), with the global surface temperature anomaly estimated at approximately +0.8 °C (Fischer et al., 2018).

While the LIG is thus an important target for testing the sensitivity and robustness of climate model simulations under a range of forcing regimes, the evolution of LIG winter climate has been especially difficult to unravel due to the model dependence on sea-ice feedbacks, which strongly affect high-latitude winter temperatures (Bakker et al., 2013). In Europe, ensembles of transient climate model simulations show poor agreement in winter temperature trends, and often opposite trends compared to available proxy data, while for summer, individual models and proxy data show a consistent early LIG temperature maximum in mid- and high latitudes (Brewer et al., 2008; Bakker et al., 2013; Salonen et al., 2018). The modeling challenges are compounded by the paucity of proxy data from the high northern latitudes due to erosion of LIG deposits during the last glaciation (Helmens, 2014).

To elucidate European winter climate evolution during the LIG, we present a data–modeling synthesis with high-resolution climate proxy archives spanning along the Atlantic seaboard from southern Iberia to the Arctic. Our data reveal strongly contrasting trends in the northern and southern European winter temperature and precipitation, consistent with a shift toward an atmospheric circulation regime with strengthened midlatitude westerlies seen in early and late LIG climate model experiments.

Proxy Data Sets

We reconstructed southwest European winter temperature and precipitation using pollen data from three marine cores drilled from R/V Marion Dufresne during the IMAGES cruises of 1995, 1999, and 2004 (Fig. 1; for details about the marine cores, see Sánchez Goñi et al. [2018]). Cores MD99–2331 from the Galician margin (42.15°N, 9.68°W) and MD04–2845 from the Bay of Biscay (45.35°N, 5.22°W) captured pollen rain from the temperate forest of the adjacent landmass, with winter temperature as a primary ecological control. We inferred winter temperature qualitatively as the sum of temperate tree pollen and quantitatively as pollen-based mean January temperature (TJan) reconstructions. Core MD95–2042 was drilled farther south at the southwest Iberian margin (37.80°N, 10.17°W), and in the Mediterranean climate zone, with pollen data primarily influenced by winter precipitation. To represent winter precipitation variations, we used the pollen sums of Mediterranean forest taxa and semiarid taxa. Several isotopic stratigraphic events were identified and dated by Shackleton et al. (2003) in core MD95–2042 and used to develop a linear-interpolation age model. The chronologies of cores MD99–2331 and MD04–2845 were based on correlating forest increases to equivalent pollen stratigraphic features identified and dated in core MD95–2042 (Sánchez Goñi et al., 2005, 2012).

For northern Europe, we used the Sokli sequence from northern Finland (67.80°N, 29.30°E), with pollen-based reconstructions of both January and July (TJul) mean temperatures (from Salonen et al., 2018) and supporting aquatic proxy data for winter conditions (Plikk et al., 2016, 2019; Kylander et al., 2018). In the Sokli age model, the interglacial onset and end identified in the pollen data are aligned with the start and end of interglacial conditions in northern and central European U/Th-dated speleothems (Salonen et al., 2018).

Paleoclimate Reconstruction

We reconstructed TJan from cores MD99–2331 and MD04–2845 using the same method as that used earlier for Sokli (Salonen et al., 2018). This method uses an ensemble of six pollen–climate calibration models, but with the calibration data set here expanded southward to include the Mediterranean climate zone and a total of 1295 surface pollen samples (see the Supplemental Material1). Each reconstruction was summarized as the median of the six models, and with a LOESS (locally estimated scatterplot smoothing) smoother (span = 0.1) fitted to the median. The performance of the calibration models was estimated by 10-fold cross-validations, suggesting coefficients of determination (R2) ranging from 0.72 to 0.91 and reconstruction errors (root-mean-square error of prediction) from 2.08 °C to 3.67 °C (for further details, see Figs. S1 and S2 and Table S1 in the Supplemental Material).

Paleoclimate Modeling

We discuss here results from two global climate model simulations for 130 ka and 120 ka, performed by Otto-Bliesner et al. (2013) using the fully coupled Community Climate System model, version 3 (CCSM3,, with dynamic components for the atmosphere–land surface–ocean–sea ice system (Collins et al., 2006). The main difference between these two simulations concerns changes in orbital parameters, leading to a marked seasonal difference in the incoming solar insolation at the top of the atmosphere at 50°N decreasing by 20 W/m2 from 130 to 120 ka in July, but increasing by 10 W/m2 in January (Berger, 1978). In addition, the atmospheric levels of greenhouse gases CO2 and CH4 declined, as CO2 decreased from 300 ppmv (parts per million volume) at 130 ka to 272 ppmv at 120 ka, and CH4 decreased from 720 ppbv at 130 ka to 570 ppbv at 120 ka. All other boundary conditions (including land-ocean distribution, vegetation, ice sheets, solar constant) were kept fixed at modern values. The applied model version has a T85 spectral resolution in the atmosphere, corresponding to ∼1.4° latitude-longitude grid spacing. The ocean model has a horizontal resolution of ∼1° latitude-longitude.

In northern Europe, the TJan reconstruction from Sokli shows a strong warming trend spanning the LIG (Fig. 2B). In the underlying pollen data, distinct plant indicators of oceanic climate like Corylus (hazel), Quercus (oak), and Osmunda (royal fern) (Birks and Paus, 1991; Salonen et al., 2012; Seppä et al., 2015) endure in the late LIG during the summer cooling trend (Fig. 2A). The mild late-LIG winters are supported by diatom data from the same sequence, indicating a reduced seasonality (Plikk et al., 2016). These multiproxy data from Sokli are further corroborated by the Mertuanoja sequence from western Finland, which, while of low resolution, shows a distinct late-LIG peak in Osmunda (Eriksson et al., 1999). Reliable precipitation or water balance proxy indicators are rare in northern Europe. At Sokli, chironomid, plant macrofossil, diatom, and geochemical data indicate increased chemical weathering and input of dissolved organic carbon during the late LIG, which may reflect increased precipitation, but which may also be related to an expansion of Picea (spruce) in the catchment (Plikk et al., 2016, 2019; Kylander et al., 2018). By contrast, in southern Europe, a winter cooling trend is indicated by the decreasing temperate forest pollen percentages and the falling reconstructed TJan in cores MD04–2845 (Figs. 2D and 2E) and MD99–2331 (Figs. 2F and 2G), while a winter drying trend is suggested by the decrease in Mediterranean forest pollen (Fig. 2H) and the increase in semiarid pollen (Fig. 2I) in core MD95–2042. In the intervening central European sector, the proxy data are more ambiguous, with the synthesis of Zagwijn (1996), which incorporated pollen and plant macrofossil data from 31 core sequences from 45°N to 56°N and used the reconstructed paleo-isotherms to extract a paleotemperature curve for one site (Amsterdam), showing a mid-LIG maximum for TJan (Fig. 2C). However, other studies from individual sites at this latitude band have indicated an early-LIG TJan maximum, followed by either moderate (Kühl and Litt, 2003) or strong (Field et al., 1994; Brewer et al., 2008) cooling.

Our data thus show strong and opposite trends in northern European (warming) and southern European (cooling and drying) winter climate over the LIG. Importantly, the reconstructed early-LIG (ca. 130–125 ka) timing of the southern European TJan maximum contradicts the modeling of Bakker et al. (2013), which suggested a TJan maximum at 120–119 ka in northern Spain and western France, but with large uncertainties (standard deviation ∼4 k.y.). No timing was established by Bakker et al. (2013) for our northern European data region, as most of the ensemble models were unable to find a TJan peak outside the errors of the LIG mean.

The latitudinal pattern of temperature and precipitation trends seen in our data is broadly consistent with the spatial signature of a shift in North Atlantic Oscillation (NAO)–type circulation, from a more negative (NAO–) state in the early LIG to a more positive (NAO+) state in the late LIG (Wanner et al., 2001). The simulated sea-level pressure (SLP) patterns for 130 ka (Fig. 3A) and 120 ka (Fig. 3B) show a more extensive Icelandic low, and also a more expressed Azores high at 120 ka, while the south-to-north SLP gradient becomes steeper at 120 ka (Fig. 3C). The NAO index values, calculated based on the difference of normalized SLP between Lisbon (Portugal) and Stykkisholmur/Reykjavik, show a shift from slightly negative (−0.245) at 130 ka to moderately positive (+0.829) at 120 ka. Zonal winds become stronger at 120 ka compared to 130 ka over most of northwest Europe (Fig. 3D), consistent with a steeper SLP gradient at 120 ka. For winter precipitation, the modeled trends (Fig. 3E) share some large-scale features of a NAO+ signature (Wanner et al., 2001; Gouveia et al., 2008), with a tendency toward increased precipitation over northern Europe, shifting to a drying trend at approximately 50°N and with reduced precipitation also modeled around southern Greenland. Our precipitation proxy data site, MD95–2042, is located within a broad latitudinal band marked by a drying trend, and here the modeled trend matches the proxy data (Figs. 2H and 2I). However, the LIG TJan patterns (Fig. 3F) deviate from a classical, mean NAO+ signature, which is characterized by a positive winter temperature anomaly over most of Europe north of 40°N (Wanner et al., 2001; Gouveia et al., 2008). Instead, a negative TJan trend is modeled over northwest Europe, similar to an east-west tilted NAO+ configuration, but lacking the cooling belt along the Mediterranean (Rousi et al., 2020). In consequence, both our northern (Sokli) and southern European (cores MD99–2331 and MD04–2845) temperature proxy sites are located in regions where the sign of the modeled trend changes. This is also true for the central European zone, which has a high density of LIG sites (Zagwijn, 1996) showing opposite temperature trends modeled across a southwest-northeast–trending line cutting across France and Germany. While some of the differences in published winter temperature reconstructions have been attributed to different reconstruction algorithms (Kühl and Litt, 2003), the discrepancies may thus also be underlain by real shifts in the temperature trend within central and western Europe.

The change to a regime with an enhanced meridional pressure gradient and strengthened midlatitude westerlies, as suggested by our data and modeling, was likely forced by the decrease in summer insolation (Fig. 2L), which resulted in expansion of the ice sheets, a more extensive sea-ice cover and snowpack, and a shift from taiga to tundra vegetation commencing in climate modeling at ca. 122 ka (Crucifix and Loutre, 2002) and supported by European proxy data sets showing a major concurrent southward displacement in the deciduous tree line (Sánchez Goñi et al., 2005). The insolation forcing and albedo feedbacks significantly amplified the north-south thermal gradient over ca. 122–120 ka (Sánchez Goñi et al., 2005), leading to a more vigorous atmospheric circulation, implying stronger winds after the LIG temperature maximum at ca. 130–125 ka (Rind, 2000). Our data are broadly consistent with the modeled 122 ka timing of the mid-LIG transition (Crucifix and Loutre, 2002), with the best breakpoints (Gurarie, 2014) in the interglacial sections of all seven winter proxy series timed at 122.2 ± 1.6 ka (Figs. 2B, 2D2I).

The strong winter warming trend seen in the northern European proxy data could have been influenced by the Atlantic Meridional Overturning Circulation (AMOC), which, based on transient modeling experiments (Bakker et al., 2013) and proxy data, continued as vigorous or even strengthened (Sánchez Goñi et al., 2012; Galaasen et al., 2014) through the late LIG (Fig. 2J), and proximally by the warm late-LIG Nordic seas (Fig. 2K) (Zhuravleva et al., 2017). However, the reverse trends of strong winter cooling and drying in southwest Europe run against these forcing mechanisms, as well as the rising winter insolation in the midlatitudes (Fig. 2M), and so they would need an additional explanatory factor. Also, the modeled changes in TJan are within ±0.5°C (Fig. 3F) compared to the reconstructed TJan decrease of up to 10°C over the LIG (Figs. 2D and 2F). The muted modeled temperature changes could be the result of missing feedbacks due to the fixed vegetation and the fixed ice sheets in Greenland and Antarctica. Moreover, the preindustrial control run of the CCSM3 model shows a cold bias in the North Atlantic region (Lunt et al., 2013), possibly leading to a damped cooling response due to restricted capacity for warming at 130 ka. By comparison, the modeling and data agree on the southwest European precipitation trend, with a strong reduction in modeled January precipitation (Fig. 3E) likely sufficient to explain the response of the southern Iberian vegetation (Figs. 2H and 2I). In modern southwest Iberia, vegetation dynamics are directly controlled by the variation in the regional influence of the westerlies, with a more northward trajectory of the westerlies seen under modern NAO+ states inducing reduced vegetation activity (Gouveia et al., 2008), a response consistent with the reduction of the southwest Iberian Mediterranean forest in the late LIG.

Salonen acknowledges funding from the Academy of Finland (project 310649). We thank Karin Helmens (Swedish Museum of Natural History) and three anonymous reviewers for comments, Anastasia Zhuravleva (GEOMAR Helmholtz Centre for Ocean Research Kiel, Germany) for help in acquiring data, and Bette Otto-Bliesner (U.S. National Center for Atmospheric Research, Colorado, USA) for making CCSM3 model output available for further analysis.

1Supplemental Material. Paleoclimate reconstructions, pollen–climate calibration data, and cross-validation results. Please visit to access the supplemental material, and contact with any questions.
Gold Open Access: This paper is published under the terms of the CC-BY license.