Summertime surface mass balance and energy balance of Urumqi Glacier No. 1, Chinese Tien Shan, modeled by linking COSIMA and in-situ measured meteorological records

To get a better overview of atmosphere-driven mass changes at Urumqi Glacier No.1, Chinese Tien Shan, the surface energy budget and mass balance is modeled by linking the COupled Snowpack and Ice surface energy and MAss balance model (COSIMA) with in-situ measured meteorological records for the ablation period in 2018. The COSIMA is calibrated by manual optimization and the modeled results agree well with the in-situ surface temperature, snow height and seasonal mass balance. Our results reveal that Urumqi Glacier No.1 experienced a significant mass loss, with an average value of − 0.77 m w.e. over the ablation period 2018. The surface energy budget components can be classified into two categories: radiation (shortwave and longwave) and turbulent fluxes. Surface melt and solid precipitation were dominated components of mass balance. The COSIMA can reproduce the glaciological mass balance compared with other models. Sensitivity analysis showed that the mass balance was more sensitive to the temperature than precipitation, and mass loss caused by temperature increase of 1 K needed to be compensated by at least 40% precipitation increase. Air temperature during the ablation period was more important than annual precipitation in controlling mass balance changes on Urumqi Glacier No. 1. These findings will enhance our understanding of the mechanisms underlying mass balance processes of ablation period and their contribution to the acceleration of glacier retreat in Tien Shan.


Introduction
The glaciers in the Tien Shan provide the sources of rivers that supply water to billions of people (Farinotti et al. 2015). Several studies have found the glaciers in Tien Shan represented an accelerated melting tendency recently Miles et al. 2021). The majority of glaciers have been receding and shrinking (Yao et al. 2012), which could endanger the water resource in the downstream region (Immerzeel et al. 2010).
Given the importance of glaciers in the Tien Shan, the quantitative understanding of the physical relationship between glaciers and climate is necessary to help predicting the responses of glaciers to climate change and its potential impacts. Most of previous studies on the response of glacier mass balance to climate change in the Tien Shan have focused on the effect of annual air temperature change (Farinotti et al. 2015;Pieczonka et al. 2015;Wang et al. 2016Wang et al. , 2017Sakai and Fujita 2017). However, the air temperature in the Tien Shan shows a large spatial discrepancy in seasonal variations (Wang et al. 2020a). Few studies on the Tien Shan quantify the effects of seasonal air temperature changes on changes in glacier mass balance. Furthermore, glaciers in the western and eastern Tien Shan exhibit strong heterogeneous changes due to their different climatic settings (Wang et al. 2015).
Due to the lack of in-situ glacio-meteorological observations in different regions of the Tien Shan, our understanding of the relationship between glaciers and climate change and glacier melt water runoff patterns are limited. The relationship between glaciers and climate change can be resolved by studying glacier energy and mass balance, which has been conducted in detail on many mountain glaciers and ice sheets around the world (Bintanja et al. 1997;Oerlemans and Klok 2002;Kaser et al. 2004). The glacier energy and mass balance modeling has been systematically studied on the individual Tibetan Plateau glaciers (Fujita and Ageta 2000;Yang et al. 2011;Mölg et al. 2012;Huintjes et al. 2015a, b;Sun et al. 2018;Li et al. 2018Li et al. , 2019Zhu et al. 2020). These studies mainly focus on the impact of meteorological factors or regional-scale atmospheric circulation on the mass balance, the reconstruction of the glacier mass balance has also accomplished using an energy and mass balance model. Few systematic studies have been conducted on Tien Shan glaciers (Kang et al. 1992;Han et al. 2005), and glacier mass balance and meteorological data are currently only available from the Urumqi Glacier No. 1.
Compared with the previous studies in the Urumqi Glacier No.1, hourly high-resolution surface mass balance modeling was only available until 2012 (Che et al. 2019;Li 2020), but the parameterizations of subsurface energy and mass fluxes were not incorporated or poorly resolved by these commonly models. Huintjes et al. (2015a, b) has applied COupled Snowpack and Ice surface energy and MAss balance model (COSIMA) with observations from Automatic Weather Station (AWS) or High Asia Refined (HAR) reanalysis data as atmospheric input to investigate glacier-wide surface energy budget and mass balance processes. Subsurface processes can be well resolved by COSIMA. In general, parameterization scheme within COSIMA makes the calculation of mass balance more reasonable. However, COSIMA has not been tested in the Tien Shan using hourly AWS observation data as input.
Taking this into account, COSIMA was applied to the Urumqi Glacier No. 1, Chinese Tien Shan, forcing by hourly meteorological data collected at an AWS on the glacier surface. The energy fluxes and mass balance at the glacier surface and in the subsurface layers were calculated and a comparison with different types of field in-situ measurements was conducted to validate model performances. Our specific aims were: (1) to understand the temporal variability of meteorological conditions and characteristics of glacierwide energy budget; (2) to quantitatively characterize the glacier-wide mass balance components and their profiles; and (3) to analyze the climate sensitivity of the Urumqi Glacier No.1 and to identify the control of atmospheric variables on the mass balance of Urumqi Glacier No.1.

Study site
Urumqi Glacier No. 1 (43° 06′ N, 86° 49′ E) is a typical continental, valley-type glacier, covering an area of 1.52 km 2 . It consists of two independent small glaciers: the east branch (EB) and the west branch (WB), and its altitude ranges between 3743 and 4267 m a.s.l. Urumqi Glacier No.1 is situated on the northern slope of Tianger Summit II in the eastern Tien Shan (Fig. 1). It is situated in the continental climate zone. According to Liu et al. (1992), this glacier is a typical summer-accumulation-type glacier because both accumulation and ablation take place simultaneously between June and September. The region is mainly controlled by three dynamical factors: A westerly jet in the upper troposphere, Siberian anticyclonic circulation, and westerly cyclonic disturbance. According to the observation data of the Daxigou Meteorological Station (3539 m a.s.l.), 3 km southeast of Urumqi Glacier No. 1, the average annual temperature was approximately − 4.9 °C, and the annual precipitation was 465 mm during 1959-2018. The average annual temperature increased by a rate of 0.29 °C (10a) -1 and the annual precipitation increased gradually at a rate of 20 mm (10a) -1 based on the linear analyses shown in Fig. 2. The main ablation period is from May to August, and the precipitation mainly comes from the water vapor carried by the west wind. The precipitation from May to August accounts for 78% of the annual amount, with the main type of solid fall, such as snow, hail and sleet.

