Validation of the SNTHERM Model Applied for Snow Depth, Grain Size, and Brightness Temperature Simulation at Meteorological Stations in China

: Validation of the snow process model is an important preliminary work for the snow parameter estimation. The snow grain growth is a continuous and accumulative process, which cannot be evaluated without comparing with the observations in snow season scale. In order to understand the snow properties in the Asian Water Tower region (including Xinjiang province and the Tibetan Plateau) and enhance the use of modeling tools, an extended snow experiment at the foot of the Altay Mountain was designed to validate and improve the coupled physical Snow Thermal Model (SNTHERM) and the Microwave Emission Model of Layered Snowpacks (MEMLS). By matching simultaneously the observed snow depth, geometric grain size, and observed brightness temperature (T B ), with an RMSE of 1.91 cm, 0.47 mm, and 4.43 K (at 36.5 GHz, vertical polarization), respectively, we ﬁnalized the important model coe ﬃ cients, which are the grain growth coe ﬃ cient and the grain size to exponential correlation length conversion coe ﬃ cients. When extended to 102 meteorological stations in the 2008–2009 winter, the SNTHERM predicted the daily snow depth with an accuracy of 2–4 cm RMSE, and the coupled SNTHERM-MEMLS model predicted the satellite-observed T B with an accuracy of 13.34 K RMSE at 36.5 GHz, vertical polarization, with the fractional snow cover considered.


Introduction
The seasonal snow cover, as an important component of the cryosphere, plays a key role in the climate system [1,2]. The high albedo of the snow surface, the latent heat generated by the internal ice/water phase transition, and the adiabatic effect of the snow layer significantly shift the energy We reached the site on October 22th, 2017, and installed the L, C, and X band radiometer from the Institute of Remote Sensing and Digital Earth (RADI) and the L, Ku, and Ka band radiometer from the Northwest Institute of Eco-Environment and Resources (NIEER). The observed field has a dimension of 50 m in the North-South direction and 30 m in the East-West direction. As can be seen from Figure 1, the grass was mowed carefully before the first snowfall to prepare a bare soil background for the microwave modeling. On the north edge of the observed field, from west to east, the RADI radiometer, a ground-based synthetic aperture radar system (SAR) and the NIEER radiometer were installed in sequence. The radiometers provided the brightness temperature measurements at 1.4, 6.925, 10.65, 18.7, and 36.5 GHz at 2 to 3 h steps. Manufactured by the RPG-XCH-DP (Radiometer Physics GMbH) company, the radiometers installed a stable noise diode that support an automatic gain calibration to keep their stability for several months. This noise injection module, together with the temperature stable reference at the other side of the Dicke switch, provides two calibration points to correct the system noise equivalent temperature (T n ) and the system gain (G). The non-linearity of the system can be assumed unchanged for its intrinsic stability. Both the RADI and the NIEER radiometers were sky-tipping calibrated when they were first installed. Later, the manually controlled sky tipping calibrations were implemented for the RADI radiometer on clear-sky conditions approximately per month. The NIEER radiometer was not further sky tipping calibrated, because it is found due to the sensitivity of 36.5 GHz to water vapor content, the sky tipping calibration introduced more errors and was no better than the automatic calibration.
On November 10th, 2017, the first snowfall covered the observed field. Starting from this day to March 11th, 2018, a total of 109 snowpits were excavated within 50 m distance to the radiometers. For each snowpit, different snow layers were recognized by observing the color and crystal size of the snow and feeling the snow hardness by fingers. For each snow layer, the thickness, geometric grain size, and temperature were measured. The thickness was measured by a woolen folding ruler. The geometric grain size was measured by post-processing photos of a large amount of snow particles on a 2 mm × 2 mm crystal card taken by a high-resolution Leica camera. The snow density was measured by a 250 mL wedge-type sampler (10 cm × 10 cm × 5 cm) at per 5 cm step. Before the first snowfall, three Decagon 5TM soil sensors were buried at −2, −5, and −10 cm depth below the soil surface to provide the soil temperature and the equivalent volumetric water content (unfrozen part) by measuring the dielectric constant of the soil. The data was collected by an Em−50 series logger.
Remote Sens. 2020, 12, 507 5 of 29 two calibration points to correct the system noise equivalent temperature (Tn) and the system gain (G). The non-linearity of the system can be assumed unchanged for its intrinsic stability. Both the RADI and the NIEER radiometers were sky-tipping calibrated when they were first installed. Later, the manually controlled sky tipping calibrations were implemented for the RADI radiometer on clear-sky conditions approximately per month. The NIEER radiometer was not further sky tipping calibrated, because it is found due to the sensitivity of 36.5 GHz to water vapor content, the sky tipping calibration introduced more errors and was no better than the automatic calibration.   Figure 2b-e shows the corresponding snow particle photos for each layer. The geometric snow grain size was post-processed using the method as follows. First, a human-computer interactive program was coded to georegistrate the photo using the 2-mm grid intersections. Then, an interactive dialog box was called to save the two ends of the maximum dimension (diameter) of one particle visually recognized and clicked by the observer. Only two observers were involved in the snow particle digitization, and the final result was checked and corrected by a single one. The X, Y coordinates of the paired points were converted to the lengths by the computer. Each photo was separated into 16 subphotos, and in each subphoto at least 10 particles were digitalized. Using the over 160 particles, the geometric grain size (D max ) was calculated using the mean of the maximum dimensions, according to the definition in Mätzler (2002) [29].
Remote Sens. 2019, 11, x FOR PEER REVIEW 5 of 29 Figure 1. The Altay experiment: (a) the snow covered field observed by the instruments on the north edge, from left to right are: the Remote Sensing and Digital Earth (RADI) radiometer, a ground-based synthetic aperture radar system (SAR), and the Northwest Institute of Eco-Environment and Resources (NIEER) radiometer; (b) the grass mowed field before the first snowfall.
On November 10th, 2017, the first snowfall covered the observed field. Starting from this day to March 11th, 2018, a total of 109 snowpits were excavated within 50 m distance to the radiometers. For each snowpit, different snow layers were recognized by observing the color and crystal size of the snow and feeling the snow hardness by fingers. For each snow layer, the thickness, geometric grain size, and temperature were measured. The thickness was measured by a woolen folding ruler. The geometric grain size was measured by post-processing photos of a large amount of snow particles on a 2 mm × 2 mm crystal card taken by a high-resolution Leica camera. The snow density was measured by a 250 mL wedge-type sampler (10 cm × 10 cm × 5 cm) at per 5 cm step. Before the first snowfall, three Decagon 5TM soil sensors were buried at −2, −5, and −10 cm depth below the soil surface to provide the soil temperature and the equivalent volumetric water content (unfrozen part) by measuring the dielectric constant of the soil. The data was collected by an Em−50 series logger. Figure 2 shows an example of the snowpit observed on 15:34 January 24th, 2018. The snow had 4 layers with a thickness of 4.5 cm, 2 cm, 2.5 cm, and 7.5 cm from top to bottom. Figure 2b-e shows the corresponding snow particle photos for each layer. The geometric snow grain size was postprocessed using the method as follows. First, a human-computer interactive program was coded to georegistrate the photo using the 2-mm grid intersections. Then, an interactive dialog box was called to save the two ends of the maximum dimension (diameter) of one particle visually recognized and clicked by the observer. Only two observers were involved in the snow particle digitization, and the final result was checked and corrected by a single one. The X, Y coordinates of the paired points were converted to the lengths by the computer. Each photo was separated into 16 subphotos, and in each subphoto at least 10 particles were digitalized. Using the over 160 particles, the geometric grain size (Dmax) was calculated using the mean of the maximum dimensions, according to the definition in Mätzler (2002) [29]. From the Altay meteorological station, we acquired the hourly air temperature, precipitation, wind speed, relative humidity, air pressure, and solar radiation to drive SNTHERM. The hourly soil temperature was also provided. We utilized the measured soil temperature at −320 cm depth as the low boundary condition, whereas the upper boundary condition is the air temperature. The required downwelling longwave radiation forcing data was from the 3-h 0.25 degree resolution Noah model Global Land Data Assimilation System (GLDAS) dataset [44]; its influence will be discussed by comparing the SNTHERM simulated and the station measured net radiation in Section 5.2.

