**3. Results**

*3.1. Air Temperature and Precipitation*

The near-surface air temperature and precipitation are essential parameters in forcing glacier mass balance models. These parameters greatly influence the ablation and expansion of mountain glaciers in the summer, with high air temperatures and liquid precipitation accelerating the ablation process and high solid precipitation accelerating the expansion process. The observed and modeled T2 and daily rainfall, and their diurnal variation over the Parlung No. 4 Glacier are shown in Figure 2. Compared with the T2 observations, both Sens1 and Sens2 experiments applying the realistic land cover yielded similar T2 values and diurnal patterns in the 2016 ablation season, while the CTL experiment applying the default land cover greatly overestimated T2. The mean T2 in the ablation season was observed to be 3.68 ◦C and was estimated to be 5.24 ◦C in the CTL experiment, 3.17 ◦C in the Sens1 experiment, and 3.23 ◦C in the Sens2 experiment. The mean T2 deviation (absolute deviation) reached 1.56 (1.88) ◦C, −0.51 (1.19) ◦C, and −0.45 (1.21) ◦C in the CTL, Sens1, and Sens2 experiments, respectively. The RMSE of T2 between the observations and CTL (Sens1) (Sens2) estimates was 2.42 (1.51) (1.54) ◦C and the CC was up to 0.7. On the whole, the model accurately recreated the diurnal cycle of T2 on the glacier. The T2 estimates were significantly correlated with the ground observations at the specified significance level of 0.01. WRF applying the real updated land surface type (Sens1 and Sens2 experiments) successfully reproduced T2 in the ablation zone of the glacier in summer 2016, confirming the significance of using realistic land surface types in model simulations of near-surface meteorological elements (Figure 2a).

**Figure 2.** Observed and modeled (**a**) 2 m air temperature (T2), (**b**) daily rainfall, and (**c**) accumulated diurnal variation of rainfall in the ablation zone over the Parlung No. 4 Glacier, with a box whisker diagram in the upper right corner.

Precipitation processes of rainfall and snowfall present different mechanisms for glaciers' mass balance. The heat of rainfall will be transferred to the glacier contributing to the melting of glaciers, while snowfall contributes to the expansion of glaciers by providing solid precipitation particles. Therefore, the precipitation differences for the contributions from liquid precipitation (rain) were investigated first, and then solid precipitation (snow). The observed rainfall mainly occurred in July (89.67 mm) and August (47.21 mm), accounting for 47% and 25%, respectively, of the total rainfall in the ablation season over the Parlung No. 4 Glacier. The modeled rainfall also mainly occurred in July and August. However, the model greatly overestimated light rainfall and underestimated moderate rainfall. For the accumulated rainfall in the ablation season, the observed total rainfall amount was 190.6 mm, the Sens2 estimate was 417.2 mm, and the Sens1 and CTL estimates were approximately 433 mm. The three experiments estimated the total rainfall to be more than double the observed total amount (Figure 2b). In terms of diurnal variation of accumulated rainfall, the observed maximum accumulated rainfall in summer 2016 occurred at nighttime (21:00–23:00 LST; LST = UTC + 8 h), reaching 48.33 mm, and rainfall from the afternoon to sunrise (15:00–08:00 LST) accounted for 92% of the daily rainfall. The model demonstrated its capacity to simulate the diurnal variation curves of rainfall but notably overestimated the rainfall amount, with the maximum accumulated rainfall in summer reaching 92.83 mm (00:00–02:00 LST), 111.06 mm (03:00–05:00 LST), and 80.95 mm (00:00–02:00 LST) for CTL, Sens1, and Sens2 modeling, respectively. For the main rainfall duration period (15:00–08:00 LST), the three experiments modeled rainfall that accounted for more than 80% of the total daily rainfall (Figure 2c).