COSIMA descriptions
COSIMA is a physical surface energy and mass balance model that assumes mass conservation in the snowpack. Thus, COSIMA couples a surface energy budget to a multilevel subsurface model, and it can calculate energy and mass balance for different components, including energy budget, meltwater percolation, refreezing and densification (Huintjes et al. 2015a). Now, COSIMA is available online as an opensource software and used for this study (https:// bitbu cket. org/ glaci ermod el/ cosima; last accessed on 17 June 2021; Huintjes et al. 2015a). It is in detail presented by Huintjes et al. (2015a, b) and Wang et al. (2020b). COSIMA calculates the mass balance of the surface melt, subsurface melt, sublimation, deposition, solid precipitation, and refreezing at hourly time step. The surface energy balance equation governs the calculation of the mass fluxes: where SW in , LW in and LW out represents incoming shortwave radiation, incoming longwave radiation and outgoing longwave radiation, respectively. SW in (1 − α) is net incoming shortwave radiation (SW net ). α is the surface albedo. Q sens and Q lat are turbulent sensible and latent heat flux, respectively. Q G is the ground heat flux which consists of fluxes of heat conduction (Q C ) and penetrating shortwave radiation (Q PS ). Q PS is equal to SW net (1 − ζ). For ζ we take the values of 0.8 for ice, and 0.9 for snow, respectively, as determined by Bintanja and van den Broeke (1995). The heat flux from liquid precipitation is neglected within COSIMA according to Huintjes et al. (2015a, b). The flux towards the surface is considered positive, while the flux from the surface has a negative sign. Sublimation occurs if Q lat is negative. Surface melt requires the surface temperature to be at the melting point (273.15 K) and the resulting energy flux F to be positive. In this case, F is equal to the heat flux for snow/ice melting (Q melt ). If surface temperature is below the melting point, F is zero. To calculate F, COSIMA computes the atmospheric energy fluxes at the glacier surface from meteorological variables and the subsurface temperature profile inside the glacier to determine Q G . The Q lat and F are converted into mass fluxes that contribute to the surface mass balance of the glacier together with calculated accumulation (Huintjes et al. 2015a).
The double critical temperature index method is used to separate solid and liquid precipitation in COSIMA. The precipitation is liquid precipitation when the air temperature is higher than 6 ℃, whilst it is solid precipitation when lower than 0 ℃ in our study region. The air temperature between the two critical air temperatures is calculated by linear interpolation (Jia et al. 2020). In COSIMA, the albedo is calculated using the parameterization of Oerlemans and Knap (1998). The ice albedo variation in the Urumqi Glacier No. 1 has a range from 0.06 to 0.44 due to topographic effects and light-absorbing impurities (Yue et al. 2017). Manual tests and inspection of albedo in Landsat-8 data showed that the albedo decreased from a fresh snow of 0.83 to a firn albedo of 0.58, with an aging factor of 1.1 day after the snowfall events. The albedo drops to 0.35 when glacier surface is without snow layer (Yue et al. 2017). The depth scaling factor of the albedo is 3 cm after manual tests ( Table 2). The vertical profiles of subsurface temperature and density are key for Q G . The subsurface temperature is solved by the thermodynamic thermal equation (Anderson 1976). And subsurface density is calculated through an empirical relation following Herron and Langway (1980). The initial temperature profile is defined by the available subsurface measurements and linear interpolation with depth. The initial snow density profile is calculated from a linear interpolation between 250 kg m −3 of the uppermost snow layer and 550 kg m −3 from the bottom snow layer.
The distributed SW in including effects from both terrain shading and cloud cover were calculated based on the radiation model within COSIMA and the detailed information about the parameters are given in Kumar et al. (1997). LW in and LW out are obtained by the Stefan-Boltzmann law. For LW in , atmospheric emission rates ε is calculated following Klok and Oerlemans (2002). To calculate the turbulent heat fluxes, COSIMA assumes that the surface roughness length varies with the duration of the fresh snow to firn snow (Mölg et al. 2009). Surface roughness length ranged from a fresh snow value of 0.24 mm (Gromke et al. 2011) increases linearly to a firn snow value of 4 mm (Brock et al. 2006). If the glacier surface is without snow, the bare ice is assumed to be 1.7 mm (Cullen et al. 2007). The correction of turbulent fluxes under stable conditions is based on the Richardson number (Braithwaite 1995). The Q C is determined by the temperature difference between the surface layer and the two top subsurface layers, depending on the thermal conductivity of the medium (ice or snow). The thermal conductivity is calculated from the subsurface density (ρ, in kg m −3 ) (Anderson 1976). The Q PS is calculated according to Bintanja and van den Broeke (1995). At the top level of the model, the proportion of SW net absorbed in surface layer reaches an exponentially decreasing flux to the bottom layer and increases the subsurface temperature. For the fraction of the SW net absorption and extinction coefficients, we take the values of 0.8 and 2.5 for ice, and 0.9 and 17.1 for snow, respectively, as determined by Bintanja and van den Broeke (1995). The detailed information of the COSIMA is presented in the Appendix A.

COSIMA configuration
The final COSIMA run was initiated from 00:00 a.m. 1 May 2018. The analyzed COSIMA output also started from this time. At this time there was a little fresh snow observed at the glacier surface facilitating the initialization process of the model. The input meteorological series of hourly air temperature, relative humidity, precipitation and incoming shortwave radiation were used for the timespan between 00:00 a.m. 1 May and 00:00 a.m. 31 August 2018. The input meteorological data and their processing were shown in Sect. 4. For the COSIMA modeling, the bottom temperature was 266.15 K as mean value for the ablation period at the AWS1. Therefore, the surface energy budget and mass balance from AWS1 is conducted for hourly values.

Mass changes assessment
In this paper, by analyzing the interannual changes of temperature (σ T ) and annual precipitation (C v ) during the ablation period, the meteorological factors that control the mass balance of glaciers can be quantitatively determined. The air temperature in ablation period and annual precipitation were used to inter-annual variability of air temperature in ablation period (σ T ) and annual precipitation (C v ), respectively. This method has been successfully applied to the Chhota Shigri glacier in the western Himalayas, Parlung Glacier No. 94 in the southeast Tibetan Plateau, Muji Glacier in the northeastern Pamir Azam et al. 2014;Zhu et al. 2020). σ T refers to the standard deviation of air temperature during ablation period, C v is the ratio of the standard deviation to the annual precipitation. σ T or C v has been used to indicate interannual changes in ablation period air temperature and annual precipitation, respectively. Then COSIMA was run by virtue of modifying air temperature or precipitation to obtain respective mass balance change. In this study, the forcing data include the air temperature and annual precipitation of Daxigou meteorological station and the detailed information has been described in Sect. 4.1.

