Coupled European and Greenland last glacial dust activity driven by North Atlantic climate.

Centennial-scale mineral dust peaks in last glacial Greenland ice cores match the timing of lowest Greenland temperatures, yet little is known of equivalent changes in dust-emitting regions, limiting our understanding of dust-climate interaction. Here, we present the most detailed and precise age model for European loess dust deposits to date, based on 125 accelerator mass spectrometry 14C ages from Dunaszekcső, Hungary. The record shows that variations in glacial dust deposition variability on centennial-millennial timescales in east central Europe and Greenland were synchronous within uncertainty. We suggest that precipitation and atmospheric circulation changes were likely the major influences on European glacial dust activity and propose that European dust emissions were modulated by dominant phases of the North Atlantic Oscillation, which had a major influence on vegetation and local climate of European dust source regions.

Atmospheric mineral dust (hereafter called dust) plays a major but poorly understood role in the climate system, both responding to and driving climate change (1). Dust concentrations in the Greenland ice cores show a close relationship to Greenland temperature at millennial to centennial timescales, although the precise causes and effects of this relationship remain unclear (2). At millennial to centennial timescale, last glacial climate variability over Greenland, and likely over the North Atlantic generally, is dominated by the Dansgaard−Oeschger (D-O) cycles (3, 4). In Greenland, each D-O cycle over the last 120 ka is characterized by (i) gradual cooling over a century to many millennia; (ii) subsequent, more rapid decline into peak stadial conditions; and (iii) very abrupt return to interstadial conditions, varying in amplitude from 5 K to 16.5 K (5, 6). Several mechanisms behind these abrupt climatic changes have been proposed, including ocean circulation (7), sea ice cover–ice sheet (8, 9), and wind pattern changes (10). The effects of D-O events extended across much of the Northern Hemisphere (11) and were identified in precisely dated terrestrial records such as speleothems in Europe and Asia (12–15).

However, changes in climate and dustiness in potential dust-emitting source regions are less well understood at millennial scales (16). Furthermore, recent isotopic and clay mineral evidence from Northern Hemisphere loess deposits reopened the debate about the origins of Greenland dust, as both East Asian and European dust sources are compatible with current Greenland dust data (17). Modeling studies highlight that changes in precipitation patterns and vegetation cover over these potential source areas profoundly affected the dust cycle in response to D-O–Heinrich events (18, 19), but well-dated, high-resolution dust records covering parts of or the whole last glacial period are scarce. Some studies on extensive windblown terrestrial dust deposits of loess that cover large areas of land surface have suggested imprints of the D-O events exist in loess grain size records in Europe (20, 21) and Asia (22). However, the chronologies of these archives are often poor or not fully independent. Furthermore, the vast majority of fully independently dated loess sequences are dated using luminescence methods (23, 24). Although these methods generally provide accurate age estimates for Late Quaternary deposits, they do not have sufficient precision to investigate the timing and underlying mechanisms of millennial–centennial-scale environmental changes, including dust accumulation variability on the continents (25, 26).

Thus, despite being globally widespread and closely tied to or acting as major global dust source regions, loess records have so far remained a largely untapped source of information on dustiness at millennial timescales and below. As other dust records tend to be geographically isolated and restricted in area, this has been a major impediment to understanding cause and effect in climate−dustiness feedbacks. Recent advances in high-precision accelerator mass spectrometry (AMS) radiocarbon dating of small gastropods in loess (26–28) provide a new opportunity to gain unprecedented insights into the temporal variations of the terrestrial dust cycle at millennial to centennial scale from loess deposits. A recent radiocarbon dating study of earthworm granules from soils in a loess site in Germany points to a link between Greenland climate and soil formation in western European loess (29). Although this highlights the potential of radiocarbon dating of loess, it is not clear whether the inferred precipitation changes in that study are accompanied by shifts in dust accumulation, crucial to understanding past dustiness. Here we present a high-resolution loess grain size and dust accumulation record from southern Hungary covering the period of 36 to 24 ka before 2000 (b2k) with a fully independent radiocarbon chronology based on a uniquely detailed age dataset consisting of 125 AMS 14C dates and an associated Bayesian age−depth model. Our time resolution allows investigating the timing of abrupt European climate shifts compared with the North Greenland Ice Core Project (NGRIP) Ca2+ record in Greenland (5). We also measured species-specific mollusc shell carbon and oxygen isotopic data as potential indicators of coupled precipitation–vegetation changes, which may be a major driver of continental dust production and deposition variations in east central Europe at centennial–millennial timescales.

