Estimation of subpixel snow sublimation from multispectral satellite observations

Abstract. Snow sublimation is an important hydrological process and its spatial and temporal variation remains largely unknown; however, few studies have been conducted to quantify its spatial variability. Our study focuses on the evaluation of two algorithms, Penman–Monteith (P–M) equation and the bulk aerodynamic (BA) parameterization of snow sublimation. The two methods were first evaluated against eddy covariance (EC) measurements of latent heat flux at towers located in the upper reaches of the Heihe River Basin (China). Both methods were in good agreement with the ground observations with high coefficient of determination (R2) and low root mean squared error (RMSE). Next, we estimated subpixel snow sublimation using remote sensing data at a 1-km×1-km spatial resolution. The results based on satellite data were evaluated against ground measurements at the two experimental sites. The P–M equation gave R2=0.75, RMSE=8.4  W m−2 for Dashalong site and R2=0.36, RMSE=9.1  W m−2 for the Dadongshu site and performed better than the BA parameterization, which gave R2=0.65, RMSE=17.5  W m−2 for the Dashalong site and R2=0.06, RMSE=21.2  W m−2 for the Dadongshu site. Overall, the results indicate that P–M is promising for estimating snow sublimation at the regional scale using satellite observations.


Introduction
Snow plays an important role in the climate system of high elevations and cold regions because of its high albedo and low heat conductivity and it is also an important water resource for these regions. [1][2][3][4] Snow sublimation has been identified as an important hydrological process and results in a significant loss of snowpack water, which involves complex mass and energy exchanges. [4][5][6][7] Extreme cases of sublimation have been shown to be very efficient at removing snowpack water with losses of up to 90% of annual snowpack on preferred alpine crests. 5 A wide range of snow sublimation rates have been reported in the literature. [7][8][9][10][11] High sublimation rates were observed in the Colorado Front Range with daily values of 2.35 mm and in the alpine region of the Sierra Nevada with daily values of 2.17 mm. 8,9 However, sublimation rates were reported with 0.39 and 0.15 mm d −1 at wind-exposed and wind-sheltered sites, respectively, at Owyhee Mountains of the USA. 7 These various results may be explained by environmental differences among study sites because snow sublimation is affected by multiple factors such *Address all correspondence to: Li Jia, E-mail: jiali@radi.ac.cn as wind speed, net radiation, air temperature, and vapor pressure deficit, 5 and the spatial variability of snow sublimation is dependent on the spatial variability of these factors. 5,10,11 Sublimation as the moisture fluxes between the snow and the atmosphere is commonly measured by an evaporation pan. The evaporation pan has been widely used with the advantage of simplicity and economy but with more manual work. 4,[12][13][14][15] Schmidt et al. 13 estimated that sublimation from snowpack was 78-mm water equivalent and represented about 20% of the normal peak water equivalent of the snowpack by a 65-cm diameter pan during a 40-day accumulation period in a Colorado subalpine forest. However, measurements by an evaporation pan were found to be invalid when precipitation or blowing snow occurred. 15 The eddy covariance (EC) method is another approach to measure snow sublimation and is widely used. 7,16-21 EC uses high-frequency (10 to 20 Hz) measurements of wind vector and scalar quantities as the basis for flux computations. A detailed description of the EC method and data processing to determine sublimation can be found in Reba et al. 17 EC sensors respond poorly during precipitation events, 7,17 however. Additionally, the so-called energy balance closure problem is a widespread occurrence in EC measurements with reported discrepancy of 10% to 30%. 22,23 Despite these challenges, EC data can be a very valuable data source for snow research when thoroughly corrected to ensure accuracy. 7,17,18,21 In addition to in situ observation, methods to estimate snow sublimation were developed in the past several decades. Initially, empirical formulas were widely adopted, based on the statistical relationship between snow sublimation and meteorological conditions (e.g., wind speed and vapor pressure deficit). 24, 25 Strasser et al. 5 used empirical formulas to estimate sublimation from the snowpack to be about 100-mm snow water equivalent (SWE) during the 2003/2004 winter period in a mountainous region in the southeast of Germany. Empirical methods are simple but have low accuracy and poor physical basis. 8,10 The aerodynamic profile method was developed later based on complex turbulent transfer theory and is more reliable than the empirical methods. 8,10,14 Hood et al. 8 used the aerodynamic profile method to estimate snow sublimation: total net sublimation was 195-mm water equivalent or 15% of annual snow accumulation at Niwot Ridge in the Colorado Front Range during the 1994/1995 snow season. The snowpack is generally assumed to be a saturated surface, which has been proved by many studies, [26][27][28][29][30] thus the bulk aerodynamic (BA) method was developed to simplify the aerodynamic profile method and remove the need for multilevel observations. The BA model needs only meteorological measurements at a single height above the snowpack and soon became the most widely used method due to a widespread lack of multilevel meteorological data. 7,[26][27][28][29][30] Reba et al. 7 found that the measured snow sublimation by EC was very similar to estimated snow sublimation by the BA model with a correlation coefficient of 0.85 and root mean square difference of 0.15 mm d −1 .
Many of the physical snow sublimation models have been based upon the energy balance approach. 21,29,31,32 The relation between the amount of energy used for snow sublimation with net radiation and subsurface heat flux can be described by a conceptual model of the energy balance of a control volume, which applies to a single, thin layer of snow placed at the air-snow interface (Fig. 1).
The snow surface energy balance equation can be written as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 7 2 3 where the net radiation flux (R n ) is the principal driver of the snow sublimation (W m −2 ), the sensible heat flux (H) is the energy transfer due to convection in the surface boundary layer (W m −2 ), the snow subsurface heat flux (G s ) is the energy transfer to the snowpack (W m −2 ), and the latent heat flux (LE) is the energy used in the phase transition of frozen/liquid water to vapor (W m −2 ). The method based on the snow surface energy balance also applies one-level meteorological measurements to estimate snow sublimation with the P-M equation, which has been widely adopted in estimating evaporation. 33,34 The P-M method for estimating latent heat flux calculation has been extensively used in hydrological, atmospheric, and environmental modeling, and can be used to estimate snow sublimation as air in the vicinity of the snowpack that is generally assumed to be vapor saturated. 29,31,32 Furthermore, many studies on snow sublimation are constrained by the unavailable or very sparse ground observation. At high latitudes, meteorological stations are sparse and located mostly at lower elevations. 2 However, snow is more abundant at high latitudes and altitudes. Although many ground-based observations and modeling studies of snow sublimation have provided valuable information (i.e., snow surface temperature and snow depth) at a particular location, it is difficult to extend such information to large scales. The spatial variability of snow sublimation over the vast regions remains largely unknown. 2 Satellite observations have long been recognized as effective in providing spatially distributed snow information to estimate snow sublimation. 35 Remotely sensed data provide temporally and spatially continuous information over snow surface, i.e., snow surface temperature and fractional snow cover (FSC). 36,37 The moderate resolution imaging spectroradiometers (MODIS) onboard NASA's Terra and Aqua satellites provide unprecedented information regarding snow and surface energies, which can be used for regional-and global-scale snow sublimation estimation in near realtime. 38,39 However, to our knowledge, very few studies have been published on the estimation of snow sublimation with satellite observations.
The BA method and P-M equation were employed in this study because both methods assume that the snow surface is vapor saturated and use atmospheric measurements at only one level above the snow surface. This is the foundation of calculating the snow sublimation by remote sensing because we cannot get gradient measurements in the surface boundary layer with satellite observations. We present an approach to estimate subpixel snow sublimation using multiple satellite observations by characterizing the subpixel energy balance of the snow cover fraction. We have evaluated our method using both tower and satellite observations as input and compared alternate parameterizations of snow sublimation. The objectives of our study were: (1) to evaluate the P-M equation and the BA parameterization using tower observations, (2) to analyze the energy balance of the subpixel snow cover fractional area, and (3) to present and validate the regional snow sublimation estimated using satellite data.

