A Test Study of an Energy and Mass Balance Model Application to a Site on Urumqi Glacier No. 1, Chinese Tian Shan

In this study, energy and mass balance is quantified using an energy balance model to represent the glacier melt of Urumqi Glacier No. 1, Chinese Tian Shan. Based on data from an Automatic Weather Station (4025 m a.s.l) and the mass balance field survey data nearby on the East Branch of the glacier, the “COupled Snowpack and Ice surface energy and Mass balance model” (COSIMA) was used to derive energy and mass balance simulations during the ablation season of 2018. Results show that the modeled cumulative mass balance (−0.67 ± 0.03 m w.e.) agrees well with the in-situ measurements (−0.64 ± 0.16 m w.e.) (r2 = 0.96) with the relative difference within 5% during the study period. The correlation coefficient between modeled and observed surface temperatures is 0.88 for daily means. The main source of melt energy at the glacier surface is net shortwave radiation (84%) and sensible heat flux (16%). The energy expenditures are from net longwave radiation (55%), heat flux for snow/ice melting (32%), latent heat flux of sublimation and evaporation (7%), and subsurface heat flux (6%). The sensitivity testing of mass balance shows that mass balance is more sensitive to temperature increase and precipitation decrease than temperature decrease and precipitation increase.


Introduction
There are 7934 glaciers comprising an area of 7179.77 km 2 in the Chinese Tian Shan, accounting for 16% and 14% in number and area of glaciers in China, respectively [1]. A number of studies show that as a result of global warming, glaciers in the Tian Shan have retreated with remarkable mass loss during recent decades, which has an important impact on local water resources, regional ecological environment and industrial and agricultural production [2][3][4][5][6][7].
Driven by solid precipitation and the surface energy budget, glacier mass balance is an immediate indicator of climate variability, and provides linkage between the glaciers and water resources [8,9]). Stake/snow pit observations are the traditional method of in-situ glacier mass balance measurement, but they are time and labor consuming. Comparison of digital elevation models (DEMs) in different periods provides a useful method to estimate the large-scale and long-term glacier mass balance, but cannot be confidently used to compute annual or seasonal surface mass balance of mountain glaciers due to the limitation of temporal and spatial resolution and accuracy [10][11][12][13]. Therefore, energy balance-based modeling of glacier mass balance has drawn great attention in order to investigate the relationship between climate and mass balance. Earlier studies in the Chinese Tian Shan deploy various empirical models, such as linear regression models and Degree-Day models, which can describe the relationship between mass balance and meteorological components [14][15][16][17]. However, they show shortcomings for glaciers with accumulation mainly in summer, and for large areas with complex terrain where radiative processes and non-linear feedback between air temperature, precipitation and albedo play a significant role. Therefore, physical models have been developed to investigate energy and mass balance of glaciers all around the world [18][19][20][21][22][23][24][25][26]. [22] presented the steady-state physically based "COupled Snowpack and Ice surface energy and Mass balance model" (COSIMA). The model solves the energy balance at the surface, driven by atmospheric variables. Energy and mass fluxes are also resolved at the glacier surface and in subsurface layers within the upper several meters of snow, firn or ice. Finally, it computes the glacier mass balance at any useful temporal resolution between minutes and several hours. Previously, it has successfully been applied on Zhadang Glacier, Puruogangri Ice Cap and two glaciers in the eastern Nyainqêntanglha Range, southeastern Tibet in China [22,23,27] as well as on outlet glaciers of the Southern Patagonia Icefield in the Andes [28].
Urumqi Glacier No. 1 (UG1) has a long monitoring history and is a reference glacier in Central Asia. Although many studies have been carried out on this glacier, few focused on glacier energy and mass balance modeling. A parameterized energy balance model of glacier melting was previously used by [29] based on the meteorological observation during the ablation seasons in the 1980s. In their study, statistic correlations between meteorological parameters and energy fluxes were established according to surface melting observations, then daily energy fluxes were calculated and discussed. In the present study, two automatic weather stations (AWSs) were installed to the side and on the east branch of UG1. Detailed meteorological and mass balance observations were conducted during the ablation period in 2018. COSIMA was then applied to estimate energy and mass balance of the east branch of UG1 to investigate the temporal patterns of energy fluxes and mass balance.
Despite the fact that one single ablation season is only a short study period, we are convinced that the findings provide novel insight on the relevant physical processes governing the surface mass balance. Especially, almost the total annual ablation and major parts of the annual accumulation occur in the ablation season, as most of the precipitation within the study region occurs in summer and the continuously very low air temperatures during the winter season. In addition, the AWS cannot be reliably operated during winter due to accessibility. Therefore, for example, the exact verticality of the mast cannot be guaranteed, and riming on the surfaces of sensors corrupts measurements.