Results

Site and Chronology.

The Dunaszekcső loess−paleosol sequence is situated in southern Hungary (Dsz, Fig. 1). This part of the Carpathian Basin is an area of low relief between the main mountain ranges of central Europe and is under Atlantic, Mediterranean, and continental climatic influence. This is expressed in the amount of annual rainfall (575 mm⋅y−1, with extremes of 276 mm⋅y−1 to 882 mm⋅y−1) and mean air temperatures during winter (∼0.3 °C, January) and summer (∼22.2 °C, July) as measured at a nearby meteorological station for the period 1998–2013. Back-trajectory analyses revealed five major sources of present-day precipitation at K-puszta, a site on the Great Hungarian Plain: Mediterranean region (57.0%), local moisture (14.8%), Atlantic region (14.2%), north Europe (7.4%), and east Europe (6.6%) (30).

Fig. 1.

Map showing the location of the study sites and paleoclimate archives mentioned in Results and Discussion. C, Crvenka in Serbia; Dsz, Dunaszekcső in southern Hungary (this study); 7H, Sieben Hängste speleothem; N, Nussloch in Germany; NGRIP, NGRIP ice core; S, Surduk in Serbia.

The exposed section of the Dsz record (ca. 17 m) is the upper part of a ∼70-m-thick Quaternary aeolian sediment sequence and covers the last glacial–interglacial cycle (26). Sediment samples were collected at 5-cm resolution between 250 and 1,045 cm in the profile for dating and proxy analyses (see ). As shown by calibrated AMS 14C ages of numerous charcoal fragments and mollusk shells (Datasets S1 and S2), the sampled profile covers a period of loess sedimentation between 36 cal y BP and 23.4 cal y BP. Radiocarbon ages of charcoals and small mollusks are consistent despite the very different origin and genesis of these materials (26, 28), and demonstrate that the age dataset is reliable and robust. The Bayesian age−depth model of the composite Dsz record is calculated at 1-cm intervals and has a minimum 95% age uncertainty of 188.7 y at 488 cm (25.796 cal y BP) and a maximum of 1,328.7 y at 1,045 cm (35.988 cal y BP) (Figs. S1–S3). In the older section, charcoal 14C ages become less precise as they approach the practical limit of radiocarbon dating.

Sedimentation–Mass Accumulation Rates and Grain Size.

Sedimentation rates have changed significantly over time, as indicated by relatively low values (mean: ∼0.4) of the posterior distribution for memory (or autocorrelation) after Bayesian age modeling using Bacon software (31) (Figs. S2 and S3). We therefore infer that the environmental conditions or internal dynamics of the aeolian system that influence loess accumulation were highly variable at the study site. Bulk mass accumulation rates (BMARs) and dust mass accumulation rates (DMARs; <10 μm dust fraction) range between 380 g⋅m−2⋅y−2 and 2,885 g⋅m−2⋅y−2 (1-cm-resolution model) and 69 g⋅m−2⋅y−2 and 629 g⋅m−2⋅y−2 (5-cm-resolution model), respectively, showing an almost 10-fold increase between the minimum at 32.6 ka b2k and the maximum at 26.1 ka b2k (1-cm-resolution model, Dataset S2). BMAR and DMAR display a longer-term increasing trend, following insolation at 45°N, and exhibit centennial- to millennial-scale variability corresponding, in general, to the stadial–interstadial pattern recorded in Greenland (Fig. 2). Significant decreases and increases in BMAR occur within short periods of time, generally coeval with the NGRIP dust record (Figs. 2 and 3).