Input meteorological data
An automatic weather station (AWS1; Table 1) was set up on a relatively flat surface (slope < 5°) near the equilibrium line altitude (ELA) of the Urumqi Glacier No.1 from 29 April 2018 (43° 06′ 22″ N, 86° 48′ 23″ E; 4025 m a.s.l.; Fig. 1a, c, d). AWS2 was installed on the glacier terminal moraine in 2011 at an altitude of 3835 m a.s.l. providing the same meteorological datasets as AWS1 and stored in the Campbell CR1000 data logger as 30-min mean values. Table 1 shows the technical specifications of all sensors in AWS1 and AWS2. Detailed information of both AWSs can be found in Wang et al. (2020b). The air temperature (T), relative humidity (RH), precipitation (P r ) and incoming shortwave radiation (SW in ) recorded by AWS1 and AWS2, are required for this study. Surface temperature recorded by AWS1 for evaluation of model performance. Unfortunately, there was a data gap spanning from 12:00 p.m. 20 August to 12:00 p.m. 22 August 2018 due to harsh environmental conditions. The following part will offer the process for filling in data gap. The annual precipitation and the mean air temperature in the ablation period of the Daxigou meteorological station (43° 06′ N, 86° 50′ E) during 1959-2018, located at an altitude of 3539 m a.s.l., were used to assess how climate factors control the change in mass balance.
Due to frequent snowfall events, the top of the radiation sensor often forms snow and frost, which may affect SW in measurements, causing defects in the data. Correction for SW in must be carried out and the method was proposed by van den Broeke et al. (2004). Since the SW in correlation coefficients of AWS1 and AWS2 reached 0.95, the gap in the SW in was filled by establishing the linear regression between AWS1 and AWS2.
Both forms of rain and snow were recorded using a Geonor T-200B weighing bucket precipitation gauge. Considering the impact of wind speed on solid precipitation, the raw precipitation data was corrected by He et al. (2009). To prepare the continuous precipitation for the 2018 ablation period, the precipitation recorded by AWS1 and AWS2 were used to fill the data gap and obtain a vertical gradient of precipitation. And the air temperature data gap has also been filled by lapse rate for air temperature. The vertical gradient of precipitation and the lapse rate of air temperature were 1.6 mm (100 m) −1 and − 1.1 ℃ (100 m) −1 , respectively. The air pressure lapse rate was − 0.067 hPa m −1 and the relative humidity lapse rate was 0.022% m −1 during the modeling period. Cloud cover (N) is a forcing variable but cannot be observed directly and is estimated using a method described by Favier et al. (2004).
where SW TOA represents the solar radiation at the top of atmosphere with the unit of W m −2 , which can be calculated according to the method of Duffie et al. (1991). Hock et al. (2005) found that the wind speed depended on the meteorological station location, and it was impossible to get a general scheme to quantitatively analyze the wind speed on the glacier. Linear regression analysis of wind speed from AWS1 and AWS2 yielded hourly regression parameters that could then be used to reconstruct wind speed during ablation period at the AWS1. This method was previously applied by Yang et al. (2016). The interpolation of relative humidity was also based on the method described by Hock et al. (2005), and the relative humidity varied relatively little between AWS1 and AWS2. Therefore, the linear regression analysis of relative humidity between AWS1 and AWS2 was also conducted and the relative humidity in the same period have similar change trend (Correlation coefficient, R = 0.85; Root Mean Square Error, RMSE = 0.23%), so that it can be extended to the overall glacier surface. Moreover, the isothermal atmospheric pressure formula was used to obtain continuous air pressure (Zhou et al. 1997). Anyway, the gap time is just a short period and the effect on modeled results adopted by the methods mentioned above is limited.

Topography and glacier mask
A Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model Version 3 (GDEM 003) with a 30 m spatial resolution, provided by the National Aeronautics and Space Administration (NASA) Land Processes Distributed Active Archive Center Table 1 Information and parameters of AWSs sensor and precipitation gauge used in this study T air temperature, P r precipitation, P air pressure, RH relative humidity, u wind speed, SW in and SW out incoming and outgoing shortwave radiation, LW in and LW out incoming and outgoing longwave radiation, T s surface temperature

Instrument
Sensors Company Parameter Accuracy Height  . Two flights were carried out with the flying height of 120 m for each flight. The overlap of each flight was 80%, and the camera pitch angle was 70°. The UAV-DEM was generated with a spatial resolution of 0.36 m based on Pix4D Mapper software. Glacier outlines were extracted from the UAV-DEM and used to split the terrain into glacier and non-glaciated area (Fig. 1a), which was used to modeling energy and mass balance for the Urumqi Glacier No. 1.

Mass balance measurement
Glacier mass balance of the Urumqi Glacier No. 1 was measured following the glaciological method using ablation stakes or snow pits. Thirty-six stakes were evenly distributed on the glacier at different elevation bands by a steam drill (Fig. 1a). The ablation stakes were measured at 29 April 2018 and 1 September 2018, respectively. The measurement items include the vertical height of stakes above the glacier surface, superimposed ice thickness and the density and thickness of each snow/firn layer at an individual ablation stake. Then all the measurement items were converted to water equivalent (w.e.) using the single-point density values from in-situ measurements. The summer mass balance at each stake was then calculated based on these data. Point values were calculated for individual elevation intervals. The mean specific mass balance, expressed as a function of elevation, was calculated by the interpolation and extrapolation of the point results to the whole glacier surface in order to compare with modeled mass balance.

Meteorological conditions during ablation period
The daily mean air temperature was approximately − 3.2 ℃ throughout the ablation period 2018. Daily mean air temperature fluctuated dramatically for the first few days of ablation period (1 May-8 June) and stayed comparatively low. The daily mean air temperature almost remained above the melting point between 9 June and 28 August (Fig. 3a). The daily mean relative humidity and air temperature exhibited similar variations. Daily mean relative humidity was about 69% during ablation period indicating a fairly humid condition (Fig. 3b). Daily mean relative humidity in May was relatively low. High moisture was particularly evident in the mature period of the rainy season (June to August).
The daily mean wind speed typically fell between 0.6 m s −1 and 5.3 m s −1 , with an average of 2.3 m s −1 (Fig. 3c). Relatively high wind speed occurred before 27 June and then it became lower. The pressure during the ablation period was relatively stable, with an average of 625 hPa (Fig. 3d). Moreover, the cloud cover declined until 22 May and then increased (Fig. 3e). The total precipitation during our observation period was approximately 692 mm, and mostly  (e), precipitation (f) and incoming shortwave radiation (g) at AWS1 site. The black dotted lines in a-d, indicate the mean value occurred from June to August (Fig. 3f). The incoming shortwave radiation (SW in ) is the important model forcing. SW in recorded in the AWS1 exhibited obvious diurnal variation ( Fig. 3g) and the average value of SW in during the ablation period was 224.3 W m −2 .