Description of Study Area and Measurements
The upper reaches of the Heihe River Basin (Fig. 2) are located between 37°41′ N to 39°05′ N and 98°34′ E to 101°11′ E, and cover an area of ∼2.85 × 10 4 km 2 . Elevation ranges from 1631 to 5245 m, the mean annual precipitation is 350 mm and the mean annual air temperature is 2.0°C. 40 Glaciers covering 421 km 2 are distributed above elevation 4500 m. At elevations above 4000 m, vegetation is very sparse and dominated by cushion plants, while meadows and shrubs occur below 3300 m. 41 Temporary snow usually exists under 2700 m, patchy snowpack from 2700 to 3400 m, and permanent snowpack exists above 3400 m. 42 Snowmelt provides most water resources in the Heihe River Basin, which supplies agriculture in the middle reaches and the arid ecosystem in the lower reaches. 43 The data were collected at two observation stations, Dadongshu mountain pass station (4101 m, 100.23°E, 38.01°N) and Dashalong station (3739 m, 98.94°E, 38.84°N), part of the Heihe Watershed Allied Telemetry Experimental Research (HiWATER) project. 44,45 These stations were selected because of the locally flat terrain and homogeneous land cover without major obstacles, such as houses and trees. The land cover at these sites is alpine meadow and bare soil. The land surface at these sites is always covered by snow most of the time in the late autumn, winter, and early spring, but snowfall is larger in spring and autumn than in winter. 46 The sensible and latent heat flux data during the period between November 1, 2014, and January 31, 2015, at Dadongshu station, and October 17, 2014, and December 31, 2014, at Dashalong station were collected. The fluxes were measured by the EC systems installed at 3 m at the Dadongshu site and 4.5 m at Dashalong site above the ground. The raw data were sampled at 10 Hz, then 30-min mean values were obtained after a series of corrections, such as data filtering, sonic temperature correction, air density correction, and coordinate rotation. 44 The auxiliary data, including the net radiation, air temperature, relative humidity, snow depth, air pressure, wind speed, and snow infrared surface temperature, were also observed by automatic weather stations (AWS) near the EC system in the same period. All the data were processed and recorded at 30-min intervals. There were no snow depth measurements at Dashalong station. The latent heat flux at this station was used to validate the estimates based on satellite data.