Fig. 2.

Comparison of insolation and paleoclimate records for the 37- to 22-ka interval. (A) Insolation for June 21 at 45°N (32), and (B) integrated insolation for the spring season at 45°N (33). (C) Median grain size of quartz grains (d50quartz) in the Dsz loess record (5-cm resolution). (D) Median grain size of bulk loess (d50bulk) in the Dsz loess record (5-cm resolution). (E) DMARs of loess at Dsz calculated from the lower (5 cm) resolution age model. (F) BMARs at Dsz calculated from the higher (1 cm) resolution Bayesian age–depth model. (G) Sieben Hängste (7H) composite stalagmite δ18O record (Western Alps) (15). (H and I) NGRIP dust and ice/water δ18O records (5). Yellow bars denote GI periods as given by Rasmussen et al. (5), while gray bars indicate phases of reduced Ca2+ (less dust) within GS periods.

Fig. 3.

Dust accumulation variations in Greenland and east central Europe and precipitation patterns in the Western Alps for the 32- to 24-ka period. (Lower) Records shown are (A) high-resolution BMAR in southern Hungary (Dunaszekcső record), (B) the 7H stalagmite δ18O record (15), and (C) NGRIP dust (Ca2+) (5). (Upper) The 2-sigma age model (or counting) uncertainties of the displayed records. Yellow and gray bars indicate decrease in dust accumulation over Greenland. Dust minima identified as interstadials in the NGRIP record are marked by yellow bars (GI-5.1, GI-4, and GI-3; ref. 5).

Median particle diameters of bulk loess (d50bulk) and quartz in loess (d50quartz) vary between 18.5 μm to 55.7 μm and 22.4 μm to 64.3 μm, respectively, and both exhibit centennial-scale variability. While quartz grain size does not show an insolation-driven trend, bulk grain size displays a general coarsening after ∼30 ka b2k to 31 ka b2k. A striking feature of the grain size records is that they are not characterized by consistent coarsening–fining trends that match Greenland Stadial–Greenland Interstadial (GS/GI) patterns, and the two records often show opposite trends to the Dsz BMAR and DMAR, and the NGRIP Ca2+ records (e.g., GI-5.1 and GI-4, Fig. 2).

Discussion

Comparison of dust accumulation in Greenland and Hungary reveals some striking patterns on orbital, millennial, and centennial timescales, which we suggest demonstrates close and potentially causative links between North Atlantic climate and European dust activity. Over orbital timescales, there is a pattern of increasing BMAR in the Dsz record after 31 ka b2k (Fig. 2), which corresponds to decreasing spring–summer insolation at 45°N latitude, culminating in the highest dust inputs during the early Last Glacial Maximum (LGM). This pattern is broadly consistent with the Greenland Ca2+ dust record, although the latter shows less of a marked increase at the LGM and this increase seems to last longer (5). The LGM was an especially cold and dry period in the North Atlantic and Europe (34, 35), accompanied by changes in European vegetation cover (36) and different atmospheric circulation patterns driven by profound changes in insolation (32), topography (ice sheets, ref. 37) and trace gas compositions (34). The orbital timescale changes in the dust cycle, as recorded in the Dsz sequence, were likely driven by insolation changes and longer term cooling–drying trends that resulted in surface conditions appropriate for more intense dust emissions in Europe during the LGM.

