Introduction

Mineral dust in the atmosphere, one of the key components of the climate system, can be used as an indicator to examine climate change. It is transported by the Earth’s wind regime, affected by both the regional and large-scale climate system1. Dust in proximal regions is mainly carried by near-surface winds via creeping and saltation2,3, while long-range transport (LRT) from distal regions is largely driven by synoptic-scale atmospheric circulation3,4. The frequency and intensity of such LRT events may be modulated by millennial-scale climate oscillations5. Additionally, atmospheric dust directly absorbs and scatters radiation and indirectly modifies the properties of clouds, affecting Earth’s climate systems6,7. Dust radiative forcing also makes an important contribution to polar amplification8. Thus, investigating provenance, distribution, and transport of dust is crucial to understanding dust-climate linkages and how atmospheric conditions respond to climate change.

During the last glacial period, dust production and deposition were greatly enhanced9,10. Dust is recorded in many deposits, particularly loess, which can incorporate both LRT dust and locally derived material9. Loess is primarily composed of silt-sized particles (2–63 μm), but often incorporates a significant fraction finer than 10 μm, potentially reflecting background dust transported hundreds to thousands of kilometers by westerlies11,12. Accordingly, loess and loess-like deposits serve as valuable archives of past dust dynamics. Located across central and northern Siberia, Yakutia is one of the major regions where loess is distributed11,13,14, and there is a special type of loess-like deposits called Yedoma (Fig. 1; Supplementary Fig. 1). Yedoma, or ice complex, is a silt-dominated ice-rich syngenetic permafrost with massive ice wedges, which developed across Beringia during the Late Pleistocene11,15,16. Although Yedoma in Alaska (i.e., eastern Beringia) has been considered to be aeolian-processed deposits since ref. 17, the depositional processes of Yedoma in Siberia (i.e., western Beringia) remain debated13,15,16. Because Siberia lies at the interface between major mid-latitude dust sources and the Arctic, understanding the provenance and transport mechanisms of Yedoma source material is crucial for reconstructing past atmospheric dust circulation and its climatic implications.

Fig. 1: Map of the study area.