Satellite and Regional Forcing Data
The satellite and regional forcing data used in this study are listed in Table 1. They are acquired from various sources and comprised of well-validated products.
Land surface temperature was retrieved using the radiometric data acquired by the MODIS on Terra and downloaded from National Aeronautics and Space Administration (NASA), which covers the entire upper reach of the Heihe River Basin.
We used the FSC retrieved with the MODIS radiometric data (MOD10A1 product) and downloaded from the National Snow and Ice Data Center (NSIDC). The MOD10A1 FSC is estimated by applying an empirical relationship with the normalized difference snow index (NDSI), established using FSC retrieved at higher spatial resolution with Landsat ETM+ data. 37 The original MOD10A1 has a 500-m spatial resolution and was resampled to a 1-km grid.
Land cover type is required to unmix the land surface temperature and estimate the radiometric surface temperature of the snow cover fraction. The land cover map at a spatial resolution of 1 km was obtained from the Cold and Arid Regions Sciences Data Center. 47,48 The land cover map was generated by combining multisource land cover/land use classification maps including a 1:1,000,000 scale vegetation map, a 1:100,000 scale land use map for the year 2000, a 1:1,000,000 scale swamp-wetland map, a glacier map, and a MODIS land cover map for China in 2001. 48 The dominant land cover types in winter in these regions were forest and bare soil. We aggregated the land cover classes to just two: (a) snow and forest and (b) snow and bare soil.
The regional atmospheric forcing data were obtained from the Cold and Arid Regions Sciences Data Center. [49][50][51] Data with a spatial resolution of 5 km were downscaled to 1 km by applying a statistical downscaling approach to hourly downward longwave and shortwave radiations, air temperature, air pressure, and specific humidity. 52 Data of wind speed and precipitation were downscaled to 1 km using the bilinear resampling method.