Beyond this first-order insolation control, centennial–millennial timescale variability in Dsz BMAR and DMAR also became more pronounced from 31 ka b2k onward. Some traces of these higher frequency oscillations can be observed before 31 ka, although the precise timing becomes more uncertain because of the lower number of radiocarbon ages from this part of the sequence and the fact that the age model is less precise at this point. In general, the high-precision 14C chronology and BMAR record at Dunaszekcső allows unprecedented insights into the timing of abrupt European dust cycles changes and a detailed comparison with Greenland climate and dust. A striking feature is that Dsz dust accumulation minima generally match GI phases (GI-5.1, GI-4, and GI-3; Fig. 3) in the NGRIP δ18O/Ca2+ records (5) within dating uncertainties, and also with shorter periods of low Ca2+ during the LGM (gray bars, Figs. 2 and 3). While the timing of Dsz BMAR decrease and subsequent increase nicely matches the onset and end of GI-5.1 in the NGRIP record, the onset of GI-4 is missing in the Dsz record. This is because BMAR stays relatively low (at interstadial levels) after a doubling from ∼600 g⋅m−2⋅y−1 to ∼1,000 g⋅m−2⋅y−1 at ∼30.6 ka b2k to 30.5 ka b2k. The DMAR minimum around 29 ka b2k as well as the subsequent abrupt increase in BMAR and DMAR are within dating uncertainty of their GI-4 counterparts in the Greenland record (Figs. 2 and 3). Similarly, a DMAR and BMAR minimum is synchronous with the Greenland GI-3 within dating uncertainty. Significant increases in dust sedimentation occur during GS phases, with BMAR peaks reaching and exceeding 2,500 g⋅m−2⋅y−2 between ∼25.8 ka b2k and 26 ka b2k. This maximum matches the timing of peak dust accumulation over central Greenland (Figs. 2 and H and 3), although the highest dust accumulations in Dsz are split into two distinct peaks. By contrast, grain size variations (d50bulk and d50quartz) do not necessarily coincide with BMAR, DMAR, and the GI/GS variability, and tend to show coarsening during GI periods (e.g., GI-5.1 and GI-4; Fig. 2 ), when the opposite would be expected. Mismatches between grain size and dust accumulation rates (AR) were also found at Crvenka in Serbia (23) and at sites on the Chinese Loess Plateau, using luminescence dating (38). This result contrasts with previous observations at Nussloch, Germany, where layers of coarse (fine) loess grain size have been suggested to correspond to GS (GI) phases in Greenland (20, 21). This may reflect the differences in aeolian sedimentation between a more humid, western European site, where grain size is largely controlled by soil formation, and those in drier east central–SE Europe where grain size can be directly linked with aeolian transport, trapping, and deposition. The general disagreement between the BMAR, DMAR, and d50bulk, quartz proxies is likely due to the fact that grain size is a much more complex parameter on short timescales. Grain size distributions (GSDs) are affected by multiple, often stochastic processes, and the clay- to sand-sized particles making up loess are mobilized, transported, and deposited in many different ways (39). Conversely, BMAR only reflects the amount of particles deposited per unit time and area, which is a function of atmospheric loading, local trapping, and preservation conditions. In this study, we will therefore refrain from using grain size to identify millennial- to centennial-scale events in loess deposits.

The match between the Dunaszekcső and Greenland dust records can be further analyzed within a wider European context. The onset of GI-5.1 in Greenland is synchronous with the reduction in dust accumulation in southern Hungary within the age scale uncertainties, and there is an excellent match between GI-5.1 and the dust minimum at ∼30.7 ka b2k (Fig. 3). While BMAR stays relatively low (still in the range of 600 g⋅m−2⋅y−1 to 1,000 g⋅m−2⋅y−1) after GI-5.1, a significant increase in dust input starts at 28.8 ka b2k, coincident with a drop in δ18O of the 7H Alpine stalagmite composite record, and the end of GI-4 within dating uncertainty (Fig. 3; note that 2σ error at 28.8 ka is 574 y for Dunaszekcső). The onset of GI-3 is coeval in Greenland and the Dsz record, and the duration of the dust flux decrease in the Dsz record is longer than that of GI-3 (240 y) (5). At the same time, the 7H record exhibits the shortest D-O event at ∼27.7 ka b2k. Further smaller decreases in dust accumulation at Dsz occurred in line with similar events in Greenland, which can also be found in the 7H record within dating uncertainty (Fig. 3, gray bars). Unfortunately, the combined uncertainty of the Dsz age model together with the IntCal13 calibration uncertainties (IntCal13 is the latest calibration curve published by the IntCal group of the international radiocarbon community, ref. 40) preclude the identification of any lead or lag between the various records. Nevertheless, the observation that there is a lot of shared variability between NGRIP and Dsz over 32 ka b2k to 24 ka b2k suggests a good general agreement between the ice core and IntCal13 chronologies and may potentially indicate a possible European dust contribution to Greenland ice cores, an issue that has not yet been satisfactorily resolved from provenance studies and has important implications for dust transport pathways (17).

