Serial Landsat images were used, including Landsat 5 Thematic Mapper (TM) images acquired in 2000, Landsat 7 Enhanced Thematic Mapper (ETM+) images acquired in 2008, and Landsat 8 Operational Land Imager (OLI) images acquired in 2019, to extract the glacier boundaries (Table 1). Several data gaps occurred in the ETM + images due to failure of the Landsat 7 Scan Line Corrector, which were dealt with using the Scan Line Corrector Gap-Fill Methodology of landsat gapfill.sav (Scaramuzza et al., 2004; Wang et al., 2020b). In addition, a DEM of the glacier surface in 2000 was derived from the Shuttle Radar Topography Mission (SRTM) DEM (https://lta.cr.usgs.gov/SRTM1Arc), which is regarded to provide the basic elevation information of the glacier surface. In this study, an SRTM DEM with a spatial resolution of 30 m was used. The Advanced Land Operation Satellite World 3-D (ALOS 3D) 30 m digital surface model (DSM) of the study area was produced by the Japan Aerospace Exploration Agency (JAXA) using several images from 2005 to 2010 (https://www.eorc.jaxa.jp/ALOS/en/aw3d30/index.htm). The time span of the DSM on YSM was denoted as 2008 to conveniently calculate the mass balance of the glacier. Therefore, the DSM was used as the DEM of the glacier in 2008 due to a lack of obstructions on the glacier. Finally, the SRTM DEM for 2000 and the ALOS DEM for 2008, which had the same spatial resolution of 30 m were used to assess the changes in the elevation of the glacier surface.

Mass balance stakes and snow pits have been widely used to measure and record in situ glacier mass balance (Østrem and Brugman, 1991; Zemp et al., 2013). This method usually has enough flexibility to resolve the details of small glaciers and allows researchers to measure short-term glacial variations (Ai et al., 2019). To further understand the melting and accumulation processes involved in the mass balance of the BRG1, we adjusted the observation network of the glacier mass balance in 2019 relying on data from the Yulong Snow Mountain Glacier and the Environment Observation and Research Station. A total of 12 stakes and seven snow pits were set up on the surface of the glacier. The stakes were used to record glacial melting, and snow pits were mainly used to observe the glacier accumulation. In addition, we made observations every 7–10 days in the summer of 2019 to obtain high-frequency glacial melt records.

Precipitation and air temperature are the primary climatic factors controlling glacier mass balance; therefore, meteorological observations in glacial areas are an important part of studying changes in glaciers (Radić and Hock, 2014). In fact, several meteorological stations have been set up near the BRG1 since 2005 to monitor and record meteorological data, including the Lijiang meteorological station, the records for which began in 1951 (Wang et al., 2020a). The Lijiang meteorological station, near Lijiang city in Yunnan Province, China, was built by the China Meteorological Administration to record the climatic conditions in the Lijiang region. It is located in the southern region of the YSM and is located nearly 26 km from Lijiang city. The dataset was obtained from the China Meteorological Data Service Center (http://data.cma.cn/).

Accumulation, c, refers to the solid water collected by glaciers, including snowfall, condensation, and refrozen rainwater, on the glacier surface, as well as any accumulation of blown snow and avalanche snow redistributed by wind and gravity. Ablation, a, refers to the expended portion of the glacial solid water, including runoff from the melting of ice and snow, evaporation (sublimation), the disintegration of ice bodies, wind-blown snow transported off of the glacier, and snow lost due to avalanches. On cold glaciers, part of the meltwater is refrozen in the granular snow layer, on the ice surfaces, and in crevasses. This process is generally referred to as internal recharge, and this part of the meltwater cannot strictly be counted as ablation because it does not leave the glacier. However, since internal recharge is not easy to observe and calculate, all of the melted snow and ice are generally counted as ablation in general mass balance observations, and calculations are conducted using meteorological methods.

Accumulation minus ablation is defined as the mass balance, that is, b in Eq. 5: b = c - a . (5)

1) Contour line method of glaciology

The mass balance contour line method was used to calculate the net mass balance (Østrem and Brugman, 1991). The net mass balance at two adjacent contour lines is defined as the specific net mass balance, and the annual net mass balance value for the entire glacier is calculated using Eq. 6: b n = ∑ 1 n s i ′ b i ′ / S , (6) where s ′ i and b ′ i are the projected areas between two adjacent contours and the average net balance, respectively; n is the number of areas between the contour lines, and b _n is the net mass balance of the entire glacier.

2) Geodetic method