Methodology
Two methods to estimate snow sublimation were evaluated: the widely used BA method and P-M combination equation adapted to snow cover. The two methods were applied assuming the air near the surface to be vapor saturated and the results were compared with the snow sublimation measured by the EC system in two experimental stations in the upper reaches of the Heihe River Basin in China. The theoretical background of these two methods is described in the following sections.

Penman-Monteith Combination Equation
The Penman 53 equation was first derived to estimate evaporation from open water and water-saturated surfaces, where equilibrium at the surface, i.e., under equilibrium with vapor-saturated air. Later, this equation was extended by Monteith 33 to nonsaturated surface, i.e., to actual evaporation. The combination equation is obtained by rewriting the surface energy balance equation using the Clausius-Clapeyron equation for (liquid and frozen) water-vapor equilibrium. In Penman and Monteith a liquid water-vapor equilibrium was assumed, while in our case the equilibrium is frozen water (ice)-vapor, i.e., we need to use the appropriate dependence of the saturated water vapor pressure on temperature. The P-M combination equation was applied in the form: 33 For full snow cover E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 1 3 9 (2) For FSC E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 1 0 0 where LE (W m −2 ) is the latent heat flux with L as latent heat of sublimation (J kg −1 ) and E (kg s −1 m −2 ) as sublimation, FSC (-) is the FSC, Δ (Pa K −1 ) is the slope of the curve relating saturated vapor pressure in the ice-water vapor equilibrium e sat (Pa) to air temperature (T a ), and for ice or snow, Δ ¼ 21.87×265.5×e sat ð265.5þT a Þ 2 , (here T a is in°C). 31 R n (W m −2 ) is the net radiation flux, G s (W m −2 ) is the snow subsurface heat flux due to the energy transfer in the snowpack, and G s was estimated using a linear relationship between G s and R n obtained from AWS and EC observations. ρ (kg m −3 ) is the air density, C p (J kg −1 K −1 ) is the specific heat capacity of air, e sat (Pa) is the saturation vapor pressure at the air temperature (T a ) in the frozen water-air equilibrium, i.e., e sat ¼ 611 × expð 21.87×T a T a þ265.5 Þ (here T a is in°C). e (Pa) is the actual water vapor pressure, r a (s m −1 ) is the aerodynamic resistance, and γ (Pa K −1 ) is the psychrometric constant. Equation (2) is a particular form of the combination equation. Menenti derived a general equation and showed that both the Monteith and Penman equations can be obtained as limiting conditions of the Menenti equation. 54 Both the Monteith and the Penman equations applied to a case, where the liquid (frozen) water-vapor phase transition occurs at a surface, where air is vapor saturated at its temperature. More in general, the liquid (frozen) water-vapor phase transition occurs at locations underneath the surface of the evaporating (sublimating) body, such as soil, stomata, or the snowpack particularly when meltwater may be present within the snowpack. Our assumptions, i.e., that the sublimation occurs at the surface of the snowpack and that air in the vicinity of the surface is vapor saturated at the temperature of the snowpack surface, may lead to differences between sublimation estimated with Eq. (2) and total latent heat flux measured by an eddy covariance device. The latter may include contribution due to subsurface sublimation and evaporation of meltwater.
The aerodynamic resistance is calculated as follows: where Ri (−) is the Richardson number given as 28 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 2 8 3 where g is the acceleration due to gravity (9.8 m s −2 ), T a (K) is the air temperature, T snow (K) is snow surface temperature, and T mean ¼ ðT a þ T snow Þ∕2 in K. The MOD11A1 land surface radiometric temperature retrieved from MODIS radiometric data is an average applying to an entire pixel. We assumed that in the case of mixed pixels, e.g., composed of snow and other land cover types such as bare soil or vegetation, the fractional cover of snow FSC is the MODIS data product, while the fractional abundance (1 − FSC) is assigned to the land cover type given by a land cover map for the entire pixel. To retrieve the surface temperature of the snow fractional area, we applied the method proposed by Sun et al. 55 The component temperatures of a mixed pixel are assumed to include just two endmembers, in our case either snow and vegetation (forest) or snow and bare soil. The end-member radiometric temperatures were estimated using MODIS land surface temperature and the spectral emissivities in the MODIS bands 31 and 32 , and G ¼ f m ε m31 ; T m and T n are the end-member temperatures within a mixed pixel, T s (K) is the MOD11A1 pixel land surface temperature, f m and f n are the fractional covers of end-members, ε m31 and ε n31 are the end-member emissivities in band 31 (0.982319 for snow, 0.9851 for forest, and 0.9832 for soil), ε m32 and ε n32 are the end-member emissivities in band 32 (0.962614 for snow, 0.9844 for forest, and 0.9731 for soil), ε 31 and ε 32 are the emissivities of the pixel in bands 31 and 32 and were estimated as fractional area weighted average of the end-member emissivities. If FSC ¼ 1, then the temperature of the pixel is the snow temperature; if 0 < FSC < 1, the fractional abundance (1 − FSC) is assigned to the land cover type given by a land cover map for the entire pixel. Once the fractional abundances f i have been determined, Eqs. (8) and (9) are applied to calculate subpixel snow temperature. Net radiation flux (R n ) is the principal driver of the snow sublimation and was expressed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 6 ; 5 3 6 where R swd is downward shortwave radiation (W m −2 ), ε is the spectrally integrated emissivity of the snowpack (-), R lwd is downward longwave radiation flux (W m −2 ), and α is the albedo of the snowpack (-) calculated by an empirical formula, 56 For melt condition∶ α t ¼ ðα t−1 − α min Þ · e − 0.24×Δt 86;400 þ α min ; where Δt is the time step after snowfall (s), τ a is the empirical coefficient (0.008), and α min is the minimum snow albedo (0.5). Extensive analyses of snow albedo and snowfall events indicate that 3-mm SWE of snowfall is a threshold for newly fallen snow to cover the old snow and resets the snow surface albedo to a high value. 58 In this work, a snowfall exceeding 3-mm SWE is defined as a "new snow" event, and the time between two consecutive "new snow" events is counted as length of elapsed time after the first snowfall. If the precipitation (snowfall here) from the forcing data is >3 mm (SWE), then the snow albedo will be the maximum value of 0.85.

