Reconstruction of paleoglacial equilibrium-line altitudes during the Last Glacial Maximum in the Diancang Massif, Northwest Yunnan Province, China

1 Introduction

The term “Tali Glacial Stage” [1,2,3], which was first proposed by Wissmann, is synonymous with the last glaciation in China. The Last Glacial Maximum (LGM) [4], which centered about 21 ± 3 ka B.P., is the latest extremely cold period in geological history and the most well-known period of the Ice Age [5,6]. The study of LGM is of great significance to understanding global paleoclimatic changes. As the named place of the “Tali Glacial Stage,” the Diancang Massif is not only the connecting point of the transition from the Tibetan Plateau to the Yunnan-Kweichow Plateau but also the breach area for the intrusion of the Indian Ocean monsoon northward and the transition belt of the Tibetan Plateau to the low relief of East and South Asia [3]. In addition, it is one of the southernmost mountains in the Asian continent that underwent glaciation during the Quaternary but without any modern glaciers, and also the first and frontal area where the southwest monsoon affects the Quaternary glaciation in the “Hengduan Mountains” region [7]. With these unique features and characteristics, the Diancang Massif has become an irreplaceable area for the study of inland paleoclimatic and paleo-environmental changes in China [3,8,9]. Accordingly, numerous glaciological studies have been conducted in this region. However, a detailed and systematic estimation of the paleoglacial equilibrium-line altitudes (ELAs) of the Diancang Massif has not been carried out thus far.

The ELA is a conceptual line, which refers to the altitude on a glacier where mass–balance equals zero, separating the zones of accumulation and ablation [10,11,12]. It is an important parameter in the study of former glaciers and can be used to estimate the ablation gradient, surface morphology, and paleo-environment (e.g., winter precipitation, summer mean temperature, annual average temperature), to name a few. ELAs are an important proxy of the fluctuations and mass–balance characteristics of former glaciers, and they have been widely used to reconstruct regional paleoclimatic conditions [13]. Various approaches have been conducted to determine glacier ELAs. Generally, we can broadly divide them into two categories. The first type is the direct field-observation method, which mainly uses modern surveying techniques such as global positioning system, remote sensing, and unmanned aerial vehicle to directly observe glacier ELAs in the field [14,15]. Its main advantages lie in its direct approach, effectiveness, and high accuracy. However, this type of method also has obvious shortcomings in that it is only suitable for modern glaciers, difficult to apply on a large scale, time consuming, and laborious. Furthermore, because glacial geomorphology is complex, field surveys and investigations are difficult to conduct in landforms [16], such as ridges, valleys, and steep slopes.

The second type is the indirect estimation method, in which glacier ELAs are estimated mainly on the basis of one or more representative indicators, such as glacial geomorphology, glacier area, glacier elevation, and glacier mass–balance. According to the information sources of the indicators used for estimation, this type of method can be further divided into two categories. The first is the estimation method based on geomorphological evidence, mainly including cirque-floor altitudes (CF), maximum altitude of lateral moraine (MALM), and toe-to-headwall altitude ratio (THAR) methods [17]. Among them, in the CF method, the attitude of the bottom of the cirque is considered to play an extremely important role in indicating glacier ELAs, approximately equivalent to ELA [18,19,20]. The MALM method considers the flow of the glacier to the center in the accumulation area and to the edge in the ablation area, the occurrence of net ablation only in the ablation area, and the deposition of marginal moraine only under ELA, resulting in the ELA being higher than the MALM [17,21]. The THAR method assumes that the glacier ELA can be approximated by some constant ratio between the altitude of the terminus and the head of the glacier [22]. These methods enable rapid and easy estimation of the position of glacier ELAs, but they are affected by geomorphological preservation or subjective factors, and their accuracy is low. Furthermore, the glacier area and altitude information are not taken into account while estimating glacier ELAs. The second is the parametric estimation approach, such as accumulation area ratio (AAR) and accumulation area balance ratio (AABR) method. This type of method is mainly based on assumed forms of the glacier mass–balance gradient and therefore can overcome the limitation associated with the nonconsideration of the glacier area and/or altitude information. For example, the AAR method [23] assumes that if glaciers are in a steady state, the ratio between the accumulation and ablation areas is constant, usually in the range of 0.5–0.8, with typical values around 0.55–0.65. However, it does not consider the mass–balance gradient [22,24] and hypsometry of glaciers [25]. Furthermore, it relies on similarity of glacier morphology and regimen over a region [26]. The AABR method [26,27,28] takes into account both the mass–balance gradient and hypsometry of glaciers and assumes that the mass–balance gradient consists of approximately two linear segments, normally with different slopes above and below the glacier ELAs. It uses a balance ratio for the estimation [12,25]. These methods are more accurate than the first type, especially in the case of sufficient data. In particular, the AABR method has the highest accuracy, and it is the most widely applied technique for glacier ELAs estimation.

