Eurasian contribution to the last glacial dust cycle

. The last 130 000 years have been marked by pro-nounced millennial-scale climate variability, which strongly impacted the terrestrial environments of the Northern Hemisphere, especially at middle latitudes. Identifying the trig-ger of these variations, which are most likely associated with strong couplings between the ocean and the atmosphere, still remains a key question. Here, we show that the analysis of δ 18 O and dust in the Greenland ice cores, and a critical study of their source variations, reconciles these records with those observed on the Eurasian continent. We demonstrate the link between European and Chinese loess sequences, dust records in Greenland, and variations in the North Atlantic sea ice extent. The sources of the emitted and transported dust material are variable and relate to different environments corresponding to present desert areas, but also hidden regions related to lower sea level stands, dry rivers, or zones close to the frontal moraines of the main Northern Hemisphere ice sheets. We anticipate our study to be at the origin of more sophisticated and elaborated investigations of millennial and sub-millennial continental climate variability in the Northern Hemisphere.


Introduction
During the last glacial interval, abrupt climate changes have been documented worldwide in different types of records, but especially in ice cores (Dansgaard et al., 1969(Dansgaard et al., , 1982Johnsen et al., 1992Johnsen et al., , 1972Johnsen et al., , 2001. Their interpretation mostly focused on the more spectacular character corresponding to abrupt warmings (named Dansgaard-Oeschger -DO events;Broecker et al., 1990) of some ten degrees to interstadial conditions in Greenland (Kindler et al., 2014). These warmings were followed by a schematic two-step return to stadial conditions. Modeling experiments are able to reconstruct this abrupt warming and the two-step return to stadial conditions, indicating a periodicity of about 1500 years (Schulz, 2002;Rahmstorf, 2003), which, however, still remains questionable (Ditlevsen et al., 2007;Thomas et al., 2011;Boers et al., 2017). More precisely, the very-high-resolution analysis of the last deglaciation in the North Greenland Ice Core Project (NGRIP) record (NGRIP, 2004) revealed that different parameters show different abruptness of the warming events. For the last two warming events (14.7 and 11.7 ka b2k (years before AD 2000)), deuterium excess increased within 1 to 3 years, alongside more gradually increasing temperatures as represented by δ 18 O. Dust concentration would require at least 15 years to decrease, preceding the deuterium excess increase by 10 ± 5 years (Steffensen et al., 2008). Care-fully considering the DO events it was noticed that during the last climatic cycle in NGRIP, some variability that occurs as the last deglaciation timing of transitions does not seem to be reproduced every time in the older parts of the record . The climate trend, cooling and increase in dustiness, within these particular events is variable as well. High-resolution studies of the temperature signal in older interstadials show the occurrence of sub-millennialscale elements like precursor events of about centennial duration before the interstadial itself, rebound events exhibiting abrupt cooling towards stadial conditions, and cooling events occurring towards the end of the interstadial (Capron et al., 2010). These sub-millennial events make the understanding of the climate variability during these interstadials even more complicated than a simple warming followed by a two-step cooling. Ice cores nevertheless provide much more information than on temperature and dust concentration only, as they release records of numerous components of the climate system, isotope values of the transported water vapor, mineral aerosols and greenhouse gas concentrations, chemical elements, etc., which show different origins and transport patterns to the high-latitude ice sheets. Such richness in proxies allows comparisons with other records of millennial-scale variability preserved in both marine (Henry et al., 2016) and terrestrial deposits, as well as in other ice cores (Barbante et al., 2006;Buizert et al., 2015). In this paper, after briefly describing the abrupt changes observed in the very-high-resolution δ 18 O and dust records from NGRIP (NGRIP, 2004), we compare the dust particle sedimentation rates over Europe and China as expressed in key loess sequences ( Fig. 1). This is a complementary study of Rousseau et al. (2017), which essentially focused on paleosol development. In the last section we provide an interpretation of the link between Greenland dust records and Eurasian loess sequence development.

Analysis
We investigate in this study the 18 Greenland interstadials (GIs), labeled from 17.1 to 2 (Rasmussen et al., 2014), which have been identified during the MIS 3 and 2 (59-14 ka b2k) interval. We performed this comparison by studying in parallel the δ 18 O (NGRIP, 2004; and the dust (Ruth et al., 2003) records from the NGRIP ice core at the highest resolution. These two indices correspond to different sources, marine and continental but also to different origins, mostly Atlantic  and East Asian (Biscaye et al., 1997;Svensson et al., 2000;Bory et al., 2002). The two indices are not directly related so that if similarities would appear, this should highlight a more global phenomenon.
The compilation of the 18 GIs in the NGRIP record  shows that DO events can be characterized using a numerical algorithm (see Supplement) by an increase in δ 18 O occurring on average in 36.4 ± 13.4 (1σ ) years, with a mean GI duration of 1048 ± 1163 (1σ ) years (Table 1). When determined visually by considering the onset of the abrupt change at the start and the return to the same initial value as the end of the event, the δ 18 O changes occur, on average, in 55.4 ± 16.1 (1σ ) years, with a mean duration of 1053 ± 1068 (1σ ) years. The larger errors than the mean values are due to the few long intervals. These two methods have been described in detail in Rousseau et al. (2017). Although both methods allow the identification of the GIs and their parameters (onset, transition, and duration), the numerical approach also allows the identification of statistically significant events that the virtual one could consider as questionable. However, we consider that these two methods are complementary, the reason why we decided to release values obtained by both of them. These compiled characteristics fit with the values generally considered from the literature (Wolff et al., 2010;Rasmussen et al., 2014).
Following the detailed correlations defined between Nussloch and NGRIP stratigraphies ( Fig. 2) (Rousseau et al., 2002(Rousseau et al., , 2007b, we then applied the dates obtained for every start and end of a GI, from both the δ 18 O and the dust NGRIP timescale, to the loess sequence. This was performed by considering that in western European loess sequences, paleosols developed from the underlying loess deposits after a stop of the eolian sedimentation (Taylor et al., 2014), and that the eolian sedimentation itself restarted on top of the developed paleosols. This makes the time evolution nonlinear and a bit more complex than the classical continuous sedimentation observed in other terrestrial, marine, and ice-core records (Kukla and Koci, 1972;Rousseau et al., 2007a). Therefore, a determined eolian interval, equivalent to the Greenland stadial (GS), includes the loess unit and the overlying paleosol (blue arrow in Fig. 3), while the paleosol development itself fits with the GI duration (red arrow in Fig. 3). Doing so, the Nussloch stratigraphy can be read as expressed in Table 2, allowing then to better estimate the sedimentation and the mass accumulation rates required for comparison with other loess records and model outputs. Conversely, in Chinese loess sequences no such paleosols developed during the last climate cycle. Therefore, own reasoning does not apply to these records.

Discussion
The mineral dust record in the NGRIP ice core is obtained from variations in dust concentration, measured in the terms of the number of particles larger than 1 µm mm −1 of melt water, which also shows abrupt changes that are quite synchronous with the DO events expressed in the δ 18 O record (Mayewski et al., 1994;Ruth et al., 2003;NGRIP, 2004;Rasmussen et al., 2014). Abrupt decreases in the dust concentration occur in 60 ± 21.2 (1σ ) years on average, with a mean GI duration of 1079 ± 1135 (1σ ) years when deter-  Table 1. NGRIP ice-core record. Estimates of NGRIP GI start, end, and duration and GI abrupt transition start, end, and duration (from Rousseau et al., 2017). Temperature reconstruction for the GI transitions as published by Kindler et al. (2014). Panels (a) and (b) show NGRIP transition dates determined visually and algorithmically, respectively.  mined with the same algorithm as for δ 18 O; when determined visually, the dust decrease occurs in 56.8 ± 19.6 (1σ ) years on average, with a mean duration of 1079 ± 1079 (1) years (Table 1). Interestingly, the dust change reaches its minimum value on average about 6 years after the δ 18 O. These abrupt changes also correspond to abrupt temperature differences of on average 11.6 ± 2.7 (1σ ) • C for the GIs 17.1 to 2 (Kindler et al., 2014; Table 1 in this paper). The amplitude of the change in the dust concentration corresponds on average to a factor of 6 in about 57 years. These values are different from previous investigations of the NGRIP dust record, which were showing much more rapid changes in a few years, but they seem more realistic than extremely abrupt shifts that would be particularly difficult to interpret from a dynamical point of view. Isotope studies of the dust recorded in the Greenland ice cores demonstrated an Asian origin for both summer and winter seasons (Biscaye et al., 1997;Svensson et al., 2000), while the analysis of modern dust preserved in modern firn indicates a different origin for summer (essentially Taklimakan Desert) and winter (Taklimakan Desert and north-ern Chinese and Mongolian deserts) (Bory et al., 2002). The transport of the modern chemical aerosols towards Greenland also partly supports an Asian origin and therefore allows interpretation of the dust record in the ice cores (Goto-Azuma and Koerner, 2001). Indeed, the major dust event recorded in China on 6 April 2001, the largest ever recorded dust storm worldwide, provides relevant information about the Asian dust transport, past and present, towards Northern Hemisphere ice sheets. The identification of the dust particles related to this event on top of Mount Logan, Alaska, indicates that the assumption of Asian dust origin in the past is also particularly adequate (Zdanowicz et al., 2006). The presence of the large Laurentide ice sheet over North America at least at the last glacial maximum caused a split of the polar jet stream, inducing two main pathways for the Asian fine dust transport eastwards (Fig. 4). A comparison of the different GI dust records indicates that the dust concentration falls by a factor of 8 in about 60 years from the beginning of the warming during DO events. These values are similar to the variations of a factor of 5 to 7 during 4 decades observed during the 15.5-11.5 ka b2k interval (Steffensen et al., 2008).  Reducing so drastically the dust concentration in the ice core implies a similar drastic reduction of the dust concentration in the atmosphere, relying on changes in the sources, in the transport, and/or in the deposition of the particles (Fischer et al., 2007).
Concerning the changes in the sources, several points can be taken into consideration. At present the amount of dust emitted from different Chinese deserts into the atmosphere occurs mostly in April, about 43.4 % of the 1451 Mt of dust emitted from 1996 to 2001 (Laurent et al., 2006). However, these sources do not behave similarly. While the average quantity of dust annually emitted is similar in the Gobi and Taklimakan deserts, the frequency of the dust storms is different, being more numerous in the Taklimakan than in the Gobi (Laurent et al., 2005). A reduction in the size of these sources does not seem to be a reliable cause, as these deserts did not vary much in the relevant past (Rittner et al., 2016). Therefore, a reduction in the availability of dust material should be sought by considering atmospheric circulation changes. Studying the electron spin resonance (ESR) intensity signal from marine core MD01-2407 from the Japan Sea, Nagashima et al. (2011) have shown that the Mongolian Gobi source was dominant during the GSs, while the Taklimakan Desert source was dominant during the GIs, with some variation between the interstadials. The analysis of Chinese speleothems from caves located southward of the Chinese Loess Plateau (CLP) indicates millennial-scale variations in the stalagmite growth, related to variations in pre-   cipitation over the cave area linked to stronger summer monsoon activity, which are synchronous to the GIs (Wang et al., 2001(Wang et al., , 2008. From two key Chinese loess sequences, Sun et al. (2012) have described abrupt millennial-scale changes expressed by detailed synchronous mean grain size variations over the 60-10 ka b2k interval in deposited wind blown material, from coarser material during the GSs to thinner particles during the GIs. The sedimentation rate estimated from the key sequences of Jingyuan and Gulang varies between 6 and 63 cm kyr −1 for Jingyuan and between 10 and 130 cm kyr −1 for Gulang (Fig. 4). Therefore, as meridional circulation prevailed at least in eastern Asia due to the monsoonal system, aridity could have been reduced during the favorable season (around April, with some variation and delay, depending on the GSs or Heinrich stadials (HS) as observed in eastern Europe by Sima et al., 2013) of dust emission of these short intervals (Sun et al., 2012). A recent study of the decline of the snow cover over northeastern Russia, southwestern Asia, northern India, and the Tibetan Plateau reports that the snow cover decline creates favorable dynamical conditions for strong winds over the Arabian Sea and thus favoring a stronger summer monsoon (Goes et al., 2005). Complementarily, a modeling experiment indicates that a strong East Asian monsoon could be related to planetary waves related to the extension of the Northern Hemisphere ice sheets, affecting the precipitation band at the latitude of the CLP (Yin et al., 2008). However, as the main interval for dust emission is April, one can hardly assume that the summer monsoon is the main driver of a reduction in dust emissions. On the contrary, modern East Asian dust storms are related to surges of cold air originating from Siberia, and such conditions should have been reduced during the GI compared to the GSs at least during April, the main dust emission season. Reduced cold outbreaks could have resulted from the reduction of the snow cover due to warming Eurasia and to a negative feedback related to dust deposition on snow as deduced from another modeling experiment (Krinner et al., 2006). Europe has been strongly impacted by these North Atlantic millennial-scale climate changes, as observed in different types of deposits (Genty et al., 2003;Müller et al., 2003;Rousseau et al., 2002Rousseau et al., , 2007b (Fig. 1). The influence of the westerlies and the position of the polar jet stream were constrained by the variation in the extension of the sea ice during the last glacial interval (Sima et al., 2009(Sima et al., , 2013. Furthermore, the presence of ice sheets over Great Britain and Scandinavia, and an ice cap over the Alps enhanced the zonal circulation that reflects the location of the thickest loess deposits (Fig. 1b). Extensive investigation of European loess series along a west-to-east transect at 50 • N (Fig. 1b) reveals that the millennial-scale climate variations observed in the North Atlantic are preserved in these particular eolian deposits (2009; Antoine et al., 2001;Rousseau et al., 2002Rousseau et al., , 2007b. They clearly show alternating loess and paleosol units, which are continental equivalents, respectively, of the GS and GIs (Moine et al., 2017;Rousseau et al., 2017) ( Figs. 2, 3). The Nussloch loess sequence along the Rhine Valley is the most detailed record for the interval of 50 to 15 ka b2k. The nature of the observed paleosols is related to the duration of the corresponding GIs themselves : GI8, the longest of the last eight GIs, is represented in Nussloch by a brown arctic soil, while the youngest GIs correspond to tundra gleys or embryonic soils (oxidized horizons) for GI3 and GI2, which are particularly short Rousseau et al., 2002Rousseau et al., , 2007b (Figs. 2, 3). An estimate of the duration of the continental equivalent of the GIs can be deduced from the Greenland dust record by considering the interval after the abrupt warming, when the dust concentration was at a minimum in the atmosphere and therefore shows low values in the Greenland ice cores. Such succession is observed over an area more than 2000 km wide from western Europe eastward to Ukraine during the last climatic cycle (Rousseau et al., 2011), showing the influence of the zonal circulation, as it is the case in modern time (Fig. 5a, b). In the Nussloch key sequence, the sedimentation rate of dust particles varies between 23 and 157 cm kyr −1 for the loess units synchronous to the GSs. These values are similar to those observed at the northern edge of the CLP by Sun et al. (2012) (Table 2, Figs. 1a, 4), supporting an apparently similar eolian characteristic within more global dynamics. Furthermore, modeling experiments show that dust emission mainly occurred in April (Sima et al., 2009) similar to modern deserts, although the European sources are not deserts at all presently. These modeling studies also indicate that the zonal circulation not only prevailed during the GI intervals but also during the GSs as recorded by loess deposits. For the coarsest material these eolian units are composed of material of local and regional (up to about 500 km) origin, mainly from dried river beds . The finest grains originate from longer distances, but still from deflation areas, mostly located in the emerged English Channel, North Sea, or northern European plain, and at the margin of the Fennoscandian frontal moraines (Sima et al., 2009;Rousseau et al., 2014) (Fig. 1b). Since we observe in Kurtak,Siberia (Figs. 1a,4) synchronous alternations between eolian loess units and paleosols similar to those observed over Europe (Haesaerts et al., 2005), we infer that a zonal circulation similar to the present one (Fig. 5) prevailed over Eurasia during both the GSs and GIs, synchronizing the loess-paleosol sequences between Europe and Asia (Haesaerts et al., 2005). The initial warming over Greenland and the North Atlantic is transported eastward by the prevailing zonal circulation. As explained above, this propagating warming should have reduced the production of cold surges contributing to the strong atmospheric dynamics (Sun et al., 2012) responsible for the dust emission in northern Chinese deserts, which are the main dust suppliers of both the CLP and Greenland ice sheet. The stronger winds during the GSs and weaker winds during the GIs permitted the record of DO-like intervals in Chinese loess sequences through grain size variations with coarser material deposition during the GSs and finer during the GIs (Sun et al., 2012).
All the mechanisms described above involve the Northern Hemisphere and do not address the origin of these abrupt changes. Recent investigations in West Antarctic Ice Sheet Divide (WAIS) ice core (Buizert et al., 2015) have demonstrated the north-to-south-directed transfer of the abrupt changes with Greenland warmings leading Antarctic coolings; the associated heat transfer is argued to be modulated by the Atlantic Meridional Overturning Circulation rather than by the atmosphere (Henry et al., 2016). Although the origin of the DO events is not yet elucidated, the North Atlantic changes still remain the drivers of those observed over Eurasia. Indeed, the modern cyclonic pattern of moisture transport towards the Greenland ice sheet  depends on the sea surface conditions over the North Atlantic Ocean and clear pathways of this transport have been identified (Chen et al., 1997). During the last climate cycle, this pattern was related to both the sea ice extent and the size and expansion of the Laurentide ice sheet in North America, which impacts the Northern Hemisphere atmospheric circulation (Kutzbach, 1987). The large changes in deuterium excess at the onset of the warming reveal large variations in the moisture source regions. These variations are most likely associated with changes in sea ice extent, which shift the source regions, combined with changes in atmospheric circulation patterns Sodemann and Zubler, 2010;Sodemann et al., 2008). Due to the changing albedo, sea ice extent itself directly impacts the atmospheric circulation. Conversely, sea ice is mainly driven by the meridional overturning oscillation, which would have thus induced the observed δ 18 O and hence temperature changes noticed in the Greenland ice cores (Henry et al., 2016). In turn, this would have also impacted the Northern Hemisphere atmospheric circulation, contributing to the abrupt millennial-scale changes, including the DO events in Eurasian terrestrial records. At middle latitude these changes constrained the emission of local and regional dust particles, their transport, and their deposition in the loess sequences, which can be characterized by similar sedimentation rates over the area corresponding to a general climate dynamics, but nevertheless different mass accumulation rates (MARs) related to the bedrock of the source areas of the deposited material. In a previous study, Kohfeld and Harrison (2003) estimated MARs over CLP varying between 21 and 809 g m −2 yr −1 for MIS 3 and between 60 and 5238 g m −2 yr −1 for MIS 2. In Europe, Frechen et al. (2003), applying the same method, indicate MARs varying between 100 and 7000 g m −2 yr −1 , with particularly high values for Nussloch (1213-6129 g m −2 yr −1 ). These estimates are higher than those determined in our study after reevaluation of the chronology of the key sequence, still following the same calculation, which yields varying MAR values between 376 and 2586 g m −2 yr −1 (376-1952 g m −2 yr −1 for MIS3 and 724-2586 g m −2 yr −1 for MIS2) or between 395 and 2515 g m −2 yr −1 (395-1952 g m −2 yr −1 for MIS3 and 723-2515 g m −2 yr −1 for MIS2) when considering the δ 18 Oand dust-related NGRIP chronologies, respectively (Table 2). Still, as Nussloch represents an exceptional record of MIS 3 and 2, we consider our results as corresponding to the highest boundary values for this time interval in European deposits.

Conclusion
Our study thus shows that a strong emphasis has to be placed on the past dust record, which is still poorly understood and weakly integrated into general circulation models. The important uncertainties associated with mineral aerosols in the recent IPCC report (Solomon et al., 2007) are a clear indication that more must be done on this particular parameter. Because this is a key factor in the climate system, understanding the past dust cycle is an important requirement, especially when estimating the dust load in the past atmospheres, and this should open new fields of investigation for a better constraint of the climate variability in different contexts.
The present study provides new insights into the analysis of the millennial-scale variability and of abrupt climate changes by proposing the links by gathering atmospheric, marine, and continental records. It shows that the complete understanding of the whole climate system requires investigations at high resolution in every domain with the support of modeling experiments. In the present study, we show that both δ 18 O increase and dust decrease take place over an interval of about 50 years on average from the start of the abrupt change. This corresponds to the 4 decades previously mentioned for the two warmings occurring during the 15.5-11.5 ka b2k interval, which are associated with strong resumptions of the meridional overturning circulation (Mc-Manus et al., 2004). This makes the potential change in the atmospheric dynamics more reliable. Our investigation provided an explanation of the record of the abrupt climate changes in the Northern Hemisphere dust records, both in the Greenland ice sheet and over Eurasia. It shows that dust should be considered one of the major players in the past millennial climate variability. Data availability. The data sets used in this study are available at http://www.icecores.dk.
Author contributions. DDR designed the research. DDR, SJJ, AS, MB, AS, and JPS performed the research. SJJ, AS, MB, and JPS performed drilling and analysis of NGRIP ice cores. AS performed the modeling experiment. DDR designed and performed the loess sequences investigation. DDR performed the new Nussloch chronological sequence and calculated the associated sedimentation rates and MARs. NB performed the algorithmic determination of the DO events in both δ 18 O and dust NGRIP records. DDR, NB, and AS wrote the paper.