Study Area
Located on the northern slope of Tianger Peak, Chinese Tian Shan, UG1 (43 • 06' N, 86 • 49' E) is a valley glacier facing northeast ( Figure 1). This glaciated region is dominated by the westerly jet stream in the upper troposphere, the Siberian anticyclonic circulation and cyclonic disturbances of the west wind circulation [4]. The annual mean temperature and the annual precipitation were approximately −4.6 • C and 460 mm, respectively, from 1959 to 2017 according to the meteorological records of the Daxigou Meteorological Station (3593 m a.s.l.) 3 km away from the glacier terminus. Precipitation mainly originates from moisture flux through westerlies, and 78% of the annual precipitation within the period 1959-2015 occurs from May to August, which predominantly falls as solid precipitation (e.g., snow, graupel, sleet) at higher elevations [30,31].
Due to continuous strong melting, the glacier was separated into two branches (east branch and west branch) in 1993. The total glacier area was 1.56 km 2 in 2015, with an area of 1.06 km 2 for the east branch ranging from 3743 to 4267 m a.s.l. According to [32], the altitude of the equilibrium line (ELA) of UG1 showed a general increasing trend in the period 1959-2008 ascending by 108 m and reaching an altitude of 4168 m a.s.l. in 2008. The mean elevation of ELA during this period was 4056 m a.s.l.
The glaciological mass balance measurement indicated that the ELA was 4180 m a.s.l. in 2018 for the east branch of UG1 marked in Figure 1b. Affected by surface albedo decrease and air temperature increase, the accelerated melting has led UG1 to undergo terminus retreat and area reduction with rates of 4.2 m yr −1 and 0.01 km 2 yr −1 , respectively, over the past 50 years [33,34].

Meteorological Data
Two AWSs were installed, one on the glacier surface (AWS1, Figure 1b) and the second one on the terminus moraine north of the east branch (AWS2, Figure 1b). AWS1 was set up near to the glacier central flowline at an elevation of 4025 m a.s.l in a relatively flat area with a slope of approximately~2 • . It started normal operation on 29 April 2018. The meteorological data observed include air temperature (T), precipitation (P r ), air pressure (P), relative humidity (RH), wind speed (u), incoming and outgoing shortwave radiation (SW in and SW out ), and incoming and outgoing longwave radiation (LW in and LW out ). The glacier surface temperature (T s ) was measured with a SI-111 sensor (Apogee Company). A borehole (8 m) was drilled on 22 August 2018 to measure the ice temperature (T ice ) in a profile at nine depths (see Table 1) with a Campbell Model 109 temperature probe. All these data were stored as 10-min mean values in a Campbell CR1000 data logger. AWS2 was installed in 2011 on the glacier terminus moraine at an elevation of 3835 m a.s.l, providing the same datasets as AWS1 except that AWS2 stores data in 30-min intervals as mean values. Two shield Geonor T200B precipitation gauges were operated close to AWS1 and AWS2, respectively, to measure precipitation in mm w.e. (water equivalent). The detailed sensor information and technical specifications of both AWSs are shown in Table 1. All sensors were checked and leveled every 10-15 days during the field survey in the ablation period. Meteorological data of AWS1 were collected from 29 April to 1 September 2018. During this period, 48 h of data from 12:00 p.m. 20 August to 12:00 p.m. 22 August were missing. The proportion of the data loss period is very small. The gap was filled using data from AWS2 according to their correlation relationship.
The COSIMA model forcing variables include air temperature, air pressure, relative humidity, precipitation, cloud cover, wind speed and downward shortwave radiation from the AWS1. The data correction and gap filling method for various variables are described separately in the following section.