As discussed earlier, no method is universally ideal for all situations, but many factors such as glacier type, catchment topography, climatic environment, and local conditions come into play, and the best approach is to use multiple methods wherever possible, with a range of values for relevant ratios [22]. In this study, several approaches, including CF, MALM, THAR, AAR, and AABR, were employed to estimate ELAs of former glaciers in the Diancang Massif during LGM on the basis of previous studies and the reconstruction of the surface of former glaciers. The results were compared and analyzed, and a weighted arithmetic average method was employed to determine the final results. The accuracy of the results was then evaluated by applying the law of propagation of measurement error [29,30]. Two critical parameters, differences between positive and negative effects and the distance ratio between the accumulation area and ablation area, are discussed in detail. The main objective of this study is to provide fundamental data for a better quantitative understanding of the paleoclimatic conditions and paleoglacial morphology of the Diancang Massif during LGM.

2 Materials and methods

2.1 Materials

2.1.1 Study area

The Diancang Massif (99°57′–100°12′E, 25°34′–26°00′N) (see Figure 1), located at the southeast end of the “Hengduan Mountains” and to the west of “Erhai Lake,” is a layered fault-block mountain with uplift in the fold-parallel ridge valley of western Yunnan Province, China. It is approximately 50 km long from north to south and 19–21 km wide from east to west. There are 19 peaks distributed in the north–south direction, all of which are more than 3,500 m above sea level. The highest is “Malong Peak,” with an elevation of 4,122 m. The lowland on the east slope and river valley on the west slope have elevations of approximately 1,980 and 1,700 m, respectively. The rock composition of the Diancang Massif mainly includes metamorphic rocks (granulite, marble, etc.), Triassic to Eocene sandstones, Permian basalts, and Devonian limestones. The geotectonic units belong to the Yanyuan-Lijiang platform margin depression belt on the western margin of the Yangtze block, and the geological evolution process is complex and has experienced multistage tectonic movement. Climatically, it is in a transitional belt between the south and middle subtropical zones, controlled by the Indian Ocean monsoon and westerlies [7] and is characterized by an annual alternation of dry and wet seasons.

Figure 1

Regional location of Diancang Massif. (a) SPOT5 satellite remote sensing image of the massive in September 2018 (2.5 m resolution) and (b) DEM (10 m resolution).

As the Diancang Massif has experienced the “Tali Glacial Stage,” it exhibits a typical glacial geomorphology, mainly featuring horns, aretes, cirques, cirque lakes, cirque thresholds, and small glaciated hollows. They are mainly distributed on both sides of the main ridge line between “Wutai Peak” in the north (3,581 m) and “Malong Peak” in the south (4,122 m). In this study, two representative glaciated sites, the north slope of “Yuju Peak” and the south slope of “Longquan Peak,” were selected as the study areas.

2.1.2 Data sources

