1.1.

Remotely Sensed Snow Distribution

Microwave snow depth (SD) and snow water equivalent (SWE) products have not been sufficiently developed for monitoring snow cover in rugged mountainous regions.^9,10 Especially in regions of shallow snow in western China, the passive microwave algorithm with gradient brightness temperature between 36.5 and 18.7 GHz has been shown to be of little value.¹¹ Snow products from the Advanced Microwave Scanning Radiometer for EOS (AMSR-E) overestimated SWE in the TP according to the previous study.¹²

Visible and near-infrared (VNIR) RS data, especially snow cover fraction (SCF) data, are more frequently used in the snow distribution simulation,^13–15 owing to their high temporal and spatial resolutions and accuracy for identification of snow areas. However, the main disadvantage of the VNIR snow data is that they are influenced heavily by clouds.¹⁶

There are studies directly combining MODIS and AMSR-E data,^17–19 that take advantage of both the high spatial resolution of VNIR data and the cloud transparency of passive microwave data.²⁰ Some studies used both Terra and Aqua MODIS snow cover products to reproduce a new composited snow cover dataset, with a user-defined cloud cover threshold.²¹

1.2.

Simulation of the Snow Distribution Using RS Data

In some studies, RS data were used as a criteria for whether snow processes exist, such as in the snow runoff model (SRM).²² In other studies, RS data were used to validate, calibrate, and adjust the model output. They are usually used in physically based models in which the snow cover area could not be regarded as direct input into the model.²³ Some studies incorporated remotely sensed SCF data with the simulated daily snowmelt to improve the accuracy of the simulated distributed snowmelt.²⁴ More generally, data assimilation (DA)—in which different RS data are combined with the snow model—has proven to be the most efficient method for estimating the snow distribution. The SCA/SCF and SD data could be used alone or combined in different assimilation experiments.²⁵

1.3.

Snow Depletion Curve (SDC) and its Application to DA

The SDC, representing the functional relation between SCF and SD/SWE,²⁶ has been widely used to map SCF to SD in assimilation experiments.

Using the concept of the SDC, Clark et al.²⁷ assimilated SCF information into a simplified snowmelt scheme and obtained results with minor improvements in the accuracy of streamflow simulations. Besides SCF, remotely sensed SWE data have also been jointly assimilated into a snow model, yielding significantly improved root-mean-square error (RMSE) and correlation values.²⁸

At present, SCF is the most sufficient RS data used in the DA of snow in mountainous regions with few meteorological stations. The key for DA using SCF data is the application of the SDC. However, it is not easy to find a suitable SDC in ungauged regions because of the absence of enough prior in situ observations, especially in the TP.

1.4.

Challenges in Simulating the Snow Distribution in the TP

In the TP, it is difficult to simulate the snow distribution because of the many limitations, especially the scarcity of station observations. A viable method is to utilize RS data and an empirical snow model. For example, Immerzeel et al.²⁹ monitored the snow distribution using MODIS snow products and an empirical degree-day method in Himalayan river basins.

The difficulties that should be accounted for in simulating snow distribution in the TP can be summarized as follows:

Accounting for these problems, we present a synthesized scheme for the simulation of the snow distribution in the TP. This scheme includes a combination method of remotely sensed SCF and SD data. The reproduced SD distribution data are assimilated into a physically based snow model, using a relevant SDC derived from reliable station observations.

2.1.

Study Region and RS Data

The study region of the TP is illustrated in Fig. 2. There are 94 stations with SD observations used in the study. Because of the data missing and short snow seasons in some regions of the TP, some stations have few observations with SD greater than zero. All the observations from 94 stations are used to compare with the assimilated results.

Fig. 2

The TP region with the meteorological stations overlaying.

