INTRODUCTION
The last glacial termination is a global climate transition phase from the last glacial maximum (LGM) to full Holocene conditions. It may represent the largest readjustments of the Earth's natural climate over the past 100,000 yr (Broecker and van Donk, Reference Broecker and van Donk1970; Denton et al., Reference Denton, Anderson, Toggweiler, Edwards, Schaefer and Putnam2010). During this interval, a warming trend was generally interrupted by rapid millennium-scale cold reversals, such as the Greenland (Isotope) Stadial 1 (GS-1, Younger Dryas) and GS-2a (Heinrich 1) events (Björck et al., Reference Björck, Walker, Cwynar, Johnsen, Knudsen, Lowe and Wohlfarth1998; Rasmussen et al., Reference Rasmussen, Andersen, Svensson, Steffensen, Vinther, Clausen and Siggaard-Andersen2006). This may have induced pauses in glaciers retreating from their LGM limits and even glacial readvances, so the period is usually called the “late-glacial” period. Currently, it is far from clear how glaciers responded to these cold climate reversals. For example, we are not sure how many glacial stillstands and/or advances occurred during this period in different North Hemisphere (NH) regions. Constraining the timing of these glacial events can therefore help us recognize the processes contributing to the last glacial termination and understand the role of millennial-scale climate events in adjusting glacial extents.
Mountain glaciers are sensitive to regional and local climate change, and thus glacial records are usually used to retrieve past climatic change (Thompson et al., Reference Thompson, Yao, Davis, Henderson, Mosley-Thompson, Lin, Beer, Synal, Cole-Dai and Bolzan1997; Thackray et al. Reference Thackray, Owen and Yi2008). Former glacial extents, marked by lateral and terminal moraines along with erosional trimlines, reflect spatial and temporal changes of glacial mass balance that are mainly influenced by climatic conditions (Cuffey and Paterson, Reference Cuffey and Paterson2010). Therefore, paleoclimatic inferences can be made by reconstructing the timing and extent of past mountain glaciers.
As the most prominent topographic feature across the middle to low latitudes in Asia, the Tibetan Plateau (TP) rises to a mean elevation of more than 4000 m above sea level (m asl), with an area of 2.5 million km2 (Zhang et al., Reference Zhang, Li and Zheng2002). Its uplift has great impact on atmospheric circulations, initiating the Asian monsoons and splitting the westerlies of the NH into two branches (Prell and Kutzbach, Reference Prell and Kutzbach1992; Raymo and Ruddiman, Reference Raymo and Ruddiman1992; Benn and Owen, Reference Benn and Owen1998). Such geographic features make the plateau a distinctive connection to global and regional climate (Molnar and England, Reference Molnar and England1990; Benn and Owen, Reference Benn and Owen1998; Zheng et al., Reference Zheng, Xu and Shen2002). Understanding the timing and features of late-glacial events on the TP is therefore important, because it can yield information on climate evolution during the last glacial termination.
Throughout the TP, cosmogenic 10Be exposure dating has been widely used to define the timing of glaciations (e.g., Owen et al., Reference Owen, Caffee, Finkel and Seong2008; Chevalier et al., Reference Chevalier, Hilley, Tapponnier, Van Der Weord, Liu-Zeng, Finkel, Ryerson, Li and Liu2011; Heyman et al., Reference Heyman, Stroeven, Harbor and Caffee2011). A few examples of late-glacial moraines have been identified using the cosmogenic 10Be dating method (Owen and Dortch, Reference Owen and Dortch2014; Heyman, Reference Heyman2014). However, relative to the LGM and other glaciations, such as Marine Isotopic Stage 4 (MIS 4) and MIS 3 (e.g., Owen et al., Reference Owen, Caffee, Finkel and Seong2008; Dortch et al., Reference Dortch, Owen and Caffee2013; Murari et al., Reference Murari, Owen, Dortch, Caffee, Dietsch, Fuchs, Haneberg, Sharma and Townsend-Small2014), reliable chronological evidence of the millennium-scale glacial events in the late-glacial period is still scarce. Moreover, little work has been done to quantitatively reconstruct glacier–climate conditions for the late-glacial period on the TP.
Our study focuses on the dating and modeling of the late-glacial glacial events in the region of Qiongmu Gangri peak, which lies at the westernmost margin of the west Nyaiqentanggulha Mountains (Fig. 1). Situated in the transition zone that is dominated by the westerlies during the winter and spring and by the Indian summer monsoon during the summer and fall (Yao et al., Reference Yao, Masson-Delmotte, Gao, Yu, Yang, Risi and Sturm2013), the region is well suited for assessing glacial response to these two climate systems. The landforms of the Qiongmu Gangri peak afford a record of glacial activities since the last glacial (Dong et al., Reference Dong, Xu, Zhou, Fu, Zhang and Li2017). Here we present a 10Be exposure chronology of moraines that documents the timing of late-glacial glacial events in the Pagele valley on the eastern slope of the Qiongmu Gangri peak. In the valley, the main glacier terminates at an elevation of 5470 m asl and has a length of 2500 m. To infer the past climate that prevailed in the last glacial termination from the glacial geomorphological record, we use a coupled glacial mass-balance and ice-flow model to generate a glacier–climate reconstruction for the late-glacial glacial events in the Pagele valley.
Figure 1. Google Earth map showing the location of the Pagele valley, Qiongmu Gangri peak in the southern Tibetan Plateau (TP). The inset shows the main atmospheric circulation systems operating on the TP. Note that the PM1, PM2, and PM3 moraine crests are delineated by dark pink, green, and yellow lines, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
METHODS
Geomorphology and 10Be sampling method
In the Qiongmu Gangri peak region, the LGM end moraines are prominent and are located at tributary valley mouths with an elevation of ~4800 m asl. The LGM lateral moraines can be traced upvalley to an elevation of ~5100 m asl (Dong et al., Reference Dong, Xu, Zhou, Fu, Zhang and Li2017). Inside the LGM moraines to the modern glaciers, several post-LGM moraine loops and remnants can be found, but the numbers of these post-LGM moraines are very different between valleys. We chose the Pagele valley as a 10Be sampling target, because eight moraine loops can be clearly found above 5100 m asl in the valley. This offers great potential to collect samples to study the millennium-scale cold reversals in the last glacial termination. In this study, we focus on the three outermost moraine sets that were most likely formed by late-glacial glacial events in the valley (Figs. 2 and 3). We named the three moraine sets PM1, PM2, and PM3 from outer (older) to inner (younger). The lithology of the clasts on the three moraines is very similar and consists of schist, gneiss, and granitic diorite. Although the three moraine sets are incised by a river, the loop shapes can be traced in the field. The PM1, terminating at an elevation of 5170 m asl, is the most discontinuous moraine, mainly due to the longtime river incision. It is covered by a 5- to 10-cm-thick soil layer with thin grass-topped turf. Glacial boulders, with diameters of up to 1.5 m, protrude on the moraine surface. These boulders commonly present slight weathering characteristics with knobs and cavernous pits on their surfaces. Upvalley to about 500 m, there is the second moraine set of PM2, rising 10–15 m above the riverbed. The lateral-frontal moraine PM2 terminates at an altitude of 5200 m asl and can be traced to an elevation of 5250 m asl. Some patches of soil, along with sparse vegetation, have also developed on this moraine. The boulders on the moraine set are larger than those on PM1, with the largest having a diameter of ~3 m. PM3 has features similar to PM2 but terminates at a slightly higher elevation of 5230 m asl. The positions of the three moraines indicate that the glacier retreated 900–1000 m from the PM1 to PM3, and the glacier's length at the PM1 position was 3500 m longer than the modern glacier.
Figure 2. Photographic views for the PM1, PM2, and PM3 moraines. (A) Field photo for the PM1 moraine with one of the sampled boulders in the foreground, as viewed looking northeast on the moraine crest; (B) the PM2 moraine with boulders on it, as viewed looking north on the moraine crest; (C) the PM3 moraine with one of the sampled boulders in the left lower corner of the photo, as viewed looking northwest on the moraine; (D) field photo showing the relative position of the PM2 and PM3 moraines, as viewed looking northwest on the PM2 moraine crest; (E) the relative position photo for the PM1, PM2, and PM3, as viewed looking southeast on the PM3 moraine; (F) a typical boulder sampled for 10Be exposure dating on the PM1. Note that the PM1, PM2, and PM3 moraine crests are delineated by dark pink, green, and yellow lines on the photos, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Figure 3. (color online) Geomorphological mapping for the PM1, PM2, and PM3 moraine from Google Earth map and cosmogenic 10Be dating results for the three moraine sets in the Pagele valley. Note that the age values in italics are outliers in each group. The inset shows plots of probability density function (PDF) for the 10Be exposure ages on the PM1, PM2, and PM3 moraines. Individual ages are plotted as a normal distribution (PDF) using the exposure age and its internal uncertainty. The cumulative PDFs are made by summing individual-age PDFs on the same moraines. The uncertainty bar (±1σ) represents the mean of internal uncertainties of the individual ages from the same moraine.
About 0.5 kg of rock samples for the 10Be exposure dating was chiseled from the upper surfaces of quartz-rich boulders on the three moraine sets (Supplementary Figs. S1–S3). Sampling locations were targeted on the moraine crests to avoid the likelihood of burial/exhumation histories and slope instability. To reduce the possibility that the boulders and rock surfaces may have been covered with snow for significant periods (several months per year) or were previously covered with sediment, the highest parts of the largest boulders were chosen for sampling. Five samples from each of the moraine sets were sampled to check the reproducibility of the dating results and the possibility of the 10Be inheritance by prior exposure. Each sample location was measured using a handheld GPS, and the size of each boulder was also recorded. The surrounding conditions were photographed to document the position and geomorphic features for each sampling site. Cosmic-ray shielding by topography was calculated using the codes from Li (Reference Li2013) on the 30-m ASTER DEM. No corrections were made for erosion or snow coverage.
10Be exposure dating method
Sample processing and 10Be measurements were conducted at Xi'an Accelerator Mass Spectrometry (Xi'an-AMS) Center, Institute of Earth Environment, Chinese Academy of Sciences. Quartz purification, Be separation, and cathode preparation were carried out following the methods of Kohl and Nishiizumi (Reference Kohl and Nishiizumi1992) and Dortch et al. (Reference Dortch, Owen, Haneberg, Caffee, Dietsch and Kamp2009). The 10Be/9Be ratios were measured by AMS with normalization to the revised10Be ICN standard with 10Be/9Be ratio of 2.851 × 10−12 (Nishiizumi et al., Reference Nishiizumi, Imamura, Caffee, Southon, Finkel and McAninch2007).
For each 10Be sampling set from the same moraines, one blank sample was used to correct the measured isotope ratios, which were then converted to 10Be concentrations for age calculations in the next step. We calculated the ages for three scaling models using CRONUS-Earth online calculators (version 3; Balco et al., Reference Balco, Stone, Lifton and Dunai2008; http://hess.ess.washington.edu, December, 2018). The sample ages are presented in Table 1, with external (analytical and production rate uncertainty) and internal (analytical uncertainty only) uncertainties at 1σ. Based on analytical approximations to modeled fluxes of main atmospheric cosmic-ray particles responsible for in situ cosmogenic nuclide production, the LSDn model by Lifton et al. (Reference Lifton, Sato and Dunai2014) predicts realistic atmospheric cosmic-ray fluxes. It works well at low-latitude and high-altitude sites like our study area. Also, the model has fewer systematic age uncertainties than the St and Lm (time-independent and time-dependent) models of Lal (Reference Lal1991) and Stone (Reference Stone2000). We therefore chose to use the exposure ages from the LSDn scaling model when discussing chronology in the paper. Furthermore, to statistically test whether there are outliers in an age group from one moraine, this study uses Peirce's criterion (Peirce, Reference Peirce1852) following the procedure of Blomdin et al. (Reference Blomdin, Stroeven, Harbor, Lifton, Heyman, Gribenski and Petrakov2016).
Table 1. Cosmogenic 10Be exposure dating results on the moraines of PM1 (NQPGL15-11 to NQPGL15-15), PM2 (NQPGL15-16 to NQPGL15-20), and PM3 (NQPGL15-21 to NQPGL15-25), with the age values in italics being outliers in each group.
a Time-independent production scaling model by Lal (Reference Lal1991) and Stone (Reference Stone2000).
b Time-dependent production scaling model by Lal (Reference Lal1991) and Stone (Reference Stone2000).
c Production scaling model by Lifton et al. (Reference Lifton, Sato and Dunai2014).
Glacial modeling method
This study used a glacial model to reconstruct the ice thickness and to invert the paleoclimate conditions from the three sets of moraines. The model couples a two-dimensional glacial mass-balance model to a shallow ice-approximation ice-flow model, which is fully described by Plummer and Phillips (Reference Plummer and Phillips2003). This model takes into account the topographic and atmospheric controls on energy and mass of the glacial surface, including explicitly incoming solar radiation, air humidity, wind speed, cloudiness, and snow avalanching on steep hillslopes. It therefore allows us to explicitly quantify the topographic shading and consider the impact of albedo on ice melting. We used the model to calculate the glacial mass balance at a monthly time step to get annual net mass balance. The calculated net mass balances were then input to the shallow ice-approximation ice-flow model, which was run iteratively until the calculated ice thickness attained a steady state.
The model domain was defined by a digital elevation model (DEM) with 30 m × 30 m grids, which was downloaded from the Geospatial Data Cloud (http://www.gscloud.cn, July, 2018). Modern glacial thickness data were adopted from Farinotti et al. (Reference Farinotti, Huss, Fürst, Landmann, Machguth, Maussion and Pandit2019) and subtracted from the DEM to generate the ice-free topography on which the model ran. The model was forced by monthly mean meteorological data over the 1981–2016 period from the Dangxiong Station (30°29′N, 91°06′E, 4200 m asl), the long-term observational station nearest to the Pagele valley. To scale the monthly temperature and precipitation with elevation, we used the regional values of temperature lapse rates and precipitation gradients from Xu et al. (Reference Xu, Pan, Dong, Yi and Glasser2017a). Monthly mean climate variables from the Dangxiong Station, as a reference for the glacier–climate modeling, are listed in Supplementary Table S1. Due to the large uncertainty in albedo and its dominant impact on energy, and thus net mass ablation calculation, we calibrated this parameter in an effort to match simulated glaciers with the observed modern glaciers in the Pagele valley. After multiple simulations with different albedo values, we confirmed that when the low and high albedos were 0.2 and 0.7, respectively, the simulated glaciers were well matched with the observed ones (Fig. 4). We therefore used the two values as optima in the following simulations. A full description of the model inputs is given in Supplementary Table S2. By altering temperature and precipitation from contemporary conditions (hereafter referred to as ΔT, temperature change relative to present; F p, precipitation as fractional value relative to modern), we ran the model iteratively until the simulated glaciers fit well with their corresponding moraine positions.
Figure 4. (color online) Comparison of the glacial distribution for observed (A) and modeled (B) glaciers under modern climate conditions in the Pagele valley.
RESULTS AND DISCUSSION
10Be exposure dating
In total we dated 15 boulder samples for the three moraine sets of PM1, PM2, and PM3. All of these 10Be ages are listed in Table 1 and Figure 3. Five samples from the PM1 moraine (NQPGL15-11 to NQPGL15-15) present ages ranging from 15,640 to 16,190 yr with an uncertainty-weighted mean age of 15,850 ± 980 yr. These ages cluster tightly and overlap with each other within 1σ internal uncertainty. The PM2 samples have five ages (NQPGL15-16 to NQPGL15-20) that constitute a wide spread from 13,790 to 31,740 yr. The two ages of 20,740 and 31,740 yr (NQPGL15-17 and NQPGL15-18) are identified as outliers by the Peirce's test, and they are too old to represent the timing when the PM2 moraine was abandoned by the Pagele glacier. The two older samples can be explained by 10Be inheritance (exposed before they were deposited on the moraine PM2). Excluding the two oldest ages (as outliers), three ages (NQPGL15-16, NQPGL15-18, and NQPGL15-20), however, are statistically indistinguishable within 1σ external uncertainty and have an uncertainty-weighted mean age of 14,140 ± 880 yr. For the PM3 moraine, five 10Be ages (NQPGL15-21 to NQPGL15-25) range from 1490 to 12,980 yr. The youngest age of 1490 yr (NQPGL15-24) is obviously an outlier of this group of ages (the boulder was most likely exhumed) and can be excluded when interpreting the age of the PM3 moraine. The remaining four ages overlap with each other within 1σ external uncertainty with an uncertainty-weighted mean age of 12,430 ± 790 yr. Collectively, excluding the outliers, the PM1, PM2, and PM3 ages are statistically distinguishable from each other based on the age probability peaks (fig. 3), and the three groups of ages are consistent with the morphostratigraphic features. This means we can confidently constrain the timing of the late-glacial glacial events using these ages from the respective PM1, PM2, and PM3 moraines.
Two scenarios have been generally applied to interpret the multiple cosmogenic 10Be ages from one moraine. One uses the mean of the multiple boulder 10Be ages to represent the timing when moraines were abandoned by glaciers (e.g., Schaefer et al., Reference Schaefer, Denton, Kaplan, Putnam, Finkel, Barrell and Andersen2009; Young et al., Reference Young, Briner, Schaefer, Zimmerman and Finkel2019). This scenario assumes that the moraine surfaces remained stable after they formed and that inheritance of 10Be is as important as postdepositional shielding in its effect on sample ages (Chevalier et al., Reference Chevalier, Hilley, Tapponnier, Van Der Weord, Liu-Zeng, Finkel, Ryerson, Li and Liu2011). The other scenario considers the oldest ages as the timing when moraine surfaces were exposed (e.g., Zech et al., Reference Zech, Zech, Kubik, Kharki and Zech2009; Heyman, Reference Heyman2014). This scenario assumes that the dated boulders are likely to be exhumed to moraine surfaces and/or rotated by postdepositional movement, considering the fact that moraine surfaces are unlikely to remain stable after their formation (Heyman et al., Reference Heyman, Stroeven, Harbor and Caffee2011). However, our 10Be ages from each of the three moraines cluster tightly at most within a ~1000 yr range after the 10Be age outliers are excluded. This allows us to correlate the moraines to the millennial-scale climate events during the last glacial termination. Therefore, regardless of the mean (15,850 ± 980, 14,140 ± 880, and 12,430 ± 790 yr) and the oldest ages (16,190 ± 1000, 14,710 ± 910, and 12,980 ± 820 yr) from each of the three moraines, they all suggest that the PM1, PM2, and PM3 moraines respectively correspond well to the GS-2a (or Heinrich 1, 16.9–14.7 ka), GI-1 (or Bølling-Allerød, 14.7–12.7 ka), and GS-1 (or Younger Dryas, 12.7/12.9–11.5/11.7 ka) events in the Greenland ice-core record (Björck et al., Reference Björck, Walker, Cwynar, Johnsen, Knudsen, Lowe and Wohlfarth1998; Rasmussen et al., Reference Rasmussen, Andersen, Svensson, Steffensen, Vinther, Clausen and Siggaard-Andersen2006; Ressen and Isarin, Reference Ressen and Isarin2001).
GS-2a and GS-1 glacial events have also been reported in other valleys of the west Nyaiqentanggulha Mountains. At the eastern margin of the west Nyaiqentanggulha Mountains, Owen et al. (Reference Owen, Finkel, Barnard, Haizhou, Asahi, Caffee and Derbyshire2005) reported a late-glacial moraine that is a latero-frontal moraine and terminates at an elevation of ~4920 m asl, 3300 m away from the modern glacier. They collected four boulder samples for cosmogenic 10Be dating from the moraine crest and dated the moraine back to 15,900 ± 900 yr (mean age). Using the same method in our study, we recalculated the ages for the four samples (Supplementary Table S3). The recalculated ages range from 15,540 ± 1050 to 17,900 ± 1220 yr (±1σ external uncertainty), with an uncertainty-weighted mean age of 16,640 ± 1110 yr using the LSDn scaling model of Lifton et al. (Reference Lifton, Sato and Dunai2014). Closer to our study site, Chevalier et al. (Reference Chevalier, Hilley, Tapponnier, Van Der Weord, Liu-Zeng, Finkel, Ryerson, Li and Liu2011) also dated a late-glacial moraine set (YanBaJian inner moraine) that is located at an elevation of ~5300 m asl, with a distance about 3000 m from the modern glacier. Five boulder 10Be samples dated the YanBaJian inner moraine back to 11,000 ± 2000 yr after one outlier was rejected (Chevalier et al., Reference Chevalier, Hilley, Tapponnier, Van Der Weord, Liu-Zeng, Finkel, Ryerson, Li and Liu2011). We also recalculated the ages for these five samples and found that the ages are between 11,850 ± 740 and 21,300 ± 1320 yr using our scaling assumptions. After outlier evaluation by Peirce's criterion (Blomdin et al., Reference Blomdin, Stroeven, Harbor, Lifton, Heyman, Gribenski and Petrakov2016), we identified two outliers (15,350 ± 960 and 21,300 ± 1320 yr). Rejecting the outliers, the remaining ages have an uncertainty-weighted mean age of 12,830 ± 810. These recalculated moraine 10Be ages are comparable to our 10Be ages for the PM1 and PM3 and thus suggest that the GS-2a and GS-1 glacial events occurred in the west Nyaiqentanggulha Mountains. We note that this is the first time that the chronology for the GI-1 glacial event in the region has been reported, making it impossible to correlate the corresponding ages to other sites.
Glacial modeling
To establish the climate scenarios that support the glacial geometries constrained by the PM1, PM2, and PM3 moraines, we set the F p to be from 0.6 to 1.4 with an incremental change of 0.1, and found the corresponding ΔTs that can force the model to produce the respective glacial extents (Fig. 5). Without precipitation change, the respective glacial extents for PM1, PM2, and PM3 required temperatures 1.9°C, 1.6°C, and 1.2°C lower than present, meaning the glacier is more sensitive to the temperature in wetter conditions than in drier conditions. For example, under precipitation of 30% more than the present (F p = 1.3), 0.5°C of warming is needed to simulate the glacial retreat from the PM1 to PM3; but under the precipitation of 30% less than the present (F p = 0.7), 1.0°C of warming is needed to simulate the same retreat. For each glacial extent, varying the precipitation by 10% in the wetter conditions can be compensated by a 0.1°C–0.2°C change in temperature; but in the drier conditions, a 0.3°C–0.4°C temperature change is required.
Figure 5. Plots of the temperature and precipitation combinations (ΔT–F p) that yield the GS-1, GI-1, and GS-2a glacial extents, which match the PM3, PM2, and PM1 moraine positions, respectively.
Figure 5 displays all ΔT–F p climate scenarios that can force the model to produce each of glacial extents constrained by the PM1, PM2, and PM3 moraines. Under these scenarios, the reproduced glaciers have areas of about 12.6, 11.6, and 10.9 km2 for the PM1, PM2, and PM3 moraines, respectively. For the PM1 glacial extent, however, the modeled glacial volumes range from 1.03 to 1.14 km3, and the maximum ice thicknesses vary from 397 to 424 m, with larger volumes and ice thicknesses in the wetter conditions. The glacial equilibrium line altitudes (ELAs) under these ΔT–F p conditions for PM1 are between ~5400 and ~5440 m, with lower values under cooler conditions. The successful ΔT–F p scenarios reproduce the glacial volumes varying from 0.79 to 0.88 km3 for the PM2 moraine, and from 0.61 to 0.68 km3 for the PM3 moraine. The modeled maximum ice thicknesses range from 386 to 402 m and from 376 to 387 m for the PM2 and PM3 moraines, respectively. The modeled ELAs for PM2 and PM3 vary from ~5460 to ~5490 m and from ~5510 to ~5540 m, respectively. Figure 6 demonstrates the modeled net annual mass balance and corresponding ice thickness for each of the three glacial events under three ΔT–F p combinations of −2.9°C–0.6, −1.6°C–1.0, and −1.5°C–0.8, respectively, and the corresponding ELAs are ~5420, ~5480, and ~5530 m under these three climate scenarios.
Figure 6. (color online) Simulations of the GS-1, GI-1, and GS-2a glacial mass balances and extents in the Pagele valley under the ΔT–F p combinations of −1.5°C–0.8 (A, D), −1.6°C–1.0 (B, E), and −2.9°C–0.6 (C, F), respectively. Note that the GS-1, GI-1, and GS-2a glacial limits (PM3, PM2, and PM1 moraine positions) are also shown for comparison.
When interpreting the paleoclimates based on the model results, we assume that the ΔT–F p combinations in Figure 5 are the only possible scenarios for the three glacial events. Uncertainties exist in these combinations due to potential uncertainties in the model input climate data and biases in the calculations. These include the altitudinal gradients of monthly mean temperature and precipitation, which are based on a small number of weather stations, and estimates of second-order climate data that are not available locally (e.g., wind speed, relative humidity). In addition, our lack of knowledge of how much these second-order climate variables differed during the last glacial termination also adds uncertainty to the modeled ΔT–F p combinations. Although we do not evaluate the uncertainties caused by those effects, previous studies have suggested that errors due to limited input data produce uncertainties in the ΔT and F p of about ± 0.5°C and ± 0.3, respectively (e.g., Plummer, Reference Plummer2002; Laabs et al., Reference Laabs, Plummer and Mickelson2006; Xu et al., Reference Xu, Hu and Qiao2013). As this investigation used a method similar to these studies, we believe that the uncertainties in our estimated ΔT and F p are not beyond these values. In addition, our study did not account for impacts of debris cover on the paleoclimate inferences due to lack of measurements of debris cover on the modern glacier. A thin supraglacial debris cover enhances glacial melting by reducing surface albedo and increasing the absorption of the incident heat, whereas a thick debris cover can protect the glacier from melting by insulating the surface. Which impact is dominant depends on a critical thickness (1–10 cm) of the debris cover (e.g., Lana et al. Reference Lana, Nakawo, Fukushima and Ageta1997; Pu et al. Reference Pu, Yao and Duan2003; Nicholson and Benn, Reference Nicholson and Benn2006; Mihalcea et al., Reference Mihalcea, Mayer, Diolaiuti, D'Agata, Smiraglia, Lambrecht, Vuiller-moz and Tartari2008). Considering the debris cover promotes glacial ablation, the inferred ΔT (temperature reduction in our study) would be overestimated by the glacial model (using a lower albedo value). In contrast, the inferred ΔT would be underestimated if the debris cover was thicker than the critical thickness. Furthermore, there are no data to quantify the highly variable distribution of debris on the glacier's surface and how this differs between modern and past glaciers. This makes it difficult to precisely evaluate the impact of debris cover on the paleoclimate inferences.
For the three glacial events, no unique paleoclimate solution can be found in our model experiment results (three ΔT–F p combination curves). The estimates of ΔT–F p combinations shown in Figure 5 are the only possible climate scenarios. To better constrain the ranges of the ΔT and F p, we have to reference other independent climate proxies. Oxygen isotopic ratios (denoted as δ18O) have been used extensively in reconstructions of long-term averaged air temperature (ice-core δ18O records) and Asian monsoon intensity (cave speleothem δ18O records), though several complicated moisture processes affect the seasonal variability in δ18O values in ice cores and cave speleothems (Dayem et al., Reference Dayem, Molnar, Battisti and Roe2010; Yao et al., Reference Yao, Masson-Delmotte, Gao, Yu, Yang, Risi and Sturm2013). Here we use the summer insolation at 30°N (Berger and Loutre, Reference Berger and Loutre1991) and δ18O records from the Guliya ice core on the TP (Thompson et al., Reference Thompson, Yao, Davis, Henderson, Mosley-Thompson, Lin, Beer, Synal, Cole-Dai and Bolzan1997) as proxies for the temperature change, and the speleothem δ18O Indian summer monsoon records from the Dongge, Yamen, Xiaobailong, Timta, and Bitto caves (Yuan et al., Reference Yuan, Cheng, Edwards, Dykoski, Kelly, Zhang and Qing2004; Sinha et al., Reference Sinha, Cannariato, Stott, Li, You, Chen, Edwards and Singh2005; Yang et al., Reference Yang, Yuan, Cheng, Zhang, Qin, Lin, Zhu and Edwards2010; Cai et al., Reference Cai, Fung, Edwards, An, Cheng, Lee and Tan2015; Kathayat et al., Reference Kathayat, Cheng, Sinha, Spötl, Edwards, Zhang and Li2016) as proxies for precipitation variability in the region (Fig. 7). The proxy chronologies are defined by using the dating methods of 14C, 36Cl, 230Th, and orbital tuning, and the proxy errors are typically less than 5%. More information about the quality and associated errors of each proxy is provided in the original publications. In Figure 7, the temperature proxy records show that the temperature was progressively increasing from GS-2a to GI-1, and then dropped in GS-1. The speleothem δ18O records indicate that the precipitation from the Indian summer monsoon was increasing from GS-2a to GI-1, and then decreasing in GS-1. Moreover, in the speleothem Indian summer monsoon records, the average δ18O value in GI-1 approximates to the past 1000 year δ18O value, and the values in GS-2a and GS-1 were higher than the past 1000 year δ18O value (Supplementary Table S4). Assuming a linear relationship between speleothem δ18O and precipitation, it can be deduced that the precipitation values in GS-2a, GI-1, and GS-1 were 60%–70%, 100%, and 80%–90%, respectively, of the value for the past 1000 years. From the modeling results, these precipitation reductions require temperature drops of 2.6°C–2.9°C, ~1.6°C, and 1.4°C–1.5°C to reproduce the PM1, PM2, and PM3 glacial extents, respectively. This suggests that a warming trend from the PM1 to PM3 glacial events could have dominated the glacier's receding during the last glacial termination, with the corresponding ELAs lower by ~300, ~250, and ~190 m relative to the modern ELA of ~5730 m. We acknowledge that a more specific climate inference is difficult until more accurate temperature or precipitation proxies are available for the three glacial events. However, in the Samdainkangsang Peak region, Xu et al. (Reference Xu, Dong, Pan, Hu, Bi, Liu and Yi2017b) also suggested that a temperature drop of 2.6°C–2.8°C and 60%–70% of modern (1981–2010) precipitation can support the late-glacial (GS-2a) glacial extents in the Barenduo valley, based on reconstructions of the late-glacial glacial extents in the valley. Due to lack of other quantitative climate constraints for the GI-1 and GS-1 glacial events, it is impossible to make a full comparison of the modeled climate scenarios in the region.
Figure 7. Time series of summer insolation at 30°N (Berger and Loutre, Reference Berger and Loutre1991) and δ18O record from Guliya ice core on the northwestern Tibetan Plateau (Thompson et al., Reference Thompson, Yao, Davis, Henderson, Mosley-Thompson, Lin, Beer, Synal, Cole-Dai and Bolzan1997) as proxies for temperature change, and Indian summer monsoon speleothem δ18O records from Dongge, Yamen, Xiaobailong, Timta, and Bitto caves (Yuan et al., Reference Yuan, Cheng, Edwards, Dykoski, Kelly, Zhang and Qing2004; Sinha et al., Reference Sinha, Cannariato, Stott, Li, You, Chen, Edwards and Singh2005; Yang et al., Reference Yang, Yuan, Cheng, Zhang, Qin, Lin, Zhu and Edwards2010; Cai et al., Reference Cai, Fung, Edwards, An, Cheng, Lee and Tan2015; Kathayat et al., Reference Kathayat, Cheng, Sinha, Spötl, Edwards, Zhang and Li2016) as proxies for precipitation variability in the region. The yellow and orange rectangles denote the GS-1 (12.9–11.7 ka) and GS-2a (16.9–14.7 ka) glacial stages, respectively, between which is the GI-1 period (14.6–13.0 ka; Björck et al., Reference Björck, Walker, Cwynar, Johnsen, Knudsen, Lowe and Wohlfarth1998; Rasmussen et al., Reference Rasmussen, Andersen, Svensson, Steffensen, Vinther, Clausen and Siggaard-Andersen2006). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Climatic relationships with oceanic and atmospheric circulations
During the last glacial termination, meltwater and iceberg outbursts from the margins of the NH ice sheets led to the glacial stadials of GS-2a (Heinrich 1) and GS-1 (Younger Dryas) in the North Atlantic region by reducing the thermohaline circulation (THC). A substantial drop in sea-surface temperature continued during the GS-2a and GS-1 periods in the region (Naughton et al., Reference Naughton, Sánchez Goñi, Kageyama, Bard, Duprat, Cortijo and Desprat2009). This resulted in an expansion of winter sea ice and introduced a highly seasonal climate in the region (Denton et al., Reference Denton, Alley, Comer and Broecker2005), because the spread of winter sea ice created Siberia-like conditions, dominated by large temperature differences between winter and summer in the North Atlantic region. This explains why the mean annual temperatures in Greenland and northern Europe decreased by 12°C to 17°C relative to today's values in GS-1, with a 22°C to 28°C decrease in winter and a 3°C to 6°C decrease in summer (Denton et al., Reference Denton, Alley, Comer and Broecker2005, Reference Denton, Anderson, Toggweiler, Edwards, Schaefer and Putnam2010). Although the HTC also controlled the Indian Ocean (where the Indian summer monsoon transports water vapor to the TP) for a major portion of the last glacial termination (Naik et al., Reference Naik, Basak, Goldstein, Naidu and Naik2019; Sun et al., Reference Sun, Zhang, Shulmeister, Bird, Chang and Shen2019), our estimated GS-1 temperature depression of 1.4°C–1.5°C for the Qiongmu Gangri peak is much less than the values from the North Atlantic region. This argument is in agreement with the view of Rehfeld et al. (Reference Rehfeld, Münch, Ho and Laepple2018) that the temperature decrease from the LGM to Holocene had a clear zonal pattern, with a progressive reduction in change from the high latitudes toward the tropics.
Using a global atmospheric circulation model, Barnett et al. (Reference Barnett, Dümenil, Schlese and Roeckner1988) related weakened Asian monsoons to long and cold high-latitude winters, similar to a consequence of the winter sea-ice cover in the North Atlantic. Chiang and Bitz (Reference Chiang and Bitz2005), using a Community Climate Model (version 3), also illustrated that the Intertropical Convergence Zone (ITCZ) shifted southward with an imposed increased ice cover anomaly like that of the North Atlantic during GS-2a. Model results therefore suggest that the expansion of sea ice across the North Atlantic, particularly in winter, was probably the key factor in spreading the impacts of the millennial-scale cold events throughout the NH and into midlatitudinal regions (Denton et al., Reference Denton, Anderson, Toggweiler, Edwards, Schaefer and Putnam2010) such as the TP. In addition, the cooling in the North Atlantic region during GS-2a might have increased the latitudinal temperature gradient in the NH. This then increased the meridional atmospheric pressure gradient and thus strengthened the westerlies, because the westerlies are formed by the air movement transition from meridional to zonal under the Coriolis force. The strengthened westerlies could have brought cold air into the TP. Consequently, during GS-2a, the cold westerlies and weak Indian summer monsoon dominated the TP, and thus the Qiongmu Gangri peak region was affected by cool and dry conditions. Moreover, by analyzing temperature proxy data from the Arabian Sea and climate model simulations, Tierney et al. (Reference Tierney, Pausata and de Menocal2016) suggested that during GS-2a, the sea-surface cooling in the Indian Ocean was a critical link between the North Atlantic event and the Indian monsoon failure.
The warmer and wetter conditions of the GI-1 period can be attributed to the enhanced THC, with the ITCZ shifting northward and thus the Indian summer monsoon being strengthened (Sinha et al., Reference Sinha, Cannariato, Stott, Li, You, Chen, Edwards and Singh2005; Liu et al., Reference Liu, Otto-Bliesner, He, Brady, Tomas, Clark and Carlson2009; Banakar et al., Reference Banakar, Baidya, Piotrowski and Shankar2017). Despite the relatively warmer conditions in the GI-1 period, it still appeared to be colder relative to the GS-1 period in the Qiongmu Gangri peak region, as evidenced by our glacier–climate modeling on the basis of the greater glacial extent of PM2 than PM3. This warming trend from GS-2a through GI-1 to GS-1, especially the less marked GS-1 signature, is also similar to the modeled temperature trend (Liu et al., Reference Liu, Otto-Bliesner, He, Brady, Tomas, Clark and Carlson2009; Tierney et al., Reference Tierney, Pausata and de Menocal2016) and the chironomid-based temperature reconstruction in southwestern China (Zhang et al., Reference Zhang, Chang, Shulmeister, Langdon, Sun, Cao, Yang and Shen2019). This pattern follows more closely the long increase in the NH summer insolation during the last glacial termination (Fig. 7), suggesting that the influence of the North Atlantic millennial-scale cold events on the TP declined as the last glacial termination progressed. Zhang et al. (Reference Zhang, Chang, Shulmeister, Langdon, Sun, Cao, Yang and Shen2019) also proposed that during the last glacial termination, there was a progressive increase in Indian summer monsoon influence on climate change in southwestern China. Furthermore, using a pollen discrimination index, Zhu et al. (Reference Zhu, Lü, Wang, Peng, Kasper, Daut and Haberzettl2015) suggested that the Lake Nam Co area, ~70 km northeast of the Qiongmu Gangri peak, was influenced by the Indian summer monsoon in the GS-1 period, although other lake-sediment proxies indicated a drier climate relative to the GI-1 period in the region. They also argued that because of the limited temperature decrease in the north Indian Ocean upwelling through the THC in the GS-1 period, the Indian summer monsoon was still dominant in the north Indian Ocean, Indian subcontinent, and even farther north, although the THC might have been weak during this period.
These patterns of oceanic and atmospheric circulations are consistent with our climate inferences using glacier–climate modeling. This suggests that the glacial events on the TP during the last glacial termination can be well correlated with the millennium-scale climate events in the North Atlantic region driven by the westerlies, and the Indian summer monsoon played a positive role in sustaining the glaciers under a warming climate trend in the region.
CONCLUSION
Using cosmogenic 10Be exposure dating and glacial modeling methods, we constrained the timing and climatic conditions for the late-glacial glacial events in the Pagele valley, Qiongmu Gangri peak, southern TP. Three lateral-frontal moraines (PM1, PM2, and PM3) were found to terminate at elevations of 5170, 5200, and 5230 m asl, respectively. These moraines indicate the glacier retreated about 1000 m from the PM1 to PM3, and the glacier's length was 3500 m at the PM1. After outliers are discarded, the mean (15,850 ± 980, 14,140 ± 880, and 12,430 ± 790 yr) and the oldest ages (16,190 ± 1000, 14,710 ± 910, and 12,980 ± 820 yr) from each of the three moraines all suggest that the PM1, PM2, and PM3 moraines respectively correspond well with the GS-2a (or Heinrich 1), GI-1 (or Bølling-Allerød), and GS-1 (or Younger Dryas) events in the Greenland ice-core record. The PM1 and PM3 moraines can be correlated to the late-glacial moraine found in the Barenduo valley and the Yanbajian inner moraine of the region, respectively (Owen et al., Reference Owen, Finkel, Barnard, Haizhou, Asahi, Caffee and Derbyshire2005; Chevalier et al., Reference Chevalier, Hilley, Tapponnier, Van Der Weord, Liu-Zeng, Finkel, Ryerson, Li and Liu2011). Using the glacier–climate model, we showed that the glacier had areas of 12.6, 11.6, and 10.9 km2 and volumes of 1.03–1.14, 0.79–0.88, and 0.61–0.68 km3 at the PM1, PM2, and PM3 moraine positions, respectively. Possible climatic scenarios of temperature and precipitation have also been established to support these three glacial extents. Considering more realistic precipitation values of 60%–80%, ~100%, and 80%–90% of present value, the model results indicated that the temperature in the region decreased by 2.6°C–2.9°C, ~1.6°C, and 1.4°C–1.5°C during the GS-2a, GI-1, and GS-1 periods, respectively. The temperature drop of 2.6°C–2.9°C in GS-2a is also compatible with the climatic reconstruction (2.6°C–2.8°C) for the late glacial in the Barenduo valley (Xu et al., Reference Xu, Dong, Pan, Hu, Bi, Liu and Yi2017b).
Taken together with information from oceanic and atmospheric circulations, these results imply that on the TP, the glacial events during the last glacial termination are well correlated with the millennium-scale climate events in the North Atlantic region by the westerlies, and the Indian summer monsoon plays a positive role in sustaining the glaciers under the warming climate trend. Therefore, any explanation for the relationship between glacier and climate during the last glacial termination on the TP must consider the interactions between the westerlies and the India summer monsoon in the region.
ACKNOWLEDGMENTS
This research was funded by the Second Tibetan Plateau Scientific Expedition and Research Program (STEP, 2019QZKK0101), National Natural Science Foundation of China (NSFC, grant no. 41771019), Strategic Priority Research Program (A) of the Chinese Academy of Sciences (XDA2007010201). We thank Jing Wei for processing the 10Be samples and Zhang Qian and Wu Yubin for their help in collecting boulder samples. Also, we would like to thank M.-L. Chevalier, the anonymous reviewer, J. Schulmeister (associate editor), and L. Owen (senior editor) for their suggestions and comments, which greatly improved the article.
SUPPLEMENTARY MATERIAL
The supplementary material for this article can be found at https://doi.org/10.1017/qua.2020.7.
