A Novel Approach to Validate Satellite Snowfall Retrievals by Ground-Based Point Measurements

: A novel method has been proposed for validating satellite radar snowfall retrievals using surface station observations over the western United States mountainous region, where the mean snowfall rate at a station depends on its elevation. First, all station data within a 1 ◦ × 1 ◦ grid are used to develop a snowfall rate versus elevation relation. This relation is then used to compute snowfall rate in other locations within the 1 ◦ × 1 ◦ grid, as if surface observations were available everywhere in the grid. Grid mean snowfall rates are then derived, which should be more representative to the mean snowfall rate of the grid than using data at any one station or from a simple mean of all stations in the grid. Comparison of the so-derived grid mean snowfall rates with CloudSat retrievals shows that the CloudSat product underestimates snowfall by about 65% when averaged over all the 768 grids in the western United States mountainous regions. The bias does not seem to have clear dependency on elevation but strongly depends on snowfall rate. As an application of the method, we further estimated the snowfall to precipitation ratio using both ground and satellite measured data. It is found that the rates of increase with elevation of the snowfall to precipitation ratio are quite similar when calculating from ground and satellite data, being about 25% per kilometer elevation up or approximately 4% per every degree Cuisses of temperature drop.


Introduction
Glacier and snowpack over high mountains are critical water resources to populations living in the lowland regions. For example, upland areas (above 2000 m elevation) in the high mountains of Asia supply the five basins of the Indus, Ganges, Yellow, Brahmaputra, and Yangtze rivers, providing water to 1.4 billion people in the downstream region [1,2]. In the western United States (U.S.), the bulk of surface water resources, as represented by the flow of the Colorado and Columbia River systems, is derived from melted winter snowpack [3]. In a warming climate, glacier is depleting over a global scale [4,5]. In the western U.S., it is observed that mountain snowpack is declining accompanied with rising surface temperature and early start of snow melting in the recent decades [6], which threatens the water resources in the western states [7]. The long-term variation of snowfall is an important indicator of climate change in both regional and global scales [8,9] and a key component of the hydrological cycle in the mid-and high latitudes [10][11][12][13].
For estimation of snowfall over a large-scale, satellite remote sensing becomes inevitable, particularly over remote mountainous regions. There are currently several satellite precipitation products available, for example, radar products [14,15], microwave radiometer products [16], radar-radiometer combined products [17], and multisensor merged products [18][19][20]. Since passive sensors measure the combined contribution from the atmosphere and surface to the upwelling radiation, the quality of their snowfall retrieval suffers greatly from the complexity of snow-covered surfaces, particularly over mountainous terrains [10,21]. Satellite radars are considered to be the most suitable sensors for

Data
In this section, we describe the datasets used in this study, which include surface station data from SNOTEL and GHCND and satellite radar retrievals from CloudSat CPR and GPM DPR.

Surface Station Data
The SNOTEL network is designed by the United States Department of Agriculture to collect snowpack and related climate data in the western U.S. and Alaska. It provides fully automated data for high snow accumulation regions including many remote areas. Snow water equivalent measurements are made using snow pillows filled with an antifreeze solution. As snow accumulates, a pressure transducer monitors the pressure of the fluid and converts the pressure to snow water equivalent [3]. In addition, snow depth, all-season precipitation accumulation, and air temperature with daily maximums, minimums, and averages are available. Although the snow pillow provides hourly data, we only use daily snow measurements in this study. The daily snow accumulation is calculated by the difference in cumulative values between two consecutive days, and a negative snow accumulation value is screened out [41]. In this study, data from 768 SNOTEL stations in the western U.S. were used to compute annual mean snowfall rate covering the period of 2006 through 2017. In addition, total precipitation measured by gauges is also reported at each station. We computed rainfall using daily total precipitation minus snowfall for the purpose of deriving snowfall to precipitation ratio.
The GHCND dataset has provided daily climate data including maximum and minimum temperature, total daily precipitation, snowfall, and snow depth etc. over global land areas [40], although we only use data over the western U.S. in this study. The U.S. collection contains the most complete daily data including some of the 19th century observations [43] as well as the 21st century measurements from the U.S. Climate Reference Network. We converted the variable "snowfall" (daily snowfall depth) to snow water equivalent to make it comparable to snowfall quantity in other datasets by assuming snowpack density of 0.1 g cm −3 [44]. Only data that passed all quality assurance check as indicated in the dataset are included in the data analysis. The period covered in this study for annual mean snowfall rate calculation is from 2006 through 2017. In addition to snowfall, total precipitation reported at the GHCND stations was used for the computation of snowfall to precipitation ratio.

