Introduction

Earth’s orbital variations drove periodic Pleistocene climate oscillations between glacial and interglacial conditions1,2,3,4. During cold glacials, continental-scale ice sheets developed in the mid- to high-latitudes of North America and Eurasia, and the Antarctic Ice Sheet expanded beyond its present-day margins5,6,7. These ice sheets “retreated” to present-day configurations, or beyond, during warm interglacials5,6. Glacial oscillations notably transitioned from a lower-amplitude 41-thousand years (kyr) cyclicity to a higher-amplitude mean ~100-kyr cyclicity between ~1.25 and ~0.6 million years ago (Ma)6,8,9,10,11,12,13. This climatic event, known as the mid-Pleistocene transition (MPT), featured substantial Northern Hemisphere ice sheet expansion during glacials6,8,9,10,11,12,13, and was accompanied by distinct ocean thermohaline circulation weakening at ~0.9 Ma14,15. The MPT was also associated with prominent inland Asian aridification and desertification, which may have enhanced global cooling through dust emission, iron fertilization of oceans16,17, and an atmospheric CO2 concentration decrease linked to increasing dust-induced oceanic productivity18,19.

After more than two decades of research, the drivers of climate change across the MPT remain debated11,12. The MPT occurred without a concomitant shift in the orbital forcing rhythm, which suggests that the transition may have been caused by nonlinear internal feedback rather than by external (insolation) forcing11,20. Several prominent hypotheses have been proposed to explain the cause of the MPT, including gradual lowering of atmospheric CO2 concentrations, sustained regolith removal from North America and Europe, long-term nonlinear feedbacks between ice sheets and global climate, and their combined effects11,20,21,22,23. Furthermore, there is no consensus about the forcing and/or pacing mechanisms for the mean ~100-kyr glacial cycles from the MPT onward: they are interpreted (i) as short eccentricity (95-kyr and 125-kyr) cycles, (ii) as 2–3 obliquity (41-kyr) cycles, (iii) as 4–6 precession (19-kyr and 23-kyr) cycles, or (iv) as semi-random fluctuations without true periodicity11,24,25,26. In addition, the timing of the MPT remains debated; an increasing number of studies appear to support a gradual transition between ~1.25 Ma and ~0.6 Ma rather than an abrupt shift at a single point around 0.9–1 Ma11,12,20,27. Finally, our current understanding of the MPT relies primarily on marine records11,12,20,21,22. However, continuous high-resolution terrestrial records through the MPT are rare, which hinders development of a comprehensive understanding of the global nature and drivers of the MPT.

The ~640,000 km2 Chinese Loess Plateau (CLP) extends across the northeastern Tibetan Plateau margin from ~100 to 115°E and from ~34 to 41°N (Fig. 1). The CLP loess-palaeosol deposits yield long-term and near-continuous records of past climate change across the MPT. Dust, which makes up the CLP loess-palaeosol sequence, was transported from the inland Gobi Desert, other nearby sandy deserts, and poorly-vegetated areas by near-surface northwesterly winter monsoon winds28 (Fig. 1a). The grain size of the Quaternary CLP loess-palaeosol sequences reflects dust transport by winter monsoon winds and is generally interpreted as a winter monsoon intensity indicator29,30,31. Stronger winter monsoon winds are associated with coarser dust accumulation on the CLP. Winter monsoon winds are caused by the outflow of cold and dry air from high-pressure cells over the cool Asian continental interior, which blows toward lower-pressure cells over the warmer western Pacific and Indian oceans. Strong winter monsoons during cold global glacials resulted in the deposition of thick loess layers; they consist of a mixture of clays, silts, and fine sands, and are largely unaltered by pedogenesis, with a yellow color32,33. The Asian summer monsoon transports heat and moisture from the western Pacific and Indian oceans toward intense low-pressure (warm) cells over South and East Asia during the boreal summer, and contributes 60–75% of annual precipitation on the CLP34,35. Enhanced summer monsoon precipitation on the CLP during warm interglacials drove more intense pedogenesis, formation of abundant iron oxides (e.g., magnetite, maghemite, and hematite), and red soil development within the yellow loess sequence31,34,36. Since the concentration of pedogenically formed iron oxides can be measured by magnetic susceptibility (χ), CLP loess-palaeosol χ is a much-used proxy for summer monsoon precipitation31,36,37. Stronger pedogenesis during periods of increased precipitation accelerates fine magnetite/maghemite and hematite formation, which causes higher χ values31,36,37. The shift from yellow to red color is caused dominantly by increased pigmentary red hematite formation during intense pedogenesis38. Cyclic stratigraphic alternations between yellow loess and red palaeosol layers are comparable across North China29,31,33,34, and provide a unique, continuous continental archive of orbital- to millennial-scale Asian monsoon and environmental variability. Furthermore, they shed light on both low- and high-latitude processes (i.e., summer and winter monsoon, respectively), and can be used to understand the relationship between global and regional Asian climate changes across the MPT27,29,31,32,33,34,35,36,38,39,40,41,42.

