Dansgaard–Oeschger-like events of the penultimate climate cycle: the loess point of view

The global character of the millennial-scale climate variability associated with the Dansgaard–Oeschger (DO) events in Greenland has been well-established for the last glacial cycle. Mainly due to the sparsity of reliable data, however, the spatial coherence of corresponding variability during the penultimate cycle is less clear. New investigations of European loess records from Marine Isotope Stage (MIS) 6 reveal the occurrence of alternating loess intervals and paleosols (incipient soil horizons), similar to those from the last climatic cycle. These paleosols are correlated, based on their stratigraphical position and numbers as well as available optically stimulated luminescence (OSL) dates, with interstadials described in various Northern Hemisphere records and in GLt_syn, the synthetic 800 kyr record of Greenland ice core δ¹⁸O. Therefore, referring to the interstadials described in the record of the last climate cycle in European loess sequences, the four MIS 6 interstadials can confidently be interpreted as DO-like events of the penultimate climate cycle. Six more interstadials are identified from proxy measurements performed on the same interval, leading to a total of 10 interstadials with a DO-like event status. The statistical similarity between the millennial-scale loess–paleosol oscillations during the last and penultimate climate cycle provides direct empirical evidence that the cycles of the penultimate cycle are indeed of the same nature as the DO cycles originally discovered for the last glacial cycle. Our results thus imply that their underlying cause and global imprint were characteristic of at least the last two climate cycles.

1 Introduction

The last climate cycle (130–12 ka) was punctuated by two different types of abrupt climate changes. The first one corresponds to abrupt warmings of up to 16.5 ^∘C (Kindler et al., 2013) in about 50 years over Greenland, followed by returns to glacial conditions in several centuries. Described as Dansgaard–Oeschger (DO) events, 26 interstadials occurred over the last 115–12 ka b2k interval at varying frequencies. A recent interpretation of their occurrence is that they result from the interactions between (i) variations in the expansion of the ice shelves surrounding the three main ice sheets covering the Northern Hemisphere, especially Greenland, (ii) the expansion of the sea ice, and (iii) alterations of the Atlantic Meridional Overturning Circulation (Petersen et al., 2013; Boers et al., 2018). The DO events have also been described in variations of planktonic counts characterizing sea surface temperatures in North Atlantic cores recording the last climate cycle (Bond et al., 1992; McManus et al., 1994; Clark et al., 1999, 2007; Pisias et al., 2010). DO cycles have also been grouped into cycles of increasing cooling magnitude named Bond cycles (Alley, 1998).

The second type of abrupt climate change, called a Heinrich event (Heinrich, 1988; Broecker, 1994), is most likely related to the DO events and appeared at the end of Bond cycles (Alley, 1998). The Heinrich events are characterized by massive iceberg discharges calved from the three main ice sheets into the North Atlantic. They are identified in North Atlantic records by the occurrence of mineral detritus eroded from North America, Greenland, and northern Scandinavia, carried out by the massive icebergs and deposited at the bottom of the ocean after the melt of the icebergs. These particular iceberg detritus deposits have been identified in the North Atlantic records down to 650 ka (MIS 16) from the Hudson Strait (Hodell et al., 2008; Hodell and Channell, 2016).

The Greenland ice core warm intervals were initially correlated with the classical terrestrial interstadials described in European pollen records (Dansgaard et al., 1993). Since then, other correlations have been proposed from numerous European pollen sequences. A more detailed overview of these terrestrial records shows that the European pollen sequences, either terrestrial or marine, did indeed record numerous warm and humid intervals during the last climate cycle, which are synchronous with the Greenland and North Atlantic interstadials (Fletcher et al., 2010). During the DO events, there was temperate forest in southern latitudes as opposed to the more open shrub–tundra grassland in northern latitudes, higher than 45^∘ N.

More recently, the high-resolution study of the last climate cycle in European loess–paleosol sequences at about 50^∘ N has demonstrated that the Greenland interstadials described in the Greenland ice cores correspond one to one with the paleosols preserved in European loess (Rousseau et al., 2002, 2007, 2017a, b). Contrary to the pollen records, wherein some dating uncertainty remains, the ¹⁴C dates obtained from earthworm granules preserved in each individual paleosol described in the loess sequences perfectly confirm this correlation (Moine et al., 2017). These new results confirmed the initial hypothesis that the length of the Greenland interstadials (GISs) corresponds to the matureness of the paleosols in Europe, from boreal brown paleosols, to tundra gleys, and even embryonic tundra gley horizons. Interestingly, this relationship between GIS and European paleosol–loess alternation is rather observed in sequences located at about 50^∘ N, while it is less or not visible further southward and more difficult to evidence eastward where the environment was too dry for the formation of soils (Rousseau et al., 2011).

In the present paper, we discuss the MIS 6 paleosol–loess succession recently revealed by multiproxy high-resolution investigation of the very well-developed loess–paleosol record of Harletz in Bulgaria (Antoine et al., 2019; Lomax et al., 2019) and compare this succession with other regional (Mediterranean and European) and more global records (Northern Hemisphere). The aim of this research is to understand whether DO-like events also prevailed prior to the last climate cycle, although there is no direct record of the penultimate climate cycle in Greenland.

