Snow Depth Estimation and Spatial and Temporal Variation Analysis in Tuha Region Based on Multi-Source Data

Abstract

In the modelling of hydrological processes on a regional scale, remote-sensing snow depth products with a high spatial and temporal resolution are essential for climate change studies and for scientific decision-making by management. The existing snow depth products have low spatial resolution and are mostly applicable to large-scale studies; however, they are insufficiently accurate for the estimation of snow depth on a regional scale, especially in shallow snow areas and mountainous regions. In this study, we coupled SSM/I, SSMIS, and AMSR2 passive microwave brightness temperature data and MODIS, TM, and Landsat 8 OLI fractional snow cover area (fSCA) data, based on Python, with 30 m spatially resolved fractional snow cover area (fSCA) data obtained by the spatio-temporal dynamic warping algorithm to invert the low-resolution passive microwave snow depths, and we developed a spatially downscaled snow depth inversion method suitable for the Turpan–Hami region. However, due to the long data-processing time and the insufficient arithmetical power of the hardware, this study had to set the spatial resolution of the result output to 250 m. As a result, a day-by-day 250 m spatial resolution snow depth dataset for 20 hydrological years (1 August 2000–31 July 2020) was generated, and the accuracy was evaluated using the measured snow depth data from the meteorological stations, with the results of r = 0.836 (p ≤ 0.01), MAE = 1.496 cm, and RMSE = 2.597 cm, which are relatively reliable and more applicable to the Turpan–Hami area. Based on the spatially downscaled snow depth data produced, this study found that the snow in the Turpan–Hami area is mainly distributed in the northern part of Turpan (Bogda Mountain), the northwestern part of Hami (Barkun Autonomous Prefecture), and the central part of the area (North Tianshan Mountain, Barkun Mountain, and Harlik Mountain). The average annual snow depth in the Turpan–Hami area is only 0.89 cm, and the average annual snow depth increases with elevation, in line with the obvious law of vertical progression. The annual mean snow depth in the Turpan–Hami area showed a “fluctuating decreasing” trend with a rate of 0.01 cm·a⁻¹ over the 20 hydrological years in the Turpan–Hami area. Overall, the spatially downscaled snow depth inversion algorithm developed in this study not only solves the problem of coarse spatial resolution of microwave brightness temperature data and the difficulty of obtaining accurate shallow snow depth but also solves the problem of estimating the shallow snow depth on a regional scale, which is of great significance for gaining a further understanding of the snow accumulation information in the Tuha region and for promoting the investigation and management of water resources in arid zones.

1. Introduction

Water is a very important natural resource in arid and semi-arid regions, and snow resources, as a frozen form of water resources, are an indispensable variable in the study of climate dynamics and the global hydrological cycle, with the high albedo and the adiabatic effect of the snow layer having far-reaching effects on hydrological and water resources, on the energy exchanges between the land surface and the atmosphere, and on the ecosystem as a whole [1,2,3]. As a major snowpack variable and a leading indicator for climate studies, hydrological applications, weather forecasting, and disaster assessment, snow depth reflects the dynamics of the snowpack and is an important parameter involved in climate and hydrological model simulations [4,5]. Snow depth data can be used not only to analyze climate change but also to estimate snow accumulation, predict snowmelt runoff, and assess flood risks [6,7,8,9]. Seasonal snowmelt is a stable and reliable freshwater recharge resource in arid and semi-arid regions, regulating river runoff and groundwater volume [10], ensuring habitats for wildlife species [11], and providing local residents with normal water for production and living. Located in the hinterland of the Asia–Europe continent and surrounded by mountains on all sides, the Turpan–Hami region is an extremely arid area, where ecosystems, spring climate, climate change, and the development of production and life benefits from the glaciers can be observed in the eastern foothills of the Tianshan Mountains and occur as a result of the seasonal snowmelt of the Tianshan water system that provides important freshwater resources and abundant groundwater [12,13,14,15]. This shows that snow resources play an important role in the Turpan–Hami region. Therefore, the accurate estimation of snow depth is essential for water resource management, understanding climate change, and maintaining ecosystem stability in the Turpan–Hami region.