The Study Area and the Datasets Utilized
The study area is the Xinjiang, Qinghai, and Xizang (Chinese name, English name is: Tibet) provinces in China, which can be considered as the extended Asian Water Tower (AWT) region. These three provinces are far from the ocean and have a dry weather. The solar radiation is stronger on the high-elevation Tibetan Plateau (TP), which is not beneficial to retain the snow cover. In the flat area, snow cover in the northern part of Xinjiang can last for over 5 months, whereas most new snow disappeared within 2 weeks on the TP unless the snowfall is heavy enough or the site is shaded, north-faced, or on high mountains. The snow disappearance is attributed from the strong sublimation, melting of the topmost snow layer under strong solar radiation (although the air temperature is below 0 • C), and from the wind blowing effect. The night on the TP is cold in the winter months, and it can refreeze the snow surface to maintain a dry snow condition at AMSR-E descending overpass time (1:30 a.m. local time approximately).
The accuracy of the snow parameter estimation is determined by the accuracy of the models themselves and the forcing data. From the National Meteorological Information Center of the China Meteorological Administration (CMA), we got the daily meteorological datasets to support our research. Figure 3 is a map of three provinces marked by the locations of the reference meteorological sites. The entire meteorological dataset given to us has 134 sites, whereas we used only 102 sites. At some sites, there was no snowfall in the 2008-2009 winter (magenta pentagon), average snow depth below 1 cm (red triangles), or incomplete data to support SNTHERM running (blue stars); these sites were not used. Most of the sites not used locate around the Taklimakan Desert in the southern Xinjiang or in the Qaidam Basin in the northern Qinghai.
The daily dataset from the meteorological stations was merged with the 3-h 0.25 degree resolution Noah model Global Land Data Assimilation System (GLDAS) dataset [44]. The SNTHERM 3-h meteorological forcing is produced as follows. For the air temperature, the measured daily mean, maximum and minimum measurements were provided by the stations. On each day, the GLDAS per 3 h simulations were translated to the daily statistics. Using the three points (mean, maximum, and minimum), a simple linear model, as T measured = a + b × T GLDAS , was used to generate the corrected GLDAS air temperature, with a and b fitted. For the relative humidity and wind speed, the GLDAS values were corrected by comparing with the daily mean measurements, using a scaling model as Y measured = a × Y GLDAS (Y represents the parameter). For the precipitation, the same scaling model is utilized to match the observed per 12 h total precipitation (8:00-20:00 and 20:00-8:00). SNTHERM requires only one average air pressure, and thus we utilized the average daily mean air pressure measured from November 2008 to March 2009.
The only tuned parameter is the radiation flux. Only 13 sites have the daily mean solar radiation measurements. Four of them locate in Tibet, 4 in Qinghai, and 5 in Xinjiang provinces. For these sites, the measured daily solar radiation was merged with GLDAS using the scaling model. The statistics showed the GLDAS daily solar radiation has a 20% relative error at most. Therefore, for the other sites, it is allowed to tune the time-series downwelling shortwave and longwave radiation on all days by a same factor between 0.8 to 1.2. There is no further artificially tuning for certain specific days.
The satellite AMSR-E T B observations used for validation are from the reprocessed MEaSUREs Calibrated Enhanced-Resolution Passive Microwave Daily EASE-Grid 2.0 Brightness Temperature (EnTB) [45]. The nominal spatial resolution is 3.125 km at 36