Model optimization and evaluation
This study applied COSIMA with an hourly time resolution and 30 m spatial resolution to present mass balance and energy budget for the Urumqi Glacier No.1 in the ablation period 2018. The surface temperature, snow height at AWS1, and the in-situ measured mass balance data of Urumqi Glacier No. 1 were used to model optimization and evaluation.
The parameter combinations of COSIMA should be carefully selected to obtain an optimal model performance. Since the modeling period is ablation period 2018, a manual trial and error method outperforms an automated Monte Carlo method for long-term documented modeled mass balance (Mölg et al. 2012;Cullen et al. 2014). A manual tuning of model parameters provides a better understanding of the simulation results in response to configuration changes. A list of parameterization options available in COSIMA is presented in Table 5 in Appendix B. The range of values tested and optimized values for COSIMA run is shown in Table 2. Input data from the AWS1 during the period 00:00 a.m. 1 May to 00:00 a.m. 31 August were used for optimization. In the next step, the model is applied to Urumqi Glacier No. 1 and the results are also evaluated. Given the possible range of COSIMA parameters in Table 2, a manual trial and error approach was used to obtain the optimal value of the COSIMA modeling. Afterwards, the model evaluation of COSIMA is conducted for the ablation period (from 1 May to 31 August 2018). The visual comparison is presented in Fig. 4 for AWS1. The comparison against the observation from AWS1 reveals a good agreement of COSIMA with the measured surface temperature (R = 0.71, RMSE = 2.61 K). In terms of snow height change, there is a significant offset in the AWS1. Especially the modeled snow height not captured after 7 August. This is partly attributed to the measurement error during the strongly ablation period. The sensor used to measure the distance between the surface and it also has a membrane, which is affected by the weather conditions of the ice surface causing this membrane to fail and produce false values for snow height. This was a generally problem for these instruments (Mölg et al. 2020). The effect of snowdrift on the snow height has been reported by Li et al. (2018) for Qiangtang Glacier No.1.
As strong ablation usually occurs on Urumqi Glacier No.1 during this period , there are few effects of snowdrift on the snow height. COSIMA overestimates the snow height at the location of AWS1. The RMSE between the daily measured and modeled snow height was 0.15 m over the observation period and R was 0.94. This study carried out the comparison between the measured and modeled summer cumulative mass balance in order to validate the robustness of COSIMA performance (Fig. 5). Negative values point to ablation and decreasing surface height, while the positive values indicate accumulation and increasing surface height. During the ablation period, summer cumulative mass balance calculated using the glaciological method and COSIMA was − 0.85 m w.e. and − 0.77 m w.e., respectively. This difference (i.e., 0.08 m w.e) is acceptable given the measured mass balance uncertainty. The modeled mass balance agreed well with the measured mass balance (R total = 0.70; RMSE total = 0.43 m w.e.). The R is 0.75 and the RMSE is 0.35 m w.e. for the WB, while the R is 0.67 and RMSE is 0.49 for the EB. Such a good performance regarding the mass balance modelling makes us confident in the COSIMA's ability to capture the mass balance on Urumqi Glacier No. 1. Therefore, COSIMA in its current configuration can offer a good estimate of surface  Figure 6 depicts the daily heat fluxes for snow/ice melting (Q melt ), net shortwave radiation (SW net ), incoming longwave radiation (LW in ), outgoing longwave radiation (LW out ), sensible heat flux (Q sens ), latent heat flux (Q lat ) and ground heat flux (Q G ). The main energy component causing the glacier mass loss were SW net (148.18 W m −2 ), LW in (241.84 W m −2 ) and Q sens (10.08 W m −2 ). The main energy expenditure were the LW out (− 283.11 W m −2 ), Q lat (− 5.93 W m −2 ) and Q G (− 3.41 W m −2 ).

Glacier-wide energy budget
Generally, SW net dominated the the heat flux for snow/ ice melting throughout the ablation period. SW net tended to increase before early June and then decrease, which may be related to increased precipitation. This can increase surface albedo, and finally result in the decreasing SW net . Due to the small air and surface temperatures, LW net was negative and acted as a heat sink. The record of LW net indicated three separate phases. Relatively larger negative LW net occurred before the early June and after the late August. Between these two periods, LW net was less negative. Affected by larger relative humidity, wind speed, and air temperature, the positive Q sens indicated that the heat transferred from the air to the glacier surface in the ablation period and then constantly provided energy gain for mass loss. Q lat was negative throughout the ablation period, mainly affected by the air-surface specific humidity difference, indicating that the sublimation process also existed throughout the ablation period and led to mass loss. This was because relatively lower temperatures and relative humidity led to specific humidity gradients signing from the glacier surface to the atmosphere. Q G provided a more minor heat sink compared with other heat fluxes depending on SW in and glacier surface temperature (Table 3). Figure 7 shows the spatial distribution of the modeled cumulative mass balance derived from COSIMA during the ablation period 2018. Spatially, the glacier mass balance ranged from 0.23 to − 2.12 m w.e. The significant mass loss was at the lower altitudes, while mass loss decreased at higher altitudes. The glacier ablation and accumulation were affected by the glacial topography. Spatial distributions of modeled mass balance were in line with the observation (Fig. 5). The spatial variation of the mass balance was influenced by slope and aspect. The cumulative mass balance was − 0.77 m w.e. during the ablation period, mass loss is mainly in the southeast of EB and WB. For the EB, the glacier melting in the western part was stronger than that in eastern part with the higher altitudes and steeper, because of the mountain shadow covered. In the same time, the cumulative mass

Seasonal specific mass balance components
Urumqi Glacier No. 1 is a typical summer-accumulation type glacier and the ablation and accumulation mainly take place simultaneously in the ablation period. The annual mass balance is significantly related to the summer mass balance (Liu and Han 1992;Li et al. 2011;Wang et al. 2016), therefore it is necessary and important to study the modeled mass balance over Urumqi Glacier No.1 during the ablation period. Moreover, the individual components of the glacier mass gain and loss were also derived from COSIMA. As formulated by COSIMA, the mass balance is composed of solid precipitation, sublimation, refreezing, surface melt, subsurface melt and deposition. Figure 8 shows the glacierwide daily mass balance components. Modeled mass balance continued to decrease until August 17. In general, the mass loss was significantly dominated by the surface melt with a value of − 0.73 m w.e., and then subsurface melt with a value of − 0.17 m w.e. (Table 4). The melting showed an increasing trend during May to July. Sublimation resulted in mass loss of 0.03 m w.e. during ablation period and mainly occurred in May. Although several studies have shown that the sublimation in the ablation period was smaller than in the cold season because of the dry environment, our result about sublimation was significantly smaller compared with other continental glaciers Zhu et al. 2018). Mass gain was derived from solid precipitation with a contribution of 0.19 m w.e., which mainly occurred from June to August. Refreezing contributed to mass gain could be neglected during the ablation period. Whether refreezing leads to mass gain during cold season provides an important internal mass storage component, needs to be further studied. Figure 9 compares the modeled and measured summer mass balance during the ablation period of 2018 as a function of altitude. The modeled mass balance well agreed with the measured mass balance, except for a discernible underestimation at the glacier terminus ranging from 3800 to 3900 m a.s.l. The change in the modeled mass balance with altitude during the ablation period can be divided into three stages (Fig. 9a): Mass balance tended to decrease and then increased in the altitude range 3800-3900 m a.s.l.; mass balance decreased gradually above 3900 m a.s.l.; mass balance was more positive in the altitude range 4150-4500 m a.s.l.