Old snow melted and small amounts of sleet occurred before 20 September, then snowfall ensued thereafter. The Sens1 and Sens2 experiments successfully estimated snowfall in late September, although Sens1 falsely simulated a large amount of snow before 20 September. The observed and modeled three-hourly snow water equivalent (SWE) in late September and the corresponding accumulated diurnal variation in snowfall and snow melt is shown in Figure 3. Compared with the snowfall observations, the Sens2 experiment reproduced snowfall on 20 September followed by rapid snow melt, but failed to reproduce snowfall on 24 September. The two snowfall events were reproduced by Sens1, but the first snowfall amount was remarkably overestimated and this was followed by slow snow melt. The total snowfall amounts from the ground observations, Sens1, and Sens2 estimates were 11.77, 28.58, and 18.47 mm, respectively, with the model greatly overestimating the snowfall amount (Figure 3a). In addition, the observed snowfall mainly occurred in the morning before 11:00 LST, with the largest amount at sunrise, presenting a single peak pattern. However, both Sens1 and Sens2 simulated the snowfall occurrence as double peaks, with the first high peak occurring at 6:00 LST and the second low peak at 17:00 LST for Sens1, and the first peak occurring at 5:00 LST and the second slightly lower peak at 11:00 LST for Sens2. Remarkably, snow melted at local noon in the estimates of both experiments, and the peak of snow melt from Sens2 lagged behind Sens1 by approximately 2 h (Figure 3b).

**Figure 3.** (**a**) Observed and modeled three-hourly SWE in late September over the Parlung No. 4 Glacier (solid line indicates snow accumulation and the dashed line indicates snow melt). (**b**) Corresponding accumulated diurnal variation of snowfall and snow melt estimates.

### *3.2. Net Radiation Flux and Albedo*

Net radiation is a direct and immediate forcing of glacier mass balance, which consists of net shortwave radiation and net longwave radiation; excess energy will melt the snow and ice and affect the near-surface air temperature. A comparison of the net radiation and its components, and the diurnal variation of net radiation between observations and model estimates, is displayed in Figure 4. It is clear to see that similar performance of WRF simulating net longwave radiation was achieved for Sens1 and Sens2 experiments with the similarly temporal pattern and the same range of values (Figure 4a,b). The mean observation of the net longwave radiation was only 13.77 W m<sup>−</sup><sup>2</sup> and the mean estimates reached 30.56 W m<sup>−</sup><sup>2</sup> for both experiments (Sens1, Sens2). Compared with the observed net longwave radiation, the net longwave radiation was rather consistently overestimated by the two experiments (Sens1, Sens2), with an average overestimation value of 17 W m<sup>−</sup>2. A similar performance of WRF in relation to the net shortwave radiation and the net radiation was achieved, revealing the variation in the net shortwave radiation to be a fundamental factor in the variation in net radiation. Owing to the distinct glacial albedo scheme that was implemented in WRF, the instantaneous net shortwave radiation and net radiation was distributed discretely on both sides of the perfect fitting line (y = x) for Sens2, and concentrated below the perfect fitting line (y = x) for Sens1, with rather large underestimations (Figure 4c–e). The average net shortwave radiation for observation, Sens1, and Sens2 estimates was 137.21, 84.2, and 172.41 W m<sup>−</sup>2, respectively. The average net radiation was calculated to be 123.44 W m<sup>−</sup><sup>2</sup> using in situ observations, 53.64 W m<sup>−</sup><sup>2</sup> for Sens1 estimates, and 141.85 W m<sup>−</sup><sup>2</sup> for Sens2 estimates. The linear regression lines for Sens1 and Sens2 estimates were y = 0.47x + 19.61 (CC = 0.81) and y = 0.92x + 46.16 (CC = 0.77), respectively, for the net shortwave radiation; and y = 0.45x − 1.95 (CC = 0.83) and y = 0.88x + 30.92 (CC = 0.78), respectively, for the net radiation. These illustrate that Sens2 apparently outperforms Sens1 in estimating the instantaneous net radiation and its dominant component (net shortwave radiation), and shows a consistent conclusion in terms of the daily mean net radiation estimates that Sens2 apparently outperforms Sens1. Owing to the notable discreteness of net shortwave radiation and the net radiation from Sens2 estimates, Sens2 yielded large RMSEs of ~140 W m<sup>−</sup><sup>2</sup> for both the net radiation and net shortwave radiation. Such large RMSEs were also simulated by the Sens1 experiment because of remarkable underestimations of the net shortwave radiation and net radiation. Compared with the observed diurnal variations in the net radiation, the two experiments achieved identical diurnal cycles, with maximum values occurring at local noon. Yet, different peaks were calculated in Sens1 (228.92 W m<sup>−</sup>2) and Sens2 (503.99 W m<sup>−</sup>2), the latter being closer to the observed value of 454.62 W m<sup>−</sup><sup>2</sup> (Figure 4f). Coincident high CCs of ~0.8 were calculated between the observations and experimental (Sens1, Sens2) estimates of net radiation and net shortwave radiation, which passed the significance test at the specified 0.01 level.