Snow Thermal Model (SNTHERM) Introduction
The SNTHERM (Snow THERmal Model) [22] was published by Rachel Jordan working in the Cold Regions Research & Engineering Laboratory (CRREL) in 1989. It is a one-dimensional model, which aimed to accurately predict snow temperature and snowmelt using the meteorological forcing data of air temperature, relative humidity, air pressure, wind speed, precipitation, and radiation fluxes. The model can calculate the snow albedo using the simulated snow grain size and the upwelling longwave radiation using the simulated snow temperature. We activated these modules, and therefore only the downwelling shortwave and longwave radiations are needed.
In SNTHERM, the snow-covered land surface is modeled as several layers of snow and soil medium, with the vegetation neglected. Its mass balance includes the precipitation (snow, rain, sleet), sublimation, and evaporation at the snow surface and the liquid water and water vapor movement in the snow layers. Its energy balance includes the radiation balance, turbulent fluxes (sensible + latent heat), and convective heat from precipitation at the snow surface, soil heat flux at the bottom of the snow, and the heat conducted and carried by mass fluxes in the snow layers. Most snow process models, complex or simple, considered the same flux terms. A the diagrammatic explanation of the SNTHERM physical processes was published on the CRREL website (https://www.erdc.usace.army.mil/Media/Fact-Sheets/Fact-Sheet-Article-View/Article/476650/sntherm/) and can be referred for a better understanding of these fluxes.
Within the snowpack, the SNTHERM considers the snow temperature and phase change, snow compaction and the snow grain growth. The grain growth module is unique, dependent on snow temperature, but simpler than the SNOWPACK and Crocus models, as:

Snow Thermal Model (SNTHERM) Introduction
The SNTHERM (Snow THERmal Model) [22] was published by Rachel Jordan working in the Cold Regions Research & Engineering Laboratory (CRREL) in 1989. It is a one-dimensional model, which aimed to accurately predict snow temperature and snowmelt using the meteorological forcing data of air temperature, relative humidity, air pressure, wind speed, precipitation, and radiation fluxes. The model can calculate the snow albedo using the simulated snow grain size and the upwelling longwave radiation using the simulated snow temperature. We activated these modules, and therefore only the downwelling shortwave and longwave radiations are needed.
In SNTHERM, the snow-covered land surface is modeled as several layers of snow and soil medium, with the vegetation neglected. Its mass balance includes the precipitation (snow, rain, sleet), sublimation, and evaporation at the snow surface and the liquid water and water vapor movement in the snow layers. Its energy balance includes the radiation balance, turbulent fluxes (sensible + latent heat), and convective heat from precipitation at the snow surface, soil heat flux at the bottom of the snow, and the heat conducted and carried by mass fluxes in the snow layers. Most snow process models, complex or simple, considered the same flux terms. A the diagrammatic explanation of the SNTHERM physical processes was published on the CRREL website (https: //www.erdc.usace.army.mil/Media/Fact-Sheets/Fact-Sheet-Article-View/Article/476650/sntherm/) and can be referred for a better understanding of these fluxes.
Within the snowpack, the SNTHERM considers the snow temperature and phase change, snow compaction and the snow grain growth. The grain growth module is unique, dependent on snow temperature, but simpler than the SNOWPACK and Crocus models, as: According to the magnitude of the grain size measurements (written as '0.04-0.2 mm for new snow, 0.2-0.6 mm for fine-grained older snow and 2.0-3.0 mm for older snow near 0 • C') and the time of the study (in 1960s to 1980s), the definition of the SNTHERM d is closer to the geometric grain size (D max ).
In [22], it suggested a default value of g 1 as 5.0 × 10 −7 m 4 /kg. Due to the semi-empirical natural of g 1 , in this paper, it is considered as a revisable parameter and is planned to be updated according to the Altay snowpit measurements.

Snow Thermal Model (SNTHERM) Improvement
The published SNTHERM code from the CRREL website was composed in 1989 and updated in 2002. Starting from 2010, we did a series of revisions. Some revisions are related to bug fixing. For example, one line was added as follows: after the snow density is calculated by snow mass divided by layer thickness, an instant check of the calculated density is carried out-if the snow density is larger than 917 kg/m 3 , reduce it to 917 kg/m 3 and correct the layer thickness to keep snow mass unchanged. A maximum snow layer thickness limit of 5 cm was added to prevent an abrupt melting rate when a very thick layer suddenly reaches 0 • C. The use of SNTHERM module to calculate snow albedo instead of a user-defined constant value was set as a mandatory option. More important, the SNTHERM89 code has a module to consider the influence of soil background on the shallow snow albedo, written in the shallow.f function. It is of critical importance for simulations in our study region and was activated.
In the original SNTHERM code, although the snowmelt flux is calculated at the bottom-most snow layer, such flux was not inputted to the topmost soil layer. The soil moisture can only be changed due to water vapor movement, which is of tiny magnitude compared to the liquid water movement. The SNTHERM can be run starting from an initial snowpack status. However, to reduce the reliance on the in-situ measurements, our choice is to start SNTHERM from before the snow season. The rainfall and the snow melted water can change the soil water mass. The later further influences the thermal properties of the soil layer, and consequently influences the snow temperature via the soil heat flux. Therefore, for the sake of physical validity, we thought the simulation of the soil water movement process is needed, and should be implemented for all soil layers for completeness. It can be added to the SNTHERM modeling framework by providing the liquid water mass flux, U l [kg m −2 s −1 ], at the bottom of each soil layer. In this case, the downward drainage of water from the snowpack into the soil layers is allowed.
According to Darcy's law, U l can be calculated as [46,47], where U l is defined as positive in the upward direction. ρ l is the liquid water density, 1000 kg/m 3 . k l is the hydraulic conductivity [m/s]. ψ is the soil matric potential [m], defined in negative values (otherwise, the sign before ∂ψ/∂z in Equation (2) should be negative). z the vertical position [m]. ∂ψ/∂z is unitless. We utilized the Van-Genuchten (1980) (VG80) soil-water retention curve model [48] to calculate k l and ∂ψ/∂z. The hydraulic conductivity is a product of saturated hydraulic conductivity (k s ) and relative hydraulic conductivity (k r ) as: where, S e is the effective liquid water saturation, defined in SNTHERM and also in [48] as: where, θ is the volumetric soil water content [m 3 /m 3 ], and θ s and θ r are the saturated and the irreducible volumetric soil water content. s is the unitless liquid water saturation, defined as θ/φ. φ is the soil porosity, calculated as 1 − ρ b /ρ s , where ρ b is the soil bulk density [kg/m 3 ] and ρ s is the soil specific density [kg/m 3 ]. S e based on θ or s is the same for soil. m is the VG80 model shape factor, and has a relationship as m = 1 − 1/n with the other shape factor, n. The soil matric potential is related to S e as [49]: where α [m −1 ] is a VG80 model parameter, defined as positive hereafter; it is the inflection point where the slope of ∂θ/∂ψ reaches its maximum value. Then, Following the SNTHERM modeling framework, ∂S e /∂z was discretized to (S e,i − S e,i−1 )/(0.5 × (∆z i + ∆z i−1 )), and inputted into a big matrix to solve S e at each layer. We also added a further adjustment of k l for frozen soil according to Equation (5) of [50].
U l calculated from the bottom-most snow layer is used as a boundary condition for the topmost soil layer. If j is the topmost soil layer, its discretized water balance equation at time t is: where ∆t is time step, ∆z is layer thickness. c7 is the conversion factor from S e to liquid water mass change rate, and M is the phase change term; see Jordan (1991) [22] (pp. 30) for more details. Originally, the SNTHERM code has only two default soil type choices: sand and clay, and only the thermal parameters were given. To check the soil model performance, we conducted three groups of calculations, as listed in Table 1. In the 1st scenario, we used the original SNTHERM89 code choosing a clay soil type. In the 2nd scenario, we added the soil water movement for clay soil, and assumed a soil texture of 22% sand, 20% silt, and 58% clay fractions, according to the clay soil setting in the CLM model [51]. Then, the soil hydraulic parameters (θ s , θ r , α, n, and m) were able to be calculated using the Carsel and Parrish (1988) model listed in [49] (eq. C3). In the 3rd scenario, the soil texture (35% sand, 46% silt, and 19% clay fractions) and porosity (0.51) at the Altay station was extracted from the GLDAS soil dataset (https://ldas.gsfc.nasa.gov/gldas/soils) to calculate the thermal coefficients using [51] and VG80 model parameters using [49]. During the thermal coefficient calculation, an equivalent organic fraction of 0.1033 was reverted from the texture and porosity, using Equation (3) in [51], and used to calculate the dry soil thermal conductivity, thermal conductivity of dry solids, and heat capacity, using Equations (7), (10), and (14) in [51], respectively. Table 1 lists the values of these parameters. In the 1st scenario when the original SNTHERM code was utilized, it was found a minimum extent of bug fixing (3 lines) is needed to produce a normal, comparable simulation result. We listed these debugging efforts in the Appendix A. To maintain the consistency with scenarios 2 and 3, in scenario 1, we also turned on the snow albedo calculation module and the shallow snow albedo correction module.
In simulation for the Altay experiment, the time-varying soil temperature at −320 cm depth is used as the boundary condition during the temperature calculation. The initial status of soil temperature was set according to the measurements at −5, −10, −15, −20, −40, −80, −160, and −320 cm depth. The Decagon 5TM sensor only measured the soil moisture of the top 10 cm, as about 10% volumetric water content. The volumetric water content in the 1st and 2nd scenarios were set homogenous as 10%. A constant volumetric water content at −320 cm depth was set as the boundary condition for the S e calculation. We have found setting the volumetric water content below −22 cm depth as 20% gave better temperature simulation for the top soil layer, and therefore this setting was used for the 3rd scenario. The reason is uncertain, but could have compensated for lots of impacting factors, like the inhomogeneous soil hydraulic parameters and time-varying soil moisture boundary condition.
The SNTHERM simulation for the 102 stations started from August 1st 2008 to August 1st 2009. the initial status of soil was set according to GLDAS data on August 1st and a constant soil temperature and moisture boundary condition was utilized. The soil parameters were set the same as the Altay site-specific soil values. An improvement in future can be utilizing the local GLDAS soil parameters or the Wei et al. (2013) Tables 2 and 3 of [53].

Brightness Temperature Model
The snow brightness temperature was simulated by the Microwave Emission Model of Layered Snowpacks (MEMLS)   [39] based on the Improved Born Approximation (IBA) [40]. It is a multiple layer snow radiative transfer model, which uses the temperature, snow micro-structure parameter (exponential correlation length, p ex , in its case), density, and thickness of all snow layers and the background soil reflectivity at vertical and horizontal polarizations to calculate the brightness temperature observed from the snow surface. The IBA considers the multiple-scattering of diffuse radiation between the snow particles. In MEMLS, when the snow radiation is emitted out from the air-snow boundary, the trapped radiation due to total internal reflection (TIR) is considered using 6-flux radiative transfer equation resolution. The MEMLS model details are explicitly presented in [39,40], and a further model explanation can also be found in [43].
The soil reflectivity model was the frequency-independent QHNfi model written in [54]. It is a semi-empirical model with the coefficients recently refitted using the field measurements covering a wide range of bare soil conditions and frequencies at 1.4 to 90 GHz.
MEMLS requires an exponential correlation length (p ex ) whereas SNTHERM simulates the geometric grain size (D max ). We utilized the equation and coefficients validated in [43] to do the D max to p ex conversion: The brightness temperature simulated by MEMLS is for fully-covered snow. In order to represent the AMSR-E brightness temperature in the pixel that nests the meteorological station, a fractional snow cover (FSC) retrieved from the 7-band MODIS reflectance [55,56] is used to calculate the pixel T B as: where e background is the background emissivity at certain frequency calculated before the first snowfall when air temperature is close to 0 • C. T background in theory should be the background physical temperature uncovered by snow. It is not the SNTHERM-simulated topmost soil temperature beneath the snowpack because the uncovered soil does not have the thermal insulation from the snow. T background should have a stronger sensitivity but less fluctuation compared to the air temperature. As a compromise, we utilized the top snow layer temperature. Figure 4 shows a comparison of the simulated snow and soil state parameters using the original SNTHERM89 code (scenario 1) and the revised SNTHERM with the soil water movement considered (scenarios 2 and 3), driving by the same meteorological forcing data. The default g 1 of 5.0 × 10 −7 m 4 /kg was utilized for simulation. Figure 4 also includes the field measurements if available. For the temperature variables, the daily-averaged mean values were plotted for a better presentation compared to strongly fluctuated hourly values.

The Revised Soil Modeling in SNTHERM and the Validation Based on the Altay Experiment
From Figure 4a, it shows the simulated snow depth from three scenarios is close, except the original SNTHERM89 code gives a slower melt for the last few days. In the snow accumulation phase, the difference in the simulated snow depth based on different soil settings is tiny, in the magnitude of 0.7 cm in maximum. The simulated snow depth matches the observed snow depth. On the other hand, during the snowmelt phase, the simulated snow by scenarios 2 and 3 (NewCode-Clay and NewCode-Altay) completely melts on 22:00 and 21:00, respectively, March 16th, whereas the measurement from the meteorological site gives a first snow-free day on March 13th. The end of snow season simulated by SNTHERM delays by 3 days. It will be improved by correcting the bias in the snow grain size estimation, as will be shown later in Figure 8a.
In Figure 4b, the simulated average snow density is also close among 3 soil parameter scenarios, and matches the observations from the snowpits. Figure 4c-e show the soil thermal properties influence the simulated topmost soil layer temperature (c), the average snow temperature (d), and the average snow temperature gradient ((T snow,bottom − T snow,top )/SD) (e). All the lines presented here are for daily average values. In Figure 4c, the daily average top soil layer temperature measured by the station and by our 5TM sensor are also added. It shows the use of the clay soil type (scenarios 1 and 2) underestimates the soil temperature for over 10 K in the coldest time. A significant improvement is observed by using the Altay site-specific soil type and a volumetric liquid water content of 20% for deep soil layers at the initial time. The simulation result using a homogeneous 10% volumetric liquid water content profile is between scenarios 1 and 3 (not shown here). Scenario 3 (NewCode-AltaySoil) gives a higher top soil temperature and because the topmost snow temperature approximates the air temperature, it gives a higher averaged snow temperature ( Figure 4d) and stronger temperature gradient (Figure 4e). According to Equation (1), a higher temperature and temperature gradient can accelerate the grain growth speed. Therefore, scenario 3 predicted a slight higher bottom geometric snow grain size (D max ) in Figure 4f, but is not always higher for all days compared to scenario 2. However, the diversity in D max is not as strong as we speculated. The resulted diversity in the predicted brightness temperature will be discussed in Section 5.1. Figure 4f shows all three scenarios underestimated the observed D max from the snowpits. In order to support the update of grain growth coefficient (g 1 ), an accurate estimation of the snow temperature and snow temperature gradient is the prerequisite condition. In Figure 5, the calculated and the observed temperature gradient from snowpits are compared. The later was calculated by the difference of snow temperature divided by the difference of the measured heights (−∂T/∂z). It shows the simulated temperature gradient has a strong diurnal variation-low in the daytime and high in the night time. The highest temperature gradient can reach 2 K/cm, whereas that at noon it is in general around 0.5 K/cm or can even be negative. Two reasons resulted in a measured smaller temperature gradient compared to the simulations. First, most measurements of snow temperature were at the daytime, and reasonably it reaches the low limit of the simulations. Second, the measured temperature of the bottommost snow layer can be biased low. The insertion of the temperature sensor easily broke the sugary depth hoar at the bottom and can mix cold air into the snow around the sensor. We think the temperature of the topmost soil layer is a more reliable validation reference than the bottom snow temperature measurements. Figure 6 compares the instant snow layer temperature and the SNTHERM snow temperature. It showed the temperature using the clay soil (scenarios 1 and 2) biased low for 2 to 3 K. The prediction using the Altay site-specific soil (scenario 3) is better. As a conclusion, the temperature simulation using scenario 3 improves over scenarios 1 and 2, and can be used for the grain growth coefficient adjustment. Figure 4g,h show the simulated total water content and liquid water content for the topmost soil layer. In Figure 4g, the total water content using SNTHERM89 nearly changes due to model defect as speculated. When the soil water movement is added, both the scenario 2 using clay soil and the scenario 3 using Altay site-specific soil can reach the saturated volumetric water content (0.36 for scenario 2 and 0.42 for scenario 3) after snowmelt or a rainfall. However, the drainage of the Altay soil is quicker than clay. In combination with a slower freezing speed, the Altay soil type gives a smaller total water content during the snow-covered days. We had some total gravimetric water content (m g ) measurements from drying the extracted soil sample by the oven. The measurements showed 9-12% m g in the 0-4 cm depth between November 7th to 16th. It was around 5% m g after November 22nd; however, after the soil was frozen, a representative sample was difficult to take. The measurements in mid-November indicate that the simulation using Altay soil is more plausible.
From Figure 4h, it shows the temperature has a strong influence on the liquid water content. Scenarios 1 and 2 give similar simulations because their simulated soil temperature is close and the soil freezing curve parameter (plasticity index) is the same. The plasticity index of Altay soil was adjusted to match the 5TM observed equivalent liquid water content.

The Grain Size Improvement and Brightness Temperature Modeling for Altay Experiment
After the snow temperature and the temperature gradient were accurately provided using the Altay site-specific soil parameters, the next step is to improve the grain growth coefficient (g1). In Figure 7a, the simulated Dmax from scenario 3 is plotted against the measured Dmax from the snowpits. The SNTHERM simulated much more snow layers (1-34 layers) compared to the in-situ measurements (1-5 layers). Here we extracted the SNTHERM snow layers within the height range of the snowpit layer and took the average. Currently, using the default g1, the simulated Dmax biased low,

The Grain Size Improvement and Brightness Temperature Modeling for Altay Experiment
After the snow temperature and the temperature gradient were accurately provided using the Altay site-specific soil parameters, the next step is to improve the grain growth coefficient (g1). In Figure 7a, the simulated Dmax from scenario 3 is plotted against the measured Dmax from the snowpits. The SNTHERM simulated much more snow layers (1-34 layers) compared to the in-situ measurements (1-5 layers). Here we extracted the SNTHERM snow layers within the height range of the snowpit layer and took the average. Currently, using the default g1, the simulated Dmax biased low, Figure 6. Comparison of simulated and observed layered snow temperature using three soil parameter settings: (a) scenario 1 using the original SNTHERM89 code with clay soil type, (b) scenario 2 using the new SNTHERM code with clay soil type, (c) scenario 3 using the new SNTHERM code with Altay site-specific soil type.

The Grain Size Improvement and Brightness Temperature Modeling for Altay Experiment
After the snow temperature and the temperature gradient were accurately provided using the Altay site-specific soil parameters, the next step is to improve the grain growth coefficient (g 1 ). In Figure 7a, the simulated D max from scenario 3 is plotted against the measured D max from the snowpits. The SNTHERM simulated much more snow layers (1-34 layers) compared to the in-situ measurements (1-5 layers). Here we extracted the SNTHERM snow layers within the height range of the snowpit layer and took the average. Currently, using the default g 1 , the simulated D max biased low, for about −0.55 mm in average (RMSE = 0.74 mm). The ratio of the simulated to the measured D max is 0.69 ± 0.22.
Assuming the geometric grain size at time (t) = 0 is d 0 , Equation (1) has an analytical solution as: Because the model was run before the first snowfall, it called an empirical function to calculate the fresh snow density and used the fresh snow density to calculate d 0 . d 0 is 0.16~0.17 mm for a fresh snow density of 50~200 kg/m 3 ; it makes d 0 negligible in Equation (10). As a result, to scale d by a factor of b approximately equals to multiplying g 1 by a factor of b 2 . B = 1/0.692 = 1.445. Therefore, an improved g 1 should be 5.0 × 10 −7 m 4 /kg * 1.445 × 1.445 = 1.044 × 10 −6 m 4 /kg. Figure 7b showed the D max comparison using the updated g 1 . After revision, the mean bias is reduced to −0.15 mm (RMSE = 0.47 mm). The ratio of the simulated to the measured D max is changed to 0.97 ± 0.33. By changing g 1 to match the observed geometric grain size, a well-defined physical meaning is introduced to the SNTHERM snow grain size, which was previously considered as a model-specific, effective grain size. Figure 7c  Assuming the geometric grain size at time (t) = 0 is d0, Equation (1) has an analytical solution as: Because the model was run before the first snowfall, it called an empirical function to calculate the fresh snow density and used the fresh snow density to calculate d0. d0 is 0.16~0.17 mm for a fresh snow density of 50~200 kg/m 3 ; it makes d0 negligible in Equation (10). As a result, to scale d by a factor of b approximately equals to multiplying g1 by a factor of b 2 . B = 1/0.692 = 1.445. Therefore, an improved g1 should be 5.0 × 10 −7 m 4 /kg * 1.445 × 1.445 = 1.044 × 10 −6 m 4 /kg. Figure 7b showed the Dmax comparison using the updated g1. After revision, the mean bias is reduced to −0.15 mm (RMSE = 0.47 mm). The ratio of the simulated to the measured Dmax is changed to 0.97±0.33. By changing g1 to match the observed geometric grain size, a well-defined physical meaning is introduced to the SNTHERM snow grain size, which was previously considered as a modelspecific, effective grain size. Figure 7c presented a comparison of the histogram of the simulated to measured Dmax ratio. Figure 7d,e showed more details of the time-series varied Dmax trend at different layers in the bottom-first order. After the grain growth coefficient (g1) was updated, the simulated snow depth and geometric grain size and density profiles are presented in Figure 8. In Figure 8a, after the g1 was improved, the simulated snow melts on 23:00 March 13th, 2018. Comparing with the measured snow-free day on March 13th (daily scale), the error is within one day. The new simulation has been significantly improved over the originally-simulated March 16th. It demonstrated that, after the grain size is accurately predicted to higher values, the snow albedo is reduced to promote a quicker snow melt. In Figure 8a, the snow depth prediction accuracy is −1.49 cm mean bias (MB) and 1.91 cm root-meansquared-error (RMSE). Figure 8b,e show the snow grain size and the snow density profiles. In Figure 8c,e, depthweighted average values and the values at the surface and the bottom of the snow are presented. The fluctuation of the surface layer geometric grain size and density is due to new snowfalls. After a snow layer is initialized from a new snowfall, in theory, the snow grain size monotonically increases. However, the combination of snow layers can reset the simulated values. For the bottom snow layer, two reasons from the model numerical implementation can result in a decreased snow grain size and density. Firstly, when the bottom layer is compressed to a thickness smaller than the limit, the bottom-most layer will be combined with the 2nd bottom layer and the snow parameters are averaged. Such fluctuations can also be found using the SNTHERM89 code. Secondly, the bottom snow mass can be reduced in SNTHERM via the water vapor movement, and this can result in a 'disappearance' of the bottom layer. After the grain growth coefficient (g 1 ) was updated, the simulated snow depth and geometric grain size and density profiles are presented in Figure 8. In Figure 8a, after the g 1 was improved, the simulated snow melts on 23:00 March 13th, 2018. Comparing with the measured snow-free day on March 13th (daily scale), the error is within one day. The new simulation has been significantly improved over the originally-simulated March 16th. It demonstrated that, after the grain size is accurately predicted to higher values, the snow albedo is reduced to promote a quicker snow melt. In Figure 8a, the snow depth prediction accuracy is −1.49 cm mean bias (MB) and 1.91 cm root-mean-squared-error (RMSE). Figure 8b,e show the snow grain size and the snow density profiles. In Figure 8c,e, depth-weighted average values and the values at the surface and the bottom of the snow are presented. The fluctuation of the surface layer geometric grain size and density is due to new snowfalls. After a snow layer is initialized from a new snowfall, in theory, the snow grain size monotonically increases. However, the combination of snow layers can reset the simulated values. For the bottom snow layer, two reasons from the model numerical implementation can result in a decreased snow grain size and density. Firstly, when the bottom layer is compressed to a thickness smaller than the limit, the bottom-most layer will be combined with the 2nd bottom layer and the snow parameters are averaged. Such fluctuations can also be found using the SNTHERM89 code. Secondly, the bottom snow mass can be reduced in SNTHERM via the water vapor movement, and this can result in a 'disappearance' of the bottom layer. After the grain growth coefficient (g1) was updated, the simulated snow depth and geometric grain size and density profiles are presented in Figure 8. In Figure 8a, after the g1 was improved, the simulated snow melts on 23:00 March 13th, 2018. Comparing with the measured snow-free day on March 13th (daily scale), the error is within one day. The new simulation has been significantly improved over the originally-simulated March 16th. It demonstrated that, after the grain size is accurately predicted to higher values, the snow albedo is reduced to promote a quicker snow melt. In Figure 8a, the snow depth prediction accuracy is −1.49 cm mean bias (MB) and 1.91 cm root-meansquared-error (RMSE). Figure 8b,e show the snow grain size and the snow density profiles. In Figure 8c,e, depthweighted average values and the values at the surface and the bottom of the snow are presented. The fluctuation of the surface layer geometric grain size and density is due to new snowfalls. After a snow layer is initialized from a new snowfall, in theory, the snow grain size monotonically increases. However, the combination of snow layers can reset the simulated values. For the bottom snow layer, two reasons from the model numerical implementation can result in a decreased snow grain size and density. Firstly, when the bottom layer is compressed to a thickness smaller than the limit, the bottom-most layer will be combined with the 2nd bottom layer and the snow parameters are averaged. Such fluctuations can also be found using the SNTHERM89 code. Secondly, the bottom snow mass can be reduced in SNTHERM via the water vapor movement, and this can result in a 'disappearance' of the bottom layer.  Figure 9 shows the comparison between the ground-based radiometer measured and the MEMLS simulated TB at 6.925, 10.65, 18.7, and 36.5 GHz, dual polarization, using the snowpit measurements (including layer thickness, geometric grain size, density, and temperature) and the soil parameters at −2 cm depth from the Decagon 5TM sensor. Figure 10 is the comparison result using all snow and soil parameters from the SNTHERM simulations. Here we utilized the Dmax to pex coefficients in Equation (8). The soil reflectivity was modeled using a fixed soil roughness of 2 cm for all cases.

Brightness Temperature Simulation for the Altay Experiment
An error statistic is shown in Table 2. It shows using the snowpit measurements, the root-mean square error (RMSE) of the brightness temperature at 18.7 V, 36.5 V, 18.7 H, and 36.5 H are 2.56 K, 5.54 K, 8.07 K, and 6.11 K, respectively. The mean bias (MB) is 1.86 K, −3.16 K, 7.51 K, and −3.51 K, respectively. Figure 10 shows the comparison using the SNTHERM simulated parameters based on both the original default g1 and the updated g1. It shows the updated g1 corrected the overestimated TB from the original g1. The difference is stronger at the end of the snow season when both the snow depth and the grain size is higher to become sensitive to the microwave volume scattering. It implies an important point that, because SNTHERM snow grain size monotonically increases with time, the simulation of snow particles needs to be very careful to prevent "over-growth" at the early stage when the TB is still not sensitive. Table 2 also shows the RMSE and MB using the updated g1 are of the same level with the simulation using the snowpit measurements.
As a summary, using the Altay experiment data, we managed to simultaneously match the observed snow depth, snow grain size and brightness temperature. Our validation reconfirmed the Dmax to pex conversion coefficients and suggested an increase of the grain growth coefficient by 2.088 times. The same parameter values will be used for simulations at 102 stations in the next section.  Figure 9 shows the comparison between the ground-based radiometer measured and the MEMLS simulated T B at 6.925, 10.65, 18.7, and 36.5 GHz, dual polarization, using the snowpit measurements (including layer thickness, geometric grain size, density, and temperature) and the soil parameters at −2 cm depth from the Decagon 5TM sensor. Figure 10 is the comparison result using all snow and soil parameters from the SNTHERM simulations. Here we utilized the D max to p ex coefficients in Equation (8). The soil reflectivity was modeled using a fixed soil roughness of 2 cm for all cases.

Brightness Temperature Simulation for the Altay Experiment
An error statistic is shown in Table 2. It shows using the snowpit measurements, the root-mean square error (RMSE) of the brightness temperature at 18.7 V, 36.5 V, 18.7 H, and 36.5 H are 2.56 K, 5.54 K, 8.07 K, and 6.11 K, respectively. The mean bias (MB) is 1.86 K, −3.16 K, 7.51 K, and −3.51 K, respectively. Figure 10 shows the comparison using the SNTHERM simulated parameters based on both the original default g 1 and the updated g 1 . It shows the updated g 1 corrected the overestimated T B from the original g 1 . The difference is stronger at the end of the snow season when both the snow depth and the grain size is higher to become sensitive to the microwave volume scattering. It implies an important point that, because SNTHERM snow grain size monotonically increases with time, the simulation of snow particles needs to be very careful to prevent "over-growth" at the early stage when the T B is still not sensitive. Table 2 also shows the RMSE and MB using the updated g 1 are of the same level with the simulation using the snowpit measurements.
As a summary, using the Altay experiment data, we managed to simultaneously match the observed snow depth, snow grain size and brightness temperature. Our validation reconfirmed the D max to p ex conversion coefficients and suggested an increase of the grain growth coefficient by 2.088 times. The same parameter values will be used for simulations at 102 stations in the next section. Remote Sens. 2019, 11, x FOR PEER REVIEW 18 of 29         Figure 13 additionally shows the statistics of daily snow depth estimation RMSE at each station. It shows that in most stations, the average snow depth is below 5 cm in the Tibet and the Qinghai provinces, but is much higher in the northern Xinjiang. The accuracy of the seasonal average snow depth is 1.27 cm RMSE and 0.98 correlation (R). The simulated snow depth is slightly biased low. The relative error is lower in northern Xinjiang where a deeper snowpack develops.  Figure 13 additionally shows the statistics of daily snow depth estimation RMSE at each station. It shows that in most stations, the average snow depth is below 5 cm in the Tibet and the Qinghai provinces, but is much higher in the northern Xinjiang. The accuracy of the seasonal average snow depth is 1.27 cm RMSE and 0.98 correlation (R). The simulated snow depth is slightly biased low. The relative error is lower in northern Xinjiang where a deeper snowpack develops. Remote Sens. 2019, 11, x FOR PEER REVIEW 21 of 29

Validation of the Brightness Temperature at the 102 Sites
Snow microstructure is the key parameter to influence the brightness temperature in snow. After the snow depth is validated, the SNTHERM simulated snow and soil parameters were inputted to MEMLS to calculate the pixel-scale brightness temperature and to compare it with the AMSR-E observations from the descending orbit (1:30 a.m. local overpass time, approximately). To reduce the influence of wet snow condition, the error statistics use the period from December to February in the next year. Figure 14 is an example of brightness temperature validation at the Altay station in the 2008-2009 winter. In this winter, the Altay station had a maximum snow depth of 60 cm. The TB at 10.65 GHz, 18.7 GHz, 23.8 GHz, 36.5 GHz, and 89 GHz, vertical polarization, was simulated with an overall RMSE of 8.8 K, from daily-based statistics. Figure 15 shows the comparison of measured and simulated daily TB at all stations. The agreement is higher for frequencies larger than 10.65 GHz. The RMSE is 6.91 K, 8.2 K, 13.34 K at 18.7, 23.8 and 36.5 GHz, vertical polarization, respectively. At 10.65 GHz, the correlation between the measured and the simulated TB is smaller, whereas the RMSE is 6.46 K and 12.93 K at vertical and horizontal polarization, respectively. The simulation at this frequency is influenced strongly by the soil modeling beneath the snow. We utilized the same soil parameters as in the Altay experiment, so the accuracy will be reduced. More important, usually the underlying background beneath the snow

Validation of the Brightness Temperature at the 102 Sites
Snow microstructure is the key parameter to influence the brightness temperature in snow. After the snow depth is validated, the SNTHERM simulated snow and soil parameters were inputted to MEMLS to calculate the pixel-scale brightness temperature and to compare it with the AMSR-E observations from the descending orbit (1:30 a.m. local overpass time, approximately). To reduce the influence of wet snow condition, the error statistics use the period from December to February in the next year. Figure 14 is an example of brightness temperature validation at the Altay station in the 2008-2009 winter. In this winter, the Altay station had a maximum snow depth of 60 cm. The TB at 10.65 GHz, 18.7 GHz, 23.8 GHz, 36.5 GHz, and 89 GHz, vertical polarization, was simulated with an overall RMSE of 8.8 K, from daily-based statistics. Figure 15 shows the comparison of measured and simulated daily TB at all stations. The agreement is higher for frequencies larger than 10.65 GHz. The RMSE is 6.91 K, 8.2 K, 13.34 K at 18.7, 23.8 and 36.5 GHz, vertical polarization, respectively. At 10.65 GHz, the correlation between the measured and the simulated TB is smaller, whereas the RMSE is 6.46 K and 12.93 K at vertical and horizontal polarization, respectively. The simulation at this frequency is influenced strongly by the soil modeling beneath the snow. We utilized the same soil parameters as in the Altay experiment, so the accuracy will be reduced. More important, usually the underlying background beneath the snow Figure 13. SNTHERM simulated and measured snow season average snow depth and the statistics of RMSE from daily snow depth for each station. The X-axis is the count number of the meteorological stations, from 1 to 102.

Validation of the Brightness Temperature at the 102 Sites
Snow microstructure is the key parameter to influence the brightness temperature in snow. After the snow depth is validated, the SNTHERM simulated snow and soil parameters were inputted to MEMLS to calculate the pixel-scale brightness temperature and to compare it with the AMSR-E observations from the descending orbit (1:30 a.m. local overpass time, approximately). To reduce the influence of wet snow condition, the error statistics use the period from December to February in the next year. Figure 14 is an example of brightness temperature validation at the Altay station in the 2008-2009 winter. In this winter, the Altay station had a maximum snow depth of 60 cm. The T B at 10.65 GHz, 18.7 GHz, 23.8 GHz, 36.5 GHz, and 89 GHz, vertical polarization, was simulated with an overall RMSE of 8.8 K, from daily-based statistics. Figure 15 shows the comparison of measured and simulated daily T B at all stations. The agreement is higher for frequencies larger than 10.65 GHz. The RMSE is 6.91 K, 8.2 K, 13.34 K at 18.7, 23.8 and 36.5 GHz, vertical polarization, respectively. At 10.65 GHz, the correlation between the measured and the simulated T B is smaller, whereas the RMSE is 6.46 K and 12.93 K at vertical and horizontal polarization, respectively. The simulation at this frequency is influenced strongly by the soil modeling beneath the snow. We utilized the same soil parameters as in the Altay experiment, so the accuracy will be reduced. More important, usually the underlying background beneath the snow cover in the pixel scale is different with the bare soil assumption. For Figure 15a,b, we did not correct the time independent offset sometimes found between the simulated and the observed T B time series caused by the soil modeling and the complex land surface conditions, which may be considered as a systematical bias like in Figures 9 and 10. The largest RMSE (of 24.95 K) is found at 89 GHz because this frequency is very sensitive to the grain size of the topmost snow layers. A good point is the simulated T B is unbiased in 160-200 K range. The T B above 200 K is biased higher, which indicates a possible underestimation of new snow grain size using the SNTHERM inner computation module instead of a user-defined value. The small 89 GHz T B measurements occurred when the snow crystals at the snow surface has grown large, before the next new snowfall.
Remote Sens. 2019, 11, x FOR PEER REVIEW 22 of 29 cover in the pixel scale is different with the bare soil assumption. For Figure 15a,b, we did not correct the time independent offset sometimes found between the simulated and the observed TB time series caused by the soil modeling and the complex land surface conditions, which may be considered as a systematical bias like in Figures 9 and 10. The largest RMSE (of 24.95 K) is found at 89 GHz because this frequency is very sensitive to the grain size of the topmost snow layers. A good point is the simulated TB is unbiased in 160-200 K range. The TB above 200 K is biased higher, which indicates a possible underestimation of new snow grain size using the SNTHERM inner computation module instead of a user-defined value. The small 89 GHz TB measurements occurred when the snow crystals at the snow surface has grown large, before the next new snowfall.  cover in the pixel scale is different with the bare soil assumption. For Figure 15a,b, we did not correct the time independent offset sometimes found between the simulated and the observed TB time series caused by the soil modeling and the complex land surface conditions, which may be considered as a systematical bias like in Figures 9 and 10. The largest RMSE (of 24.95 K) is found at 89 GHz because this frequency is very sensitive to the grain size of the topmost snow layers. A good point is the simulated TB is unbiased in 160-200 K range. The TB above 200 K is biased higher, which indicates a possible underestimation of new snow grain size using the SNTHERM inner computation module instead of a user-defined value. The small 89 GHz TB measurements occurred when the snow crystals at the snow surface has grown large, before the next new snowfall.

Analysis of the Influence of Soil Parameter Setting on the Brightness Temperature Simulation
Section 4.1 demonstrated that the setting of soil parameters influences the snow temperature, temperature gradient, and therefore influences the calculated grain size, according to the SNTHERM model theory. In Figures 4c and 6, it shows the use of a site-specific Altay soil type predicts the temperature variables better. In Figure 4f, the simulated bottom D max is in average slightly larger using the Altay site-specific soil (dark purple line) compared to clay (light blue line). The influence will be enlarged after the grain grow rate is increased. Therefore, here we added a discussion to evaluate the soil parameter influence on the T B simulation at 36.5 GHz. The T B is a more significant impact indicator compared to layered D max in small values. Figure 16 shows the comparison of T B using two different SNTHERM simulated soil properties (a) and the same property from 5TM soil sensor measurements (b). In Figure 16a, it shows surprisingly, the simulated T B is better using the original clay soil setting. The error is most likely from the soil dielectric constant modelling. By comparing with Figure 16b, it shows the reason is from the difference in the simulated soil properties caused by different soil parameter settings, not from the simulated D max difference in Figure 16c. In Figure 16c, it showed the average D max using Altay soil is still larger than the clay soil, which agrees with Figure 4f and the modeling theory; however, in the mid snow season, the bottom layer D max using clay soil can be reversely larger, due to the randomness in layer combination. V-pol., (d) 23.8 GHz V-pol., (e) 36.5 GHz V-pol., (f) 89 GHz V-pol. The color bar represents the number of points at each 1 K measured and simulated TB interval. The unit of the labeled root-mean-squared error (RMSE) and mean bias (MB) statistics is K.

Analysis of the Influence of Soil Parameter Setting on the Brightness Temperature Simulation
Section 4.1 demonstrated that the setting of soil parameters influences the snow temperature, temperature gradient, and therefore influences the calculated grain size, according to the SNTHERM model theory. In Figures 4c and 6, it shows the use of a site-specific Altay soil type predicts the temperature variables better. In Figure 4f, the simulated bottom Dmax is in average slightly larger using the Altay site-specific soil (dark purple line) compared to clay (light blue line). The influence will be enlarged after the grain grow rate is increased. Therefore, here we added a discussion to evaluate the soil parameter influence on the TB simulation at 36.5 GHz. The TB is a more significant impact indicator compared to layered Dmax in small values. Figure 16 shows the comparison of TB using two different SNTHERM simulated soil properties (a) and the same property from 5TM soil sensor measurements (b). In Figure 16a, it shows surprisingly, the simulated TB is better using the original clay soil setting. The error is most likely from the soil dielectric constant modelling. By comparing with Figure 16b, it shows the reason is from the difference in the simulated soil properties caused by different soil parameter settings, not from the simulated Dmax difference in Figure 16c. In Figure 16c, it showed the average Dmax using Altay soil is still larger than the clay soil, which agrees with Figure 4f and the modeling theory; however, in the mid snow season, the bottom layer Dmax using clay soil can be reversely larger, due to the randomness in layer combination.  Table 3 compares the TB difference caused by different soil parameter settings, by the influence of soil parameter settings on the simulated snow grain size (with soil property difference excluded) and by the different grain growth coefficients (the original and the revised g1). It shows the influence of grain growth coefficient correction is strongest (of 6.93 K in average), indicating its importance.  Table 3 compares the T B difference caused by different soil parameter settings, by the influence of soil parameter settings on the simulated snow grain size (with soil property difference excluded) and by the different grain growth coefficients (the original and the revised g 1 ). It shows the influence of grain growth coefficient correction is strongest (of 6.93 K in average), indicating its importance. The influence of soil property simulation accuracy is in the second order, with a magnitude of 5.18 K in average, but its maximum influence is within 10 K. The influence of soil parameter setting on the geometric grain size prediction is in the order of 0.1 mm in average, and the resulted T B difference at 36.5 GHz, vertical polarization is within 1 K in average. It implies the SNTHERM simulated snow properties are trustable regardless of the soil parameter setting.

Evaluation the Use of GLDAS Downwelling Longwave Radiation
In Section 2.1, we mentioned the downwelling longwave radiation inputted to the SNTHERM model is from the GLDAS dataset. The station does not have the downwelling longwave radiation measurement. We didn't force the simulated net radiation to match the observed net radiation, because we considered the difference between the SNTHERM model pure snow assumption and the in-situ condition influenced by the dry high grass. A natural condition is required by the meteorological station observation; however, as seen from Figure 17, the measured net radiation has a smaller peak at noon, influenced by the low albedo and the shade from the dry high grass.
Averaged from the November 10th, 2017 to March 12th, 2018 snow-covered period, the SNTHERM simulated net radiation on pure snow is 9.84 W/m 2 higher compared to the measurements from the meteorological site; the hourly standard deviation of difference is 31.15 W/m 2 . The influence of soil property simulation accuracy is in the second order, with a magnitude of 5.18 K in average, but its maximum influence is within 10 K. The influence of soil parameter setting on the geometric grain size prediction is in the order of 0.1 mm in average, and the resulted TB difference at 36.5 GHz, vertical polarization is within 1 K in average. It implies the SNTHERM simulated snow properties are trustable regardless of the soil parameter setting.

Evaluation the Use of GLDAS Downwelling Longwave Radiation
In Section 2.1, we mentioned the downwelling longwave radiation inputted to the SNTHERM model is from the GLDAS dataset. The station does not have the downwelling longwave radiation measurement. We didn't force the simulated net radiation to match the observed net radiation, because we considered the difference between the SNTHERM model pure snow assumption and the in-situ condition influenced by the dry high grass. A natural condition is required by the meteorological station observation; however, as seen from Figure 17, the measured net radiation has a smaller peak at noon, influenced by the low albedo and the shade from the dry high grass.
Averaged from the November 10th, 2017 to March 12th, 2018 snow-covered period, the SNTHERM simulated net radiation on pure snow is 9.84 W/m 2 higher compared to the measurements from the meteorological site; the hourly standard deviation of difference is 31.15 W/m 2 .

Error Analysis for the Brightness Temperature Simulation at 102 Stations
The accuracy of the brightness temperature simulation at the 102 stations is influenced by the following reasons. Firstly, the snow grain size was not directly validated on the Tibetan Plateau, its prediction accuracy still needs more exploration. Secondly, the soil hydraulic and thermal parameters

Error Analysis for the Brightness Temperature Simulation at 102 Stations
The accuracy of the brightness temperature simulation at the 102 stations is influenced by the following reasons. Firstly, the snow grain size was not directly validated on the Tibetan Plateau, its prediction accuracy still needs more exploration. Secondly, the soil hydraulic and thermal parameters were not set using local datasets. Although we have shown that the influence of soil setting on snow parameter estimation is limited, the soil parameters can still change the background emissivity beneath the snow and influence the T B simulation; it was about 10 K in Altay experiment. Thirdly, the background condition of each pixel, which is in most cases mixed pixel at the AMSR-E resolution, is much complicated compared with the SNTHERM assumption. We had partly considered its influence by introducing the background emissivity and the fractional snow cover inputs. However, several factors were still not included, like the forest over the snow, the accelerated snowmelt from the warm grass under the snow, the different freeze-thaw evolution process of other land cover types in the mixed pixel, terrain effect on microwave modeling, terrain effect on snowdrift by winds, error in the satellite passive microwave observations, etc. These factors are the potential error sources for our simulation result in Section 4.5.

Conclusions
In this paper, we validated the use of the SNTHERM combined with MEMLS to simulate the snow depth, snow grain size, and brightness temperature in the Xinjiang, Qinghai, and Tibet provinces in China. Results showed that the SNTHERM has a good daily snow depth estimation accuracy, in the order of l-4 cm, when the accurate meteorological forcing is provided. During the Altay experiment, comparison with the in-situ measurements shows that the snow grain size can be predicted by an RMSE of 0.47 mm and mean bias (MB) of −0.15 mm. The simulated daily brightness temperature at 36.5 GHz, vertical polarization, has an accuracy of 4.43 K RMSE using the coupled SNTHERM and MEMLS model at the point-scale by comparing with the ground-based radiometer measurements from the Altay field experiment and 13.34 K using the AMSR-E observed T B from the 102 meteorological stations of varied snow, soil property and mixed-pixel background characteristics. The statistics is based on the dry snow period from December to February.
Using the field experiment of the entire winter, we updated the grain growth coefficient (g 1 ) of the SNTHERM, and suggested a new value of 1.044 × 10 −6 m 4 /kg. After this revision, we linked the SNTHERM predicted snow grain size to the geometric grain size definition that can be observed in the field.
The snow is very different on the Qinghai-Tibetan Plateau (QTP), compared to the northern Xinjiang. Driven by high solar radiation, the snow is easier to melt on the QTP and the snow depth is shallower in flat areas. The brightness temperature for shallow snow is easier to be influenced by other nuisance parameters, especially the background emission and the signals from other land cover types. Therefore, an extended field experiment is suggested in this region to validate the geometric grain size directly. The SNTHERM snow grain size was found to be biased in previous studies [34], whereas in this paper, this bias was corrected before the snow grain size was inputted to the snow brightness temperature modelling. Used together with the snow temperature and the snow temperature gradient, the grain growth rate g 1 of SNTHERM reduces its dependency on the other snow variables and has the potential to become a relatively stable and consistent parameter. Our study also found the correction of snow grain size can improve the accuracy in End of Snow Season Day (ESSD) by 3 days. We are looking forward to more studies and validations for g 1 , as well as the similar coefficient in other physical snow grain growth models. It will enhance the physic basis in the snow parameter modeling and the microwave signal modeling supported by the snow process simulations.  are grateful to the National Meteorological Information Center of the China Meteorological Administration (CMA) for providing the daily meteorological observations at the stations, and the scientific teams and the NASA GES DISC (Goddard Earth Sciences Data and Information Services Center) and LP DAAC (Land Processes Distributed Active Archive Center) for producing and distributing the GLDAS and MODIS reflectance datasets.

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

Appendix A. The Minimum Bug Fixing Required for a Normal Operation of the SNTHERM89 Code
As found from our simulation test using the Altay experiment meteorological forcing data, three lines in the SNTHERM89 code are need to be corrected, with contents and reasons listed as follows.
(1) In sntherm.f, line 475, 'bw(i) = wmass/dz(i)': it calculates the snow water density (bw(i)) as the snow mass (wmass) divided by snow layer thickness (dz(i)), after dz(i) was updated using the snow compaction rate. It is suggested to add an upper limit on the computed density as: if(bw(i).gt.917d0)then bw(i) = 917d0 dz(i) = wmass/bw(i) endif Otherwise, an abnormal snow density may be predicted if the snow compaction is not in phase with the snow mass change. This error usually occurs for shallow or melting snow.
(3) In gemet.f, line 111, 'if (ibasestep .eq. 1 .or. ihour .eq. 4) ido = 1' is better revised as 'ido=1'. The previous code calculates the solar zenith angle only at 4:00 a.m. It will wrongly give a zero solar zenith angle for other time steps and jump the shallow snow correction module for the albedo calculation. The original code may result in a slow snowmelt because the snow albedo is overestimated.
Additionally, the activation of the shallow snow correction module requires inputting a user-defined constant snow albedo larger than 1.0 (to defunction its use) and inputting the latitude, longitude, slope, aspect and offset of the meteorological forcing data time compared to the UTC time. Important reminder: the longitude and timezone offset in the SNTHEM89 code are positive for the West Hemisphere. For example. 'dlongt' and 'itimezone' in getinput.f expect inputs as −88.0, −8 for 88 • E, UTC+8 local time offset, respectively.