Snow depths measured directly by automatic weather stations and on the ground (snow surveys) are usually considered to be the true values with the highest degree of accuracy. However, it is difficult to obtain snow cover and snow depth on a global or regional scale because of the sparse distribution of weather stations and the difficulties in obtaining manual measurements, especially in high mountains or sparsely populated areas [16,17]. In contrast, remote-sensing satellites can provide the required collection of large-scale snowpack information [18,19]. Passive microwave remote sensing, which can penetrate clouds and is unaffected by weather variations, is considered to be the most effective method for measuring snow depth on both global and regional scales because it not only provides stable snow information with day-by-day resolution but also calculates snow depth via receiving radiant energy from the snow cover and the surface and converting it into brightness temperature values [20,21]. There are three snow depth products covering the Turpan–Hami area: the AMSR2 standardized snow depth product (10 km) obtained from the website of the Japan Aerospace Exploration Agency (https://gportal.jaxa.jp/gpr/), the long time series snow depth data (25 km) obtained from the National Tibetan Plateau Data Centre (https://doi.org/10.11888/Geogra.tpdc.270194), and the daily values of snow depth in Xinjiang (500 m) from the National Glacial Permafrost Desert Science Data Centre (http://www.ncdc.ac.cn). The first two methods were developed on the basis of Chang’s algorithm [22], which is the earliest snow depth estimation method, that estimates snow depth through the linear fitting of the bright temperature difference between 18 GHz (or 19 GHz) and 36 GHz (or 37 GHz) to the measured snow depth. Due to the coarse spatial resolution of the bright temperature difference data and the difficulty in capturing the fine-scale features of snow depth, especially in mountainous areas with complex terrain, the data are insufficient to accurately estimate the snow depth in the Turpan–Hami area. The third is a random forest algorithm that introduces the vegetation index (NDVI), snow cover days (SCDs), topographic data, and spatial location data on the basis of long time series snow depth data products and snow depth measurements. However, while the algorithm increases the spatial resolution of the long time series snow depth data products and improves the estimation accuracy to a certain extent, machine learning algorithms demonstrate a high degree of dependence on their own structures, which can lead to overfitting the training data and poor interpretability [23].

Optical remote-sensing techniques cannot specifically quantify snow attributes, such as snow depth and snow water equivalent, but can provide fine-grained monitoring of the snow area coverage information. Walters et al. [24] downscaled MOD10A1, the fractional snow cover area (fSCA, 500 m) data, to a higher spatial resolution (30 m) snow-covered area (SCA) estimate by constructing a linear model. In order to reduce the uncertainty associated with the low spatial resolution of passive microwave remote sensing snow products, Dai et al. [25] established a new spatial dynamics method to improve the spatially resolved snow depth by combining passive microwave bright temperature data (AMSR2), the snow cover fraction (SCF), and surface emissivity (LST).

In order to create a high spatiotemporal resolution snow depth product for the Tibetan Plateau, Yan et al. [26] developed a spatiotemporal downscaling method that combined high spatial resolution 8-day cloud-free snow cover (SCP) data with passive microwave data (a long time series snow depth dataset) that were already available. A number of scholars have combined optical snow products (e.g., dichotomized snowpack, snow cover, and annual snow duration) with resampled passive microwave snow depth products to obtain snow depth and snow water equivalent information, which is an important step in the development of accurate snow depth products [27,28,29]. Berman et al. [30] developed a dynamic time warping (DTW) algorithm based on optical data (MODIS data and Landsat series data) to create 30 m spatial resolution daily fractional snow cover area (fSCA) data. Subsequently, Laurin et al. [11] extended their method to a global scale based on the Google Earth Engine and the R programming language’s SnowWarp tool, which opens up more possibilities for producing refined snow depth products in the Turpan–Hami region. Additionally, a few scholars have used microwave remote sensing to invert snow depth and then estimate regional snow accumulation and its spatiotemporal distribution characteristics. Yang et al. [31] conducted a long-term snowfall simulation in the Tianshan region using the WRF/Noah-MP model, based on ERA5 reanalysis data and real-time updated leaf area index and green vegetation coverage data, quantifying the snowfall in the Tianshan area. Other scholars have accurately estimated snowfall in mountainous areas by correctly reproducing topographic precipitation at high spatial resolutions based on regional climate models with land–atmosphere coupling [32,33,34]. Zhu et al. combined long-term microwave remote sensing snow depth data in China and used microwave remote sensing FY-3B/MWRI data to estimate snow density, thereby calculating the snow accumulation in the Tianshan region [35].

The Turpan–Hami region is located in the eastern part of Xinjiang, which is an important transportation hub and the core zone of the Silk Road Economic Belt in China, an area that is very rich in coal and petroleum resources, with good prospects for the development of specialty melon and fruit agriculture and tourism but with an extreme scarcity of water resources, sparse vegetation, and extremely fragile ecosystems [36]. Seasonal snowmelt is the most important freshwater resource in the region, but most of the existing studies on the snowpack in Xinjiang have focused on the snow stabilization zone in northern Xinjiang, with relatively little attention having been paid to the Turpan–Hami region. The existing snowpack products are insufficient to meet the urgent needs of snowmelt runoff simulations, hydrological water resource management, and the maintenance and improvements to ecosystem stability in the region.

In the current study, 250 m spatial resolution products for the daily fractional snow cover area (fSCA) data from 2000 to 2020 are produced through the DTW algorithm created by Berman et al. [30] using the SnowWarp tool of Laurin et al. [11]. Next, a linear fitting relationship between the measured snow depth data of the station and the passive microwave brightness temperature difference data are found using the Chang algorithm [22]. Finally, a spatially downscaled snow depth inversion algorithm applicable to the Turpan–Hami area was developed, and fractional snow cover area (fSCA) data were established. The algorithm proposed in this study not only more accurately estimates the snow depth in mountainous areas, to a certain extent, but also solves the problem of inaccurate snow depth inversion for shallow snow. In addition, the spatially downscaled snow depth inversion algorithm developed in this study generated a daily snow depth dataset with 250 m spatial resolution for the Turpan–Hami region for the period of 2000–2020, from which the temporal and spatial variations in snow depth in the Turpan–Hami region were further analyzed, aiming to provide data and technological support for the Turpan–Hami region’s development.

2. Study Area and Data

2.1. Overview of the Study Area

The Turpan–Hami region (hereafter referred to as the Tuha region) is located in the eastern part of the Tianshan mountain system in eastern Xinjiang, China, (Figure 1) and has a continental warm–temperate arid climate. It is surrounded by a low mountainous waterless desert running east–west and is separated from the Tarim Basin and Junggar Basin to the north and south, respectively. The region has long sunshine hours, strong radiation, scarce precipitation, vigorous evaporation, uneven spatial and temporal distribution of water resources impacting water scarcity, extreme scarcity of vegetation, severe land sanding, and extremely fragile ecosystems [36,37].

The glaciers and melting snow in the mountainous areas are the main water sources for the production and lives of the people in Turpan, Shanshan, and Hami, supplementing the rich groundwater resources; there are very few rivers with water flowing all year round, which means that the region belongs to the water resource scarcity area, and the contradiction between the supply and demand of water resources in the Tuha area is prominent. As one of the key areas of the national-level project, the Third Scientific Expedition to Xinjiang, Tuha region, presents certain challenges to water resources investigation and research due to its rugged terrain and relatively sparse meteorological stations [38]. Therefore, it is of practical significance to accurately understand the snow information in the Tuha region and assess the local snow water resources, water cycle processes, water resources management, and ecological environmental protection.

2.2. Data

Table 1. Information on data and datasets.

Data	The Purpose of the Data	Spatial Resolution	Time Resolution	Data Sources
MODIS	Used for the process of snow depth inversion.	500 m	1 day	https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MOD10A1 (Accessed on 7 November 2023).
Landsat		30 m	16 day	https://developers.google.com/earth-engine/datasets/catalog/landsat (Accessed on 7 November 2023).
SSM/I and SSMIS		25 km and 12.5 km	1 day	http://www.ncdc.ac.cn (Accessed on 7 November 2023).
AMSR2		10 km	1 day	https://gportal.jaxa.jp/gpr/ (Accessed on 7 November 2023).
SRTM DEM	Analysis used for downscaling snow depth products.	90 m		https://earthexplorer.usgs.gov/ (Accessed on 7 November 2023).
RFSD product	Evaluation of accuracy used for downscaling snow depth products.	500 m	1 day	http://www.ncdc.ac.cn (Accessed on 7 November 2023).
PMSD product		25 km	1 day	https://doi.org/10.11888/Geogra.tpdc.270194 (Accessed on 7 November 2023).

summarizes the spatial and temporal resolution and availability for all data sources used in this study. Detailed information for each source is provided in the following subsections.

2.2.1. Optical Data

The optical data used in this study are the NDSI snowpack data from the MODIS10A1 V6 Snow Cover Daily Global 500 m product (daily, 500 m) and the NDSI snowpack data from the Landsat series of remotely sensed imagery (16 days, 30 m). These products are used for the dynamic time warping algorithm, which was completed using the SnowWarp tool [11]. The time range was 1 August 2000–31 July 2020, and in this study, 1 August of the current year–31 July of the following year is defined as a hydrological year, in which, due to the limitations of the algorithm, the Rui Nian day (29 February) is excluded, for a total of 20 hydrological years.

2.2.2. Passive Microwave Brightness Temperature Data

This paper used SSM/I, SSMIS, and AMSR2 passive microwave brightness temperature data from the National Glacial Tundra Desert Science Data Centre and the GCOM-W1 data service network. In Xinjiang, the transit time of the AMSR-2 ascending orbit is day time, and the transit time of descending orbit is night time, taking into account that the snow surface temperature rises due to sunlight irradiation and temperature increases, which cause the snow to melt faster during the transit of ascending orbit than that during the transit of descending orbit, and it is easy to affect the inversion accuracy of the downscaling algorithm [27,39,40]. It has been shown that the horizontal polarization approach is more sensitive to snow identification [39,41]. Therefore, in this paper, this study mainly used brightness temperature data at different frequencies for the descending orbit, in horizontal polarization.

• DEM Data

The SRTM DEM data were obtained from the USGS with a spatial resolution of 90 m and a data format of GeoTIFF. Preliminary preprocessing operations were performed on the SRTM DEM data, which mainly consisted of splicing, cropping, and resampling, in order to obtain DEM data with a resolution of 250 m within the Tuha area for analyzing the inter-annual variability of elevation on snow depth.

2.

Xinjiang Area 500 m Daily Snow Depth Dataset (2010–2020)

The Xinjiang 500 m daily snow depth dataset (2010–2020), hereafter referred to as the RFSD product, is based on the Chinese long time series snow depth dataset, longitude, latitude, SRTM elevation, slope, slope direction, surface roughness, MOD13A1 Normalized Vegetation Index (NDVI) data, and the number of days of snow cover. This dataset, established using the random forest algorithm for snow depth downscaling, was obtained from the National Glacial Tundra and Desert Science Data Centre. The snow depth dataset before and after the downscaling was evaluated and compared in terms of accuracy using the measured snow depth data from meteorological stations in Xinjiang; then, a day-by-day 500 m downscaled snow depth product was produced in Xinjiang. The dataset was stored in TIF file format, and this study used this dataset to analyze and compare the accuracy of the products.

3.

Chinese snow depth long time series dataset

The passive microwave snow depth dataset (hereafter referred to as the PMSD product) used in this study comes from the National Tibetan Plateau Data Centre [42]. The dataset has a day-by-day temporal resolution and a spatial resolution of 25 km, and its raw data are daily passive microwave brightness temperature data from SMMR (1979–1987), SSM/I (1987–2007), and SSMI/S (2008–2018) archived at the National Snow and Ice Data Centre (NSIDC).

4.

Ground Meteorological Station Data

The snow depth data from ground-based stations are sourced from the National Meteorological Information Center, China’s national-level ground weather station dataset of daily values of basic meteorological elements; these data include station number, observation date, longitude, latitude, altitude, daily snow depth, and daily snow pressure. The experiment collected daily snow depth data from ten meteorological stations in the Tuha region for the hydrological years spanning from 2000 to 2020. Among these, the daily snow depth data from five hydrological years (2010–2015) at the ten meteorological stations were independently used as a validation dataset for the downscaling algorithm. Furthermore, empty and invalid data were eliminated during the pre-processing of algorithm simulation and product accuracy assessment.

Additionally, all data sources utilized in this study, along with their spatial and temporal resolutions, are summarized in Table 1.

3. Research Methodology

In this study, a spatially downscaled snow depth inversion algorithm (as shown in Figure 2) was used to combine multi-source remote sensing data to obtain snow depth data in the Tuha region. Firstly, this study obtained the day-by-day fSCA data according to the dynamic time warping algorithm; then, this study improved the appropriateness based on the Chang algorithm [22] coupled with the SSM/I, SSMIS, and AMSR2 passive microwave brightness temperature data. Finally, this study developed a set of spatially downscaled snow depth inversion algorithms for the Tuha region based on the Python construction of the linear fitting relationship between the brightness temperature difference, fSCA data, and the site-measured snow depth data and the establishment of the standardized product production system. The specific methods are described below.

3.1. Preparation of the Fractional Snow Cover Area (fSCA) Data

In this study, the SnowWarp tool (https://github.com/bermane/snowwarp/) proposed by Laurin et al. [11], based on Google Earth Engine (GEE), was used to obtain the fSCAs data for the years 2000–2020 in the study area. The algorithmic core of the tool is the DTW algorithm, which was originally developed for spoken language recognition [43] and was later developed as a data fusion algorithm for temporal normalization between two datasets [44].

The DTW algorithm used in this study was developed by Berman et al. [30], who used the DTW algorithm to match MODIS NDSI daily data and rearrange the Landsat fSCAs for historical periods to create a denser time series to obtain FSC estimates with a spatial resolution of 250 m. The original intention of this study was to use the DTW algorithm to combine MODIS fSCAs (500 m, 1 day revisit time) and Landsat fSCAs (30 m, 16 day) to obtain a time series of the daily fSCA data with a resolution of 30 m. However, due to the long data processing time and insufficient hardware to support the computational needs, the resultant output setting in the current study was changed to 250 m to obtain the fSCA data at 250 m spatial resolution. In addition, considering the relatively large scope of the study area, the large quantity of data, and the limitations of various aspects such as running memory in this study, the study region was broken into smaller subareas for data processing. Figure 3 illustrates the approach of block processing used to optimize data handling and reduce computational demands.

The specific methods were as follows:

• Upload the shapefile of the study area (consisting of a total of 124 small subregions) to the GEE platform.
• Preprocessing of MODIS and Landsat data based on GEE. Since GEE can directly acquire daily cloud-free MOD10A1 data, this study only needed to perform cloud masking on Landsat data, where the CFMASK algorithm was used and the cloudiness threshold was set to <70%.
• Obtain the fSCA data were using the SnowWarp tool [11]. The fSCA data were obtained based on the regression formula (Equation (1)) developed by Salomonson et al. [45] based on the correlation between NDSI and fSCA and adjusting the fSCA_adj data, which were obtained based on the method of Berman et al. [30] (Equation (2)) for forest canopy cover:

fSCA = 1.45 × NDSI − 0.01

(1)

fSCA_adj = fSCA/(1 − fCC)

(2)

where NDSI is the normalized snow deposit index; fSCA and fSCA_adj are the fractional snow cover area and the adjusted snow cover area, respectively; and fCC is the proportion of canopy cover data, which indicates the proportion of the pixel area covered by trees.

d.

Extract and process SnowWarp annual key statistics to obtain the day-by-day fSCA products at a spatial resolution of 250 m over each of the 124 small subregions.

e.

Mosaicking using the Seamless Mosaic tool of the Environment for Visualizing Images 5.3(ENVI) software to obtain fSCA data at the scale of the entire study area to use in the spatially downscaled snow depth inversion process.

3.2. Spatial Downscaled Snow Depth Inversion

According to previous studies and simulation analyses of the relationship between the brightness temperature of microwave radiation and snow depth, the brightness temperature difference between 18 GHz (or 19 GHz) and 36 GHz (or 37 GHz) is considered to be linearly related to the snow depth when the snow depth is small (<50 cm) [22,41,46,47], as shown in the Chang algorithm depicted in Equation (3). In particular, establishing a linear inversion algorithm from a statistical point of view is currently a feasible approach in the absence of more information about snow cover parameters [39].

$$\mathrm{S}\mathrm{D}=1.59\times ({\mathrm{T}\mathrm{B}}_{18\mathrm{H}}-{\mathrm{T}\mathrm{B}}_{37\mathrm{H}})$$

(3)

where SD is the snow depth; TB_18Hand TB_37H are the 18 GHz and 37 GHz horizontally polarized bright temperatures, respectively.

Herein, snow depth data from 10 ground-based snow depth measurement stations in the Tuha region for the years 2000–2011 and 2016–2020 are utilized as ground truth (totaling 14,247 valid snow depth measurements). Considering the low spatial resolution of the passive microwave data and transient changes in the snow accumulation, this study incorporated the daily spatial resolution of 250 m obtained by the above algorithms into the fSCA data with a daily spatial resolution of 250 m for spatial downscaling. This led to the development of an empirical model-based spatially downscaled snow depth inversion algorithm applicable to the Tuha region, and thus this study obtained the daily snow depth dataset of the Tuha region for the period of 2000–2020. In addition, before performing the spatially downscaled inversion, the Grody decision tree method was utilized to remove the above scatterers considering that the surface features such as rainfall, wind and sand, tundra, and cold desert would produce scattering features that were similar to the snow layer [48].

Before downscaling, this study considered two cases. The first case is if the fSCA image element value is 0, where regardless of whether there is snow on the passive microwave image element, the value of the image element is judged to be 0; the second case is if the fSCA image element value is between 0 and 100%, and then the downscaling inversion formula is utilized to calculate the snow depth value:

SD = (0.25 × (TB_18h − TB_36h) − 12.28) × fSCA

(4)

SD is the snow depth value of the snow-covered area on each image, TB_18h and TB_36h represent the brightness temperature value of each band on the corresponding image, and fSCA is the adjusted fSCA obtained from each pixel using the SnowWarp tool.

3.3. Snow Accumulation Estimation

This study utilized the snowpack estimation method (as shown in Equation (5)) proposed by Zhu et al. [35] to calculate the multi-year average snow mass for each raster of the snow period, accumulation period, stable period, and melt period in the Tuha region. The snow period is defined as November to March, with November–December as the accumulation period, January–February as the stabilization period, and March as the ablation period [49]. Snow density was used as 0.18 g/cm³ as proposed by Che et al. [50]:

M = A × SD × ρ

(5)

where M is the multi-year average snow mass per raster (in tons); A is the raster area per image (in km²); SD is the multi-year average snow depth within each period (in cm); and ρ is the snow density (in g/cm³).

3.4. Single-Factor Linear Regression Analysis Method

A linear trend analysis method was used to analyze the interannual trends in snow depth and snow accumulation at the spatial precipitation scale in the Tuha region. The one-way linear regression model was formed using y_i to denote the snow depth in year i and x to denote the year number:

y_i = a + bx_i

(6)

where a is the intercept in the unidirectional linear regression model, which represents the expected snow depth within the hydrological year and b is the trend in snow depth over time; if b > 0, snow depth increases; b < 0, snow depth decreases; and b = 0, snow depth does not change.

3.5. Sen + Mann–Kendall Trend Analysis

In this study, Sen trend analysis and the Mann–Kendall test were used to calculate the time series trend in snow depth in the Tuha region for the hydrological years of 2000–2020 to obtain the interannual rate of change in snow depth and to test its significance. Sen, as a nonparametric statistical trend analysis method, does not require the samples to obey a specific distribution and also shows good robustness in dealing with missing values and outliers in the series (the interference of abnormal year values can be excluded) [51].

The Sen trend analysis equation is

β = median((X_j − X_i)/(j − i)),∀_j > i

(7)

where β is the trend degree of snow depth change, the median () indicates the median function; i, j is the number of time series; X_i, X_j is the snow depth value of time series i, j. When β > 0, it means that the snow depth has an upward trend; when β < 0, it means that the snow depth has a downward trend.

The Mann–Kendall test is a nonparametric test for determining the significance of the trend in snow depth [52]. It is calculated as

$$\mathrm{Z}=\left\{\begin{array}{c}\frac{\mathrm{S}-1}{\sqrt{\mathrm{V}\mathrm{a}\mathrm{r}\left(\mathrm{S}\right)}},S>0\\ 0,S=0\\ \frac{\mathrm{S}+1}{\sqrt{\mathrm{V}\mathrm{a}\mathrm{r}\left(\mathrm{S}\right)}},S<0\end{array}\right.$$

(8)

$$S=\sum _{\mathrm{i}=1}^{\mathrm{n}-1}\sum _{\mathrm{j}=\mathrm{i}+1}^{\mathrm{n}}\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left({\mathrm{X}}_{\mathrm{j}}-{\mathrm{X}}_{\mathrm{i}}\right)$$

(9)

$$\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left({\mathrm{X}}_{\mathrm{j}}-{\mathrm{X}}_{\mathrm{i}}\right)=\left\{\begin{array}{c}+1,{\mathrm{X}}_{\mathrm{j}}-{\mathrm{X}}_{\mathrm{i}}>0\\ 0,{\mathrm{X}}_{\mathrm{j}}-{\mathrm{X}}_{\mathrm{i}}=0\\ -1,{\mathrm{X}}_{\mathrm{j}}-{\mathrm{X}}_{\mathrm{i}}<0\end{array}\right.$$

(10)

$$\mathrm{V}\mathrm{a}\mathrm{r}\left(\mathrm{S}\right)=\frac{\mathrm{n}\left(\mathrm{n}-1\right)\left(2\mathrm{n}+5\right)}{18}$$

(11)

where Z is the standardized test statistic; S is the test statistic; and n is the length of the dataset. In this study, the significance of the interannual trend in snow depth change in the Tuha region was tested at the confidence level α = 0.05. When the test result |Z| > 1.96, the change trend was considered significant; when |Z| ≤ 1.96, the change trend was considered non-significant.

4. Results

4.1. Comparative Analysis of Overall Accuracy of Snow Depth Products

The accuracy for each of the existing PMSD products, RFSD products, and the spatially downscaled snow depth produced in this study (SD) was verified using the daily snow depth values from ground stations at 10 meteorological stations over 5 hydrological years from 2010 to 2015. The overall accuracy of the three products was evaluated using five evaluation indicators: coefficient of determination (R²), bias, Pearson’s correlation coefficient (R), root-mean-square error (RMSE), and mean absolute error (MAE), as shown in

Table 2. Overall accuracy of three snow depth products.

Data	R2	Bias	R	RMSE	MAE
SD	0.69	−0.36	0.84 **	2.60	1.50
PMSD	0.51	−0.21	0.72 **	3.29	2.45
RFSD	0.61	−0.83	0.83 **	2.95	1.74

**: p < 0.01, significant correlation.

.

According to Table 2, all products passed the significance test (p ≤ 0.01). The spatially downscaled snow depth products produced in this study exhibited the highest R² and the lowest RMSE and MAE (R² = 0.69, R = 0.84, RMSE = 2.60, MAE = 1.50, Bias = −0.36). Compared to the PMSD products (R² = 0.51, R = 0.72, RMSE = 3.29, MAE = 2.45, Bias = −0.21), the RFSD products (R² = 0.61, R = 0.83, RMSE = 2.95, MAE = 1.74, Bias = −0.83) exhibited better overall accuracy. Therefore, the above analysis indicates that the accuracy ranking of the three snow depth products is as follows: spatially downscaled snow depth products produced in this study > RFSD products > PMSD products.

Furthermore, to further validate the spatially downscaled snow depth products produced in this study, accuracy evaluations were conducted according to different hydrological years. Scatter plots of computed versus measured snow depths for all five analysis years (2010–2015) are shown in Figure 4a. Additionally, scatter plots for each hydrologic year are shown in Figure 4b–f. Overall, the accuracy and performance results within the five hydrological years are satisfactory, and all have passed the significance test.

To further validate the reliability of the spatially downscaled snow depth retrieval algorithm developed in this study, this study conducted accuracy verifications of the downscaled snow depth product produced in this study based on the PMSD product (as shown in Figure 5a) and RFSD product (as shown in Figure 5b), respectively. The accuracy between the downscaled snow depth product produced in this study and the RFSD product is relatively high (R² = 0.65, R = 0.82, Bias = −0.47, RMSE = 2.34, MAE = 1.48). The PMSD product is one of the most widely used snow depth products in the field of snow cover. Therefore, this study also conducted accuracy evaluations between the PMSD product and the RFSD product and found that the accuracy is comparable to that between the PMSD product and the downscaled snow depth product produced in this study. Overall, the verification results among the downscaled snow depth product produced in this study, the RFSD product, and the PMSD product are relatively close, and the test results of the downscaled snow depth product produced in this study, the PMSD product, and the RFSD product are all significant at the 0.01 significance level. Hence, the downscaled snow depth product produced in this study is relatively reliable.

Figure 6 shows the accuracy of the three snow depth products compared to the field snow depth true values over a range of snow depth intervals to evaluate the performance of each the products over the different intervals. In this study, the field snow depth true values were categorized into six intervals (0 ≤ SD < 3 cm; 3 ≤ SD < 6 cm; 6 ≤ SD < 9 cm; 9 ≤ SD < 12 cm; 12 ≤ SD < 15 cm; and SD ≥ 15 cm). For different snow depth zones, the three snow depth products performed best in the shallow snow zone (0 ≤ SD < 3 cm), with RMSEs of 1.24 cm, 2.23 cm, and 1.84 cm, respectively, but the spatially downscaled snow depth product had the smallest RMSE (RMSE = 1.23) and the best performance in the shallow snow zone (0 ≤ SD < 3 cm), when compared to the PMSD product and the RFSD product. The PMSD product had the largest RMSE when the snow depth was in the interval range (0 ≤ SD < 3 cm, 3 ≤ SD < 6 cm, 6 ≤ SD < 9 cm, 9 ≤ SD < 12 cm). This was due to the coarser spatial resolution of the PMSD product, which is approximately 25 km and is more suitable for exploring large-scale snow depth studies, and it is unable to accurately estimate snow depths in small areas of shallow snow. Due to the limitations of the algorithm itself, the RFSD product overestimates the snow depth to a certain extent, resulting in the largest RMSE of the RFSD product in the interval range (12 ≤ SD < 15 cm and SD ≥ 15 cm). The SD product is based on the dynamic time warping algorithm, which can obtain more accurate snow distribution information; thus, the spatially downscaled snow depth inversion algorithm developed for the Tuha region reflects the snow depth information in the Tuha region more realistically through the statistical relationship between the bright temperature difference, the true value of the snow depth, and the fSCA data. As a result, the spatially downscaled snow depth product always has the smallest root-mean-square error and the highest accuracy within the six snow depth intervals mentioned above.

The spatially downscaled snow depth product assigns the snow depth to each pixel linearly based on the true value of the ground snow depth and the fSCA data, which not only simulates the true value of the ground snow depth to the maximum extent but also accurately grasps the proportion of the snow distribution in each pixel, so as to produce the high spatial resolution snow depth product, and therefore, the performance of the spatially downscaled snow depth product developed in this study is relatively good.

Although the spatially downscaled snow depth improves the accuracy of the snow depth algorithm in the Tuha area, with the spatial resolution increased to 250 m, the spatial heterogeneity of the snowpack is strong, and the representative error of the inversion process based on the snow depths of the existing meteorological stations is unavoidable. When the snow depth was ≥15 cm (the high value of snow depth in the Tuha area), the analysis demonstrated that the accuracy became relatively worse. There are multiple reasons for the poor snow depth accuracy, including elevation difference, slope direction difference, climate change, subsurface difference, wind speed and direction, and the physical properties of snow. In addition, snow depths of 15 cm and above in the Tuha region are basically distributed in mountainous areas with higher elevations, while the snow depth measurements based on the spatially downscaled snow depth algorithm developed in this study were based on station data. However, the stations are mostly distributed on plains with lower elevations, and thus, the stations are unable to provide snow depth measurements in individual higher elevation areas, which may be the main reason for the lower accuracy of the spatially downscaled snow depth products in high snow depth areas. According to previous studies, ground snow depth data are indispensable for further snow depth detection. In order to improve the accuracy of snow depth estimation in the Tuha region and to ensure the reliability of the validation results, it is essential to continuously and effectively conduct a large number of snow surveys in future work. Therefore, in subsequent research, this study can address this limitation by establishing additional stations on different surfaces and at different elevations or by conducting more extensive snow surveys.

4.2. Spatial Comparative Analysis of Different Snow Depth Products

In order to further reflect the spatial distribution of different snow depth products in Tuha, this study, based on Figure 7, examined the spatial distribution of the snow depth on 7 February 2011 as an example to show the local details of the snow depth distribution. The PMSD product estimated the maximum snow depth of that day to be 18 cm, which is more in line with the meteorological stations. In fact, the snowpack in Tuha is mostly distributed in mountainous areas with a complicated topography and subsurface conditions. Figure 7b displays the PMSD products, with a spatial resolution of 25 km. However, this resolution greatly restricts the detailed information of the spatial distribution of snow depth and only roughly reflects the spatial distribution of snow depth in the Tuha region; therefore, it is difficult to define the boundaries of the snowpack distribution. The RFSD product was developed based on the PMSD product by combining the topography, vegetation index, location factor, and the number of snow days, forming a 500 m spatial resolution product of the daily snow depth in Xinjiang. This product shows the spatial distribution of snow depth in detail, and the deep snow and shallow snow areas can be clearly distinguished in Figure 7c. From the whole region of Tuha, the snow depth values were generally larger, and the maximum snow depth reached approximately 43 cm on the same day, which may be overestimated compared with the true value of snow depth at the site. Figure 7a originates from the SD product, providing a more detailed presentation of the spatial distribution of the Tuha region. The figure shows that the deep snow areas are mainly distributed in the Bogda Mountains, western Barkun, and central Hami (North Tian Shan, Barkun Shan, and Harik Shan); the shallow snow areas are mainly distributed in other areas except the mountainous regions. The snow depth distribution of the spatially downscaled snow depth product is more consistent with that of the PMSD product, which shows that the spatially downscaled snow depth inversion algorithm developed by introducing the fSCA data obtained via the dynamic time warping algorithm not only provides a good basis for estimating the snow depth of shallow snow but also provides more detailed information about the snow distribution in the mountainous areas.

4.3. Trends in Snow Depth at Different Altitudes

Elevation is an important topographic factor that affects the distribution and change in snow depth; thus, it is to be expected that the Bogda Mountain, North Tianshan Mountain, Barkun Mountain, and Harlik Mountain in the Tuha region have deep snow cover which lasts for a large number of days, and the Turpan–Hami Basin has a lower elevation, high temperatures, and less precipitation, resulting in shallow snow cover. Finally, the snow depth in the mountainous areas of the Tuha region is generally higher than that in other regions.

Figure 8 illustrates the interannual variation in snow depth in different elevation zones in Tuha in the hydrological years 2000–2020. This study divided the elevation of the Tuha region into five elevation zones: less than 500 m, 500–1500 m, 1500–2500 m, 2500–3500 m, and more than 3500 m, with each zone separated by 1000 m increments in elevation. The results show that the average annual snow depth increases with elevation, with an obvious vertical incremental pattern. Where the elevation of the Turpan–Hami Basin is below 500 m, the mean annual snow depth is extremely low—below 0.5 cm—due to the high temperatures, which are unfavourable for the formation of snow. Where the elevation is below 2500 m, the mean annual snow depth is not more than 4 cm, and most of the snow is concentrated in the range of 2500 m. Finally, where the elevation is more than 3500 m, the mean annual snow depth is a maximum of about 8 cm. The increasing trend in snow depth (S = 0.02, p > 0.05) is due to the sudden increase in snow depth in this elevation range during the 2012–2017 hydrologic years. Overall, the inter-annual snow depths at different elevation ranges in the Tuha region have shown “fluctuating changes” since 2000.

4.4. Spatial and Temporal Variability of Snow Depth

4.4.1. Spatial Distribution of Snow Depths and Mass

The average annual snow depth in Tuha is only 0.89 cm, and the maximum annual snow depth does not exceed 17 cm. From the multi-year average snow depth spatial distribution depicted in Figure 9, it can be seen that the snow in most areas of Tuha exists in the form of shallow snow; the deep snow in Tuha is mainly concentrated in the mountainous areas, with obvious spatial heterogeneity in the distribution of snow, and there are two main areas of high-value snow depths, which are located in the north part of the Turfan area and the northwest and central parts of the Hami area, respectively. Compared with topographic maps, these two high-value snow depth zones correspond to the perennial glacier-covered areas of the Bogda, Northern Tianshan, Barkun, and Halik mountain ranges in the region.

Figure 10 illustrates the spatial distribution of snow depth in the four snow accumulation periods in the Tuha area: for many years, the average snow depth in the snow accumulation period was below 1 cm in the area accounting for 88.79% of the Tuha area, with the largest distribution range, 1–3 cm and 3–5 cm each, accounting for 6.53% and 2.35%, distributed in most of the areas in the western part of the Barkol, the low mountainous areas, and the surrounding areas, and a more than 5 cm average snow depth for many years. The average snow depth of more than 5 cm accounted for 2.32% of the Tuha area and was mainly distributed in the mountainous areas, where the maximum snow depth of many years was 14.13 cm, and the average snow depth of many years was approximately 0.83 cm. The spatial pattern of the multi-year mean snow depth during the snow period in the Tuha region is consistent with that of the stabilization period, but the value is slightly smaller, being approximately 0.71 cm. The multi-year mean snow depth during the stabilization period was 1.19 cm, which is more consistent with the study of Zhu Shuzhen et al. [35] focusing on the snow depth in the Tianshan region. The multi-year average snow depth in the ablation period was 0.67 cm, and the value of snow depth became smaller in the whole region. The distribution of the snow depth in the western part of Barkun became shallower, but the snow depth in the mountainous areas changed less, and there was almost no snow distribution in the Turpan–Hami Basin. Topography and climatic conditions are important factors affecting snow formation and snow depth distribution; usually, the region will have a high altitude, low temperatures, and a large amount of snow, while the slope direction will also profoundly affect the snow pattern and snow depth distribution on mountain slopes.

The snow mass in this study was computed using Equation (5), but since the contours are the same as in Figure 10 (with a different scale), they are not displayed here. The spatial distribution of the multi-year average snow volume during the snow accumulation period in the Tuha region is more consistent with the spatial distribution pattern of the snow depth, generally presenting the characteristics of abundant snow volumes in the mountainous areas and scarce snow volumes in the basin area. During the snow accumulation period, the multi-year average snow accumulation of each grid in the Tuha region was 55.04 t, mainly concentrated in the Bogda Mountain in the north of Turpan, the northern Tianshan Mountain in the central part of Hami, the Barkun Mountain, and the Harlik Mountain. The western part of the Barkun was also relatively abundant in snow accumulation, which is related to the lower temperatures of Barkun in the whole year. The multi-year average snow mass of each raster in the Tuha region plummeted from 78.64 t during the stabilization period to 24.80 t during the ablation period, a reduction of approximately two-thirds. The spatial distribution of the snowpack in the accumulation, stabilization, and ablation periods showed a tendency to expand and then shrink, centred on the mountainous areas and western Barkun, which was consistent with the trends in snowfall and temperature changes.

4.4.2. Inter-Annual Changes in Snow Depth

Figure 11 depicts the interannual changes in snow depth during the hydrological years 2000–2020. The average annual snow depth in the Tuha region fluctuates between 0.58 and 1.58 cm, with a multi-year average of 0.89 cm, and the average annual snow depth in the Tuha region during the hydrological years 2000–2020 shows a “fluctuating decrease” at a rate of 0.01 cm·a⁻¹. Annual average snow depths reached 1 cm or more during the 2002–2003, 2005–2006, 2009–2010, 2012–2013, and 2015–2016 hydrological years, with a peak of 1.58 cm during the 2005–2006 hydrological year; the three hydrological years from 2006 to 2009 showed a significant decrease in the snow depth “trough pattern” and an annual average snow depth of 0.65 cm, and in the hydrological year of 2008–2009 there was the lowest value of snow depth across these 20 hydrological years with a snow depth of 0.58 cm. Subsequently, the annual average snow depth rebounded within the range of the multi-year average snow depth and exhibited a trend of reduced fluctuations.

In this study, the results of the analysis of snow depth in the Tuha area during four periods (snow accumulation, accumulation, stabilization, and ablation) in 20 hydrological years were as follows, as shown in Figure 12: Since 2000, the annual mean snow depth in the Tuha area during the snow accumulation period has shown a slightly increasing trend (the trend is not significant); the changes in snow depth during the snow accumulation period and the ablation period have also shown an increasing trend, especially during the ablation period, with a large inter-annual variation, which is probably related to climate change [15]. During the stabilization period, the annual mean snow depth in the Tuha region showed a decreasing trend, which was consistent with the trend in the annual mean snow depth change in the Tuha region. The order of the size of the annual mean snow depth in the three periods of snow accumulation was stable period > ablation period > snow accumulation period. The correlation coefficient between the change in annual mean snow depth in the stabilization period and the trend in the annual mean snow depth change in the whole year is R = 0.976 (this passed the significance test of 0.01), which indicates that the snowfall in the stabilization period in the Tuha region is crucial for the Tuha region.

4.4.3. Temporal and Spatial Trends in Snow Depths

Figure 13 presents the spatial and temporal trends in the Tuha region obtained through Sen + Mann–Kendall trend analysis. From Figure 13a, it can be observed that over the past 20 hydrological years, there has been a spatiotemporal trend in snow depth variation in the Tuha region. Since 2000, the snow depths in the Northern Tian Shan, the Barkun Mountains, and the Harlik Mountains have shown an increasing trend, whereas their surroundings, the south of the Bogda Mountains in the northern part of the Turfan region and the western part of the Barkun region, have shown a decreasing trend in varying degrees.

Figure 13b overlays the trend analysis results with the Mann–Kendall test results at a confidence level of α = 0.05, categorizing the snow depth change trends into five categories: significantly increasing, non-significantly increasing, non- significant reduction, significant reduction, and stable and unchanged. This study statistically analyzed the proportions of different types in relation to the total area. The results show that the significantly increasing area accounts for 9.03%, those with a non-significant increase account for 24.54%, those with a non-significant reduction account for 36.49%, those with a significant reduction account for 27%, and those with no change in stability account for 2.93%. Overall, the snow depth in the Tuha region decreased more than it increased in the 20 hydrological years; thus, the overall snow depth change showed a decreasing trend.

5. Discussion

From the spatial distribution of the average snow depth obtained from the inversion, most of the snow in the Tuha region is mainly concentrated in the mountainous areas, with an average snow depth of approximately 0.89 cm and a maximum snow depth of no more than 17 cm, indicating that most of the snow in the Tuha region exists in the form of shallow snow. These results put forward higher requirements and challenges for the development of snow depth inversion algorithms with high temporal and spatial resolutions. The spatially downscaled snow depth inversion algorithm developed by introducing the fSCA data obtained via the DTW algorithm provides a good basis for estimating the snow depth of the shallow snow, and at the same time, it provides more detailed information on the snow distribution in mountainous areas. In the absence of additional snowpack parameters, spatial downscaling of passive microwave bright temperature data using the fSCA data obtained from the DTW algorithm is feasible and provides a more accurate and refined estimate than other products that include the depth of the snowpack in the Tuha region.

Despite the spatial downscaling of snow depth improving the accuracy of snow depth algorithms in the Tuha region, the spatial heterogeneity of the snowpack is strong. Therefore, representative errors in the inversion process based on snow depth data from existing meteorological stations are unavoidable. During the ablation period, our team found that a large amount of snow at a low altitude had melted and formed a deep frost layer. This was due to the increase in temperature, the formation of a deep frost layer in the snow layer under a certain temperature gradient, and the diffusion of water vapour, which led to changes in the snow grain size and snow density. In our future research, we not only plan to obtain physical properties such as snow temperature, snow grain size, snow density, and snow water content on the basis of snow surveys but also to use physical models such as MEMLS and HUT [53,54] to establish a snow depth model to more accurately simulate the radiative transfer process within the snow layer; additionally, this study will attempt to use multi-source sensors to estimate the snow depth changes in the Tuha area.

Further exploration of the above algorithms cannot be separated from the ground truth data of snow depth, and in order to improve the accuracy of snow depth estimation in the Tuha area and to ensure the reliability of the validation results, it will be necessary to carry out a large number of snow surveys regularly throughout the year in the future work.

6. Conclusions

This study developed a spatially downscaled snow depth inversion method for the Tuha region by coupling MODIS, TM, and Landsat8 OLI (area-likeness of snow) data with SSMI, SSMIS, and AMSR2 passive microwave brightness temperature data and produced a daily 250 m spatially resolved snow depth dataset for the Tuha region over 20 hydrologic years (2000–2020). The following conclusions were drawn from the analysis:

(1) The new spatial downscaling snow depth inversion algorithm, compared to other algorithms covering the Tuha region, not only improves the validation accuracy of the inverted snow depth data (R = 0.84 and significant at 0.01 level, R2 = 0.69, RMSE = 2.60, MAE = 1.50, Bias = −0.36) but also achieves more refined improvements at spatial scales. This, to a certain extent, enables accurate estimation of shallow snow depths in the Tuha region.

(2) The snow depth in the Tuha area has two high value areas in the spatial distribution, distributed in the north of Turpan (Bogda Mountain) and the northwestern part of Hami (Barkun Autonomous County) and in the central part of the area (North Tianshan Mountain, Barkun Mountain, and Harlik Mountain), with a maximum annual snow depth of 16.66 cm and an average snow depth of only 0.89 cm for many years. Since 2000, the snow depth in the North Tianshan Mountains, Barkun Mountains, and Harlik Mountains has demonstrated an increasing trend; meanwhile, its surroundings, south of the Bogda Mountains in the northern part of the Turfan area and the western part of the Barkun area, have shown a decreasing trend to varying degrees. Additionally, in the Tuha region, the annual average snow depth increases with elevation, showing a vertical incremental pattern. Meanwhile, the snow depth in the Tuha region from the inter-annual changes showed a “fluctuation decrease” trend, with fluctuations in the range of 0.58–1.58 cm. The average snow depths of the three periods in the snow accumulation period are in the order of stabilization period > ablation period > accumulation period, among which the snowfall amount in the stabilization period of the Tuha region is crucial for the Tuha region.

(3) The Tuha region snow accumulation period, stabilization period, and ablation period of the average annual snow amount of the mountainous areas show abundant snow; the basin area of snow demonstrates sparse characteristics; and the spatial distribution pattern of the snow depth is more consistent. Over a period of 2000 years, in the Tuha region, the average annual snow amount of each grid in the stabilization period was 78.64 t; then, it plummeted in the ablation period to 24.80 t, a reduction of about two-thirds.

(4) Additionally, in the future, when we have more computational resources available, there is a possibility of increasing the spatial resolution of the snow depth product to 30 m.

In general, the snow depth in the Tuha region decreased more than it increased in 20 hydrological years, and the snow depth showed a “fluctuating decrease” trend with a rate of 0.01 cm·a⁻¹. Comparatively speaking, glacial snowmelt is an important source of water resources in the Tuha region, and clarification of the distribution of snowpack, changes in snow depth, and the magnitude of the snowpack is of great significance to the understanding of hydrological processes, the simulation of snowmelt runoff, water resource management, and the analysis of resource decision-making in the region.