Downward Shortwave Radiation and Cloud Cover
The components of radiation were measured directly by the AWSs. Due to snowfall and resulting snow cover on the top of the radiation sensor, as well as freezing, the radiation measurements might include defective data. Therefore, a method by [36] is used to adjust radiation components effectively. This method considers the problem of radiation sensor riming and suggests that when an abnormal measured value occurs due to riming, the two-hour mean value should be adopted. Afterwards, the correlation coefficient between the downward shortwave radiation of AWS1 and AWS2 reached 0.95, and thus the AWS1 downward shortwave radiation data gaps can be filled with data from AWS2 using a linear transfer function.
Cloud cover (N) is a forcing variable but cannot be observed directly and is estimated using a method described by [37,38]: 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 [39].

Air Temperature and Precipitation
Critical values of temperature are generally used to distinguish the liquid from the solid fraction of precipitation. The precipitation is considered to be liquid when the air temperature is higher than 5 • C, and to be solid when the air temperature is lower than 1 • C. A linear interpolation was adopted in the solid-liquid precipitation transition zone, which is specified by a 1 to 5 • C temperature range [40]. Then, the original precipitation data (P _orig ) observed by T200B were corrected (P _corr ) using a method suggested by [41]: where CE is the capture rate (%) of solid precipitation and u is wind speed (m s −1 ). CE is generally considered to be 90% during calm conditions (u < 1 m s −1 ). Altitudinal gradients between AWS1 and AWS2 were found to be −0.011 • C m −1 for air temperature and 0.016 mm m −1 for precipitation. Therefore, temperature and precipitation data missing in the data record of AWS1 were substituted with data from AWS2 using these lapse rates. The larger lapse rate of temperature is probably because of the difference in the underlying surface between the two AWSs in the study period.

Wind Speed, Humidity and Pressure
Any wind speed data missing in the record of AWS1 were filled with unchanged AWS2 data, although the correlation coefficient of AWS1 and AWS2 wind speed from 24 June to 20 August 2018, was only 0. 41. Hock (2005) pointed out that relative humidity varies little among AWSs in close vicinity. However, the two AWSs in this study are located on different surfaces. Hence, the RH data voids of AWS1 was also filled with AWS2 data in the same period using a linear regression function based on the correlation between the two stations (r 2 = 0.87).
Any data gaps of air pressure are filled using two methods. The first method is to fill the gap using the AWS2 data because the correlation coefficient of air pressure between AWS1 and AWS2 during overlapping periods was 0.82. The other method is to use an isothermal atmospheric pressure equation to fill the air pressure data gaps of AWS1 by applying the correlation between air pressure and altitude of AWS2 [42]: where P AWS1 , P AWS2 , H AWS1 , H AWS2 denote air pressure and altitude of AWS1 and AWS2, respectively, and t is the average air temperature difference between H AWS1 and H AWS2 , and α is a constant (α = 0.036). Since the results from the two methods are consistent, only the second method has been used for gap filling.

Glacier Mass Balance Data
To validate the simulation, glacier mass balance was observed using the stake/snow pit method. Five stakes were deployed near AWS1 by using a steam drill (Figure 1d). The stakes were measured at an interval of approximately every two weeks from 29 April to 1 September 2018. The measured variables include the vertical height of stakes above the glacier surface, superimposed ice thickness and the density and thickness of each snow/firn layer at every individual stake. Then all the variables were converted to mm w.e. Ice density is set to 900 kg m −3 . Snow density was taken from measurements. Superimposed ice density was assumed to be the same as the density of ice since direct measurements were not feasible. The single-point mass balance for an individual stake (b n ) is calculated using where b s , b ice and b si are mass balance of snow, glacier ice and superimposed ice, respectively, in mm w.e. [43,44]. Finally, the averaged mass balance of the five stakes was used to serve as validation for the model.