We compiled existing data for former glaciers in the study area, including high-resolution satellite remote sensing images (SPOT5, 2.5 m, September 2018), high-resolution digital elevation models (DEMs; 10 m vertical distance), glacier dating data [3], concise glacier geomorphology maps [from ref. [8]], and other field investigation data such as the glacial geomorphic evidence position, glacial geologic maps, well-preserved terminal and lateral moraines, trimlines, and ice-marginal meltwater channels.

2.2 Methodology

2.2.1 Glacial geomorphic mapping

Accurate glacial geomorphological mapping [22,31,32,33,34,35] is essential for geometry reconstruction of former glaciers, and its results are usually taken as an important input in the reconstruction of paleoglacial ELAs. In traditional glacier geomorphological mapping, data are mainly collected through field investigation activities, such as fixed-point survey and recognition of geomorphic evidence (e.g., lateral-terminal moraines and trimlines), and indoor mapping is usually conducted using remote sensing images, DEMs, and other data [36]. However, the identification of historic landforms and sediments is difficult because diagnostic morphological features of glaciers are often destroyed by intense fluvial and glacial erosion [37], which poses great challenges to the accurate glacial geomorphological mapping and reconstruction of the morphology of former glaciers. To overcome these challenges, we propose the integration of existing research results, such as concise glacier geomorphology map and glacier dating data (data taken from ref. [3]), high-resolution satellite remote sensing images, high-resolution DEM, and field-investigation data, supported by relevant ArcGIS terrain analysis tools to complete the glacial geomorphological mapping of the study area. The basic process is shown in Figure 2.

Figure 2

Flowchart of glacier geomorphology mapping.

The basic idea of glacial geomorphological mapping is as follows: First, the ArcGIS spatial analysis toolbox is adopted to extract topographic features based on DEM data, such as ridgelines, valley lines, slope, aspect, slope variability, contours, and so forth. And then the concise glacier geomorphology map is rectified and re-projected according to the extracted topographic features. Among which, terrain features were obtained by applying a comprehensive method of identifying positive and negative terrain. A method combining plane curvature and slope shape was used to extract the glacier range of the LGM of the “Tali Glacial Stage.” In this method, large values of plane curvature on the positive terrain correspond to ridges, and large values of plane curvature on the negative terrain correspond to valleys. The width of ridges and valleys extracted by the combination of plane curvature and slope shape can be adjusted by selecting the size of plane curvature, and the corresponding study area can further indicate the upper edge line of the cirque wall. Considering the practicality and capability of slope variability (SOA) to represent plane curvature, plane curvature was replaced by SOA in the extraction process.

Second, taking the rectified and re-projected concise glacial geomorphology map and the superimposed high-resolution satellite remote sensing images (SPOT5) as the base map, combined with the positioning and dating data and field-investigation data, the glacial geomorphological features are interpreted and vectorized indoors, and the glacial period of each feature is determined.

2.2.2 Reconstruction of paleoglacial surface

The reconstruction of former glaciers is crucial for the estimation of ELAs and should be based on a rigorous interpretation of glacial geomorphology, as well as an assessment and testing of the reconstruction technique of the glacier surface and its physical plausibility [38]. At present, all reconstruction models of former glaciers satisfy the shallow ice approximation, which basically ignores abrupt changes in the vertical direction of any part of the glacier. Similarly, because changes in the paleoglacial surface cannot be observed, restored glaciers are assumed to be in a stable state and have a perfect-plasticity rheology, that is, the surface morphology of the former glacier remains unchanged. In this study, the python-based GlaRe toolkit developed by ref. [38] was employed for the reconstruction of former glaciers in the study area.

2.2.3 Reconstruction of paleoglacial ELAs

The range of methods for reconstructing paleoglacial ELAs has been extensively reviewed in previous research works, such as Meierding [18,27,39,40,41]. Common methods for determining the position of paleoglacial ELAs include CF, MALM, THAR, AAR, and AABR. As mentioned earlier, these are indirect estimation methods and each offers different advantages and disadvantages. The use of any of these methods may lead to inaccurate or even incorrect results. Therefore, this study adopted several approaches to comprehensively reconstruct ELAs. In this article, a python-based GIS tool developed by ref. [24] is employed for paleo-ELA reconstruction.

