Last Glacial climate reconstruction by exploring glacier sensitivity to climate on the southeastern slope of the western Nyaiqentanglha Shan, Tibetan Plateau

Abstract

Improvements in understanding glacial extents and chronologies for the southeastern slope of the western Nyaiqentanglha Shan on the Tibetan Plateau are required to understand regional climate changes during the Last Glacial cycle. A two-dimensional numerical model of mass balance, based on snow–ice melting factors, and of ice flow for mountain glaciers is used to assess the glacier sensitivity to climatic change in a catchment of the region. The model can reproduce valley glaciers, wide-tongued glaciers and a coalescing glacier within step temperature lowering and precipitation increasing experiments. The model sensitivity experiments also indicate that the dependence of glacier growth on temperature and/or precipitation is nonlinear. The model results suggest that the valley glaciers respond more sensitively to an imposed climate change than wide-tongued and coalescing glaciers. Guided by field geological evidence of former glacier extent and other independent paleoclimate reconstructions, the model is also used to constrain the most realistic multi-year mean temperatures to be 2.9–4.6°C and 1.8–2.5°C lower than present in the glacial stages of the Last Glacial Maximum and middle marine oxygen isotope stage 3, respectively.

1. INTRODUCTION

The Tibetan Plateau and its surrounding mountains are among the most prominent topographic features on the Earth. Through time, the uplift of the Tibetan Plateau has caused environmental changes, brought about by the initiation of the South Asian monsoon and the splitting of the Northern Hemisphere westerlies into two branches (Prell and Kutzbach, 1992; Raymo and Ruddiman, 1992; Benn and Owen, 1998). Moreover, the variations of the Northern Hemisphere ice sheet are expected to exert influences on the climate change of the Tibetan Plateau, through storm tracks of cold air and westerly circulations (Porter and An, 1995; Benn and Owen, 1998; Zhou and others, 1999). Therefore, the Tibetan Plateau has a profound connection to global and regional climate (Molnar and England, 1990; Owen and Benn, 2005). Understanding the Quaternary glaciations of the Tibetan Plateau can therefore yield valuable information on the evolution of Earth's climate system.

Previous studies have extensively focused on the timing and style of the Quaternary glaciations on the Tibetan Plateau, based on glacial landform mapping, relative age controls from landform weathering features and morphostratigraphy, and numerical dating technologies, including radiocarbon ¹⁴C, optically stimulated luminescence and cosmogenic radionuclides (CRN). Of importance are studies that have used CRN methods to constrain the timing of glaciations. About 1855 individual exposure ages from 485 glacial deposits have been reported across the Tibetan Plateau (Heyman, 2014). Several glacial periods across the Tibetan Plateau have been identified based on the CRN methods: the Little Ice Age, Neoglacial, mid-Holocene, early Holocene, late glacial, global Last Glacial Maximum (LGM; marine oxygen isotope stage 2 (MIS 2)), middle last glacial (MIS 3) and early last glacial (MIS 4) (Owen and others, 2008; Owen and Dortch, 2014). Furthermore, using a robust statistical method, Dortch and others (2013) and Murari and others (2014) respectively recognized 19 and 27 regional glacial stages across the semi-arid western regions and the monsoon-influenced regions of the Himalayan-Tibetan orogeny. These studies allow us to compare the timings of glacial events at different sites and thus to construct a chronological framework for regional glaciations. They also make it possible to understand the climatic drivers for the glaciations of the Tibetan Plateau by correlation with other independent climate-proxy records. Indeed, current studies of Tibetan Plateau glaciations are focused on these two points, enabling a reasonably complete picture of late Quaternary glaciations to emerge. However, quantitative reconstructions of the Tibetan Plateau glaciations, such as the glacier mass balance, ice volume and possible climate conditions during each glacial episode, still require further study on the basis of the chronological framework of the Tibetan Plateau glaciations.

Numerical modeling of former glaciers can be an effective method of inferring paleoclimate, particularly when constrained by glacier-related geomorphologic features (Hostetler and Clark, 2000; Plummer and Phillips, 2003; Blard and others, 2007). Paleo-glacier simulations therefore add understanding of glacial climatic conditions, while also providing an independent method by which to allow paleo-glaciers to be reconstructed based on geological or geomorphological field evidence (Le Meur and Vincent, 2003; Plummer and Phillips, 2003). During recent years glacier–climate modeling work has been carried out on parts of the Tibetan Plateau (Xu and others, 2013; Xu, 2014; Wang and others, 2015). However, because the vast Tibetan Plateau is influenced by different climate systems, such modeling work needs to be carried out in different parts of the region, in order to fully understand the relationships between glacier fluctuations and climatic changes. In this study, we focus our glacier–climate modeling efforts on the southeastern slope of the western Nyaiqentanglha Shan (Fig. 1). Specifically, this study aims to test glacier sensitivity to climate change and infer possible climate conditions that could have supported the glacier extents during the Last Glacial period, through simulating the glacier extents at a catchment scale under different climate scenarios.