Glacier Energy and Mass Balance Model
COSIMA consists of a surface energy model and a multi-layer subsurface snow and ice model to compute glacier mass balance accounting for meltwater percolation, retention and refreezing typically at the desired time step. Surface melt serves as a linkage between the surface energy balance and subsurface module. COSIMA is described in detail by [22]. Therefore, only a brief description is provided here.
The surface energy balance within COSIMA can be written as [45] where F is the energy flux. SW in , LW in and LW out represent downward shortwave radiation, incoming longwave radiation and outgoing longwave radiation, respectively. αis the surface albedo. Q sens and Q lat are the turbulent sensible and latent heat flux, respectively, and Q G is the ground heat flux which consists of fluxes of heat conduction (Q c ) and penetrating shortwave radiation (Q ps ). Heat flux from liquid precipitation is neglected. Energy fluxes toward the surface have a positive sign, fluxes from the surface towards the subsurface or atmosphere have a negative sign. The resulting overall energy flux F is equal to surface melt energy (Q melt ), when the surface temperature reaches the melting point (273.15 K). Otherwise Q melt is zero if the surface temperature is below the melting point. The bulk aerodynamic method was adopted after [45] to calculate turbulent heat fluxes (Q sens and Q lat ) between the surface and 2 m above the surface. Turbulent heat fluxes functions can be described as follows: ρ air = (p × 100)/(287.058 × (T air (1 + 0.608 × q air ))) (8) q air/s = (RH air/s × 0.622(E /s /(p − E /s )))/100 (9) where C p is specific heat capacity of air (1004.67 J kg −1 K −1 ), ρ air is air density (kg m −3 ), k is the von Karman constant, and h z and z 0 are the reference height (2 m) and the surface roughness length, respectively. The roughness lengths of fresh snow, firn and ice are 0.24 ± 0.05 mm, 4 ± 2.5 mm and 1.7 ± 1 mm, respectively [46][47][48]. u is wind speed (m s −1 ); T air and T s are the air temperature at 2 m above the surface and at the glacier surface (K), 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 and q s are specific humidity at 2 m above the surface and at the surface (kg kg −1 ); RH is assumed to be 100% at the surface; p is air pressure (hPa); E and E S are saturation water vapor pressure at the instrument height (2 m) and the surface, respectively. A correction of the turbulent fluxes Q sens and Q lat is required for the aerodynamically more stable atmosphere over a melting glacier [49]. The stability is defined by the bulk Richardson number (R i ): where g is the gravity acceleration (9.81 m s −2 ), T air(c) corresponds to air temperature at 2 m above the surface ( • C) and h z is the reference height (2 m), T air is air temperature at 2 m above the surface (K) and u is the wind speed (m s −1 ). Considering of R i , the correction of turbulent fluxes for stable atmospheric conditions is carried out after [50]: Q sens = Q sens and (11) The subsurface model within COSIMA divides, in our case, the 10-m snow and ice pack below the glacier surface into 50 layers with an average thickness of 0.1 m for each individual layer. Due to the effects of the physical condition (e.g., temperature, density and liquid water content) each layer is different. The parameters for the model on UG1 are shown in Table 2. The initial snow depth for this study is set to 0.32 m based on field observations. The initial temperature profile is linearly interpolated between the surface temperature and bottom temperature. The measured temperature at a depth of 8 m was between −8.5 and −5.5 • C from August 2018 to April 2019, so the 10-m temperature was set to −7 • C. The initial density profile is interpolated based on the measured density of snow, firn and ice.   [58] According to [22], the uncertainty of modeled mass balance per day (U day ) by COSIMA can be determined by the following.
where X j,max and X j,min are the maximum and minimum mass balance for the five ablation stakes by running the model with the calibrated precipitation gradient (0.016 ± 0.003 mm m −1 ), and by further averaging the results of the six measurements. The variable n j is the number of days of the measurement interval j. Finally, the uncertainty of the modeled mass balance for the ablation period is determined to be ±0.03 m w.e. by multiplying U day and the days of the ablation period.