Surface albedo is a key factor in net radiation and is the main determinant of the energy and mass balance of a glacier. Rapid changes in the albedo arise because of snowfall and the subsequent evolution of the snowpack, with large values associated with fresh snow, and decreasing values as snow melt progresses. Our experiments applied the snow albedo scheme and a bare ground albedo scheme for Sens1 and Sens2, and additionally involved a modified two-stream radiation transfer scheme considering the three-dimensional structure of the canopy in the control experiment. The albedo was calculated from a combination of the reflected radiation of snow, bare ground, and vegetation canopy in the model, and was compared with ground observations and MODIS product (Figure 5). The observed daily albedo declined sharply from 0.75 to 0.24 before 29 June, slightly increased to 0.47 on 3 July, then remained at a rather low, variable value of approximately 0.23 until 19 September. Subsequently, the observed albedo rapidly increased to 0.75 on 20 September, then to 0.78 on 24 September (Figure 5a). The variation in the observed daily albedo indicates that the seasonal snowpack in the ablation zone was melting and did not disappear until late June, when the ice became exposed on the surface. Subsequently, the observed variable

low albedo value before 19 September not only results from small amounts of sleet and subsequent snow melt increasing the albedo, but also from rainfall increasing the liquid water content above the glacier's surface and in the meantime decreasing the albedo. In the meanwhile, the observed bare ice albedo fluctuates around 0.23 which was affected by the liquid water content. Therefore, the low ice albedo of 0.23 instead of the minimum albedo is set as the background albedo of the glacier in the Sens2 experiment. The large albedo value in late September resulted from two apparent snowfall events. Owing to the very different treatments of snow, canopy, and bare-ground albedo in the WRF model, the three experiments yielded substantially different surface albedo estimates. The albedo maintained a value of 0.12 in the ablation season in the CTL experiment because the surface type was wrongly classified as open shrub-land. The Sens1 experiment applied the default albedo scheme on the glacier, presenting a very high albedo value of ~0.8 on 3 July and 9, 20, and 21 September, while retaining a constantly high albedo value of 0.68 for the rest of the ablation season. The Sens2 experiment, using the modified albedo scheme on the glacier, yielded a constantly low albedo value of 0.35, which was close to the in situ observation, before increasing to a maximum of 0.78 during snowfall in late September (Figure 5a). Apparently, the albedo error of Sens2 in the ablation season is significantly smaller than that of Sens1, while the error of Sens1 before the middle of June and after the middle of September is significantly lower than that of Sens1. This reveals that the modified glacial albedo scheme for Sens2 is only applicable to the ablation season.

**Figure 4.** Comparison of the net longwave radiation (**a**) daily mean NetLW; (**b**) instantaneous NetLW, net shortwave radiation; (**c**) instantaneous NetSW; and net radiation (**d**) instantaneous Rn; (**e**) daily mean Rn), and the diurnal variation in instantaneous net radiation (**f**) between the observations (OBS: grey) and model estimates (Sens1: green; Sens2: blue). The grey dashed line denotes the perfect linear fit line of y = x, the green solid line denotes the linear fit between the observation and Sens1 estimates, and the blue solid line denotes the linear fit between the observation and Sens2 estimates in (**c**) and (**d**).

**Figure 5.** Comparison of model-estimated albedo at 14:00 LST against observations (**a**) and MODIS retrieved product on 21 August 2016 (**b**). The solid black circle indicates the location of the AWS.