The match between the Dsz BMAR and 7H δ18O records indicates that large-scale atmospheric reorganizations are the likely cause of dust cycle changes in east central Europe. We here propose that the predominant modes of the North Atlantic Oscillation (NAO) during the LGM and resultant changes in storm tracks may explain these patterns. The 7H speleothem δ18O dataset reflects precipitation patterns and reveals intense southern moisture advection from the Mediterranean toward the Alps during the LGM, related to Rossby wave breaking of the jet stream over western Europe (15). A related southward shift in the position of the North Atlantic storm track was also found in recent modeling studies of LGM atmospheric circulation (37, 41, 42), explained by a tendency toward more cyclonic wave breaking (CWB) and less anticyclonic wave breaking (AWB) events in the LGM runs (41, 42). Importantly, under current climate, teleconnections such as the NAO (characterized by a meridional displacement of the upper-tropospheric Atlantic jet; ref. 43) have a dynamic link to Rossby wave breaking events (44). The negative (positive) phase of the present-day NAO relates to a more southward (northward) position of the eddy-driven Atlantic jet, and exhibits more CWB (AWB) than the long-term average (42). The effect of this dynamic link on climate is stark, as it affects seasonal and, in particular, wintertime temperature and precipitation patterns in Europe. For example, during the prolonged Maunder Minimum cold period from 1640 to 1715, associated with an NAO‒ phase, winters were characterized by more-frequent blocking situations connected with cold air outbreaks toward central and eastern Europe (45). Springs were cold and characterized by a southward shift of the midlatitude storm tracks, and summers in central Europe were wetter and slightly cooler than they are today. Simulations performed using a coupled atmosphere–ocean model for the Maunder Minimum revealed higher cyclone density in central–east central Europe for both the winter and spring periods (46), which were also found to be the major dust emission periods in east central Europe under cold climate states during the last glaciation (February–June; ref. 18). In contrast to the Maunder Minimum cold period, a persistent NAO+ mode was reconstructed for the Medieval Warm Period from ∼800 to 1300 (47). Modeling studies suggest that the NAO during the LGM differed from the modern one, as it was characterized by four centers of action (48), and weaker latitudinal fluctuation of the eddy-driven jet (42). However, some simulations demonstrate that latitudinal wobbling remains the primary type of North Atlantic jet variability and also dominates jet pulsing during glacial times (49). As such, we suggest that the dust minima in the Dsz BMAR and DMAR records during Greenland GI events could be explained by a prolonged NAO+ mode with a more northerly jet position (NW trajectory in the 7H composite record, Fig. 3). Conversely, dust maxima during GS events would be linked to a persistent NAO‒ phase with a more southerly storm track (S trajectory; 7H δ18O, Fig. 3), due also to nontopographic forcing such as centennial-scale sea ice or sea surface temperature changes in the North Atlantic (49). These conditions, when both winter and spring were colder and characterized by higher cyclone density (45, 46), may have favored enhanced European dust emission and deposition.

Mollusk shell stable isotopic composition measurements obtained from Dsz suggest shifts from dominantly C3 vegetation and wet summers during GS periods to mixed C3/C4 vegetation and dry summers during GI phases (see and Fig. S6). Such changes strengthen our conceptual model: During periods of dominant NAO‒ phase, eastern European summers would have been cooler, with available water for longer periods enabling C3 plant dominance, but winters and springs would have been dry and stormy, which would have increased dustiness in dust-emitting regions such as floodplains and outwash fans. In contrast, during periods of dominant NAO+ phase, eastern European summers would have been drier, benefiting C4 vegetation. NAO‒ and NAO+ phase dominant periods coincide with Greenland GS and GI, respectively.