Meteorological Observations
Meteorological data from AWS1 for the observation period are shown in Figure 2. The air temperature fluctuated strongly with a mean daily temperature of −3.2 • C in the early ablation season from 29 April to 8 June 2018. Daily mean air temperatures above 0 • C occurred primarily between early June and until the end of August. It then decreases again after 19 August until early September. Relative humidity fluctuated even more than air temperature with a daily mean relative humidity of 69% during the ablation period. Relative humidity is low in May and relatively high with less fluctuation in June and early July.
Daily wind speed over the ablation season ranged from 0.   Table 3 and Figure 3 show the monthly and daily surface energy balance components in the ablation season, respectively. Net shortwave radiation (SWnet) (84%) and sensible heat flux (16%)  Table 3 and Figure 3 show the monthly and daily surface energy balance components in the ablation season, respectively. Net shortwave radiation (SW net ) (84%) and sensible heat flux (16%) contribute the most to glacier melting. The minimum value of monthly SW net was +30.48 W m −2 and occurred in June. The maximum value of +108.43 W m −2 occurred in August. Besides the minor influence of cloud cover, the large range is mainly related to the variation of surface albedo over the study period with typically lower albedo later in the ablation season. In fact, the estimated cloud cover was least in June while still SW net was comparatively low. Nonetheless despite large solar radiation during that period, the surface melting is less since the surface was covered with snow cover resulting in high albedo and little absorbed solar radiation. Later, as a response to both diminishing snow cover and increasing liquid precipitation, the albedo on average decreases until reaching the minimal albedo in late August resulting in increased melting at that time. The sensible heat flux was positive during the ablation period (+10.89 W m −2 ) due to higher air temperature at 2 m than at the glacier surface. This results in a downward heat flux towards the glacier surface. LW net , Q lat , Q melt (melt energy), and Q G accounted for 55%, 7%, 32%, and 6%, respectively as energy expenditures. The monthly mean of net longwave radiation (LW net ) was −39.3 W m −2 during the ablation period. The relative humidity (mean of 69%) was high during June and July. Since LW out was always larger than LW in at monthly average, LW net was negative accordingly. During the observation period, air temperatures were higher than the surface temperature and the air vapor pressure was unsaturated on average. Consequently, the sublimation and evaporation were more than freezing and condensation so that the latent heat flux was mostly negative. The average value of Q G is −4.85 W m −2 . Ground heat flux is probably biased by assumptions of constant ice temperatures at the bottom of the surface layer (i.e., at 10-m depth) in the model.

Mean Diurnal Cycle of Meteorological Variables and Surface Energy Balance Components
The mean diurnal cycles of the meteorological variables and surface energy balance components in the ablation period are shown in Figure 4. Air temperature is generally higher than surface temperature throughout the ablation period. The difference between them during the day is less than during the night due to strong radiative cooling of the surface during the night. Wind speed shows small fluctuations between day and night for the ablation period. The maximum values of relative humidity and wind speed are possibly due to enhanced catabatic flow down the glacier occurring during the night. Further, the air temperature is lower at night obviously causing the relative humidity to increase if absolute humidity does not change. The mean diurnal cycle of sensible heat flux is similar to the cycle of the latent heat flux. However, the sensible heat flux shows significantly higher values during the night than during the day. The mean latent heat flux was always negative throughout the diurnal cycles but shows significant differences between day and night. The minimum value of latent heat flux was −8.57 W m −2 and appeared at 8 a.m. because of larger relative humidity. R i which provides a measure of stability is used to adjust turbulent fluxes under stable conditions. Its value is always below 0.2 indicating a neutral-to-slightly-stable boundary layer that only partly inhibits turbulent exchange. Along with larger wind speed during the night this causes the latent heat flux to become positive between 4 p.m. and 5 p.m. triggering some condensation at the surface.
Penetrating shortwave radiation as part of incoming shortwave radiation decays exponentially with depth. It shows similar changes to net shortwave radiation and incoming shortwave radiation, but amounts to much smaller absolute values. The flux of heat conduction was positive indicating that the upper ice layer obtained energy in the night while it released energy during the day when the ice temperature was equal to melting point [59]. In general, the ground heat flux consists of heat conduction and penetrating shortwave radiation. Ground heat flux was positive and directed towards the surface to compensate for the release of energy due to radiative surface cooling in the night when the subsurface snow/ice is at 0 • C. However, because ground heat flux mainly is used for the release of cold storage in the day, the mean value was negative. As expected, net radiation was positive during the day but negative in the night. Maximum surface melt occurs later in the day than the maximum of net radiation. This indicates that the interplay of the various energy fluxes shifts during the day from sublimation/evaporation and surface as well as air warming early in the day to surface melting later in the day.

Mass Balance at the Site of the AWS1
The mass balance at the glacier surface is the sum of accumulation and ablation, determined by snowfall, refreezing, condensation, and deposition as accumulation terms and surface melt, subsurface melt, sublimation and evaporation as ablation terms. The modeled mass balance change at the site of the AWS1 on UG1 during the study period is shown in Figure 5. The final net mass balance during the ablation period is −0.671 m w.e. Glacier surface melt accounted for the main mass loss (−0.742 m w.e), followed by subsurface melt (−0.114 m w.e). Sublimation and evaporation also contribute a small mass loss (−0.022 m w.e). Mass is primarily gained by snowfall (0.196 m w.e) and refreezing in subsurface layers (0.007 m w.e). This indicates that refreezing is not relevant in the study period and that almost all the surface and subsurface meltwater is removed from the glacier as glacial runoff. Since, on average, water vapor pressure at 2 m is less than at the surface, and wind speed is generally low, a small but almost negligible condensation of 0.005 m w.e is produced. Figure 5. Mass balance components and cumulative mass balance at the site of AWS1 on UG1 during the study period. Glacier mass balance is the sum of solid precipitation, refreezing-surface melt, condensation (positive) and surface and subsurface melt and sublimation (negative).