In addition, the performance of WRF in relation to the albedo estimates when applying different glacial albedo schemes was also evaluated against the MODIS-retrieved high pixel resolution and the daily albedo product over the whole glacier. MODIS retrieval of albedo on 21 August 2016 was selected to assess the model performance on simulating the spatial distribution of albedo (Figure 5b), because the albedo in the pixels of the glacier was measured in the highest quality on 21 August 2016. It revealed that the glacial albedo was different across the pixels of the glacier in the summer with a low value of 0.2–0.4 in the ablation zone and a high value of 0.5–0.7 in the accumulated zone. The model applying the default glacial albedo scheme greatly overestimated the albedo of the whole glacier particularly in the accumulated zone where the albedo kept to 0.7–0.8. The model applying the modified glacial albedo scheme showed the remarkable advantages in relation to the albedo estimates in the ablation zone, but underestimated the albedo in the accumulated zone. Therefore, in the ablation zone, using the observed background albedo had a large effect on the simulated albedo (Figure 5b). Besides, we used the MODIS-retrieved albedo in pixels of the glacier, where the ideal quality of the spectral bands, no cloud effect, and low aerosol quantity were recognized in the ablation season. The albedo RMSE and the mean error between Sens1 and the MODIS product were 0.22 and 0.18, respectively, indicating an overestimation of the albedo of the glacier when the default glacial albedo scheme was applied in WRF. In contrast, the albedo RMSE and the mean error between Sens2 and the MODIS product were 0.16 and −0.08, respectively, illustrating a slight underestimation of the albedo of the glacier when the modified glacial albedo scheme was applied in WRF. The albedo estimates from WRF were significantly correlated with the MODIS retrievals at the specified significance level of 0.01, with similar CC values of ~0.5.

On the whole, the albedo was largely underestimated with a rather low value of 0.12 kept in the CTL experiment due to wrong land cover that was used in the model. The Sens1 and Sens2 experiments using the updated land cover successfully simulated snowfall events in late September. However, Sens1 gave a rather high albedo, with the

average value of 0.68 leading to an underestimation of the net shortwave radiation and net radiation. Sens2 gave a low albedo, with the average value of 0.35 being closer to the in situ observations when little snow appeared, causing an overestimation of the net shortwave radiation and net radiation in early June, followed by underestimation in July and August (Figures 4 and 5a). On average, surface shortwave radiative forcing reached to −141 ( −60) W m<sup>−</sup><sup>2</sup> due to albedo changes that were induced by land cover in CTL updated to the one in Sens1 (Sens2), which shows the cooling effect. Surface shortwave radiative forcing reached to 85 W m<sup>−</sup><sup>2</sup> due to glacial albedo changes that were induced by albedo the scheme in Sens1 that was updated to the scheme in Sens2, which shows an apparent heating effect and may cause a considerable ablation of the glacier.

### *3.3. Turbulent Heat Flux*

The surface net radiation determines the energy redistribution between turbulent water vapor and heat fluxes. Albedo is a key factor of net radiation. In our previous analysis, albedo schemes present large impacts on shortwave radiation estimates, but similar T2 is simulated by both Sens1 and Sens2 experiments. Maybe compensating the effects of energy fluxes to/away from the ice surface results in the similar T2. In order to explain this finding reasonably, the modeled daily mean sensible and latent heat fluxes above the Parlung No. 4 Glacier, the ground heat flux beneath the glacier, as well as the daily mean surface heat budget in the ablation season are shown in Figure 6. It shows that the Sens1 and Sens2 experiments simulated similarly low values and consistent daily variations of sensible ( −60 to 0 W m<sup>−</sup>2) and latent ( −51 to 36 W m<sup>−</sup>2) heat fluxes over the entire time. The excess surface energy enters the glacier in the form of ground heat flux. High ground heat flux from Sens1 (15–184 W m<sup>−</sup>2) and Sens2 (30–295 W m<sup>−</sup>2) heats the glacier, resulting in surface and subsurface ice melt (Figure 6a). From the boxplots analysis of the daily mean surface heat budget (Figure 6b), the large difference of the net radiation between Sens1 and Sens2 estimates attributes to the large difference of net shortwave radiation that is caused by different glacial albedo schemes. The simulated net radiation, sensible and a proportion of the latent heat flux by both Sens1 and Sens2 experiments are used to heat the surface glacier first and then transmit downwards to heat the subsurface glacier. The main energy contributor to surface and subsurface ice melt from both experiments is the net radiation, supplemented by sensible heating. Similar sensible heat flux was reproduced by both Sens1 and Sens2 experiments, which is highly linked to the temperature difference between the ground and air, and the ground temperature of the glacier kept to 0 ◦C. These account for the similar T2 that was estimated by WRF applying different glacial albedo schemes. The average ground heat flux from Sens1 was 88 W m<sup>−</sup>2, which was much lower than that from Sens2 (182 W m<sup>−</sup>2). This is because the albedo from Sens1 was much higher that from Sens2, causing apparently lower net radiation and reduced energy entering the glacier.