Fig. 1. Location of the model domain on the southeastern slope of the western Nyaiqentanglha Shan and southern Tibetan Plateau. The red line delineates the model domain.

2. STUDY AREA

The western Nyaiqentanglha Shan stretches ~250 km from SW to NE on the southern Tibetan Plateau and includes more than 30 peaks over 6000 m asl, with the highest at 7117 m asl (Nyaiqentanglha Peak, Yu and others, 2013). The South Asian monsoon, carrying moisture from the Indian Ocean, dominates in summer; and the westerly circulations derived from the Mediterranean Sea and/or the Atlantic Ocean dominate in winter (Yatagai and Yasunari, 1998; Shi, 2002). This produces a distinct seasonal oscillation of the climate with warm–wet summers (May–September) and cold–dry winters (October–April). Observations at Dangxiong Weather Station (30°29′N, 91°06′E, 4200 m asl), indicate that in the period 1981–2010 the mean annual temperature was 2.06°C and the mean annual precipitation was 478 mm. More than 90% of the precipitation is delivered between May and September (Fig. 2).

Fig. 2. The monthly averaged temperature and precipitation for the period of 1981–2010 at the Dangxiong weather station (30°29′N, 91°06′E, 4200 m asl).

During the Late Quaternary, the western Nyaiqentanglha Shan was intensively glaciated, and many moraines were formed and preserved in the glaciated valleys (Li and others, 1986; Owen and others, 2005; Chevalier and others, 2011; Dong and others, 2014). Two sets of moraines are present near valley mouths on the southeastern slope of the western Nyaiqentanglha Shan (Fig. 3a). The inner moraine set consists of latero-frontal moraines that are typically perched within the valleys and extend down to the valley mouths. The moraines of this set were deposited by individual valley glaciers. The outer moraine set is characterized by subdued ridges radiating from the mountain front. These ridges extend ~2 km southeast of the range front and terminate at altitudes of 4760–4850 m asl. We interpret these as moraines formed by a coalescing glacier that advanced out of the mountains into the range foreland. The two sets of moraines indicate two glacial expansions relative to the contemporary status, corresponding to different climate conditions. Accordingly, this paper focuses on simulating the respective glacier extents constrained by the two moraine sets and inferring the possible climate conditions that supported these glacier advances.

Fig. 3. Comparison of the glacial distribution for (a) observed and (b) modeled glaciers under modern climate conditions, based on an ice DDF _i of 11.8 mm °C⁻¹ d⁻¹. Note that the LGM and mid-MIS 3 glacier limits are delineated by field geological evidence.

Using cosmogenic ¹⁰Be exposure dating, Chevalier and others (2011) dated five and 22 boulder samples for the inner and outer moraine sets, respectively. The ages on the inner moraine set range from 10.4 ± 0.9 to 20.3 ± 1.8 ka, calculated using the time-dependent ¹⁰Be production rate scaling model of Lal (1991) and Stone (2000), and the ages on the outer set spread even wider (from 12.6 ± 1.1 to 49.8 ± 4.4 ka). The wide spread of the ¹⁰Be ages makes it difficult to assign the moraines to specific glaciations. On the northwestern slope of the western Nyaiqentanglha Shan, Dong and others (2014) also dated two moraine sets near the mouth of the Payuwang Valley, which have similar morphostratigraphic settings to the moraines of this study. Dong and others (2014) assigned the respective ages of 23.8 ± 4.0 and 39.9 ± 3.7 ka to the two moraine sets and argued that the glaciers expanded during LGM and MIS 3 in the western Nyaiqentanglha Shan. In addition, by separating Gaussians from cumulative probability frequency distributions to obtain best fits on the ¹⁰Be data, Murari and others (2014) recalculated the ages for the two moraine sets. They assigned an age of 12.5 ± 3.1 ka to the inner moraine set and three ages of 58.1 ± 7, 35.3 ± 7.4 and 25.4 ± 8.6 ka to the outer set. However, using an updated global reference ¹⁰Be production rate and a more robust statistical method, Heyman (2014, in Appendix file) showed that the inner moraine set has an age of 22.2 ± 9.0 ka, and the outer set has ages of 56.0 ± 27.7, 55.4 ± 16.9 and 26.2 ± 12.2 ka. According to these ¹⁰Be ages, we ascribe the inner moraine set to the LGM. Although it is difficult to rule out the possibility that the outer moraine formed during MIS 4, there are other ¹⁰Be dates indicating that the glaciers advanced during MIS 3 (60–30 ka) influenced by the Asian summer monsoon in the Himalaya, Anyemaqen, Nianbaoyeze and Qilian Mountains (e.g. Finkel and others, 2003; Owen and others, 2003; Wang and others, 2013; Murari and others, 2014; Owen and Dortch, 2014). In addition, Shi and others (2000) investigated 23 glaciated sites distributed in 12 areas in Asia, Europe, America and Australia and concluded that during MIS 3, glacial expansion happened in those regions. Based on the isotopic record preserved in the Guliya ice core, Shi and Yao (2002) identified three sub-stages in MIS 3 (MIS 3a, b and c) and argued that the cold-and-wet climate of MIS 3b (44–54 ka) favored glacier expansion. Accounting for this evidence, combined with the recalculated ages by Murari and others (2014) and Heyman (2014), we speculate that the outer moraine set most likely corresponds to a glacier advance that occurred during mid-MIS 3, but acknowledge here that such age estimate is a working hypothesis until more data are available. Such understanding of the glacial chronologies provides a framework for modeling glacier extents and inferring related climatic changes.