Uncertainty of Measured Mass Balance
The assessment of measurement uncertainties is crucial during the field survey. Measurements are prone to errors in stake readings, snow/firn density measurements, stakes sinking, snow/firn or superimposed ice misidentification, and repeated measurements [60][61][62][63]. According to the estimation method from [61], the uncertainty of single-point mass balance during the ablation period of 2018 is estimated to be ±0.16 m w.e., taking the possible erroneous judgments of snow/firn/ice layers and repeated measurements into account. Then the relative uncertainty of observed mass balance is 25% over the ablation period.
The standard deviation of the mass balance measurement at the five stakes was also calculated over the seven observation periods (Table 4). Since strong melting occurs in late July and August, measurement intervals are denser after 20 July, while the remaining two intervals are longer. The most intensive melt occurs around 10 August and so the largest standard deviation arises in the period of 1-10 August. The mean standard deviation of five stakes is 0.029 mm w.e. for the whole observation period.

Uncertainty of Modeled Mass Balance
To validate the model results, the observed and modeled values of mass balance during the ablation season were compared and a regression analysis was performed. The mass balance was modeled according to seven observed dates. Figure 6 indicates that modeled mass balances are in agreement with observed values. The difference is within 5%. The correlation coefficient is 0.96 and the root-mean-square error (RMSE) is 0.22 m w.e. From Figure 6 it can be seen that the modeled mass balances are mostly more negative than measured values. This is hard to explain but may be related to some extent to the fact that this version of the COSIMA model uses a fixed albedo value for bare ice which might be a little too high compared to the actual real albedo at the sites of the ablation stakes. The glacier surface temperature is an important variable to judge the model performance, especially regarding daily cycles. The daily observed and modeled surface temperatures are highly correlated (r 2 = 0.88), but with a relevant RMSE of 3.24 ( Figure 7).  Moreover, different precipitation types including rain, snow and sleet have a significant impact on the glacier energy and mass balance. There are several approaches to quantify this. For example, [64] found that the precipitation types are highly dependent on surface elevation and [65] used a typical phase partitioning method to specify the types of precipitation based on near-surface air temperature. This assumption on precipitation type can cause considerable uncertainties. According to previous studies by [40,52,66], the upper threshold of precipitation phase (all liquid above the threshold) and the lower threshold for precipitation phase (all solid below the threshold) was 6.5 ± 0.5 • C and 1 ± 1 • C, respectively.

Sensitivity of Mass Balance to Climatic Factors
A model sensitivity analysis of mass balance to various scenarios of climate change was carried out (Table 5 and Figure 8). The result shows that response of mass balance to air temperature increase and precipitation decrease is more sensitive than to air temperature decrease and precipitation increase indicating a non-linear system response. When precipitation remains unchanged, the mass balance would decrease by 0.91 m w.e if air temperature increases by 1 K which is much higher than the increase in mass balance of 0.60 m w.e if air temperature decreases by 1 K. When temperature is unchanged the mass balance would decrease by 1.07 m w.e. if precipitation decreases by 10%, but only a 0.37 m w.e. mass balance increase occurs as a response to a 10% positive forcing of precipitation. The scenario with 2 K temperature increase confirms this trend. The probable reason for such a non-linear response is that the air temperature change (increase/decrease) will directly influence the share of snow of total precipitation, and precipitation decrease will cause significant reduction of surface albedo. This indicates that the effect of shortwave radiation absorption is very prominent when air temperature increases and precipitation decreases. The glacier surface albedo determines shortwave radiation absorption and is influenced by air temperature as well as precipitation. As the air temperature increases, exposure of bare ice and surface debris, increased meltwater at the surface and accelerated snow melt cause the reduction of surface albedo. Vice versa, the glacier surface albedo is sensitive to snowfall and it will increase rapidly when snowfall events occur more often. Moreover, it can be seen that a 2 K temperature increase would not cause linearly doubling effect on mass loss compared to a 1-K temperature increase effect if precipitation is unchanged. The effect of changing precipitation exhibits similar trends, but not as distinct as temperature changes. This is very interesting but needs to be investigated further, because, if true, it means that along with further climate warming, the rate of glacier mass loss would decrease.

