Evaluation of Albedo Schemes in WRF Coupled with Noah-MP on the Parlung No. 4 Glacier

Abstract

Meteorological variables (e.g., air temperature (T2), radiation flux, and precipitation) determine the evolution of glacier mass and characteristics. Observations of these variables are not available with adequate spatial coverage and spatiotemporal resolution over the Tibetan Plateau. Albedo is the key factor of net radiation and is determined by the land cover and snow-related variables. This study focuses on evaluating the performance of the albedo parameterization scheme in WRF coupled with Noah-MP in terms of glacio-meteorological variables, by conducting experiments applying the standard surface albedo scheme with the default vegetation and corrected to ice cover and the modified glacial albedo scheme to the Parlung No. 4 Glacier in the 2016 ablation season. In situ glacio-meteorological element observations and MODIS-retrieved albedo are selected to assess the performance of the model. The key results are as follows. First, compared to the air temperature bias of 1.56 °C in WRF applying the standard surface albedo scheme and the default vegetation cover, realistic land-use categories considerably reduce the model warm bias on the glacier. The model using realistic land-use categories yields similar T2 diurnal patterns to the observations, with a mean bias of only −0.5 °C, no matter which glacial albedo scheme is implemented. Second, the default glacial albedo scheme gives a rather high albedo value of 0.68, causing an apparent underestimation of the net shortwave radiation and net radiation; the modified glacial albedo scheme gives a mean albedo value of 0.35, close to the in situ observations, helping to relieve underestimations of net shortwave radiation and net radiation. Compared with the MODIS albedo of the glacier, WRF applying the default glacial albedo scheme apparently overestimates the albedo with a mean error of 0.18, while WRF applying the modified glacial albedo scheme slightly underestimates the albedo with a mean error of only −0.08. Third, the mean net radiation flux (142 W m⁻²) and high ground heat flux (182 W m⁻²) values that were estimated by WRF applying the corrected land cover and the modified glacial albedo scheme result in the heating of the glacier surface and subsurface, causing ice melt and the liquid water content to increase more quickly and preferentially, equating to an estimated ice thickness decrease of 1 m by mid-June in the ablation region. Our study confirms the ability of the WRF model to reproduce glacio-meteorological variables as long as a reasonable glacial albedo scheme and the corrected land cover is applied and provides a theoretical reference for researchers that are committed to further improvement of the glacial albedo scheme.

1. Introduction