Satellite Data
The satellite retrievals to be validated are snowfall from CloudSat CPR. The CPR operates at W band with the minimum detectability of around −30 dBZ. Its high sensitivity to light snow and broad global coverage makes it a good candidate for snowfall retrievals. The radar reflectivity is sampled in the vertical with a bin size of 240 m. The footprint size of radar reflectivity profiles is 1.4 km cross track by 2.5 km along track. The standard snowfall product, 2C-SNOW-PROFILE (Version R05), was created from CPR radar reflectivity profiles based on [15].
To estimate the snowfall to total precipitation ratio, GPM DPR retrievals are also used for the estimation of rainfall. The DPR operates at Ku and Ka bands. The Ku radar swath has a width of 245 km and consists of 49 beams with a vertical resolution of 250 m in the normal scan mode. In this study, the normal scan mode data from the DPR level 2A precipitation retrievals at version 6 [45] are used. Since the minimum detectability of DPR is around 13 dBZ, it misses a large portion of light precipitation events [24]. Therefore, we used DPR only for rainfall estimates.
Though both CPR and DPR products have their intrinsic flag for precipitation type, the classification method they use is different. To conduct consistent phase determination across the CPR and DPR datasets, we applied the phase classification scheme developed by [46]. Meteorological variables collocated from the fifth generation European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis (ERA5, [47]) hourly reanalysis data at a resolution of 0.25 • (latitude) × 0.25 • (longitude) were used as inputs, including the 2-m temperature and dew point, surface pressure, and low-level lapse rate.
The key input variable into the phase classification scheme is the 2-m temperature. According to [48], the 30-km spatial resolution of ERA5 is inadequate to provide accurate 2-m temperature at the radar pixel scale of a few kilometers over the complex terrains. To mitigate this problem, we followed the topography correction method proposed by [48] to derive the 2-m temperature at the radar footprint. The Shuttle Radar Topography Mission dataset with 15 arcseconds (~417 m) resolution, SRTM15 [49], was used as topography reference. We first calculated the average elevation from SRTM15 at each ERA5 grid, and then computed the difference between the elevation at the radar pixels and corresponding ERA5 grids. Using the low-level lapse rate, we interpolated the 2-m temperature at surface level. The temperature difference before and after correction ranges from −10 to 10 K, which would substantially impact the precipitation phase determination. We used the 2-m temperature, relative humidity, and low-level lapse rate as inputs to the phase classification scheme of [46], which provides the conditional probability of solid precipitation. A probability of 50% was regarded as the threshold to separate rain and snow. Once the precipitation is classified, the snowfall rate in CPR and the rainfall rate in DPR products were selected for analysis.
The CloudSat snow product ranging from 2006 to 2017 was used for annual mean snowfall estimates. The snowfall near the surface cannot be reliably measured due to contamination by the ground clutters. So, after we decided the precipitation phase at surface to be solid, we used the lowest available snow retrieval as the surface snowfall rate. For surface rainfall estimates, we used the level-2A product, 2ADPR, "DPR precipitation" for the surface rainfall retrievals [45] covering 2015 to 2019 to compute the annual mean rainfall rate. Similar to CPR snowfall estimation, we first decided on the surface precipitation type at the time and location of the DPR radar pixel based on the phase classification algorithm [46]. Once rain was recognized, we used the precipitation retrieval at the lowest uncontaminated level as the surface rainfall estimate.
The precipitation observations were then sorted into 1 • × 1 • grids. The annual mean value is the total snowfall or rainfall divided by the total number of observations within each grid. To avoid uneven sampling among seasons, we averaged the rainfall or snowfall rates in every 5-day window for all available years, then computed the annual mean. The annual mean snowfall to total precipitation ratio is calculated as the annual mean snowfall rate divided by the annual mean precipitation rate, which is snowfall plus rainfall rates.

