Evaluation of Ten Fresh Snow Density Parameterization Schemes for Simulating Snow Depth and Surface Energy Fluxes on the Eastern Tibetan Plateau

Abstract

Snow cover on the Tibetan Plateau has a shallow depth, plaque distribution, and repeated ablation. The applicability of the snow parameterization scheme in the current land surface process model on the TP needs to be further tested using observational data. In this paper, using the land surface process model CLM4.5 and ten fresh snow density parameterization schemes characterized by temperature, wind speed, and relative humidity, three discontinuous snow processes in Maqu, Madoi, and Yakou and two continuous snow processes in Madoi and Yakou were simulated. By comparing the simulated snow depth with the observed, it was found that this model can clearly describe repeated snow accumulation and ablation processes for the discontinuous snow cover process. The K_W scheme, compared with the original Anderson scheme, performed the best regarding snow depth simulation. However, all schemes overestimated the melting rate of snow, and were not able to simulate continuous snow accumulation. The simulation effect of the Schmucki scheme on radiation and energy flux under discontinuous snow cover was significantly improved compared with other scheme. None of schemes performed perfectly, so future studies that focus on the simulations of snow depth, radiation flux, and energy flux under continuous snow cover for accurate and wide applications are recommended.

1. Introduction

Snow cover has been identified as the most variable parameter of land surface condition in both spatial and temporal areas of a region locally via the surface energy balance budget, or remotely via snow–atmosphere coupling and the subsequent implications for circulation [1,2,3]. It is considered to play a significant role in the climate, including its high surface albedo, low thermal conductivity, and the energy required to melt the pack during the snowmelt season [1,4].

The Tibetan Plateau (TP) has an average elevation of approximately 4000 m and is a major snow-covered area in China and the northern hemisphere [5,6,7]. The changes in snow cover on the TP affect the interaction between the TP surface and atmosphere in terms of its heat and dynamics, eventually affecting the regional climate [8,9]. Many studies have shown that snow during the winter and spring on the TP is an important predictor of precipitation in Asia [10,11]. Excessive snow weakens the surface heating of the plateau and the Asian summer monsoon, leading to less summer precipitation in India and Southern and Northern China and more summer rainfall in the middle and lower reaches of the Yangtze River [8,12,13]. Some scholars have also believed that melted snow penetrates the soil and extends the climatic effect of increasing snow [14,15]. In addition, the energy gained on the ground increased as the snow area decreases, so one of the main reasons for global warming in spring is the cooling effect of snow decreases in the northern half, as the cryosphere has shrunk in recent years [16,17], and atmospheric warming could lead to glaciers melting and even ice avalanches, causing a disaster [18]. Otherwise, snow cover also affects hydrological processes, vegetation phenology, and the carbon cycle [19,20,21].

The snow cover process on the TP is different compared with other high latitudes in the northern hemisphere. Except for some mountainous areas, the snow cover in most areas of the TP cannot last for the whole winter. This characteristic is mainly manifested, as it cannot accumulate in the autumn and continue into the spring of the following year. Accumulation and ablation processes occur in cycles for one year, and the coupling between snow cover and the atmosphere is very complex. Observational data have shown that the spatial distribution of snow cover on the TP is extremely uneven [22], and has a shallow depth, plaque distribution, and repeated ablation [23]. Remote sensing and observational data have shown that snow is mainly distributed in the Himalayas, the western side of the Hengduan Mountains, the Nyenchen Tanglha Mountains, and the Bajankla Mountains, while the depth of snow in the hinterland of the TP is relatively shallow [24,25,26]. According to observations from meteorological stations over the last 50 years, 70% of the stations’ annual average snow cover days were less than 30 d, and 77% of the stations’ mean daily snow depth of approximately was less than 0.5 cm, with the same spatial distribution [25,27], with the majority of them being discontinuous snow cover process. In discontinuous or sporadic permafrost areas, snow changes rapidly at an intraseasonal time scale, and the variability is relatively high in the central and eastern regions [8,28]. Due to this unique snow distribution and change characteristics (rapidly changing) on the TP, most snow processes begin with no snow cover, which means the snowfall depth is equal to the snow depth, so fresh snow density and its influenced elements are to become the focus of this research.

Fresh snow density is a critical parameter for simulating snow depth, which influences albedo and surface energy fluxes in land surface processes and hydrological models. As fresh snow density varies widely, roughly 0.1~0.35 g/cm³ [29], it is challenging to parameterize it in the land surface processes models. A study has shown that fresh snow density is easily influenced by the temperature, relative humidity, wind speed, solar radiation, and other meteorological factors around the world [30]. Higher surface temperatures cause snow particles to melt and metamorphose. At higher wind speeds (of more than 9 m/s), ice crystals can move on the surface, which causes surface compaction. Humidity is one of the necessary conditions for the growth of snow particles, and solar radiation also is an influential factor for fresh snow density. Therefore, snow compaction, crystal deterioration, and ice water content porosity caused by changed meteorological factors will further lead to the change in snow density within a certain range. At present, the widely used fresh snow density parameterization schemes in land surface processes and hydrological models include the following: (1) The constant method; for example, 0.08 g/cm³ (temperature is less than −15 $\mathrm{℃}$) or 0.05 cm⁻³ (temperature is greater than −15 $\mathrm{℃}$) [31]; (2) the empirical formula method; for example, choosing the influence of temperature and other meteorological factors forms the empirical formula, such as an exponential function [32,33,34,35,36]; and (3) the post-treatment diagnostic method; the fresh snow density was diagnosed under two conditions according to the temperature vertical profile [37]. The constant method and empirical formula method are used to determine the fresh snow density to simulate snow depth and surface energy budget for further exploring the changing of the land–atmosphere process after snow cover in the land surface processes models.

