Numerical Reconstruction of Three Holocene Glacial Events in Qiangyong Valley, Southern Tibetan Plateau and Their Implication for Holocene Climate Changes

: The dating of well-preserved Holocene moraines in the Qiangyong Valley, southern Tibetan Plateau (TP), o ﬀ ers great potential for reconstructing Holocene glacier extents and examining climate changes in the region. Guided by Holocene moraine features, this study used Geographic Information System (GIS) model tools to reconstruct paleo-glacier surfaces and glacier equilibrium line altitude (ELA) depressions for three Holocene glacial stages in the valley. The GIS-based models showed that the Qiangyong Valley contained ice volumes of 8.1 × 10 8 , 6.2 × 10 8 , and 4.6 × 10 8 m 3 during the early Holocene, Neoglacial, and Little Ice Age (LIA) glacial stages, and that the ELA was decreased by ~230 ± 25, ~210 ± 25, and ~165 ± 25 m, respectively, compared to modern conditions. Furthermore, the summer temperatures were estimated to be 1.56–1.79, 1.37–1.64, and 1.29–1.32 ◦ C cooler than present to support the three Holocene glacier extents, based on the evidence that the respective precipitation increased by 20–98, 13–109, and 0.9–11 mm relative to the present, which were derived from the lacustrine pollen data for the southern TP. By comparison, this study found that the amplitudes of the ELA-based summer temperature depressions were much larger than the pollen-based counterparts for the three glacial stages, although the two proxies both showed increasing trends in the reconstructed summer temperatures.


Introduction
Glacier changes over time provide clear evidence of climate variations [1]. Specifically, former glacier extents indicated by moraines and trimlines reflect temporal changes in glacier mass balance that were mainly controlled by past annual solid precipitation (snow) and summer temperature [2]. Therefore, paleoclimatic information can be inferred from glacial geomorphologic features by reconstructing the age of past mountain glacier extents [3,4]. The geomorphology-derived paleoclimate can provide an independent check on fidelity of other nearby terrestrial climate reconstructions and can help to assess the quality of climate model simulations. This is particularly true for the areas

Study Area
The Qiangyong Valley originates from the Mt. Kaluxung (28.8 • N, 90.3 • E), 6674 m above sea level (m asl), and drains northward to the Kaluxiong River on the northern slope of the Himalaya (Figure 1). The valley occupies an area of 28 km 2 , with elevation ranging from 4800 m asl at the confluence with the Kaluxiong River to 6674 m at the summit. In the valley, the contemporary Qiangyong glacier, with an area of~6.25 km 2 , flows down along a~4.9 km-long flow path and terminates at an elevation of 5000 m asl. Based on the field observation at the end of Qiangyong Glacier, the annual precipitation in 1975 was~489 mm [20,21]. Downstream of the modern glacier, three sets of latero-frontal moraines can be seen near the mouth of the valley ( Figure 2) [18,22]. The outermost moraine complex extends 4 km from the valley mouth and terminates at an altitude of 4810 m asl where some hummocks are enclosed by the latero-frontal moraine. The second moraine complex comprises distinct sharp-crested lateral ridges which lie 20-30 m above the present river on both sides of the valley. The moraine encloses the Daqiangyong Lake at 4900 m asl and stretches up for~3 km from the valley mouth to an elevation of 5000 m asl. The innermost latero-frontal moraine is relatively small, with its ridges lying 10-15 m above the valley floor. This moraine is about 1 km downstream of the modern glacier front and encloses the Xiaoqiangyong Lake at 4950 m asl. Using the cosmogenic 10 Be exposure dating method,  Owen et al. (2005) [18], while moraines in Gangbu Valley and Karola West Valley were dated by Liu et al. (2017) [19]. The Kaluxiong River and moraine crests were delineated from Google Earth image.  [19]. The Kaluxiong River and moraine crests were delineated from Google Earth image.