2 Harletz paleosol succession (Fig. 1)

The Harletz loess sequence is a 20 m thick terrestrial record located in Bulgaria on the left bank of the Ogosta River, a tributary of the Danube, at 43^∘41^′52.78^′′ N, 23^∘49^′42.27^′′ E (alt.+40 m a.s.l.; Fig. 1).

Figure 1 Harletz loess sequence. (a) Location of the western sites discussed in the text in the MIS 6 maximum ice sheet extent by Florence Colleoni from Dendy et al. (2017); the black box refers to (b). (b) Location of Harletz loess sequence (red dot) in Bulgaria and nearby loess sequences (black dots) cited in the text (modified from Antoine et al., 2019; 1: eolian sands, 2: sandy loess, 3: loess >5 m, 4: loess <5 m). (c) Pedostratigraphy and variation of some parameters (grain size, magnetism, color reflectance; modified from Antoine et al., 2019). Interstadials deduced from the measured parameters (ISm) and corresponding to incipient weathered horizons and from paleosol observation (ISp) (this study).

The multidisciplinary high-resolution investigation – including grain size, organic carbon (%), magnetic susceptibility, and spectro-colorimetry – of the loess deposits and inter-stratified soil units has revealed the superposition of two interglacial–glacial cycles below the topsoil (Antoine, 2019b). Combined with luminescence dating (Lomax et al., 2019), this approach allowed us to correlate these pedosedimentary cycles with loess cycles C and B of Kukla (1977), which are themselves correlated with Marine Isotope Stage (MIS) 7–6 and MIS 5–2, respectively.

The lower cycle is composed of a 4 m thick basal interglacial pedocomplex (Ogosta soil complex) corresponding to an upbuilding soil sequence superimposed on overbank fluvial sandy silts from the Ogosta River overlain by a first eolian series, 10 m thick, showing a detailed record of alternating loess and paleosol units. The upper cycle starts with a 2 m thick interglacial paleosol (Harletz soil complex) overlain by a loess layer, about 4 m thick, capped by a modern-day truncated soil (20 cm) (Antoine et al., 2019). The sedimentological study shows that the two eolian units are built from regional windblown silts and fine sands. The general stratigraphy identified in Harletz fits the general frame described for nearby sequences from Serbia, Romania, and Hungary by Markovic et al. (2015), including the identification of a cryptotephra layer within cycle C, mainly identified through magnetic susceptibility but also with the occurrence of a few glass shards (Sébastien Nomade, personal communication, 2019). Nevertheless, the identification of the paleosol succession at the base of the penultimate eolian sequence is particularly unique (Antoine et al., 2019) for European loess sequences. This soil sequence exhibits a succession of four soil horizons separated by many loess deposits, in which the thickness and the intensity of each pedogenesis is less and less developed towards the top (Antoine et al., 2019). Moreover, while most of the paleosols (even some ≤10 cm thick) have been observed when cleaning the outcrop, others are very incipient paleosols and/or incipient weathered horizons and could only be identified through magnetic susceptibility and spectro-colorimetry measurements supported by grain size variations (Fig. 1).

The luminescence dating obtained for the upper part of the stratigraphy above the cryptotephra layer is consistent with the nearby sequences and supports our interpretation of the MIS 7–6 and MIS 5–2 assignment for the lower and the upper part of the section, respectively (Lomax et al., 2019). A tephra layer has also been described in other loess sequences of the region (Lomax et al., 2019). However, no date has been directly measured on this particular marker, but all existing estimates indicate a minimum age of ca. −170 ka (Antoine et al., 2019). This age can be compared with two volcanic eruptions identified south of Rome, Italy, the Vico B ignimbrite and the Pitigliano tuff (Giraudi and Giaccio, 2017), as well as in nearby Lake Ohrid, where the tephra have been dated to $-\mathrm{162}\pm \mathrm{6}$ and $-\mathrm{163}\pm \mathrm{22}$ ka, respectively (Leicher et al., 2016). In Harletz, the luminescence date of the sample 0.5 m above the cryptotephra indicates an age of $-\mathrm{171}\pm \mathrm{14}$ ka, consistent with estimates of the tephra layer from Mostiştea to the east and Ruma and Stalać to the northwest of Harletz. At Harletz, the lower record of MIS 6 indicates five paleosols to incipient weathered horizons of various maturation (ISp). They are distributed along a decreasing trend in the clay content, the color reflectance, and the magnetic susceptibility but along an increasing trend in the coarse sand and the grain size ratio (coarse + fine sand–clay) (Fig. 1).

Just above the cryptotephra layer, there are two more thin incipient weathered horizons only observed through the analytical data (cryptopaleosols, ISm). Moreover, ongoing magnetic studies (isothermal remanent magnetization and coercivity of remanence) identified three more incipient weathered horizons on top of the loess unit at 6.2, 7.4, and 8.4 m (Lagroix et al., 2016). Therefore, considering the Lake Ohrid recalculated dates of the tephra layer and the luminescence age available, there are 10 horizons including paleosols and incipient weathered horizons identified above the MIS 7 Ogosta interglacial pedocomplex. In the MIS 6 loess section, there are seven and three horizons in the lower and upper parts, respectively, assigned to the interval between −190 and −130 ka. These paleosols and incipient horizons are interpreted to correspond to interstadials of various duration in the continental record by reference to those that have been described in the last climate cycle (Antoine et al., 2009; Rousseau et al., 2002, 2007, 2017a) (Fig. 1).