Except in the summer, the snow density on the TP is the lowest among the three stable snow areas in China, at less than 0.14 g/cm³ [38]. This value is significantly different from that obtained for North America (0.5 g/cm³, sometimes even higher) [39]. At present, fresh snow density parameterization schemes are mainly developed for snow cover characteristics in North America. Due to the snow on the TP being dry and cold, studies have revealed that the simulation snow depths of a single point and region are all too big [40,41]. The inaccuracy of the new snow density is related to the variance of the snow depth simulation on the TP. To improve the simulation of the land surface process model in the snow cover area on the TP, the applicability of various parameterization schemes on the TP needs to be further tested in detail.

This work addressed the above issues by comparing ten fresh snow parameterization schemes using the Community Land Model 4.5 (CLM4.5) with three observation sites on the TP. Our paper was organized as follows: Section 2 described the datasets, model, and the 10 fresh snow parameterization schemes in this paper; the results of the simulated snow depths via CLM4.5 and the 10 schemes, and fluxes with the K_W and Schmucki schemes were presented in Section 3; the discussions and conclusions were presented in Section 4.

2. Data and Methods

2.1. Study Area

This study used three field observation stations located in different climatic regions in the eastern part of the plateau. As shown in Figure 1, the three field observation stations include Maqu, at an elevation of 3423 m in the cold-temperate semi-humid zone, Madoi, at an elevation of 4280 m in the cold and semiarid zones, and Yakou, at an elevation of 4148 m in the arid sub-frigid zone. The vegetation of all three stations consists of alpine meadows [42,43,44], which are located in seasonally frozen soil areas and seasonal snow-covered areas [45].

2.2. Data