Yedoma deposits (yellow-filled areas; from refs. 65,[66](https://www.nature.com/articles/s41467-025-65772-2#ref-CR66 “Strauss, J. et al. Database of Ice-Rich Yedoma Permafrost Version 2 (IRYP v2). https://doi.org/10.1594/PANGAEA.940078

(2022).“)), possible source areas (PSA; red-shaded and blue-outlined), and schematic dust transport paths (orange arrows) during the Late Pleistocene are shown. Arrows represent episodic events likely more frequent during their respective periods and do not depict mean climatology. Red, blue, and purple stars indicate the sampling locations of the Batagay, Cyuie, and Churapcha ice wedges, respectively. Reconstructed ice-sheet extents (Eurasian Ice Sheet Complex, EISC; Greenland Ice Sheet, GIS; Laurentide ice sheet, LIS; Cordilleran Ice Sheet, CIS) are shown for Marine Isotope Stages (MIS) 3 and the Last Glacial Maximum (LGM; MIS 2) as derived from refs. 60,[61](https://www.nature.com/articles/s41467-025-65772-2#ref-CR61 “Batchelor, C. L., Krapp, M., Manica, A. & Murton, D. K. The configuration of Northern Hemisphere ice sheets through the Quaternary. https://doi.org/10.17605/OSF.IO/7JEN3

(2025).“); see Methods for details. Major PSAs include the Taklimakan Desert (TK), Qaidam Basin (QD), Badain Jaran Desert (BJ), and Tengger Desert (TG) along the northern margin of the Tibetan Plateau (NMTP); the Yellow River (YR), Hobq Desert (HQ), and Mu Us Desert (MU) in the Ordos Plateau (OP); and the Mongolian Gobi (MG). Additional PSAs include the Gurbantunggut Desert (GT), Onqin Daga sandy land (OD), Hunlun Buir sandy land (HB), and Horqin sandy land (HO) along the northern boundary of China (NBC), as well as the Chinese Loess Plateau (CLP), Eastern Qinling Mountains (EQLM), Alaskan-Yukon loess and dust, Greenland glacial dust, and Western (W), Central (C), and Eastern (E) European loess. For comparison with potential local sources, we also include suspended particulate matter (SPM) from the Lena and Yana Rivers, Laptev Sea sediments (LSS), northeastern (NE) Siberian loess, Aldan Shield lamproites, and mafic magmatic rocks from the eastern margin of the Siberian Craton (Kharaulakh and Suordakh events). Base map includes rivers from Natural Earth (public domain). National boundaries are based on the TM World Borders Dataset 0.3 (available at https://larmarange.github.io/prevR/reference/TMWorldBorders.html), licensed under CC BY-SA 3.0 (https://creativecommons.org/licenses/by-sa/3.0/).

Syngenetic ice wedges, widely distributed in northeastern Siberia, grew during permafrost formation and accommodated significant amounts of sediment18, recording paleoclimate and paleoenvironment conditions during permafrost development. A substantial portion of the Yedoma deposits has been lost during post-depositional permafrost degradation, changing the permafrost environment since the Pleistocene-Holocene transition19,20. Nevertheless, well-preserved syngenetic ice wedges without thaw evidence are considered to retain original materials and geochemical signatures, which were implemented during Yedoma formation21,22,23. To elucidate the origin and transport of Yedoma sediments during the Late Pleistocene, we sampled such well-preserved syngenetic ice wedges without any thaw-related disturbance from Batagay (upper ice complex, UIC) and Central Yakutia (Cyuie and Churapcha) (Supplementary Notes; Supplementary Figs. 2, 3; Supplementary Table 1). Formation ages of the ice wedges studied in the middle part of the Batagay UIC (~60–27 ka; ref. 24) fall within the middle to late (mid-late) Marine Isotope Stages (MIS) 3 (~45–30 ka), based on newly available radiocarbon ages (see Methods; Supplementary Notes; Supplementary Table 2), which corresponds to the period before the Last Glacial Maximum (LGM). In contrast, the Central Yakutia ice wedges formed in early to middle MIS 2 (~21–18 ka; refs. 21,22), representing the LGM period.

In this study, we constrain the provenance and transport pathways of inorganic solid particles preserved in syngenetic ice wedges from Batagay and Central Yakutia by measuring trace element concentrations, particularly rare earth elements (REE), and isotope compositions of Sr, Nd, and O, which have been widely used to trace provenances for various types of materials (e.g., rock, sediment, and dust; refs. 1,25,26,27). Formed mainly by snow compaction and hoar frost accretion22,23, rather than surface water infiltration, our ice wedge samples effectively trap atmospheric input and are expected to preserve atmospheric dust signals. Unlike host sediments, well-preserved ice wedge infill is less affected by redeposition and sediment reworking, and thus can more directly record atmospheric deposition during ice wedge growth. Together, these datasets allow us to evaluate the relative contributions of local versus LRT sources during the mid-late MIS 3 (Batagay) and the LGM (MIS 2; Central Yakutia).

Results

Grain-size distributions and trace element systematics

The possibility of containing the LRT dust in both regions is suggested by the bimodal to polymodal grain-size distributions, indicative of a polygenetic origin (Supplementary Notes; Supplementary Fig. 4; Supplementary Data 1; refs. 12,28). Alongside a dominant silt-sized peak (10–80 μm range), all samples display a secondary mode or shoulder around 1 μm, consistent with the fine-grained LRT dust. Because REE, Sr isotope, and quartz oxygen isotope compositions can be grain-size dependent25,26,29, we therefore analyze size-separated fractions (see Methods).

REE concentrations varied across grain-size fractions but did not exhibit any systematic correlation with particle size (Supplementary Fig. 5; Supplementary Data 2). Differences between sample groups are more prominent than those among grain-size fractions, indicating that the REE concentrations are mainly controlled by provenance rather than grain-size effects. To compare the geochemical compositions of all samples and possible source areas (PSAs), Fig. 2 and Supplementary Fig. 6 present selected element ratios, including REE and other trace elements. These data are normalized to post-Archean Australian shale (PAAS; ref. 30). The element ratios used in these binary plots have been suggested as useful particle size-independent provenance tracers25.

Fig. 2: Bivariate plots of selected element ratios.

a (Gd/Er)PAAS vs (La/Gd)PAAS. b (Gd/Er)PAAS vs (La/Yb)PAAS. c Y/TbPAAS vs (Eu/Eu*)PAAS. d Y/∑REEPAAS vs (Eu/Eu*)PAAS. The subscript PAAS indicates that element concentrations were normalized to post-Archean Australian shale30. Red, blue, and purple symbols denote ice wedge sediments from Batagay, Cyuie, and Churapcha, respectively, covering both bulk and grain-size fractions. Possible source areas (PSAs) are shown as shaded or outlined fields, including surface mud crusts from the northern margin of the Tibetan Plateau (NMTP), sands from the Taklimakan Desert (TK), sands from the Ordos Plateau (OP), loess from the Eastern Qinling Mountains (EQLM), loess and paleosols from the Chinese Loess Plateau (CLP), suspended particulate matter (SPM) from the Lena and Yana Rivers, and local lithologies such as Aldan Shield lamproites and Kharaulakh-Suordakh magmatic rocks, together with comparable datasets from Yukon loess and St. Elias Mountain dust of probable Asian origin. Data sources and information for PSAs and comparable regions are provided in Supplementary Data 2.

The binary plots suggest that our samples represent mixtures of LRT dust from the Chinese sources and locally derived material, with varying proportions (Fig. 2; Supplementary Fig. 6). The Batagay samples overlap closely with the Ordos Plateau (OP) field, suggesting a relatively stronger input from Chinese dust. In contrast, the Central Yakutia samples fall between the fields of the Lena and Yana River suspended particulate matter (SPM) and the northern China aeolian deposits, probably the OP, suggesting a mixture of regional input and distal dust sources. Europium anomalies (Eu/Eu*) also distinguish the two groups (Supplementary Notes; Supplementary Fig. 5): the Batagay samples overlap with the OP to the northern margin of the Tibetan Plateau (NMTP) fields, whereas the Central Yakutia samples align more closely with the Lena and Yana SPM (Fig. 2c, d; Supplementary Fig. 6a–e). In detail, the Lena River SPM systematically shifts toward the Aldan Shield lamproites across all element ratios. The Central Yakutia samples show a similar deviation, trending toward both the Lena River SPM and the Aldan Shield lamproites, likely demonstrating a substantial contribution from this crystalline lithology of the Aldan Shield.

Strontium and neodymium isotopes

The measured 87Sr/86Sr ratios of bulk ice wedge sediments range from 0.715784 to 0.717100 in Batagay and from 0.716610 to 0.717291 in Central Yakutia (Fig. 3; Supplementary Data 3). The fine fractions (<10 and <5 μm) show higher ratios, ranging from 0.719567 to 0.722099 in Batagay and from 0.720124 to 0.721181 in Central Yakutia, indicating a grain-size effect with increasing 87Sr/86Sr toward the fine fractions. In contrast, εNd(0) values show little evidence of grain-size effects but display clear regional differences. Bulk εNd(0) values from the Batagay samples (−8.8 to −8.5) are higher than those from the Central Yakutia samples (−18.3 to −16.3), and the fine fractions show the same contrast (−8.2 to −7.9 in Batagay and −14.7 to −14.3 in Central Yakutia), highlighting the distinct provenance signals between the two sites.

Fig. 3: Sr-Nd isotope compositions.

a ɛNd vs 87Sr/86Sr for ice wedge sediments from Batagay, Cyuie, and Churapcha, shown for bulk samples (dominant grain-size modes of 50−80 μm, 20−70 μm, and 10−60 μm, respectively; see Supplementary Notes and Supplementary Fig. 4) as well as fine fractions (<10 and <5 μm), compared with possible source areas (PSAs) and reference datasets shown as shaded or outlined fields. εNd denotes present-day values, εNd(0), calculated as [(143Nd/144Nd)sample/(143Nd/144Nd)CHUR − 1] × 104, with (143Nd/144Nd)CHUR = 0.512638 (CHUR = chondritic uniform reservoir). PSA fields include sands from the Mongolian Gobi (MG), northern boundary of China (NBC), Ordos Plateau (OP) and northern margin of the Tibetan Plateau (NMTP); loess from the Chinese Loess Plateau (CLP), Eastern Qinling Mountains (EQLM) and northeastern (NE) Siberia; Laptev Sea sediments (LSS); and local lithologies (Aldan Shield and Kharaulakh-Suordakh volcanic rocks, whole-rock data). Additional datasets comprise European loess, Greenland dust, and Alaskan and Yukon loess. Data sources and information are given in Supplementary Data 3. b Same plot showing mixing relationships among PSAs. Dark-gray shading indicates the isotopic mixing range among the Chinese PSAs, whereas light-gray shading represents the extended mixing field between the Chinese PSAs and local lithologies. Red and blue dashed lines denote mixing trajectories between representative PSAs. Gray labels indicate the mixing ratio used to define a hypothetical endmember (e.g., 2:8 denotes an endmember obtained by mixing NMTP and OP in that proportion). Mixing trajectories were calculated using Sr and Nd isotope ratios and elemental concentrations of the respective endmembers (see Methods). For the NMTP, Sr and Nd concentrations were taken from the Taklimakan Desert (TK) sand data (ref. 67).

The Batagay values appear near the isotope range of the Greenland glacial dust and are plotted between the ranges of Mongolian Gobi Desert (MG), northern boundary of China (NBC), Chinese Loess Plateau (CLP), and the NMTP (Fig. 3). This similarity in radiogenic isotopes suggests that our samples may share a common East Asian provenance with the Greenland glacial dust. Numerous studies have reported that dust in the Greenland ice cores originated from major dust-emitting sources in East Asia (e.g., CLP, MG, and the Taklimakan Desert, TK) during the last glacial period31,32,33. Although dust of European origin has also been reported in Greenland, and Sr-Nd isotopes alone are insufficient to distinguish between European and Asian sources34, geographical distance makes it more plausible that atmospheric dust deposited in our study area originated from East Asia.

Bulk Central Yakutia samples overlap with the isotopic range of Alaskan loess, and the εNd(0) values overlap with those of Aldan Shield lamproites (Fig. 3). Although direct comparison of Sr isotopes is complicated by grain-size effects and by the use of whole-rock lamproite values26,35, the εNd(0) overlap suggests a contribution from the Aldan Shield-derived material. The more radiogenic 87Sr/86Sr in the Central Yakutia samples relative to the lamproites can be attributed to grain-size partitioning26. However, the fine fractions (<10 and <5 μm) shift for both Sr and Nd isotopes toward the fields of the Eastern Qinling Mountains (EQLM) and OP, suggesting additional contributions of Chinese dust, consistent with mixing signatures inferred from REE ratios (Fig. 2).

Although the reported 87Sr/86Sr values for the Lena and Yana River SPM (0.716624 to 0.717728 for the Lena; 0.713647 to 0.716892 for the Yana and its tributary; ref. 36) are similar to those of our samples, Nd isotopes for these rivers are scarce. We alternatively compare our dataset with the Laptev Sea sediments (LSS) and northeastern (NE) Siberian loess as more integrative local references (Fig. 3). The LSS cannot be regarded as a single definitive local endmember, because the LSS represents a mixture of sediments derived from diverse sources, including coastal erosion and river discharge, and because riverine SPM draining limited watersheds can have different chemical compositions37,38. Nonetheless, the isotopic similarity between LSS and all samples suggests that local sources likely contributed to the ice wedge sediments. Similarly, NE Siberian loess exhibits Sr isotopes comparable to our samples at equivalent size fractions (<5 μm).

However, NE Siberian loess shows Nd isotope compositions distinct from the samples. This εNd(0) discrepancy indicates that the Nd systematics require additional components, consistent with mixing between local lithologies and LRT dust. Indeed, all samples fall within the mixing field defined by several PSAs (Fig. 3b). The Batagay samples cannot be explained solely by local lithologies such as the Kharaulakh-Suordakh volcanic rocks or Aldan Shield lamproites. Because those reference values are based on whole-rock measurements and therefore not directly comparable to the fine sediment fractions, we treat them as bounding constraints on the local compositional range. The Batagay values are consistent with mixing between a set of hypothetical local endmembers (reflecting Kharaulakh-Suordakh and Aldan Shield contributions) and the OP-NMTP mixing endmembers (red dotted lines in Fig. 3b). For Central Yakutia, mixing trajectories between the OP-NMTP mixing endmembers and a local endmember defined by the bulk samples can explain the observed shift toward finer fractions (blue dotted lines in Fig. 3b). The bulk Central Yakutia compositions are taken as an approximation of the local source, showing Aldan Shield-like εNd(0) values and elevated 87Sr/86Sr compared to lamproites. Taken together, these results suggest that both study sites may contain LRT dust from the OP and NMTP in addition to local material. This interpretation is broadly consistent with the REE-based evidence (Fig. 2), which also points to contributions from Chinese PSAs. We note, however, that REE and Sr-Nd isotope signatures can be influenced by distinct grain-size and mineralogical fractions, which may account for the minor discrepancies observed between the two provenance signals. The REE data suggest OP as the dominant source (Fig. 2), whereas the Sr-Nd isotopes imply mixing with additional components (Fig. 3).

Quartz oxygen isotopes

Oxygen isotope compositions of quartz grains are effective in differentiating distinct provenances in well-mixed aeolian dust1. Since δ18O values of aeolian quartz are influenced by grain-size effects with low-temperature microcrystalline quartz29, we report grain-by-grain isotopic compositions subdivided by grain-size to avoid such effects (see Methods; Fig. 4; Supplementary Data 4). The δ18O (vs V-SMOW) distributions of all grain-size fractions are generally characterized by two dominant isotope ranges, around +10‰ and +15‰ (Fig. 4). The overall distributions differ between Batagay and Central Yakutia, and these differences are confirmed by statistical analyses (see Fig. 4 caption). For three coarser fractions (>32 μm), the Batagay samples show relatively wide peak areas from +10 to +20‰, while the Central Yakutia samples have primary peaks at +10‰ and secondary peaks around +15‰. The finest fractions (<32 µm) exhibit broad distributions at all sites, suggesting the mixing of various grain provenances. In particular, the dominant peak areas for the finest fraction of the Batagay samples shift toward values around +15‰, indicating distinct source contributions compared with the Central Yakutia samples.

Fig. 4: Distribution of δ18O values for quartz grains across four size fractions at Batagay, Cyuie, and Churapcha.

a > 125 µm. b 63−125 µm. c 32−63 µm. d <32 µm. Raw data are provided in Supplementary Data 4. Distributions are shown as Gaussian kernel density estimates (KDEs; see Methods). Statistical tests (Welch t-test, Mann–Whitney U test, permutation test) indicate highly significant differences between the Batagay and Central Yakutia (Cyuie and Churapcha) distributions across all size fractions (p < 0.001), whereas no significant differences were detected between the Cyuie and Churapcha samples (p ≈ 0.8). The blue-shaded band marks a peak range primarily attributable to local sources, whereas the yellow-shaded band represents a peak range that may reflect contributions from both local sources and long-range transported (LRT) dust originating from the Chinese arid regions. Outlined fields indicate the δ18O ranges from the Ordos Plateau to the northern margin of the Tibetan Plateau (OP-NMTP; orange dashed) and those of the Mongolian Gobi (MG; brown dashed), from ref. 27, illustrating overlaps with peaks observed in the samples.

The main +10‰ peaks of the measured δ18O values in Fig. 4 are attributed to the local influx of granitic or metamorphic components. The Batagay sampling site is located in the Yana River basin along the Verkhoyansk Range, where granitoid intrusions are common in the central part of the Verkhoyansk-Kolyma orogen (Supplementary Fig. 1; ref. 39). Through glacial erosion or river discharge, sediments from the Verkhoyansk Mountains could contain significant amounts of quartz grains with δ18O values of approximately +10‰, which may have been generated by weathering of granitoid intrusions. Similarly, in the Lena and Aldan River basins, where the Cyuie and Churapcha sampling sites are located, glaciogenic sediments and morainic materials from the Verkhoyansk Mountains during the Pleistocene are common40. Besides, the Lena River flows from the southern part of the Siberian Platform consisting of Proterozoic metamorphic rocks and granitoids41, and the Aldan River originates in the Aldan Shield which comprises Archaean metamorphic series40. Therefore, the primary sources of Lena and Aldan Rivers are crystalline rocks of metamorphic origins, may comprising the δ18O values of +10‰, which is also consistent with the radiogenic 87Sr/86Sr and εNd(0) signatures observed in the Central Yakutia samples.

The different peaks of δ18O values in the sedimentary component range of approximately +15‰ seen in Fig. 4 could be contributed by both local sources and LRT dust. The basement of central and northern Yakutia is mostly Mesozoic sedimentary rocks and Late Paleozoic and Mesozoic folded sediments, respectively (Supplementary Fig. 1; ref. 40). Either glacial erosion of local bedrock or aeolian transport like saltation from local sediments can provide a basis for δ18O values near +15‰. The sedimentary components also likely contain LRT quartz from the Chinese arid regions, as previously suggested by REE and Sr-Nd isotope systematics. In this regard, the ranges of δ18O values taking into account particle sizes from the MG (+15.1 to 17.3‰ for <16 µm; +12.6 to 15.1‰ for 16−63 µm) and the OP and NMTP (+16.1 to 18.5‰ for <16 µm; +14.0 to 17.0‰ for 16−63 µm) as reported in ref. 27 overlap with the range in which peaks for sedimentary components appear in the samples. Notably, the δ18O distribution of the Batagay samples has a prominent peak in this range (Fig. 4), suggesting that LRT dust may have contributed particularly significantly from the Chinese arid regions to Batagay during the mid-late MIS 3. In addition, for the Batagay samples, the dominant δ18O peak near +15‰ is more pronounced in the finer fractions of quartz particles than in the coarser fractions (Fig. 4d). This may reflect a larger proportion of LRT dust in the finer fractions. The contrasting proportions of +15‰ peaks between the Batagay and Central Yakutia samples likely reflect different contributions of LRT versus local material, with the Batagay distribution showing a stronger LRT component and the Central Yakutia distributions retaining more local input, consistent with the stronger local signatures in the Central Yakutia samples inferred from REE ratios (Fig. 2) and Sr-Nd isotopes (Fig. 3).

Discussion

Local sediment accumulation in Yedoma deposits primarily occurred when the surface was unfrozen, relatively dry, and exposed to aeolian processes11,13. Accordingly, a large portion of the Yedoma sediments in this region is interpreted as aeolian in origin, although there is no consensus on this and other interpretations exist regarding their origin and deposition processes, especially for sites further north16,42. At the Batagay site, sediment transport was mainly through saltation from the surrounding Verkhoyansk Mountains43, whereas the Central Yakutia site received sediment primarily via deflation from the Lena-Aldan river plain40. Compared to the relatively warm mid-late MIS 3 period, the extremely cold and dry climate of the LGM (MIS 2) likely enhanced deflation, thereby promoting greater local loess accumulation. This likely resulted in a stronger influx of local material in Central Yakutia during the LGM, thereby diminishing the relative contribution of LRT dust in the geochemical signals.

The presence of LRT dust within ice wedges reflects the synoptic-scale atmospheric circulation. Modern climatology offers useful analogues: episodic transport of East Asian dust to the Arctic across Siberia has been observed, particularly via meridional pathways associated with low-pressure systems over northern Japan in spring44,45. In such cases, dust originating from the Gobi-Taklimakan region is lofted by spring cold-surge fronts into the mid-troposphere and couples with the enhanced cyclonic flow over Northeast Asia. Case studies integrating satellite observations, ground-based lidar, trajectory modeling, and reanalysis data demonstrate that dust storms can follow a meridional shortcut over the Korean Peninsula and Japan, traverse eastern Siberia within days, and deposit even coarse (10−50 μm) particles in the Arctic44,45. A similar mechanism may have been operated during the Late Pleistocene, particularly under the relatively warmer and climatically more variable conditions of the mid-late MIS 3, when the Batagay ice wedges formed.

The enhanced LRT dust signal observed in the Batagay ice wedge samples (mid-late MIS 3) likely reflects synoptic-scale atmospheric circulation. Previous studies have shown that this period was characterized by frequent millennial-scale climate oscillations, which likely influenced atmospheric circulation and dust fluxes46,47,48,49. However, the temporal resolution of our samples does not allow us to resolve individual millennial-scale events or seasonal variations. Nevertheless, the consistent geochemical similarities with the Chinese aeolian deposits across samples of different ages suggest that episodic meridional dust transport toward East Siberia occurred repeatedly throughout the mid-late MIS 3, regardless of millennial-scale variability. We therefore suggest that the Batagay ice wedges may record evidence of such recurrent transport events facilitated by synoptic-scale circulation before the LGM (Fig. 1; Supplementary Fig. 7a).

In contrast, the Central Yakutia samples (LGM; MIS 2) exhibit a stronger local geochemical signature, despite model simulati