A Novel Approach for Satellite Snowfall Validation
As mentioned earlier, to obtain a statistically meaningful comparison between satellite radar and surface station measurements, we need to average data over a sufficiently long time and over a large enough area surrounding the surface station. In several tests of try and error, we found that we need to average the CloudSat data over an area of 1 • (latitude) by 1 • (longitude) for the entire observation record from mid-2006 to mid-2017 for the comparison statistics to be stable. Therefore, in the following discussions in this section we will use CloudSat data averaged over 1 • × 1 • for 2006 through 2017 to compare station mean data averaged during 2006-2017.
We first demonstrate the problem arising from using either SNOTEL or GHCND data alone to validate CloudSat snowfall retrievals. Figure 1 shows the annual means of (a) SNOTEL and (b) GHCND snowfall rates and the CloudSat minus (c) SNOTEL and (d) GHCND annual mean snowfall rates. The mean snowfall rates are generally larger at SNOTEL than at GHCND stations. As a result, CloudSat CPR mostly underestimates snowfall against SNOTEL while having a mixture of over-and underestimation against Remote Sens. 2022, 14, 434 5 of 17 GHCND observations. One may suspect that this systematic difference is due to the inaccuracy of either SNOTEL or GHCND or both snowfall estimates. However, we argue in the following that it is largely caused by the elevation difference of these two types of stations.
We first demonstrate the problem arising from using either SNOTEL or GHCND data alone to validate CloudSat snowfall retrievals. Figure 1 shows the annual means of (a) SNOTEL and (b) GHCND snowfall rates and the CloudSat minus (c) SNOTEL and (d) GHCND annual mean snowfall rates. The mean snowfall rates are generally larger at SNOTEL than at GHCND stations. As a result, CloudSat CPR mostly underestimates snowfall against SNOTEL while having a mixture of over-and underestimation against GHCND observations. One may suspect that this systematic difference is due to the inaccuracy of either SNOTEL or GHCND or both snowfall estimates. However, we argue in the following that it is largely caused by the elevation difference of these two types of stations. Since SNOTEL collects data automatically, the stations can be placed in remote areas where frequent human access is not needed. However, GHCND stations are mostly conventional stations where routine manual measurements are required, and they are mostly located in the valleys for easy human access. This systematic difference of station elevations is shown in Figure 2. If snowfall intensity depends on station elevation, a systematic difference in annual mean snowfall between the two types of stations will occur. To test this hypothesis, we plot annual mean snowfall rates versus station elevation in Figure 3 using both SNOTEL and GHCND data. A clear trend of snowfall rate increasing with elevation can be observed, and the SNOTEL and GHCND data seem to follow the same increasing trend. Therefore, we argue that the systematic difference between SNOTEL and Since SNOTEL collects data automatically, the stations can be placed in remote areas where frequent human access is not needed. However, GHCND stations are mostly conventional stations where routine manual measurements are required, and they are mostly located in the valleys for easy human access. This systematic difference of station elevations is shown in Figure 2. If snowfall intensity depends on station elevation, a systematic difference in annual mean snowfall between the two types of stations will occur. To test this hypothesis, we plot annual mean snowfall rates versus station elevation in Figure 3 using both SNOTEL and GHCND data. A clear trend of snowfall rate increasing with elevation can be observed, and the SNOTEL and GHCND data seem to follow the same increasing trend. Therefore, we argue that the systematic difference between SNOTEL and GHCND mean snowfall rates is not a result of observation error at the stations but rather caused by the systematic difference of station elevations. To further verify this hypothesis, we calculated the difference of mean snowfall rates for all SNOTEL-GHCND station pairs that are within 25 km horizontal distance and plotted the results in Figure 4 as a function of the elevation difference between the pair and SNOTEL station elevation. It is shown that the snowfall rate difference is near zero when the pair of stations are at the same elevation. The value of SNOTEL minus GHCND snowfall rates even becomes negative when the SNOTEL station is at a lower elevation than the GHCND station; the large positive values occur when the SNOTEL station is higher than the GHCND station.
GHCND mean snowfall rates is not a result of observation error at the stations but rather caused by the systematic difference of station elevations. To further verify this hypothesis, we calculated the difference of mean snowfall rates for all SNOTEL-GHCND station pairs that are within 25 km horizontal distance and plotted the results in Figure 4 as a function of the elevation difference between the pair and SNOTEL station elevation. It is shown that the snowfall rate difference is near zero when the pair of stations are at the same elevation. The value of SNOTEL minus GHCND snowfall rates even becomes negative when the SNOTEL station is at a lower elevation than the GHCND station; the large positive values occur when the SNOTEL station is higher than the GHCND station.   GHCND mean snowfall rates is not a result of observation error at the stations but rather caused by the systematic difference of station elevations. To further verify this hypothesis, we calculated the difference of mean snowfall rates for all SNOTEL-GHCND station pairs that are within 25 km horizontal distance and plotted the results in Figure 4 as a function of the elevation difference between the pair and SNOTEL station elevation. It is shown that the snowfall rate difference is near zero when the pair of stations are at the same elevation. The value of SNOTEL minus GHCND snowfall rates even becomes negative when the SNOTEL station is at a lower elevation than the GHCND station; the large positive values occur when the SNOTEL station is higher than the GHCND station.   To account for this snowfall dependency on elevation, a new validation dataset that combines the SNOTEL and GHCND data is created using the method as schematically shown in Figure 5. This is one of the 768 1 • × 1 • grids centered at SNOTEL stations examined in this study. First, we divide each 1 • × 1 • grid centered at a SNOTEL station into 240 × 240 sub-grids (Note that the number of sub-grids shown in Figure 5a is much less than 240 × 240 for figure clarity); the mean elevation of each sub-grid is determined by SRTM15 topography data. The relation between annual mean snowfall rate and elevation is derived from combined SNOTEL and GHCND data as shown in Figure 5b. This relation is used to fill the values at those sub-grids where there is neither SNOTEL nor GHCND stations. The 1 • × 1 • grid mean snowfall rate is then computed by averaging over values of all the sub-grids. This areal mean value will be used to compare with CloudSat snowfall observations within the 1 • × 1 • area centered at the SNOTEL location. To account for this snowfall dependency on elevation, a new validation dataset that combines the SNOTEL and GHCND data is created using the method as schematically shown in Figure 5. This is one of the 768 1° × 1° grids centered at SNOTEL stations examined in this study. First, we divide each 1° × 1° grid centered at a SNOTEL station into 240 × 240 sub-grids (Note that the number of sub-grids shown in Figure 5a is much less than 240 × 240 for figure clarity); the mean elevation of each sub-grid is determined by SRTM15 topography data. The relation between annual mean snowfall rate and elevation is derived from combined SNOTEL and GHCND data as shown in Figure 5b. This relation is used to fill the values at those sub-grids where there is neither SNOTEL nor GHCND stations. The 1° × 1° grid mean snowfall rate is then computed by averaging over values of all the sub-grids. This areal mean value will be used to compare with CloudSat snowfall observations within the 1° × 1° area centered at the SNOTEL location.
A key procedure in this method is to determine the mean snowfall rate versus elevation relation, which varies depending on locations. We carefully examined this relation for all the grids corresponding to the 768 SNOTEL stations; each grid was assigned a unique relation although some of them are quite similar. For most of the grids, we found that an exponential function fits the relation well. In Figure 6   The exponential fitting is applied with the fitting function shown in the legend. The 1° × 1° area is divided into 240 × 240 sub-grids. Note that the number of sub-grids shown in (a) is much less than 240 × 240 for figure clarity. The derived snowfall rate versus elevation relation will be used to compute snowfall rate at sub-grid without a SNOTEL or GHCND station, shown as black empty circles in (a). The exponential fitting is applied with the fitting function shown in the legend. The 1 • × 1 • area is divided into 240 × 240 sub-grids. Note that the number of sub-grids shown in (a) is much less than 240 × 240 for figure clarity. The derived snowfall rate versus elevation relation will be used to compute snowfall rate at sub-grid without a SNOTEL or GHCND station, shown as black empty circles in (a). A key procedure in this method is to determine the mean snowfall rate versus elevation relation, which varies depending on locations. We carefully examined this relation for all the grids corresponding to the 768 SNOTEL stations; each grid was assigned a unique relation although some of them are quite similar. For most of the grids, we found that an exponential function fits the relation well. In Figure 6 we show several examples illustrating the typical shapes of the snowfall-elevation relation. A fitting curve (sometimes stepwised for a better fitting) is also shown in each figure, which is used to compute the snowfall rate in the sub-grids that are not occupied by either SNOTEL or GHCND stations.   Essentially, this method is to mimic the grid mean as if ground stations are densely positioned in every sub-grid. It is noted that the gap-filling operation is not archived by spatial interpolation but rather by relying on the elevation dependence of climatological mean of snowfall. This elevation dependence seems to be a characteristic of snowfall precipitation as it is also reported over the Tibetan Plateau [48]. The mean snowfall rate so derived for a 1° × 1° grid should be more representative of that for the area measured by CloudSat, thus validation against this value is more meaningful. In the following discussions, this new value is referred to as "grid mean" snowfall rate. Essentially, this method is to mimic the grid mean as if ground stations are densely positioned in every sub-grid. It is noted that the gap-filling operation is not archived by spatial interpolation but rather by relying on the elevation dependence of climatological mean of snowfall. This elevation dependence seems to be a characteristic of snowfall precipitation as it is also reported over the Tibetan Plateau [48]. The mean snowfall rate so derived for a 1 • × 1 • grid should be more representative of that for the area measured by CloudSat, thus validation against this value is more meaningful. In the following discussions, this new value is referred to as "grid mean" snowfall rate.
In Figure 7 we show how much difference the above method makes in annual mean snowfall rates against original SNOTEL values. Overall, the grid mean snowfall rate is lowered by 0.43 mm day −1 (or 25%), averaged over all stations in the region. However, the sign of the changes is not uniformly distributed as it generally adjusts the value downward where the SNOTEL's elevation is higher than its surrounding areas and vice versa.
In Figure 7 we show how much difference the above method makes in annual mean snowfall rates against original SNOTEL values. Overall, the grid mean snowfall rate is lowered by 0.43 mm day −1 (or 25%), averaged over all stations in the region. However, the sign of the changes is not uniformly distributed as it generally adjusts the value downward where the SNOTEL's elevation is higher than its surrounding areas and vice versa. CloudSat CPR retrievals are then compared with the grid mean snowfall rates, and the results are shown in Figures 8 and 9. Recall that the CloudSat annual mean values are derived by averaging data from 2006 through 2017 over a 1° × 1° area centered at SNOTEL stations. In Figure 8, the comparison indicates that CloudSat underestimates snowfall for most of the grids in the west U.S. mountainous region. In Figure 9, the scatterplots show the difference and ratio between CloudSat retrieval and grid mean snowfall rate over a 1° × 1° area centered at SNOTEL station as a function of grid mean elevation or snowfall rates. It is seen that the bias or ratio does not seem to have a dependency on elevation. However, they strongly depend on snowfall rate-the higher the snowfall rate is, the larger the underestimation by CloudSat will be. When averaged over the grids around all the 768 SNOTEL stations, the bias of CloudSat estimates is about −0.85 mm day −1 , which is −65% with respect to the grid mean snowfall rate estimated by the ground measurements. While further studies are needed, we suspect that this underestimation is largely a result of the space radar's inability to obtain valid measurements close to the surface. The lowest level without ground contamination is often 1 km or higher above the ground, and snowfall rate often decreases with height in this region. CloudSat CPR retrievals are then compared with the grid mean snowfall rates, and the results are shown in Figures 8 and 9. Recall that the CloudSat annual mean values are derived by averaging data from 2006 through 2017 over a 1 • × 1 • area centered at SNOTEL stations. In Figure 8, the comparison indicates that CloudSat underestimates snowfall for most of the grids in the west U.S. mountainous region. In Figure 9, the scatterplots show the difference and ratio between CloudSat retrieval and grid mean snowfall rate over a 1 • × 1 • area centered at SNOTEL station as a function of grid mean elevation or snowfall rates. It is seen that the bias or ratio does not seem to have a dependency on elevation. However, they strongly depend on snowfall rate-the higher the snowfall rate is, the larger the underestimation by CloudSat will be. When averaged over the grids around all the 768 SNOTEL stations, the bias of CloudSat estimates is about −0.85 mm day −1 , which is −65% with respect to the grid mean snowfall rate estimated by the ground measurements. While further studies are needed, we suspect that this underestimation is largely a result of the space radar's inability to obtain valid measurements close to the surface. The lowest level without ground contamination is often 1 km or higher above the ground, and snowfall rate often decreases with height in this region.