Fig. 1: Site location map.
figure 1

a Map of Asian dust sources (Tibetan Plateau, Gobi Desert, sandy deserts, and wind-eroded land) and aeolian loess deposits in North China28. The distributions of Gobi Desert, sandy deserts, wind-eroded land, and loess were from a recent study28, which refers to the Environmental and Ecological Science Data Center for western China, National Natural Science Foundation of China (http://westdc.westgis.ac.cn). b Map of the Chinese Loess Plateau with the locations of our studied Chaona and Luochuan sections (red stars), and other loess-palaeosol sections discussed here (black stars). The Yellow River system in North China is indicated by the dark blue line. We created these maps with ArcGIS (version 10.7) and Illustrator 2020 software.

Full size image

Here we present new continuous millennial-resolution grain size records for the last 1.6 Myr from two parallel loess-palaeosol sections on the Central CLP to improve understanding of orbital-scale Asian climate variability and dynamics during the transition to the “~100-kyr world”. We combine them with existing χ records from the same sections42,43 and with grain size and χ records from other loess-palaeosol sections29,33,44,45,46 to investigate glacial-interglacial changes in the winter monsoon, summer monsoon, dust and moisture transport, and Asian interior climatic conditions across the MPT, with particular focus on two extreme pulses. We consider these CLP palaeoclimate records within a broader context of existing terrestrial and global palaeoclimate records. We also perform new Earth System model simulations to assess the influence of Northern Hemisphere ice sheet expansion in driving the observed continental-scale Asian glacial climate anomalies across the MPT.

Results and discussion

Insights from CLP palaeoclimate records

To investigate winter monsoon dust transport evolution over the last 1.6 Myr, 982 and 1115 samples were collected for grain size analysis from the Luochuan (109°24′E, 35°48′N; 98.3-m thick) and Chaona (107°12′E, 35°7′N; 110-m thick) loess-palaeosol sections on the Central CLP, respectively (Fig. 1b). The 10-cm sampling interval is equivalent to a temporal spacing of ~1–2 kyr. Analogous to the global chronology for benthic foraminiferal δ18O records from marine sediment cores47, the loess-palaeosol chronology has been established based on different sections/cores across the CLP using orbital tuning, land-sea correlation, and/or grain-size age models, which result in similar ages and the same cycle-to-cycle correlations of loess and palaeosol layers to glacial and interglacial periods defined by the marine benthic δ18O record27,29,31,33,43. Our CLP loess-palaeosol chronology was established by correlating both new median grain size and existing χ records42,43 from the Luochuan and Chaona sections to the marine benthic foraminiferal δ18O record47, with additional support from palaeomagnetic constraints, including boundaries between the Brunhes and Matuyama polarity chrons and the Jaramillo subchron (see Methods and Supplementary Figs. 1, 2). Other CLP loess-palaeosol sections used for comparison were also synchronized to this chronology (Figs. 2 and 3). Our new median grain size records from both the Luochuan and Chaona sections have consistent glacial-interglacial variations that also correlate with those observed in other sections across the CLP (Fig. 2). In addition to median grain size records, we calculated the U-ratio (ratio of 16–44 μm versus 5.5–16 μm particle concentrations)48 and grain-size index (GSI, the ratio of 26–52 μm versus <16 μm particle concentrations)49 for the Luochuan and Chaona sections to assess Asian winter monsoon dust transport variability (Fig. 2). Median grain size, U-ratio, and GSI records from both Luochuan and Chaona have almost identical variability patterns (Fig. 2a–f). To reduce the effects of local changes, and to better reveal large-scale glacial-interglacial winter monsoon dust transport across the CLP, we compiled a median grain size stack based on our new records from Chaona and Luochuan and existing records from the Lingtai33, Jingchuan33, and Baicaoyuan44 sections (Fig. 2; Methods). Our median grain size stack varies consistently with the previous stack of mean grain size of quartz particles (MGSQ) based on data from the Lingtai and Zhaojiachuan sections29 (Fig. 2j, k). Here we use our new multiple grain size records to offer a new, more global, perspective on the dynamics of extreme Asian climate events across the MPT, in contrast to previous CLP grain size studies that focused on the general orbital variability and/or land-sea correlations29,30,31,33,43,44.

Fig. 2: Grain size variations of loess-palaeosol sequences across the Chinese Loess Plateau.
figure 2

Median grain size records (black lines) from the a Luochuan (this study), d Chaona (this study), g Lingtai33, h Jingchuan33, and i Baicaoyuan44 sections on our newly refined loess-palaeosol chronology. b, c U-ratio (ratio of 16–44 μm versus 5.5–16 μm particle concentrations) and grain-size index (GSI, a ratio of 26–52 μm versus <16 μm particle concentrations) from the Luochuan section. e, f U-ratio and GSI from the Chaona section. In both sections, U-ratio records are depicted with magenta lines and GSI records in light blue. j Chinese Loess Plateau median grain size stack. k Chinese Loess Plateau mean grain size of quartz particles (MGSQ) stack based on data from the Lingtai and Zhaojiachuan sections29. l LR04 marine benthic δ18O record47. L-numbers refer to consecutive loess horizons counting back from the present day. Numbers on the benthic δ18O record refer to Marine Isotope Stages, counting back from the present day. Pink bars indicate the correlation of the median grain size of marker loess layers (L15 and L9-1) to glacial stages (MIS 38 and MIS 22).

Full size image
Fig. 3: Magnetic susceptibility variations of loess-palaeosol sequences across the Chinese Loess Plateau.
figure 3

Magnetic susceptibility (χ) records from the a Luochuan43, b Chaona42, c Jingchuan45, d Zhaojiachuan29, e Lantian46, and f Lingtai29 loess-palaeosol sections across the CLP on our refined loess-palaeosol chronology. g Chinese Loess Plateau χ stack (this study). h LR04 marine benthic δ18O record47. S-numbers refer to consecutive palaeosol horizons counting back from the present day. Numbers on the benthic δ18O record refer to Marine Isotope Stages, counting back from the present day. Pink bars indicate the correlation of χ of marker loess layers (L15 and L9-1) to glacial stages (marine isotope stages 38 and 22).

Full size image

To assess summer monsoon variability on the CLP coeval to our grain-size-based winter monsoon record, we developed a new loess-palaeosol χ stack by compiling existing χ records from Luochuan43, Chaona42, Jingchuan45, Zhaojiachuan29, Lantian46, and Lingtai29 (Fig. 3). Both our grain size and χ stacks involve more CLP loess-palaeosol sections than previous stacks29,33, which further minimizes the impacts of local changes. As is the case for grain size records, temporal χ variability matches well among sections and generally correlates cycle-by-cycle with glacial-interglacial cycles in the benthic foraminiferal δ18O record (Figs. 2, 3). Nevertheless, the records contain subtle differences in detail, largely because each represents distinct climate features. In short, χ reflects summer monsoon precipitation, whereas grain size is affected by dust transport, and benthic foraminiferal δ18O is a function of both deep-sea temperature and global ice volume.

Consistent with previous studies27,29,31,33,39, glacial loess layers have overall larger median grain sizes, U-ratios, and GSI values (stronger winter monsoons), and lower χ values (weaker summer monsoons and lower precipitation) than interglacial palaeosol layers (Figs. 2–4). Loess layers L15 (correlating to MIS 38 at ~1.25 Ma) and L9-1 (correlating to MIS 22 at ~0.9 Ma) have notably lower χ values and exceptionally larger median grain sizes, U-ratios, and GSI values than other loess layers (Figs. 2–4). The distinct grain size increases in L15 and L9-1 are observed in all loess-palaeosol sequences from both the eastern and western CLP (Fig. 2), which indicates dominant and widespread dust transport changes during MIS 38 and MIS 22 across the entire CLP. We infer that winter monsoon conditions over Asia during these periods were amplified (i.e., cooler, drier, and windier) compared to preceding and succeeding glacials.

Fig. 4: Terrestrial and global climate variability on glacial-interglacial timescales across the mid-Pleistocene transition.
figure 4

Our newly established Chinese Loess Plateau a χ and b median grain size stacks. S-numbers and L-numbers refer to consecutive palaeosol and loess horizons counting back from the present day, respectively. c Mediterranean sea level record8, d Global sea level reconstruction based on Pacific benthic foraminiferal δ18O and Mg/Ca records10. e Seawater δ18O from ODP site 1123, South Pacific Ocean9. f Global sea level reconstruction based on the LR04 benthic foraminiferal δ18O stack6. g LR04 benthic foraminiferal δ18O record47. Numbers refer to Marine Isotope Stages, counting back from the present day. h global average sea surface temperature (SST) stack70. i Enhanced Asian aridification, new desert formation, and desert expansion from ~1.25 Ma onward inferred from geological records19,28,32,56,58,59. Pink bars indicate the correlation of coarse marker loess layers (L15 and L9-1) to glacial stages (marine isotope stages 38 and 22).

Full size image

To investigate the main (orbital) periodicities of the Asian summer and winter monsoon, and to compare these with global (i.e., high-latitude) climatic variations, we present time-evolutive spectral analyses of the χ, grain size, and marine benthic δ18O stacks. Spectral analyses of the median grain size and χ stacks suggest a major transition from a predominant 41-kyr to a mean ~100-kyr periodicity across the MPT, albeit with subtle differences in the exact expression of the transition (Fig. 5a, b). These differences likely reflect more nuanced regional responses to MPT climate change. CLP precipitation (indicated by χ) is dominated by moisture transport from the western Pacific and Indian oceans to inland Asia by the summer monsoon, which is a low-latitude process, whereas transport of cold and dry air from Eurasia toward the tropical oceans by the winter monsoon (indicated by median grain size) represents a high- to mid-latitude process. Spectral analyses of the benthic δ18O record47 reveal a prominent switch from predominant 41-kyr to mean ~100-kyr cycles across the MPT, with co-occurring 41-kyr and ~100-kyr cycles between ~1.2 Ma and ~0.6 Ma (Fig. 5c). A weak obliquity (41-kyr) cyclicity until ~0.6 Ma suggests that the shift in MPT periodicity was more gradual and delayed in the global mean glacial cycle pattern reflected in the benthic foraminiferal δ18O record than in the CLP precipitation (χ) and winter monsoon (median grain size) records (Fig. 5). The low-frequency spectral power in all three records is not centered exactly on the 100-kyr band before ~0.65 Ma and after ~0.2 Ma (Fig. 5). Spectral power splits broadly into two branches close to 80-kyr and 120-kyr, respectively, during the MPT in all three records, although the timing and subtle features of the split vary. These two branches converge toward a band close to 100-kyr by the MPT termination in all three records, which gradually evolves into a ~120-kyr band. The dominance of mean ~100-kyr power after the MPT coincides with a substantial decline in the 41-kyr band. These observations are broadly consistent with previous findings for the MPT11,25. In addition, we find a 50- to 60-kyr band in the χ stack around 1.3 Ma and the median grain size stack between ~1.3 Ma and ~1.4 Ma, which is not distinct in the benthic δ18O stack (Fig. 5). Overall, an apparent shift from a predominantly 41-kyr to mean ~100-kyr orbital periodicity across the MPT in the CLP χ, median grain size, and marine benthic δ18O stacks suggests that orbital-scale Asian summer and winter monsoon variations are closely linked to glacial-interglacial pacing of Northern Hemisphere ice sheets.

Fig. 5: Orbital-scale climate variability across the mid-Pleistocene transition.
figure 5

Evolutive spectra of the a Chinese Loess Plateau χ stack, b Chinese Loess Plateau median grain size stack, and c LR04 benthic foraminiferal δ18O record47 with a 320-kyr sliding window and 2.6-kyr step. Red dashed lines indicate 41-kyr and 100-kyr periodicities.

Full size image

Mechanisms for CLP loess coarsening across the MPT

Our multiple new grain size records reveal both orbital-scale variability and extreme pulses across the MPT, in agreement with previous records29,33,39,40. Two prominent loess grain size anomalies at ~1.25 Ma (L15) and ~0.9 Ma (L9-1) have been explained previously in terms of phased Tibetan Plateau uplift39. However, evidence for major plateau uplift during the MPT is tenuous. The Tibetan Plateau was already close to its present-day elevation and configuration by at least the late Miocene, with only limited and more regional Quaternary adjustments50,51,52,53,54,55. Thus, plateau uplift cannot explain the distinct coarsening of L15 and L9-1 across the CLP, nor their astronomical pacing. This leaves the cause(s) of L15 and L9-1, and their palaeoclimatic significance, open to further investigation33.

Fundamentally, these exceptionally coarse loess layers must reflect a combination of (i) widespread wind strength increase, (ii) transport pathway shortening due to enhanced and expanded Central Asian aridity, (iii) enhanced coarse dust production through increased aridity and sediment availability, and/or (iv) reduced vegetation cover with lower soil stability and greater soil erosion by wind during glacials MIS 38 and MIS 22. These terrestrial changes probably became prominent at the onset of (~1.25 Ma) and halfway through (~0.9 Ma) the MPT. For example, new sandy deserts (e.g., Badain Jaran Desert, Tengger Desert) are known to have formed at ~1.2–0.9 Ma to the north of the CLP28,56, while existing sandy deserts (e.g., Mu Us Desert) expanded southward at ~1.25 Ma32 (Fig. 4i). Sandy desert environments first appeared in the Hobq Desert at ~1.3–1.2 Ma, replacing preceding fluvio-lacustrine environments57. The Tarim and Qaidam basins (northwestern China) also experienced increased aridity from ~1.25 Ma onward19,58,59 (Fig. 4i). Increasing carbonate δ13C values of the Kazakhstan loess sequence suggest that Central Asia aridified across the MPT60. Glacial loess deposits in Tarim Basin, which were deposited under the control of mid-latitude Westerlies, became generally coarser across the MPT18. The loess coarsening and parallel dust flux increases in the Tarim Basin are consistent with the expansion of Central Asian arid regions across the MPT18. Moreover, Tibetan Plateau cooling across the MPT increased physical weathering intensity, which in turn produced more sand-sized material available for transport to western China through riverine and rain erosion19. These Asian arid regions, especially the neighbouring Badain Jaran, Tengger, Mu Us, and Hobq deserts (see Fig. 1a), provided coarse dust sources for the CLP28. We propose that desert expansion to the west and north of the CLP together with winter monsoon wind strengthening and potential dust source changes led to loess coarsening on the CLP. This agrees with previous results of Sr and Nd isotopic compositions and detrital zircon U-Pb dating of CLP loess-palaeosol sediments, which track dust provenance changes during the MPT, with increased contributions from the NE Tibetan Plateau and Gobi Altay Mountains61,62,63,64,65. The inferred Asian aridification and winter monsoon strengthening is also supported by the extension of aeolian loess deposits to the south of the Qinling Mountains and Yangtze River catchment in southeast China across the MPT66.

We note that intensification of glacial Asian climate and environmental conditions coincided with Northern Hemisphere ice sheet expansion at the onset and middle of the MPT, when expression of mean ~100-kyr glacial cyclicity initiated and enhanced, respectively6,8,9,10,13,67. Both marine and terrestrial data suggest that glacial Northern Hemisphere ice sheets expanded substantially at the beginning of, and halfway through, the MPT. For example, various sea level reconstructions suggest notable low stands during MIS 38 and MIS 22 relative to preceding glacials, albeit with subtle amplitude differences among reconstructions that relate to variable uncertainties in different methods6,8,9,10 (Fig. 4c–f). In addition, 26Al-10Be burial dating of tills suggests that the Laurentide Ice Sheet advanced to its extreme southern limit (~40°N) by ~1.3 Ma67. The ODP Site 887 magnetic susceptibility and Deep Sea Drilling Project (DSDP) Site 607 carbonate concentration records suggest that Northern Hemisphere ice sheets expanded and shed more ice-rafted debris into the Gulf of Alaska at ~1.3 Ma and Central North Atlantic Ocean at ~0.9 Ma, respectively68,69. Associated with ice sheet expansion, glacial global average sea surface temperatures decreased markedly after ~1.25 Ma70 (Fig. 4h). We infer that Asian glacial climate intensification at ~1.25 and ~0.9 Ma as indicated by CLP loess coarsening may be linked to concomitant shifts to greater glacial Northern Hemisphere ice sheet expansion in the context of enhancing global cooling across the MPT.

To assess the effect of Northern Hemisphere ice sheet expansion on Asia’s interior glacial conditions, aridification and desertification, we performed climate model sensitivity experiments. We used the Community Earth System Model (CESM 1.2) to perform two experiments to examine the Asian (hydro-) climate response to Northern Hemisphere ice sheet expansion (see Methods). One experiment (Large Ice Sheet (LIS) experiment) was conducted with Northern Hemisphere ice sheet distributions, orography, vegetation, lakes, aerosol conditions, and orbital parameters, set at Last Glacial Maximum (LGM, ~20 ka) values5,71, and with 220 ppm CO2 and 450 ppb CH4 concentrations that are similar to values during glacial MIS 50 (~1.5 Ma) as reconstructed from Antarctic ice cores72. The second experiment (Small Ice Sheet (SIS) experiment) was conducted with the same boundary conditions as the LIS experiment, but with a smaller Northern Hemisphere ice sheet distribution, for which we employ the early Holocene/post-LGM configuration at 13 ka5 (Fig. 6a). Northern Hemisphere glacial ice sheet expansion across the MPT was similar in amplitude to the change between 13 ka and the LGM (~20 ka) according to sea level and benthic foraminiferal δ18O records6,8,9,10,47, although exact ice sheet configurations were different between the ~20–13 ka and MPT periods67. Hence, our experiments are not representative of the full range of changing MPT boundary conditions, which remain poorly constrained, but are instead sensitivity experiments that provide insights into the likely climate variations during the intense glacial events, similar to those that coincided with L15 and L9-1.

Fig. 6: Simulated Asian climate and atmospheric circulation responses to ice volume increase.
figure 6

a Northern Hemisphere ice sheet distribution used in the Small Ice Sheet (SIS, upper) and Large Ice Sheet (LIS, lower) experiments. The yellow dot represents the Chinese Loess Plateau (CLP). Simulated changes in b annual temperature, c annual precipitation, d annual net moisture (precipitation minus evaporation), e Asian high-pressure cell during winter, and f winter monsoon (700 hPa winds) due to ice sheet expansion (SIS – LIS differences). LIS:SIS ratios of atmospheric g annual dust emission flux and h annual dust loading. Solid green contours in b–h denote the 3000 m topographic contour, which includes the Tibetan Plateau. Red stars represent the Chaona and Luochuan loess-palaeosol sections. Small red dots in b–h denote regions with statistical significance above the 95% confidence level (Student’s t-test). Winter in the model is represented by December to February.

Full size image

Differences between our SIS and LIS experiments suggest that Northern Hemisphere glacial ice sheet expansion led to a lowering of Asian mean annual temperature, precipitation, and net surface moisture (precipitation minus evaporation) (Fig. 6b–d), which facilitated intensification and expansion of Central Asian aridity and increased dust production. Mean annual precipitation in the LIS experiment decreased by ~14% in arid inland regions (60–100°E, 30–60°N) and by ~9% in East Asian monsoon regions (100–120°E, 20–40°N) relative to the SIS experiment. Furthermore, ice sheet expansion strengthened Asian high-pressure cells and winter monsoon circulation (Fig. 6e, f), which enhanced winter monsoon dust transport capacity toward the CLP. Overall, annual dust fluxes sourced from arid regions to the north and east of the Tibetan Plateau increased by up to an order of magnitude from the SIS to LIS experiments (Fig. 6g). This is associated with an approximate doubling of the annual atmospheric dust loading over East Asian down-wind regions (Fig. 6h), which is comparable to the largely doubled median grain sizes of L15 and L9-1 relative to adjacent loess layers observed across the CLP (Figs. 2, 4b). Dust changes are substantially larger than precipitation changes between the SIS and LIS experiments (Fig. 6), consistent with the markedly larger grain size changes compared to χ changes at ~1.25 Ma and ~0.9 Ma (Figs. 2, 3).

Our adopted CESM 1.2 experiments do not include dynamic vegetation, nor glaciogenic dust responses. Hence, our dust inferences represent minimum responses because (i) Northern Hemisphere ice sheet expansion was associated with increased glaciogenic dust73 and (ii) accompanying temperature and precipitation lowering (with arid zone expansion) may have further decreased vegetation cover46,74, which in turn would have reduced soil stability, thereby facilitating erosion and dust production and availability. Other potential dust-producing processes not included in the model, such as enhanced physical weathering and rock fracturing through intensified frost wedging and/or glacial grinding, would also produce more dust material that becomes available for transport under colder MPT conditions19. Regardless, our model estimates, together with our land-sea palaeoclimate correlation, strongly support the hypothesis that intensification and southeastward expansion of Asian aridity increased coarse dust availability, and that strengthening of winter monsoon winds increased coarse dust transport, leading to an overall loess coarsening across the CLP, in response to phases of Northern Hemisphere ice sheet expansion during the MPT.

The CLP loess grain size records reveal different regional response intensities to boundary condition changes. Notably, loess L15 at ~1.25 Ma is coarser than L9-1 at ~0.9 Ma, which contrasts with the larger ice volume and cooler climate during the latter period (Fig. 4). This suggests that the CLP grain size response was more pronounced at the MPT onset when the global climate cooled sufficiently for Northern Hemisphere ice sheets to reach a critical size for the first time, which allowed abundant coarser particles to be transported by stronger winds. For example, the initial glaciation-induced transition from fluvio-lacustrine to desert environments in the arid Asian interior would offer greater dust material availability for transport to the CLP than the later, more sustained sandy desert conditions without fluvial-lacustrine processes. Fluvial-lacustrine conditions generally produce abundant dust for aeolian transport once water bodies desiccate and sediment becomes exposed75,76. We argue that the distinct coarsening of L15 and L9-1 across the CLP was related to a combination of greater winter monsoon intensity, shorter transport pathways, decreased vegetation cover, and increased availability of freshly eroded and transportable material in source regions. These changes responded to Northern Hemisphere ice sheet expansion to critical sizes across the MPT. These terrestrial responses were unique to the MPT, when global and Asian climates shifted from a lower-amplitude 41-kyr cyclicity to a higher-amplitude mean ~100-kyr cyclicity (Figs. 4, 5), which constituted a major global climatic reorganization.

Integrating observations, land-sea correlations, and model simulations, we propose that Northern Hemisphere ice sheet expansion drove large-scale amplification of Asian glacial conditions at the onset of (~1.25 Ma) and halfway through (~0.9 Ma) the MPT, when expression of mean ~100-kyr glacial cyclicity initiated and enhanced, respectively. These two glacial intensifications were marked by a combination of intensified and expanded Asian aridity, winter monsoon strengthening, and summer monsoon weakening, with distinct coarsening of loess layers L15 and L9-1 across the CLP. The shift from a predominantly 41-kyr to mean ~100-kyr orbital periodicity across the MPT is also apparent in our winter and summer monsoon records, which, more generally, reflect Northern Hemisphere ice sheet control on orbital-scale Asian climate variability, not just on extreme glacial Asian climate events at ~1.25 Ma and ~0.9 Ma. Our study supports a close relationship between the Pleistocene Asian-interior and global climate changes.

Methods

Following surface outcrop removal, 1115 and 982 fresh samples, respectively, were collected from the Chaona and Luochuan sections, Central CLP, from Holocene palaeosol layer S1 (corresponding to MIS 1) to palaeosol layer S22 (corresponding to MIS 55) at 10 cm intervals. These intervals are equivalent to an averaged temporal spacing of ~1–2 kyr that is of higher sampling resolution than in many previous studies. All samples were used for grain size analyses. In preparation for grain size analyses, about 0.2–0.3 g of bulk sample was first pretreated with 30% hydrogen peroxide (H2O2) to remove organic matter and subsequently with 10% hydrochloric acid (HCl) to remove carbonates and iron oxides. After dispersing with 10 mL 10% sodium hexametaphosphate ((NaPO3)6) solution in an ultrasonic bath for ~10 min, the grain size distribution was measured using a Malvern 2000 Laser Instrument at the Institute of Earth Environment, Chinese Academy of Sciences, Xi’an (China). The relative standard deviation of grain size measurement was <3%.

CLP loess-palaeosol age models were established using different approaches, including orbital tuning, land-sea correlation, and grain-size age models, that are all similar and match well with the marine benthic δ18O records27,29,31,33. In particular, the extremely coarse loess layers L15 and L9-1 are consistently correlated to MIS 38 and MIS 22 in these age models27,29,33. In Luochuan, Chaona, and other CLP sections, palaeosol layers, which developed during interglacials, have higher χ values and smaller grain sizes than loess layers that accumulated during glacials27,29,31,33. Based on the correlation of loess (palaeosol) layers to glacial (interglacial) periods27,29,31,33,43, we constructed an age model for the Luochuan and Chaona loess-palaeosol sections by correlating median grain size and χ to the marine benthic δ18O stack47. We correlated iteratively until visually satisfactory cycle-to-cycle correlations were obtained. Final age models for the Luochuan and Chaona sections were established using 37 and 38 age correlation points, respectively. These age correlation points are mostly adjacent to loess–palaeosol boundaries (Supplementary Fig. 1), which facilitated consistent correlation point selection throughout the entire section. The age model is supported by an established magnetostratigraphy for each section77,78, which provides age tie points for the Brunhes and Matuyama polarity chrons, including the Jaramillo subchron, that are consistent with the geomagnetic polarity time scale (GPTS)47 (Supplementary Fig. 1), taking into account uncertainties in the post-depositional natural remanent magnetization (NRM) lock-in depth (generally <50 cm) in CLP loess-palaeosol sediments44. The age model results in a smooth and linear relationship between age and depth for both sections, without abrupt shifts, which provides additional support for our established CLP chronology (Supplementary Fig. 2).

Other CLP loess-palaeosol sections, including Lingtai29,33, Jingchuan33,45, Baicaoyuan44, Zhaojiachuan29, and Lantian46, were synchronized to this chronology by cycle-to-cycle correlation of median grain size or χ data to their counterparts at Luochuan and Chaona. Thus, all sections considered here were placed on the same, newly refined chronology (Figs. 2, 3). After synchronization, we used the interpolating function in the Acycle software79 to conservatively resample the grain size and χ records at 0.5-kyr intervals to obtain an evenly-spaced data series. We constructed the CLP median grain size stack by averaging the evenly-spaced median grain size time series of the Chaona, Luochuan, Lingtai, Jingchuan, and Baicaoyuan sections with equal weight. Similarly, we constructed the CLP χ stack by averaging the evenly-spaced χ time series of Luochuan, Chaona, Jingchuan, Zhaojiachuan, Lantian, and Lingtai sections with equal weight. The χ, median grain size, and LR04 benthic δ18O stacks47 were subjected to spectral analysis to evaluate the robustness of their potential orbital signature. Evolutionary power spectra were calculated using the Acycle software79 with a 320-kyr sliding window and 2.6-kyr step. To improve the expression of the orbital transition across the MPT from 41-kyr to ~100-kyr cycles in the χ and median grain size records, their longer trends were removed with a low-band-pass filter; the residual records were used for evolutionary power spectral analysis.

The Community Earth System Model (CESM 1.2) was used to test the sensitivity and underlying Asian climate response dynamics to Northern Hemisphere ice sheet expansion across the MPT. The CESM consists of coupled dynamic atmosphere, ocean, land, sea-ice, and land-ice components80. We used the Community Atmosphere Model (CAM), Community Land Model (CLM), Parallel Ocean Program (POP), Community Sea-Ice Component (CICE), and Coupler modules in the CESM. The CAM has a Bulk Aerosol Model (BAM) parameterization of dust emission, transport, and deposition81, which has been shown to simulate well dust emission flux and loading in Asia82. The atmosphere has 26 vertical layers and ~0.9° (latitude) × 1.25° (longitude) horizontal resolution. The land model has 15 soil layers and the same horizontal resolution as the atmosphere. Ocean and sea-ice components have 60 vertical layers with 0.5° horizontal resolution. To evaluate Asian (hydro-) climate response to MPT Northern Hemisphere ice sheet expansion, we performed two sensitivity experiments: the SIS and LIS experiments (see main text for boundary condition details). In contrast to LGM values (185 ppm CO2 and 350 ppb CH4), which are often used in simulations82,83,84,85,86, we used 220 ppm CO2 and 450 ppb CH4 concentrations that are similar to their values during glacial MIS 50 (~1.5 Ma) as reconstructed from Antarctic ice cores, which may better reflect early Pleistocene glacial greenhouse gas concentrations during the MPT72. In our experiments, only Northern Hemisphere ice sheet distributions were variable while keeping other boundary conditions fixed. This model design enables an understanding of how Asian (hydro-) climate responds exactly to Northern Hemisphere ice sheet expansion. The ice volume difference between SIS and LIS experiments is broadly in agreement with the amplitude of change across the MPT according to sea level and benthic foraminiferal δ18O records6,8,9,10,47. Both experiments were integrated for 150 model years from the equilibrated LGM initial conditions based on the Palaeoclimate Modeling Intercomparison Project (PMIP3) LGM experiment (http://pmip3.lsce.ipsl.fr). Climatological means of the last 50 model years were used here.

Our modeling differs from previous efforts in the following respects. Previous simulations generally used the full range of LGM boundary conditions82,83,84,85,86. Although our experiments employ LGM orography, vegetation, lakes, aerosol conditions, and orbital parameters, we adjusted CO2 and CH4 concentrations from LGM to early Pleistocene glacial values to better reflect MPT glacial greenhouse gas conditions72. Moreover, previous simulations examined Asian climate responses to the combined influences of global ice volume, orbital forcing, and greenhouse gas concentration changes from the LGM to Holocene, pre-industrial, or present day83,84,86,87,88, yet the response to ice volume change by itself was poorly constrained82,89. In our experiments, only the Northern Hemisphere ice sheet extent is varied and other parameters (e.g., greenhouse gas and orbital forcing) are kept constant, which helps to improve understanding of how Asian (hydro-) climate responds exactly to Northern Hemisphere ice sheet expansion across the MPT. Northern Hemisphere ice sheet expansions from the Holocene, pre-industrial, or present day to the LGM used in previous simulations are useful for understanding the impacts of large-amplitude ice sheet changes between full glacial and interglacial conditions. However, these ice volume changes overestimate the glacial changes across the MPT. The expansion of Northern Hemisphere ice sheets between our SIS and LIS experiments better represents the ice volume changes across the MPT. Despite these model differences, both previous82,84,86,89,90,91 and our new simulations indicate that Northern Hemisphere ice sheet expansion was associated with Asian cooling, winter monsoon strengthening, and aridification. Simulations suggest that large ice sheets can also drive cooling and aridification in North China under the warmer Pliocene boundary conditions92. However, our models do not simulate the wetter-than-present conditions in Central Asia when ice volume was larger, for example, under the LGM maximum Northern Hemisphere ice volume conditions that were observed in some previous simulations93; this response may have arisen from a boundary condition other than ice volume in previous simulations.