Altitudinal distribution of modeled mass balance
The ELA was at 4150 m a.s.l. The absolute values of surface melt decreased with increasing altitude (Fig. 9a), and ranged from − 1.6 m w.e. at the glacier terminus (3800 m a.s.l.) to 0.1 m w.e. at the top of the glacier (4250 m a.s.l.). Mass balance of EB and WB as a function of elevation in-depth analysis should be helpful to shed light on the details of mass balance. Glacier mass loss mainly occurred below 4150 m a.s.l. and the mass balance of EB was more negative than that of WB (Fig. 9b, c). Figure 9d shows the altitudinal distribution of the averaged mass balance components. The solid precipitation and glacier melt (including surface melt and subsurface melt) during the ablation period were prominent components of mass balance, whereas refreezing, deposition and sublimation were less dominant. In generally, surface melt significantly dominated the glacier mass loss, and then subsurface melt. Mass loss by way of sublimation in this study could be negligible. On the other hand, the mass gain was mainly due

Error analysis for COSIMA
The uncertainties of input parameters in the process of energy mass balance calculation can influence COSIMA results, which need to be considered. At Urumqi Glacier No. 1, uncertainties of input parameters are associated with the albedo, surface roughness length, vertical gradient for precipitation and lapse rate for air temperature (Table 6 in Appendix B). 10 scenarios for parameter sensitivity tests are set by perturbing one parameter and keeping other parameters unchanged. Figure 10 shows the evaluation of parameter uncertainty in COSIMA. Mass balance is less sensitive to the lapse rate of air temperature than the vertical gradient of precipitation. When the lapse rate for air temperature increased (or decreased) by 10%, a mass change of 0.001 (or 0.108) m w.e. was achieved, and the vertical gradient for precipitation ± 10% was produced less than 0.001 m w.e. of mass loss. The highest parameter sensitivity was related to ice albedo. The sensitivity for fresh snow albedo was higher than the firn snow albedo. When ice albedo increased (or decreased) by 10%, a mass change of − 0.140 (or 0.150) m w.e. occurred on Urumqi Glacier No. 1. Changes in the albedo depth scale also have significant effects on mass balance compared with albedo time scale. When the albedo depth scale increased (or decreased) by 10%, the modeled mass change was 0.013 (or − 0.015) m w.e. occurred on Urumqi Glacier No. 1. When albedo time scale increased (or decreased) by 10%, a mass change of − 0.003 (or 0.004) m w.e. was achieved. Moreover, we also evaluated sensitivity of surface roughness length, such as surface roughness length ice, firn snow and fresh snow, indicating that the effect of surface roughness length on mass balance can be neglected in this study, because the sensitivity of surface roughness length of ice, firn snow and fresh snow was very low.

Parametrization of the surface energy fluxes
The SW net is an important source of energy for glaciers. It is determined by SW in and albedo of the glacier surface within the COSIMA. The albedo values for snow and ice are variable due to grain size and form, liquid water content, topographic effects, impurities, etc. (Cuffey and Paterson 2010;Yue et al. 2017). Albedo parameterization scheme in COSIMA tries to reproduce albedo by introducing albedo of fresh snow, firn snow and ice, albedo time scale and albedo depth scale, which can solve exponentially decreases from fresh snow albedo to ice albedo (Oerlemans and Knap 1998). In Fig. 11 we show the comparison between the measured and modeled hourly albedo during the ablation period 2018. In generally, this albedo parametrizations in COSIMA were able to capture increases except early ablation period, and several studies had also used similar albedo parametrizations (e.g. Mölg et al. 2012;Huintjes et al. 2015b), nevertheless the measured albedo was much more variable than modeled albedo. The modeled albedo as shown in Fig. 11 was small in the early ablation period, which was probably related to the small snowfall events in the same period. However, measured albedo has a slight fluctuation compared with modeled albedo. Measured albedo values were calculated as a ratio of the outcoming shortwave radiation to the incoming shortwave radiation. Because the CNR4 radiation sensor has poor cosine response quality at solar zenith angles larger than 80°, the outcoming shortwave radiation flux was greater than the incoming shortwave radiation at certain moments early in the ablation period, Therefore, we deduced that the measured albedo at AWS1 site may be wrong, and the possible reason was the failure of the CNR4 radiation sensor. Most of the measured albedo increases are associated with snowfall events (as shown in Fig. 3f), but the double critical temperature index method is adopted to deduce it from total precipitation measured at AWS1 (see Sect. 3.1). Additionally, the faster the albedo decreases after snowfall events and the lower albedo time scale (1.1 day), indicating the faster snow metamorphism during the ablation period and may be associated with increased ambient temperature. Usually the LW net makes an important contribution to the energy exchange on the glacier surface. The LW net is often negative, this is because glacier surface is like as blackbody within COSIMA and atmosphere emissivity is often smaller than 1. Schaefer et al. (2020) has reported that the variability in emissivity cannot only be explained by the variability in the cloudiness and the relative humidity may influence emissivity. The uncertainties in the change of the cloud Fig. 9 Comparison between measured (black circle) and modeled (blue line) mass-balance elevation distribution for the ablation period at 50 m altitude intervals on the Urumqi Glacier No. 1 (a), East Branch (b), and West Branch (c). Grey horizontal bars illustrate the area-altitude distribution of Urumqi Glacier No. 1. d Altitudedependent profiles of the averaged mass balance and its individual components for the ablation period, including surface melt, subsurface melt, sublimation, solid precipitation, refreezing and deposition cover might result in the large emissivity. Cloud cover as input data in COSIMA was obtained from a parametrization described by Favier et al. (2004). However, this parametrization is not unique (for instant Oerlemans 2001). Additionally, temperature of atmosphere often emits LW in , however, in this study, whether the 2 m air temperature measured at AWS1over the glacier surface can represent the temperature of atmosphere is still to be proved.
The turbulent fluxes are often affected by local meteorological conditions. Q lat peaked in the months prior to 8 June (Fig. 6), indicating sublimation of glacier surface (Fig. 8). But the air temperature remained rather low, which favored a large surface-air vapor pressure gradient, and the lower relative humidity and higher wind speeds also drive turbulence (Fig. 3a-c). Generally, monthly means of Q sens and Q lat were of opposite sign, but absolute values of Q sens were larger than Q lat when air temperature rose, especially after 8 June, increasing the importance of Q sens for surface melt (Table 3).