2.2.4 Determination of the final paleoglacial ELA’s value and accuracy evaluation

In this study, a weighted arithmetic average method was implemented to determine the final paleoglacial ELA’s value, and the weighted root mean square error was applied to evaluate the accuracy of the results, which is a good metric of weighted average model performance and can represent the error more precisely. The relevant calculation formulas are as follows:

(1) ELA final = ∑ i = 1 n ( p i × ELA i ) ∑ i = 1 n p i ,

(2) m final = ± ∑ ( p i v i v i ) n − 1 ,

(3) v i = ELA i − ELA final ,

where ELA final is the final weighted arithmetic average value of paleoglacial ELAs; ELA and p are the paleoglacial ELA value and weight of the corresponding estimation method, respectively; m final is the median error of the weighted arithmetic mean of the paleoglacial ELA.

For the determination of weight, the factor comparison method was adopted, and each estimation method was regarded as a factor affecting the final ELA result. These factors were then compared in pairs. If factor A was more important than factor B, A and B were assigned 1 and 0, respectively. If they were equally important, both were assigned 0.5. In this manner, a diagonal matrix could be constructed; the value of each element in each row is the weight of the factor (Table 1).

Table 1

Results of weight calculation using the factor comparison method

	CF	MALM	THAR	AAR	AABR	Total of the comparison	Weight
CF	0.5	0.5	0	0	0	1	0.08
MALM	0.5	0.5	0	0	0	1	0.08
THAR	1	1	0.5	0	0	2.5	0.20
AAR	1	1	1	0.5	0	3.5	0.28
AABR	1	1	1	1	0.5	4.5	0.36

3 Results and analysis

3.1 Reconstruction results

According to the method described in Section 2.2.2, paleoglacial flowlines and a raster DEM of the paleoglacial bed are required as mandatory input parameters. Numerous studies have generated paleoglacial flowlines [42,43,44,45]. This study employed a least-cost-path approach based on the graph-theory combined with the morphology of glacial valleys [43], which can better reflect the characteristics of geomorphology and achieve the best effect in the selection of glacier flowlines. The extracted results are shown in Figures 3 and 4.

Figure 3

Graphical estimation of ELA using the THAR method for grouped glaciers in the Rwenzori (from Kaser and Osmaston, 2002).

Figure 4

Distribution map of the flowlines of former glaciers. (a) Flowlines on the south slope of “Longquan Peak” and (b) flowlines on the north slope of “Yuju Peak.”

After creating the flowlines, the GlaRe toolkit was applied, and the ice thickness was initially calculated along the glacier flowlines, whose parameters for calculation are the interval distance (10 m) and the default shear stress (∼54 kPa) and then interpolated across the full extent of the glacier surface. The results are shown in Figure 5. Figure 6 shows a comparison of the subglacial geomorphology with the reconstructed former glacier surface, selecting one flowline from the north slope of “Yuju Peak” and other from the south slope of “Longquan Peak.”

Figure 5

Reconstructed surface of former glaciers. (a) Reconstructed glacier surface on the south slope of “Longquan Peak” and (b) reconstructed glacier surface on the north slope of “Yuju Peak.”

Figure 6

Comparison of profiles between the reconstructed surface of former glaciers and the subglacial geomorphology along the main flowline. (a) Profile of the south slope of “Longquan Peak” and (b) profile of the north slope of “Yuju Peak.”

According to the principle of glacier dynamics and from the overall morphology of the reconstruction, the ice surfaces of these two sites are similar, which agrees with the actual situation of the distribution of glacial relics. Comparing the subglacial topography and the ice thickness curve, elevation was found to increase suddenly and ice thickness was found to decrease abruptly when the subglacial topography encountered a steep slope. In reality, ice velocity at the ice ridge increases, resulting in the accumulation of ice less than that at the front and back of the glacier ridge.