Bulk Aerodynamic Model
According to the BA theory, the latent heat flux can be computed as: 28,30 For full snow cover where P a is the atmospheric pressure (Pa), e s (Pa) is the saturation vapor pressure at the snow surface temperature (T snow ), e s ¼ 611 × expð 21.87×T snow T snow þ265.5 Þ, (here T snow is in°C). C e (m s −1 ) is the bulk transfer coefficient for vapor exchange, which can be expressed as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 1 1 6 ; 1 2 0

Metrics for Accuracy Assessment
We applied multiple metrics to evaluate the accuracy of our estimated sublimation. The mean relative error (MRE) describes the ratio between the deviation of the estimate from the observation and the observed value. The root mean square error (RMSE) measures the absolute difference between estimation and observation. The coefficient of determination (R 2 ) is a measure of the consistency of estimated with measured values, e.g., in time. These metrics were calculated as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 1 1 6 ; 6 4 4 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 1 1 6 ; 5 9 6 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 8 ; 1 1 6 ; 5 5 6 where P i and O i are the estimates and observations, respectively, with P and O being the averages and n is the number of estimate-observation pairs.

Meteorological Conditions
At Dadongshu station, there were snowfall events from October 10, 2014, as shown in Fig. 3   the Dadongshu site was 45.2%, ranging from 9.1% to 87.2%. High relative humidity, which means large amount of water vapor in the air, will inhibit snow sublimation.

Snow Subsurface Heat Flux
The snow subsurface heat flux (G s ) was not measured and we estimated as the residual of the energy balance at the surface of the snowpack. This was done only in the Dadongshu site, when it was covered by a homogeneous snowpack during the snow season, for all available complete daytime (11:00 to 17:00) measurements between November 1, 2014, and January 31, 2015. Next, we determined a linear correlation between G s and R n , obtaining a regression coefficient of 0.575 and R 2 ¼ 0.63 (Fig. 5) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 9 ; 1 1 6 ; 2 8 9 This relationship in Eq. (19) is further applied to estimate regional snow sublimation using satellite data.