In addition, ice melt leads to an increase in the liquid water content in the glacier. The liquid water content estimates in the different layers in the ablation region of the glacier are displayed in Figure 7. This shows that excess energy is first used to melt the surface ice (the upper layer, 0–0.1 m), then the second layer (0.1–0.4 m), followed by the third layer (0.4–1.0 m), and finally the bottom layer (1.0–2.0 m). Apparently, the ice in the upper layer melts more quickly than that in the bottom layer of the glacier. Compared with Sens1 estimates, more energy from Sens2 estimates enters the glacier, causing ice to melt more quickly and preferentially, which in turn causes the liquid water content to increase more rapidly and preferentially. Both Sens1 and Sens2 experiments illustrated considerable glacier melting in the ablation region, with an ice thickness reduction of 0.4 m by 22 June for the Sens1 estimates, and an ice thickness reduction of 1.0 m by mid-June for the Sens2 estimates (Figure 7).

**Figure 6.** Comparison between the Sens1 and Sens2 estimates of the daily mean near-surface sensible (HFX, + to atmosphere), latent (LH, + to atmosphere), ground (GRDFLX, + to soil) heat fluxes (**a**) and each component of the surface heat budget (**b**), NetSW: net shortwave radiation; NetLW: net longwave radiation; NetR: net radiation).

**Figure 7.** Comparison between the Sens1 and Sens2 estimates of the liquid water content (SH2O) in the first (L1), second (L2), third (L3), and fourth layer (L4) in June 2016 in the ablation region of the Parlung No. 4 Glacier.

### **4. Discussion and Conclusions**

The glacier mass balance has a grea<sup>t</sup> effect on the physical processes in the cryosphere and hydrosphere, and exhibits strong feedback with current climate change [2,3]. However, high quality observations of glacier mass balance and glacio-meteorological variables are sparse across the TP, while the WRF model can provide useful estimates of meteorological conditions (e.g., T2, precipitation, and radiation fluxes). Previous numerical studies of glacier mass balance estimates usually combine the standard WRF model and the distributed glacier energy and mass balance models through dynamically-downscaled, offline, or interactive coupling simulations [26,27,31,32]. The land-atmosphere coupling WRF involves the mass supplement from precipitation processes and ice melt that is related to the energy budget, which provides a possibility to use WRF alone to study glacier mass changes. However, rarely do studies attempt the related research. In addition, as the key factor of surface net radiation, albedo highly links to the energy balance of glaciers. Previous studies often implement the default glacial albedo scheme when using the standard WRF

to force the distributed glacier energy and mass balance models, seldom improving the glacial albedo scheme in WRF. In our study, control and sensitivity experiments focusing on the land surface type and glacial albedo scheme were conducted using the WRF model coupled with Noah-MP to simulate glacio-meteorological variables in the 2016 ablation season above the Parlung No. 4 Glacier on the TP. The in situ meteorological element observations and the MODIS-retrieved albedo product were selected to evaluate the model performance. Our study initially tries to modify the glacial albedo scheme and evaluates the capacity of WRF alone to the ablation of a marine glacier over the TP. It is one step forward in understanding the glacier mass balance over the TP.

Surfaces with a high albedo reflect a large fraction of solar irradiance, resulting in only a small amount of shortwave radiation being absorbed by the land surface, thereby reducing the surface net radiation heating, and leading to a lower T2 through turbulent heat exchange between the glacier surface and the near-surface atmosphere. In contrast, low albedo has the opposite effect and promotes a higher T2. In this study, the three experiments were able to obtain not only daily variations, but also the mean diurnal pattern of T2. However, the model overestimated the daily mean T2 throughout the ablation season when using the unrealistic open shrub-land surface type, owing to the very low estimated surface albedo of ~0.12. Vegetation decreases the surface albedo mainly because of the low background albedo of the underlying vegetation cover in the CTL experiment. Such a low albedo contributes to a high surface net radiation through reducing the reflected irradiance, and the high net radiation heats the land surface with open shrub-land covered enhancing the turbulent sensible heat exchange among the land surface, vegetation, and near-surface atmosphere. This leads to the CTL estimated mean air temperature of 1.56 ◦C higher than the in situ observations. In addition, the vegetation increases not only the evapotranspiration, changing the surface water balance but also the roughness length and friction velocity enhancing the momentum exchange between the land surface and near-surface atmosphere. The apparent warm bias of 1.56 ◦C was significantly alleviated by sensitivity experiments using a larger glacier albedo, with the mean bias decreasing to −0.5 ◦C and a 37% improvement in RMSE. Consequently, the realistic geostatic land surface parameters such as land surface type are strongly recommended in the model simulations in the following research.