Although unique among the other loess sequences from the area, can these paleosols and incipient weathered horizons and interstadials be related to more global events as interpreted for last climate cycle paleosols?

3 Proxy records in the Mediterranean region

3.1 Loess records (Fig. 2)

Numerous loess sequences within the Carpathian and lower Danube basins show the identified penultimate loess unit assigned to MIS 6 or early MIS 5e, either based on luminescence ages or on magnetic susceptibility variations, in Bulgaria (Jordanova and Petersen, 1999; Jordanova et al., 2008), in Serbia (Fuchs et al., 2008; Markovic et al., 2009; Murray et al., 2014), in Ukraine (Buggle et al., 2009), in Romania (Balescu et al., 2010; Timar-Gabor et al., 2011; Vasiliniuc et al., 2012), and in Hungary (Novothny et al., 2010, 2011; Ujvari et al., 2014). Some of them show the occurrence of a tephra layer, e.g. in Batajnica, Ruma, and Stalać in Serbia (Buggle et al., 2009; Markovic et al., 2009; Novothny et al., 2010, 2011; Obreht et al., 2016; Ujvari et al., 2014; Vandenberghe et al., 2014) and at Mostiştea in Romania (Balescu et al., 2010; Panaiotu et al., 2001). Markovic et al. (2015) did not assign this tephra to any particular eruption, but considering the age model proposed, it rather corresponds to either the Vico B ignimbrite or the Pitigliano tuff (Giraudi and Giaccio, 2017) at about −170 ka. Furthermore, in the review of Danube loess sequences, Markovic et al. (2015) also refer to an embryonic pedogenic layer older than the identified −170 ka tephra but younger than the first loess subunit identified at the start of the L2 loess unit correlated with MIS 6, above the interglacial paleosol S2SS1. Because of its position in the stratigraphy scheme, this embryonic paleosol is correlated with the first paleosol identified in Harletz between about −15 and −14.5 m of depth (Antoine et al., 2019; Fig. S1 in the Supplement). The other nine paleosols identified in Harletz have not been mentioned or described westwards in the Serbian and eastwards in the lower Danube loess sequences, making this highly developed loess record of MIS 6 a remarkable reference sequence for this region (Antoine et al., 2019).

Figure 2 Correlation between the last two climate cycles from the southern Europe loess sequences along the Danube, modified after Antoine et al. (2019). Map of the locations (see Fig. 1b); the colored rectangles on the map (left) refer to the pedostratigraphic sequences correlated on the correlation sketch (right). Interstadials ISp and ISm (this study; see Fig. 1c) are reported on the Harletz record.

3.2 Pollen records (Fig. 3, Table 1)

Located at 40^∘54^′–41^∘10^′ N and 20^∘38^′–20^∘48^′ E in the Balkan Peninsula within the Dinaride–Albanide–Hellenide mountain belt, Lake Ohrid is a transboundary (Macedonia and Albania) lake. This site is the closest pollen record to the Harletz sequence, where two different open vegetation zones represent MIS 6, namely OD-5 (−190 to −160 ka) and OD-4 (−160 to −129 ka) (Sadori et al., 2016). The former corresponds to a grassland-dominated environment and the latter to a rather steppe-dominated environment. A tephra layer also marks the boundary between the two vegetal formations, indicating that OD-5 and OD-4 could be correlated with the lower and the upper parts of the Harletz MIS 6 record, respectively (Antoine et al., 2019). Considering the variations in both the percentage of total arboreal (AP) or total arboreal minus pine + juniper + birch pollen grains in OD-5, only four tree expansion phases can be noticed, which could correspond to four out of the five interstadials identified in the lower part of the Harletz MIS 6 record. Their onsets are dated, according to the age model used (Sadori et al., 2016), at about −185, −178, −175, and −169 ka below the tephra layer. In the following OD-4, four other interstadial onsets are dated at about −159, −151, −145, and −139 ka.

Figure 3 Western (ODP 977 SST – alkenone; Martrat et al., 2004a, b) to eastern (Soreq speleothem δ¹⁸O; Bar-Matthews et al., 2003a, b) Mediterranean records of MIS 6, including a pollen record from Ioannina I-284 (AP, AP − pine + birch + juniper; Roucoux et al., 2011), compared with summer insolation with an indication of the time occurrence of Sapropel 6 (Bard et al., 2002). The pedostratigraphy of the Harletz MIS 6 record is added with the newly described paleosols and pedogenic horizons for comparison. Labels of the interstadials are as identified in the Alboran Sea. Complementary labels assign warm events to the ODP 977 (AI-5'a, b, AI-9'a, b), Ioannina (I1–12), and Soreq (S1–9) Mediterranean MIS 6 records used in Table 1. MIS 6 ice sheet extent by Florence Colleoni in Dendy et al. (2017).

