**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<sup>2</sup> and estimated water reserves of 5600 km<sup>3</sup> (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–10]. Under global warming, glaciers on the TP have been retreating and shrinking for decades, a trend that has accelerated in recent years [11–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 grea<sup>t</sup> 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–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–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 dynamicallydownscaled, 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–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<sup>−</sup><sup>2</sup> 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<sup>2</sup> 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–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.

**Figure 1.** (**a**) WRF domains and model topography. (**b**) Terrain elevation from the WRF model, shaded in units of meters with the green line denoting the glacier boundary and the green solid circle denoting the observation site in the ablation zone of 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 glaciometeorological 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 agreemen<sup>t</sup> 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 threehourly 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.


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

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 + (\alpha\_{old} - 0.55)e^{\frac{-0.01dt}{\text{score}}} \tag{1}$$

$$f\_{su} = \tanh(\frac{h\_{su}}{2.5z\_0(\frac{\rho\_{sn}}{\rho\_{nuv}})^{f\_m}}) \tag{2}$$

$$a\_{\\$} = a\_1 + f\_{\text{sil}}(0.84 - a\_1) \tag{3}$$

$$
\mathfrak{a}\_{sd1} = \mathfrak{a}\_{sd2} = \mathfrak{a}\_{s1} = \mathfrak{a}\_{s2} = \mathfrak{a}\_s \tag{4}
$$

where *αold* is the albedo of the last time step (*dt*). *fsn* is fractional snow cover. *hsn* is snow depth in unit of m. *ρsn* is the bulk density of snow in unit of kg m<sup>−</sup>3. *ρnew* is the fresh snow density with the value of 100 kg m<sup>−</sup>3. *z*0 is the snow surface roughness length with the value of 0.002 m. *fm* is melting factor determining the curves in melting season which is adjustable and sets to 1.0 in Noah-MP. *αs* is the albedo of snow. *αsd*1 and *αsd*2 denote the direct albedo of snow for visible and near infrared bands, respectively, and *αsi*1 and *αsi*2 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\tag{5}$$

$$\mathfrak{a}\_{\rm si1} = 0.95(1 - 0.2A\_{\rm \varepsilon}) \tag{6}$$

$$a\_{\rm si2} = 0.65(1 - 0.5A\_{\rm c})\tag{7}$$

$$
\alpha\_{\rm sd1} = \alpha\_{\rm si1} + 0.4Z\_{\rm c}(1 - \alpha\_{\rm si1}) \tag{8}
$$

$$
\alpha\_{\text{sd2}} = \alpha\_{\text{s2}} + 0.4Z\_{\text{c}}(1 - \alpha\_{\text{s2}}) \tag{9}
$$

where *Z* is the solar zenith angle and *Ac* 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 *αs* × 1.2 in Equation (6) and 0.65 to *α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.


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

### *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]:

$$\mathbf{a}\_{\text{short}} = 0.160\mathbf{a}\_1 + 0.291\mathbf{a}\_2 + 0.243\mathbf{a}\_3 + 0.116\mathbf{a}\_4 + 0.112\mathbf{a}\_5 + 0.081\mathbf{a}\_7 - 0.0015 \tag{10}$$

where *αshort* is the surface broadband albedo, and *α*1 to *α*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.