Albedo is a key factor in surface energy balance and affects the redistribution of surface net radiation between turbulent water vapor and heat fluxes. Noah-MP imposes a minimum ice surface albedo of 0.55 which does not adequately represent the surface conditions in the ablation zone of the glacier. This was the main cause of the large albedo value of 0.68 that was simulated in the model that applied the default glacial albedo scheme. A rather high albedo was retained in the default glacial albedo scheme resulting in significant underestimations of the net shortwave radiation and the net radiation. However, an albedo value of 0.35 was held in the modified glacial albedo scheme; this contributed to relieving the underestimation of not only the net shortwave radiation and net radiation, but also the peak net radiation value at local noon. In addition, our simulation results confirm that the net shortwave radiation is the dominant energy source during the ablation season on the southeast TP [20]. Notably, the albedo error of Sens2 in the ablation season is significantly smaller than that of Sens1, but the error of Sens1 before the middle of June and after the middle of September is significantly lower than that of Sens1, revealing that the modified glacial albedo scheme for Sens2 is only applicable to the ablation season. This is because the observed bare ice albedo in the ablation season is set as the background albedo and used in the modified glacial albedo scheme. Besides, the WRF model cannot precisely grasp the variation of albedo during snowfall processes in early June and late September, even though the model applies the observed background albedo to the modified glacial albedo schemes. This is highly linked to the improper snow albedo scheme and the inappropriate description of the relationship between snow albedo and the underlying ice surface.

In addition, the net longwave radiation was apparently overestimated by 17 W m<sup>−</sup><sup>2</sup> in Sens1 and Sens2 experiments. The overestimation of net longwave radiation was highly related to the cloud effect on the downward longwave radiation. The glacier emitted the same intensity of longwave radiation due to the unified observed and simulated surface temperature of 0 ◦C. In realistic cloudy conditions, the cloud increased the downward longwave radiation causing a decrease of the net longwave radiation. However, there is more rainfall in both sensitive experiments (Sens1 and Sens2) on days when no rainfall was observed, and clouds may be present in both sensitive experiments but less matching the timing of the observed clouds. More rainfall in WRF sensitive experiments was caused by cloud particles gradually growing into raindrops through complex microphysical processes and falling to the ground. Afterwards, the cloud dissipated accompanied by a decrease of downward longwave radiation. This would result in an increase of the net longwave radiation for both sensitivity experiments. What is more, the complex topography is a potential factor causing large RMSEs of net radiation and its components for both sensitive experiments; our sensitive experiments that were reported herein did not consider the effects of complex local topography on the direct solar irradiance, and also neglected the reflected and diffuse radiation from surrounding mountainous topography.

Both sensitivity experiments (Sens1 and Sens2) generated relatively small and consistent turbulent water vapor and heat exchange values, showing the similarity of model performance in relation to surface sensible and latent heat, despite two different glacial albedo schemes being applied. The sensible and latent heat fluxes were estimated to be a little higher than Yang's results using the energy balance model [22]. For example, the mean downward sensible heat flux was calculated to be 28 W m<sup>−</sup><sup>2</sup> from our experiments and 17 W m<sup>−</sup><sup>2</sup> from Yang's estimate [22]. Compared with WRF using the default glacial albedo scheme, the excess surface energy that was obtained from WRF using the modified glacial albedo scheme heats surface ice and enters the glacier, heating subsurface ice in the form of ground heat flux, which is estimated to be 182 W m<sup>−</sup>2. This causes surface and subsurface ice melt and the liquid water content to increase more quickly and preferentially in WRF with the modified glacial albedo scheme. Our study confirms net radiation to be the main contributor to surface ice melt, supplemented by sensible heating, and illustrates considerable glacier melting in the ablation region during the ablation season. The ice thickness decreased by 1.0 m by mid-June when applying the modified glacial albedo scheme in WRF, which closely matches the ablation that was measured by stakes and calculated using the surface energy-balance model from Yang's results [22]. However, the ice thickness of glaciers is defined to be 2 m in the current Noah-MP, and the maximum ice thickness reduction is limited to 2 m. Therefore, it is reasonable to model the ablation of glaciers under the condition of ice thickness more than 2 m but the cumulative reduction less than 2 m. In our study, the thickness of the Parlung No. 4 Glacier and the observed cumulative reduction are more than 2 m in the ablation zone in summer 2016 [22]. In order to accurately simulate the large ablation, the ice thickness should be adjusted in the Noah-MP, which was not considered in our study but should be considered further through conducting more numerical tests to ensure the stable calculation during the model's integration. Understanding glacier mass balance anomalies requires quantification and insight into subtle shifts in the energy balance at high altitude glaciers [30]. In order to understand the accumulated and melt regimes of glaciers, further work will focus on the investigation of the roles of each component of glacial surface and subsurface energy balance.