3.2 Estimated results of paleoglacial ELAs

In this study, paleoglacial ELAs were estimated using five different methods, CF, MALM, THAR, AAR, and AABR. Among them, AAR and AABR were directly employed using a Python-based toolkit developed by ref. [24]. AABR values ranged between 0.9 and 3.0, and the interval of iteration for AABR was 0.1. For AAR, the suggested ratios ranged from 0.4 for maritime glaciers to 0.8 for polar glaciers, and the interval of iteration was 0.05. Considering that the resolution of DEM is 10 m, the altitude interval for surface volume calculation was set to 10 m. The value of THAR was 0.46. The iterative calculation results of the AAR and AABR methods are shown in Figure 7.

Figure 7

Distribution chart of paleoglacial ELA values calculated by AAR and AABR.

As shown in Figure 7, the paleoglacial ELAs of the south slope of “Longquan Peak” are larger than those of the north slope of “Yuju Peak”; moreover, the values calculated using AABR are more concentrated, whereas the values calculated using AAR show an uneven distribution. The estimation results of all methods adopted in this article are summarized in Table 2.

Table 2

Summary of estimation results of paleoglacial ELAs by various methods (in m a.s.l.)

	CF	MALM	THAR	AAR	AABR
South slope of “Longquan Peak”	3,830	3,850	3,837	3,931	3,891
North slope of “Yuju Peak”	3,780	3,785	3,611	3,873	3,783

Overall, the results from these five estimation approaches (Table 2) show that AAR provides higher values, whereas THAR provides lower values, especially for the north slope of “Yuju Peak.” According to the method described in Section 2.2.4, we finally determined the final paleoglacial ELA of the south slope of “Longquan Peak” to be 3,883 m a.s.l. and that of the north slope of “Yuju Peak” to be 3,773 m a.s.l., with median errors of 18.73 and 45.00 m, respectively.

4 Discussion

The Diancang Massif formed by an uplift that occurred during MIS 3b (54–44 ka B.P.), with numerous summits attaining an altitude of more than 3,500 m a.s.l. The geomorphological evolution of this area is well documented. However, a detailed and systematic estimation of paleoglacial ELAs during LGM is lacking. In this study, various approaches were employed to estimate paleoglacial ELAs based on the reconstructed glacial surfaces of the south slope of “Longquan Peak” and the north slope of “Yuju Peak.” The final paleoglacial ELAs of these two areas were then determined by the weighted arithmetic average method, and the median error of the final results were assessed.

“Yuju Peak” and “Longquan Peak” are located in the central and southern part of the Diancang Massif, respectively, which feature the most well-developed glacial geomorphology. The arithmetic mean of their paleoglacial ELAs is 3,828 m a.s.l., which can represent the paleoglacial equilibrium-line characteristics of typical glacier development in the Diancang Massif. This final result is in good agreement with the simple assumption (3,800  m a.s.l.) of Cui et al. (2004) and the inferential result (3,833 m a.s.l.) of ref. [46]. In order to better understand this final estimated result of the Diancang Massif, it was comprehensively compared with those of six other regions near the study area (Table 3).

As shown in Table 3, the paleoglacial ELA value of the Diancang Massif is relatively close to those of “Laojun Mountains” and “Xuebang Mountains,” indicating similar paleoclimatic conditions among them during LGM. At the same time, the paleoglacial ELA value of Diancang Massif is the second lowest among the seven regions, which agrees with the Diancang Massif being located at the southwest end of the “Hengduan Mountains.” Specifically, one of the main reasons is the difference in the periods of mountain uplift. For example, the “Qianhu Mountains,” “Laojun Mountains,” “Yulong Snow Mountains,” and other sites all retained complete glacier sequences from early to late “Tali Glacial Stage” [46], and the inverted secondary glaciation remained in the “Yulong Snow Mountains” [9,48], whereas the Diancang Massif was uplifted during MIS 3b. Therefore, if calculated on the basis of the uplift rate of 4.7 mm/annum since the late Pleistocene, the top of the Diancang Massif exactly corresponds to the lower boundary of glacier development over 3,700 m during that period. Accordingly, the later uplift and low-relief environments resulted in the lower altitudes of the equilibrium-line of the Diancang Massif.