Table 1 Synchronism of the various interstadials discussed in this study. Indication of the predicted DO event occurrence by Barker et al. (2011) for the various assigned interstadials from Sanbao11 and Soreq speleothem δ¹⁸O, ODP 984, ODP 983, U1308 remaining percentage in N. pachyderma s., ODP 977 and MD01-2443 sea surface temperature (SST), and Ioannina I-284 arboreal pollen variations. Marine Isotope Stage (MIS) 6 is decomposed into an early part corresponding to the occurrence of the Harletz paleosols (192–170 ka), and a later part (170–130 ka) with incipient weathered horizons.

At lower latitude, the Tenaghi Philippon site is located in Greece at 41^∘10^′ N, 24^∘20^′ E; 40 m a.s.l. This is a pollen record famous for its long continental climate record covering about eight climate cycles (Tzedakis et al., 2006). This record indicates woody taxa expansions, which can be considered interstadials, during the interval equivalent to MIS 6. Following the age model used (Tzedakis et al., 2006), they occurred at about −185, −180, −171, −163, −153, −151, and −138 ka. Therefore, some of them seem potentially synchronous with the tree expansions described at higher resolution in Lake Ohrid and by extension to Harletz interstadials. Another long pollen core retrieved in Greece, the Ioannina I-284 core (39^∘45^′ N, 20^∘51^′ E; 472.69 m a.s.l.), which covers several climate cycles as well, reveals the record of MIS 6 in the section −138 to −99 m. Roucoux et al. (2011) have described the high-resolution variations of the vegetation during this time interval. Phases of temperate tree expansion, interpreted as interstadials, and phases of tree contraction, interpreted as stadials, have been identified. Based on the employed age model (Roucoux et al., 2011), these temperate tree expansions occurred at about −186, −181, −177, −173, −170, −165, −159, −156, −145, −138, and −135 ka that we labeled I11 to I1 (Fig. 3). They appear more numerous than in Tenaghi Philippon and in Lake Ohrid, probably due to the higher resolution of the pollen analysis of the Ioannina core, and similarly some of these tree expansions could be correlated with the Harletz MIS 6 interstadials, five expansions being recorded older than −165 ka (Table 1).

3.3 Marine records (Fig. 3, Table 1)

On both sides of the Strait of Gibraltar, Martrat et al. (2004a, 2007a) have described the climate variations expressed during at least the past two climate cycles from cores retrieved from the Alboran (ODP 977A; 36^∘1.9^′ N, 1^∘57.3^′ W; 1984 m below sea level) and the Iberian seas (MD01-2443; 37^∘52.85^′ N, 10^∘10.57^′ W; 2925 m below sea level). The analysis of δ¹⁸O from the planktonic foraminifera Globigerina bulloides shows variations with abrupt changes mimicking – for the last climate cycle (MIS 5 to 2) – those described in the Greenland ice core records. During the penultimate cycle (MIS 7–6, −245 to −130 ka) the foraminifera δ¹⁸O record also shows abrupt warmings similar to the DO interstadials, which have also been observed in the variations of the ${\mathrm{U}}_{\mathrm{37}}^{{\mathrm{k}}^{\prime }}$ alkenone index, a proxy for the sea surface temperature. Nine interstadials are therefore identified, named Alboran or Iberian margin interstadials, although the older one, AI-9', seems to be a composite in which three different events can be discriminated. The onsets of these penultimate cycle interstadials are dated at about −186 (−186, −182, −179), −176, −170, −165, −159, −152, −150, −141, and −133 ka. They were labeled AI-9', AI-8', AI-7', AI-6', AI-5', AI-4', AI-3', AI-2', and AI-1', respectively. Similarities in terms of the number of warming events appear with the long Mediterranean pollen records mentioned previously. Although no tephra is identified, five events can once more be identified prior to −165 ka, as also observed in Harletz (Fig. 3, Table 1). Martrat et al. (2007a) stated that the observed variations in sea surface temperature, expressed by ${\mathrm{U}}_{\mathrm{37}}^{{\mathrm{k}}^{\prime }}$ for the last four climate cycles in the western Mediterranean, show a nonlinear response to external triggers of climate: obliquity and precession. However, some interstadials are synchronous with Mediterranean sapropels, which are a direct response to orbital forcings. Rohling et al. (2015) presented an updated review of the present and past Mediterranean climate and oceanography with a clear differentiation between the western and eastern basins, which may have impacted the adjacent terrestrial regions.

3.4 Speleothem records (Fig. 3, Table 1)