3. METHODS

The paper simulates the glacier extents during the LGM and mid-MIS 3 that are constrained by the glacial geomorphologic features and cosmogenic ¹⁰Be exposure ages on the southeastern slope of the western Nyaiqentanglha Shan. To do this, we use a coupled two-dimensional (2-D) mass-balance and shallow ice-flow model that was previously used to model mountain glacier extents on the northwestern and eastern parts of the Tibetan Plateau (Plummer and Phillips, 2003; Xu and others, 2014). The basis and strategy for the modeling are presented in these previous studies, and here we focus on the model input and model parameterization that are necessary to run the model in the study area.

3.1. Input data

The input data for the glacier–climate model include a DEM to represent the subglacial bed topography and monthly averaged climate records of temperature and precipitation. A DEM for the study area, with 30 × 30 m² spatial grids was downloaded from the Geospatial Data Cloud (http://www.gscloud.cn/). From the DEM we extract a watershed as a model domain with an area of 116.3 km² on the southeastern slope of the western Nyaiqentanglha Shan. The catchment consists of four valleys, which are termed Langkaku, Zuolai, Qiongmuda and Ranbuqu from west to east (Fig. 3a). Within the catchment, we catalog 16 modern glaciers occupying an area of 18.2 km², nine of which are larger than 0.5 km² with the largest glacier being 5.87 km² in the Langkaku valley (Table 1). In the catchment, elevations rise from 4500 to 6300 m asl and down-valley floor slopes are <10°. This makes the shallow ice approximation equations an appropriate model choice (Le Meur and Vincent, 2003; Leysinger Vieli and Gudmundsson, 2004). No ice thickness data exist within the catchment. To obtain the subglacial bed topography, we estimate ice thickness on the trimmed DEM using a simple procedure that assumes that glacier ice behaves as a perfectly plastic material such that ice thickness (H) can be determined by surface slope (α) and basal shear stress (τ = 100 kPa) from H = τ/(ρgsinα). The subglacial bed topography is described by subtracting the estimated ice thickness from the 30 × 30 m² DEM grids.

Table 1. Cataloged modern glaciers with an area larger than 0.5 km² in the model domain

* Denotes glacier terminal.

We compile 30 years (1981–2010) of daily climate records of temperature and precipitation measured at the Dangxiong Station (30°29′N, 91°06′E, 4200 m asl), the nearest long-term meteorological station to our model domain. The daily temperature and precipitation records are then converted into monthly means and used as ‘reference’ variables for the interpolation on the 30 × 30 m² DEM grids. The station is located in the same basin as our model domain (southeast of the mountains), so the station can be representative of climate at the same elevation in the basin. Xie and others (2010) calculated the monthly temperature lapse rates based on 1 year (August 2006–July 2007) of temperature record from nine automatic weather stations (AWSs) set up on the southeastern slope of the western Nyainqentanglha Shan. The nine AWSs are located in the range of 30°28′–30°32′N, 91°02′–90°03′E and distributed at elevations between 4300 and 5500 m asl. Using the monthly averaged temperatures at the Dangxiong Station and the linear monthly temperature lapse rates (Table 2), we calculate the monthly temperature grids in the model domain. When modeling the glacier mass balance on the western Nyainqentanglha Shan, both Caidong and Sorteberg (2010) and Zhao and others (2014) adopted a precipitation–elevation gradient of 5%/100 m from the observations of Li and others (1986) in the region. Here we use the precipitation–elevation gradient of 5%/100 m and the monthly precipitations at Dangxiong station to calculate the monthly precipitation grids for the model domain. The glacier accumulation (monthly proportion of precipitation that falls as snow) is determined by the fraction of the month falling below the snowfall threshold (T _s) of 1°C (Anderson and others, 2006). Glacier ablation is calculated by the degree-day approach that considers the fraction of the month with temperature above 0°C. The climatic input (temperature, precipitation) used to define the modern glacier distribution can then be systematically changed (ΔT: temperature change relative to present; F _p: precipitation as fractional value relative to modern) to simulate the glacier extents under hypothesized paleoclimatic conditions for the specific glacial stages.

Table 2. Monthly temperature lapse rates used in the model

3.2. Model parameterization

The mass-balance model calculates the net annual mass balance of each DEM cell and determines the accumulation and ablation areas for the modeled region. The transfer of mass between the accumulation and ablation areas is calculated based on the ice-flow model used in previous studies (e.g. Plummer and Phillips, 2003; Kessler and others, 2006). The coupled mass-balance and ice-flow model has parameters including snow and ice degree-day factors (DDF _s, DDF _i), snowfall threshold (T _s), ice deformation and sliding coefficients (A, B), and ice deformation versus sliding factor (f). Based on the observational mass-balance data in the melt seasons of the year 2007–2008 on the Zhadang Glacier, western Nyainqentanglha Shan, Wu and others (2010) estimated the DDF _s to be 5.3 mm °C⁻¹ d⁻¹, and the DDF _i ranges from 4.0 to 14 mm °C⁻¹ d⁻¹. Because the DDF _i is highly uncertain relative to other parameters and has a direct impact on the mass-balance calculation, we choose to calibrate this parameter following the method developed by Xu (2014). Forced by the climate of the past 30 years (1981–2010), the model runs under different DDF _i (4.0–14 mm °C⁻¹ d⁻¹) to match the observed present glacier distribution. When the model best reproduces the observed glacier distribution in the model domain, we adopt the corresponding DDF _i value in the following modeling experiments. Moreover, we assume that the glaciers deposited these two sets of moraines when they were in or near equilibrium with the climate during the LGM and MIS 3 glacial stages, and thus a steady-state simulation can be made. For such steady-state glacier simulation, the net mass-balance gradient effectively constrains the glacier extent, and variations in the two ice transport processes of ice deformation and sliding lead to only a small modification in ice extent (Kessler and others, 2006). Therefore, we do not attempt to adjust the coefficients of ice deformation (A) and sliding (B) and prescribe values of 1.0 × 10⁻⁷ a⁻¹ kPa⁻³, 1.5 × 10⁻³ m a⁻¹ kPa⁻² and 0.5 for the parameters A, B and f, respectively, which are widely used in previous studies (Plummer and Phillips, 2003; Laabs and others, 2006; Refsnider and others, 2008).

3.3. Model experiments

Four sets of sensitivity tests are conducted by changing the climate variables from modern conditions and by evaluating location, extent and growth pattern of the resultant glaciers. The first experiment is to establish the degree of influence exerted by temperature change on glacier extent. In the experiment, incremental changes in temperature relative to present day (ΔT) from −0.5 to −5.5°C (with a step of 0.5°C) were applied uniformly across the model domain. As the simulated glacier extents approached the respective positions of the two moraine sets, the ΔT was changed by a 0.1°C step in subsequent simulations to obtain the best-fitting glacier extents with the geological field evidence. The second experiment is designed to test the glacier sensitivity to precipitation change. Using incremental changes in precipitation from 1.0 to 10.0 times present-day amount and no change in temperature, we will see how the glacier evolves in the model domain. Guided by the first two sensitivity tests, the third experiment sets a suite of ΔT − F _p combinations to force the model. By fitting the simulated glaciers under these ΔT − F _p combinations to the glacier extents constrained by the inner and outer moraine sets, we will identify the climate scenarios that support the LGM and mid-MIS 3 glacier extents. The model parameters that dominantly influence our paleoclimate estimates are DDF _i and T _s in the coupled mass-balance and ice-flow model. To evaluate the impacts of the two parameters on the paleoclimate inferences, the fourth experiment uses different parameter values of DDF _i (10.8 and 12.8 mm °C⁻¹ d⁻¹) and T _s (0 and 2°C) to estimate the LGM ΔT − F _p combinations by reproducing the LGM glacier extent in the model domain.

4. RESULTS

By running multiple simulations under current climate conditions with different DDF _i values, we found that when the DDF _i is 11.8 mm °C⁻¹ d⁻¹, the model can best reproduce the observed glacier distribution in the model domain (Fig. 3b). The modeled modern glaciers have a total area of 19.7 km², comparable to the total area (18.2 km²) of the cataloged 16 modern glaciers. The total glacier volume was estimated to be 0.83 km³, with the maximum ice thickness of 167 m in the Langkaku valley where the largest modern glacier is located. For the four sensitivity experiments, we present key results as following.

4.1. Step cooling with present precipitation

Figure 4 shows steady-state ice conditions under the step cooling scenario in the first experiment. With temperature reductions between 0 and −2.0°C, the glaciers are constrained in the four individual valleys; between −2.0 and −3.0°C lowering, the glaciers are characterized by wide-tongued shape; and beyond −3.0°C the glaciers coalesce on the foreland of the mountain range. Figure 5 shows the changes in ice area and volume under incremental cooling from −0.5 to −5.0°C. For one specific cooling, the model demonstrates that the valley glaciers change somewhat more intensively than the wide-tongued and coalescing glaciers. When the temperature decrease is between −0.5 and −2.5°C, the glacier volume will rise from 2.45 to 6.41 km³ with the maximum ice thickness and area respectively increasing from 269 to 443 m and 32.2 to 55.5 km². With imposed temperature reductions between −3.0 and −5.0°C, the glacier volume rises from 7.28 to 9.14 km³ with the maximum ice thickness and area increasing only from 452 to 468 m and 60.8 to 78.5 km², respectively. These results demonstrate the sensitivity of the model to temperature changes without precipitation change relative to the present day. We find that the model can reproduce the glacier extents constrained by the inner and outer moraine sets under temperature decreases of −2.3 and −3.2°C, respectively.

Fig. 4. Glacier growth under imposed step coolings of 0.5–5.0°C and present-day precipitation. Note that the glaciers vary from valley glaciers, through wide-tongued glaciers, to a coalescing glacier.

Fig. 5. Modeled glacier area, volume and their respective change rates as a function of the imposed step coolings of −0.5 to −5.0°C (a, b) and step precipitation increases of 1–10 times present values (c, d).

4.2. Step precipitation increasing with present temperature condition

When an incremental, spatially uniform precipitation increase is imposed across the region with no change in temperature, the general pattern of glacier growth is essentially the same as that in the situation of progressive temperature cooling (not shown). With precipitation increases up to 3.0 times present-day precipitation, the glaciers are distributed in the four individual valleys; and beyond 4.5 times, the glaciers expand and coalesce on the foreland of the mountain range. When the precipitation is from 2.0 to 4.0 times present, the glacier volume rises from 3.66 to 7.80 km³ with the maximum ice thickness and area respectively increasing from 343 to 495 m and 38.9 to 53.4 km². With precipitation from 5.0 to 7.0 times present, the glacier volume rises from 9.80 to 11.6 km³ with the maximum ice thickness and area increasing only from 527 to 569 m and 59.7 to 66.6 km², respectively (Fig. 5). The model can reproduce the glacier extents constrained by the inner and outer moraine sets under the precipitation scenarios 3.5 and 5.5 times present, respectively. We also note that for one specific modeled glacier area, the modeled volume by precipitation increasing alone is larger than that by temperature cooling alone.

4.3. Climate scenarios for the LGM and mid-MIS 3 glacier extents

Guided by the first two sensitivity tests, incremental changes in ΔT from −0.5 to −5.0°C and in F _p from 0.2 to 2.1 are applied uniformly across the model domain. By fitting the simulated glaciers under these ΔT − F _p combinations to the glacier extents constrained by the inner and outer moraine sets, we identified the climate scenarios that support the LGM and mid-MIS 3 glacier extents (Fig. 6a). In this experiment, for one specific glacier extent the modeled glaciers have larger volumes in the wetter conditions and smaller volumes in the colder conditions. This is typically the case for the model when ice-flux parameters are kept constant (Plummer and Phillips, 2003).

Fig. 6. Plots of the temperature and precipitation combinations (ΔT − F _p) (a) that yield LGM and mid-MIS 3 glacial extents, respectively and (b) glacier-parameter sensitivity test for the LGM climate conditions.

All ΔT − F _p combinations that can reproduce the LGM glacier extent are displayed on the upper curve in Figure 6a. Under these combinations, the modeled glaciers have an area of ~50.9 km²; the modeled glacier volumes however vary from 5.58 to 5.72 km³ with the maximum ice thicknesses ranging from 411 to 442 m. For the mid-MIS 3 glacier extent, all the successful ΔT − F _p combinations are shown on the lower curve in Figure 6a. Under these combinations, the modeled glacier areas are ~63.2 km²; the modeled glacier volumes however vary from 7.41 to 7.64 km³ with the maximum ice thickness ranging from 457 to 491 m. Figure 7 shows the modeled net annual mass balances and corresponding ice extents for the LGM and mid-MIS 3 glacial stages that match well with the field geological data under the two ΔT − F _p combinations of (−3.6°C, 0.5) and (−2.4°C, 1.5), respectively.

Fig. 7. Simulated glacier mass balances and extents for the LGM and mid-MIS 3 glacial stages under the ΔT − F _p combinations of (−3.6°C, 0.5) (a, b) and (−2.4°C, 1.5) (c, d), respectively. Note that the LGM and mid-MIS 3 glacier limits are also shown for comparison.

4.4. Impact of DDF _i and T _s on the climate inferences

To evaluate the impacts of the two parameters on the paleoclimate inferences, we performed the modeling experiments with different parameter values of DDF _i (10.8 and 12.8 mm °C⁻¹ d⁻¹) and T _s (0 and 2°C) in attempts to reproduce the LGM glacier extent in the model domain. Figure 6b shows the results of the model sensitivity test by varying the two parameters. The difference in modeled temperature lowering between using a DDF _i value of 10.8 and 12.8 mm °C⁻¹ d⁻¹ is −2.0 to −2.9°C, and the difference between a T _s value of 0 and 2°C is 1.1–1.4°C. If F _p is set to 1 (under modern precipitation), varying the DDF _i by ±1.0 mm °C⁻¹ d⁻¹ alone results in a ±1.3°C difference in ΔT, and a ±1°C difference in T _s modifies ΔT by ±0.7°C. Likewise, if we set the ΔT to be −3°C (3°C temperature depression relative to modern), varying the DDF _i by ±1.0 mm °C⁻¹ d⁻¹ alone results in a ±0.5 difference in F _p, and a ±1°C difference in T _s modifies F _p by ±0.25. Therefore, a difference of ±0.5 mm °C⁻¹ d⁻¹ in DDF _i and variation of ±1°C in T _s have approximately equal impact on our temperature and precipitation estimates. Similar results can be expected for the mid-MIS 3 glacial climate inferences, although we did not test for this glacial period.

5. DISCUSSION

Because glacier extent depends strongly on both temperature and precipitation, the results in this study do not provide a single paleoclimate solution for each of the two modeled glacial stages. However, the results do have implications for understanding patterns of glacier growth in response to climate change in the region, and the approximate magnitudes of climate change necessary to produce the LGM and mid-MIS 3 glacier extents. The differences in the inferred magnitude of climate change between LGM and mid-MIS 3 have further implications for understanding the likely driving mechanism for the glacier advances during the Last Glacial period.

5.1. Glacier sensitivity to climate change

In the step cooling and precipitation increasing experiments, the model reproduced the valley glaciers, wide-tongued glaciers and a larger coalescing glacier. The results of these experiments suggest that the valley glaciers respond more sensitively to a specific climate change than the larger wide-tongued and coalescing glaciers. For example, under 1°C cooling alone, the total volume of the modeled valley glaciers increases by 2.18 km³ on average, while the volume of the modeled coalescing glacier only increases by 1.17 km³. Moreover, the experiment results highlight that the dependence of glacier growth on temperature and/or precipitation is nonlinear, although the glacier area and volume always increase under increasingly cold or wet climates (Fig. 5). By calculating change rates of the glacier area and volume curves, we are able to identify how the change rates of glacier area and volume evolve with respect to these equal and incremental cold or wet climates. Figures 5a, b show that the change rates of glacier area (dA/dΔT) and volume (dV/dΔT) decrease steadily with cooling between −0.5 and −2.5°C as glaciers are constrained in the four individual valleys type, and once glaciers become wide-tongued glaciers with cooling between −2.5 and −3.0°C, the change rates slightly increase, probably reflecting the topographic difference between high-altitude valleys and low-altitude range forelands. Cooling between −3 and −3.5°C, however, leads to a decrease in the change rates as glaciers coalesce on the foreland, with temperatures colder than this leading to another increasing trend of change rates. Similar patterns can be seen when the increasing precipitation is imposed to force glacier growth (Figs 5c, d). The coupled mass-balance and ice-flow model directly uses climate (precipitation and temperature) data to determine net annual mass balance, and it runs on the DEM that describes topography of the model domain. Therefore, the varying characteristics of the modeled glacier area and volume change rates imply that some key climate–topography thresholds appear to exist in our model domain, which may control the formation process of different glaciers (Payne and Sugden, 1990). Although we did not simulate the recession of glaciers in the region, we predict the reverse to be true for recession scenarios (such as during deglaciation).

The model is a 2-D coupled approach where calculated grids of mass balance on a GIS are used as input to a glacier flow model. The algorithm solves for gridded mass balance first, and then calculates ice flow to an equilibrium state, indicating that in the model, mass balance is not dynamically coupled to flow regime. This prevents the exploration of how a specific paleoclimate scenario would have been reflected in glacier dynamics, erosion, or geomorphological formations. Also the model does not examine the timescale at which one specific glaciation responds to climate change. This therefore limits the interpretation for proposed climate shift, for example, from warmer/wetter conditions of more extensive mid-MIS 3 glacier extent to colder/drier LGM conditions. However, the model indeed provides a means to hypothesize reasonable combinations of temperature and precipitation for the two sets of moraines in the region, although the ¹⁰Be ages are not robust to constrain the mid-MIS 3 moraine set.

5.2. Climate conditions during the Last Glacial

For glacial stages of LGM and mid-MIS 3, our model simulations defined two specific sets of climatic ΔT − F _p combinations that can reproduce respective glacial extents (Fig. 6a). The curves denoted by the successful ΔT − F _p combinations are essentially parallel between the two glacial stages. As Xu (2014) suggested, this characteristic is a result of the assumption that the seasonal patterns of temperature lapse rate and precipitation for each glacial stage modeling are the same as present. The estimates of ΔT − F _p combinations shown in Figure 6a are the only possible climatic scenarios. To better constrain the ranges of the ΔT and F _p, we have to recognize the general climate features for each of the two glacial stages.

The western Nyaiqentanglha Shan receives moisture, which is dominantly transported by the South Asian summer monsoon from the Indian Ocean (Yatagai and Yasunari, 1998; Shi, 2002). As suggested by numerical simulations of atmospheric circulation patterns, the South Asian monsoon correlates well with the distribution of solar radiation as determined by the cyclical variability in the precession of the Earth's orbit (Clemens and others, 1991; Prell and Kutzbach, 1992). The precessional cycle of 30°N reveals strong summer insolation during MIS 3, which sometimes exceeds the present-day conditions and nearly equals early Holocene values, and subsequent to that the insolation steadily decreased and reached its minimum at the LGM (18–24 ka) (Berger and Loutre, 1991; Laskar and others, 2004). Furthermore, during the middle part of MIS 3 (MIS 3b, 44–55 ka), the summer insolation is less than that during the early and late MIS 3, but is characterized by a northward increasing trend with latitude (Berger, 1978; Shi, 2002). This would be favorable for the establishment of low pressure over the Tibetan Plateau and high pressure over the southern subtropical Indian Ocean, and thus enhanced pressure gradient, strengthened monsoon and increased precipitation (Zonneveld and others, 1997; Shi and Yao, 2002). By considering the Northern Hemisphere summer radiation and glacial–interglacial ice-sheet boundary conditions, Prell and Kutzbach (1987) modeled the precipitation in South Asia and argued that it was drier than today during MIS 4 and LGM, and wetter in the interglacial period of MIS 3. Compiling 75 paleoclimatic records (mainly lake and loess records), Herzschuh (2006) reconstructed the moisture evolution history in monsoonal Central Asia during the last 50 ka. Most records from the area also yield a wet climate during the middle and late MIS 3, after which there were dry conditions (until the LGM). Furthermore, from multiple lines of evidence including pollen, periglacial phenomena and ice cores from the Tibetan Plateau, Shi and others (1997, 2001) proposed that during MIS 3 precipitation was 40–100% higher than today, while during the LGM precipitation was 30–70% lower than today. By referencing these precipitation-related reconstructions, this study tentatively constrains F _p values to be 0.3–0.7 and 1.4–2.0 for the LGM and mid-MIS 3, respectively. Considering these precipitation amounts, our model results suggest that a temperature lowering between −2.9 and −4.6°C is required to reproduce the LGM glacial extent, and a temperature lowering between −1.8 and −2.5°C can support the mid-MIS 3 glacial extent in the model domain (Fig. 6a). This suggests that temperature depression is necessary to support the LGM and mid-MIS 3 glacier extents and that any reasonable scenario of precipitation change alone would not suffice. As shown in our experiments of step precipitation increases in the absence of temperature change, the LGM and mid-MIS 3 glacier extents require precipitation amounts of 3.5 and 5.5 times present, respectively. Such increases in precipitation are not consistent with the climate proxy data, indicating that temperature lowering is a necessary trigger for regional glaciation in the region. During the mid-MIS 3, the temperature lowering could be associated with less insolation relative to that during the early and late MIS 3 (Shi, 2002; Shi and Yao, 2002). Below we first consider our results in the context of other climate reconstructions in the southern Tibetan Plateau, and subsequently within the wider framework of the whole Tibetan Plateau.

The magnitudes of temperature changes necessary to produce the LGM and mid-MIS 3 glacier extents determined in this study are generally comparable with those found in other paleoclimate studies on the southern Tibetan Plateau. For example, in the Payuwang valley on northwestern slope of the western Nyaiqentanglha Shan, Xu and Glasser (2015) studied the glacier sensitivity to the equilibrium-line altitude (ELA) change. Based on the ELA depression, they inferred that during the LGM, the temperature was 3.3–4.4°C lower than the period of 2005–2008, and that the temperature during MIS 3 was slightly lower (<1.5°C lower) than the period of 2005–2008. Using a similar modeling method to this study, Xu (2014) estimated the LGM temperature to be 4.0–5.9°C lower than present in the Yingpu Valley, Queer Shan and southeastern Tibetan Plateau. In addition, based on biogeographical data from ground beetles of the genus Trechus, and from juniper tree haplotypes, Schmidt and others (2011) suggested a maximum LGM temperature depression of 3–4°C in July for the southern Tibetan Plateau. Miehe and others (2011) also argued that Alpine steppes of Tibet persisted during the LGM with 3–4°C lower summer temperature, by hypothesizing that the Alpine steppes of the Tibetan highlands remained ecologically stable during the LGM.

In the wider context of the whole Tibetan Plateau, it appears that during the glacial periods of the LGM and mid-MIS 3, the magnitudes of temperature depression are larger on the northern Tibetan Plateau than on the southern part. For example, from the δ ¹⁸O values of the Guliya ice core in the western Kunlun Shan, northwestern Tibetan Plateau, it was suggested that during the mid-MIS 3, the temperature was ~5°C lower than today; and during the LGM, the mean annual temperature was 6–9°C lower than today (Shi and others, 1997, 2001). Using a climate–glacier modeling method in the Kuzigun Valley, northwestern Tibetan Plateau, Xu and others (2014) showed that the LGM climate was 4–9°C colder than today. Moreover, based on the ELA depression in the Dalijia Shan, northeastern Tibetan Plateau, Wang and others (2015) argued that during the mid-MIS 3 temperature depression was 3.3–6.8°C. These estimates of temperature lowering are larger than in the southern Tibetan Plateau. This may reflect different degrees to which the Northern Hemisphere ice sheets influenced the northern and southern parts of the Tibetan Plateau during the glacial periods, as the climate of the northern Tibetan Plateau tended to be tuned to the climatic rhythm of the high northern latitudes by storm tracks of cold air and westerly circulations (Porter and An, 1995; Benn and Owen, 1998; Zhou and others, 1999; Seong and others, 2009).

6. CONCLUSIONS

In this paper, we tested glacier sensitivity to climate changes and inferred the climate conditions that could have supported the LGM and mid-MIS 3 glacier extents on the southeastern slope of the western Nyaiqentanglha Shan, using a coupled mass-balance and ice-flow model and guiding constraints from glacial geological data. Our sensitivity experiments explored glacier growth in response to temperature depression and precipitation increases. In the sensitivity experiments, the model reproduced the valley glaciers, wide-tongued glaciers and a coalescing glacier. The experiment results showed that the dependence of glacier growth on temperature and/or precipitation is nonlinear. It can also be concluded that on the southeastern slope of the western Nyaiqentanglha Shan, temperature depression is necessary to support the LGM and MIS 3 glacier extents and that any reasonable scenario of precipitation change alone would not suffice. Guided by field geological evidence and other independent paleoclimate reconstructions, the study constrained the more realistic temperatures to be 2.9–4.6°C and 1.8–2.5°C lower than present to support the glacier extents of the LGM and mid-MIS 3 glacial stages, respectively. Although the model did not dynamically couple the mass-balance to ice-flow regime and examined the timescale at which one specific glaciation responds to climate change, the temperature estimates for the two glacial periods are generally comparable with those found in other paleoclimate studies on the southern Tibetan Plateau. However, the last glacial temperature depressions estimated by the glacier–climate model are smaller than those reported from on the northern Tibetan Plateau.