The comparison between the morphologies of the main flowlines and the profiles of reconstructed glacier surfaces along the main flowlines in Figure 8 shows the occurrence of an obvious slope change in the paleoglacial ELA of the north slope of “Yuju Peak,” but not on the south slope of “Longquan Peak.” Elevation differences of the accumulation area on the north slope of “Yuju Peak” and the south slope of “Longquan Peak” are 198 and 225 m, respectively, and height differences of the ablation area are 225 and 275 m, respectively. Accordingly, differences in the positive and negative effects are 1:1.14 and 1:1.28, respectively. The distances of the accumulation area and ablation area on the north slope of “Yuju Peak” were 758 and 326 m, respectively, and the distance ratio between the accumulation area and the ablation area was 2.33:1. The distances of the accumulation area and ablation area on the south slope of “Longquan Peak” were approximately 548 and 660 m, respectively, and the distance ratio was 0.83:1. The aforementioned analysis shows that the distance ratio of the north slope of “Yuju Peak” agrees with the development characteristics of the cirque-floor in this area, but it was not obvious for the south slope of “Longquan Peak.” Moreover, the small change of the subglacial slope on the site also explains this situation. Combined with the subglacial ladder topography on the south slope of “Longquan Peak,” the ice at this site can be considered to have accumulated in many depressions.

Figure 8

Comparison between the morphology of the main flowlines and reconstructed glacier surface profiles of the study area. The left and right charts correspond to the south slope of “Longquan Peak” and the north slope of “Yuju Peak,” respectively. The x-axis is the length of the flowlines, and the y-axis is the height above sea level.

Given its position at the key interface facing the southwest monsoon, the reconstruction of paleoglacial ELAs in this area during LGM is of great significance and has scientific value for quantitatively understanding paleoclimatic conditions, which would guide more accurate modeling of paleoclimatic changes in this area.

5 Conclusions

In this article, various approaches were adopted to estimate the paleoglacial ELAs of the north slope of “Yuju Peak” and the south slope of “Longquan Peak” during LGM. The corresponding estimated results were compared and analyzed. The main findings of this study can be summarized as follows:

1. The values estimated using CF and THAR were relatively small, and the AABR method could provide a more reliable distribution of paleoglacial ELAs. From the comprehensive results, the final estimated paleoglacial ELAs of the north slope of “Yuju Peak” and the south slope of “Longquan Peak” were 3,773 and 3,883 m a.s.l., respectively. Compared with the mass–balance characteristics, the horizontal distance ratio between the accumulation area and the ablation area of “Yuju Peak” was 2.33:1, reflecting the glacier characteristics of cirques, and that of the south slope of “Longquan Peak” was 0.83:1. This could be attributable to the formation of a fault by the geological action of the east slope promoting the development of glacial geomorphology. The difference in the positive and negative effects of glacier action on the north slope of “Yuju Peak” could be quantified as 1:1.14 and that of the south slope of “Longquan Peak” as 1:1.28. The lower boundary of glaciation was almost the same, which could be attributed to the influence of geological structures on the slope.

2. The glacier surface on the north slope of “Yuju Peak” agrees with the theoretical glacial surface morphology and reflects the power function characteristics; the south slope of “Longquan Peak” has lower ice thickness and a flat glacier surface, reflecting a subglacial stepped geomorphology to a certain extent.

In summary, this study provides fundamental data for better quantitative understanding of the regional paleoclimatic conditions during LGM, thus contributing to more accurate modeling of paleoclimatic changes in this area.