Geodetic v.s. modeled melt rates
In this study, glacier mass balance of Urumqi Glacier No. 1 was modeled using the AWS1-driven COSIMA during the ablation period in 2018. The mean surface velocity in the same investigative period was 0.026 m day −1 corresponds to 3.3 m year −1 , which was derived by the comparison of two high-resolution UAV photogrammetries . Assuming no significant speed up during the ablation period and considering the 30 m spatial resolution of COSIMA, the dynamical change can be neglected for modeling at seasonal timescales due to derived small surface velocities. Figure 12a illustrated the spatial differences of surface elevation changes between COSIMA and UAV. Modeled results using COSIMA within this study and the geodetic results referred to Wang et al. (2021) based on repeated high-resolution UAV photogrammetries, respectively. Since different densities for snow, firn, and ice were used within COSIMA, we employed an average ice density of 900 ± 17 kg m -3 for the conversion from mass changes to surface elevation changes in COSIMA. The in-situ measured densities of firnsnow data (change in ice thickness) was used to estimate the single-point density conversion of 752 ± 34 kg m -3 during the ablation period of 2018. Differences between both datasets at the tongues are small, while high differences occur at the middle part. The latter is caused by a constant ice flow into this branch over time. Surface velocity is larger in the middle-lower part than in the upper stream . By comparison, the geodetic results of Wang et al. (2021) indicated a stronger thickening in some upper parts,  Table 2) which might be compensated by stronger thinning in the lowest regions compared to this study. The accuracy was within decimeter accuracy, with a mean value of 0.14 m for the ablation period (Fig. 12b). Such decimeter-scale uncertainty supports the acquisition of the glacier elevation changes derived from COSIMA. Figure 12c shows that profiles of differences between surface elevation changes derived using COSIMA and repeated UAV surveys. Overall, both agreed well with each other, but the difference still existed (R = 0.56; Standard error = 0.54). Repeated UAV surveys observed glacier thinning even in the upper-elevation areas, while COSIMA estimated a gain of mass in the upperelevation areas and a loss of mass in the ablation area. The modeled melt rates in this study ranged from 0.2 to 1.7 cm w.e. day −1 in the ablation period of 2018 with the mean value of 0.6 cm w.e. day −1 . According to Xu et al. (2019), the melt rates of 0.9 and 0.8 cm w.e. day −1 were derived from long-range terrestrial laser scanning measurements in the ablation periods of 2015 and 2016 on Urumqi Glacier No. 1, respectively, indicating the slightly reduced mass loss in recent year together with our modeled results. This phenomenon was also founded in Li et al. (2021).

Sensitivity of mass balance to air temperature and precipitation
To assess the sensitivity of the mass balance of the Urumqi Glacier No.1 to climatic factors, various air temperature or precipitation changes as input data were applied to run COSIMA over the ablation period of 2018. Eight independent temperature change scenarios were designed with air temperatures adjusted in 0.5 K steps from − 2 to 2 K while other variables and model parameters were held constant. Eight independent precipitation change scenarios were also established in the same way by perturbing the precipitation in 10% steps from -40 to 40%. The COSIMA was run under the background of these sixteen scenarios as a sensitivity analysis and the results are presented in Fig. 13. The sensitivity of mass balance to increasing air temperature was higher than that to increasing precipitation on Urumqi Glacier No.1, and the dependence of mass balance on changes in air temperature and precipitation was close to linear. However, this does not mean that air temperature is more important than precipitation for controlling changes in mass balance. It only shows that the mass balance will change accordingly when the air temperature changes by 1 K or the precipitation changes by 10%. To roughly keep the mass balance on the Urumqi Glacier No. 1, 1 K increase in air temperature would have to be compensated by at least 40% precipitation change in our study. For Urumqi Glacier No. 1, Che et al. (2019) found that 23% increase of the precipitation could justly compensate to the mass loss caused by 1 K increasing in air temperature (Che et al. 2019). By comparison, the sensitivity result in this study with higher accuracy is larger than Che et al. (2019). There are two probable reasons. Firstly, the forcing data is from AWS1 on the glacier surface, which can accurately show glacio-meteorological conditions. Secondly, COSIMA coupled surface and subsurface mass balance process together, which can account for meltwater percolation, retention, and refreezing. For the continental Haxilegen Glacier No. 51 located close to Urumqi Glacier No. 1, the mass balance was more sensitive to 1 K air temperature change than to a 65% precipitation change . The mass balance of the Guliya ice cap was more sensitive to the changes in moisture related variables than that in temperature . When the air temperature of Qiyi Glacier increased by 1 K, the ELA increased by 172 m, while the precipitation increased by 10%, the ELA decreased by 62 m ). However, the mass loss after a 1 K change temperature in Muji glacier needs to be compensated for by increasing precipitation by approximately 39% (Zhu et al. 2020). To roughly maintain the mass balance of the Shiyi Glacier, a 1 K increase in air temperature must be compensated by at least 35% of the precipitation changes . As one of the maritime glaciers, the mass balance of Parlung Glacier No. 94 was approximately 2 ~ 3 times more sensitive to 1 K temperature change than to 30% precipitation change ). The sensitivities of glacier mass balance in response to climate change were different in various mountain ranges, which were mainly resulted from the discrepancies in the ratio of snowfall to precipitation during ablation period, the amount of melt energy during ablation period, and precipitation seasonality in the different local regions. Although there are significant differences in the sensitivity of different types of glaciers to air temperature and precipitation, extreme continental glaciers have a lower percentage increase than precipitation required for maritime glaciers in order to balance the effects of a 1 K temperature increase (Fig. 13).