Comparison to Other Glaciers in Northwest China
Although the direct comparison with other continental glaciers is very difficult because of the differences in glacier mass balance, surface energy and mass balance models, and study periods, it can still reveal overall characteristics and variability of energy fluxes and mass balance for glaciers in the Chinese Tian Shan and Qilian Mountains (Table 6) since they are located in similar climate conditions and not very far from each other. Further, it has to be taken into account whether measurements and modeling refer to a specific site in the ablation or accumulation zone on a glacier, or if areal estimates for the whole glacier are considered. The comparison shows that the net radiation (R net ) accounts for 61% to 93% of the total energy intake on these glaciers. On UG1, R net is 61% and 85% in this study and in earlier years, respectively; it is 81% and 92% for the different altitude of Keqikar Glacier in the Chinese Tian Shan, 73% and 93% for the different altitude of Laohugou Glacier No. 12, and 82% for Qiyi Glacier in Qilian Mountains (Table 6). This indicates that R net is significantly different between different glaciers and even between different elevations and periods on a same glacier, meaning the surface characteristics are very important to shortwave radiation absorption.
The sensible heat flux rates of Qiyi Glacier (+14.2 W m −2 ) [19] and Keqikar Glacier (+14.4 W m −2 ) [67] are larger than that of the site of AWS1 on UG1 (+10.89 W m −2 ) and Laohugou Glacier No. 12 (+5.7 W m −2 ) [25]. Nonetheless these numbers indicate a fairly proportional relationship between overall energy balance and sensible heat flux driven by similar ranges of the main driving variables: air temperature and wind speed. Overall negative values of the latent heat flux at all study sites indicate that the latent heat exchange removes energy from the surface. The latent heat exchange of Keqikar Glacier (−23 W m −2 ) was more intense during the ablation period compared to other glaciers, indicating a larger vapor pressure gradient and stronger evaporation. Further studies ought to be carried out at similar locations on glaciers, a similar time in the annual cycle, and with the same modeling approach in order to better understand how different meteorological conditions impact the climatic surface mass balance of glaciers in the Chinese Tian Shan.

Conclusions and Outlook
The coupled snow and ice surface energy and mass balance model COSIMA was used to simulate the glacier surface energy and mass balance at site AWS1 on UG1 during the summer season 2018 using air temperature, air pressure, relative humidity, precipitation, cloud cover, wind speed and downward shortwave as forcing variables. The modeled cumulative mass balance of UG1 was −0.67 ± 0.03 m w.e. and the modeled mass balance was in very good agreement with the observed value of −0.64 ± 0.16 m w.e. The difference is within 5%, with a correlation coefficient of 0.96 and RMSE of 0.22 m w.e. The overall good agreement indicates that COSIMA can be efficiently applied to calculate the energy and mass balance for UG1. The energy components causing glacier ablation are mainly from net shortwave radiation (84%) and then from sensible heat flux (16%). The energy expenditures include the net longwave radiation (55%), the heat flux for snow/ice melting (32%), latent heat flux (7%) and ground heat flux (6%). Affected by the energy budget, the modeled cumulative mass balance is mainly dependent on surface melt and snowfall.
The sensitivity test of mass balance to air temperature and precipitation shows that mass balance is more sensitive to temperature increase and precipitation decrease than temperature decrease and precipitation increase, probably due to significant negative surface albedo feedback to increasing temperature and decreasing precipitation. Unsurprisingly, comparisons with other glaciers in the Chinese Tian Shan and Qilian mountains indicate that there are large differences between different elevations on the same glacier, indicating that surface characteristics that change with altitude have an important influence on surface energy balance.
In order to further study and understand the glacier processes and mechanisms, as well as the precise energy and mass balance characteristics of UG1 and its response to climate change, distributed energy and mass balance modeling should be applied in the future. To do this, more observations need to be carried out at various altitudes, especially further down in the ablation zone where the strongest melting occurs. In addition to the conventional observation data, remote sensing data, such as lidar data and time-lapse photography might be incorporated to assist in constraining surface albedo in space and time.