The Tibetan Plateau (TP) has a complex topography and unique geographical environment, with a mean elevation of approximately 4000 m above sea level (a.s.l.). It is known as the Earth’s ‘third pole’ [1], containing many of the world’s middle- and low-latitude glaciers. According to the Second Chinese Glacier Inventory, there are 48,571 glaciers in China, with a total area of 51,480 km² and estimated water reserves of 5600 km³ (http://news.sciencenet.cn/htmlnews/2015/1/310736.shtm, accessed on 1 June 2020), of which approximately 80% are found on the TP. Glaciers represent an important land surface type, and their glacier–atmosphere interaction affects the exchange of water and energy in the land–atmosphere system. Very complex physical feedback mechanisms link glaciers and the climate system [2,3]. Therefore, glaciers are considered to be natural indicators and sensitive recorders of climatic and environmental changes [4,5].

As a significant component of the cryosphere, mountain glaciers have attracted unprecedented attention, in particular in regard to their mass balance [6,7,8,9,10]. Under global warming, glaciers on the TP have been retreating and shrinking for decades, a trend that has accelerated in recent years [11,12,13,14]. It has been noted that glaciers that are located in the southeastern TP and central Himalayas have retreated rapidly, while those that are located in the Karakoram and Eastern Pamirs have retreated slowly, revealing the great spatial variability in glacier mass balance across the whole TP [15]. Glacier mass balance change has an important impact on the availability of glacial meltwater to recharge the surrounding rivers and lakes of the Yangtze River basin. The retreat of glaciers has contributed to rising lake levels in regions with extensive glacier coverage, such as the Nam Co Lake and Selin Co Lake [16], and has contributed to global mean sea level rise [17,18].

Glacier mass balance has been observed sparsely and far from comprehensively over the topographically-complex TP. Previous research has been limited to a small number of glaciers, including the Qiyi, Xiaodongkemadi, and Parlung glaciers [15,16]. Most investigations of glacier mass balance have depended on energy-based models [6,10] and remote sensing retrievals [9,19]. Glacio-meteorological variables (i.e., near-surface air temperature, precipitation, wind speed, relative humidity, and radiation fluxes) greatly affect the glacier mass balance and are essential factors in mass balance models. The glacial meteorology, point energy, and mass balance of Parlung No. 4 Glacier has previously been investigated [20,21,22,23], revealing that net radiation fluxes (especially net shortwave radiation) govern the surface melt of the glacier, with net shortwave radiation contributing 98% of the surface melt. The temperature index model has been proven to be applicable for mass balance and ablation modeling when incorporating solar irradiance [20]. Modeling the spatial distribution of glacier mass balance requires distributed glacio-meteorological forcing, but this is difficult to implement owing to the sparse and uneven distribution of in situ observation stations across the TP. Also, collecting valid in situ measurements of glacio-meteorological variables and energy balance is difficult owing to logistical problems that are associated with the harsh, high-elevation environment of the TP.

Land–atmosphere interactions are evident at the interface of glaciers and the lower atmosphere and drive the rapid response of glaciers to surrounding environmental changes. Temperature, precipitation, and general atmospheric circulation are essential factors influencing changes in glacier mass [15,24]. As one of several coupled atmosphere–land surface models (LSMs), the advanced Weather Research and Forecasting (WRF) model [25] is a good candidate for estimating the glacio-meteorological variables that are required to force glacier mass balance models (e.g., a distributed energy and mass balance model). Numerous studies have evaluated the ability of WRF to produce forcing data for glaciological studies with a correct representation of the glacierized area [26,27,28,29,30]. Great efforts have been made to estimate glacio-meteorological variables using WRF coupled with an LSM (e.g., Noah, multi-physics Noah (Noah-MP), and Rapid Update Cycle (RUC)). These schemes have been used to drive the physically-based, distributed glacier energy and mass balance models that were developed to estimate mountain glacier mass variability in dynamically-downscaled, offline, or interactive coupling simulations [26,27,31]. Provisional results have indicated that WRF-modelled meteorological variables at high spatial resolution (i.e., 1 km) can be used to force distributed simulations of Kersten Glacier mass balance with acceptable accuracy [26]. This method has also been successful in simulating the Zhadang Glacier, a small alpine glacier, although feedbacks from the glacier surface mass change into the regional atmospheric forcing were neglected [31]. Models of the interactive coupling between WRF and glacier mass balance have shown promise in studying glacier mass variability [27,32]. However, previous studies of glacier mass balance are seldom based only on the WRF model although the mass accumulation and ice melt as well as energy budget had been involved in this model, leaving more possibility in future research.

The land surface type is a significant factor affecting land surface properties (e.g., emissivity, albedo, and roughness length), and in WRF has an important influence on the modeled surface and near-surface meteorological variables (e.g., temperature, radiation, albedo, wind speed, and snowmelt). However, the statistical land-use product in WRF is wrong to match the Parlung No. 4 Glacier land cover. In addition, snow albedo determining the surface energy budget and influences the glacier mass balance, undergoes large variations during the snow melting and accumulation periods, which is essential for the ice surface energy and mass balance because of its strong controls on the length of the accumulation and ablation seasons. It is significantly affected by many parameters, i.e., snow depth and age, snow cover, surface temperature, cloud cover fraction, wind speed, positive accumulated degrees days, solar zenith angle, and impurities [33,34,35,36,37]. Generally, the snow albedo schemes depend on the observation data and involve the empiric parameters with the most important to be the maximum and minimum albedo. From the review of the currently existing snow albedo schemes [33], many glacial albedo schemes use the minimum albedo about 0.5, which is mostly suitable for the thick ice but not suitable for relatively thin ice. The maximum prescribed snow albedo is usually 0.8–0.85, but the fresh snow albedo is observed up to 0.95. What’s more, the simplest snow albedo schemes apply constant values of albedo for different land cover. Other schemes depend on temperature to account for snow metamorphism and snow thinning. More sophisticated schemes consider the snow-related variables and solar zenith angle (for example Biosphere-Atmosphere Transfer Scheme (BATS) [38] and LSM [39]) and impurities [40,41]. It is revealed that the deposition of absorbing aerosols decreases the snow albedo of the Himalayan region by 0.15 ± 0.13, causing a positive radiative effect of 14 ± 13 W m⁻² and an increase of the surface temperature by 1.33 ± 1.2 °C as well as the reduction of the snow cover fraction by 7 ± 11% [42]. Therefore, the choice of the snow albedo scheme has a considerable impact on the simulations of both weather and climate and the glacier energy and mass balance.

In addition, glaciers in the Himalayas are mostly sensitive to monsoon-related precipitation perturbations and summer air temperature, which are highly linked to albedo owing to the crucial snow–albedo feedback in summer [43]. In this study, we conducted three numerical experiments using WRF, with one applying the default land surface type (open shrub-land) and the other two adjusted to the snow and ice type and modified albedo scheme on Parlung No. 4 Glacier in the ablation season. Glacio-meteorological variables (near-surface air temperature, wind speed, precipitation, albedo, and radiation fluxes) were modeled and evaluated based on in situ observations and satellite-retrieved albedo at the glacier. This preliminary work is helpful in assessing the two-way coupling of WRF and glacier mass balance models when estimating the mass change of maritime mountain glaciers.

2. Data and Methodology

2.1. Study Area and In-Situ Measurements

Parlung No. 4 Glacier is located at the southeast margin of the TP (Figure 1). This geographic region is strongly dominated by the South Asian monsoon and receives frequent precipitation after the onset of the Indian summer monsoon; thus, glaciers here are of the maritime type. It is a debris-free glacier with an area of 11.7 km² and a length of 8 km [22]. Many high-quality glacio-meteorological and mass balance observations are available for Parlung No. 4 Glacier, and these have allowed detailed study of the glacier [20,21,22,23]. The glacial mass change in the ablation zone is a crucial component of the whole glacier mass balance. Therefore, this study focuses on the assessment of the glacio-meteorological variables simulations in the ablation zone of the Parlung No. 4 Glacier.

An automatic weather station (AWS) is installed at 29.25°N, 96.93°E, at an elevation of 4800 m a.s.l. in the ablation zone of Parlung No. 4 Glacier (Figure 1). Specific meteorological variables including air temperature at 2 m height (T2), and upward and downward shortwave and longwave radiation, are collected hourly by a CR1000 Campbell Scientific data logger. The hourly precipitation is measured by a Geonor T-200B weighing bucket gauge. The observed T2, components of the radiation fluxes, and precipitation in summer were used to evaluate the numerical estimates in our three experiments. These observational data were obtained from the TP scientific data center website and are freely available from the National Tibetan Plateau Data Center (https://data.tpdc.ac.cn/en/, accessed on 1 July 2020). The observed surface albedo was calculated as the ratio of upward shortwave radiation flux to solar irradiance. Detailed instrumental information and descriptions of the effects of the observed meteorology and surface energy fluxes on the glacier in the ablation season have already been provided [20,22].

2.2. Model Configuration and Experimental Design

The WRF model was developed through a partnership of the National Center for Atmospheric Research (NCAR), the National Oceanic and Atmospheric Administration (NOAA), the United States Air Force, the Naval Research Laboratory, the University of Oklahoma, and the Federal Aviation Administration (The model can be downloaded from https://www2.mmm.ucar.edu/wrf/users/, accessed on 1 July 2020). It is a state-of-the-art atmospheric modeling system, comprising of a fully compressible and non-hydrostatic model with a terrain-following pressure vertical coordinate and Arakawa C-grid horizontal coordinates [25]. It uses Runge–Kutta second- and third-order integration in the time schemes, and second- to sixth-order integration in the advection schemes. The options for atmospheric and land surface processes can be chosen with a broad range of grid sizes, from tens of meters to thousands of kilometers.

WRF has the capacity to estimate glacio-meteorological variables in the low boundary layer and can successfully force distributed energy and mass balance models of mountain glaciers [27,31]. In order to investigate the performance of WRF with respect to glacio-meteorological variables above maritime mountain glaciers in the ablation season, the more recent WRF version 4.3.1 was applied to Parlung No. 4 Glacier in summer 2016. We configured 3 nested domains, with the inner-most domain covering the glacier and its surroundings (Figure 1). The model was centered on 29.23°N, 96.92°E, with spatial resolutions of 12.5, 2.5, and 0.5 km. All the domains were set to 50 terrain-following vertical levels, stretching from the surface to 50 hPa. A dataset from the Interim European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-Interim) [44], which is found to have the closest agreement with in situ measurements of air temperature in the Tibetan Plateau [45], was chosen to produce the initial and boundary meteorological conditions with a 0.25° × 0.25° horizontal resolution and a six-hourly intervals. The Noah-MP LSM includes a separate glacier treatment and improved snow physics, with up to three layers in the snowpack, representing improvements over the original Noah scheme [46,47]; it also features a modified two-stream radiation transfer scheme, which considers the three-dimensional canopy structure to calculate radiation fluxes that are reflected, absorbed, and transmitted by vegetation. This LSM uses a ‘tile’ approach to calculate albedo, a key factor in the energy budget, considering bare ground, vegetation canopy, and snow cover [48]. The Noah-MP coupled with WRF has been shown to provide suitable robust precipitation estimates across the TP [49]; hence, this scheme was chosen for our current study. The model was run from 1 May to 1 October 2016, producing three-hourly output meteorological variables. The first month of simulation was regarded as the model spin-up. The physics schemes that were selected when using WRF and the multi-physical parameterization schemes from Noah-MP are detailed in

Table 1. Detailed options that were selected in WRF coupled with the Noah-MP LSM.

Simulations Time Period	From 1 May to 1 October 2016
Nest	Three nested domains (two-way)
Projection	Mercator
Center of domain	29.23°N, 96.92°E
Resolution	12.5 km, 2.5 km, 0.5 km
Microphysics	Thompson scheme
Longwave radiation	RRTM scheme
Shortwave radiation	Dudhia scheme
Surface layer	Revised MM5 Monin-Obukhov scheme
Planetary boundary layer	YSU scheme
Cumulus parameterization	Grell-Devenyi ensemble scheme in domain 1
Land surface	Noah-MP
Canopy stomatal resistance	Ball-Berry scheme
Soil moisture factor for stomatal resistance	Noah scheme
Runoff and groundwater	TOPMODEL with groundwater
Surface layer drag coefficient	Monin–Obukhov
Soil permeability option	Linear effect, more permeable
Radiation transfer	Two-stream applied to vegetated fraction
Ground surface albedo option	CLASS (Canadian Land Surface Scheme)
Precipitation partitioning between snow and rain	Jordan scheme
Snow/soil temperature time scheme	Semi-implicit
Noah-MP glacier treatment	Includes phase change

.

Land-use type strongly influences radiation fluxes and near-surface air temperatures [50]. The default static land-use in the current WRF version 4.3.1 is the Moderate Resolution Imaging Spectroradiometer (MODIS) land-use product with 30 arc-seconds spatial resolution. This land-use product incorrectly classifies land cover within the ablation zone of Parlung No. 4 Glacier as open shrub-land. To evaluate the importance of land surface type and albedo-related parameters in the accurate estimation of meteorological elements, numerical experiments were conducted: the first experiment used the default land-use product and default albedo scheme of CLASS [51] as the control experiment (CTL); the second experiment (Sens1) used the true land-use type of snow and ice for the extent of the glacier, default bare ice albedo (visible = 0.8, near infrared = 0.55, background albedo = 0.55), and snow albedo scheme of CLASS; the third experiment (Sens2) was similar to Sens1, but corrected bare ice albedo according to previous results (visible = 0.5, near infrared = 0.2 [52], background albedo = 0.23 [22]), and additionally included snow age and solar zenith angle in the CLASS scheme according to the principle of the BATS snow albedo scheme [53]. The CLASS (Equations (1)–(4)) snow albedo scheme and snow cover fraction that were used in the CLASS scheme are described in the Sections 2 and 3.4 of the technical description of Noah-MP [53] in the following equations:

$${\alpha }_1=0.55+\left({\alpha }_{old}-0.55\right)e^{\frac{-0.01dt}{3600}}$$

(1)

$$f_{sn}=\tanh (\frac{h_{sn}}{2.5z_0{(\frac{{\rho }_{sn}}{{\rho }_{new}})}^{f_m}})$$

(2)

$${\alpha }_s={\alpha }_1+f_{sn}\left(0.84-{\alpha }_1\right)$$

(3)

$${\alpha }_{sd1}={\alpha }_{sd2}={\alpha }_{si1}={\alpha }_{si2}={\alpha }_s$$

(4)

where ${\alpha }_{old}$ is the albedo of the last time step ($dt$). $f_{sn}$ is fractional snow cover. $h_{sn}$ is snow depth in unit of m. ${\rho }_{sn}$ is the bulk density of snow in unit of kg m⁻³. ${\rho }_{new}$ is the fresh snow density with the value of 100 kg m⁻³. $z_0$ is the snow surface roughness length with the value of 0.002 m. $f_m$ is melting factor determining the curves in melting season which is adjustable and sets to 1.0 in Noah-MP. ${\alpha }_s$ is the albedo of snow. ${\alpha }_{sd1}$ and ${\alpha }_{sd2}$ denote the direct albedo of snow for visible and near infrared bands, respectively, and ${\alpha }_{si1}$ and ${\alpha }_{si2}$ denote the diffuse albedo of snow for visible and near infrared bands, respectively.

The BATS (Equations (5)–(9)) snow albedo scheme is described in the Section 3.3 of the technical description of Noah-MP [53] in the following equations:

$$Z_c=\frac{1.5}{1+4\cos Z}-0.5$$

(5)

$${\alpha }_{si1}=0.95\left(1-0.2A_c\right)$$

(6)

$${\alpha }_{si2}=0.65\left(1-0.5A_c\right)$$

(7)

$${\alpha }_{sd1}={\alpha }_{si1}+0.4Z_c\left(1-{\alpha }_{si1}\right)$$

(8)

$${\alpha }_{sd2}={\alpha }_{si2}+0.4Z_c\left(1-{\alpha }_{si2}\right)$$

(9)

where Z is the solar zenith angle and $A_c$ is the snow age.

The default snow albedo scheme in the model is developed based on the deep snow with slow melting, which shows a large snow-related simulation deviation on the TP where the snow is shallow and melts rapidly. The insufficient consideration of snow age leading to the lag of melting is the potential reason. In order to more accurately account for the impact of snow age on snow melting on the TP, we attempted to simultaneously consider the snow age in both CLASS and BATS schemes in the Sens2 experiment. However, the albedo of snow for visible and near infrared bands is parameterized to the same value in the CLASS scheme. In reality, the spectral albedo is the different values according to the spectral albedo measurements with a higher albedo for visible band and a lower for near infrared band and with different spectral albedo curves for fresh snow and old snow [54]. In terms of spectral albedo that is related to snow age in the BATS scheme, multiplicative factors of 0.95 in Equation (6) and 0.65 in Equation (7) represent the diffuse fresh snow albedo for visible and near infrared bands, respectively, which corresponds to about 1.2 and 0.8 times the prescribed fresh snow albedo value of 0.8 for broadband in the model. Therefore, we boldly modified multiplicative factors of 0.95 to ${\alpha }_s$ × 1.2 in Equation (6) and 0.65 to ${\alpha }_s$ × 0.8 in Equation (7) when integrating the CLASS and BATS snow albedo schemes in Sens2 experiment. Eventually, the Equations (1)–(3) and (5), the modified Equations (6)–(9) were used ordinally to calculate the spectral snow albedo for the direct and diffuse irradiance in the Sens2 experiment.

The different initial surface conditions and the applied snow albedo schemes in the ablation zone of the glacier in three experiments are outlined in

Table 2. Initial surface conditions in the ablation zone of the glacier and snow albedo schemes that were used in our experiments (vis = visible, nir = near infrared).

Experiment	Land-Use Categories	Vegetation Fraction/%	Soil Categories	Background Albedo	Bare Ice Albedo	Snow Albedo Scheme
CTL	open shrub-land	30.5	loam	0.16	-	CLASS
Sens1	snow and ice	0	land-ice	0.55	vis = 0.8, nir = 0.55	CLASS
Sens2	snow and ice	0	land-ice	0.23	vis = 0.5, nir = 0.2	combined CLASS and BATS

.

2.3. Evaluation of Model Performance

In order to assess the performance of the model in the ablation region of the glacier, the root-mean-square error (RMSE), mean absolute deviation, and the correlation coefficient (CC) between the in situ observations and model estimates were used to evaluate the model performance in terms of the glacio-meteorological variables (T2 and shortwave/longwave radiation). Also, we applied linear regression to the observed and modeled net radiation. Pearson linear cross-correlation was chosen to calculate the CC, and the t-test was chosen to test the significance of the correlation. A significance level of 0.01 was specified in this study. In addition, space remote sensing instrument developed by National Aeronautics and Space Administration, MODIS is used to monitoring global climate change. The product of MOD09GA Version 6 (available from: https://lpdaac.usgs.gov/products/mod09gav006/, accessed on 1 July 2020) provides an estimate of the daily surface spectral reflectance of MODIS/Terra in bands 1 to 7 with a spatial resolution of 500 m. This product was used to assess the performance of the model across the whole glacier, including the main body of the glacier. Quality assurance information regarding the quality control code and atmospheric condition flag of the MOD09GA product were considered to achieve the ideal quality (ideal quality of bands, no cloud, low aerosol quantity) broadband albedo product, which can be estimated using Liang’s algorithm [55]:

$${\alpha }_{short}=0.160{\alpha }_1+0.291{\alpha }_2+0.243{\alpha }_3+0.116{\alpha }_4+0.112{\alpha }_5+0.081{\alpha }_7-0.0015$$

(10)

where ${\alpha }_{short}$ is the surface broadband albedo, and ${\alpha }_1$ to ${\alpha }_7$ represent the surface reflectance in MODIS bands 1 to 7, respectively. The spectral coverage for MODIS bands 1 to 7 is 0.62–0.67, 0.84–0.87, 0.46–0.48, 0.54–0.56, 1.23–1.25, 1.63–1.65, and 2.11–2.15 μm, respectively.

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).

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).

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⁻² and the mean estimates reached 30.56 W m⁻² 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⁻². 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⁻², respectively. The average net radiation was calculated to be 123.44 W m⁻² using in situ observations, 53.64 W m⁻² for Sens1 estimates, and 141.85 W m⁻² 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⁻² 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⁻²) and Sens2 (503.99 W m⁻²), the latter being closer to the observed value of 454.62 W m⁻² (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.

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 (Figure 4 and Figure 5a). On average, surface shortwave radiative forcing reached to −141 (−60) W m⁻² 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⁻² 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⁻²) and latent (−51 to 36 W m⁻²) 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⁻²) and Sens2 (30–295 W m⁻²) 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⁻², which was much lower than that from Sens2 (182 W m⁻²). 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).

4. Discussion and Conclusions

The glacier mass balance has a great 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⁻² 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⁻² from our experiments and 17 W m⁻² 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⁻². 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.