In the study, we used the calibrated AMSR-E SD data from the long time series of passive microwave satellite dataset in China, which is provided by the Environmental and Ecological Science Data Center for West China, National Natural Science Foundation of China ( http://westdc.westgis.ac.cn). In this dataset, the brightness temperatures from different sensors (SMMR for 1978 to 1987, SSM/I for 1987 to 2008, and AMSR-E for 2002 to 2010) were cross-calibrated.³² SD data were retrieved using the modified Chang algorithm,³³ which is dynamically adjusted based on the seasonal variation of grain size and snow density. It has been validated by considering the influences from vegetation, wet snow, precipitation, cold desert, and frozen ground. The extent of the snow distribution from this dataset was indirectly validated by MODIS snow cover products.¹²

The continuous SCF data with 500-m resolution were reproduced based on the cubic spline interpolation algorithm, and the accuracy of this dataset has been demonstrated.³⁴ The used dataset for reproduction are from MODIS/Terra Snow Cover Daily L3 Global 500-m grid products (MOD10A1), which is developed from a linear fit of Thematic Mapper snow cover to the normalized difference snow index of MODIS.³⁵

2.2.

SDC for Patchy Snow Distribution

Since the SCF values are $<100\%$ where snowpacks are shallow in most regions of the TP, the SDC, which represents the empirical corresponding relation between SCF and SD, is more useful for snow modeling. However, in situ observations of SD at a single station are not abundant for calibrating the regional SDC. Here, we present an alternative method by collecting a set of station observations.

The basic idea is that, since it is difficult to calibrate the SDC in a specific station, we choose all stations and their adjacent regions with a resolution of 0.1 deg as subjunctive regions. Then the in situ observations become sufficient and more reliable for calibrating a general SDC in the subjunctive region. In this subjunctive region, there are $n$ grids and $n$ SD observations of stations. Each grid has the same area and the SD was measured at each station in each grid. It is assumed that the SD observation at a station can represent the mean SD of the snow-covered area of the grid. Then for the total subjunctive regions, the mean SD can be evaluated from

The corresponding SCF of the subjunctive region can be obtained from

Using Eqs. (1) and (2), we can obtain a pair of $h_s$ and ${\mathrm{scf}}_s$ for 1 day. For a set of days, an SDC is expected to be regressed by pairs of $h_s$ and ${\mathrm{scf}}_s$ of the subjunctive region.

Of course, for a reliable calculation, $n$ should be large enough and the stations should be randomly distributed throughout each grid. Moreover, the SCF of the region is known because we have the daily MODIS SCF dataset. Then, the SDC of the subjunctive region can be evaluated.

Using the above idea, we collected in situ SD data of 94 stations in the TP and used the SCF data of each grid from the improved MODIS dataset. This gives a subjunctive region with the 94 grids in which each has an SD observation. The scale of the subjunctive region is about 0.1 deg × 0.1 deg.

An exponential function is used to fit the SDC between SCF and SD (Fig. 3). The coefficient of determination ($R^2$) of the regression is 0.67. The statistical SDC is used here for transferring SCF to SD. It is assumed that the calibrated SDC is applicable in constructing a relation between SCF and SD at the 0.1-deg scale. The SDC can be formulated as

Fig. 3

The gray curve is the fitted SDC using data from 94 stations in the TP. The years for the calibration are 2006, 2007, and 2008.

As shown in Fig. 3, there are no SCF values greater than 50%, which means that Eq. (3) cannot be used in cases where SCF is greater than 50%. When SCF is greater than 50%, we adopt the AMSD combining the microwave remotely sensed SD data and the VNIR SCF data. We investigate the reliability of microwave SD data when SCF is greater than 50%. The $R^2$ is about 0.27 between them. The accuracy is not high, however, the microwave SD may be the only available data for the assimilation scheme when SCF is greater than 50%, because there is no prior information for constructing an SDC in these cases. On the other hand, the cases when SCF is greater than 50% are few. Our investigation indicates that these cases only account for 8–18% in snow covered grids in different years in the TP.

2.3.

Downscaling of AMSR-E SD Data in Combination with MODIS SCF Data

In this section, we define a combination rule for reproducing a new SD dataset (with a spatial resolution of 0.1 deg), using the improved AMSR-E (and SSM/I) SD data (with a spatial resolution of 0.25 deg),¹² and the improved MODIS SCF data (with a spatial resolution of 500 m).^34,36 The new synthetic SD dataset combining AMSR-E (and SSM/I) and MODIS datasets is referred to as AMSD in the following text.

There are two assumptions needed for this rule. First, although SD is inhomogeneous in a grid of microwave data with a spatial resolution of 0.25 deg, we assume that the SD heterogeneity can be effectively described by using MODIS SCF information (500 m). Second, the SD value from AMSR-E (and SSM/I) is assumed to be the average value of the grid.

Then in a 0.25-deg grid (which is denoted as $A$),

The MODIS SCF data are resampled to 0.05 deg first. Then, the SD data are redistributed from a resolution of 0.25 deg to a resolution of 0.05 deg using Eq. (4). For avoiding obvious data errors, three special cases are accounted for:

An example is illustrated in Fig. 4.

Fig. 4

Illustration of the filtered SCF data at 0.1 deg (a), AMSR-E SWE of 0.25 deg (b) and their combination of 0.1 deg (c) on March 25, 2008.

2.4.

Distributed Meteorological Data

There are few meteorological stations in the TP. It is not enough to provide a reliable spatially distributed meteorological dataset for the vast region. For reliable meteorological forcing data in simulating the SD distribution, we adopt the downscaled meteorological data from the WRF model.³⁷ The WRF model is a fully compressible, Euler nonhydrostatic mesoscale forecast model with a runtime hydrostatic option. This model is effective in downscaling weather and climate data at scales from 1 to 100 km.

By using meteorological data with a spatial resolution of 1 deg and a temporal resolution of 6 h from the global forecast system (GFS), downscaled daily meteorological data with a spatial resolution of 0.1 deg were simulated for the TP. This dataset has been validated in the Heihe River Basin, which is located in the northeastern TP, using daily meteorological observation data from 15 stations of the China Meteorological Administration and hourly meteorological observation data from Watershed Airborne Telemetry Experimental Research project sites.^38,39 Comparisons indicate that the WRF model simulations of air temperature, surface pressure, and relative humidity are very reliable and that the precipitation, wind speed, and downward shortwave radiation are relatively good.⁴⁰

2.5.

SnowModel

We choose SnowModel³¹ as our model operator. SnowModel is a distributed snow evolution system which has been used and validated in many snow regions. This model system includes a distributed meteorological forcing conditions module, a surface energy exchange module (EnBal), a snow mass evolution module (SnowPack), and a snow redistribution module. Because the meteorological forcing has been generated from the WRF model, and the snow redistribution at the 0.1-deg scale is ignored in this study, we choose EnBal and SnowPack modules to simulate the snow distribution evolution. It is assumed that the snow evolution is homogeneous in the snow-covered area in a 0.1-deg grid. The vegetation distribution data are from the Environmental and Ecological Science Data Center for West China, National Natural Science Foundation of China ( http://westdc.westgis.ac.cn).

2.6.

DA Strategy

In the study, the EnKF method⁴¹ is used to assimilate the improved SD estimates into SnowModel. We use 50 ensemble members in the ensemble forecast. Selected forcing and state variables are perturbed; these include the four forcing datasets (downward shortwave radiation, downward long-wave radiation, precipitation, and air temperature) and the model state SD and the state variable (combined 0.1 deg AMSD data). The perturbation rule used here follows that of the previous studies.^28,42 The SD is initially perturbed. The four forcing data are perturbed during loops. Since we use a one-dimensional EnKF in the study, the perturbations of forcing data are not spatially correlated. The shortwave radiation and precipitation suffer from multiplicative perturbations, whereas perturbations of air temperature and longwave radiation are additive. The mean values of all these perturbations are zero. The standard deviations and cross-correlations are listed in Table 1.

Table 1

Perturbation parameters in the EnKF estimation of SD.

3.1.

Simulation Results of the SD Distribution

The time series of the assimilated SD dataset is calculated from 2000 to 2008. The monthly mean SD values are presented in Fig. 5.

Fig. 5

Comparison of different evaluated snow depths from 2000 to 2008 in the TP.

When the mean SD ranges from 0.005 to 0.010 m, the downscaled microwave SD is close to the evaluated SD using the SDC method. When the mean SD is $<0.005\mathrm{m}$ and $>0.010\mathrm{m}$, there is a large difference between the evaluations. Because the SDC method is ineffective when the SCF is equal to 100%, the disagreement in the high-value range is reasonable. However, when the snow is shallow, the microwave SD values have a greater accuracy. This is the main reason why the two results have a greater difference in the lower value range. The assimilated results are slightly less than the downscaled microwave SD values and the evaluated SDs using the SDC method. Some studies have indicated an overestimation of SD from AMSR-E data;⁴³ this can be attributed to the limitation of the current AMSR-E algorithm to account for the large grains that typically develop in snowpacks.¹⁰ Lower SD results may be more reasonable.

3.2.

Spatial Validation Using Distributed Snow Stations

One of the dominant constraints in simulating snow distribution using models is the lack of “real” validation data. Of course, remotely sensed SD and SCF data can be regarded as possible validation data sources. However, they are indirect. Unknown system and algorithm errors may be involved when using remotely sensed data. In contrast, in situ station observations are relatively reliable and are continuous over multiple years, but the scale effect would affect their representative capacity. Because complex topography dominates the TP, in different study grids with complex terrain the mean SD on the grid scale can be significantly different from the point SD.

Following a method similar to that in Sec. 2.2, we design a random sampling comparison method for a reliable validation. We construct a subjunctive region with $n$ stations, but we do not use a single point observation to validate the assimilated results. By using Eq. (3), the mean SD value of the subjunctive region can be evaluated using observations at stations. If the assimilated results are correct, these should be approximately equal to the assimilated mean SD of the selected subjunctive region. Of course, the station locations should be randomly distributed within each grid. If these stations are all randomly distributed in all grids, then Eq. (3) is more reliable, since we use a random sampling method.

First, those grids with meteorological stations are selected as comparison samples. For validating whether these stations are randomly distributed within each grid, the grid is divided into two parts by the station elevation: (1) the area in which the elevation is greater than that of the station and (2) the area in which the elevation is lower than that of the station. We illustrate the histogram of the proportions of the first part in all these grids (Fig. 6). It is obvious that the elevation areas divided by these stations are randomly distributed about 0.5 in all grids, thus confirming that these stations are randomly distributed in the selected region.

Fig. 6

Histogram of the area proportions of the partial regions in which the elevation is higher than the station elevation.

Then, the mean SD value at these grids in the subjunctive region is regarded as the “real” SD. For comparison, the evaluated SD values in each grid from the assimilated results are averaged. The comparison results are illustrated in Fig. 7.

Fig. 7

Comparison between simulated and observed results. The observation is the mean observed SD of select stations and the values on the $y$ axis are the mean simulated SD of corresponding grids. The red cycles denote the assimilated results and the gray circles denote the AMSD results. From years 2000 to 2008, $R^2$ is 0.10, 0.12, 0.49, 0.11, 0.35, 0.54, 0.53, 0.32, and 0.72, respectively, Pearson $r$ is 0.73, 0.66, 0.81, 0.61, 0.79, 0.79, 0.82, 0.76, and 0.87, respectively.

The assimilated SD and AMSD data are compared to the in situ observations. The RMSE is used to quantify the accuracy. The RMSE of the assimilated results are in the range [0.002 m, 0.008 m], and that of the AMSD results are in the range [0.004 m, 0.021 m]. In all the years, the accuracy of the simulation is improved, most especially in 2005, 2006, 2007, and 2008. This suggests that the assimilated results become closer to the in situ observations and that the evaluating accuracy is improved by using the synthesis method suggested in this study. We investigated the determination coefficients and the Pearson product-moment correlation coefficients ($r$) for these nine years. The range of the determination coefficients is in the range [0.10, 0.72], and the range of Pearson correlation coefficients is in the range [0.61, 0.87]. These coefficients indicate there are stronger correlations between the assimilated results and the in situ observations. It should be noted there are some low values in the determination coefficients. Such as, the $R^2$ is 0.10 in 2000. This suggests there are some incorrect estimations in the results, although the Pearson $r$ is 0.73.

There are noticeable cases in which the assimilated SD and AMSD are close to zero whereas the observations are much greater than zero. The reason for this discrepancy is that remotely sensed data cannot provide the snow information in these cases. Shallow snow cover and wet snow hardly attenuate microwaves and do not give the volume scattering necessary for snow volume determination.⁴⁴ This explains why snow cover is not monitored by RS while there is snow or even deep snow at the stations. When the assimilated SD is close to zero, this means that SCF or AMSD is close to zero also, because the SD observation used in the assimilation strategy is from SDC transferring or AMSD. Although better simulation results are obtained by using an assimilation strategy, combining information from different RS data and meteorological data, the assimilated results are also apparently hampered by the accuracy of observations. Improvement of retrieval algorithms to obtain microwave SD in rugged mountain areas is needed for a better simulation of the snow distribution.

3.3.

Factors Influencing Evaluation Results

Snowpack in the TP is patchy and shallow and changes rapidly in most cases. There are complex terrain conditions and a varied distribution of vegetation in the TP. With the consideration of data availability and the features of the snow distribution in the TP, a series of empirical methods are used in the study. Topographic features and snow distribution characteristics, which affect the RS data availability and snow processes, are the most important reasons why we use these methods. There are three important steps in our study: downscaling of microwave SD data, transferring between SCF and SD using the SDC, and simulation of snow processes. Here, we discuss how characteristics of the snow distribution and topography influence the evaluating results in these steps.

When the microwave SD data are downscaled, SD at the subgrid scale of 0.1 deg is evaluated by the proportional weight of SCF to the whole snow-covered area in the grid with a resolution of 0.25 deg. The implicit assumption here is that the SCF of a grid is directly proportional to the SD. Obviously, the assumption varies over different terrains. However, we have not adopted further improvement in the combination method. We did attempt to construct an SDC relation between SCF and microwave SD, but found that the goal could not be achieved. There is only a very weak correlation between the two. We suppose the most important reason for this failure is that the snowpack is shallow and that the terrain is mountainous. In some of the literature, it has been suggested that microwave SD values have large uncertainties in patchy and shallow snowpack regions.⁴⁴ However, there are not enough data to support reliable prior relations between SD and SCF in different terrain conditions. In fact, the SCF is patchy in most cases. In these cases, we use an empirical SDC method to transfer SCF to SD. The combined SD is regarded as an observation operator in the assimilation framework only in cases in which SCF is greater than 50%. Although the different terrain conditions could affect the downscaling accuracy of SD, we believe the present method to be reasonable given the limitation of data and available information.

In the cases in which the snow distribution is patchy, we construct a general SDC relation between SD and SCF. Of course, the most reliable method is to construct different regional SDCs under different topographic conditions. However, there are few data to construct these relations. Although there are 94 stations with snow observations in the study region, there are almost no continuous observations nor enough data to support the regression of SDC under different conditions because of the short snow period and missing observations. Therefore, we suggest a compromise in this study, in which the SDC is constructed based on the mean in situ SD of all stations. There are good correlations between the averaged in situ SD and SCF, indicating that the SDC concept, used in many other studies,^45,46 is also suitable to the TP region. Although there are some calibrated SDCs in other study regions with good data,⁴⁷ those calibrated parameters could not be used directly in the TP, because the snow in this region has features that differ from those of other regions. Liston et al.⁴⁸ also indicated that the key parameters will be dependent upon different topographic scales. With only SCF assimilation, Lannoy et al.²⁸ indicated that the simulated timing of the onset of the snow season was improved, but without an improvement of SWE amounts. There are many reasons for the result, but mostly the complex relationship between SCF and SWE is responsible. Since we could not calibrate regional SDC in the TP based on the existing observations, a general SDC is reasonable as a substitute. The results show relevant simulation results that are close to the truth.

Modeling of snow processes is also influenced by topography. In our study, the forcing data have a spatial resolution of 0.1 deg, so the corresponding modeling unit is also 0.1 deg. In the TP region at this scale, the terrain is very complex, with features such as meadows, forests, and so on. In this study, snowmelt is regarded as uniform because of the limitation of computation speed and forcing data. Accurately examining the snow simulation at the subgrid scale is very important. With different parameterized considerations of inhomogeneous snow processes at the subgrid scale, the simulations are expected to lead to very different results. This subject is of interest for further study. A method to take into account computation performance and inhomogeneous snow processes at the subgrid scale is also needed in the future.