-
A severe snow event took place over a large fraction of the TP in March 2017 with a maximum SD deeper than 70 cm in Nyalam County of the Tibetan Autonomous Region, China. This extreme event featured ground observations of SD never recorded before in the eastern TP (i.e., in Minyang), and within the thickest 4th percentile in the eastern and southern TP —Nyalam and Songpan — by historical statistics of SD observations for more than 50 years. Variable snowfall intensity was observed, from light snow in some regions to heavy snow in the eastern, central, and southern edges of the TP. The sensitivity of this snow event to the applied LSM and initial and boundary conditions was investigated in our previous study (Liu et al., 2019).
This study focuses on the same snow event but includes an improved albedo parameterization and an updated set of static land surface parameters in WRF + Noah LSM, such as land cover and fractional vegetation cover (FVC). Observations of near-surface air temperature and snow water equivalent (SWE) were taken from 32 World Meteorological Organization surface synoptic observation (SYNOP) stations run by the China Meteorological Administration (CMA), and albedo observations were taken from six observatories operated by the Chinese Academy of Sciences (CAS). The geographic information and distribution of in situ stations are shown in Table 1 and Fig. 1.
Station Type Station number Name ID Latitude (°N) Longitude (°E) Elevation (m) CMA 1 Mangya 51886 38.25 90.85 2945.0 2 Daqaidam 52713 37.85 95.37 3174.0 3 Gangca 52754 37.33 100.13 3302.0 4 Golmud 52818 36.42 94.90 2807.6 5 Doulan 52836 36.30 98.10 3189.0 6 Xining 52866 36.62 101.77 2295.2 7 Shiquanhe 55228 32.50 80.08 4278.6 8 Baingoin 55279 31.37 90.02 4700.0 9 Nagqu 55299 31.48 92.07 4507.0 10 Xainza 55472 30.95 88.63 4672.0 11 Xigaze 55578 29.25 88.88 3836.0 12 Lhasa 55591 29.67 91.13 3648.9 13 Tingri 55664 28.63 87.08 4300.0 14 Lhunze 55696 28.42 92.47 3860.0 15 Pagri 55773 27.73 89.08 4300.0 16 Tuotuohe 56004 34.22 92.43 4533.1 17 Zadoi 56018 32.90 95.30 4066.4 18 Qumarleb 56021 34.13 95.78 4175.0 19 Yushu 56029 33.00 96.97 3716.9 20 Madoi 56033 34.92 98.22 4272.3 21 Darlag 56046 33.75 99.65 3967.5 22 Zoige 56079 33.58 102.97 3441.4 23 sog 56106 31.88 93.78 4024.0 24 Dengqen 56116 31.42 95.60 3873.1 25 Qamdo 56137 31.15 97.17 3315.0 26 Sertar 56152 32.28 100.33 3894.0 27 Barkam 56172 31.90 102.23 2664.4 28 Songpan 56182 32.67 103.60 2850.7 29 Batang 56247 30.00 99.10 2589.2 30 Nyingchi 56312 29.57 94.47 2991.8 31 Deqen 56444 28.45 98.88 3319.0 32 Jiulong 56462 29.00 101.50 2925.0 CAS 33 QOMS 28.36 86.95 4276.0 34 SETS 29.77 94.73 3326.0 35 NASDE 33.39 79.70 4264.0 36 MASWE 38.41 75.05 3668.0 37 SHSEX 33.22 88.83 4947.0 38 MAQU 33.92 102.10 3440.0 Note: Qomolangma station (QOMS), Southeast Tibet station (SETS), Ngari station (NASDE), Muztagh Ata station (MASWE), ShuangHu station (SHSEX), MAQU station (MAQU) Table 1. Information of in situ China Meteorological Administration (CMA) and Chinese Academy of Sciences (CAS) stations.
Figure 1. Land cover types and in situ stations over the Tibetan Plateau and the surrounding regions; the legend number indicates the land cover type index, with the associated specific description to the right of each index; solid red circles represent the location of in situ stations with the station number near its position. For a description of station numbers, see Table 1.
Due to the sparse and uneven distribution of the few CAS observatories, the in situ observations of albedo are poorly representative of the TP. Spaceborne observations of radiance reflected by the Earth-Atmosphere system provides an effective means to monitor land surface albedo both regionally and globally (Li and Garand, 1994), and MODIS albedo products are more accurate than other satellite data products when evaluated against in situ observations (Liang et al., 2002, 2005). The products MOD09A1, MOD09CMG, and MYD09CMG were used to calculate surface albedo by converting narrowband reflectance to broadband albedo following Liang (2000). The product MOD09A1 is a MODIS/Terra eight-day surface reflectance product with a spatial resolution of 500 m, provided on a sinusoidal coordinate grid. The product MOD09CMG (MYD09CMG) is a MODIS/Terra (MODIS/Aqua) daily surface reflectance product with a spatial resolution of 0.05°, provided on the Climate Modeling Grid (CMG). The three surface reflectance products are L3 data products and provide surface reflectance in bands 1 to 7, from which broadband albedo can be estimated (Liang, 2000):
where
${\alpha _{{\rm{short}}}} $ is the surface broadband albedo;$ \alpha_1 -\alpha_7$ is the surface reflectance in MODIS bands 1 to 7. 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.The dataset MCD12Q1 is the MODIS Terra + Aqua combined land cover product. It is a yearly L3 level product with a spatial resolution of 500 m, provided on a sinusoidal coordinate grid. It contains five land cover classification schemes and is calculated using a supervised decision tree classification algorithm. The primary land cover scheme identifies 17 land cover classes defined by the International Geosphere-Biosphere Program (IGBP). Due to large differences in reflected radiation for different land cover types, the land cover 2017 product from MODIS was used to calculate land surface albedo using the IGBP classification in this study.
-
The default surface albedo parameterization scheme in the Noah LSM considers snow cover and age (Ek et al., 2003; Livneh et al., 2010). The Noah LSM default albedo scheme is described by:
where
$ L_v$ is a prescribed empirical parameter;$ \alpha_{s0}$ is a prescribed maximum albedo for each land cover type; t is the number of days since the last snowfall; A and B are constants set to 0.94 and 0.58, respectively for the snow accumulation period, i.e., from October to February, and set to 0.82 and 0.46 otherwise;$ \alpha_{\rm{bg}}$ is the background albedo, and$ f_{\rm c} $ is the fractional snow cover. Snow cover is set equal to 1.0 when SWE meets or exceeds the prescribed threshold SWE. Otherwise, snow cover is parameterized using the following equations:where
$S_{{\rm{up}}}$ is a prescribed threshold SWE in meters that implies 100 percent snow cover for each land cover type;$S_{{\rm{eqv}}}$ is SWE in meters; and$ P_0$ is a constant shape parameter of the distribution function of snow cover (set to 2.6 in the Noah default scheme).The Noah LSM default albedo scheme (Livneh et al., 2010) imposes a seasonally-variable decay in albedo from its initial fresh snow value. The decay is slower (faster) during the accumulation (ablation) season [Eq. (3)], with an initial maximum snow albedo value that is usually very high [Eq. (2)]. Apart from snow surface albedo, this scheme also includes the influence of ground surface albedo through the fractional snow cover [Eq. (4)], which leads to some shortcomings. For example, snow cover in the Noah LSM is calculated using an areal depletion threshold (
$S_{{\rm{up}}}$ ), which assumes partial snow cover when SWE is below the threshold [Eqs. (5) and (6)]. Otherwise, snow cover is set equal to 1.0. The parameter$S_{{\rm{up}}}$ for non-forest land cover is fixed at 0.02 m. In summary, the areal snow cover fraction remains at 1.0, and the albedo is only determined by snow age while the snow is thicker than 0.1 m (Ek et al., 2003; Livneh et al., 2010). This approach is independent of the influence of ground albedo. Therefore, regarding the greater sensitivity to SD than snow cover for the snow accumulation and ablation periods, considering SD explicitly in the Noah LSM albedo scheme may be a better choice.To improve the Noah LSM albedo scheme, we used the snow albedo parameterization proposed by Oerlemans and Knap (1998), which takes into account SD and snow age explicitly, and includes an application specific to glaciers:
where s is the last day with snowfall prior to day i, and s–i, thus giving the snow age in days. The parameters
$\alpha_{\rm{firn}} $ ,$\alpha _{\rm{freshsnow}}$ , and$ \alpha_{\rm{ice}} $ are the albedos of firn, fresh snow, and bare ice, respectively; d is SD in meters; t* and d* are scaling parameters for the snow age and SD, respectively;$ \alpha^{(i)}$ is actual albedo on day i. Equation (9) determines whether or not snowfall has occurred.Land cover in the TP is fragmented (Fig. 1) and the snow-free surface albedo depends on land cover and location. Our first change to the albedo scheme was to replace
$\alpha_{\rm{ice}} $ in Eq. (8) with the snow-free albedo appropriate to the land cover at a given location. This change accounts for the large differences in reflected radiance between different land cover types (Dong and Li, 1994; Hu et al., 2019). For example, the reflectance of a water body is much lower than the reflectance of bare soil or a sparsely vegetated surface (Menenti et al., 1989). -
The non-hydrostatic WRF model (Skamarock et al., 2008), version 3.7.1 (released August 2015), was used for this study with a horizontal resolution of 25 km. The WRF experiments were configured for a single large domain, with an upper-right boundary at 46°N and 110°E, and a lower-left boundary at 20°N and 60°E, to fully include regions important to the Indian monsoon and westerlies. We ran the model for 10 days and 18 hours starting at 0800 Beijing Standard Time (LST) on 5 March 2017, which was at least one day before snowfall. This temporal gap allowed the model ample time to stabilize and rapidly converge against a physically dynamic balanced state. The model initial and boundary conditions were provided by the European Centre for Medium Range Weather Forecasts (ECMWF) reanalysis dataset (ERA-Interim) provided at a 0.25° spatial resolution and six-hour temporal resolution. The model was configured to use the Noah LSM to describe all land-atmosphere interactions; the Lin scheme to represent microphysical processes; the RRTM scheme to describe longwave radiation; the Dudhia scheme to represent shortwave radiation; the YSU scheme to describe the planetary boundary layer, and the Kain-Fritsch cumulus parameterization scheme for convective clouds.
Four experiments were carried out with WRF, each implementing a different albedo parameterization scheme (Table 2). The default Noah LSM albedo scheme [Eqs. (2)−(4)] was implemented in the control experiment (EXP1), with hourly model output. In the second experiment (EXP2), the improved albedo scheme [Eqs. (7)−(9)] was implemented in WRF + Noah LSM, using remote sensing products (MOD09A1, MOD09CMG, and MYD09CMG) and sparse in situ observations of SD to accurately estimate
$ \alpha_{\rm{firn}}$ ,$ \alpha_{\rm{freshsnow}}$ ,$ t^*$ ,$ d^*$ , and bare ground albedo in the albedo parameterization. Although the WRF model simulation has model errors including systematic deviation and mode integration error accumulation, the WRF + Noah LSM appears to be able to verify not only the spatial pattern but also the observed values of SD at the sparse meteorological stations, confirming the model’s performance of regional snow estimates (Liu et al., 2019). Therefore, to improve WRF model albedo estimates during the snow event, the modeled SD was considered to generate the improved albedo scheme, which was implemented in the third experiment (EXP3). The third experiment (EXP3) was similar to EXP2 but used hourly estimates of SD from EXP1. Model results for EXP2 and EXP3 were output at six-hour intervals. To evaluate the performance of the WRF + Noah LSM using the improved albedo scheme, i.e., EXP2 and EXP3, we also considered results from our previous numerical experiments with WRF + CLM, in which the more complex WRF + CLM albedo parameterization was used (Liu et al., 2019), referred to here as EXP4. The first day was used for model spin-up.Experiment Land surface physics Albedo parameterization Data used to estimate
parameters in improved
albedo schemeLand cover type Fractional vegetation
coverEXP1 Noah Default Noah none default default EXP2 Noah Improved albedo scheme Observed snow depth, MOD09CMG and MYD09CMG default default EXP3 Noah Improved albedo scheme Model snow depth and MOD09CMG default default EXP4 CLM Default CLM none default default EXP5 Noah Improved albedo scheme Model snow depth and MOD09CMG MCD12Q1 default EXP6 Noah Improved albedo scheme Model snow depth and MOD09CMG MCD12Q1 CR algorithm (Carlson and Ripley, 1997) EXP7 Noah Improved albedo scheme Model snow depth and MOD09CMG MCD12Q1 GI algorithm (Gutman and Ignatov, 1998) Table 2. Overview of the design of numerical experiments with the WRF.
Land cover and vegetation are important factors in the snow event simulation through their influences on the estimation of surface albedo and the interception of snowfall. In this study, the four experiments used the same default static land surface products to define land cover and green vegetation coverage. The default land cover dataset in the model was produced by the United States Geological Survey (USGS), using multispectral image data acquired by the Advanced Very High Resolution Radiometer (AVHRR) from April 1992 to March 1993, and adopted the 24 classification categories from the USGS. The default green vegetation coverage dataset in the model has a spatial resolution 0.144°, and was derived by Gutman and Ignatov (1998) using the AVHRR Normalized Difference Vegetation Index (NDVI) from 1985 to 1990.
-
The albedo parameterization scheme is expanded to represent albedo as a function of SD, snow age, snow-free albedo, fresh snow albedo, firn albedo [Eqs. (7)−(9)]. Snow depth (SD) and snow age were taken from ground observations and WRF outputs. Fresh snow albedo and snow-free albedo were taken from data retrieved from MODIS. Firn albedo and the scales for snow age and SD were estimated from nonlinear fitting. Based on the WRF + Noah LSM, the expanded albedo scheme was used in EXP2 and EXP3.
-
The MOD09A1 and MCD12Q1 products from 2017 were used to determine snow-free albedo and fresh snow albedo. More specifically, these data were used:
(1) Broadband albedo was calculated from the MOD09A1 dataset and Eq. (1) for each pixel.
(2) Pixels with MODIS retrieved albedo greater than 0.7 and lower than 0.3 were classed as fresh snow and snow free, respectively. Thus distributions of fresh snow and snow-free albedo were assigned to each specific land cover type, assuming land cover types as mapped in the MCD12Q1 product.
(3) The third quartile of the distribution of MODIS retrieved fresh snow albedo, selected as described in (2), was used as fresh snow albedo [
$\alpha_{{\rm{freshsnow}}}$ in Eq. (7)] for each specific land cover type. Similarly, the third quartile of the distributions of snow-free albedo found for each land cover type in (b) was used for land cover specific$\alpha_{{\rm{snowfree}}}$ , which we replaced with$\alpha_{{\rm{ice}}}$ in Eq. (8).Taking all land cover types together, the first and the third quartiles of the fresh snow albedo are 0.71 and 0.84, respectively (Fig. 2), and the first and third quartiles of the snow-free albedo are 0.02 and 0.27, respectively. For simplicity, the third quartile values for the two subsets of different land cover types were averaged to give the fresh snow albedo (0.79) and snow-free albedo (0.19), which were then used in the improved albedo parameterization when it was implemented in the Noah LSM scheme.
-
Surface spectral reflectance from the MOD09CMG product and hourly SD estimates from EXP1 were used to calculate the parameters needed for the improved albedo scheme [Eqs. (7) and (8)] in EXP3. The in-band reflectances from the MODIS product were combined to estimate the broadband albedo using Eq. (1). We combined the WRF SD values with the albedo values calculated from the retrievals corresponding to the time closest to the snowfall. The snowfall time was extracted from the hourly WRF output for EXP1. As SD increases, the albedo increases for shallow snow but remains constant for deep snow (where the albedo is already high). We define deep snow as snow where SD is at least 20 cm. Therefore, snow albedo is considered equal to the albedo retrieved from the MODIS data, where the SD is at least 20 cm. The MODIS albedo for locations where the SD values from WRF (EXP1) are at least 20 cm was used to determine and , for use in the improved albedo scheme in EXP3 [Eq. (7)].
To determine the relationship between the ground observations of SD and albedo, MODIS surface spectral reflectance products from the Terra (MOD09CMG) and Aqua (MYD09CMG) platforms were used in the improved albedo scheme for EXP2. Similar to EXP3, where ground observations of SD are at least 20 cm, the albedo from MODIS was used for
$\alpha_{\rm{firn}} $ and t* in the improved albedo scheme [Eq. (7)] for EXP2. -
Spatial variability in SD and fragmentation of snow cover immediately after snowfall are both high for this snow event. The modeled SD and MODIS retrieved albedo during the modeled snowfall are shown in Fig. 3. It is clear that albedo is less than 0.7 for a large range of SD estimates (values from WRF), with a wide range of values. For example, albedo corresponding to a SD greater than 100 cm varies between 0.1 and 0.75 at higher elevations in the Himalayas, which illustrates that the modeled snowfall is actually a snowmelt process in this complex terrain region (Fig. 3). In these conditions, the difference may be attributable to the difference in spatial resolution between the MODIS (500 m × 500 m) and WRF (25 km × 25 km) data, from which the albedo and SD were taken. The different resolutions mean that the highly variable surface types and elevation are captured to different degrees in the two datasets. It is concluded that the coarse grid resolution WRF generally overestimates SD for the high elevations in the Himalayas. We, therefore, used SD values from WRF only where the estimates were less than 100 cm to determine d* for use in the improved albedo scheme [Eq. (8)] in EXP3. We considered snow albedo equal to the fresh snow albedo during the snowfall event. The ground observed SD values were used to estimate d* in the improved albedo scheme [Eq. (8)] in EXP2 during snowfall.
Figure 3. (a) Scatterplot of snow depth (SD) from EXP1 estimates and Terra MODIS albedo and (b) spatial distribution of SD estimates from EXP1 during snowfall.
The steps to estimate the albedo-dependent parameters (i.e.,
$ \alpha_{\rm{freshsnow}}$ ,$ \alpha_{\rm{snowfree}}$ ,$ \alpha_{\rm{firn}}$ , t*, and d*) in the improved albedo parameterization scheme are shown in Fig. 4. The values used for these parameters in EXP2 and EXP3 are listed in Table 3. Our value for firn albedo in EXP2 is 0.51, which is similar to values from other studies, e.g., 0.53 in Oerlemans and Knap (1998) and 0.5 in Yang et al. (2013), but is different than that used in EXP3 (Table 3).Figure 4. Flowchart illustrating the steps to estimate the parameters for the improved albedo parameterization scheme in EXP2 and EXP3 using MODIS albedo products and snow depth (SD).
Experiments Fresh snow albedo Snow-free albedo Firn albedo d* (m) t* (d) EXP2 0.79 0.19 0.51 0.023 0.42 EXP3 0.79 0.19 0.13 0.11 1.38 Table 3. Optimal values for parameters in the improved albedo parameterization scheme.
-
To assess the performance of the WRF when applying various albedo parameterization schemes, we compared field observations of near-surface air temperatures and SWE from 32 CMA stations, albedo values from 6 CAS stations, and MODIS retrieval products with the model’s estimated values. At local solar noon in Lhasa (1400 LST, LST = UTC + 8 h), the observed albedo value is closer to the Lambertian albedo described by the WRF model when coupled with LSMs. Thus, we used albedo observations at 1400 LST to evaluate the model estimated albedo. The root mean square error (RMSE), mean absolute deviation (MAD), temporal correlation coefficient (CC), spatial correlation coefficient (SCC), and relative difference (
$\delta $ ) were used to evaluate the model’s performance [Eqs. (10)−(13)]. The calculation of the SCC follows that of the CC, but areas are weighted by the sine of the latitude. Weighting is also performed appropriately for unequally spaced grids. The SCC was calculated in the region from 65°E to 106.8°E, and from 25°N to 40.8°N, for 10752 grids in total. Specifically:where N is the total number of observed or simulated data (N = 10 when evaluating the performance of albedo estimates at one station, and N = 40 when evaluating the performance of near-surface air temperature estimates at one station);
$X_{\rm{p}}^{(j)}$ and$X_{\rm{o}}^{(j)}$ are the modeled and observed values at timestep j, respectively;${\overline {X_{\rm{p}}}}$ and${\overline {X_{\rm{o}}}}$ are the mean of the modeled and observed values, respectively;${\delta ^{\left( j \right)}} $ is the relative difference at timestep j in units of percentage.
Station Type | Station number | Name | ID | Latitude (°N) | Longitude (°E) | Elevation (m) |
CMA | 1 | Mangya | 51886 | 38.25 | 90.85 | 2945.0 |
2 | Daqaidam | 52713 | 37.85 | 95.37 | 3174.0 | |
3 | Gangca | 52754 | 37.33 | 100.13 | 3302.0 | |
4 | Golmud | 52818 | 36.42 | 94.90 | 2807.6 | |
5 | Doulan | 52836 | 36.30 | 98.10 | 3189.0 | |
6 | Xining | 52866 | 36.62 | 101.77 | 2295.2 | |
7 | Shiquanhe | 55228 | 32.50 | 80.08 | 4278.6 | |
8 | Baingoin | 55279 | 31.37 | 90.02 | 4700.0 | |
9 | Nagqu | 55299 | 31.48 | 92.07 | 4507.0 | |
10 | Xainza | 55472 | 30.95 | 88.63 | 4672.0 | |
11 | Xigaze | 55578 | 29.25 | 88.88 | 3836.0 | |
12 | Lhasa | 55591 | 29.67 | 91.13 | 3648.9 | |
13 | Tingri | 55664 | 28.63 | 87.08 | 4300.0 | |
14 | Lhunze | 55696 | 28.42 | 92.47 | 3860.0 | |
15 | Pagri | 55773 | 27.73 | 89.08 | 4300.0 | |
16 | Tuotuohe | 56004 | 34.22 | 92.43 | 4533.1 | |
17 | Zadoi | 56018 | 32.90 | 95.30 | 4066.4 | |
18 | Qumarleb | 56021 | 34.13 | 95.78 | 4175.0 | |
19 | Yushu | 56029 | 33.00 | 96.97 | 3716.9 | |
20 | Madoi | 56033 | 34.92 | 98.22 | 4272.3 | |
21 | Darlag | 56046 | 33.75 | 99.65 | 3967.5 | |
22 | Zoige | 56079 | 33.58 | 102.97 | 3441.4 | |
23 | sog | 56106 | 31.88 | 93.78 | 4024.0 | |
24 | Dengqen | 56116 | 31.42 | 95.60 | 3873.1 | |
25 | Qamdo | 56137 | 31.15 | 97.17 | 3315.0 | |
26 | Sertar | 56152 | 32.28 | 100.33 | 3894.0 | |
27 | Barkam | 56172 | 31.90 | 102.23 | 2664.4 | |
28 | Songpan | 56182 | 32.67 | 103.60 | 2850.7 | |
29 | Batang | 56247 | 30.00 | 99.10 | 2589.2 | |
30 | Nyingchi | 56312 | 29.57 | 94.47 | 2991.8 | |
31 | Deqen | 56444 | 28.45 | 98.88 | 3319.0 | |
32 | Jiulong | 56462 | 29.00 | 101.50 | 2925.0 | |
CAS | 33 | QOMS | 28.36 | 86.95 | 4276.0 | |
34 | SETS | 29.77 | 94.73 | 3326.0 | ||
35 | NASDE | 33.39 | 79.70 | 4264.0 | ||
36 | MASWE | 38.41 | 75.05 | 3668.0 | ||
37 | SHSEX | 33.22 | 88.83 | 4947.0 | ||
38 | MAQU | 33.92 | 102.10 | 3440.0 | ||
Note: Qomolangma station (QOMS), Southeast Tibet station (SETS), Ngari station (NASDE), Muztagh Ata station (MASWE), ShuangHu station (SHSEX), MAQU station (MAQU) |