The geodetic method has traditionally been used to calculate the glacial mass balance (Andreassen et al., 2002; Cogley, 2009; Andreassen et al., 2016). In this study, DEMs from SRTM and JAXA, with spatial resolutions of 30 m, were used in the geodetic mass balance estimation. Because these DEMs were obtained from two different sources using different image acquisition techniques, horizontal or vertical offsets may exist (Kumar et al., 2017). Thus, it was necessary to compare the two phases of the glacial DEMs (Nuth and Kääb, 2011). We used the geodetic method to indirectly calculate the mass balance of the BRG1, in which two of the DEMs of the glacial surface were subtracted to calculate the changes in the volume, which was then converted to the mass balance using a density conversion (Cogley et al., 2011). The change in the total volume Δ V was determined by summing the change in the elevation Δ h i at an individual pixel r during a specified time period: Δ V = ∑ i = 1 N r Δ h i , (7) where N is the number of pixels covering the glacier at its maximum extent and r is the pixel size (30 m × 30 m).

The calculated change in the volume can be converted to a geodetic mass balance (m w.e.) using Eq. 8. B = Δ V S ¯ ⋅   ρ ρ w a t e r = Δ V S t 1 + S t 2 2   ⋅ ρ ρ w a t e r , (8) where ρ is the average density of 850 ± 60 kg m⁻³, S is the mean glacier area of the BRG1, and t ₁ and t ₂ are the two acquisition dates.

(3) Profile curve method

The principle of the profile curve method is similar to that of the contour line method. Due to objective reasons, sometimes the ground mass balance stakes and snow pits collapse or are missing, so less observation information than is desirable is available. For example, only a few limited points were documented within a certain elevation range. Therefore, the curve pattern of the mass balance can be plotted using only a few observation sites, and then, the mass balance of the corresponding elevation zone can be calculated by combining the results with the glacier area of the elevation zone (Andreassen et al., 2016). That is, the regional mass balance between adjacent contour lines within each elevation interval (50 m or 100 m) was integrated using the profile line method to calculate the total amount of accumulation and ablation for each elevation interval in order to obtain the glacier mass balance. Therefore, in this study, the profile curve method for the mass balance calculation was tested on the BRG1 for the first time.

The SRTM and ALOS DEMs were used to extract the elevation of the ice surface of the BRG1 during two phases in order to analyze the changes in the elevation of the glacier surface in this region. Glaciers tend to be in constant motion due to their own physical properties and gravity, but the bare rock surrounding a glacier remains comparatively unchanged. To ensure a good spatial match between the two phases of the DEMs, 13 fixed points were selected, mainly in the mountainous areas or peaks without glaciers, to validate the consistency of the two phases in the DEM space (Figure 1). As is shown in Figure 2, the elevations of the selected points that did not change position were statistically tested. The results show that r and R ² values were 0.98 and 0.96, respectively, and they passed the significance test (p < 0.0001). We also found that some differences in the fixed points were occurred, which was mainly due to the geographic coordinate system and resampling of the elevation raster. Therefore, the elevations of these selected validation points were consistent. That is, the SRTM and ALOS DEMs can be used to assess the changes in the elevation of the glacier surface due to their good spatial matching.

To understand the mass balance of the BRG1, the annual mass balance estimates from 2000 to 2010 were obtained (Wang et al., 2020b). As is shown in Figure 3, the glacier mass balance fluctuated and changed during the period of 2000–2010. The average mass balance was -0.63 m w.e., and it exhibited a significant decreasing trend of 0.16 m w.e./a during this period. The mass balance of the BRG1 was initially positive and reached a maximum of 0.92 m w.e. in 2002. The mass balance of the BRG1 was negative from 2003 to 2010, that is, the glacier underwent mass loss after 2002. The largest mass loss occurred in 2003, with a mass balance of -1.77 m w.e.

To deepen our understanding of the glacial melting and accumulation processes, an enhanced observation network was established on 27 July 2019, and in situ data were recorded weekly until 15 September 2019. During this period, the glacier underwent strong melting and mass loss. The total glacial melting at the observation positions ranged from 1.60 m w.e. to 2.00 m w.e., with an average glacier mass balance of −1.75 m w.e. (Figure 4A). The mean daily glacial melting was at 35 mm w.e., and the strongest glacial melting occurred in August. The mass balance of the entire glacier was −1.41 m w.e. In addition, the annual mass balance of the glacier was also calculated based on the observational records until September 2020. From September 2019 to the beginning of September 2020, the winter mass balance ranged from 1.09 m w.e. to 3.51 m w.e., with an average accumulation of 2.06 m w.e. (Figure 4B). The annual net mass balance at the observation points ranged from −0.88 m w.e. to −3.44 m w.e., with an annual mass balance of −1.95 m w.e. The annual mass balance increased significantly (i.e., decreased melting) with increasing elevation above 4,700 m a.s.l. Finally, using the contour line method, the annual mass balance of the entire glacier was determined to be −1.31 m w.e. in 2019/2020.