More realistic land-use parameters in the model can reduce the modelled T2 and precipitation biases, and more accurate surface albedo will provide further model bias reductions [56]. This study used the satellite-retrieved spectral albedo and the in situ observed albedo in the ablation zone of the glacier to modify the glacial albedo scheme, which achieved better performance in relation to the glacio-meteorological elements in the focused ablation region. However, the modified albedo scheme is not reasonable in the main area of the glacier because the parameters in the ablation zone cannot represent the parameters in the main area of the glacier. Further improvements in glacial albedo schemes

in land–atmosphere coupling models are urgently required in order to obtain better model performance in the estimation of glacio-meteorological variables across the entire glaciers.

The ablation of glaciers has contributed to rising not only the lake levels in the glacier concentration areas [16], but also the global mean sea level [17,18]. The lake and sea level rise will threaten the lives of people who reside by lakes and in the coastal cities. Therefore, the accurate prediction of glacier ablation contributions to the lake and sea level rise is essential in regions with extensive glacier coverage. Our study assesses the ability of WRF itself to estimate the ablation of a mountainous glacier, which is the first step before WRF itself can be widely used in the accurate simulation and prediction of glacier ablation. Our study confirms the potential of the WRF model to generate reasonable glacio-meteorological variables, such as T2, radiation, and ground heat flux except precipitation, as long as a reasonable glacial albedo scheme is applied. These variables are essential to estimate glacier energy and mass changes. Precipitation is another decisive factor in glacier mass changes and different phases of precipitation (i.e., rain and snow) show different mechanisms of glacier mass balance; however, the simulated precipitation was twice that of the in situ observations during the ablation season over the glacier. This is not directly influenced by the glacial albedo scheme but was affected not only by the microphysical and cumulus parameterization schemes in the WRF model, but also by the complex mountainous topography. At present, we have not investigated complex parameters in the microphysics and cumulus schemes, nor the role of terrain and the terrain-induced water vapor transport, which relate to the very high precipitation duration, intensity, and amount in WRF estimates. Therefore, it is important to explore how to improve the parameterizations of precipitation in regional models over complex topography. In the future, we will focus on evaluating and improving the performance of the model with regard to the total amount of summer precipitation over the southeastern TP, through investigation of the model microphysics and cumulus schemes. Also, we will further investigate the roles of mountainous topography and the terrain-induced water vapor from the Indian summer monsoon in increasing precipitation over the southeastern TP.

**Author Contributions:** Conceptualization, L.L. and M.M.; methodology, L.L.; writing original draft preparation, L.L., M.M. and Y.M.; funding acquisition, Y.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Strategic Priority Research Program of Chinese Academy of Sciences (XDA20060101), the Second Tibetan Plateau Scientific Expedition and Research program (STEP) (2019QZKK0103), and the National Natural Science Foundation of China (91837208).

**Data Availability Statement:** The ERA-Interim dataset is available from http://apps.ecmwf.int/ datasets/data/interim-full-daily/, accessed on 1 July 2020. Meteorological observations on Parlung No. 4 Glacier are available from National Tibetan Plateau Data Center (https://data.tpdc.ac.cn/en/, accessed on 1 July 2020).

**Acknowledgments:** The first author would like to acknowledge all group members for their help in completing this paper.

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