Methods

Site and Sampling.

Sampling for grain size analyses was performed at 5-cm resolution on two overlapping profiles (Fig. S1) of a carefully cleaned loess–paleosol sequence at Dunaszekcső (southern Hungary, 46°05′25″N, 18°45′45″E, 135 m above sea level), freshly exposed by a landslide in 2008. AMS radiocarbon dating was conducted on samples collected at 5-cm-depth resolution between the depths of 485 and 850 cm and at 15- to 30-cm-depth resolution for the rest of the sequence. For 14C-dating, loess cuboids with dimensions of 15 × 5 × 10 cm (width–height–length) were prepared and cut from the sediment. Sample blocks were subsequently disintegrated in the laboratory by soaking them in distilled water. Charcoal fragments and mollusk shells were extracted by washing the sediments through a 1-mm-mesh sieve, then dried at 50 °C and handpicked using gloves and precleaned forceps to avoid modern carbon contamination. After being identified at the species (or family) level, shells were wrapped in Al foil and put in closed plastic bags. The charcoal fragments were handled and packed in a similar way, but separately from mollusk shells. The general validity of the AMS 14C-based chronology was confirmed by optically and infrared-stimulated luminescence (OSL-IRSL) dating (26).

AMS Radiocarbon Dating.

Before the 14C measurements, charcoal fragments were treated using the standard acid–base–acid (ABA) method. After the final acid wash, the samples were washed again with distilled water to neutral pH and freeze-dried overnight. First, dried charcoal fragments were combusted in an on-line combustion system using CuO in one step at 1,000 °C (ABA-OSC1000), while, during subsequent runs, they were subjected to stepped combustion in a pure O2 gas atmosphere, first at 400 °C and then at 800 °C (ABA-TSC400 and ABA-TSC800) (28). Two samples with low carbon yields were measured via a gas ion source attached to the AMS.

Mollusk shells were ultrasonically washed and etched using weak acid (2% HCl) to remove surface contaminations (20 to 30% of the original mass removed). Subsequently, acid-cleaned shells were dried and put into vacuum-tight two-finger reaction ampoules (∼100 cm3 inner volume) and dissolved using phosphoric acid. CO2 was produced by acid hydrolysis of shells at 75 °C, further purified cryogenically, and then graphitized (50).

All of the 14C measurements were undertaken using a compact radiocarbon AMS system (MICADAS) at the Hertelendi Laboratory of Environmental Studies, Institute for Nuclear Research, Debrecen (Hungary). Conventional radiocarbon ages were converted to calendar ages using OxCal online (version 4.2; ref. 51) and the IntCal13 calibration curve (40).

Composite Profile.

The composite profile was obtained using the AMS 14C ages and bulk median diameter (d50bulk) values measured in profiles 1 and 2 (Fig. S1). The validity of this approach was subsequently checked by Monte Carlo simulation with the “Intra-site correlation age modelling” code (52).

Age−Depth Modeling.