Two Mediterranean speleothems provide records for MIS 6. The Argentarola Cave, located (42^∘23^′8.79^′′ N, 11^∘3^′59.17^′′ E) by the Tyrrhenian Sea, shows an interval of low δ¹⁸O values between −180 and −165 ka with minima at −178, −172, and −168 ka, the latter being in phase with high-latitude insolation. The −180 to −165 ka interval is interpreted as corresponding to the Mediterranean Sapropel S6 (Bard et al., 2002) (Fig. 3). The recorded variations are correlated with total organic content variations recorded in the two Mediterranean cores MD84641 and KC19C, indicating rather warm and pluvial conditions during this interval over the western Mediterranean. Eastward, in Soreq Cave in Israel (31.458^∘ N, 35.038^∘ E; 400 m a.s.l.), a speleothem record covering the last 180 kyr was obtained. The measured δ¹⁸O record of MIS 6 shows variations at about −178, −170, −165, −160, −152, −149, −144, −138, and −134 ka (Fig. 2). The peaks at −178 and −152 ka show very low values, interpreted as corresponding to intense pluvial intervals over southern Israel, with the former correlated with the start of Sapropel S6, while no particular sapropel is associated with the latter (Ayalon et al., 2002; Bar-Matthews et al., 2003a), although it is referred to as a monsoon index maximum (Melières et al., 1997). These records indicate high precipitation rates at least during the −180 to −165 ka interval, which also corresponds to the most clearly expressed interstadials and paleosols in Harletz. Describing Soreq and another record in northern Israel, Bar-Matthews et al. (2003a) concluded that these speleothems recorded the proxy signal of global and regional eastern Mediterranean climate over the last 250 kyr. As Fig. 3 shows, the age model of Soreq is based on U−Th dating and agrees with the Bard et al. (2002) results, while the marine and continental cores are partly orbitally tuned.

Although the interstadials identified in the Harletz MIS 6 sequence have not been described in the closest loess sections from the Carpathian Basin and lower Danube region, the overview of the various and diversified MIS 6 records presented above clearly indicates that these paleosols and incipient weathered horizons correspond to climatic events that are recognized all along the northern Mediterranean and therefore have more than a local to regional significance. This shows that the climate mechanism proposed to explain the paleosol–loess alternations of the last climate cycle (Antoine et al., 2009; Boers et al., 2017, 2018; Rousseau et al., 2002, 2007, 2017a, b), in northern European sequences at about 50^∘ N, seems to have already prevailed during the penultimate climate cycle. This suggests that these alternations are part of more global dynamics. Therefore, in the next step of our study we will look for similar events in other records outside the Mediterranean region at a more global scale.

4 DO-like events during MIS 6 in different regions (Fig. 4, Table 1)

The present investigation of the Harletz MIS 6 loess records reveals the occurrence of alternating loess intervals and paleosols, similar to those from the last climatic cycle in Europe. These paleosols are correlated with interstadials described in various Mediterranean records and correspond to warmer and moister events, which may also correlate with other Northern Hemisphere interstadials. In northwestern Europe, MIS 6 tundra gleys have also been described in loess sequences. Locht et al. (2016) report three tundra gleys and one paleosol in northern France. In Belgium, Juvigne et al. (1996) describe five tundra gleys in Kesselt, while Pirson et al. (2018) report two tundra gleys in Albert Canal, and Haesaerts et al. (2016) mention one and two tundra gleys in Harmignies and Remicourt, respectively. Finally, in the lower Rhine valley in Germany, Schirmer (2010), defining a synthetic Rhine pedostratigraphical record, refers to eight tundra gleys over two thin MIS 6 paleosols (a basal humus zone and a calcareous cambisol). However, the precise dating of all these tundra gleys throughout MIS 6 remains unclear without any particular investigation and therefore prevents correlations as defined for the last climate cycle.

Figure 4 Comparison of Northern Hemisphere records of MIS 6. Greenland synthetic δ¹⁸O (Barker et al., 2011), North Atlantic SST (N. pachyderma s.) (cores ODP 984, Mokeddem and McManus, 2016; ODP 983, Barker et al., 2015; U1308, Obrochta et al., 2014a, b), Iberian margin and Alboran Sea SST (alkenone MD01-2443, ODP 977; Martrat et al., 2004a, b, 2007a, b), and Chinese speleothem δ¹⁸O (Sanbao11; Wang et al., 2008). Labels of the interstadials are as identified in the marine and speleothem records. Green dots denote predicted DO event occurrence from Barker et al. (2011). The pedostratigraphy of the Harletz MIS 6 record is added with the newly described paleosols and pedogenic horizons for comparison. Complementary labels assign warm events to the ODP 984 (IS6-3a, b), U1308 (1–8), and ODP 977/MD01-2443 (AI-5'a, b; AI-9'a, b) North Atlantic and Mediterranean MIS 6 records used in Table 1. MIS 6 ice sheet extent by Florence Colleoni in Dendy et al. (2017).

The Greenland ice cores are the key paleoclimatic references used to interpret the paleosol–loess unit alternations identified during the last climate cycle in the European loess sequences at about 50^∘ N. They all show a very precise imprint of the abrupt changes over the last 130 kyr. Among them, the NGRIP δ¹⁸O record has an extremely high resolution, although layer-counting errors accumulate down the core (Boers et al., 2017), yielding suitable dates for the onset and end of the abrupt warmings as observed in various parameters (Fischer et al., 2015; Schupbach et al., 2018), and therefore allows for precise correlations with the terrestrial sequences (Moine et al., 2017; Rousseau et al., 2002, 2007, 2017a, b).