Investigation of the Snowfall to Precipitation Ratio
The change in snowfall to total precipitation ratio can serve as a useful indicator on how global warming impacts on hydrological cycle. In a recent study using present weather reports, Shi and Liu found that the ratio of snowfall to precipitation occurrence has been decreasing globally in the past 40 years while at high latitudes shows an increasing trend [50]. The ratio of snowfall to precipitation amount cannot be examined in their study because of the lack of snowfall data. Over the contiguous United States, Feng and Hu studied the ratio of snowfall to precipitation amount using the U.S. Historical Climatology Network (USHCN) data and found a decreasing trend over the past five decades in most of the regions [51]. However, similar to the GHCND data, the USHCN data used in their study are skewed toward the lowland areas where snowfall amounts are generally lower than those over high mountains. As an application of the new validation dataset created in this study, in this section we study the snowfall to precipitation ratio (hereafter, referred to as S/P ratio) over the west U.S. mountainous region and examine how well this ratio can be estimated from satellite radar observations. Remote Sens. 2022, 14, x FOR PEER REVIEW 10 of 18 Figure 8. The difference between CloudSat (SCloudSat) and 1° × 1° grid mean (SGrid_mean) snowfall rates. Figure 9. Scatterplots of the difference (a,b) and the ratio (c,d) between CloudSat retrieval (SCloudSat) and grid mean snowfall rates (SGrid_mean) over 1° × 1° grids centered at SNOTEL locations.   To calculate S/P ratio, we first examine whether rainfall rates also have elevation dependence in this region. Similar to Figures 3 and 4, in Figure 10a,b we show the scatterplot of annual mean rainfall rate observed at SNOTEL and GHCND stations versus elevation and the difference of snowfall rates between SNOTEL and GHCND station pairs that are within 25 km as a function of their station's elevation difference. Unlike the case of snowfall, the annual mean rainfall rate does not show a clear dependency on elevation. We also examined SNOTEL and GHCND rainfall data in the 1 • × 1 • area surrounding the 768 SNOTEL stations. No clear elevation dependency is found either. Based on the above findings, we calculated the ground truth 1 • × 1 • grid mean rainfall simply by averaging all ground observations within the grid. The S/P ratio based on ground data was then computed by the grid mean snowfall rate and the average rainfall rate over a 1 • × 1 • area surrounding SNOTEL stations.
tology Network (USHCN) data and found a decreasing trend over the past five decades in most of the regions [51]. However, similar to the GHCND data, the USHCN data used in their study are skewed toward the lowland areas where snowfall amounts are generally lower than those over high mountains. As an application of the new validation dataset created in this study, in this section we study the snowfall to precipitation ratio (hereafter, referred to as S/P ratio) over the west U.S. mountainous region and examine how well this ratio can be estimated from satellite radar observations.
To calculate S/P ratio, we first examine whether rainfall rates also have elevation dependence in this region. Similar to Figures 3 and 4, in Figure 10a,b we show the scatterplot of annual mean rainfall rate observed at SNOTEL and GHCND stations versus elevation and the difference of snowfall rates between SNOTEL and GHCND station pairs that are within 25 km as a function of their station's elevation difference. Unlike the case of snowfall, the annual mean rainfall rate does not show a clear dependency on elevation. We also examined SNOTEL and GHCND rainfall data in the 1° × 1° area surrounding the 768 SNOTEL stations. No clear elevation dependency is found either. Based on the above findings, we calculated the ground truth 1° × 1° grid mean rainfall simply by averaging all ground observations within the grid. The S/P ratio based on ground data was then computed by the grid mean snowfall rate and the average rainfall rate over a 1° × 1° area surrounding SNOTEL stations. The regional distribution and elevation dependency of S/P ratio estimated by ground data are shown in Figure 11. S/P ratio is larger than 20% for most of the grids and there is a clear positive correlation between elevation and S/P ratio. Roughly, S/P ratio increases 25% with an increase of 1 km elevation. If the mean temperature lapse rate is 6.5 °C km −1 , it translates to a roughly 4% S/P ratio change per every degree Celsius temperature The regional distribution and elevation dependency of S/P ratio estimated by ground data are shown in Figure 11. S/P ratio is larger than 20% for most of the grids and there is a clear positive correlation between elevation and S/P ratio. Roughly, S/P ratio increases 25% with an increase of 1 km elevation. If the mean temperature lapse rate is 6.5 • C km −1 , it translates to a roughly 4% S/P ratio change per every degree Celsius temperature change. Albeit rather crudely, this may be considered to be an observation-based estimate of the sensitivity of precipitation phase change to a warming environment.
Similar to the approach proposed by [48], we may estimate the S/P ratio based on the satellite radar data using CloudSat snowfall averaged from 2006 to 2017 and DPR rainfall averaged from 2014 to 2020. The S/P ratio maps estimated by satellite data and its dependency on elevation are shown in Figure 12. Compared to the results shown in Figure 11, although the satellite radar data derived elevation dependency of the S/P ratio is noisier in the scatterplot, the relation seems to follow a similar trend, i.e., about 20~25% increase in S/P ratio per every kilometer increase in elevation. The noisy relation may be improved by averaging data over a large area, for example, 1 • × 2 • as done by [48], which is beyond the scope of this study and will be examined in the future. The S/P ratio estimated from the satellite is somewhat lower (by 20% on average) than that estimated by ground data. Recall that CloudSat underestimates snowfall in this region by about 65% against ground measured grid mean values. The lower S/P ratio by satellite data is partially caused by the CloudSat snowfall underestimation. However, based on our evaluation, the GPM DPR retrievals also underestimate rainfall in this region by about 50%, possibly because the GPM rainfall, too, is derived from DPR reflectivity at levels of 1 km or above over actual surface. As a result, the S/P ratio underestimation becomes less severe than either the CloudSat snowfall or the GPM DPR rainfall underestimation. change. Albeit rather crudely, this may be considered to be an observation-based estimate of the sensitivity of precipitation phase change to a warming environment. Figure 11. (a) Distribution and (b) elevation dependency of snowfall to precipitation (S/P) ratio estimated using ground-based SNOTEL and GHCND data.
Similar to the approach proposed by [48], we may estimate the S/P ratio based on the satellite radar data using CloudSat snowfall averaged from 2006 to 2017 and DPR rainfall averaged from 2014 to 2020. The S/P ratio maps estimated by satellite data and its dependency on elevation are shown in Figure 12. Compared to the results shown in Figure 11, although the satellite radar data derived elevation dependency of the S/P ratio is noisier in the scatterplot, the relation seems to follow a similar trend, i.e., about 20~25% increase in S/P ratio per every kilometer increase in elevation. The noisy relation may be improved by averaging data over a large area, for example, 1° × 2° as done by [48], which is beyond the scope of this study and will be examined in the future. The S/P ratio estimated from the satellite is somewhat lower (by 20% on average) than that estimated by ground data. Recall that CloudSat underestimates snowfall in this region by about 65% against ground measured grid mean values. The lower S/P ratio by satellite data is partially caused by the CloudSat snowfall underestimation. However, based on our evaluation, the GPM DPR retrievals also underestimate rainfall in this region by about 50%, possibly because the GPM rainfall, too, is derived from DPR reflectivity at levels of 1 km or above over actual surface. As a result, the S/P ratio underestimation becomes less severe than either the CloudSat snowfall or the GPM DPR rainfall underestimation.   Similar to the approach proposed by [48], we may estimate the S/P ratio based on the satellite radar data using CloudSat snowfall averaged from 2006 to 2017 and DPR rainfall averaged from 2014 to 2020. The S/P ratio maps estimated by satellite data and its dependency on elevation are shown in Figure 12. Compared to the results shown in Figure 11, although the satellite radar data derived elevation dependency of the S/P ratio is noisier in the scatterplot, the relation seems to follow a similar trend, i.e., about 20~25% increase in S/P ratio per every kilometer increase in elevation. The noisy relation may be improved by averaging data over a large area, for example, 1° × 2° as done by [48], which is beyond the scope of this study and will be examined in the future. The S/P ratio estimated from the satellite is somewhat lower (by 20% on average) than that estimated by ground data. Recall that CloudSat underestimates snowfall in this region by about 65% against ground measured grid mean values. The lower S/P ratio by satellite data is partially caused by the CloudSat snowfall underestimation. However, based on our evaluation, the GPM DPR retrievals also underestimate rainfall in this region by about 50%, possibly because the GPM rainfall, too, is derived from DPR reflectivity at levels of 1 km or above over actual surface. As a result, the S/P ratio underestimation becomes less severe than either the CloudSat snowfall or the GPM DPR rainfall underestimation.