Bayesian age−depth modeling was performed using the Bacon code (31), based on 125 radiocarbon data points (Dataset S1). Inverse AR (sedimentation times, years per centimeter) were estimated from 35.6 million Markov Chain Monte Carlo (MCMC) iterations, and these rates form the age–depth model. Inverse AR were first constrained by default prior information: accumulation (acc.)shape = 1.5 and acc.mean = 20 for the gamma distribution, and memory (mem.)mean = 0.7 and mem.strength = 4 for the beta distribution describing the memory effects (or autocorrelation) of inverse AR. All input age data were provided as 14C y BP, and Bacon used the IntCal13 calibration curve to convert conventional radiocarbon ages to calendar ages. Age modeling was run to achieve a 5-cm final resolution. In a second attempt, as part of a sensitivity test, the parameters were set as acc.shape = 2 and acc.mean = 20 for the gamma distribution and mem.mean = 0.4 and mem.strength = 17 based on the posterior gamma and beta distributions from the first modeling attempt (Fig. S2). Posterior gamma and beta distributions better fitted the priors for the second model (Fig. S3), and the mean of inverse AR after modeling was found to be ∼25 y⋅cm−1, matching very well the regional average (26.55 y⋅cm−1) of east central European loess records (Table S1). Although only insignificant differences were found between these age models (Figs. S4 and S5), the second age–depth model was used to interpret proxy data in this study (Dataset S2). Finally, the age modeling was run again using the second parameter set to obtain a 1-cm-resolution age–depth model (Dataset S2). To make the age scale comparable to the ice core timescale (GICC05), 50 y were added. All ages are therefore presented in figures in ice core years b2k (before year AD 2000).

Mass Accumulation Rate Calculations.

AR were calculated using the age–depth model as AR (meters per year) = d2 − d1/a2 − a1, where d1 and a1 are consecutive depths and weighted mean model ages in the profile. BMARs have been computed using the equation BMAR (grams per square meter per year) = AR × ρ × f, where ρ is the dry bulk density (grams per cubic meter), and f is the sediment fraction that is aeolian in origin. Since loess is entirely aeolian, then f = 1 for all of the calculations. A value of 1.5 g⋅cm−3 was used for dry bulk density as the best estimate for Hungarian loess (53). For the purpose of model–paleodata comparison, DMARs were calculated as DMAR (grams per square meter per year) = AR × ρ × f, where f is the <10-μm fraction of the total/bulk aeolian mass (16, 53). This fraction was determined from the laser particle size measurements. DMAR is considered a minimum estimate of dust deposition, as DMAR may be lower than the real dust deposition flux due to occasional erosion events (16).

Grain Size Analyses.

Before laser diffraction measurements, 3 g of loess samples were pretreated with 10 mL of 20% H2O2 and 10 mL of 10% HCl to remove organic matter and carbonates. Subsequently, 10 mL Na(PO3)6 was added to the samples, which were finally ultrasonicated for ∼1 min. These are chemically fully dispersed samples and the presented bulk loess grain size data are based on fully dispersed GSDs.

To obtain pure quartz separates, air-dried bulk sediment samples (2 g each) were treated with 30% H2O2 (10 mL/sample) for 24 h to remove organic matter. Subsequently, 50 mL of 6 M HCl was added, and the solution was boiled at 90 °C to 100 °C for 1 h to remove carbonates and iron oxides. Quartz was isolated by the sodium pyrosulfate fusion–hydrofluorosilicic acid method. This procedure resulted in pure quartz separates without affecting the grain size, grain shape, and surface textures of the quartz crystals, as has been proven by SEM imaging (39).

Grain size of bulk loess samples and quartz separates was analyzed using a Malvern Instruments Mastersizer 3000 laser diffractometer with a Hydro LV wet dispersion unit having a measurement range of 0.01 μm to 2,100 μm divided into 100 size bins. Two light sources were utilized, a red He–Ne laser at a wavelength of 0.633 mm and a blue LED at 0.470 mm. Diffracted light intensity was measured by 50 sensors over a wide range of angles. Constants of 1.33 for the refractive index of water, 1.544 for the refractive index of solid phases (valid for quartz, and most clay minerals and feldspars), and an absorption index of 0.1 were applied for both bulk loess and quartz. Bulk grain size analyses reported in this paper are the average of seven successive laser diffraction runs (total of 70,000 snaps), while those of quartz are the mean of three repeat measurements. The recorded data were processed using the Malvern’s Mastersizer 3000 software (version 3.10), which transformed the scattered light data to particle size information based on the Mie Scattering Theory. Median grain size of bulk loess samples (d50bulk) and quartz separates (d50quartz) was calculated from the Mastersizer 3000 software outputs.