Although initially described in the GRIP record (Dansgaard et al., 1993), no high-resolution record of the penultimate climate cycle is available from the Greenland Ice Sheet. Moreover, the Antarctic ice core records are not directly usable in a first step to compare them with the European loess sequences (see below). A potential solution is therefore given by considering the long high-resolution records from Chinese speleothems. Indeed, the Hulu Cave δ¹⁸O record has been demonstrated to perfectly correlate with the NGRIP record of the last climate cycle (Wang et al., 2001; Cheng et al., 2006) through precise ²³⁰Th ages. Moreover, this precise and detailed climate record has been extended back in time until −224 ka b2k through the addition of the measured δ¹⁸O performed on the Sanbao cave speleothem SB11, which has a resolution similarly high as that of the Hulu Cave but for the penultimate climate cycle. This is the reason why we use this particular record as a reference for our comparison.

4.1 Sanbao11 (Figs. 4, S1)

The Sanbao cave is located at 110^∘26^′ E–31^∘40^′ N, 1900 m a.s.l., in central China, south of the Chinese Loess Plateau. It is strongly influenced by the East Asian monsoon variations and correlates with the summer insolation gradient between 65 and 15^∘ N. Moreover, similarly to the last climate cycle record measured from the Hulu Cave, the δ¹⁸O record of Sanbao speleothem B11 indicates abrupt changes at a millennial scale during MIS 7 and 6, which have been interpreted as interstadials by Wang et al. (2008). Considering the age boundary between MIS 7 and 6 (Lisiecki and Raymo, 2005) at 191 ka, the Sanbao speleothem preserved five interstadials, named B24 to B20, in MIS 7 (although missing the lower 19 ka; limit at −243 ka) and 13, named B19 to B1, in MIS 6. The δ¹⁸O minima in MIS 6 are dated at −189 ka for B19, −181 ka for B18, −175.5 ka for B17, −172 ka for B16, −168 ka for B15, −166 ka for B14, −162 ka for B13, −160 ka for B12, −157 ka for B11, −151 ka for B10, −149 ka for B9, −147.5 ka for B8, −145 ka for B7, −136 ka for B2, and −134 ka for B1 (Fig. 3); note that B3 to B6 are better recorded in the Hulu Cave MSP speleothem (Cheng et al., 2006). According to the prevailing interpretation, the lower the δ¹⁸O values from the speleothem record, the more intense the East Asian summer monsoon intensity. Interestingly, the minimum δ¹⁸O values for B17, B16, and B15 in MIS 6 are much lower than the other minimum values identified in the past 224 kyr record, including those characterizing the past two interglacials (MIS 7 and 5e) and the Holocene one. Hence, the time interval between −178 and −165 ka must have corresponded to extremely high precipitation rates over the East Asian region. This observation is in accordance with the high precipitation interval derived from the speleothems from Soreq and Argentarola caves in the Mediterranean region during Sapropel 6. Interestingly, the Sanbao cave record yields five interstadials between −190 and −167 ka and hence the same number of interstadials as the number of paleosols observed in Harletz below the tephra layer. Moreover, the B21, B16, and B8 δ¹⁸O minima correspond to peaks of decreasing magnitude in summer insolation at 15^∘ N. However, although summer insolation at both 15^∘ N (tropics) and 65^∘ N (arctic circle) show maximum values corresponding to the Sanbao δ¹⁸O main minima, the insolation gradient (15–65^∘), which induces a potentially strong meridional heat transport in the atmosphere and a possible excess in precipitation over the continent, indicates the weakest maximum for the time interval B17–B15 (about −177 to −166 ka) (Fig. 4).

4.2 Greenland GL_t_syn (Figs. 1, 3, Table 1)

Using the argument of the bipolar seesaw interpretation of the antiphased climate variation prevailing between Antarctica and Greenland, Barker et al. (2011) have developed a synthetic δ¹⁸O record of Greenland for the last 800 kyr based on the EPICA results and timescale, named GL_t_syn. This synthetic record nicely replicates the millennial-scale variability observed during the last climate cycle in the NGRIP ice core and therefore makes the authors confident that it reliably reproduces the last 800 kyr of Greenland δ¹⁸O abrupt variability. Two timescales have been proposed: one based on EPICA DC3 and another based on the Sanbao speleothem, with some discrepancies. Because of the more precise dates obtained from the Chinese speleothems, we will use the GL_t_syn with the Sanbao timescale for our comparisons. Plotted versus the Sanbao δ¹⁸O record, the GL_t_syn record indicates intervals showing peak maximum values aligned with most of the Chinese speleothem interstadials. Among these peak values, Barker et al. (2011) predicted the occurrence of 11 DO-like events, which should correspond to as many interstadials (Fig. 4, Table 1). Five of them can be clearly identified between −190 and −166 ka, which align with the Sanbao ones, therefore supporting the global value of Harletz paleosols interpreted as interstadials.

4.3 North Atlantic (Figs. 1, 4, S1, Table 1)

South of Greenland, the North Atlantic Ocean is the classical region in which to observe the Greenland warmings and interstadials, as this was demonstrated for the last climate cycle (Bond et al., 1992; McManus et al., 1999; Henry et al., 2016; Hodell and Channell, 2016).