Climatic factors controlling mass changes
Summer mass balance on Urumqi Glacier No. 1 plays a dominant role for the annual mass balance Wang et al. 2016). So the modeled summer mass balance was used to assess how climatic factors controlling the mass changes. However, it should be noted, this result is inevitably subject to uncertainties when compared with other studies, which referred to annual mass balance change. Based on the method used in this study by analyzing interannual Fig. 13 Mass balance sensitivity by perturbing the mean air temperature at a step of 0.5 K and the precipitation at percentage intervals of 10% variability of ablation period air temperature and annual precipitation, mass balance on Urumqi Glacier No. 1 is considered to be more strongly controlled by ablation period air temperature than by annual precipitation, because the sensitivity of air temperature on mass balance change (− 0.149 m w.e.) is larger than annual precipitation on mass balance (0.910 m w.e.) (Fig. 14). Similar results found that the mass loss from increasing in air temperature was significantly higher than that from compensating in precipitation in Urumqi Glacier No. 1 based on air temperature and annual precipitation during the period of 1958-2015 (Che et al. 2019). Therefore, it is deduced that the glacier mass loss in Urumqi Glacier No.1 was mostly resulted from increasing in air temperature.
We present the past studies about the control of air temperature or precipitation on mass balance in the western China together with our results in Fig. 14. The difference between mass balance changes at a single glacier strongly underlines the controlling of climate on mass balance, and presents the response of glaciers to climate change. The Urumqi Glacier No. 1 and Haxilegen Glacier No. 51 located in the eastern Tien Shan are continental glaciers. The mass loss of the Urumqi Glacier No.1 was similar to that of the Haxilegen Glacier No. 51, and a statistically significant relationship in mass balance has also been identified . Mass loss of them could be attributed to air temperature rise during the ablation period. For Urumqi Glacier No.1, the former combined with ice temperature increase and albedo reduction on the glacier surface must be considered together. The physical mechanisms could also be suitable for the Haxilegen Glacier No. 51, because both of the glaciers are situated in the eastern Tien Shan with a relatively dry continental climate. More negative mass balances occurred in Zhadang Glacier as a continental glacier, located in the Nam Co basin, south Tibetan Plateau, air temperature can contribute to this mass loss. In summer, positive incoming shortwave radiation as the most significant heat flux, together with sensible heat flux, contributed to positive heat flux for melting for glacier surface melting. As a result, strong summertime melting occurred on the glacier surface (Zhang et al. 2013). The Muztag Ata Glacier No. 15 and Muji Glacier are extreme continental glaciers. The smallest mass loss is observed at the Muztag Ata Glacier No. 15 and Muji Glacier in the eastern and northeastern Pamir regions due to the strengthening westerlies with controlling of precipitation. The mass balance change of the Palung Glacier No. 94, which is a maritime glacier, is also controlled by air temperature. In conclusion, the intensity of air temperature and precipitation controlling mass balance is various in different regions. The mass loss of these glaciers is mainly controlled by air temperature except the Muztag Ata No. 15 Glacier and Muji Glacier, the mass balance of which are mainly dominated by the annual precipitation. Figure 15 shows the comparison of our results with the glaciological mass balance, Degree-day model and Geodetic method results on Urumqi Glacier No. 1. The mass balance obtained by different methods compares with glaciological mass balance to clearly present the optimal model performance. The energy balance model is usually regarded as reference model in calculating mass balance. Che et al. (2019) conducted an energy balance modeling experiment forcing by AWS2 datasets and the result was in line with actual observation. The relative coefficient between the modeled and measured cumulative mass balance was 0.86, and the coefficient of determination was 0.75 (Che et al. 2019). The Degree-day model was more suitable for long-term mass balance estimates, because overall mass balance estimates were in a good agreement with glaciological mass balance for the long term (Wu et al. 2011). In term of the enhanced Degree-day model, the spatial distribution of mass balance in Urumqi Glacier No. 1 showed that the performance was less performed compared with the glaciological mass balance (Huintjes et al. 2010).