Validation of Estimated Snow Sublimation
We evaluated first estimates of sublimation obtained with both Eqs. (2) and (13) using in situ measurements at the Dadongshu site of all the energy balance components, i.e., net radiation, latent heat flux, and sensible heat flux for the period with complete snow cover (FSC ¼ 1). The  (Fig. 6), while RMSE was 10.48 and 10.21 W m −2 , respectively. The MRE for P-M was 28.1%, while it was 1.12% for BA. The relatively small MRE of BA is due to the offset effect of the positive and negative biases. It is clear that the P-M equation overestimates the sublimation of snow more than the BA approach (Fig. 7). Overall, the evaluation against in situ measurements showed that the performance of both P-M and BA is satisfactory, especially taking into account the low RMSE values.

Sensitivity Analysis
The sensitivity of estimated sublimation to input data and parameters was evaluated to identify the most influential ones and understand how to deal with data of different accuracies. The sensitivity of estimated snow sublimation to net shortwave radiation, snow surface temperature, air temperature, wind speed, and relative humidity was evaluated. The sensitivity analysis (SA) was performed on one variable at a time by assuming a range between AE20% and AE10% of error for each input variable, while keeping the others constant. 59 A nondimensional sensitivity coefficient was calculated as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 0 ; 1 1 6 ; 3 3 9 where SC is the dimensionless sensitivity coefficient, LE and CV are the reference values of latent heat flux and input variable, respectively, ΔLE is the change of estimated sublimation when ΔCV is changed for AE20%. The reference value of each variable was the mean value between 11:00 and 17:00 during the period November 1, 2014, till January 31, 2015, at the Dadongshu site. The reference values of both input variables and estimated sublimation are given in Table 2. The sensitivity results are shown in Fig. 9. The SA indicates that the wind speed is the most influential variable when using either P-M or BA. The P-M equation is less sensitive, i.e., SC ¼ 0.77, to wind speed, however, than BA, i.e., SC ¼ 1.2. This is a clear advantage of the P-M equation for large area estimates because the regional wind field is less accurate than in situ measurements. 52 Another influential variable is relative humidity with SC ¼ −0.81 and −0.52 for BA and P-M respectively. Relative humidity has a negative impact on LE for both methods. The influence of net shortwave radiation on P-M estimates of sublimation is comparable to relative humidity, i.e., SC ¼ 0.52. Both approaches showed the least sensitivity to air temperature and surface temperature. These agreed with former studies on the sensitivity of snow sublimation calculation to the error of input variables, which highlighted that snow sublimation was highly sensitive to error of winds speed and relative humidity, while not sensitive to the error of temperature. 31,60,61

Regional Estimates of Snow Sublimation with Remote Sensing Data
The evaluation of the P-M and BA methods with in situ measurements shows that the P-M estimation is more accurate. Since it is derived from the surface energy balance equation, the P-M equation is also most suitable for regional analyses based on the multispectral radiometric data collected by spaceborne imaging spectroradiometers. This notwithstanding, we applied both methods for our regional case study. Regional snow sublimation was estimated for the entire upper reaches of the Heihe River Basin by applying retrievals from satellite data, i.e., subpixel snow surface temperature and FSC, in combination with atmospheric forcing data. The regional analyses were carried out under clear sky conditions within the period from November 1, 2014, to January 31, 2015. The snow component surface temperature was obtained using the method described in Sec. 3.1.
The map of estimated sublimation obtained on a fully clear day (i.e., November 8, 2014) over the entire upper reaches of the Heihe River Basin was chosen to demonstrate our method. The FSC was rather high and the subpixel snow surface temperature was in the range of 260 to 273 K (Fig. 10). The spatial distribution of snow sublimation estimated by the P-M model and BA model is shown in Fig. 10. The estimated sublimation over the study area ranges between 0 and 80 W m −2 with the highest frequency values around 0 to 20 W m −2 when applying the P-M equation [Figs. 11(a) and 11(b)]. The BA estimates were similar, although the highest frequency values were slightly lower, i.e., around 0 to 10 W m −2 [Figs. 11(c) and 11(d)].
The regional estimates of snow sublimation under clear sky conditions were evaluated against the in situ EC measurements at the Dadongshu and Dashalong sites. In the period November 1, 2014, till January 31, 2015, there were 18 usable remotely sensed measurements at Dadongshu and 17 measurements at Dashalong. The P-M estimated sublimation at Dashalong was in good agreement with the measurements (Fig. 11)  It is important to note here that the EC measurements are essentially local-scale measurements representing only very small spatial scales when compared with the 1-km pixel size of the snow sublimation estimated with satellite data. These observations were considered; however, the best option available for evaluating the model performance because we lacked the consistent and long-term snow sublimation data measured by other methods.