For the penultimate cycle, several records also indicate warming events interpreted by the respective authors as interstadials. The warming events observed in the sea surface temperature (SST), deduced from alkenone studies from the Iberian margin, are of similar magnitude as those observed in the Alboran Sea (ODP 977/MD12443) (Fig. 3). Furthermore, the Iberian margin core MD01-2444 shows temperate tree pollen percentages, which vary in line with the planktonic δ¹⁸O measured in the same core (Margari et al., 2014). Such high percentages of arboreal pollen during the base of MIS 6 (−185 to −160 ka) are interpreted in terms of reduced summer aridity and increased pluvial conditions. Back to the sea surface temperature, which has been interpreted as mimicking the Greenland δ¹⁸O variation (Shackleton et al., 2000), plotting ODP 977 alkenone SST reconstructions of MIS 7 and 6 against GL_t_syn, one can notice that numerous discrepancies appear in the number of interstadials and for those identical in their occurrences, which seems to be due to the age model used for ODP 977. Further north, in core U1308 (49^∘53^′ N–24^∘14^′ W, 3871 m n.s.l – near sea level) Obrochta et al. (2014a) report important variations in Neogloboquadrina pachyderma sinistral, an indicator of cold surface water, during MIS 6. The minimum counts are much stronger at the base of the MIS 6 record, between −190 and −160 ka, than in its upper part. Although dating these warming events can be proposed from the published material, at about −190, −175.5, −167.5, −160, −166.5, −143, and −137 ka, these dates remain preliminary as the time resolution is not high enough compared to GL_t_syn or ODP 977/MD01-2443 marine cores (Fig. 4). The M23414-9 gravity core (53.537^∘ N, 20.288^∘ W; water depth 2199 m), selected to study variations in the North Atlantic drift, indicates the highest values in N. pachyderma s. percentages between −165 and −130 ka. However, minima correspond to warmings in both summer and winter sea surface temperature estimates (Kandiano et al., 2004). A total of 11 of such events can be dated at about −192, −187.5, −180.5, −175.5, −172, −163.5, −158.5, −141, −140, −137, and −133.5 ka BP, but they show very little overlap with the other marine cores previously mentioned, probably because of a lower resolution. Still, the five lowest occurred before −165 ka. A similar pattern is noticed further north in ODP 983 (60.4^∘ N–23.6^∘ W, 1984 m n.s.l), with stronger minimum counts of N. pachyderma s. between 190 and 160 ka than between −160 and −130 ka. Using the GICC05/NALPS/China timescale (Barker et al., 2015), warming events are observed at about −189.9, −178.7, −173.2, −169.5, −168.2, −163, −160.7, −151, −149.1, −147.4, −140.4, and −134.5 ka (Fig. 4). These dates vary slightly when using the other proposed timescale based on the EDC3 age model. Closer to Greenland and south of Iceland, ODP 984 (drilled at 61.25^∘ N and 24.04^∘ W, 1648 m) yielded a record of the past 220 kyr (Mokeddem and McManus, 2016). It includes the complete MIS 6, showing millennial-scale variability that is best expressed during the interval −170 to −130 ka. The foraminifer composition indicates alternation between episodes of northward polar-front retreats (low values of N. pachyderma s.), synchronous with warm and salty water inflows (high values of Turborotalia quinqueloba), and episodes of southward polar-front advance (high values of N. pachyderma s.), with fresh and cold water inflows (low T. quinqueloba) and high ice-rafted debris (IRD) content. The cold inflow episodes are interpreted as stadials, while the warm ones are considered interstadials equivalent to those observed during the last climate cycle. They occurred at about −165 to −163, −156 to −152, −148 to −146, −139 to −137, and −135 to −132 ka, and they are named IS6-5, IS6-4, IS6-3, IS6-2, and IS6-1, respectively (Fig. 4). Interestingly, these interstadials are the nearest geographically described to Greenland during MIS 6, but considering the resolution, more events could be considered when looking at the percentage of N. pachyderma s. (see Fig. 4 and Table 1).

When comparing the MIS 6 records described previously, e.g., the reconstructed Greenland δ¹⁸O, the Sanbao speleothem δ¹⁸O, and the foraminifer counts and δ¹⁸O measurements from Sites 984 and 983, one can notice that the interstadials observed in the North Atlantic records can have an equivalent in both the reconstructed Greenland synthetic and Sanbao δ¹⁸O interstadials (Fig. 4). Plotting every record used in the present study on its individual timescale, a global synchronism is far from being evident, except for the early MIS 6. When reconstructing the Greenland δ¹⁸O variations, Barker et al. (2011) predicted DO event occurrences for the last climate cycle, which fit with the DO events described from the Greenland ice cores. They expanded the occurrence prediction to the previous climate cycles covering the past 800 kyr. Plotted versus the Sanbao interstadials, although these predicted DO-like events fit for the early part of MIS 6, between −192 and −166 ka, fewer abrupt changes are assigned DO-like labels than in the Sanbao record of interstadials for the −166 to −130 ka interval, especially in the upper (more recent) MIS 6. This is demonstrated when comparing the reconstructed Greenland δ¹⁸O, the variations in N. pachyderma s. percentages in the North Atlantic ODP 984, the closest to Greenland, and the U1308 records. Reporting all the interstadials in Table 1 supports the previous interpretation that these interstadials are hemispheric and can be interpreted as DO-like events, without any doubt at least for the early MIS 6. However, the question remains as to why these marine records did not show all of the interstadials described in the other domains, contrary to what has been observed for the last climatic cycle. Is this related to the time resolution of the studied cores, which would not be high enough? This is an ongoing problem that this paper cannot address, as it requires further investigation and additional data with higher resolution.