The spatial variation in the glacier mass balance can clearly reflect the spatial variation trend of a glacier. To understand the spatial mass balance of the BRG1, the spatial mass balance was calculated using the geodetic method based on the SRTM and ALOS DEMs. A contour line map of the mean annual mass balance was drawn for the 2000–2009 period (Figure 5A). The pattern of the glacier mass balance was not the same as that of other glaciers. In general, the largest glacier accumulation occurred in the firn zone on the upper surface of the glacier at the highest elevation. However, the largest accumulation on the BRG1 occurred on the glacier surface in the central glacier zone, which was significantly affected by hill shadows. The mass balance at the glacier front was less than -3.0 m w.e., indicating very significant mass loss. To further validate and understand the pattern of the glacier mass balance, the contour lines of the glacier mass balance in summer 2019 were also drawn using mass balance records with high-frequency observations (Figure 5B). The spatial characteristics of the contour line map were similar to those of the annual mass balance during the period of 2000–2009. Therefore, we found that the glacier mass balance increased with increasing elevation, while it decreased at elevations of greater than ∼4,800 m a.s.l. Although the spatial patterns of the glacier mass balance were very similar, a significant difference was observed. For example, accumulation did not occur on the glacier surface during the period of 2000–2009, and it was completely absent during the summer of 2019, that is, the glacier mass loss was more significant than before. In addition, the spatial pattern did not remain stable indefinitely, which was predominantly dependent upon the glacier size and surface features, the surrounding terrain, and the local climate.

The mass balance at different elevations on the glacier surface was related to the elevation and exhibited a pattern containing specific elevation zones. To explore the spatial pattern of the mass balance and elevation on the glacier surface, the surface was divided into 11 zones using an elevation interval of 100 m. The extracted glacier area was 1.47 km² based on Landsat 7 ETM + images acquired in 2008. Then, this area was divided into the corresponding intervals (Figure 6). We found that the glacier area was mainly concentrated in the elevation range of 4,600 m a.s.l. to 5,000 m a.s.l., and the largest glacier area was located at 4,800 m a.s.l. First, the mass balance of the individual grid cells was statistically analyzed. In addition, the mass balance was statistically counted using an interval of 100 m. Next, the average mass balance in each elevation band was calculated. A profile curve of the change in the mass balance with elevation was obtained (Figure 6). In addition, we also found that the pattern of the glacier mass balance was similar to the pattern of the contour line map. The mass balance of the entire glacier was obtained according to the weighted sum of the average mass balance along the profile curve multiplied by the area in each corresponding elevation band. According to the final calculations, the mean annual mass balance of the BRG1 from 2000 to 2010 calculated using the contour method was -0.99 m w.e., and that calculated using the profile curve method was -1.01 m w.e. The mean annual mass balance of the glacier was calculated to be -1.18 m w.e. during the period of 2000–2010 using the geodetic method. The results obtained using the three methods were consistent. In addition, the equilibrium line elevation was ∼4,800 m a.s.l. according to the profile curve of the mass balance of the BRG1 during the period of 2000–2010.

To understand the local climate change on Yulong Snow Mountain, air temperature and precipitation were analyzed using the records from the Lijiang meteorological station. The annual air temperature ranged from 11.8 to 14.2°C during the period of 1951–2019 (Figure 7A). The mean annual air temperature was 12.9°C in the study area. A significant increasing trend (0.20°C/decade, p < 0.0001) occurred during this period. The annual precipitation ranged from 648 to 1,283 mm during the period of 1951–2019, exhibiting strong fluctuations with time (Figure 7B). Although the annual precipitation exhibited an increasing trend of 4 mm/decade, the trend was not significant (p > 0.1). Therefore, the amount of precipitation did not change significantly or fluctuate significantly during the study period, with a mean annual precipitation of 954 mm. An increase in air temperature tends to accelerate glacial melting.