Discussions
In this study, we proposed a new approach to validate satellite measured areal mean quantity by point station measurements. The novelty of this approach is that it utilizes the characteristic of climatological mean snowfall rate increasing with elevation, so that we can downscale the snowfall distribution with measurements at very sparsely distributed stations. Although the rate of increase of snowfall with elevation is different from location to location as shown in Figure 6, the general trend seems to be universal; it is true for the 768 SNOTEL locations we examined in this study and for the areal mean results in Tibetan Plateau as reported by [48]. It is noted that this trend is generally not true for rainfall, as shown in Figure 10. Therefore, the approach is not applicable to rainfall validations. It is also noted that the approach proposed here is conceptually different from spatial interpolation in filling measurement gaps, which uses the horizontal spatial pattern of measured quantities and needs the stations with known observables to be distributed relatively uniformly in the horizontal plane for a better interpolation. Our approach, however, uses the unique property of mean snowfall rate varying with elevation regardless of the horizontal distribution of existing stations. To the authors' best knowledge, we are the first ones proposing this approach.
One other issue may be raised related to the time mismatch for snowfall and rainfall estimations in the computation of satellite derived S/P ratio in Section 3.2. CloudSat snowfall retrievals are available only from 2006 to 2017 and GPM DPR rainfall retrievals are available from March 2014 onward. To compute S/P ratio, we assumed that the mean of each retrieval in a 1 • × 1 • box represents the climatological mean of snowfall or rainfall, implying that the roughly 10-year mean of snowfall or the 6-year mean of rainfall is much greater than the difference of them between the two time periods. To test this assumption, we analyzed the relative difference of mean snowfall rates between the two time periods. The results are shown in Figure 13. The differences are generally less than 20% with both positive and negative values, resulting in the difference averaged for all stations almost canceling out (3% for SNOTEL and 0.4% for GHCND stations). Considering that we use 1 • × 1 • averaged satellite estimates to compute S/P ratio, the uncertainty of mean snowfall rate caused by the time difference should most likely be only a few percent, a value significantly smaller than the uncertainty of the retrievals. Song and Liu performed a similar analysis using ERA5 reanalysis data [48] for the Tibetan Plateau and reached a similar conclusion, i.e., the variation of multi-year mean of precipitation being far less than either the mean values themselves or the observational uncertainties. Therefore, the S/P ratio derived from the satellite data is still physically meaningful.
posing this approach.
One other issue may be raised related to the time mismatch for snowfall and rainfall estimations in the computation of satellite derived S/P ratio in Section 3.2. CloudSat snowfall retrievals are available only from 2006 to 2017 and GPM DPR rainfall retrievals are available from March 2014 onward. To compute S/P ratio, we assumed that the mean of each retrieval in a 1° × 1° box represents the climatological mean of snowfall or rainfall, implying that the roughly 10-year mean of snowfall or the 6-year mean of rainfall is much greater than the difference of them between the two time periods. To test this assumption, we analyzed the relative difference of mean snowfall rates between the two time periods. The results are shown in Figure 13. The differences are generally less than 20% with both positive and negative values, resulting in the difference averaged for all stations almost canceling out (3% for SNOTEL and 0.4% for GHCND stations). Considering that we use 1° × 1° averaged satellite estimates to compute S/P ratio, the uncertainty of mean snowfall rate caused by the time difference should most likely be only a few percent, a value significantly smaller than the uncertainty of the retrievals. Song and Liu performed a similar analysis using ERA5 reanalysis data [48] for the Tibetan Plateau and reached a similar conclusion, i.e., the variation of multi-year mean of precipitation being far less than either the mean values themselves or the observational uncertainties. Therefore, the S/P ratio derived from the satellite data is still physically meaningful. There are a few studies in the literature on comparing CloudSat snowfall product with surface station measurements [52][53][54], ground-based radar measurements [55][56][57][58], and model reanalysis [59][60][61]. Most of these studies are conducted over high latitudes. There are a few studies in the literature on comparing CloudSat snowfall product with surface station measurements [52][53][54], ground-based radar measurements [55][56][57][58], and model reanalysis [59][60][61]. Most of these studies are conducted over high latitudes. Among them, King [53] evaluated CloudSat snowfall product against four station measurements and three gridded snow water equivalent products throughout the Canadian Arctic and found that CloudSat has better performance north of 70 • N, with underestimation compared to measurements at most of the stations and the reanalysis. Ryan et al. [61] derived a 15-km resolution snowfall climatology from CloudSat snowfall retrievals over Greenland ice sheet and concluded that CloudSat accumulation climatology has an uncertainty of ± 28% with respect to accumulation rates derived from ice cores. Edel et al. [60] compared CloudSat snowfall climatology with several reanalysis datasets and found that similar general geographical patterns are observed in all datasets, although there are significant mean snowfall rate differences over the Arctic between 58 • and 82 • N. Using conventional surface weather station data over Canada, Hiley et al. [52] found that CloudSat snowfall retrieval does not correlate well with surface station measurements, except for at some high latitude stations where CloudSat has more frequent sampling and mixed phase precipitation is less of an issue. While these studies indicate that CloudSat snowfall product is of great value in climate research, its accuracy is still unclear, and thus a quantitative assessment is still needed even in the high latitudes.
Validation of CloudSat snowfall over the midlatitudes, which is more relevant to this study, has been conducted using ground-based weather radar observations over the U.S. Of particular interest is the work done by Cao et al. [55], in which it is found that CloudSat underestimates heavier snowfall and the underestimation is well correlated to the snowfall intensity, a result that is consistent with the conclusion of this study. In a subsequent study, Chen et al. [56] found that there is a large discrepancy between CloudSat and surface radar snowfall estimates with a low correlation coefficient of only about 0.41, and the discrepancy seems to be related to snowfall intensity and the bin height where CloudSat surface snowfall is determined. On the other hand, Matrosov [57] compared the same surface weather measurements with CloudSat retrievals but resampled the surface radar data to match the bin height at which CloudSat snowfall is derived and found that the above correlation coefficient is much higher, being around 0.8 to 0.85 depending on the method used for the reflectivity to snowfall rate conversion. Therefore, it is argued that the discrepancy between CloudSat and surface radar retrievals is largely due to the CloudSat CPR's inability to measure near surface snowfall because of the contamination of ground clutter. We argued the same point in this study and pointed out that this problem is particularly severe in the mountainous regions with the combination of orographic precipitation enhancement and the existence of satellite radar blind zone. Clearly, future verification of this speculation is needed with joint analysis of surface radar, weather station, SNOTEL, and CloudSat data.