4.4 Chinese loess records (Fig. S1)

One potential record of all these events has been proposed from sequences from the Chinese Loess Plateau, located north of the high-resolution speleothems described previously.

For the last glacial cycle, the study of the grain size variations from temporally highly resolved loess sequences indicates changes in the size of the deposited particles ,which are aligned, using a timescale based on luminescence dates, with the DO events (Sun et al., 2012). As the deposited material originates from northern Chinese deserts, during interstadials, when the wind velocity is reduced due to stronger East Asian monsoon blowing from the south, finer material is deposited. In contrast, during stadials, the wind velocity is stronger and the East Asian summer monsoon weakens, allowing coarser material to deposit. Following this idea, Yang and Ding (2014) have proposed a stack of several loess sequences covering the last two climate cycles (Fig. S1). They show that for the last 130 kyr, the median grain size variations match the former results by Sun et al. (2012) and correlate with the Hulu Cave speleothem δ¹⁸O record (Cheng et al., 2006; Wang et al., 2008). Expanding to the penultimate cycle, Yang and Ding (2014) indicate that the median grain size variations during MIS 7 and 6 also mimic the Sanbao cave δ¹⁸O variations within the uncertainties related to the dating of these sequences. Interestingly, the authors refer to interstadial paleosols that do not correspond to the identified individual interstadials but rather to a group of them (L2-4 groups B17, B16, and B15, while L2-2 gathers B11 to B5). This observation is similar to what was described for the last climate cycle record by Sun et al. (2012): not every interstadial is individually resolved in the stratigraphy, contrary to the various paleosols in the European loess sequences at about 50^∘ N, which correspond one to one with the Greenland interstadials.

5 Concluding comments

The identification of interstadials and paleosols in the Harletz MIS 6 record, which were not known previously in the loess records of the nearby sequences, led us to evaluate the significance of these warm and moist episodes. Such events appear to have equivalents in the northern Mediterranean region and are expressed through different parameters. Interestingly, some of these events correspond to very humid episodes as recorded in speleothems, corresponding to the deposition of Sapropel 6 in the Mediterranean Sea, associated with a northward shift in the position of the summer Intertropical Convergence Zone (ITCZ). North Atlantic cores also yield evidence of these interstadials, although the available records are not complete except in the vicinity of the Strait of Gibraltar. In the Iberian Sea, similar events in the SST reconstructions based on alkenones as in the Alboran Sea can be observed. These North Atlantic MIS 6 interstadials are also associated with lower benthic δ¹³C values, inducing the continental uptake of ¹²C through the development of vegetation increase, in agreement with synchronous arboreal peaks in the pollen records and synchronous soil development.

At a broader geographic scale, these warm events are also well-recorded synchronously in Chinese speleothems and the Greenland synthetic δ¹⁸O variations of MIS 6. Therefore, contrary to the initial interpretation as regional events, the ISm and ISp interstadials mostly expressed in Harletz provide evidence of events at least at the Northern Hemisphere. In Greenland, they correspond to strong warmings, and in the North Atlantic they correspond to abrupt SST increases, partly due to the development of incipient weathered horizons and paleosols in European loess, arboreal vegetation in European (mostly Mediterranean) pollen records, minimum values in the δ¹⁸O from Mediterranean and Chinese speleothems, and maxima in Chinese stack median grain size. There is therefore a similarity between these interstadials described from MIS 6 records in the Atlantic, Asian, and European regions, including the Mediterranean area, and those described for the last climate cycle in the same regions.

In summary, we argue that there are similarities between the well-established, globally synchronous interstadials of the last glacial cycles and corresponding episodes of the penultimate glacial cycle. This is true for the local character of these episodes in the different kinds of records, ranging from marine and lake sediments to ice cores and loess sequences, but there is rising evidence that it is also true for the global synchronicity of the interstadials identified in the different records. Clearly, the available data for the penultimate cycle are substantially sparser, have much coarser resolution, and are overall less reliable. Nevertheless, taken together, the available empirical evidence suggests that the MIS 6 interstadials described in the various environments, and correlated in a similar way as for the last 130 kyr, i.e., the various paleosols described at the base of MIS 6 in Harletz and corresponding to the interval −192 to −166 ka, can confidently be interpreted as DO-like events of the penultimate climate cycle. Moreover, in addition to the four DO-like events deduced from direct paleosol observations, six more were deduced from reliable proxy measurements, which also find equivalents in other Northern Hemisphere records of MIS 6. These are therefore 10 DO-like events that we identified in our record of MIS 6. With improved resolution of the records for the penultimate glacial, particularly more detailed dating, a thorough statistical analysis of the synchronicity of interstadials in the different records may become possible in the near future.

At least the last two climate cycles, and probably further previous ones, are likely to have experienced millennial-scale variability around the globe, which cannot be directly explained by any astronomical forcing. Such climate variability may have been induced by forcings superimposed to the classical orbital parameters inducing the climate cycles and is most likely caused by self-sustained oscillations induced by interactions between ice cover and ocean circulation changes.