Reconstruction of Snow Depth Data at Moderate Spatial Resolution (1km) from Remotely Sensed Snow Data and Multiple Optimized Environmental Factors: A Case Study over the Qinghai ‐ Tibetan Plateau

: Snow depth distribution in the Qinghai ‐ Tibetan plateau is important for atmospheric circulation and surface water resources. In ‐ situ observations at meteorological stations and remote observation by passive microwave remote sensing technique are two main approaches for monitoring snow depth at regional or global levels. However, the meteorological stations are often scarce and unevenly distributed in mountainous regions because of inaccessibility, so are the in ‐ situ snow depth measurements. Passive microwave remote sensing data can alleviate the unevenness issue, but accuracy and spatial (e.g., 25 km) and temporal resolutions are low; spatial heterogeneity in snow depth is thus hard to capture. On the other hand, moderate imaging at moderate spatial resolution (1 and high temporal resolution but area Fusing passive microwave snow depth data with optical snow area extent provides an unprecedented opportunity for generating snow depth data at moderate spatial resolution and high temporal resolution. In this article, a linear multivariate snow depth reconstruction (LMSDR) model was developed by fusing multisource snow depth data, optical area extent data, and environmental factors (e.g., spatial distribution, terrain features, and characteristics), to reconstruct daily at a We found that a significant impact on the of depth over the QTP. Relatively high accuracy (root mean square error (RMSE) = 2.26 cm) was observed in the snow depth when compared with in situ data. Compared with the passive microwave remote sensing snow depth product, constructing a nonlinear snow depletion curve product with an empirical formula and fusion snow depth product, the LMSDR model demonstrated a significant improvement in accuracy of snow depth reconstruction. The overall spatial accuracy of the reconstructed snow depth was 92%. Compared with in situ observations, the LMSDR product performed well regarding different snow depth intervals, land use, elevation intervals, slope intervals, and SCD and performed best, especially when the snow depth was less than 3 cm. At the same time, a long time snow depth series reconstructed based on the LMSDR model reflected interannual variations of snow depth well over the QTP.


Introduction
As an important constituent of the cryosphere and an indispensable variable used in hydrological science research, snow cover has a significant impact on global climate change and the hydrological cycle [1]. The relatively high levels of albedo and thermal insulation of snow cover directly affect the energy exchange between the atmosphere and the earth surface [1,2]. In the global hydrological cycle, snow accumulation and melting contribute greatly to water resource redistribution. It also provides the most important resources of freshwater in arid and semi-arid regions [3,4]. As the basic physical characteristics of snowpacks, snow depth is broadly applied in climatic and hydrological simulations [5][6][7]. Meanwhile, snow depth plays an important role in weather forecast [8], snowmelt runoff simulation [1], and drought and flood prediction [9,10]. Since the Qinghai-Tibetan Plateau (QTP) forms the target port of the space in Asia, its massive distribution of ice and snow considerably impact the ecosystem of the region and those of Asia and the Northern Hemisphere [11]. As a major open water storage in cold regions, snow cover buildup and decay are closely related to the hydrological cycle [12], the atmosphere-land interaction [13], the green-up date of vegetation in QTP [14], and even the frequency of summer monsoon in East Asian and heatwaves in northern China [15][16][17][18]. Knowledge of snow depth is essential to obtain accurate spatio-temporal variations of snow water equivalent for management of water resources and understanding plateau climate change.
Although real-time snow depth data obtained at conventional meteorological stations and automatic stations make it possible to investigate long-time series snow depth variations [19,20], in-situ measurements reflect merely snow depth variations in a limited geographical area. Besides, meteorological stations in QTP were located primarily in its eastern and low elevation areas; their applications in hydrology are thus limited [21]. However, remote sensing technologies provide long-time series measurements of snow on large scales [22][23][24]. In the visible and near-infrared (VNIR) spectral regions, satellite sensors such as the moderate resolution imaging spectroradiometer (MODIS), advanced very high-resolution radiometer (AVHRR), sentinel-2, etc. can obtain snow cover extent and snow cover fraction (SCF) with relatively high spatial resolution [25][26][27]. Snow cover extent products have been widely applied to hydrological models for streamflow forecast in many catchments in the world [28][29][30][31][32]. However, optical sensors only measure snow coverage but not snow depth; also, they are susceptible to cloud interference. Cloud and snow have similar spectral signatures in VNIR, leading to misclassification between snow cover and cloud [33,34]. Passive microwave sensors such as SMMR (scanning multichannel microwave radiometer), SSM/I (special sensor microwave/imager), SSMI/S (special sensor microwave imager/sounder), and AMSR2 (advanced microwave scanning radiometer 2) were based on the difference between the depth and shallow of snow cover and the microwave brightness temperature. The snow depth can be inferred by the brightness-temperature difference of GHz 18-37 [35][36][37] have the capacity to map snow depth on the global scale because microwaves can penetrate through cloud [38,39], but it is difficult to identify areas of shallow snow for their low spatial resolution and snow depth is often overestimated in arid and frigid areas [21].
Accuracy or temporal-spatial resolution is often limited for a single data set by a single technology. Reconstruction of accurate snow depth at high resolution by fusing multiple data sets developed through multiple remote sensing technologies becomes attractive. Assimilation of in-situ observations and snow depth data of remote sensing had proven to be effective in improving hydrological models and snow depth models [40,41]. Although the ensemble Kalman filter method (EnKF) is widely used for snow depth data assimilation, snow depth data obtained through EnKF has hardly been used effectively, given the few and uneven distribution of meteorological stations and large fluctuation in surface elevation [42,43]. Snow depth simulation based on the remotely sensed SCF data involves building a nonlinear snow depletion curve, which is hard to accomplish because of its high nonlinearity and spatial heterogeneity of QTP, the accuracy of the simulation is often compromised [44]. A downscaling algorithm of multiple data sources was developed based on passive microwave remote sensing data and optical remote sensing data [45]. New, improved snow depth downscaling methods have been developed and can generate good simulations of snow depth in Alaska, QTP, and northern China [24,46,47]. However, terrain features can affect surface spectral reflectance and scattering of microwave signals, leading to low accuracy in the microwave snow depth dataset; meanwhile, MODIS data are severely subject to cloud interference, leading to errors in the final snow depth product [21,47]. Snow depth downscaling modeling has been used in simulating snow depth variations in QTP with some success [48]. However, using a combination of multichannel brightness temperature data with surface information to estimate the spatial distribution of snow depth can either overestimate or underestimate snow cover extent to a certain degree.
On the other hand, environmental factors such as longitude, latitude, slope, surface roughness, and surface elevation fluctuation often affect snow accumulation at a specific location. The objective of this study is to develop a linear multivariate snow depth reconstruction (LMSDR) model, incorporating the fused snow depth from multiple remote sensed data sets, in-situ snow depth observations, and relevant environmental factors controlling snow accumulation at a specific location to reconstruct snow depth in QTP with high accuracy and high spatial resolution.

Study Area
QTP is located in the west of China between 26°00′-39°47′N and 73°19′-104°47′E. It starts from the southern foot of the Himalayas in the south and ends at the northern side of Kunlun Mountains and Qilian Mountains in the north, with a span of about 1532 km. The distance from Pamir Plateau in the west to Hengduan Mountains in the east is about 2945 km. The QTP is high in the northwest and low in the southeast (Figure 1), with an average elevation of about 4000 m and an area of 2,620,100 km 2. QTP is the highest and largest plateau in the world. The area features long and frigid winter and short and warm summer. There is a great diurnal temperature variation. The mean annual temperature ranges from −15 °C to 10 °C. Precipitation mainly occurs in May through September but is extremely spatially uneven [49,50]. Snow exists primarily from September to March the second year. Major rivers in QTP include the Yangtze River, the Yellow River, the Lantsang (Mekong) River, the Brahmaputra, the Ganges, and the Salween that all originate from the plateau. The QTP also has numerous lakes, snowpacks, and glaciers.

Long-Term Passive Microwave Snow Depth Dataset (PMSD)
The passive microwave snow depth dataset (PMSD) used for this study was developed by Che et al. [51][52][53] based on satellite passive microwave measurements and provided by the national Tibetan plateau data center (TPDC) (https://data.tpdc.ac.cn). The daily snow depth data has a 25 km spatial resolution. The raw data used in generating the PMSD data set were the daily passive microwave brightness temperature data from SMMR (1979SMMR ( -1987, SSM/I , and SSMI/S (2008-2018) archived by the US national snow and ice data center (NSIDC). The PMSD data set has been widely applied to climate change analysis, hydrological simulation, water resource management, vegetation variation, etc. [54][55][56]. In this study, we used the nearest-neighbor interpolation method for the raw data resampling [48,57]. Then through preprocessing such as projection transformation and cropping, the snow depth data with a 1 km spatial resolution were generated as inputs for snow depth reconstruction at high accuracy.

In-Situ Snow Depth Observations
In-situ snow depth observations include those at conventional and automatic meteorological stations. In-situ data from conventional meteorological stations was provided by the national meteorological information center (http://data.cma.cn/site). Data included the date of observation, longitude, latitude, elevation, and snow depth (minimum value: 1 cm). In-situ data from the automatic meteorological stations were provided by TPDC (https://data.tpdc.ac.cn). Data included the date, snow depth, soil moisture, and land surface temperature [63].

Other Data Sources
To improve the accuracy of the snow depth reconstruction (SDR) model, other data sets used for model development include digital elevation model (DEM) data, land-use and land-cover change (LUCC), and soil particle size distribution data. National aeronautics and space administration (NASA) SRTM DEM at 90 m spatial resolution of version 4 (SRTM 90 m V4) [64] were processed to generate slope, aspect, surface roughness, and surface elevation fluctuation datasets. LUCC data provided by resource and environment science and data center (http://www.resdc.cn) include cultivated land, forest, grassland, waterbody, construction land, and unused land. These six types of LUCC have 25 subtypes. TPDC (https://data.tpdc.ac.cn) also provided soil particle size distribution data that was a multi-layer soil particle size distribution dataset developed based on 1:1,000,000 soil maps and 8595 soil profiles obtained during the second national soil survey of China [65,66]. The contents of sand and clay in the surface layer (0-30 cm) were used to evaluate the result of snow depth reconstruction.

Development of the LMSDR Model
Given the limitation of a single data source by a single remote sensing technology and the success of data fusion of multiple remote sensed data sets, we proposed developing a linear multivariate snow depth reconstruction (LMSDR) model to reconstruct snow depth in QTP at high accuracy and spatial resolution (1 km). The strategy for the model development was as follows: firstly, the algorithms of data fusion and downscaling were used to produce fused snow depth (FSD) data based on long-term passive microwave snow depth dataset (PMSD) and MODIS-SCF data. Secondly, the FSD was combined with in-situ observations to build a multifactor linear model to reconstruct snow depth data using relevant environmental factors such as longitude, latitude, slope, surface roughness, surface fluctuation, and SCD. Finally, the LMSDR model was used to simulate a long-time series of snow depth (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) in the QTP, and simulation accuracy was then evaluated.

Multisource Snow Depth Data Fusion
The algorithm of multisource data downscaling was adopted to downscale the PMSD product from a 25 km spatial resolution to 1 km. Three cases were considered: (1) when the value of a PMSD pixel (25 km × 25 km) was greater than 0, and at least one of its SCF sub-pixels (1 km × 1 km) had pixel values greater than 0 (SCFi > 0), the snow depth values of these sub-pixels were assigned based on the fusion-downscaling algorithm (see Equation (1a)); (2) when the value of a PMSD pixel was greater than 0, and all its SCF subpixels had a pixel value of 0 (SCFi = 0), the snow depth values of its sub-pixels were assigned to 0 (see Equation (1b); and (3) when the snow depth value of the PMSD pixel was 0, and at least one of the SCF sub-pixels had a value greater than 0 (SCFi > 0), the snow depth values of the sub-pixels were assigned according to the relationship between the perennial snow depth observations and SCF [47] (see Equation (1c)). These three cases were summarized in Equation (1) if where SCFi is the SCF value of the i-th SCF sub-pixel within a PMSD pixel; PMSD is the snow depth value of the PMSD pixel; FSDi is the fused (assigned) snow depth value for the i-th SCF sub-pixel; n is the total number of the sub-pixels within a co-registered PMSD pixel (here the value of n is 625 when 25 km is downscaled to 1 km).

Snow Depth Reconstruction (SDR) Model Development and Evaluation
We utilized 14,520 samples selected from eight even hydrological years of in-situ observations, snow cover characteristics, and environmental factors (e.g., spatial distribution, terrain features, etc.) to develop an SDR model. The tentative candidate models were as follows where y is the in-situ observation at the meteorological stations; X1, X2, X3… Xi are snow characteristics and environmental factors; β1, β2, β3 … βi are the model coefficients corresponding to the factors; ui is the error. Equations (2-5) represent linear, logarithmic, power, and inverse models, respectively. We utilized the 14,520 in-situ observations selected from eight even hydrological years for model construction. We evaluated the models using the sample leave-one-out cross-validation (LOOCV) method [67,68]. Firstly, one sample from the in-situ observations was chosen as the verification data, and the remaining samples were then used for model development. Secondly, the accuracy of the model was evaluated using the sample that was excluded in advance. This process was repeated for all data samples (14,520 in this case). Finally, the performance of all models (Equations (2-5)) was evaluated using the mean coefficient of determination (R 2 ) and the root mean square error (RMSE). The best model was then selected.
Ground truth data used for model validation included 148,274 snow depth data from the 100 meteorological stations within the eight odd hydrological years and the 10 automatic meteorological stations for the period 2015-2016. The accuracy analysis for snow depth products of PMSD, SDC, FSD, and SDR was performed by evaluating model accuracy matrices such as RMSE, mean absolute error (MAE), positive mean error (PME), and negative mean error (NME) that are defined as below where yi is the in-situ observation at meteorological stations; E (yi) is the i-th estimated snow depth value; y is the mean observed value; n is the total number of surface observations; and p and r are samples with an observed value more and less than the estimated values, respectively. Figure 2 shows the procedure and methods used for the SD model development.

Selection of Environmental Factors
The selection of environmental factors that affect the multivariate SDR models is a very important step to optimize the model sensitivity. The LOOCV approach and linear regression method were used to evaluate factors that affect snow depth. Candidate factors for optimization include location, terrain features, LUCC, soil particle size proportion, and SCD. Table 1 shows the regression analysis results between in-situ snow depth observations and each candidate factor. We can see that a single factor generally has a poor correlation with the observations, with R between 0.07 and 0.69. It showed that the highest correlation (R = 0.69) occurred between SCD and in-situ observations, and the lowest correlation (R = 0.07) occurred between aspect and in-situ observations. We chose six environmental factors with a correlation coefficient larger than 0.3 to develop the final SDR model: latitude, longitude, slope, surface roughness, surface elevation fluctuation, and SCD.

Model Optimization
The six optimized environmental factors were thus used in combination with FSD to build a final multivariate reconstruction model. As shown in Table 2, the LMSDR model (R 2 = 0.63, RMSE = 2.28 cm) has higher accuracy and stability than the power model (R 2 = 0.63, RMSE = 6.94 cm), the logarithmic model (R 2 = 0.59, RMSE = 6.89 cm), and the inverse model (R 2 = 0.58, RMSE= 6.84 cm). Previous studies verified the presence of a highly linear relationship between microwave brightness temperature and snow depth over homogeneous snow-covered areas [69]. However, microwave brightness temperature data have a poor spatial resolution (e.g., 25 km); a single brightness temperature pixel can contain multiple types of land cover. Topographical complexity enhances the heterogeneity in the spatial distribution of snow depth. These factors greatly limit the spatial expression of microwave brightness temperature data [70,71]. The SDR data at higher spatial resolution (1 km) can alleviate this weakness effectively by introducing MODIS-SCF data of moderate high resolution (1 km) and environmental factor data at high resolution (90 m). Figure 3 shows the spatial distributions of snow depth in QTP generated by the SDR models. The LMSDR model generated better snow depth spatial distribution than the logarithmic model, the exponential model, and the inverse model. Meanwhile, we tried a similar data analysis by introducing all other factors without SCD to simulate snow depth and found that the R 2 = 0.56 and RMSE = 2.44 cm (the LMSDR model, R 2 = 0.63, RMSE = 2.28cm.) and that the introduction of SCD improves the accuracy of the reconstructed snow depth. In summary, the LMSDR model better reflected the snow depth simulation effect and its spatial distribution over the QTP. Based on the LMSDR model, the daily snow depth data of 16 hydrological years (2002-2018) were reconstructed.

Overall Accuracy Comparison
The daily snow depth samples (148,274) of in-situ observations in the eight odd hydrological years of 100 meteorological stations, ten automatic meteorological stations from 2015 to 2016, and the products of LMSDR, PMSD, SDC, and FSD extracted from the matched 110 stations were used for accuracy analysis. Based on the modeling accuracy matrices of RMSE, MAE, PME, and NME (Equations (7-10)), the overall accuracy of the four snow depth products (PMSD, SDC, FSD, and LMSDR) were compared and analyzed using the in-situ observations as the reference (Figure. 4a). The LMSDR model predicted snow depth with RMSE and MAE being 2.26 cm and 0.46 cm, respectively, followed by the SDC product (RMSE = 2.28 cm, MAE = 0.46 cm), then by the products of PMSD and FSD. The values of PME (4.17 cm, 3.91 cm, and 4.00 cm) were smaller than the absolute values of NME (−4.33 cm, −7.75 cm, and −11.63 cm) for the products of PMSD, SDC, and FSD, respectively. These accuracy matrix data indicated that the products of PMSD, SDC, and FSD underestimated the observed snow depth. In comparison, the absolute values of PME (3.19 cm) and NME (-3.90 cm) of snow depth data generated by the LMSDR model were close in magnitude, both smaller than those of the other three snow depth products (PMSD, SDC, and FSD), indicating that the LMSDR product outperformed the PMSD, SDC and FSD products in accuracy. Figure 4b shows the accuracy comparison of the four snow depth products against in-situ observation for different snow depths to investigate the performance of each product for snowpack of different depths. About 95% of the in-situ observation samples had a snow depth of less than 3 cm, and about 1% of the samples had a snow depth of more than 30 cm. Overall, the RMSE values of all snow depth products increase with increasing snow depths. Products of PMSD, SDC, FSD, and LMSDR performed relatively well in areas of shallow snow (depth less than 3 cm), with RMSE values of 3.19 cm, 1.57 cm, 5.55 cm, and 1.52 cm, respectively. The newly developed LMSDR product performed the best in accuracy for all snow depth products. As shown in Figure 4, although the LMSDR product outperformed all other products for snowpack of different depths, it was noteworthy that poorer accuracy occurred when the value of snow depth was more than 30 cm. The reasons for this are twofold: (1) the scarcity of in-situ snow depth data in western QTP with deep snow was available for model development, and (2) the spatial heterogeneity of snow cover distribution in the topographically complex regions. The accuracy of reconstructed snow depth can be subject to many factors in the study area, e.g., patchiness prevalence, snowdrift frequency, wide distribution of frozen soil, the undulation of terrain, and the difference in the underlying surface of snow cover. The uneven distribution of in-situ observations in the deep snow areas can contribute greatly to the low accuracy in the deep snow areas [72], which can be alleviated with more meteorological stations being installed in western QTP. At the same time, it was found that the accuracy of the FSD product by the fusion algorithm performed not so well as the PMSD and SDC products. There may be two causes: (1) the fusion-downscaling algorithm assigned snow depth to each subpixel linearly using the SCF value as a weight and the PMSD value as invariable to get a good downscaling effect. However, the complex local terrain resulted in a discontinuous distribution of snow, hence the uncertainty of the fusion algorithm [47,73,74]; and (2) extreme heterogeneity of snow cover distribution with prominent patchiness and low spatial resolution (25 km)of the PMSD product resulted in relatively great errors in FSD in areas of high altitudes.

Snow Depth Accuracy Comparison of Different Snow Characteristics and Environmental Factors
Further performance analysis of the four data products was performed for snowpack over different LUCCs, elevation intervals, slope intervals, and SCDs. Table 3 shows that all products compared better with in-situ measurements over cultivated land than other LUCC types except PSD, with the LMSDR product performing the best (RMSE = 1.58 cm). For the construction land and bare land, the LMSDR product performed better than the other three products in matching the in-situ observations, with the RMSE being 1.86 cm and 1.90 cm, respectively. The LMSDR model performed worst over grassland (RMSE = 2.28 cm), but still better than the other three products. The SDC product outperformed the LMSDR product for snow over the forest, probably because of a lack of enough Elevation distribution strongly affected the amount and distribution of precipitation and snow accumulation at specific locations. As shown in Table 3, the PMSD, SDC, FSD, and LMSDR products performed relatively well in matching in-situ observations at the elevation zone of 1000-2000 m, with RMSE being 1.48 cm, 1.16 cm, 3.11 cm, and 1.13 cm, respectively. The LMSDR product performed the best in accuracy. All four products have relatively lower accuracy in the region of elevation above 4000 m. All products are positively correlated with elevation since the probability of precipitation as snowfall increases with increasing elevation [75,76]. The SDC product was found to perform the best in the region of elevation between 2000-2500 m (RMSE = 0.87 cm). The FSD product was found to perform the worst in regions of elevation between 3500-4000 m (RMSE = 13.24 cm). The differences in accuracy among the products were closely related to sample distribution at different elevation zones as well as the interaction between local climate and the complex terrain [77,78].
Slope was an important environmental factor of terrain features. The local incident angle and the polarization direction of microwave radiation depends on the slope. Slope has an impact on snow cover distribution to some extent by changing solar radiation and water vapor channel [48,74,79]. As shown in Table 3, the RMSE of the four snow depth products increased with increasing slope. The PMSD, SDC, FSD, and LMSDR products showed relatively higher accuracy in relatively flat areas with a slope less than 5°, with the RMSE values being 2.98 cm, 2.30 cm, 7.07 cm, and 1.80 cm, respectively.
As an important indicator of snow cover longevity, SCD had a prominent positive correlation with snow depth. Normally, areas with longer snow cover duration had a persistently lower temperature, and snow could be accumulated for a longer period [80]; it takes a longer time for a deeper snowpack to melt completely. However, as SCD increases, the accuracy of the PMSD, SDC, FSD, and LMSDR products exhibited a decreasing trend, as shown in Table 3. When SCD was less than 30 days, all snow depth products performed relatively well in accuracy, with the LMSDR product being the most accurate (RMSE = 1.18 cm). When SCD was more than 60 days, the PMSD, SDC, FSD, and LMSDR products performed poorly in accuracy, with RMSE being 3.88 cm, 2.71 cm, 8.38 cm, and 2.50 cm, respectively.

Spatial Accuracy Comparison of Different Snow Depth Products
The PMSD product fused with the moderate-resolution MODIS-SCF data can better capture snow depth heterogeneity than PMSD itself, reducing errors caused by the low resolution-related mixing problem [47]. The FSD data produced by the fusiondownscaling algorithm using the PMSD data and the cloud-reduced MODIS-SCF datasets effectively improved the accuracy of snow depth distribution in QTP while preserved snow depth information in the PMSD products.
Taking the spatial distribution of snow depth on 2 February 2013 as an example ( Figure 5), parts of the typical snow-covered area were presented in detail to reflect the spatial distribution of the four snow depth data products in QTP. As shown in Figure 5a, the PMSD product reflected roughly the spatial distribution of snow depth, but it is hard to effectively discriminate the snow-covered from snow-free areas. Also, the maximum snow depth of the day observed by the PMSD product was 23 cm, while the in-situ observed maximum depth was 93 cm. The spatial resolution of 25 km limited the details of the snow depth spatial distribution. As shown in Figure 5b, the SDC product presented the spatial distribution of snow depth variation in more detail at a resolution of 1 km. However, as the number of pixels with a 100% SCF increased in the winter season, snow depth data generated by the empirical formulas failed to discriminate the area of deep snow from that of shallow snow. Compared with the mean values of the in-situ observations, the snow depth in the SDC product was obviously overestimated in the central QTP. In the perennially snow-covered areas such as the Himalayas and Nyainqentanglha Mountains, the SDC value was significantly underestimated compared to the in-situ observations, with the maximum value being only 26.9 cm. Figure 5c shows the spatial distribution of snow depth of the FSD product, from which area of deep snow can be discriminated effectively from the area of shallow snow. Nevertheless, because of the complex local terrain, the linear mixing approach based on SCF as a weight had some uncertainty, which resulted in serious heterogeneity and fragmentation of pixel distribution in the snow-covered area in the mid-eastern QTP. Figure 5d shows the snow depth distribution of the LMSDR product. Compared with the PMSD, SDC, and FSD products, the LMSDR product presented snow depth spatial distribution in detail by introducing impacts on snow depth by multiple environmental factors. The areas of deep snow shown in the LMSDR product were located primarily in the Karakoram Range, the west of the Kunlun Mountains, the Himalayas, the Kailas Range, and the Nyainqentanglha Mountains, with snow depth exceeding 100 cm. The areas of shallow snow were located primarily in the northern QTP, Ngari, the eastern mountains of the Western Sichuan Plateau, and the northern Hengduan Mountains. Only a little snow was seen in the eastern part of Qinghai Province, the southern part of western Sichuan Plateau, the southern part of Hengduan Mountains, the Qaidam Basin, and the lower Brahmaputra in southeastern Tibet. These observations derived from the LMSDR product were consistent with previous studies [81][82][83]. Taking the in-situ observations from the 100 meteorological stations as a reference to estimate the spatial accuracy of the four snow depth products quantitatively. When the snow depth from each product and the in-situ observations at the corresponding station within the same time period match, the error matrix was defined as true; otherwise, the error matrix was defined as false. By calculating the ratio of the counts of the samples that each product matches the in-situ observation to the total sample counts, the spatial accuracy for the product was derived in the study area. Results show that the PMSD, SDC, FSD, and LMSDR products had a spatial accuracy of 66.8%, 92.5%, 92.1%, and 92%, respectively. Overall, the SDC, FSD, and LMSDR products had high spatial accuracy, while PMSD has poor spatial accuracy. PMSD products had a better spatial accuracy of 70%-90% in the northeast and southwest of QTP, but less than 50% in the southeast and hinterland of QTP (Figure 6a). The overall spatial accuracies in the SDC, FSD, and LMSDR distributions were similar. Most of the spatial accuracy in the northeast, southeast, and southwest of QTP was in the 70%-95% range. The regions with lower spatial accuracy are mainly distributed in the central and eastern regions of QTP (Figure 6b-d). In order to show spatial accuracy at different geographic zones, the meteorological stations 52,869,52,908,55,228,55,655,56,128, and 56,173 (see Figure 6) were selected for comparison (see Table 4). Results show that the LMSDR, SDC, and FSD products performed relatively well in accuracy in the shallow snow areas; for example, spatial accuracy exceeded 93% at station 55,228 and station 56,128 (Table 4). In the areas with deep snow, spatial accuracy was lower for the LMSDR, SDC, and FSD products. For example, at station 55,655, the accuracy was less than 70% ( Table 4). The PMSD product at station 52,908 and station 52869, located in the northern QTP (Figure 6), performed the best with accuracy higher than 80%, followed by station 55,228 and station 56,173 with an accuracy of 65% and 61%, respectively. The PMSD product performed the worst at station 55,655 and station 56,128, with accuracy being 36% (Table 4). Through the spatial accuracy analysis of the snow depth in QTP, a significant correlation between the spatial accuracy of snow depth and the values of snow depth were revealed. The products performed prominently better in spatial accuracy in areas of shallow snow than in areas of deep snow, which may be due to a relatively large difference between the in-situ observations at the meteorological stations and those by the remote sensing satellite. For the LMSDR model, the spatial distribution of the meteorological stations also affected the accuracy of the snow depth products since the in-situ observations have been dominantly located in eastern QTP, where snow was generally shallower. For instance, snow measurements at stations 55,655 and 56,128 in the valleys at the southern foot of the Himalayas and the Brahmaputra Basin were affected by the warm and humid air from the Indian Ocean and the Bay of Bengal [81,82], widespread occurrence of cloud and frequent rainfall caused a relatively low spatial accuracy in snow depth. At station 52,908 in the central QTP and station 55,228 in Ngari in the north of the Himalayas ( Figure 5, Figure 6), these products performed better in snow depth because of less rainfall and cloud interference. In summary, compared with the PMSD, SDC, and FSD products, the LMSDR product generated by the LMRSD model developed in this study performed better in the overall spatial accuracy when compared to in-situ observations.

Interannual Variation of the LMSDR Model
Based on the daily snow depth over the QTP in the 16 hydrological years from 2002 to 2018, the performance of the LMSDR product was further evaluated by analyzing its interannual and seasonal variation of the annual and seasonal mean snow depths. As shown in Figure 7a, the LMSDR product could well capture the variation trend of snow depth over the QTP when compared to the in-situ measurements. Figure 7b-d showed the annual variation between the mean reconstructed snow depth and the in-situ observations for autumn, winter, and spring, respectively. The LMSDR product underestimated the snow depth in winter and spring, especially in spring, with a difference near 0.2 cm. For some years, the deepest snow from the LMSDR product occurred in winter, while the deepest snow from the in-situ observation occurred in spring. In autumn, the LMSDR product was generally consistent with the in-situ observation. In some years, such as the autumn of 2012 and winter of 2014, the LMSDR product showed an overestimate. Generally, the LMSDR product worked best for snow depth reconstruction in autumn, followed by spring, and had relatively low accuracy in winter. In particular, in the years 2014 and 2015, the differences in the mean snow depth between the reconstructed and insitu observed were more than 1.2 cm. The reconstruction of a long-time series of snow depth data using the LMSDR model developed in this study could effectively reduce some errors caused by excessive fluctuations in snowfall in some years over the QTP. The interannual variation of the annual mean value of the snow depth in the hydrological years was generally consistent with the in-situ observations. However, during the analysis of the interannual variations in different seasons, it was found that the accuracy in snow depth reconstruction in winter by the LMSDR product was relatively low. Moreover, the reconstruction in spring showed underestimation. As a key variable for snow depth reconstruction, the PMSD product used for reconstruction in winter and spring when snow was deep reflected the cumulative amount of all previous snowfall, while the in-situ observed snow depth provided by the meteorological stations represented a daily accumulation. The two values of snow depth were quite different, and the results of the reconstruction model showed obvious errors when snow was deep. With continuous snow accumulation in winter and spring, the possibility of bias and underestimation of the reconstructed values of snow depth was relatively high. On the other hand, seasonal winds significantly affected snow cover extent in higher elevation zones during winter and spring [84,85]. At the same time, the difference in the measurement timing between satellite sensors and in-situ observations, and the influence of various surface environmental factors increased the likelihood of such deviations.

Conclusions
Accurate measurement of the dynamic snow depth over mountainous areas such as QTP is important to understanding alpine climate dynamics, hydrological cycle, production and life, and ecological changes. In this study, multiple data sets of long-time series of snow depth such as the PMSD and snow cover data such as MODIS-SCF were used for an optimized multivariate reconstruction model development and validation. We took the in-situ observations of snow depth as "true value" for model development and validation. Such a model can effectively reduce accidental errors caused by excessive or low snowfall in some years. At the same time, the effective fusing of PMSD of low resolution (25km) and MODIS-SCF of moderate resolution (1 km) improved the spatial accuracy and resolution (1km) of snow depth distribution and avoided the spatial loss caused by simple data fusion or single data source identification in most studies. Based on this, the introduction of snow cover day and optimized influencing environmental factors such as longitude, latitude, slope, surface roughness, and surface elevation undulation made the LMSDR product outperform several other snow depth products when in-situ snow depth data were taken as the reference. As a result, the RMSE and MAE of the reconstructed snow depth by the model were only 2.26 cm and 0.46 cm, respectively; when compared to in-situ measurements, snow spatial accuracy exceeded 92%. The RMSE of the LMSDR product was 1.52 cm in areas of shallow snow where snow depth was less than 3 cm. The mean annual snow depth in hydrological years provided a sound representation of the trend of snow depth over QTP, with the best reconstruction for autumn in the years. Despite the success of the LMSDR model in reconstructing the longterm snow depth data in QTP, there were still some deficiencies that require improvement. For instance, due to the uneven spatial distribution of the meteorological stations in QTP, the snow depth reconstructed in areas of deep snow that occurs more frequently at elevation areas, low valleys, and forests had high uncertainty and were hard to evaluate because of lack of in-situ observations. Long-term in-situ observations in these areas are important for future model improvement and accuracy of reconstructed snow depth.
Author Contributions: All authors have made significant contributions to the manuscript. P.W. conceived the study. T.Z. designed the study. P.W., T.Z., and X.Z. wrote the manuscript. X.Z. and G.Y. helped to analyze the results and improved the manuscript. P.W. and N.W. processed the data. J.L. and B.W. helped revise some figures and the manuscript. All authors read and approved the final manuscript.

Data Availability Statement:
The data used to support the findings of this study are available from the corresponding author upon request.

Acknowledgments:
The author thanks anonymous reviewers for providing invaluable comments on the original manuscript.

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