Maqu and Madoi are two observation points located at the Zoige Ecosystem Research Station of the Chinese Academy of Sciences. Maqu station is located at the Hequ Horse Farm in Maqu County, Gansu Province, and Madoi station is located at the Ngoring Lake in Madoi County, Qinghai Province. These observation data all have strict quality control, and they include meteorological elements, soil temperature, soil humidity, and eddy flux. The observation data from Yakou come from the Heihe River Basin Surface Process Observation Network of the Chinese Academy of Sciences, Dadongshu-Yakou, Qilian County, Qinghai Province. They include two automatic meteorological observation points (10 m apart). The time steps of the meteorological elements and soil temperatures are 30 min and 10 min, respectively; the time step of the eddy data is 30 min. These data were downloaded from the National Tibetan Plateau Science Data Center (https://data.tpdc.ac.cn/en/, accessed on 24 March 2022). The time resolution of the snow depth measured via ultrasonic ranging was one day. The observation items and instrument heights at the three stations are shown in the Supplementary Materials [46,47,48,49], detailed in Table S1.

The daily snow depth data at the Maqu and Madoi stations were obtained from the TP snow dataset from CMA (http://data.cma.cn, accessed on 24 March 2022). The geography range was 25–40° N, 75–105° E, and the time range was from 1951 to 2018.

2.3. Methods

2.3.1. The Representatives of Stations and Snow Processes

According to previous research results [25,26,27], over the last 50 years, meteorological observation stations have shown that the annual average snow cover days on the TP are relatively short, with less than 25 d, and the largest season is during the winter, with approximately 10 d; the mean daily snow depth is relatively shallow, and does not exceed 0.5 cm, mostly being discontinuous. Compared with the TP from the past 50 years, Maqu exhibits similar snow cover days and snow depths, while Madoi shows slightly larger values calculated for the period of 2014–2017 [50]. Yakou recorded an annual average of 211 snow cover days, which aligns with previous research results [51]. The mean daily snow cover depth was 4.05 cm, which is more consistent with the distribution characteristics of snow cover in some mountainous areas on the plateau [52,53]. Therefore, the selected stations are representative in their respective contexts. Detailed data are presented in Table S2.

Due to changes in the sun’s altitude and air humidity within a day, the daily variation of the surface albedo exhibits a U-shape. Surface albedo increases significantly after a snowfall, with fresh snow cover ranging from 0.8 to 0.9. As the snow ages, it becomes uneven and wet, and albedo gradually decreases. Generally speaking, dense, dry, and clean snow has an albedo of 0.9; damp snow ranges from 0.5 to 0.6; porous dirty snow has an albedo of 0.3–0.4 [54]. The daily average of albedo in the snow-free period of alpine grassland is around 0.2–0.27 [55]. As shown in Figure 2, the mean daily surface albedo at Maqu from 2014 to 2015 was 0.26, and it increased significantly to about 0.6 after snowfalls (snow depth missed in December 2014–February 2015 and April 2015). The albedo characteristics of the other two stations were the same as that detailed above. Therefore, the snow cover period was defined when the daily average albedo is larger than 0.4 and the daily snow depth is more than 0 cm.

2.3.2. Model and Setup

The CLM is a series of land models developed through the Community Earth System Model (CESM) project. As a relatively well-developed land surface model, which includes biogeophysical and biogeochemical processes, it has been widely used and developed around the world. With the update of the version, this model has been improved and perfected in all aspects, and parameterization relevant to snow cover is also being improved. CLM4 included more sophisticated representations of soil hydrology and snow processes based on its previous version. In particular, new treatments of soil column–groundwater interactions, soil evaporation, vertical burial of vegetation by snow, snow cover fraction and aging, and vertical distribution of solar energy for snow were implemented, and a new snow cover fraction parameterization was incorporated that reflects the hysteresis in fractional snow cover for a given snow depth between the accumulation and melt phases [56]. Later, CLM4.5 uses a nested sub-grid hierarchy, in which each grid cell is composed of multiple land units, snow or soil columns, and plant functional types [57]. Snowpack evolution, such as snow accumulation, depletion, densification, metamorphism, and percolation, along with the refreezing of water were simulated with up to five unevenly spaced vertical layers, and the number of vertical layers was determined based on the snow depth [58]. The thickness of each layer was prescribed using predetermined minimum and maximum thickness values that were described by Jordan [59], which were determined in a Eulerian fashion. When the snow depth is less than a specified threshold value (0.01 m), snow exists in the model without being represented by explicit snow layers [57], whereas when it is greater than or equal to 0.01 m, a snow layer is initialized. This model handles combinations and subdivisions by adjusting the snow depth in each layer and the number of snow layers, depending on whether there is an accumulation or a depletion of snow.

In CLM4.5, the parameterizations for snow are primarily based on Anderson [60], Jordan [59], and Dai and Zeng [61]. Among them, fresh snow density is calculated separately for snow accumulation and depletion. When snowfall occurs, fresh snow density is calculated following the Anderson scheme (the detailed formula will be mentioned in Section 3.1.1).

The annual snow cover days and mean daily snow depth at Maqu are 24.75 d and 0.18 cm, respectively. Therefore, the Maqu station only used a discontinuous snow process in the winter for simulation. The simulation period (Period A) spanned from 12 December 2014 to 20 January 2015 (40 d), with a time step of 30 min in the model. For this experiment, the model’s surface parameters and initial soil temperature and humidity were set as follows: an alpine meadow land cover type, with a vegetation coverage of 100% and a canopy height of 0.2 m. Since the soil was frozen when the simulation was started in December, it is crucial to enhance the stability and accuracy of the model’s results by adding the initial ice content into the soil (Table S3). In Table S3, the initial ice content was calculated for each soil layer by converting the ice volume and water volume based on differences between the current soil water content and the average liquid water content during the summer months, along with other data obtained from field observations and sampling experiments.

The annual average number of snow days at Madoi from 2014 to 2017 was 49.25 d, and the mean daily snow depth was 0.24 cm. Compared with Maqu, there were more snow days, and Madoi had a snowy year in 2018. Therefore, a continuous snow process over a long period of time and a discontinuous snow process over a short period of time were selected for our simulation tests. The two time periods were as follows:

(1)

From 1 November 2014 to 20 January 2015, 81 d in total (Period B), among which the spin-up period was 39 d, from 1 November to 9 December, and atmospheric forcing was also updated every 30 min.

(2)

From 21 October 2018 to 17 April 2019, a total of 179 d (Period D). The initial soil temperature and humidity settings were added to make the model run stably in this period, due to data limitations caused by spin up (Table S4). Since the simulation starts in October, the initial ice content of the soil is zero.

The surface parameters and soil mechanical composition were set under the grassland land cover type, with a vegetation coverage of 55% and a canopy height of 0.05 m.

According to the statistical results of Yakou, with 211 d and 4.05 cm from 2014 to 2017, a relatively thick level of snow existed for nearly 2/3 of the year from the autumn of the previous year to the spring or even the summer of the following year. Therefore, two snow processes were also selected for the simulation test:

(1)

Another discontinuous snow process in a short period of time, from 11 August 2015 to 30 September 2015, with a spin-up time period from 1 August to 10 August (Period C).

(2)

One continuous snow process similar to that at Yakou, from 1 October 2014 to 20 March 2015, with a spin-up time period from 1 August to 30 September (Period E).

As the snow accumulation period of Yakou has a spin-up time, we added the soil mechanical composition data gained via soil sample collection in the field, pretreatment (weighted, dyed, and sift), and detection (Table S5). In addition, the ground cover type was also the alpine meadow, with a vegetation coverage of 70% and a canopy height of 0.02 m.

2.3.3. Evaluation Method of Simulation Effects

The evaluation of model simulation values was mainly based on the following three statistical methods:

(1)

The correlation coefficient (R): indicates the degree of similarity between the simulated value and the observed value change trend.

(2)

Mean deviation (ME): represents the size of the overall deviation between the simulated value and the observed value.

(3)

Root mean square error (RMSE): represents the simulated value. The magnitude of the deviation between the value and the observed value is the superposition of the simulation effect of the simulation value at each moment in the entire simulation period.

3. Results

3.1. Simulation of Snow Depth in Schemes

3.1.1. Air Temperature Schemes

This part mainly used four fresh snow density parameterization schemes that were only related to the air temperature: the Anderson (1976) scheme (CLM4.5 origin scheme), the Kampenhout (2017) scheme (K_T), the Pomeroy (1998) scheme, and the test scheme designed by the author; specific calculations have been detailed in Formulas (4)–(7). The test scheme chose the formula of maximum fresh snow density to be calculated for each temperature interval from the other three air temperature schemes, so that it can reflect the effects of density and snow depth.

The calculated values of fresh snow density for each scheme are shown in Figure 3. The range of calculated density values using the Anderson and K_T schemes was 50~169 kg/m³. The Pomeroy scheme had a relatively high density when the temperature was lower than −10 $\mathrm{℃}$, but had the lowest values when the temperature was between −10 $\mathrm{℃}$ and 2.59 $\mathrm{℃}$. The test scheme had the highest density in all temperature ranges, above 70 kg/m³.

Anderson (1976):

$${\rho }_s=\left\{\begin{array}{c}{50+1.7\left(17\right)}^{1.5} T-T_{frz}>2\\ 50+1.7(T-T_{frz}+15{)}^{1.5} -15<T-T_{frz}\le 2\\ 50 T-T_{frz}\le -15\end{array}\right.$$

(1)

Kampenhout (2017):

$${\rho }_s=\left\{\begin{array}{c}{50+1.7\left(17\right)}^{1.5} T-T_{frz}>2\\ 50+1.7(T-T_{frz}+15{)}^{1.5} -15<T-T_{frz}\le 2\\ -3.8328\left(T-T_{frz}\right)-0.0333(T-T_{frz}{)}^2 T-T_{frz}\le -15\end{array}\right.$$

(2)

Pomeroy (1998):

$${\rho }_s=\left\{\begin{array}{c}67.92+51.52e^{T-T_{frz}/2.59} T-T_{frz}\le 0 ℃\\ 119.2+20\left(T-T_{frz}\right) T-T_{frz}>0 ℃\end{array}\right.$$

(3)

Test:

$${\rho }_s=\left\{\begin{array}{c}119.2+20(T-T_{frz}) T-T_{frz}>2\\ 50+1.7(T-T_{frz}+15{)}^{1.5} -10<T-T_{frz}\le 2\\ 67.92+51.52e^{T-T_{frz}/2.59} -20<T-T_{frz}\le -10\\ -3.8328\left(T-T_{frz}\right)-0.0333(T-T_{frz}{)}^2 T-T_{frz}\le -20\end{array}\right.$$

(4)

Among them, T denotes the air temperature, with its unit being K, and T_frz is the ice–water phase transition temperature, with a value of 273.15 K.

The snow depth simulated using the four schemes solely based on temperature in the five simulation periods is illustrated in Figure 4. The black dotted line represents the observed snow depth. All four schemes were able to simulate three distinct snow processes during Period A when compared to the observations. However, the results show some deviation between the observations and simulations in snow depth, snowfall days, and end time. In detail, the simulated snowfall day occurred one day later than the observed, while it ended one day earlier, and the duration of continuous snow lasted for only one day from 12 to 14 December. Among these results, the Anderson scheme was the largest (9 cm), the Pomeroy scheme was the smallest (6 cm), and the other two schemes formed results of 7 cm, which all exceeded the observed results.

For the snowfall on 18 December, all schemes exhibited identical values for both the number of snowfall days and snow depth; however, ablation occurred two days later than observed. During the snow cover in January, although the simulated number of snowy days aligned with observations’, ablation took place two days earlier than what was recorded. Overall, the variation tendency of the snow processes in the simulated and observed results were basically consistent, and the maximum snow depth days were the same (6 January). On that day, the simulated snow depth of the Anderson, KT, Pomeroy, and test schemes were 15 cm, 15 cm, 19 cm, and 14 cm, respectively, and the value of the test scheme was the closest to the observed (11 cm).

Figure 4b shows that the four schemes can simulate the discontinuous snow process in Period B basically; however, there was a clear distinction of the four schemes when compared to observed snow depth, which may be caused by the precipitation forcing variable inputs from the Madoi station over time. However, under this hypothetical precipitation condition, the test scheme outperforms.

For period C of the Yakou station, the results of each scheme in Figure 4c show slight differences. This is due to the daily average air temperature being 1.94 $\mathrm{℃}$ during the autumn, accompanied with the daily average air temperature on snowfall days being above −1 $\mathrm{℃}$. When the temperature exceeded 0 $\mathrm{℃}$, the calculated fresh snow densities were the same in the Anderson, K_T, and test schemes, and the density of the Pomeroy scheme within this temperature range was slightly different from the other schemes (Figure 3). Two snowfalls were predicted in August, but the lack of observed data in this period may be attributed to wind migration or melting caused by a higher surface temperature. In the subsequent three snow processes, the simulation performance was superior compared to the other periods.

The following are the simulation results of the long-term continuous snow processes at Madoi station (Period D) and Yakou station (Period E) using the four abovementioned schemes (Figure 4d,e). The snow process accumulates from the fall of the previous year to the spring of the following year.

The simulated snow depth in Figure 4d exhibits a relatively minor deviation compared to the observed values prior to November. Subsequently, both the observed and simulated snow depths demonstrate a consistent trend of accumulation and melting characteristics. Notably, with the decreasing trend, the simulated values quickly decreased to zero. This difference was noticeable, and the simulated value was generally larger in the accumulation period, and the longer the snow accumulated, the more obvious this phenomenon became. As shown in Figure 4e, the overall simulation trend at Yakou station is similar to the previous process at Madoi. Therefore, a common issue in simulating temporal variations of continuous snow covers lies in accurately describing prolonged snow accumulation; the model significantly overestimated the rate of snowmelt, and this was clearly related to the unique characteristic of the snow cover on the plateau. However, among all four schemes tested for simulating the snow depth magnitude, it is worth mentioning that the test scheme remained superior.

Combining the calculation results of the statistics in Figure 5, the simulated values of these four schemes have a good linear correlation with the observed in Period A (Figure S1a1–a4), all being above 0.85. From the ME and RMSE, the deviation of the Pomeroy scheme was the largest and the test scheme was smallest (0.9 cm). In Period B, the simulated value was approximately five to thirty-five times as much as the observed; the linear correlation coefficient was less than 0.5, and the ME and RMSE were 4 cm and 9 cm, respectively. In Periods A–B, data inspection results indicated that the test scheme was the best.

According to the consistency of the ME and RMSE, the calculation methods have slight differences among the four schemes in Period C.

The correlations in Periods D–E, compared with the discontinuous snow cover process (A–C), were very poor, when all the values were below 0.2. In Period D, the high-value areas of the four schemes were distributed near 75 cm, while that of the observed were all below 40 cm, indicating that the simulated snow depth is generally overestimated. Since none of the schemes simulate continuous accumulation, the large simulated values throughout offset the zero values during complete ablation, resulting in a negative mean error for the total sample. In period E, the MEs were approximately 5 cm. Based on the RMSE, the minimum value of the test scheme still shows its optimal effect.

To sum up, the four schemes in this section within the five time periods have better simulation effects on the discontinuous snow cover process at Maqu (Period A) and the continuous snow cover process at Madoi (Period D). Compared with the various schemes, the test scheme has the best effect on the snow depth in each period.

3.1.2. Wind Speed and Air Temperature Schemes

In this part, four schemes related to air temperature and wind speed were used to conduct simulations in every period. The four schemes used were the Liston (2007) scheme, the Jordan (1999) scheme, the Kampenhout (2017) scheme (hereinafter referred to as K_W), and the Vionnet (2012) scheme (specific calculations have been detailed in Formulas (8)–(11).

Among them, the Liston scheme adds wind speed items, while the wind speed is greater than 5 m/s at a height of 2 m. The K_W scheme adds a smooth calculation term at 5 m/s based on the K_T scheme. The air temperature input variable in the Jordan scheme is Kelvin, and 260.15 K is the critical temperature of the two calculation methods. The Vionnet scheme only uses one calculation method under all conditions. Except for the Liston scheme, the height of the input instantaneous wind speed of all other schemes mentioned is 10 m. Therefore, the wind speed has already corrected the required height (Formula (11)) for more accurate simulation effects.

Figure 6 compares the fresh snow density calculated for the four schemes in this section when the temperature was between −20 and 1 $\mathrm{℃}$, and when the wind speed was between 1 and 15 m/s. As shown in the Liston scheme in Figure 6a, when the wind speed was less than 5 m/s, the density of fresh snow was less than 150 kg/m³; when the temperature was lower than −15 $\mathrm{℃}$, the density changed inconspicuously; when the temperature was higher than −15 $\mathrm{℃}$, the density gradually increased to 350 kg/m³. The fresh snow density calculated using the Vionnet scheme had no obvious boundary between the air temperature and wind speed, and the maximum density was approximately 200 kg/m³. Relative to the Liston scheme, the K_W scheme had a smooth calculation at a wind speed of 5 m/s. The interval of fresh snow densities greater than 350 kg/m³ was larger than that of the Liston scheme. The Jordan scheme was still able to reach a maximum of 350 kg/m³, but its interval was smaller than that of the Liston and K_W schemes.

Liston (2007): $T_{wb}\ge 258.16$

$${\rho }_s={\rho }_T+{\rho }_{ws}$$

$${\rho }_T={50+1.7\left(T_{wb}-258.16\right)}^{1.5}$$

$${\rho }_{ws}=25+250\left\{1.0-exp\right.\left[-0.2(W_S-5.0)\left]\right\} W_S\left(2m\right)>5 \mathrm{m}/\mathrm{s}\right.$$

(5)

Jordan (1999):

$${\rho }_s=500[1-0.951\mathit{exp}\left(-1.4\left(278.15-T_a{)}^{-1.15}-0.008{W_s}^{1.7}\right)\right] 260.15 \mathrm{K}<T\le 275.65 \mathrm{K}$$

$${\rho }_s=500[1-0.904\mathit{exp}(-0.008{W_s}^{1.7}\left)\right] T\le 260.15 \mathrm{K}$$

(6)

Kampenhout (2017):

$${\rho }_s={\rho }_T+{\rho }_{ws}$$

$${\rho }_T=\left\{\begin{array}{c}{50+1.7\left(17\right)}^{1.5} T-T_{frz}>2\\ 50+1.7(T-T_{frz}+15{)}^{1.5} -15<T-T_{frz}\le 2\\ -3.8328\left(T-T_{frz}\right)-0.0333(T-T_{frz}{)}^2 T-T_{frz}\le -15\end{array}\right.$$

$${\rho }_w=266.861\left(0.5\right(1+\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}(W_s/5)){)}^{8.8}$$

(7)

Vionnet (2012):

$${\rho }_s=109+6({T-T_{frz})+26(W_s)}^{0.5}$$

(8)

$$V(z2)=V(z1){\left(\frac{z2}{z1}\right)}^{0.14}$$

(9)

In Period A, all schemes were able to effectively simulate the occurrence of three discontinuous snow processes. However, due to low wind speeds during the two snow processes in December, each snow depth of the simulation was approached. The wind speed was 5.73 m/s, and the maximum instantaneous value reached 9.67 m/s on 5 January; the simulation effects of the four schemes were significantly different. Although the K_W scheme exhibited a one-day lag compared to the peak snow depth simulated via the K_T scheme, the peak value was significantly reduced to 9 cm, which is closer to the observed value than that of the K_T scheme (15 cm). The Liston and Vionnet schemes showed thicker snow depths than the K_W and K_T schemes and 1~2 day earlier melting days than the observations. The snowfall day simulated via the Jordan scheme lags one day behind the observations and other schemes, and has an unclear description of the snowfall, accumulation, and gradual melting process.

Figure 7b shows the result of snow depth in Period B. For the snow process starting from 5 January, both the Vionnet and Jordan schemes indicate higher snow cover days compared to the observations. Regarding snow depth, the K_W scheme was still more suitable for the discontinuous snow cover process.

From the results of the snow depth in Periods B and C (Figure 7b,c), the simulation effects of the four schemes are better than those in Period B. The main reason for this is that the precipitation (input) and snow depth (output) were both collected from the Yakou field observation station, and the matching degree of the two series of data was relatively high. The simulated snow depth on the day of each scheme was larger than the observed depth on the snow cover process at the end of August, with the largest value formed via the Vionnet scheme (8 cm). In September, there were four discontinuous snow processes with short intervals. In each phase, the K_W scheme value was the closest to the actual, and the rest were larger than the observed values to varying degrees.

During the simulations of Periods D and E (Figure 7d,e), the parameterization schemes, after adding wind speed, displayed a greater improvement in the snow depth simulation than the previous temperature schemes. However, the evident character of faster snow melting did not ameliorate in this part.

Figure S2 shows a scatter plot comparison of the simulated and observed values of the four schemes related to air temperature and wind speed in the five simulation periods. In Period A, combined with the calculation results of the statistics in Figure 8, the R of the four schemes were more than 0.5 at 0.77, 0.58, 0.81, and 0.70, respectively, which were deemed to be similar to the temperature-related schemes. From the ME results, except for the Vionnet scheme, which was larger than 1 cm, the K_W, Jordan, and Liston schemes were less than 1 cm, which were smaller than the MEs of the temperature schemes. In summary, the K_W and Liston schemes have better simulation effects on snow depth during this period. Regarding Period B, the simulation effect of the K_W scheme was the best, having the smallest ME and RMSE.

The linear correlation coefficients (R) of the four schemes were reached in Period C. Combined with the characteristics shown in Figure S2, and the values computed in Figure 8, the K_W scheme was still the most suitable for simulating the discontinuous snow process.

In Period D, the R values were poor, at 0.17, 0.15, 0.16, and 0.20, respectively. The MEs of the K_W, Jordan, and Liston schemes were all negative, while the Vionnet scheme was positive. Considering the small ME, due to the simulated no-snow period, the largest absolute value of the K_W scheme should represent a relatively well-simulated effect. Thus, the largest ME of the Vionnet scheme was consistent with the discontinuous snow process, which all show unsatisfactory results. Although the snow depth and precipitation data have a high degree of matching at Yakou, the parameterization schemes with added wind speed terms are still not ideal for the simulation effect of snow depth in Period E. Therefore, the simulation of the snow depth for the continuous snow accumulation process still needs further improvement.

From this point of view, the four schemes in this section also have a better simulation effect for the discontinuous snow cover process at Maqu station (Period A), and the continuous snow cover process at Madoi station (Period D). In each simulation period, the Jordan scheme simulation was smaller than the observed value, while the Vionnet scheme was larger, and the K_W scheme was found to be the best scheme for the TP’s snow cover.

3.1.3. Relativity Humidity, Wind Speed, and Air Temperature Schemes

Two parameterization schemes of fresh snow density were used in this section for comparative research. These schemes add the calculation of relative humidity based on air temperature and wind speed; that is, fresh snow density is a function of air temperature, wind speed, and relative humidity. The calculation methods of the Lehning (2002) scheme and the Schmucki (2014) scheme are shown in Formulas (10) and (11), respectively. The conditions of the Lehning scheme are that the relative humidity is greater than 70%, and that the snow depth must be output in a time step of 30 min or 60 min. Therefore, the model output step of the snow depth was set to 30 min first; then, the daily snow depth was processed and compared with the observed data. The conditions for the Schmucki scheme for calculation are wind speeds greater than 2 m/s at a 10 m height and a relative humidity greater than 80%. Two methods were used to calculate the critical temperature in the Schmucki scheme, which was found to be −14 $\mathrm{℃}$.

Lehning (2002):

$${\rho }_s=70+6.5T+7.5T_s+0.26RH+13W_s-4.5T{T}_s-0.65T{W}_s-0.17RH{W}_s+0.06T{T}_{s}RH$$

(10)

Schmucki (2014):

$${\rho }_s={10}^c$$

$$\mathrm{c}=3.28+0.03T-0.36-0.75\arcsin \left(\sqrt{RH}\right)+0.3{log}_{10}\left(W_s\right)$$

$$\mathrm{c}=3.28+0.03T-0.75\arcsin \left(\sqrt{RH}\right)+0.3{log}_{10}\left(W_s\right) T<-14 \mathrm{℃}$$

(11)

Figure 9 shows the results of the snow process in Periods A–E using the two schemes detailed in this section. In Period A, the Lenhing scheme did not reflect the occurrence of the two snowfalls in December and displayed a smaller snow depth than the observed for the third snowfall. The simulation results of the Schmucki scheme in the first snow cover process were similar to the previous schemes and do not show the second snowfall. However, the snow depth was thick during January, and the peak value (5 cm) was 2 cm greater than the observed value.

Comparing the results of the two schemes in Figure 9b, the depth of the Schmucki scheme was larger than that of the Lehning scheme in mid-December, with a peak value of 37 cm. The simulated and observed maximum snow depths of the Lehning scheme have tremendous differences of 50 cm and 2 cm in mid-January, respectively.

From the aspects of the snow accumulation-melt process and maximum snow depth, the Schmucki scheme was generally more accurate in the snow process of Period C (Figure 9c). Compared with the Liston scheme, which had the best effects, the Schmucki scheme lagged behind the snowfall on 27 August by one day, and the magnitude of the snow depth simulation was similar to that of the Liston scheme.

As shown in Figure 9d, these two schemes showed similar performances in most cases. During the snow cover around mid-November and early March, the Schmucki scheme values were larger than those of the Lehning scheme. Compared to the observed, the maximum accumulation snow depth values were more accurate than those of the previous schemes. However, it is worth mentioning that all schemes do not reproduce the continuous snow accumulation.

Next, combined with the statistical calculation results in Figure 10 and Figure S3, we evaluated the simulation effects of two schemes for calculating the fresh snow density after adding humidity. It can be seen from the scatter plots that the Lehning scheme was smaller than that observed with a negative ME. The ME of the Schmucki scheme was close to 1 cm, and the R of the two was above 0.7. Combined with the trend of the snow process simulated in Figure 9, the simulation effect of the Schmucki scheme was better than that of the Lehning scheme in Period A. However, larger MEs and RMSEs were produced using these two schemes than those of the previous scheme, and the K_W scheme was still the best in Period C.

In simulation Periods D and E, the correlation values of the two schemes significantly decreased. Nevertheless, the magnitude of the simulated snow depth was closer to the observed depth than those of the previous schemes, with reduced RMSEs. For the continuous snow process during the D–E periods, the simulation effect of the Lehning scheme was better than that of the Schmucki scheme.

3.2. Improved Radiation and Energy Simulation with Better Fresh Snow Density Schemes

The surface radiation balance equation and energy distribution equation are as follows: Rn = Hs + LE + G₀ = DLR + DSR − ULR − USR. Within this equation, downward longwave radiation (DLR) and downward shortwave radiation (DSR) function as forcing variables that were entered into the model, while the output variables are heat flux (Hs), latent flux (LE), soil heat flux (G₀), upward longwave radiation (ULR), upward shortwave radiation (USR), and net radiation flux (Rn). In the previous results, the K_W and Schmucki schemes showed significant improvements in the snow depth simulation of the discontinuous snow process in short periods. Therefore, this section will research the changes in these variable simulations with the better fresh snow density schemes to further reveal the land–atmosphere energy exchange during discontinuous snow cover processes on the TP.

Figure S4 shows the simulated snow depths of the Schmucki, Kw, and Anderson schemes observed at Maqu. The Anderson scheme had the largest simulation, while the other two schemes displayed similar trends and magnitudes. Compared with the observed value, the amelioration of the snow depth simulation was mainly reflected in January.

Figure 11 shows the difference of the radiation and energy flux between the observed and simulated results after improving the snow depth via the Schmucki, Kw, and Anderson schemes in Maqu. For two radiation components, the error of three schemes was roughly in the range of 50 W m⁻², and the Schmucki scheme showed a more obvious difference; perhaps the lowest snow depth in the second snow period (the shortest shadow area) affects higher land surface temperatures and reflect less solar radiation, causing a distinct error. Though the earlier end of the snow period was revealed via the K_W and Schmucki schemes (17 January), a lower level of snow depth leads to a higher albedo and more reflected solar radiation, meaning that the upward shortwave radiation is related to the snow coverage, land surface coverage, and other reasons alongside snow depth in the model.

Compared with the radiation components, energy flux showed the more distinct deviation in the three schemes, especially soil heat flux. For sensible heat flux, the Anderson scheme showed the closest simulated result than the other schemes, and the lowest RMSE also represented this character (

Table 1. Calculation of related statistics between the simulated and observed radiation and surface energy fluxes after improved snow depths in the snow period of Maqu.

Scheme	Statistic	ULR	USR	Rn	Hs	LE	G0
Schmucki	ME	16.65	−24.72	8.07	4.68	8.09	−10.39
	RMSE	22.31	66.28	59.74	37.58	21.19	90.67
Kw	ME	14.30	−24.51	10.21	−0.10	16.00	−11.32
	RMSE	20.53	68.75	63.45	37.57	28.27	87.51
Anderson	ME	13.91	−24.54	10.63	−0.36	15.49	−10.06
	RMSE	20.08	68.52	63.33	37.07	27.66	87.38

). For latent heat flux, the Schmucki scheme was obviously superior to the others. The reasons for the more significant positive deviations (80~90 W m⁻² of RMSE) for soil heat flux may be that these models do not simulate the level of thinner snow cover, which can weaken the energy exchange on the TP, accompanied with error observations due to surface soil heat flux retrieved from the calculation via observations in 5 cm under the ground.

In general, the improvement of the Schmucki scheme for radiation and energy flux was better than the other schemes.

Figure S5 shows that the two schemes have a desired simulation effect on the snow depth at Madoi station. It can be seen that the Anderson scheme had the largest value, although the two schemes were slightly improved compared with the original scheme. Nevertheless, the relevant statistics have not changed, and are dozens of times larger than the observed statistics; therefore, they have almost no effect on the improvement of radiation and energy. This section of the simulation process has not been discussed further.

For the two schemes (the K_W and Schmucki schemes) with better snow depth improvement effects and the original scheme in the short period of time (Figure S6), the value of the Schmucki scheme was slightly larger than that of the Anderson scheme, and the K_W scheme was obviously closer to the actual observed values.

In comparison with the observations, since the ME of the snow depth between the K_W and Anderson schemes was within 1 cm, the simulated values of radiation and energy for these two schemes basically showed no change from the statistical results during the snow-covered period. The Schmucki scheme was slightly larger than the Anderson scheme in the snow depth simulation, and the RMSE of the radiation flux and energy flux obtained via the simulation were the lowest (

Table 2. Calculation of related statistics between the simulated and observed radiation and surface energy fluxes after improved snow depth over short periods of time at Yakou.

Scheme	Statistic	ULR	USR	Rn	Hs	LE	G0
Schmucki	ME	19.84	−71.29	51.45	44.17	53.06	16.62
	RMSE	25.66	138.73	117.55	65.60	81.73	90.41
Kw	ME	20.10	−72.30	52.20	44.92	53.57	16.33
	RMSE	25.84	139.62	118.30	66.02	81.97	90.64
Anderson	ME	20.10	−72.30	52.20	44.92	53.57	16.33
	RMSE	25.84	139.62	118.30	66.02	81.97	90.64

); that is, as the snow depth increased, the surface temperature and upward longwave radiation decreased, and the upward shortwave radiation increased. Additionally, Figure 12b shows a distinct bias in the snow period; a rather large snow depth still reflects the obvious negative bias in upward shortwave radiation and is more obvious in the higher elevation region, so it explains the particularity of snow cover on the TP; the weight of factors affecting snow albedo on the plateau should be further studied.

In summary, Maqu showed that the change in upward longwave radiation has a negative correlation with the snow depth, which indicates that when the snow depth is less than 20 cm, the snow cover has a cooling effect on the soil. The upward shortwave radiation may either decrease or increase with increasing snow depth, indicating that the influence of snow depth on snow albedo may depend on other snow parameters, such as snow particle size or a snow water equivalent [45].

As none of the schemes were able to simulate the main feature of the continuous accumulation, it is no different from the discontinuous snow cover; thus, we have no further research on the simulation of the radiation flux and energy flux under the continuous snow process.

4. Discussion

(1)

Although we have found several schemes that improved simulation effects on the discontinuous snow depth, all schemes simulated larger snow depths during the accumulation process and earlier ablation characteristics of continuous snow cover process. The Tibetan Plateau (TP) has a dry and cold climate, resulting in a low snow density under low temperature and humidity conditions. Consequently, the model produces exaggerated snow depths due to the same rainfall forcing data. Additionally, during snow melting, the melted water penetrates the snow layer and freezes into ice, leading to a significant increase in snow density. However, due to climate dryness on the TP, it causes the model to overestimate the rates of snowmelt and snow density caused by the refreezing of meltwater [62].

(2)

The snow cover fraction (SCF) is usually less than 100%; yet, for local scale-like site observations, the SCF should be 100% once the snow depth is more than a certain threshold. Otherwise, a less than 100% SCF will remarkably reduce surface albedo and enhance snow melting, leading to smaller SCFs; this positive feedback leads to quick snow ablation in continuous snow cover processes. Furthermore, solar radiation, air temperature, precipitation, and the orientation of the slope are crucial factors in the snowmelt modeling process [63,64,65]. Therefore, other parameterization schemes in the model need to be further evaluated to improve the land surface process model for continuous snow cover on the TP in the future.

(3)

Due to the limited availability of snow depth and precipitation observation data in the field at the Maqu and Madoi stations, this study utilized the products of snowfall and precipitation at the station where the China Meteorological Administration was located, or where the precipitation was observed in the field. The analyzed data were strictly controlled, but it was still impossible to completely avoid the deviation of the model from the snow simulation due to the insufficient accuracy of the precipitation-forced input. At present, various remote sensing snow monitoring data can compensate for the shortage of observation stations and snow parameters on the plateau [66,67,68]. We hope to fully utilize these two types of data in the future and carry out more comprehensive and accurate studies on the plateau snow.

(4)

Further improvements are required in the calculation method used in this study to compare surface soil heat flux with the model outputs. From the comparison, the calculation results of soil heat flux were very different from the simulation, and the energy non-closure rate during the snow accumulation process was large [69]. A detailed consideration of energy dissipation during transmission processes will enable a better verification and evaluation of land surface process model capabilities.

5. Conclusions

This study compares the simulated snow depth of ten parameterized schemes of fresh snow density with CLM4.5 during a total of five simulation periods at three field observation stations on the eastern TP and radiation and energy flux with the observed results. The following are the primary conclusions:

(1)

All schemes can simulate the accumulation and ablation process of discontinuous snow cover in a short period at three stations, and the K_W scheme (adding the wind speed calculation component) was the best-performing scheme regarding snow depth simulation. However, no scheme was able to adequately describe the continuous accumulation of snow cover with the characters of repeated accumulation and ablation at Madoi and Yakou for a long period of time. Under the circumstances, the Lehning scheme performed marginally better than other schemes in terms of snow depth during continuous snow.

(2)

Compared with other schemes, the simulation effect of the Schmucki scheme on radiation flux and energy flux under discontinuous snow cover was significantly improved.

(3)

According to the simulation effects of the improved fresh snow density scheme on the radiation and energy flux in the discontinuous snow process, the change in upward longwave radiation has a negative correlation with the snow depth; that is, the snow accumulation has an overall effect on cooling when the snow depth is less than 20 cm, and the upward shortwave radiation may either decrease or increase with increasing snow depth. The simulated depth of snow has a considerable effect on the heat exchange (sensible heat flux and latent heat flux) between the ground and air, demonstrating that the effective thermal conductivity in the model is sensitive to varying snow density responses.