Derivation of Glacier Bed Topography
Reconstruction of the paleo-glacier surface requires knowledge of past glacier basal topography. Thus, evidence of post-glacial geomorphologic activities, which would change the glacier bed topography, should be considered when tailoring present-day topography to the real paleo-glacier bed conditions. The ASTER Global Digital Elevation Model V2 (GDEM V2), with 30 × 30 m spatial resolution (downloaded from Geospatial Data Cloud, http://www.gscloud.cn/, and then trimmed slightly beyond the watershed of the Qiangyong Valley), was corrected to represent the paleo-glacier basal topography. In the study area, the thickness of the contemporary glacier and the depth of the two lakes, enclosed by the Neoglacial and LIA moraines, had to be subtracted from the GDEM V2 data to obtain the basal topography of reconstructed paleo-glaciers.
We calculated the contemporary glacier bed topography starting from the current surface topography given by the GDEM V2, and subtracting from it the current glacier thickness, as calculated applying the VOLTA model developed by James and Carrivick (2016) [23]. We used the glacier boundary of the glacier reported in the Randolph Glacier Inventory (RGI) version 6.0 [24]. In the calculation, we used the VOLTA default method to assign the basal shear stress ( ) value on the flow line. The was calculated using the method proposed by Driedger and Kennard (1986) [25],

Derivation of Glacier Bed Topography
Reconstruction of the paleo-glacier surface requires knowledge of past glacier basal topography. Thus, evidence of post-glacial geomorphologic activities, which would change the glacier bed topography, should be considered when tailoring present-day topography to the real paleo-glacier bed conditions. The ASTER Global Digital Elevation Model V2 (GDEM V2), with 30 × 30 m spatial resolution (downloaded from Geospatial Data Cloud, http://www.gscloud.cn/, and then trimmed slightly beyond the watershed of the Qiangyong Valley), was corrected to represent the paleo-glacier basal topography. In the study area, the thickness of the contemporary glacier and the depth of the two lakes, enclosed by the Neoglacial and LIA moraines, had to be subtracted from the GDEM V2 data to obtain the basal topography of reconstructed paleo-glaciers.
We calculated the contemporary glacier bed topography starting from the current surface topography given by the GDEM V2, and subtracting from it the current glacier thickness, as calculated applying the VOLTA model developed by James and Carrivick (2016) [23]. We used the glacier boundary of the glacier reported in the Randolph Glacier Inventory (RGI) version 6.0 [24]. In the calculation, we used the VOLTA default method to assign the basal shear stress (τ b ) value on the flow line. The τ b was calculated using the method proposed by Driedger and Kennard (1986) [25], who suggested separating the glacier in elevation bands with areas (A i ) and mean slopes (α i ) as follows: where the A i is in m 2 and the τ b is in Pa. The basal shear stress obtained for the modern glacier is 162 kPa.
To measure the depth of the two lakes, we used a sonar detection system with Global Navigation Satellite System that was carried on an unmanned surface vehicle. The measuring precision is of ±1 cm, and the retrieved average depth of the two lakes was 3.90 m and 14.67 m for the Xiaoqiangyong and Daqiangyong lakes, respectively.

Reconstruction of Paleo-Glacier Surface
Paleo-glacier geometry reconstructions rely highly on morphologic evidence of former glacier extent (e.g., lateral and terminal moraines, glacier trimlines on bedrock). However, in reality, these glacial morphologic features are fragmentary and even totally missing, and often become degraded with time. In addition, these landforms often provide information exclusively along the margins of the paleo-glaciers, and the accuracy in reconstructing the surface topography sharply decreases with increasing distance from the margins, especially on the accumulation area. The best alternative is to use the available geomorphic evidence as constraints for numerical reconstructions of glacier thickness based on the mechanics of glacier flow [26,27]. The GlaRe, a GIS tool based on an iterative solution to the perfect plasticity assumption for glacier rheology [28], was utilized to numerically estimate the ice thickness distribution and paleo-glacier geometry for the three Holocene glacial events in the Qiangyong Valley. This tool first generates ice thicknesses along flow lines and then interpolate them towards the limits of the glacier. Specifically, the numerical approach iteratively solves the Equation (2) as a quadratic equation along the central glacier flowline.
where h is ice surface elevation, b is glacier bed elevation, H is ice thickness (in meters), ∆x is step length (in meters), τ av is average basal shear stress over an interval (in Pa), F is shape factor of cross-section, g is gravity acceleration (9.81 ms −2 ), and i is iteration step number. This is a derivation from the perfect plasticity assumption formula (τ = ρgHsinα) for the calculation of basal shear stress (τ) at the glacier bed [27], with ρ and α being ice density (~900 kgm −3 ) and ice surface slope, respectively. Details and explanations on the derivation of Equation (2) can be found in Benn and Hulton (2010) [29] and Van Der Veen (2013) [30]. For the reconstruction of paleo-glacier surfaces, we manually drew the flowlines according to the convex and concave directions on topographic contours of the valley. To calculate the F, we made cross-sections orthogonal to the flowlines at an interval of 50 m. For each paleo-glacier surface reconstruction, we applied the GlaRe tool using the catchment divide of the Qiangyong Valley, the lower point of the glacier as constrained by the frontal moraine, and the bedrock topography calculated by the VOLTA model as inputs. The glacier basal shear stress (τ av ), as a primary input to the GlaRe, makes the first-order control on the glacier surface elevation. Its values usually lie in the 50-150 kPa range for valley glaciers, indicated by field observational and experimental data [2,27]. Because the latero-frontal moraines of the three Holocene events can be clearly traced up to the contemporary glacier tongue, we tuned the τ av on cross sections (with a 50 m step length) along the flow line, until the modeled ice thicknesses matched the heights of the lateral and frontal moraines. For the middle and upper parts of the paleo-glacier during each glacial stage, we used the basal shear stress value obtained at the upper end of the lateral moraine [29]. The shear stress values that have been assigned along the flowlines of the three paleo-glaciers are shown in Figure 3. For topographically constrained glaciers, such as valley glaciers, not all the glacier driving stress is supported by the basal shear stress, because significant resistance to flow is also provided by lateral drag effect. To consider this effect in the basal shear stress formula (τ = ρgHsinα), the F factor (shape friction factor) was incorporated into the formula to obtain the adjusted basal shear stress (τa) where A is area of glacier cross-section (m 2 ) and p is perimeter of the cross-section (m) [2]. The GlaRe calculated the F factor on defined cross-sections and propagated values up the glacier until the next cross-section or the end of the flow line was reached. This tool automatically generated additional ice thickness points along the defined cross-sections, and thus helped to improve the interpolation accuracy, as errors in interpolation naturally increase with distance from an input point (on the flow line).

Reconstruction of Paleoglacier ELA and Paleoclimate
Reconstruction of ELA is a simple technique to reconstruct past climatic conditions from geomorphologic evidence of paleo-glacier extent, which is based on the fact that accumulation above For topographically constrained glaciers, such as valley glaciers, not all the glacier driving stress is supported by the basal shear stress, because significant resistance to flow is also provided by lateral drag effect. To consider this effect in the basal shear stress formula (τ = ρgHsinα), the F factor (shape friction factor) was incorporated into the formula to obtain the adjusted basal shear stress (τ a ) where A is area of glacier cross-section (m 2 ) and p is perimeter of the cross-section (m) [2]. The GlaRe calculated the F factor on defined cross-sections and propagated values up the glacier until the next cross-section or the end of the flow line was reached. This tool automatically generated additional ice thickness points along the defined cross-sections, and thus helped to improve the interpolation accuracy, as errors in interpolation naturally increase with distance from an input point (on the flow line).

Reconstruction of Paleoglacier ELA and Paleoclimate
Reconstruction of ELA is a simple technique to reconstruct past climatic conditions from geomorphologic evidence of paleo-glacier extent, which is based on the fact that accumulation above the ELA must be compensated by equal ablation below the ELA in a steady state glacier condition. Since the study of Porter (1975) [31], the method has often been used and became a traditional method to assess climate during Quaternary glaciations. Benn et al. (2005) [32] reviewed the methods of ELA reconstruction for former glaciers. In this study, we use the accumulation area ratio (AAR) and area-altitude balance ratio (AABR) methods on the reconstructed glacier surfaces to estimate the ELAs for the three Holocene glacier events in the Qiangyong Valley. For the AAR method, we assumed that the zero-balance AAR ranges between 0.5 and 0.8, as suggested by Benn et al. (2005) [32], Osmaston (2005) [33], and Dyurgerov et al. (2009) [34]. For the AABR method, we used a balance ration of 1.75 ± 0.71, which is the average value provided by Rea (2009) [35]. These methods have been implemented for the paleo-glacier in the Qiangyong Valley using the GIS-based tool developed by Pellitero et al. (2015) [36].
The differences in ELA (∆ELA) between these Holocene and modern glaciers represent the climate shifts associated with the changes in the Holocene glacier geometry. Specifically, Ohmura et al. (1992) [37], using linear perturbation analytic mathematics, derived a formula to relate the ∆ELA to changes in temperature (∆T) and precipitation (∆P) at the ELA where (∂P/∂T) ELA is the gradient of the empirical function f (P, T) = 0 that formulates the relationship between annual precipitation and summer temperature at the glacier ELA, ∂P/∂z is the gradient of precipitation with elevation, and ∂T/∂z is the temperature lapse rate in a specific area. Ohmura et al. (1992) [37] suggested the (∂P/∂T) −1 ELA ranges from 2.5 × 10 −3 to 3.3 × 10 −3 • C/mm on the basis of the best-fit polynomial regression curve P = 645 + 296T + 9T 2 for 70 globally-distributed glaciers. Using weather station records in the Kaluxiong River catchment, Zhang et al. (2015) [21] estimated the annual ∂T/∂z and ∂P/∂z to be −0.72 • C/100 m and 28 mm/year/100 m, respectively, using one year observation data on the glacier in 1975. We assumed that these parameters have not changed during the Holocene, and used the values reported by Zhang et al. (2015) [21] for calculating air temperature changes relative to the year of 1975 for the three Holocene paleo-glaciers, using Equation (5).

Glacier Extents of the Three Holocene Glacial Periods
The modern glacier thickness derived from the VOLTA model is shown in Figure 4. The VOLTA model estimated the modern glacier volume to be 3.0 × 10 8 m 3 , with maximum and average ice thickness being~155 and~47 m, respectively (Table 1).     To evaluate the validity of the AAR and AABR methods in calculating ELAs, we compared the modern glacier ELAs, calculated by using the GIS-based tool, with the ELAs estimated by using mass balance field observations. By locating at which elevation belt the observed annual mass balances was zero, we found that the zero-balance  [38,39], indicating that the two years did certainly not show specific balanced-budget conditions. We believe that the modern balanced-budget ELA lies somewhere between the zero-balance elevations of the two years (5800 and 6250 m asl). Using the AAR and AABR methods, we estimated the modern ELAs to be between 5700-6100 m ( Table 2), overlapping with the range of 5800-6250 m asl. This suggest that the AAR and AABR methods can be used to estimate the ELA changes in the Qiangyong Valley.   Table 1. The early Holocene glacier occupied an area of 10.0 km 2 and had a volume of 8.1 × 10 8 m 3 . Relative to the early Holocene glacier, the Neoglacial glacier decreased to an area of 8.45 km 2 , and the volume of this period decreased to 6.2 × 10 8 m 3 . The LIA glacier area continued to decrease to 7.68 km 2 , and the ice volume was 4.6 × 10 8 m 3 . The estimated average ice thickness values range from 47 to 83 m between the reconstructed glaciers, while the maximum values vary widely from 223 to 268 m.  The estimated ELAs using different AAR and AABR values are listed in Table 2. Obviously, the ELA obtained with the AAR method for AAR varying between 0.5-0.8 ranges between a large elevation interval, 500-600 m on average. Considering the AAR of 0.65, the ELAs of the three paleoglaciers are 5669, 5677, and 5735 m asl, respectively. Under the worldwide average AABR value of 1.75 [35], the ELAs for the three glacier periods are 5669, 5687, and 5735 m asl, respectively, which are comparable to the values from the AAR of 0.65. Compared to the modern glacier ELA of 5900 m asl (AAR = 0.65 or AABR = 1.75), the ELAs of the three paleo-glaciers were 230 ± 25, 210 ± 25, and 165 ± 25 m lower, respectively. The AABR method explicitly accounts for glacier surface hypsometry and uses glacier mass balance gradients with altitude below and above glacier ELA (BR value), thus yielding more reliable ELA reconstruction than the AAR method [33]. Considering an uncertainty (±0.71) of the global-wide AABR value [35], the ELAs differ by 150-200 m, less than the ELAs difference (500-600 m) derived by the AAR (0.5-0.8) for each of the three glacier periods (Table 2). Therefore, we chose the AABR-derived ELA values to discuss the climate implications.

Climate Implications of the Three Glacial Events
To establish the climate scenarios for the three glacial events using Equation (5), we set the ΔP to be from −100 to 100 mm relatively to the modern annual mean precipitation (with an incremental change of 25 mm), and found the corresponding ΔT values for the respective paleo-glacier ELAs ( Figure 6). Without precipitation change, the respective ELA changes required temperatures of 1.82-1.87, 1.66-1.71, and 1.30-1.34 °C lower than present for the three glacial events, with a range that depends on the ( ⁄ ) −1 range of 2.5 × 10 −3 -3.3 × 10 −3 °C /mm. If the annual precipitation were set to increase (decrease) by 100 mm, the summer temperature would be 1.54-1.57 (2.07-2.20), 1.38-1. 41 (1.91-2.04), and 1.01-1.05 (1.55-1.67) °C lower than present to sustain the three glacial events. Therefore, the ELA shift is highly controlled by the summer temperature, as can be expected by the different coefficients for ΔT and ΔP in Equation (5).
The uncertainties of calculated ΔT and ΔP for the three Holocene glacial stages are determined by the errors in ELA, temperature lapse rate and precipitation gradient (Equation (5)). According to Pellitero et al. (2015) [36], the uncertainty in the estimated ELA mainly derives from the delineation of the paleo-glacier boundaries. Pellitero et al. (2015) [36] also suggested that the calculation error (i.e., the minimum error) is equal to the chosen contour interval. In calculation of the ELA, we used the contour interval of 50 m, so the uncertainty of the reconstructed ELA is about ±25 m. Using Equation (5), we tested the sensitivity of the inferred paleoclimate to the uncertainties related to the ELA, temperature lapse rate and the precipitation gradient. The ±25 m error in the ELA results in an The estimated ELAs using different AAR and AABR values are listed in Table 2. Obviously, the ELA obtained with the AAR method for AAR varying between 0.5-0.8 ranges between a large elevation interval, 500-600 m on average. Considering the AAR of 0.65, the ELAs of the three paleo-glaciers are 5669, 5677, and 5735 m asl, respectively. Under the worldwide average AABR value of 1.75 [35], the ELAs for the three glacier periods are 5669, 5687, and 5735 m asl, respectively, which are comparable to the values from the AAR of 0.65. Compared to the modern glacier ELA of 5900 m asl (AAR = 0.65 or AABR = 1.75), the ELAs of the three paleo-glaciers were 230 ± 25, 210 ± 25, and 165 ± 25 m lower, respectively. The AABR method explicitly accounts for glacier surface hypsometry and uses glacier mass balance gradients with altitude below and above glacier ELA (BR value), thus yielding more reliable ELA reconstruction than the AAR method [33]. Considering an uncertainty (±0.71) of the global-wide AABR value [35], the ELAs differ by 150-200 m, less than the ELAs difference (500-600 m) derived by the AAR (0.5-0.8) for each of the three glacier periods (Table 2). Therefore, we chose the AABR-derived ELA values to discuss the climate implications.

Climate Implications of the Three Glacial Events
To establish the climate scenarios for the three glacial events using Equation (5), we set the ∆P to be from −100 to 100 mm relatively to the modern annual mean precipitation (with an incremental change of 25 mm), and found the corresponding ∆T values for the respective paleo-glacier ELAs ( Figure 6). Without precipitation change, the respective ELA changes required temperatures of 1.82-1.87, 1.66-1.71, and 1.30-1.34 • C lower than present for the three glacial events, with a range that depends on the (∂P/∂T) −1 ELA range of 2.5 × 10 −3 -3.3 × 10 −3 • C/mm. If the annual precipitation were set to increase (decrease) by 100 mm, the summer temperature would be 1.54-1.57 (2.07-2.20), 1.38-1. 41 (1.91-2.04), and 1.01-1.05 (1.55-1.67) • C lower than present to sustain the three glacial events. Therefore, the ELA shift is highly controlled by the summer temperature, as can be expected by the different coefficients for ∆T and ∆P in Equation (5). uncertainty of ±0.2 °C for the inferred temperature. Varying the temperature lapse rate of ±0.10 °C/100 m makes the inferred temperature change by ±0.23 °C; similarly, varying the precipitation gradient by ±10 mm/year/100 m, leads to variations in the inferred temperature of ±0.06 °C. We therefore suggest that the uncertainty in the inferred summer temperature change (ΔT) for each of the Holocene glacial stages is of about ±0.53 °C . Equation (5) gives no unique paleoclimate solution for each of the three glacial events. The estimates of ΔT − ΔP combinations shown in Figure 6 are just possible climatic scenarios. However, these solutions for ΔT − ΔP combinations are valuable in assessing other quantitative climatic-proxybased reconstructions. Using pollen assemblage data from lake and peat sediments, Chen et al. (2020) [40] quantitatively reconstructed the summer temperature changes and annual precipitation variations during the Holocene period at different parts of the TP. For the southern TP and nearest to our study area, the Lake Hidden and Qongjiamong Co were chosen by Chen et al. (2020) [40] to reconstruct the Holocene summer temperature and annual precipitation changes. The pollen-based climate reconstructions from Lake Hidden showed that the average annual precipitation was 98, 109, and 11 mm higher than the present during the early Holocene (10.8-8.9 ka), Neoglacial (3.7-1.9 ka) and LIA (0.4-0.1 ka), respectively; from Qongjiamong Co, the corresponding values were 20, 13, 0.9 mm higher than the present, respectively. If these values were correct, the model of Equation (5) would predict that the summer temperature was 1.56-1.79, 1.37-1.64, 1.29-1.32 °C lower than the present during the three Holocene glacial stages, respectively. This suggests an increasing temperature trend between the three Holocene glacial periods. The summer temperature in the Hidden Lake area was averagely 0.90, 0.38, 0.13 °C lower than present in the three Holocene glacial periods, respectively; and in the Qongjiamong Co area, it was averagely 0.11, 0.59 °C higher than present in the early Holocene and Neoglacial periods, but slightly lower (0.03 °C ) than present in the LIA. Therefore, the lake pollen-based summer temperature reconstructions, especially from Lake Hidden, also indicate an increasing temperature trend between the three Holocene glacial periods (Figure 7). Despite this increasing trend in both of the reconstructions, the amplitudes of the ELAbased summer temperature depressions relative to the present are much larger than the pollen-based counterparts (0.90, 0.38, 0.13 °C at Lake Hidden site) during the three glacial stages. Such difference between the ELA-based and pollen-based reconstructions may partly result from the accuracy of the The uncertainties of calculated ∆T and ∆P for the three Holocene glacial stages are determined by the errors in ELA, temperature lapse rate and precipitation gradient (Equation (5)). According to Pellitero et al. (2015) [36], the uncertainty in the estimated ELA mainly derives from the delineation of the paleo-glacier boundaries. Pellitero et al. (2015) [36] also suggested that the calculation error (i.e., the minimum error) is equal to the chosen contour interval. In calculation of the ELA, we used the contour interval of 50 m, so the uncertainty of the reconstructed ELA is about ±25 m. Using Equation (5), we tested the sensitivity of the inferred paleoclimate to the uncertainties related to the ELA, temperature lapse rate and the precipitation gradient. The ±25 m error in the ELA results in an uncertainty of ±0.2 • C for the inferred temperature. Varying the temperature lapse rate of ±0.10 • C/100 m makes the inferred temperature change by ±0.23 • C; similarly, varying the precipitation gradient by ±10 mm/year/100 m, leads to variations in the inferred temperature of ±0.06 • C. We therefore suggest that the uncertainty in the inferred summer temperature change (∆T) for each of the Holocene glacial stages is of about ±0.53 • C.
Equation (5) gives no unique paleoclimate solution for each of the three glacial events. The estimates of ∆T − ∆P combinations shown in Figure 6 are just possible climatic scenarios. However, these solutions for ∆T − ∆P combinations are valuable in assessing other quantitative climatic-proxy-based reconstructions. Using pollen assemblage data from lake and peat sediments, Chen et al. (2020) [40] quantitatively reconstructed the summer temperature changes and annual precipitation variations during the Holocene period at different parts of the TP. For the southern TP and nearest to our study area, the Lake Hidden and Qongjiamong Co were chosen by Chen et al. (2020) [40] to reconstruct the Holocene summer temperature and annual precipitation changes. The pollen-based climate reconstructions from Lake Hidden showed that the average annual precipitation was 98, 109, and 11 mm higher than the present during the early Holocene (10.8-8.9 ka), Neoglacial (3.7-1.9 ka) and LIA (0.4-0.1 ka), respectively; from Qongjiamong Co, the corresponding values were 20, 13, 0.9 mm higher than the present, respectively. If these values were correct, the model of Equation (5) would predict that the summer temperature was 1.56-1.79, 1.37-1.64, 1.29-1.32 • C lower than the present during the three Holocene glacial stages, respectively. This suggests an increasing temperature trend between the three Holocene glacial periods. The summer temperature in the Hidden Lake area was averagely 0.90, 0.38, 0.13 • C lower than present in the three Holocene glacial periods, respectively; and in the Qongjiamong Co area, it was averagely 0.11, 0.59 • C higher than present in the early Holocene and Neoglacial periods, but slightly lower (0.03 • C) than present in the LIA. Therefore, the lake pollen-based summer temperature reconstructions, especially from Lake Hidden, also indicate an increasing temperature trend between the three Holocene glacial periods (Figure 7). Despite this increasing trend in both of the reconstructions, the amplitudes of the ELA-based summer temperature depressions relative to the present are much larger than the pollen-based counterparts (0.90, 0.38, 0.13 • C at Lake Hidden site) during the three glacial stages. Such difference between the ELA-based and pollen-based reconstructions may partly result from the accuracy of the 10 Be exposure ages for the three glacial stages. For example, the 10 Be ages frame, provided by Owen et al. (2005) [18] and Liu et al. (2017) [19], can be estimated to be systematically older than published ones (Table S1), when recalculated with the updated 10 Be production scaling models in the online calculator of CRONUS-Earth V3 [41]. Specifically, the recalculated sample ages for the LIA, Neoglacial and early Holocene stages respectively become a few decades, about 1000 years and 3000 years older than the corresponding ages published by Owen et al. (2005) [18] and Liu et al. (2017) [19]. This is caused by different 10 Be production scaling models. Until now, there is no consensus about which scaling model is the best [42]. Moreover, when establishing the pollen-climate functions, the inherent "edge effects" within weighted averaging partial least squares (WA-PLS) method would overestimate the target climatic variable (temperature or precipitation) at its low-gradient end and underestimate it at the high-gradient end [40,43]. These factors may all contribute to the difference between the ELA-based and pollen-based reconstructions. Therefore, more investigations are encouraged to understand the possible reasons for the detected discrepancies between the glacier-model based and pollen-based temperature reconstructions.  [19], can be estimated to be systematically older than published ones (Table S1), when recalculated with the updated 10 Be production scaling models in the online calculator of CRONUS-Earth V3 [41]. Specifically, the recalculated sample ages for the LIA, Neoglacial and early Holocene stages respectively become a few decades, about 1000 years and 3000 years older than the corresponding ages published by Owen [19]. This is caused by different 10 Be production scaling models. Until now, there is no consensus about which scaling model is the best [42]. Moreover, when establishing the pollen-climate functions, the inherent "edge effects" within weighted averaging partial least squares (WA-PLS) method would overestimate the target climatic variable (temperature or precipitation) at its low-gradient end and underestimate it at the high-gradient end [40,43]. These factors may all contribute to the difference between the ELA-based and pollen-based reconstructions. Therefore, more investigations are encouraged to understand the possible reasons for the detected discrepancies between the glaciermodel based and pollen-based temperature reconstructions.

Conclusions
On the basis of 10 Be exposure age frame for Holocene moraines [18,19], we have reconstructed paleo-glacier surfaces for three Holocene glacial stages in the Qiangyong Valley, southern TP, using the GIS tool GlaRe [28]. On the reconstructed paleo-glacier surfaces, the ELAs were derived using the  [19]), Neo-Neoglacial, EH-Early Holocene, LIA-Little Ice Age; comparison of the ELA-based with the pollen-based summer temperature changes relative to the modern conditions for the three Holocene glacial periods, with the red dots and black triangulars being the respective average values of these changes in Lake Hidden and Qongjiamang Co areas (c); and the pollen-based summer temperature (d) and annual precipitation (e) during the Holocene period from Lake Hidden and Qongjiamong Co (Data from Chen et al. (2020) [40]).

Conclusions
On the basis of 10 Be exposure age frame for Holocene moraines [18,19], we have reconstructed paleo-glacier surfaces for three Holocene glacial stages in the Qiangyong Valley, southern TP, using the GIS tool GlaRe [28]. On the reconstructed paleo-glacier surfaces, the ELAs were derived using the AAR (0.5-0.8) and AABR (1.75 ± 0.71) methods for the three glacial stages. Constrained by the Holocene moraine positions, the model results suggested that the Qiangyong Valley contained ice volumes of 8.1 × 10 8 , 6.2 × 10 8 , and 4.6 × 10 8 m 3 with glacier areas of 10.0, 8.45, and 7.68 km 2 during the early Holocene, Neoglacial, and LIA glacial stages, respectively. The ELA depressions compared to current glacier geometry were 230 ± 25, 210 ± 25, and 165 ± 25 m during the three Holocene glacial stages, respectively. Without precipitation changes relative to the present, these ELA depressions would have required temperature of 1.82-1.87, 1.66-1.71, and 1.30-1.34 • C lower than at present for the three glacial stages. Lower temperature depressions (1.56-1.79, 1.37-1.64, and 1.29-1.32 • C) would have been required to support the three reconstructed glacial stages when using precipitation anomalies inferred from pollen data, which suggests higher precipitation in the past compared to current conditions. By testing the sensitivity of the inferred temperature change (∆T) to the ELA, temperature lapse rate and precipitation gradient, we suggested that the uncertainty in the ∆T for each of the Holocene glacial stages is of about ±0.53 • C. Moreover, we found that the amplitudes of the ELA-based summer temperature depressions were much larger than the pollen-based counterparts during the three glacial stages, although the summer temperature reconstructions from the two methods both showed an increasing trend between the three Holocene glacial stages. Therefore, we highlight that the glacier ELA-climate model is valuable for quantitatively assessing other independent proxy-based climate reconstructions.

Conflicts of Interest:
The authors declare no conflict of interest.