INTRODUCTION
Sand seas in the arid and semiarid regions of northern China cover ~600,000 km2 (Zhu et al., Reference Zhu, Wu, Liu and Di1980), accounting for a large portion of the desert landscape in the mid-latitudes of the Northern Hemisphere (Goudie, Reference Goudie2002; Williams, Reference Williams2014). These regions are important in the context of global climate change as potential sources of substantial dust flux, and atmospheric dust exerts important influences on climate state by regulating formation of clouds, energy and water budgets, and biogeochemical cycles (Martin and Fitzwater, Reference Martin and Fitzwater1988; Biscaye et al., Reference Biscaye, Grousset, Revel, Van der Gaast, Zielinski, Vaars and Kukla1997; Miller and Tegen, Reference Miller and Tegen1998; Jickells et al., Reference Jickells, An, Andersen, Baker, Bergametti, Brooks and Cao2005; Uno et al., Reference Uno, Eguchi, Yumimoto, Takemura, Shimizu, Uematsu, Liu, Wang, Hara and Sugimoto2009; Ravi et al., Reference Ravi, D'Odorico, Breshears, Field, Goudie, Huxman and Li2011; Shao et al., Reference Shao, Wyrwoll, Chappell, Huang, Lin, McTainsh, Mikami, Tanaka, Wang and Yoon2011; Bullard et al., Reference Bullard, Baddock, Bradwell, Crusius, Darlington, Gaiero and Gassó2016). Consequently, paleoenvironmental reconstructions in Chinese sand seas have been the subjects of numerous studies over the last decades. Despite a general trend towards increasing aridity during the Quaternary, several studies report the presence of lacustrine deposits in the aeolian sequences in these deserts, suggesting the occurrence of more humid periods in northern China during the Late Quaternary (Kelts, Reference Kelts1992; Hövermann, Reference Hövermann1998; Chen et al., Reference Chen, Wang and Tang2000; X.P. Yang et al., Reference Yang, Liu and Xiao2003, Reference Yang, Ding and Ding2006, Reference Yang, Forman, Hu, Zhang, Liu and Li2016; Yang and Scuderi, Reference Yang and Scuderi2010). For example, lacustrine and fluvial deposits in the southern Taklamakan Desert reveal that a humid period began at ca. 2 ka and ended at ca. 1.5 ka, which is in line with the timing of the ancient cities nearby (Yang et al., Reference Yang, Du, Liang, Zhang, Chen, Rioual, Zhang, Li and Wang2021). Thus, unrefuted evidence of paleoenvironmental changes in these deserts is essential for understanding the climatic and human-activity changes in Asian mid-latitudes (Yang et al., Reference Yang, Ma, Dong, Zhu, Xu, Ma and Liu2010, Reference Yang, Li and Conacher2012; G.Q. Li et al., Reference Li, Jin, Wen, Zhao, Madsen, Liu, Wu and Chen2014; Z.J. Li et al., Reference Li, Sun, Chen, Wang, Zhang, Guo, Wang and Li2014; Lancaster et al., Reference Lancaster, Wolfe, Thomas, Bristow, Bubenzer, Burrough and Duller2016).
The Kumtag Desert (KMD) is located at the westernmost end of the Hexi Corridor, northwestern China (Fig. 1). Due to its remoteness and harsh natural conditions, after the first scientific report made by Hedin (Reference Hedin1903), studies in the KMD region have remained relatively scarce. Xia and Fan (Reference Xia, Fan and Xia1987) first argued that this desert began to form in the middle Pleistocene and it continued to expand northward from the late Pleistocene to Holocene. Tang et al. (Reference Tang, Su, Ding, Zhu, Zhai, Yi, Liu, Zhang and Li2011) and Dong et al. (Reference Dong, Lü, Lu, Qian, Zhang and Luo2012) suggested that the desert became dry during the Late Miocene–Pliocene, and that it had experienced a continuous drying process during the entire Quaternary. Zhao et al. (Reference Zhao, Lu, Zhang, Wang, Yi, Chen, Zhang and Wu2015), however, reported that the KMD experienced two significant climate changes since 23 ka, one being the shifting from cold-dry to warm-humid, and the other being a gradual change from warm-humid to warm-dry conditions. Ding (Reference Ding2017) assumed, based on the studies about sandy loess, that climate in the KMD had been affected by multiple cycles of cold-dry and warm-dry conditions since the early Holocene. Although previous work established a basic framework of paleoenvironmental changes in the KMD, two key issues remain unsolved. First, the current reconstruction of environmental changes in the KMD lacks evidence from sedimentary sequences in the desert itself and is solely based on interpretation of grain-size data. Second, there is a lack of high-quality physical dating for the loess sections near or surrounding the desert, despite their close association with current reconstructions of past environmental changes in the KMD.
Figure 1. Overview of the study area. (A) General map of arid Central Asia (ACA) and climatic records mentioned in the text. MLW, ISM, and EAM are abbreviations for mid-latitude Westerlies, Indian summer monsoon, and East Asian summer monsoon, respectively. (B) Topographical map of the Kumtag Desert (KMD), northwestern China and locations of the KM, YG, and ASK sections.
The eastern edge of the Kumtag Desert consists of rivers and oases interspersed with dunes. In this paper, we present data from several sections located along the eastern edge of the KMD, with chronologies established by optically stimulated luminescence (OSL) dating. We use granular and geochemical data as paleoenvironmental proxies. Based on these results, we aim to reconstruct the paleoenvironment since the late Pleistocene, to discuss possible driving mechanisms and to explore the interaction between environmental change and human activities in the eastern KMD.
REGIONAL SETTING
The KMD is located north of the Altyn-Tagh Mountains, south of the Aqike Valley, east of Lop Nur, and west of the Dang River (Fig. 1). The sand sea has an area of ~19,500 km2 (Zhu et al., Reference Zhu, Wu, Liu and Di1980). The word Kumtag means “sand mountain” in the local language. Various dune types are found in the KMD sand sea, including barchans, linear dunes, net-shaped dunes, star dunes, sand sheets, and distinctive feather-like dunes (Dong et al., Reference Dong, Qian, Yan and Su2010a, Reference Dong, Wei, Qian, Zhang, Luo and Hub, Reference Dong, Qian, Luo, Zhang and Lü2013a, Reference Dong, Zhang, Qian, Luo, Lv and Lub; Qian et al., Reference Qian, Dong, Zhang, Luo and Lu2012, Reference Qian, Dong, Zhang, Luo, Lu and Yang2015; Lü et al., Reference Lü, Narteau, Dong, Rozier and du Pont2017). The KMD is characterized by a hyper-arid, continental climate because the East Asian summer monsoon hardly reaches the area due to the long distance and terrain barriers. The mean annual precipitation recorded at the Dunhuang weather station on the eastern margin of the desert was 42.2 mm during 1971–2000 (Lu et al., Reference Lu, Wu, Dong, Lu, Xiao and Wang2012). Several seasonal streams rising from the alluvial fans of the Altyn-Tagh Mountains in the south flow northwards and vanish within the KMD. The monthly runoff of these rivers varies dramatically due to the concentration of rainfall from May to August in the KMD and Altyn-Tagh Mountains (Kang et al., Reference Kang, Chen, Zhang, Wang, Shang and Cheng2015). For example, the runoff of the Xitugou Stream in May 2015 was 350.9 m3/second and only 53.4 m3/second in April of that same year. Generally speaking, runoff is relatively high between May and September (Chen et al., Reference Chen, Niu, Huang, Zhang, Wang and Zhang2017).
A branch of low mountains occurs on the northern edge of the Altyn-Tagh Mountains, while the southern parts of the KMD overlie these mountainous forelands. Although the eastern KMD is similar to other parts of the desert with regard to dune types and climate, the drainage network extends deep into the sand sea and the paleo-meanders remain continuous and reach the Xihu Wetland in the northeast. This resulted in an extensive occurrence of lacustrine deposits (Fig. 1). In addition, several oases lie along the seasonal streams, one of the largest being the Yangguan Oasis, which was the site of one of the main towns along the ancient Silk Road. The earliest records of human occupation in the eastern KMD has been dated to ca. 4 ka (Liu et al., Reference Liu, Wang, Zhao and Ding2004).
MATERIAL AND METHODS
Sampling
We carried out fieldwork in the eastern KMD, mainly along the margin of the sand sea and along the gullies incising the dune landscape (Fig. 1). Seventy-five samples were collected from three sections for measuring grain size and geochemical compositions. Samples of the aeolian sand were collected at ~10-cm intervals from the sedimentary sections, but with a resolution of 5 cm for the lacustrine sediments intercalated in these sequences. Thirteen samples for OSL dating were collected from the sections to establish a chronological framework. All OSL samples were taken using stainless steel tubes with a 4-cm diameter and a 30-cm length. The sampling position for every OSL sample is shown in Figure 2.
Figure 2. Field photographs of the KM, YG, and ASK sections and OSL sampling sites in the eastern KMD (see Fig. 1 for geographical locations).
Luminescence dating
OSL dating was carried out at the OSL laboratory of the University of Freiburg, Germany. Bulk samples were first wet-sieved and grain-size fractions of 100–250 μm were selected, then treated with 20% HCl to remove carbonates, followed by 15% H2O2 to remove organic matter. The quartz fraction was isolated by two steps of density separation using heavy liquids with densities of 2.58 g/cm3 and 2.70 g/cm3, respectively. The samples were subsequently etched with 40% HF for 60 minutes to remove any remaining feldspar grains and the outer part of the quartz grains affected by external alpha irradiation. Then, all samples were treated with 1 mol/L HCl for 60 minutes to remove fluorites. The quartz grains were mounted on aluminum discs with silicone oil to make up aliquots for equivalent dose (De) measurement. Both 2-mm and 4-mm diameter aliquots were measured for all samples, because the OSL signal is partly rather dim and the larger aliquot size provided much better signal-to-noise ratios.
All OSL measurements were conducted on a Freiberg Instruments Lexsyg-Smart TL/OSL reader equipped with blue LEDs (458 ± 5 nm) and IR LEDs (850 ± 20 nm). Laboratory irradiation was carried out using Sr-90 beta sources mounted on the reader. Equivalent dose (De) was measured using the single-aliquot regenerative-dose protocol (SAR) procedure (Murray and Wintle, Reference Murray and Wintle2003). OSL decay curves decrease to background, indicating that the quartz OSL signal is dominated by the fast component (Fig. 3A). The De values of the investigated samples did not reach saturation (see Fig. 3A, inset).
Figure 3. OSL characteristics of quartz grains from 4-mm aliquots sample KM2-1. (A) OSL decay curves of the natural and regeneration dose of 19.14 Gy rapidly decrease to background, indicating that the quartz OSL signal shows a similar shape and is dominated by the fast component (inset: example growth curve fitted using a single saturation exponential function). (B) Preheat plateau and (C) dose recovery tests. (D) Example of a De distribution. Open diamonds = individual aliquot values; solid diamonds = mean aliquot values; open circles = each aliquot; error bars indicate ± 0.68 Gy to ca. 1.24 Gy.
Preheat plateau and dose recovery tests were conducted on aeolian sample KM2-1 before the De measurement. A preheat plateau test was carried out with six aliquots for preheat temperatures of 160, 180, 200, 220, and 240°C. The result of the preheat plateau test shows a De plateau at 200–240°C (Fig. 3B). Dose recovery ratios for the sample at different preheat temperatures range from 0.9 to 1.1, and decrease slightly with the temperature (Fig. 3C). Thus, we chose to preheat at 220°C for 10 seconds for all samples to avoid possible thermally unstable luminescence components potentially remaining after low-temperature preheating. Possible feldspar contamination of samples was checked at the end of each SAR cycle by one infrared stimulated luminescence (IRSL) measurement and no aliquot showed any measurable IR signal, indicating the purity of the quartz extracts.
Figure 3D displays the De distribution of sample KM2-1, which shows basically a normal distribution. While most samples have results that were similar to sample KM2-1, we found that some samples had relatively widely scattered and positively skewed distributions with high overdispersion values. This is likely caused by partial bleaching of the OSL signal in the different aliquots (Olley et al., Reference Olley, Pietsch and Roberts2004). The average De was calculated using the Central Age Model (CAM) and the Minimum Age Model (MAM) for all samples for comparison, and the overdispersion and distribution shape of De distributions were used for age model selection (see Supplemental Information for details).
The concentration of U and Th of each sample was measured using inductively coupled plasma mass spectrometry (ICP-MS, NexION-300D) and the content of K was determined by a wavelength dispersive X-ray fluorescence spectrometry (XRF, Axios-mAX) at the analytical laboratory of the Beijing Research Institute of Uranium Geology, China National Nuclear Corporation (CNNC). Dose rates were calculated assuming all samples had an average water content of 10 ± 5%. Cosmic ray contributions were calculated based on the geographical latitudes, altitudes, and present-day burial depths of each sample (Prescott and Hutton, Reference Prescott and Hutton1994). Quartz OSL ages were calculated using the LDAC calculator v1.0 (Liang and Forman, Reference Liang and Forman2019; results presented in Table 1).
Table 1. Summary of optically stimulated luminescence (OSL) data for 13 samples from the KM, YG and ASK sections in the eastern KMD (a mean water content of 10 ± 5% is assumed for all samples). Od. = Overdispersion.
Grain size
Grain size measurements were carried out using the methods of Lu and An (Reference Lu and An1998) with minor modifications at the Key Laboratory of Geoscience Big Data and Deep Resource of Zhejiang Province in the School of Earth Sciences, Zhejiang University, China. First, ~2 g of each bulk sample were heated at 70°C with 10 ml of 30% H2O2 and 10% HCl to remove organic matter and carbonates, respectively. Second, the samples were put in deionized water for 24 hours to remove HCI. The samples then were added to 10 ml of 0.05 mol/L (NaPO3)6 and ultrasonicated for 10 minutes to accelerate dispersion. Grain size distributions were measured using an OMEC LS-909E laser particle-size analyzer, which has a measurement range of 0.02–2000 μm with an accuracy that is better than 1%, and better than 1% variation in terms of reproducibility.
The grain size parameters of all samples were calculated following the method of Folk and Ward (Reference Folk and Ward1957). However, previous research in the KMD reported that it was difficult to interpret these grain size parameters because the grain sizes of the sediments often had multi-modal distributions (X.L. Liang et al., Reference Liang, Niu, Qu, Liu, Liu, Zhai and Niu2019; A.M. Liang et al., Reference Liang, Dong, Su, Qu, Zhang, Qian, Wu, Gao, Yang and Zhang2020). Thus, we used the end-member modeling analysis (EMMA) to distinguish the components of multi-modal grain size distributions. The EMMA was expected to provide valuable information on the provenance of sediments and characteristics of the paleoenvironment (Weltje, Reference Weltje1997; Dietze et al., Reference Dietze, Hartmann, Diekmann, IJmker, Lehmkuhl, Opitz, Stauch, Wünnemann and Borchers2012; Y. Li et al., Reference Li, Song, Fitzsimmons, Chen, Wang, Sun and Zhang2018; Zhang et al., Reference Zhang, Zhou, Wang, Song, Zhao, Xie, Russell and Chen2018; Liu et al., Reference Liu, Wang, Zhao and Yang2019), while multiple algorithm models have been developed for unmixing end-members (EMs) (Yu et al., Reference Yu, Colman and Li2016; Zhang et al., Reference Zhang, Zhou, Wang, Song, Zhao, Xie, Russell and Chen2018). We chose the AnalySize software (version 1.2.2) to unmix the grain size distribution data (Paterson and Heslop, Reference Paterson and Heslop2015). Unmixing was conducted through non-parametric analyses, with the best number of EMs being identified by the software. The optimal number of EMs needs to be most strongly correlated with the dataset (e.g., high R 2), having both the lowest correlation between EMs themselves and the lower angular deviation, and minimizing the number of EMs (Paterson and Heslop, Reference Paterson and Heslop2015; Duan et al., Reference Duan, An, Wang, Herzschuh, Zhang, Zhang, Liu, Zhao and Li2020).
Major, REE, and trace elements
Geochemical measurements were made on 31 samples from the KM section. A Retsch mixer mill MM-400 grinder (vibrational frequency at 20 Hz, 60 minutes grinding time) with an agate cup was used to grind the samples into a homogenized powder with diameters <75 μm (200 mesh according to Tyler standards). Major elements were determined using an XRF (Axios-mAX), with an analytical precision better than ± 3% of the value for all elements except for MnO (up to ± 10%). All major element concentrations (Al, Mg, Ca, Na, K, Ti) are expressed as weight percentages (wt. %), with iron expressed as Fe2O3. Loss on ignition (LOI) was determined by weighting pre-dried samples after heating at 1000°C for 1 hour. Trace elements (Zr, Rb, Ba), including rare earth elements (REE), were measured by inductively coupled plasma mass spectrometry (ICP-MS, NexION-300D). Rh and Re were used as the internal standards, and the standard W-2a of the United States Geological Survey (USGS) was used as the reference. The relative standard deviation is less than ± 1% and ± 5% for the REE and trace elements, respectively. Both major- and trace-element measurements were conducted at the analytical laboratory of the Beijing Research Institute of Uranium Geology, CNNC.
In this study we used the chemical proxy of alteration defined as CPA (= [Al2O3/(Al2O3 + Na2O)] × 100, molar contents) (Buggle et al., Reference Buggle, Glaser, Hambach, Gerasimenko and Marković2011) as the standard to estimate the weathering intensity of samples. We also used ratios of mobile to immobile elements such as (CaO + Na2O + MgO)/TiO2 as an index of the degree of weathering in order to reduce the grain size effect on chemical composition (S.L. Yang et al., Reference Yang, Ding and Ding2006) because these ratios vary little with grain size changes.
RESULTS
Stratigraphy
The KM section, which consists of three sub-sections (39.81°N, 93.45°E, 1554 m asl), is exposed in a dry valley on the west side of the Duobagou gully (Fig. 2). KM-1 is located at the top with a thickness of 1.7 m, facing NE (65°). KM-2 is located just below KM-1, with a thickness of 1.6 m and the same orientation as that of KM-1. At a horizontal distance of ~5.1 m below KM-2 is the KM-3 sub-section with a thickness of 1.75 m, facing NE (55°). The composite KM section consists of four layers of lacustrine sediments intercalated by aeolian sands and fluvial sediments. The sedimentary facies were identified in the field on the basis of texture and bedding. The stratigraphy of this section can be divided into five units upwards (Fig. 4A). Unit 1: 10.15–8.40 m, medium- to coarse-grained brown sand (10YR 5/3; Munsell Color, 2000) with inclusion of some small pebbles; the section did not reach the lower boundary of this unit (depositional interpretation = fluvial). Unit 2: 3.30–2.70 m, well-sorted, coarse-grained yellow-brown (10YR 5/4) sand with cross-bedding (depositional interpretation = aeolian). Unit 3: 2.70–2.06 m, alternations of clay and sand showing multiple facies changes with several shifts between aeolian and lacustrine deposition (depositional interpretation = shifting between aeolian and lacustrine processes). The upper 0.51 m (Unit 3b) is medium- to coarse-grained, yellow-brown (10YR 5/4) sand with wavy bedding (depositional interpretation = fluvial). The lower portion of the section at 2.70–2.57 m (Unit 3a) consists of fine-grained, yellow-brown (10YR 5/4) sand, and a 3-cm thick silty clay layer at the bottom (depositional interpretation = lacustrine). Unit 4: 2.06–0.50 m, one distinct gray-white (10YR 6/3) silty clay units occur within the poorly sorted medium- to coarse-grained yellow-brown (10YR 5/4) aeolian sands at 2.06–1.96 m (depositional interpretation = fluvial). Unit 5: 0.50–0 m, a strongly cemented gray-white (10YR 6/3) clay layer (depositional interpretation = lacustrine; Fig. 5).
Figure 4. (A) Lithology, OSL ages, and proxy indexes of mean grain size (Mz), EM components, LOI, CPA, and (CaO + Na2O + MgO)/TiO2 of the KM section. (B, C) Lithology, OSL ages, and mean grain size and EM components of the (B) YG section and (C) the ASK section. The light gray bars indicate possible humid events.
Figure 5. Field photo of the upper part of the KM section (units 4 and 5 in Fig. 4A) showing cemented lacustrine deposits in the upper 50 cm (12.5–17.5 on the ruler) and loose aeolian sand in the lower part.
The YG section (39.98°N, 93.99°E, 1196 m asl) is located on the west bank of Xitugou gully, ~1 km north of Yangguan Oasis. The section faces NE and is 3.65 m thick. The stratigraphy of this section can be divided into three units upwards (Fig. 4B). Unit 1: 3.65–1.80 m, three light gray-white (10YR 6/3) lacustrine silty clay units occurring at 3.65–3.55, 3.35–3.21 (strongly cemented), and 2.05–1.80 m, within the overwhelmingly brown (10YR 5/3) coarse-grained sands (depositional interpretation = fluvial). Unit 2: 1.80–1.65 m, well-sorted yellow-brown (10YR 5/4) sand without bedding (depositional interpretation = aeolian). Unit 3: 1.65–0 m, two white (10YR 6/3) silty clay layers occurring at 1.65–1.25 and 0.50–0.30 m, within the yellow-brown (10YR 5/4) coarse-grained sand (depositional interpretation = fluvial and lacustrine).
The ASK section is located at the southeastern margin of the KMD (39.69°N, 94.27°E, 1577 m asl). The section is 1.10 m thick, facing NE (25°). The sedimentary facies can be divided into two units (Fig. 4C). Unit 1: 1.10–0.45 m, yellow-brown (10YR 5/6) sand with oblique bedding, and with a dip angle of 15°, while at the depth of 0.75 m the color of the sand changes to gray-white (7.5YR 5/6) (depositional interpretation = aeolian). This unit belongs to typical aeolian sands. Unit 2: 0.45–0 m, gray-white (7.5YR 5/6) clay, with distinct horizontal bedding at a depth of 0.10 m (depositional interpretation = lacustrine deposits).
Chronology
Detailed discussion of which aliquot size and age model were chosen is given for each sample individually in the Supplementary Information. In general, while 4-mm aliquots may mask partial bleaching due to averaging effects (Wallinga, Reference Wallinga2002), the higher number of grains produces more robust results due to larger signal-to-noise ratios. The finally selected results are summarized in Table 1. The oldest sample is from the lower part of the KM section, dated to 16.8 ± 3.0 ka, while the youngest one is from the upper part of the ASK section, dated to 0.4 ± 0.1 ka.
Combined with stratigraphic units, the sedimentary history of the three sections can be reconstructed on the basis of these OSL ages. For the KM section, two OSL samples from the basal unit are dated to 16.8 ± 3.0 ka (sample KM3-1) and 15.2 ± 2.6 ka (KM3-2), suggesting that the KM-3 section probably formed during the last glacial maximum (LGM). We did not dig downwards to reach the lower boundary of this unit. The upper boundary of stratigraphical unit 2 yields an OSL age of 10.3 ± 1.2 ka (KM2-2). The lower part of the fluvially formed unit 4 is dated to 13.1 ± 1.9 ka (KM1-1), while the upper part is dated to 7.6 ± 0.8 ka (KM1-2). Thus, the overlying lacustrine unit should have been deposited after ca. 7 ka (Fig. 4A).
Chronology of the YG section is established by five OSL ages. The samples underlying and overlying the lacustrine layer at a depth of 3.35–3.21 m are dated to 4.4 ± 0.7 ka (YG1-1) and 2.2 ± 0.3 ka (YG1-3), respectively. Similarly, the samples from the sandy layer under- and overlying the lacustrine sediments at a depth of 2.05–1.80 m are dated to 2.3 ± 0.5 ka (YG1-4) and 2.4 ± 0.3 ka (YG2-1), respectively. Hence, the dating results place these lacustrine layers between ca. 4.4–2.2 ka. The upper part of sedimentary unit 3 yields an OSL age of 1.2 ± 0.4 ka (YG2-3), indicating that the uppermost lacustrine layer should have been deposited before ca. 1 ka (Fig. 4B).
For the ASK section, two OSL samples are dated to 0.7 ± 0.2 (ASK-1) and 0.4 ± 0.1 ka (ASK-2), indicating that the overlying lacustrine sediments were deposited not more than 400 years ago in this currently hyper-arid environment (Fig. 4C).
Grain size
There is a clear variation of grain size in the KM section (Fig. 4A), but due to instrument limitation the gravels in the samples were not measured. At the bottom of the section (i.e., stratigraphical unit 1) the sediments include fine-grained sands (mean diameter 140–170 μm) and gravels. The sediments then change to coarse-grained sands at a depth of 3.30–2.70 m (stratigraphical unit 2), with a mean grain size of 150–200 μm. The grain size changes frequently between 2.70 and 2.06 m (stratigraphical unit 3), from clayey silt at 2.67 m to coarse-grained sand at 2.10 m. The sediment consists of sand at a depth of 1.96–0.50 m, and sandy silt at 2.06–1.96 m. The upper 0.50 m of the section (stratigraphical unit 5) consists of clayey silt.
The grain sizes of the YG and ASK sections are shown in Figure 4B and 4C, respectively. For the YG section, the mean grain size is smaller in lacustrine facies than in fluvial deposits (gravels were not measured here, either). In unit 1, clayey silt occurs at depths of 3.65–3.55, 3.35–3.21, and 2.05–1.80 m. Unit 2 consists of aeolian sands with a mean grain size of 104–141 μm. Unit 3 consists of fluvial sands intercalated with clayey silts at depths of 1.65–1.25 and 0.50–0.30 m. Medium- to fine-grained sand occurs at a depth of 1.25–0.50 m, with a grain size range of 174–431 μm. Sandy silt dominates the uppermost 0.30 m.
In the ASK section, the mean grain size shows a clear trend of decreasing upwards. Aeolian sand with a mean diameter of 204–240 μm occurs in unit 1, and silty sediments with mean grain sizes of 36–161 μm occur in unit 2 (Fig. 4C).
For the end-member modeling analysis (EMMA), we calculated the correlation within the dataset (R 2) and between the internal end members and angular deviation to identify the optimal number of end members (see detailed discussion about the optimal number of end members selected summarized in Supplementary Information). Two EMs were identified for the KM section, with modal grain sizes of 44.8 μm (KM-EM1) and 183.3 μm (KM-EM2). Figure 4A shows that the values of KM-EM1 generally remain low, except for three abrupt increasing events. KM-EM2 shows an inverse trend compared to KM-EM1.
In the sediments of YG section, four EMs were identified with modal grain sizes of 17.5 μm (YG-EM1), 114.7 μm (YG-EM2), 183.4 μm (YG-EM3), and 370.9 μm (YG-EM4). The percentage of YG-EM 1 is generally low from ca. 3.9 ka to ca. 2.0 ka; two high values occurred after ca. 2 ka. In Figure 4B only the sum of YG-EM2 and YG-EM3 is shown, indicating an inverse trend of variation compared to YG-EM1. The percentage of YG-EM4 is mostly quite low, except for the middle portion of unit 3.
Two EMs were identified for the ASK section with modal grain sizes of 90.7 μm (ASK-EM1) and 231.9 μm (ASK-EM2). The percentage of ASK-EM 1 increases after ca. 0.3 ka while that of ASK-EM2 shows an opposite trend of changes (Fig. 4C).
Weathering intensity indexes and LOI
Geochemical measurements were done only for samples from the KM section, which is the thickest of the three sections in this study. The weathering intensity index CPA shows no significant change in unit 1 (mean value 73.9), but is characterized by two abrupt increases at the bottom and top of unit 2 (Fig. 4A). The (CaO + Na2O + MgO)/TiO2 ratio (mean value = 23.5) is generally large in the entire section, although frequent variations can be recognized (Fig. 4A).
The LOI value shows the opposite trend from the mean grain size and is generally low (~7%) in the KM section. The LOI of samples from unit 1 is ~7.5%. In unit 2, LOI varies between 5.6 and 6.8%. The LOI of unit 3 is generally around 6% and reaches 13.9% at the bottom of this unit (2.67 m). The LOI of samples from unit 4 varies between 5.9% and 7.7%, with a peak (8.5%) at a depth of 1.96 m. The values LOI of samples from unit 5 are ~9.8%.
DISCUSSION
Climatic implications of changes in paleoenvironmental proxies
Grain size and EMs
The grain-size distributions of unconsolidated sediments may provide information on changes in the source and transport of clastic materials and the paleoenvironment (Pye, Reference Pye1987; Ding et al., Reference Ding, Derbyshire, Yang, Yu, Xiong and Liu2002; Qin et al., Reference Qin, Cai and Liu2005; Yin et al., Reference Yin, Qin, Wu and Ning2009; Xiao et al., Reference Xiao, Fan, Zhou, Zhai, Wen and Qin2013). KM-EM1 and ASK-EM1 exhibited an asymmetric unimodal distribution pattern, with peak values of ~45 μm (coarse-grained silt) and 91 μm (very fine-grained sand), respectively (Fig. 6). YG-EM1 exhibited an asymmetric bimodal distribution pattern, with main peak values of ~18 μm (medium-grained silt) and secondary peak values of ~183 μm (fine-grained sand) (Fig. 6). Compared to the grain size distribution curves of the YA section (40°29′26″N, 93°14′28″E) in the southern Beishan Mountains, northeast of the KMD (Liang et al., Reference Liang, Niu, Qu, Liu, Liu, Zhai and Niu2019), the end members of KM-EM1, YG-EM1, and ASK-EM1 are likely associated with lake sediments, which is consistent with our field observations (Fig. 4).
Figure 6. Grain-size frequency distributions and end members for samples of the three sections. (A, D) KM section, (B, E) YG section, and (C, F) ASK section.
KM-EM2, YG-EM2, YG-EM3, and ASK-EM2 are characterized by narrow unimodal shapes with peak values of ~183 μm (fine-grained sand), 115 μm (very fine-grained sand), 184 μm (fine-grained sand), and 232 μm (fine-grained sand), respectively (Fig. 6D–F). Surface aeolian sand from the KMD exhibits a unimodal distribution, with a mode value range of 200–450 μm (Liang et al., Reference Liang, Niu, Qu, Liu, Liu, Zhai and Niu2019). Therefore, YG-EM3 and ASK-EM2 are characteristic of aeolian sediments, although they are slightly finer grained than the value reported by Liang et al. (Reference Liang, Niu, Qu, Liu, Liu, Zhai and Niu2019). YG-EM2 has a fine tail with coarse- and medium-sized silt, and the content of fine components is significantly lower than lake sediments. The fine tail of YG-EM2 could represent a suspension component (Pye, Reference Pye1987; Tsoar and Pye, Reference Tsoar and Pye1987). Thus, YG-EM2 is also likely characteristic of an aeolian sedimentary environment. KM-EM2 has the same grain-size frequency distributions as YG-EM3 and ASK-EM2. However, since the stratigraphy indicates that the KM section is a complex section that includes aeolian, fluvial, and lacustrine sediments, it is inappropriate to infer that KM-EM2 is characteristic of aeolian sediment. Considering that the gravels in the samples were not measured, KM-EM2 probably represents variations of both aeolian and fluvial sands and would not be characteristic of any particular sedimentary process.
YG-EM4 exhibits an asymmetric bimodal distribution pattern, with main peak values of 371 μm (medium-grained sand) and secondary peak values of 35 μm (coarse-grained silt) (Fig. 6E). Fluvial sediments in the eastern KMD exhibit discontinuous bimodal distributions, with main peak grain sizes of 150–400 μm and sub-peak grain sizes of 10–20 μm. Mode values of grain sizes from the coarse and fine components in sediment from the Tarim River are 120 and 10 μm, respectively (Li et al., Reference Li, Mu and Xu2012), and mode values of Niya River sediments range from 80 to 450 μm and from 5 to 12 μm (Yang et al., Reference Yang, Du, Liang, Zhang, Chen, Rioual, Zhang, Li and Wang2021), which is consistent with features of YG-EM4. Thus, YG-EM4 is most likely characteristic of a fluvial sedimentary environment.
Fluvial deposits in deserts represent flowing waters in the past, and lacustrine sediments indicate that the section, along with its surrounding areas, had once been covered by water. Thus, both fluvial and lacustrine processes indicate relatively wetter conditions in deserts (X.P. Yang et al., Reference Yang, Preusser and Radtke2006). Based on this interpretation, we infer that KM-EM1, YG-EM1, and ASK-EM1, along with YG-EM4, likely are linked to relatively wetter environments, while YG-EM2, YG-EM3, and ASK-EM2 probably indicate drier environments with stronger aeolian processes.
Chemical paleoenvironmental proxies
Two parameters that mainly influence the chemical composition of sediments are variations in source material and sink-source distance. The immobile elements, Al, Ti, and Zr, and their ratios reflect the weighted average composition of the source material. In order to evaluate variations in the composition of the source rock, we use an Al-Ti-Zr ratio diagram following the method proposed by Fralick and Kronberg (Reference Fralick and Kronberg1997). If the data exhibit a linear trend in the ratio plot rather than a scattered distribution, it can be inferred that the composition of the source rock is invariable. The Al-Ti-Zr ratio diagram of the KM section (Fig. 7A) shows that the data have a linear trend, indicating that the source was relatively constant. Therefore, we conclude that variations in the CPA and the (CaO+Na2O+MgO)/TiO2 ratio in the KM section are not caused by changes in the source materials.
Figure 7. (A) Plot of TiO2/Zr versus (Zr × 100,000)/Al2O3 for sediments from the KM section. (B) Plot of K/Ba versus K/Rb for sediments from the KM section and potential source region. Data for alluvial sediments in the West Altyn-Tagh and Beishan Mountains are from Liang et al. (Reference Liang, Dong, Su, Qu, Zhang, Qian, Wu, Gao, Yang and Zhang2020). Data for the fluvial sediments in the East Altyn-Tagh Mountains and lake sediments in Yangguan Oasis are from the final report of the State Scientific Survey Project of China (2017FY101000) (unpublished). The numbers in parentheses are the number of samples.
K-feldspar has a unique geochemical signature of K/Rb and K/Ba values because of the specific mineral crystallization histories from magma. In addition, unlike other K-bearing minerals (e.g., mica), K-feldspar can survive after abrasion during mechanical transport processes (Muhs, Reference Muhs2017). Therefore, we chose the K/Rb versus K/Ba ratio as a sediment-source indicator to determine the provenance of the KM section in the eastern KMD, because the previous mineral analysis results show that both aeolian sands and their possible source sediments contain K-feldspar (Xu et al., Reference Xu, Lu, Zhao, Wang, Su, Wang, Liu, Wang and Lu2011). The K/Rb versus K/Ba ratio diagram of the KM section (Fig. 7B) reveals that none of the sediment samples of the section really matched the samples from potential source regions. Geographically, the section is located near the fluvial fans in the forelands of the East Altyn-Tagh. However, fine-grained sediments could have been blown to the site of sampling by northly winds. Thus, we assume that the deposits of the KM section are a mixture of materials from the mountains located both north and south of the KMD (as indicated Fig. 7B). Hence, the sediments of the KM section are probably from nearby areas. With a relative short distance between the source and site of depositional, the provenance of the sediments would not have caused a major change to chemical weathering indexes in the KM section. Therefore, we conclude that the CPA and the (CaO + Na2O + MgO)/TiO2 ratio can be used to indicate the weathering intensity in the eastern KMD.
Although Liang et al. (Reference Liang, Dong, Su, Qu, Zhang, Qian, Wu, Gao, Yang and Zhang2020) inferred that the sorting process probably did not entirely overprint the provenance signatures of the sediment in the KMD, results from the Loess Plateau (S.L. Yang et al., Reference Yang, Preusser and Radtke2006) and the Taklamakan Desert (Jiang and Yang, Reference Jiang and Yang2019) indicate that elemental distribution of the sediment depends on grain size to a considerable extent. However, previous studies also have demonstrated that the CPA and the (CaO + Na2O + MgO)/TiO2 ratio are grain-size independent (S.L. Yang et al., Reference Yang, Preusser and Radtke2006; Chen et al., Reference Chen, Li, Dong, Yu, Zhang and Yu2021). Thus, we infer that the CPA and the (CaO + Na2O + MgO)/TiO2 ratio of bulk samples also could reflect the chemical weathering intensity in the eastern KMD.
Higher chemical weathering intensity values commonly indicate a warmer and more humid environment, while decreased weathering intensity reflects increased dust input, presumably from cool and dry environments (S.L. Yang et al., Reference Yang, Preusser and Radtke2006; Buggle et al., Reference Buggle, Glaser, Hambach, Gerasimenko and Marković2011; G.Q. Li et al., Reference Li, Jin, Wen, Zhao, Madsen, Liu, Wu and Chen2014). Therefore, the CPA and LOI values confirm that the eastern KMD had been constantly under a generally dry climate from the LGM to the Early Holocene, except for two slightly warm and humid but short-lasting intervals between 13 ka and 9 ka (Fig. 4A). A similar interpretation could be derived from the variations of the (CaO + Na2O + MgO)/TiO2 ratios, which show lower values while CPA and LOI are higher.
Paleoenvironmental changes since the late Pleistocene in the east part of KMD
On the basis of the stratigraphical facies and their chronology and paleoenvironmental proxies in the KM, YG, and ASK sections, we tentatively reconstructed four main paleoenvironment stages since the late Pleistocene in the eastern KMD.
In stage 1 (ca. 17–15 ka, 10.15–8.40 m depth in the KM section) (Fig. 4), the stratigraphical facies are dominated by fluvial sediments, which indicate an increase in effective moisture in the hyper-arid environment. Even though the CPA value remains at a low level, the (CaO + Na2O + MgO)/TiO2 ratio fluctuates slightly.
In stage 2 (ca. 13–7 ka, 2.70–0.40 m depth in the KM section), the interval from ca. 13 ka to 7 ka corresponds to sedimentary units 2–4 in the KM section. The percentages of KM-EM1 show three abrupt-increase events, and the same happens to LOI and CPA (Fig. 4A). These events suggest that the eastern KMD had increases in water availability and potentially also in temperature three times between 13 ka and 9 ka, although precise chronologies for these events still are not available.
In stage 3 (ca. 4.4–2.2 ka, 3.65–1.65 m depth in the YG section), the percentage of YG-EM2+3 is significantly high, while the YG-EM1 component is high only in two short intervals (Fig. 4B). We interpret this stage to have been generally wetter due to continuous fluvial sedimentation and the peaks of YG-EM1, which could reflect short-lived increases in water availability.
In stage 4 (ca. 2.2 ka–present, 1.65–0 m depth in the YG section and 1.1–0 m depth in the ASK section), the percentages of the YG-EM1 show two increase events, with the higher value occurring immediately after ca. 2.4 ka, followed by a strong fluvial process as shown by the increase in YG-EM4 component. ASK-EM1 reveals a lacustrine environment during the last 400 years, probably associated with an increase in effective moisture (Fig. 4C).
Although we aim to reconstruct paleoclimatic changes in the KMD since the LGM, the three sections of our study allow us to recognize only the characteristics of several short epochs during the last 17,000 years. It is obvious that three sections cannot be sufficient for reconstructing paleoenvironmental history in such a vast area. Because all three sections are located in the downstream parts of the desert gullies, no flooding can come that far during modern times. Thus, our interpretation of wetter environments should be viewed as typical evidence of paleoenvironmental changes.
All three sections that we studied contain fluvial and lacustrine facies indicating availability of water in this usually hyper-arid region. The fluvial sediments in the lower part of KM section suggest an increase in effective moisture between ca. 17 and 15 ka. A distinctly drier event occurred during ca. 13–7 ka, but probably not continuously during that interval because fluvial facies did not appear until ca. 10 ka. For the Late Holocene, the fluvial and lacustrine facies, which are dated 4.4–2.2 ka in the YG section and the last 400 years in the ASK section, indicate increases in moisture. Although many studies show a combination of increased moisture and temperature, the slight variation of CPA does not enable us to assess temperature changes.
Driving mechanisms of past environmental changes in the eastern KMD
Variation of the East Asian monsoon system has been regarded as one of the most important aspects that drives climate change in Chinese drylands (Zhou et al., Reference Zhou, Dodson, Head, Li, Hou, Lu, Donahue and Jull2002; Ding et al., Reference Ding, Derbyshire, Yang, Sun and Liu2005; Lancaster et al., Reference Lancaster, Wolfe, Thomas, Bristow, Bubenzer, Burrough and Duller2016). However, this does not apply for the KMD because the East Asian summer monsoon can barely reach our study sites. Evidence from fluvial and lacustrine sediments in the southern part of the Taklamakan Desert (mainly along the Niya River) indicates that the Late Holocene was the wettest period of the past 12 ka due to increase in runoff caused by increased melting of glaciers and potential increase of rainfall (Yang et al., Reference Yang, Du, Liang, Zhang, Chen, Rioual, Zhang, Li and Wang2021), which is similar to our reconstruction suggesting two more humid events—between 4.4 and 2.2 ka and during the last 400 years. Hence, paleoenvironmental change in the KMD probably reflects the same driving mechanism affecting the Taklamakan Desert (mid-latitude Westerlies, MLW). Ding et al. (Reference Ding, Tang, Su, Zhang and Lu2017) suggested that the sandy loess BL section started to form at ca. 8.3 ka, which is supposed to reflect the melting of glaciers on the northern slope of the Altyn-Tagh Mountains. Thus, we speculate that the increase in water availability due to melting glaciers in the Altyn-Tagh Mountains should be viewed as one of the main drivers causing environmental changes in the KMD. But convincing interpretation requires chronological data about glacial variations, which remain as a large knowledge gap.
Previous research in the arid central Asia (ACA) demonstrates that climate condition is mainly affected by the MLW (Li, Reference Li1990; Wu and Guo, Reference Wu and Guo2000; Yu et al., Reference Yu, Chen, Liu and Wang2001; Vandenberghe et al., Reference Vandenberghe, Renssen, van Huissteden, Nugteren, Konert, Lu, Dodonov and Buylaert2006; Chen et al., Reference Chen, Yu, Yang, Ito, Wang, Madsen and Huang2008). During the warm period, the MLW moved northward, and the ACA received less precipitation and a dry climatic condition prevailed. During the cold period, the MLW became stronger and moved southward, bringing more precipitation towards the ACA (Thompson et al., Reference Thompson, Mosley-Thompson, Davis, Lin, Dai, Bolzan and Yao1995; B. Yang et al., Reference Yang, Braeuning, Johnson and Yafeng2002; X. Yang et al., Reference Yang, Rost, Lehmkuhl, Zhu and Dodson2004; Chen et al., Reference Chen, Chen, Holmes, Boomer, Austin, Gates, Wang, Brooks and Zhang2010; Lauterbach et al., Reference Lauterbach, Witt, Plessen, Dulski, Prasad, Mingram and Gleixner2014).
Our paleoenvironmental reconstruction is generally consistent with several studies in the ACA. For example, loess-paleosol sequences and their magnetic susceptibility records from the northern slope of the Tianshan Mountains (Chen et al., Reference Chen, Jia, Chen, Li, Zhang, Xie, Xia, Huang and An2016) and Kyrgyzstan (Li et al., Reference Li, Song, Orozbaev, Dong, Li and Zhou2020) show a wetter climate during the Late Holocene (Fig. 8). Tang et al. (Reference Tang, Ding and Zhang2017) found multiple variations between dry-windy conditions and opposite wetter conditions during the Holocene, with some of the wetter periods reconstructed from the YG section (Fig. 8). Several teams of scientists have reconstructed Holocene climate changes on the basis of lacustrine sediments in Bosten Lake. Among them, Huang et al. (Reference Huang, Chen, Fan and Yang2009) proposed that the driest climatic conditions occurred between 17–8 ka, while a moderately dry period took place between 8–6 ka, based on pollen assemblages (A/C ratio and relative Ephedra abundances) and grain size characteristics in the sediments. In a different core from the same lake, grain size data of detrital sediments, C/N and C/S ratios of organic matter, as well as δ18O and δ13C values of lake carbonate indicate that Bosten Lake had a high salinity between 5.2 ka and 3.0 ka, but low salinity 3.0–2.2 ka (Zhang et al., Reference Zhang, Feng, Yang, Gou and Sun2010), meaning a later beginning of wetter conditions in the Late Holocene.
Figure 8. Holocene paleoenvironmental records of the KMD and adjacent areas reflected by different proxies. (A) Moisture index from core YKD0301 in Lop Lur (Liu et al., Reference Liu, Zhang, Jiao and Mischke2016). (B) Salinity index from core BSTC2000 in Bosten Lake (Zhang et al., Reference Zhang, Feng, Yang, Gou and Sun2010). (C) Moisture proxy ΔL* record from BSK loess section in the western Tian Shan, Kyrgyzstan (Li et al., Reference Li, Song, Orozbaev, Dong, Li and Zhou2020); the dash line represents loose soil, which is affected by modern plants. (D) Holocene moisture index in the core area of the Westerlies-dominated region, based on the magnetic susceptibility record of the LJW10 section (Chen et al., Reference Chen, Jia, Chen, Li, Zhang, Xie, Xia, Huang and An2016). (E) Percentage Ephedra abundances record from core BST04H in the Bosten Lake (Huang et al., Reference Huang, Chen, Fan and Yang2009); the dash line indicates no measurable pollen existed. (F) Aridity index in the KMD region, based on the record of the >110 μm grain size fraction in the BL section (Tang et al., Reference Tang, Ding and Zhang2017). (G) Moisture intensity and (H) human activity intensity in the eastern KMD (this study). The gray-shaded intervals represent the occurrences of relatively wet conditions since the Middle Holocene in these records.
Based on variations of grain size, pollen, ostracode species, and the composition of the soluble salts, Liu et al. (Reference Liu, Zhang, Jiao and Mischke2016) reported the occurrence of high salt content during 5.1–2.4 ka, suggesting saline conditions in Lop Nur and a possibly dry environment (Fig. 8). The difference between our results and those from Liu et al. (Reference Liu, Zhang, Jiao and Mischke2016) might reflect regional variation, but more importantly the difference may reflect the complexities between the open landscape where we investigated and the closed lake system of Lop Nur.
Influence of paleoenvironmental change on human activity in the eastern KMD
Although the KMD has a relatively harsh environment, human activity has a long history in its fringe areas, especially in the eastern KMD. The earliest known human remains in the eastern KMD are from the Xitugou site, which is located on the edge of the Yangguan Oasis. Pottery fragments and smelting slag indicate that the site could be dated to 3900–3400 yr BP (Liu et al., Reference Liu, Wang, Zhao and Ding2004; Y.X. Li et al., Reference Li, Chen, Qian, Chen and Wang2018). During the western Han Dynasty (206 BC through 25 AD), human activity in the eastern KMD increased rapidly, with large towns appearing in the oasis and several military facilities, such as beacon towers, built on their outskirts. These human settlements gradually expanded and became strategic strongholds along the ancient Silk Road (Zhai, Reference Zhai2017). Based on the absence of human activity remains in these historic sites after the mid-Tang Dynasty, Li (Reference Li1998) suggested that the towns of the eastern KMD probably have been abandoned since then.
Fluvial sediments and their grain size changes in the YG section show frequent fluvial/alluvial events in a probably more humid period starting at ca. 4.4 ka. Modern precipitation and hydrological data indicate occurrences of short-term rainstorms and heavy rainfalls in the eastern KMD and on the northern slope of the Altyn-Tagh Mountains from June to September (Hu et al., Reference Hu, Ning, Kang, Wang, Shang and Yang2017), and the runoff of rivers in the region occasionally increases abruptly between May and September (Chen et al., Reference Chen, Niu, Huang, Zhang, Wang and Zhang2017). Related fluvial/alluvial runoff transports large amounts of sediment from the Altyn-Tagh Mountains to the desert region. In light of the modern processes in this mountain–desert system (Fig. 9), we may assume that beyond human factors, the abandonment of the ancient Silk Road in the eastern KMD was probably caused by a combination of floods and desertification, because large floods could have destroyed the towns in the desert margins on one side, and potential desert encroachments due to abundant source material transported by the floods might have challenged human survival on the other side.
Figure 9. Schematic representation illustrating the possible influences of past environmental changes on human activity in the eastern KMD.
CONCLUSIONS
Although the KMD is one of the largest sand seas in China, field-based paleoenvironmental studies on this sand sea are still rare. Here, we presented a newly established chronologic framework using OSL dating for three aeolian sections with intercalated fluvial and lacustrine deposits in the eastern KMD. After careful assessment regarding the climatic implications of proxy data such as grain size, LOI, CPA, and (CaO + Na2O + MgO)/TiO2 ratios in these three sections, we used stratigraphy and its OSL chronology to reconstruct the history of paleoenvironmental changes in the eastern KMD since the late Pleistocene.
Considering the generally hyper-arid nature of the region, we conclude that the paleoenvironment was characterized by aridity in the eastern KMD despite some fluctuations in effective moisture that have occurred since the late Pleistocene. Several short-lived, more-humid periods or phases can be recognized in the sections we studied. Although our data do not allow us to draw a continuous curve of paleoclimatic variation, we are still able to conclude that changes in water availability and potentially also in temperature have occurred during the last 17 ka. Our data can confirm that in the interval ca. 17–15 ka, the eastern KMD experienced an increase in effective moisture. During the interval ca. 13–7 ka, the eastern KMD was likely still drier and windy but with one short-lived wetter event at ca. 9 ka. A possibly wetter period occurred ca. 4.4–2.2 ka, which is consistent with the occurrence of human activity in this region. We invoke melting of glaciers on the Altyn-Tagh Mountains as a significant contributor to paleoenvironmental change in the eastern KMD.
Supplementary material
The supplementary information for this article can be found online at https://doi.org/10.1017/qua.2023.38.
Acknowledgments
This study was jointly supported by the National Natural Science Foundation of China (grant No. 42230502), the Ministry of Science & Technology of China (grant No. 2017FY101001) and the China Scholarship Council (grant No. 202106320129). We are also very grateful to Dr. Mackenzie Day, an anonymous reviewer, associate editor John Dodson, and senior editor Lewis Owen for their helpful comments and suggestions on the initial draft of this paper.
Competing interests
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Open access