Discussion
The daily mean value of snow sublimation measured by EC from November 1, 2014, till January 31, 2015, at Dadongshu was 0.15 mm d −1 , compared with 0.41 mm d −1 in the Rocky Mountains of Colorado by Molotch et al. 62 and 0.16 mm d −1 in Mongolia by Zhang et al. 10 The magnitude  of sublimation has been shown to vary widely across different land surface environments and elevations. 5,10,11,62 The different values of snow sublimation may be explained by environmental differences among study sites. 5 Our estimates of daily sublimation are well within the bounds of previous research.
Both the P-M and the BA methods performed well at the Dadongshu site. The lower sensitivity to the wind speed gives the P-M method an enormous advantage when applied to regional scales considering the scarce knowledge of the large area wind field at the 1 km × 1 km spatial resolution. Despite the ability of the P-M method to capture the spatial and temporal pattern of snow sublimation distribution, there are multiple uncertainties to deal with. First, we did not account for sublimation from blowing snow. Pomeroy and Essery 63 found that blowing snow sublimation has been estimated to account for as much as 10% to 50% of snowfall in North American Prairies and Arctic environments. The occurrence of blowing snow in the upper reaches of the Heihe River Basin has most likely significant impacts on the estimates of sublimation during the snow season. The actual amount of sublimation in this catchment could reasonably be higher than what has been reported here.
The second source of error is the accuracy of the spatial input data. These errors can result from both inaccurate forcing data and remote sensing data. Data on hourly 2-m air temperature, 2-m relative humidity, and snow surface temperature were in good agreement with tower measurements, while estimates of 10-m wind speed correlated poorly with in situ measurements (Fig. 14). The latter is a major challenge in regional analyses of sublimation with either method, given the poor accuracy of the large area wind field. 64 Errors on relative humidity, air temperature, and snow surface temperature are smaller but not negligible.

Conclusions
Snow sublimation plays a significant role in hydrological processes in cold regions and seasons, and it varies with atmospheric and surface conditions. In this study, we developed a method to  estimate snow surface sublimation with satellite observations and validated the results with EC measurements of latent heat flux. We applied the P-M equation and the BA parameterization, which use satellite data products to estimate regional snow sublimation at a 1-km spatial resolution. Model estimates have been compared with ground EC measurements. The P-M estimates were in good agreement with measurements at two sites (Dashalong and Dadongshu), with RMSE ¼ 8.4 and 9.1 W m −2 , respectively. The BA estimates were less accurate. The encouraging performance of the P-M suggested the extension of this approach to large area estimates of snow sublimation based on multiple satellite data. Remote sensing data at moderate spatial resolution enable us to estimate and map snow sublimation over large areas. The applicability and accuracy of the P-M method to estimate snow sublimation relies on the availability and accuracy of a suite of remotely sensed input data products, i.e., snow surface temperature, FSC, and atmospheric forcing data. Anyway, remote sensing has a large advantage to get the snow information with large-scale supervising ability. However, there are still some challenges. For example, clouds will greatly influence the quality of remote sensing data such as snow surface temperature and FSC. It is important to realize that the present study is limited by the unavailability of clear sky remote sensing data. As such, higher quality and accuracy of these input data are vital and expected.

Disclosures
The authors declare no conflict of interest.