Fig. 14
The spatial pattern of mass balance change of Urumqi Glacier No. 1 with other glaciers in western China when ablation period air temperature or annual precipitation is changed by the interannual variability in ablation period air temperature (σ T ) or annual precipitation (C v ) at Chinese national meteorological stations near the glaciers. The mass balance changes are derived using data from previous studies Zhu et al. 2017Zhu et al. , 2018Zhu et al. , 2020Zhang et al. 2018). The red circles with olive green portion show mass balance changes corresponding with different interannual variabilities. Glaciers are shown in light blue, and the background is the Global Multi-resolution Terrain Elevation Data (GMTED2010) dataset sourced from the USGS (http:// topot ools. cr. usgs. gov/ GMTED_ viewer/ viewer. htm). The standard map of the standard map is from the service website of the Ministry of Natural Resources of China (GS (2019)1822) There was actual phenomenon that annual mass balance was mainly related to summer mass balance (mass balance in the ablation period) in Urumqi Glacier No.1 Wang et al. 2016). As shown in Fig. 15, our result was more consistent with annual glaciological mass balance compared with Che et al. (2019). The hourly meteorological records and COSIMA combined with surface and subsurface processes together produce this optimal result. It is therefore regarded as more accurately reflect energy and mass balance process on glacier surface. Geodetic mass balance in ablation period of 2015 and 2016 was more negative compared with annual glaciological mass balance, while our estimate with the value of − 0.77 m w.e. was also more consistent with the annual glaciological mass balance. The performance of the simplified energy balance model is better than that of the Degree-day model in the short time, but the Degree-day model performed better than the simplified energy balance model in the same zone (e.g. the zone around the ELA) (Li 2020). We collect the model comparison studies during ablation period and assess their difference in space and time. Some studies reveal the enhanced Degree-day model offering significant improvements over the classical Degree-day model at the point scale, nevertheless the improvement was limited (Pellicciotti et al. 2005). However, at the glacier scale, the result is less clear. Pellicciotti et al. (2013) showed that there were obvious differences in performance between the enhanced Degree-day model and an energy balance model. MacDougall et al. (2011) also applied an energy balance model and four empirical models, and similar conclusions were obtained. The two models can be demonstrated to be clearly superior to others, and their performance strongly depends on input data and temporal and spatial resolution of the application. For the enhanced Degree-day model, the input meteorological variables need to be extrapolated from point observations to the grid cells of the glacier, as the energy mass balance require a number of input meteorological variables. Some meteorological data, such as wind speed or shortwave radiation, are difficult to model on the glacier surface. On the other hand, interpolation methods also cannot accurately obtain these meteorological data changes on the glacier surface, especially the wind speed, because there is no clear elevation or spatial dependency. It is not clear which is superior at whole glacier and larger scales between the two models. At Urumqi Glacier No. 1, the both performances have not compared with glaciological mass balance as yet. The results presented here are important, since some studies have shown that the modeled mass balance affects runoff projections (Gabbi et al. 2014).
Compared with previous studies, we have clarified the connection between glacier change and the atmospheric conditions. Through the forcing of hourly-scale value for each meteorological variable, this paper has revealed the melting process and mechanism of the Urumqi Glacier No. 1 during the ablation period. The spatial and temporal changes of the energy fluxes and mass balance components are well demonstrated, and the meteorological factors controlling the mass balance of the glacier can also be quantitatively determined. Our study is the first attempt to evaluate the performance of mass balance at Urumqi Glacier No. 1 by linking COSIMA with the in-situ measured meteorological records and to understand glacier energy and mass balance process. Our estimate result is consistent with glaciological mass balance (Fig. 4), and similar with annual glaciological mass balance (Fig. 15). However, our study is limited in time scale and the insight into for model performance, such as parameter instability. In future, we plan to extend the input data time series using the ERA-5 reanalysis data, particularly with regard to glacier projections.

Conclusions
Based on the meteorological measurements in AWS1 on the surface of Urumqi Glacier No. 1 in the Chinese Tien Shan during the ablation period 2018, the temporal variations in energy and mass balance and the response of mass balance to climate variations, were investigated using COSIMA. The optimization of our modeled results is not suitable for traditional optimization method (e.g. Monte-Carlo approach) because of the short-term investigated period. Therefore, we attempted to find the best parameter selection for AWS1. The manual model optimization conducted to meet the climatic observations from AWS1, afterwards, the model evaluation of COSIMA was conducted using surface temperature, snow height and glacier-wide mass balance measurements.
Relative humidity was small while wind speed was high prior to 8 June, but air temperature remained rather low and surface temperature rose. These meteorological conditions will drive turbulence. The main components in energy budget were radiation components (SW net and LW net ) and turbulent fluxes during ablation period. Albedo as a key parameter could affect SW net , and most of the measured increases in albedo can be associated with snowfall events. Surface melt and solid precipitation were two dominant components of mass balance. Most surface melt and solid precipitation occurred simultaneously in June-August with the strongest mass losses in July. Urumqi Glacier No. 1 experienced a significant mass loss, with an average value of − 0.77 m w.e. corresponds to 0.6 cm w.e. day −1 over ablation period 2018. Comparing with the UAV results at Urumqi Glacier No. 1 during the same period showed that the modeled surface elevation change agreed well with surface elevation change determined by repeated UAV surveys. By comparing this result with the geodetic mass balance observations of Urumqi Glacier No. 1 for ablation periods of 2015 and 2016, it is indicated the slightly reduced mass loss was experienced in recent year.
Sensitivity analysis showed that the sensitivity of temperature is greater than precipitation, and mass loss caused by temperature increase of 1 K needed to be compensated by at least 40% precipitation. The air temperature during ablation period was more important than annual precipitation in controlling mass balance changes on Urumqi Glacier No. 1. The key of performance of the energy balance model depended on the accurate of meteorological records and optimization of parameters. In general, the COSIMA can reproduce the glaciological mass balance compared with Geodetic method and other energy balance model. The Degree-day model was more suitable for long-term mass balance estimates, and the performance of the simplified energy balance model was better than that of the Degree-day model in the short time. However, at the glacier-wide scale, it is difficult to identify which is better between the enhanced Degree-day model and energy balance models, that is because it strongly depends on input data.

Incoming shortwave radiation (SW in )
The radiation model within COSIMA considers clear sky and scattered shortwave radiation, a matrix for incoming shortwave radiation generated by radiation model includes effects from both terrain shading and cloud cover (Huintjes et al. 2015a, b).The preparation of distributed incoming shortwave radiation is divided into two steps. First, the potential shortwave radiation value (SW in,x,y,pot ) was calculated according to the radiation model within the COSIMA. The radiation model mainly considers the altitude, latitude, and the spatial resolution of DEM; secondly, according to the ratio of shortwave radiation (SWin, x_AWS, y_AWS, pot) observed by the meteorological station to potential shortwave radiation (rix, y), based on the above ratio, by observing the time change of shortwave radiation during the ablation period in 2018 the incoming shortwave radiation (Qcorrx, y) (Kumar et al. 1997) on the surface is detailed in the following formula.

Albedo
The parameterization of surface albedo (α) follows the scheme of Oerlemans and Knap (1998) where α is determined as a function of snowfall frequency and snow depth α firn , α frsnow , and αice is the firn albedo, fresh snow albedo and ice albedo, respectively. t snow is the time since the last snowfall, t * is a constant for the effect of ageing on snow albedo, h is the snow depth and d * is a constant for the effect of snow depth on albedo.

Incoming and outcoming longwave radiation (LW in and LW out )
LW in , LW out is the incoming and outcoming longwave radiation, respectively. T s is the surface temperature, ε cs is the atmospheric emissivity. N is the cloud cover factor. ε cs is cloud emissivity. E is water vapour pressure (hPa). T air is air temperature at 2 m above the surface (K). RH air is relative humidity. E is the saturation vapour pressure (hPa). T air(C) is air temperature (°C). For a, b and ε cl we take the values of 2, 0.433 and 0.984, respectively, as determined by Klok and Oerlemans (2002).

Turbulent heat exchange
Turbulent heat fluxes Qsens and Qlat are calculated through the bulk aerodynamic method after Oerlemans (2001) between the surface and 2 m, using Tair, relative humidity (RH) and wind speed (u) data. ρ air is the air density (kg m -3 ), c p is specific heat capacity of air (1 004.67 J kg -1 ℃ -1 ), K is the van Karman constant (0.41), h z is the instrument height (2 m), z o is the surface roughness length that changes depending on time from fresh snow (0.24 mm), firn snow (4 mm) and ice (1.7 mm) (A13) air = (p ⋅ 100) 287.058 ⋅ T air 1 + q air ⋅ 0.608 (A14) qair/s = RHair/s ⋅ 0.622 E/s(p − E/s) 100 (Brock et al 2006;Cullen et al 2007;Gromke et al 2011); T s is air temperature, surface temperature, respectively. L E is latent heat of evaporation (2.514×10 6 J kg −1 ); L s is latent heat of sublimation (2.849×10 6 J kg −1 ), q air , q s are specific humidity at 2 m and at the surface (kg kg −1 ). RH air/s is the relative humidity at 2 m or surface. RHair is assumed to be 100% at the surface. p is air pressure (hPa). E and E s are saturation water vapour pressure at 2 m and at the surface.

Ground heat flux
The ground heat flux (Q G ) is the sum of the conductive heat flux (Q C ) and the energy flux from penetrating shortwave radiation (Q PS ). Q PS is calculated following Bintanja and van den Broeke (1995) the extinction of net shortwave radiation (SW net ) in the snow or ice layers is parameterized. Si is the remaining fraction of shortwave radiation reaching down to depth zs. In the top model layer, a fraction ζ is absorbed and an exponentially decreasing flux with constant extinction coefficient β reaches the layers at depth -zs and increases the subsurface temperature. Thus, Q PS is equal to SW net (1 − ζ). For ζ and β we take the values of 0.8 and 2.5 for ice, and 0.9 and 17.1 for snow, respectively, as determined by Bintanja and van den Broeke (1995). Q C is determined from the temperature difference between the surface and the two uppermost subsurface layers of COSIMA and depends on the thermal conductivity (λ) of the medium (ice or snow). λ is calculated from the subsurface density (ρ, in kg m −3 ) after Anderson (1976): Thus, in order to calculate Q G vertical profiles of subsurface temperature (T sub ) and density need to be known. These are calculated within the subsurface model within COSIMA. The model uses a vertical structure that consists of layers with an equal thickness (dz) of 0.2 m. Each subsurface layer is characterized by a temperature, density, liquid water content and depth (z s ). In this layer structure with a total depth z of 10 m, T sub is solved from the thermodynamic heat equation (Anderson 1976).   Data availability Glaciological mass balance data related to this study are submitted to the WGMS and be available at the following website: https:// doi. org/ 10. 5904/ wgms-fog-2019-12 (WGMS, 2019). COSIMA is available online as an open-source software and used for this study (https:// bitbu cket. org/ glaci ermod el/ cosima; Huintjes et al. 2015a). The meteorological data are available upon request by email to the corresponding author (wangpuyu@lzb.ac.cn). Table 6 The sensitivity of input parameters

Declarations
The base value (V) and range (R) were taken from the stated references.