Conclusions
This study aims at assessing the accuracy of CloudSat snowfall retrievals in the mountainous regions in the western U.S. using ground-based snowfall measurements. When validating the satellite radar snowfall retrievals, we found that the two sets of ground truth data have systematic differences with each other: the SNOTEL has significantly higher values of annual mean snowfall rate than the GHCND does. Further investigation indicates that this systematic difference is a result of SNOTEL stations being mostly placed at higher elevations than the GHCND stations, and snowfall is generally heavier at high than low elevations in the mountainous regions. Therefore, we conclude that neither SNOTEL nor GHCND dataset alone can correctly represent an areal mean of snowfall rate, thus it cannot be used to compare with the areal mean values of CloudSat retrieved snowfall.
To solve this problem, a novel approach is proposed in this study, in which all SNOTEL and GHCND station data within a 1 • × 1 • area are used to develop a snowfall rate versus elevation relation. This relation is then used to compute snowfall rate in other locations within the 1 • × 1 • area, mimicking that surface observations are available everywhere in the grid. Grid mean snowfall rates surrounding all SNOTEL stations are then derived, which should be more representative of the mean snowfall rate than those derived by either SNOTEL or GHCND dataset alone. The so-derived grid mean snowfall rates are compared with CloudSat retrieved mean snowfall rates in corresponding grids. The results show that the CloudSat product underestimates snowfall by about 65% when averaged over all the 768 grids in the west U.S. mountainous region. The bias occurs regardless of elevation but strongly depends on mean snowfall rates in the grid. That is, the heavier the mean snowfall is in a grid, the more severe the underestimation will be.
As an application of the so-derived grid mean snowfall, we further estimated the snowfall to precipitation ratio from both ground and satellite measured data. The satellite estimation is based on an approach proposed by [48] using CloudSat to estimate snowfall and GPM DPR to estimate rainfall. The general distributions of the surface and the satellite-based estimation of S/P ratio have similar horizontal pattern, although the satellite estimated S/P ratio is somewhat lower than the ground-based estimation. The increasing rates of S/P ratio with elevation derived from ground and satellite-based data are quite similar, being about 20-25% per kilometer up, which translates to approximately 4% per degree Celsius of temperature drop.
This study introduced a new approach to validate satellite precipitation retrievals by surface point measurements where the surface precipitation intensity around the stations is not randomly distributed but depends on some topographic features, such as elevation as shown in this study. Although we only demonstrated the usage of the method for snowfall validation in the western U.S., the same strategy can be applied elsewhere, should the precipitation pattern there show a clear dependency on a topographic feature. Data Availability Statement: SNOTEL data are available from the USDA: https://www.wcc.nrcs. usda.gov/snow/snotel-data.html (accessed on 10 January 2022). GHCND data are available from NOAA/NCDC: ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/daily, (accessed on 10 January 2022). CloudSat data are available from CloudSat Data Processing Center at Colorado State University: https: //www.cloudsat.cira.colostate.edu, (accessed on 5 August 2019). GPM data are available from NASA data center: https://gpm.nasa.gov/data (accessed on 18 June 2020). SRTM data are available from https://topex.ucsd.edu/WWW_html/srtm15_plus.html (accessed on 6 March 2021). ERA5 data are available from ECMWF: https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5, (accessed on 10 October 2021).