The glacial melting was mainly controlled by the local air temperature. Based on the daily records for Lijiang meteorological station from July 27 to 15 September 2019, the change in the daily air temperature at an elevation of 4,800 m a.s.l. was calculated, with a rate of decrease of 0.6°C/100 m. As is shown in Figure 8, the daily air temperature ranged from 0.8 to 7.3°C during the observation period in the summer of 2019, with an average daily air temperature of 4.0°C. Precipitation occurred on 35 of 51 days. Rainy weather also significantly affected the glacial melting process (Figure 8). In particular, glacial melting was very significant on days with a warm air temperature and without precipitation, such as during the period of August 16–26 and September 5–15. However, from August 26 to September 5, precipitation occurred frequently. The daily air temperature was low on rainy days, and the lowest daily air temperature was 0.8°C. In addition, fresh snow was observed during the fieldwork conducted in this specific time period. Accordingly, the glacial melting was relatively weak on rainy days compared with that on other days. In general, the melting of the BRG1 mainly resulted from high air temperatures, but it was also significantly affected by precipitation.

Glacial melting varied across the Qinghai–Tibet Plateau based on variations in the local climatic circulation and the terrain conditions. To understand the different impacts of climate change on the behavior of glaciers, several glaciers have been observed since 2000, including Parlung River Glaciers Nos. 04, 10, 12, 94, and 390, as well as the Demula Glacier (Yang et al., 2011; Yang et al., 2013; Zhu et al., 2018; Wang et al., 2021). The average annual mass balances of the Demula Glacier, Parlung River Glacier Nos. 04, 10, 12, and 94, and the BRG1 were -1.02 m w.e. for 2007–2010, -0.37 m w.e. for 2006–2007, -0.78 m w.e. for 2006–2009, -1.70 m w.e. for 2006–2010, -0.92 m w.e. for 2006–2010, and -1.16 m w.e. for 2008–2010 (Figure 9). The average annual mass balance of the BRG1 during the period of 2000–2010 was calculated to be -1.18 m w.e., -0.99 m w.e., and -1.01 m w.e. using the geodetic, profile curves, and contour line methods, respectively. The mean annual mass–balance of the BRG1 was -1.56 m w.e. from 2011 to 2019, and that of Parlung River Glacier No. 94 was -0.95 m w.e. from 2011 to 2018. In addition, the mass balance of the Demula Glacier ranged from -1.67 to 0.17 m w.e. during the period of 2006–2010, that of the Parlung River Glacier No. 94 ranged from -1.93 m w.e. to 0.16 m w.e. during the period of 2006–2018, and that of the BRG1 ranged from -1.87 m w.e. to -0.91 m w.e. during the period of 2008–2019. Thus, the mass balance results for the BRG1 obtained in this study using the different calculation methods are reliable. In addition, the mass balance fluctuations of these glaciers were similar and exhibited an increased melting trend with time.

We acknowledge that some differences exist among the results obtained using the contour line, geodetic, and profile curve methods due to their uncertainties (Zemp et al., 2013). The contour line method is a glaciological method, and its uncertainty is mainly derived from the uncertainties of the point measurements, spatial integration, and glacier reference area. The uncertainty of the geodetic method is mainly due to the uncertainties of the digital terrain models (DTMs) and the density conversion used to convert the glacier surface area to the water equivalent. The results obtained using the profile curve method were determined from the profile line of the mass balance with a change in elevation. In general, continued geodetic surveys every 10 years are needed to measure the overall changes and provide data for calibrating the glaciological mass balance (Zemp et al., 2013; Andreassen et al., 2016). For example, Oerlemans (2013) estimated an even higher dissipative melting rate of −0.23 m w.e. per year in Nigardsbreen, which indicated that the accumulative mass balance was very sensitive to systematic biases and generic differences between the geodetic and glaciological methods. Andreassen et al. (2016) discussed the differences and uncertainties of the geodetic and glaciological methods, and they reanalyzed and calibrated 10 glaciers with long-term mass balance series using the profile method. The glaciological method measures the glacier surface mass balance, while the geodetic method measures the glacier surface, internal, and basal mass balances. In addition, relatively small mountain glaciers with negative cumulative balances are easier to measure correctly than maritime glaciers (Andreassen et al., 2016). The accumulation and ablation patterns are needed for use in the statistical analysis to reduce the amount of field work, expenses, and personnel required (Rasmussen and Andreassen, 2005). The profile curves provided some important information about the mass changes related to ignoring the areas of the glacier surface for which no observation data were available, for example, ice falls, crevasses, and unusual phenomena in the firn/accumulation zone. The combination of the profile curve and geodetic methods provides more accurate results and can be used to calculate the mass balances of maritime glaciers (Andreassen et al., 2015; Andreassen et al., 2016). These methods can also be used to calculate the mass balance of the BRG1 in the future.