An Integrated Shadow-Adjusted Snow-Aging Index for Alpine Regions

Abstract

Detecting the variations in snow cover aging over undulating alpine regions is challenging owing to the complex snow-aging process and shadow effect from steep slopes. This study proposes a novel snow-cover status index, namely shadow-adjusted snow-aging index (SASAI), portraying the integrated aging process within the Manas River Basin in northwest China. The Environment Satellites HJ-1A/B optical images and in-field measurements were used during the snow ablation and accumulation periods. The in-field measurements provide a reference for building a candidate library of snow-aging indicators. The representative aging samples for training and validation were obtained using the proposed time-gap searching method combined with the target zones established based on the altitude of snowline. An analytic hierarchy process was used to determine the snow-aging index (SAI) using multiple optimal snow-aging indicators. After correction by the extreme value optimization algorithm, the SASAI was finally corrected for the effects of shading and assessed. This study provides both a flexible algorithm that indicates the characteristics of snow aging and speculation on the causes of the aging process. The separability of the SAI/SASAI and adaptability of this algorithm on multiperiod remote sensing images further demonstrates the applicability of the SASAI to all the alpine regions.

1. Introduction

The seasonal variability of snow cover influences the global hydrological balance and direct feedbacks in the climate system [1]. Snow cover recession is crucial in alpine areas where ablation affects the local water cycle and energy balance, water resource management, and disaster prevention and control [2,3]. During the ablation period, the snow cover area (SCA) shows a stepwise retreat until it melts completely or becomes permanent snow or a glacier while the snow particles change morphologically from fresh to grain snow and granular ice to glaciers [4]. Such seasonal morphological snow cover changing is summarized as the snow-aging process (SAP). Nowadays, the use of satellite imagery has become increasingly common and reliable for acquiring snow-cover information. Fine spatial and temporal resolutions offered by satellites allow real-time monitoring of SAP related parameters, which is reflected by the abundance of available algorithms and models from various perspectives in the literature.

At a macro level, SAP could be described as the change of geometric parameters such as snow-cover extent, snowline altitude, and snow depth [5,6]. Many snow cover extent products have been generated by combining classical algorithms such as normalized difference snow index (NDSI) [7], updated cloud removal techniques, and spatiotemporal data fusion techniques [8,9,10,11]. For example, moderate resolution imaging spectroradiometer (MODIS) data are broadly employed as a reliable data source for detecting the extent of snow cover owing to its high spatial, time, and spectral resolutions [12,13]. Its sensitivity to factors such as aerosol optical properties are also explored in alpine areas [14]. Further, the “subpixel snow-cover information”, “empirical relationship assumptions” [15,16], and “spectral unmixing” [17] models have been extensively applied for the inversion of the fractional snow cover. Transient snowline altitude and glacier elevation can be extracted by combining optical and synthetic aperture radar (SAR) imagery as well as DEM data [18]. Passive microwave-based models can penetrate clouds and provide measurements in shadowed regions and hence, are useful for inversion of the snow depth [19,20], snow water equivalent (SWE) [21], and snow cover storage [22]. Based on this, the composite snow cover products ESA GlobSnow SWE dataset [23] and snow data assimilation system (SNODAS) [24] were produced.

Microscopically, SAP mainly refers to physical metamorphism as a variation in the grain size and particle structure [1,2,3,4]. The accumulation of pollution in snow, increasing liquid water content, and surface roughness are also typical symptoms of SAP (these also influence snow metamorphism) [5,6]. Extensive research has been conducted on retrieving the snow grain size using “scaled band area” algorithm and the MODIS snow-covered area grain size (MODSCAG) model [15]. On the basis of the quasi-crystalline approximation (QCA) theory, the vertical structure of snow cover could be detected by the “interferometric synthetic aperture radar (InSAR) scatter model” [25]. The pollutant concentration was estimated by a contamination index based on snow reflectance [26,27,28]. The liquid water content could be obtained by water flux models incorporated within the 1-D snow cover model SNOWPACK [29]. By applying joint inversion using multispectral and hyperspectral data, the snowpack roughness can also be extracted [30]. Using image processing techniques, the extracted snow cover and glaciers or a fine classification of different snow cover types could also reflect the physical characteristics of the snowpack [31,32,33].

Indirectly, other external parameters also provide indicators for SAP, such as surface temperature of snow cover, snowpack duration, cloudiness, precipitation, solar zenith angle, and elevation [34,35,36,37]. Many researches have explored using correlation analysis, sensitivity and response models of snow melting and data pertaining to other parameters mentioned above [38,39]. The classical degree-day model establishes a linear correlation between accumulated temperature and snow melting [40].

Theoretically, from these parameters, the continuous monitoring of SAP can be achieved. However, extracting single-dimensional and single-index SAP parameters is not enough to understand the snow ablation process, which is complex and rapidly changing because of multiple factors. As this process requires multiperiod monitoring and normalization, it is time-consuming and inefficient. Besides, the relationship between the external SAP parameters and snow aging is not robust that results in highly biased snow-aging estimates. Multiple factors such as accumulation of pollution, liquid water content, temperature, solar zenith angle, cloud cover, precipitation, and snow cover duration influence the rate of snow aging [41,42]. Complex topography and heterogeneous land cover redistribute the snow cover over the alpine terrain, which adds difficulty in understanding the conditions of snow aging and its ablation process [43]. Because of these factors, the snow-aging variability and its sensitivity to multiple recession factors over the alpine areas has not been thoroughly investigated.

Snow heterogeneity affects microwave radiation and scattering signals. Surface features can produce scattering signals that increase uncertainties [44]. The bias in microwave radiation and scattering signals is caused by snow heterogeneity [45], scale inconsistency, and mixed pixel interference needs consideration by combining the optical and microwave data [46,47]. For these reasons, we prefer the use of optical data in the present work. The optical satellite data, acquired reflectance from the visible (VIS) to shortwave infrared (SWIR) spectrum with a high temporal resolution, were selected by constructing this snow-aging index (SAI). Therefore, the objectives are to: (i) obtain and analyze different snow-aging factors for the typical alpine area during the snowmelt period; (ii) develop an optical-based algorithm to integrally characterize these factors and to quantitatively reflect their contributions; (iii) optimize the effect of shadows on the integrated SAI especially in rugged terrains; (iv) evaluate the accuracy of the final result to characterize the spatiotemporal patterns of SAP in the alpine areas. This study will provide useful information for application in snow hydrological studies, regional climate change, and future snow projection, especially, during the ablation season.

2. Study Area and Data

2.1. Study Area

The Manas River Basin (MRB) is located between the northern foothills of the Tian Shan Mountains and south of the Junggar basin in China (43°05′ N–44°10′ N and 85°00′ E–86°20′ E) (Figure 1). Subjected to a continental drought climate, the MRB has a large altitudinal gradient. From north to south, the altitude increases from 850 to 5242 m and areas above 3900 m are covered with permanent snow. With the abundant snow and significant vertical gradient of temperature during ablation, MRB makes a suitable natural laboratory for this research, focusing on the spring–early summer period when the snowmelt and runoff processes are most active.

2.2. Satellite and Auxiliary Data

With two optical satellites in the same orbit plane, the HJ-1A/B have four-days revisit capability providing 30 m spatial resolution and four bands of CCD/IRS observations from visible to near infrared. The HJ-1A/B data are freely obtained from the China Centre for Resources Satellite Data and Application (http://www.cresda.com/CN/). The image preprocessing includes sensor calibration, orthorectification, atmospheric correction by the FLAASH model, and approved cosine terrain correction [48].

The auxiliary data include the field work and elevation data. The former is collected within the MRB for different periods to provide the basis for channel combinations that might be useful for detecting the snow-aging status. Twenty-eight sampling sites along the western and eastern branches were set up (see Figure 1). Grain size, water extent, density, reflectance, and extent of contamination were measured and recorded. The corresponding meteorological and satellite parameters such as snow surface temperature, air temperature, sensor zenith angle, and relative humidity were also collected. The Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM) Version 2 was employed for extracting the snowline and evaluating the shading correction, and for optimization.

3. Methods

Figure 2 shows the main datasets and methods applied for the shadow-adjusted SAI (SASAI). The in-field work collected the observational SAP parameters and spectra of different snow samples across the MRB. The obtained snow-aging spectral features were used as reference for constructing the candidate snow-aging indexed library. The “time-gap searching method” and target zones obtained from DEM were combined to capture the representative snow-aging training and validation samples from HJ-1A/B CCD/IRS reflectance data. A criterion was constructed to screen out the optimal snow-aging indexes with confidence. By quantifying the contributions of the indexes through the analytic hierarchy process (AHP), snow-aging index (SAI) was constructed. Within the classified areas (shaded/hot spot/flat), a terrain condition adjusted factor determined by an “extreme value optimization” algorithm was applied to adjust the SAI to eliminate its shaded and hot spot effect. Finally, the shading correction and separability of the SAI/SASAI were assessed by the local solar incidence cosine (LSIC) and validation samples.

3.1. Field Work

Two field work campaigns were carried out in the northern and southern branches of the MRB during the snow accumulation (10–17 December 2013) and ablation (15–22 March 2014) periods. The field work included the acquisition of reflectance spectra and collection of morphological and physical parameters, and other environmental variables [49,50]. The purpose of the field work was to acquire measurements that can help in choosing the channel combinations or bands that are sensitive for the detection of the snow-aging status.

Figure 3 displays the photographs of the sampling sites, and tight shots of fresh snow, clean aging snow, and contaminated aging snow. The corresponding reflectance spectra were measured with the ASD FieldSpec4 instrument, a portable and battery powered spectroradiometer that operates in three detector ranges (350–1050, 1000–1800, and 1800–2500 nm). The spectral resolution at 700 nm is 3 nm and the spectral resolution at 1400 nm and 2100 nm are both 8 nm. The wavelength accuracy is 0.5 nm and the field of view of the probe is 25°. When observing, the probe is vertically downward at a height of 1 m from the snow surface. The reflectance spectra were obtained automatically as the ratio of the target radiance and a standard white reference (Spectralon), multiplied by the Spectralon absolute spectral reflectance. Dark current and white reference were taken before each reflectance measurement. Each spectrum consists of an average of 40 individual spectra that were collected between 11:00 and 14:00 local time under clear-sky conditions and calm wind. The maximum standard deviation throughout the spectrum is 0.05.

A crystal gauge with a 40× handheld magnifier was used to measure the grain size and particle structures. At each observation point, the narrowest and longest widths of 80 snow particles were measured, and the average value was taken as the average effective optical-snow particle size [51]. The average grain size calculated using the digital snow particle property (DSPP) method was recorded as the average effective optical-snow particle size that can be simplified as the half of the long axis of the snow particle [52]. The snow density was calculated with a Snow Fork analyzer. The surface water content of snow cover was estimated with a snow cylinder. Figure 4 displays reflectance observations made using an ASD spectroradiometer, grain size measurement by crystal gauge with a magnifying hand lens, snow density estimation with a Snow Fork analyzer, and measurement of snow surface water content with the snow cylinder.

3.2. “Time-Gap Searching” for Snow-Aging Samples

The “time-gap searching” method is proposed in this study for collecting the snow-aging samples. This method can search for valid snow-aging samples through fixed locations from dual images with a short time gap. Thus, by sampling the pixels with the same locations and short time intervals, the influence from external time-affected variables (e.g., sensor zenith angle, solar elevation angle, etc.) and spatially-affected variables (e.g., undulating topography, shade effect, different slope/aspect) are reduced, thereby highlighting the temporal changes in snow reflectance associated with snow aging. Such a method is particularly suitable for the conditions of this study, for example, the research subject changes rapidly in the short term in this research area, and conditions for in-field sampling are not conducive.

In Figure 5, a pair of HJ-1A/B images indicating a sudden rise in snow-cover extent within a three-day interval are displayed. Considering that this short time interval was during the melting period, the sudden rise in snow area was caused by the fall of new snow that covered the aging snow, which was declining in the previous ablation period. Thus, 6 April 2014 was at the late stage of an ablation period, and 9 April 2014 was at the onset of the next period. By using the images of these two days, the difference of aging degrees of samples can be distinguished to a large extent. Therefore, the pixels measured at the same locations on 6 and 9 April 2014 were used to represent aging and fresh snows, respectively.

The snow samples analyzed on 6 April 2014 were graded as medium and severe aging by considering the average snowline altitude of 3900 m as the threshold [50,53]. Therefore, it can be considered that the snow samples below the snowline altitude on 6 April 2014 are classified as severe aging samples, while the others are classified as medium aging samples. Likewise, the snow covers higher than 3900 m on 9 April 2014 was defined as the target zone for fresh snow. Therefore, a total of 40,710 pixels within 306 regions of interest in three categories were collected and equally separated into two groups for training and validation, respectively.

3.3. Analytic Hierarchy Process for Construction of the SAI

The AHP synthesized multiple well-behaved aging indexes to construct an integrated SAI, which can calibrate the numeric scale for the quantitative and qualitative analyses [54]. Figure 6 shows the hierarchical structure in terms of goal, criteria, and alternatives. The goal is to construct the SAI B. The Jeffries–Matusita (J–M) separability criterion was set for assessing the statuses (P1, P2, and P3) of the samples. The optimal snow-aging indexes were considered alternatives. The intensity of importance was determined by pairwise comparisons according to the Saaty 1–9-scale method [54]. The criteria of P1, P2, and P3 was defined as equally important for the B (SAI). The rating of alternatives of S1 and S2 were pair-compared and defined with reference to the J–M separability of indexes for samples, and the weight factors of optimal snow-aging indexes were determined. The formula for SAI is given in Equation (1), where, w represents the weight factors of each optimal index, and the n represents the number of the selected index.

SAI= w₁ × Snow-aging index₁+ w₂ × Snow-aging index₂+…+ w_n × Snow-aging index_n

(1)

3.4. Elimination of the Terrain-Induced Shadow Effect

The SAI, which is based on optical images, is likely to be affected by an undulating terrain, and hence, biased in severely shaded and bright hot spot areas. Therefore, the principle of the shadow-adjusted snow-aging index is based on the implementation of a nonlinear combination of optical data to compensate for the SAI values of the severely shaded areas and suppression of the SAI of bright hot spot areas without affecting them in the other areas.

Therefore, the shaded SAI (SSAI) is established (see Equation (2)) based on the reflectivity differences between the red band of multispectral images for snow-cover indices under severely shaded areas, bright hot spot areas, and other flatter areas. This index highlights the difference of SAI values between the severely shaded and hot spot areas, and is insensitive within other flatter areas. Then, by employing a terrain condition adjusted factor, $f(\Delta )$, and constructing a nonlinear combination of SAI, the shaded and hot spot effects can be controlled. SSAI is defined as:

$$SSAI=\frac{M_r-B_r}{B_r}$$

(2)

where, for SAI calculated from every period of multispectral images, $B_r$ represents the reflectance data of the red band within the snow-covered area and $M_r$ the maximum reflectance value of the red-band data. To obtain SASAI, the terrain condition adjusted factor, $f(\Delta )$, is used to control the effect of shadows.

SASAI is defined as:

$$SASAI=SAI+f(\Delta )·SSAI$$

(3)

The iterative algorithm as an “extreme value optimization” method is introduced to determine $f(\Delta )$ as follows:

(i) Selection of the experimental area where the terrain shadow effect is significant. In the study area, the area with snow cover was extracted and considered as the experimental area.

(ii) Image classification of the experimental area for the classes: severely shaded, hot spot, and flat areas. The support Vector Machine supervised classification method was used for this purpose.

(iii) $f(\Delta )$ was determined using the iterative algorithm as follows: Let $f(\Delta )$ start from 0, increase in increments of 0.001, and calculate the maximum values of the SASAIs of the severely shaded and hot spot areas according to Equation (3). Exit the loop when these two maximums are equal or approximately equal, and determine the value of $f(\Delta )$. If no optimal value is found when the $f(\Delta )$ value is accumulated to 5, the loop is terminated. Reclassify the severely shaded and hot spot areas and re-iterate the search for the $f(\Delta )$ within the reclassified areas.

The premise of shadow adjustment is that the SAI values for severely shaded and hot spot regions are equivalent or close to each other. Thus, $f(\Delta )$ is iteratively determined until the maximum number of SASAIs are close to each other in the shaded and hot spot areas. Localized difference of energy fluxes or precipitation will cause an uneven aging of the snow cover. To preserve this factual spatiotemporal heterogeneity, the subsequent adjusted parameters within the same melting period were set as that of the initial value.

To evaluate the success of SASAI in correcting the effects of shading, the Pearson correlation coefficients of SAI and SASAI with respect to the solar incidence cosine were compared under various terrain conditions. The solar incidence cosine is defined by Equation (4):

$$\cos i=\cos \sigma \times \cos \theta +\sin \sigma \times \sin \theta \times \cos \left(\beta -\omega \right)$$

(4)

where i represents the solar angle of incidence of each pixel, $\sigma$ represents the slope angle, $\beta$ represents the aspect angle, and $\theta$ and $\omega$ represent the solar zenith angle and solar azimuth angle.

4. Results

4.1. Field Measurements and Analyses

Representative samples of new, uncontaminated aging, and contaminated aging snows were compared in terms of their morphological appearance, measured field spectra, and environmental parameters. As shown in Figure 7, a visual comparison of the photographs of magnified snow particles revealed that the colors of the new and uncontaminated aging snows are bright and clear, while the contaminated aging snow has a dim color. The fresh snow particles have morphological features such as hexagonal plates, dendritic shapes, and solid prisms. The aging snow particles show larger spherical aggregation in both the contaminated and uncontaminated snows. The detailed appearance characteristics and environmental parameters of snow particles are listed in

Table 1. Comparison of properties of new, aging, and contaminated snows.

	Fresh Snow	Aging Snow	Contaminated Snow
Color	Bright and clear	Bright and clear	Dull and impure
Shape (as shown in Figure 7)	No crystallization, snowflake, hexagon, columnar, or a mixture	Spherical and agglomerated	Spherical and agglomerated
Equivalent Grain size	200 nm–300 nm	600 nm–1000 nm	750 nm–1200 nm
Average layer water content	1.2–1.6%	2–4%	6–9%
Average density	176–221 kg/m3	320–90 kg/m3	340–420 kg/m3
Average Air temperature	−9.8 ℃	−2.5 ℃	7.8 ℃
Sampling date	10–17 December 2013	15–22 March 2014	15–22 March 2014
Time lag from snow fall	~1 day	> 7 days	> 7 days

. With increase in time after snowfall, the aging snows have a larger particle size, water content, and density. Furthermore, the aging of snow corresponds with rising temperature and time lag.

The field spectra of new, uncontaminated aging, and contaminated aging snows were compared (Figure 8a). A distinct spectral feature in all the shapes were found in the wavelength range of 350–1300 nm. Compared with the high reflectivity of clean fresh snow (maximum ranging between 0.5–1.0), the drop in the reflectance of clean aging snow might be caused by the scattering characteristics owing to an increase in particle size. A further significant drop in the reflectance of contaminated aging snow might be because of absorption of radiation by impurities that accelerate aging.

The reflectivity of uncontaminated aging snow with three observed particle sizes was compared with that of fresh snow (Figure 8b). Lower absolute reflectance values (decline ranging between 0.15 and 0.25) were measured in accordance with increasing grain size. The acquired reflectance shows the typical downward trend when the grain size increases within a wavelength range from 350 to 1050 nm.

Figure 8c displays the reflectivity of the contaminated aging snow with increasing distances from the highway, a source of pollution. Distinct lower absolute reflectance values (max ranging between 0.3 and 0.4) were measured in the polluted snow that is closest to the highway. The reflectivity gradually rises with increasing distance from the highway.

Comparison of snow reflectivity with water content (Figure 8d) shows that reflectance decreased with an increase in the water content of the aging snow. Therefore, the reflection measurements at wavelengths sensitive to the water content will help characterize the aging of snow cover.

Generally, snow surface reflectance can effectively distinguish between fresh, aging, and contaminated aging snows, especially, within the wavelength range of 430–1300 nm. An increase in particle size and particle aggregation are significant features of snow aging. The pollutants will significantly accelerate the melting progress. Similarly, an increase in water content and density are associated with snow aging. Air temperature is also an indicator of snow aging. Thus, these indicators were cataloged and a candidate SAP library (

Table 2. The candidate SAP indexes library with respective calculation formulas.

Index Types	Index Name	Calculation Formula	N/A	N/M	M/A
Snow cover indicator	Normalized Difference Snow Index [55]	NDSI =$\frac{\mathrm{Green}-\mathrm{NIR}}{\mathrm{Green}+\mathrm{NIR}}$	0.65	0.95	0.66
	Melt Area Detection Index [56]	MADI =$\frac{\mathrm{Green}}{\mathrm{SWIR}}$	1.10	0.98	0.93
Grain size indicator	Snow Grain Index [57]	SGI =$\frac{\mathrm{Green}-\mathrm{IR}}{\mathrm{Green}+\mathrm{IR}}$	1.29	1.41	1.41
Contamination indicator	Snow Contamination Index [58]	SCI = $\frac{\mathrm{Blue}-\mathrm{Green}}{\mathrm{Blue}+\mathrm{Green}}$	0.81	0.75	0.67
Vegetation indicator	S3 [59]	S3 =$\frac{\mathrm{NIR}\left(\mathrm{RED}-\mathrm{SWIR}\right)}{\left(\mathrm{NIR}+\mathrm{Red}\right)\left(\mathrm{NIR}+\mathrm{SWIR}\right)}$	0.69	1.06	0.61
Temperature indicator	Normalized Difference Temperature Snow Index [60]	NDTSI =$\frac{\mathrm{NDSI}-\mathrm{TIR}}{\mathrm{NDSI}+\mathrm{TIR}}$	1.08	0.98	0.70
	Normalized Difference Snow Temperature Index [60]	NDSTI =$\frac{\mathrm{Blue}-\mathrm{TIR}}{\mathrm{Blue}+\mathrm{TIR}}$	1.41	1.05	1.11
	S3 Temperature Snow Index [60]	S3TSI =$\frac{\mathrm{S}3-\mathrm{TIR}}{\mathrm{S}3+\mathrm{TIR}}$	0.78	1.23	0.80
Single band reflectance	Reflectance of HJ-1A/B CCD band 1	Ref Blue	1.41	1.18	1.15
	Reflectance of HJ-1A/B CCD band 2	Ref Green	1.41	0.97	1.08
	Reflectance of HJ-1A/B CCD band 3	Ref Red	0.73	1.30	1.33
	Reflectance of HJ-1A/B CCD band 4	Ref NIR	1.41	1.11	0.96
HJ-1A/B CCD with spatial resolution of 30 m		HJ-1B IRS with spatial resolution of 120 m
Blue: Band 1 (0.43–0.52 μm)		NIR: Band 5 (0.75–1.10 μm)
Green: Band 2 (0.52–0.60 μm)		SWIR: Band 6 (1.55–1.75 μm)
Red: Band 3 (0.63–0.69 μm)		MIR: Band 7 (3.50–3.90 μm)
NIR: Band 4 (0.76–0.90 μm)		TIR: Band 8 (10.5–12.5 μm)
N/A: New & Aging snow		M/A: Medium & Aging snow
N/M: New & Medium aging snow		Grade of separability: Excellent/Good/Medium/None

) was constructed.

4.2. Optimized Snow Cover Aging Indexes

By combining the snow-aging training samples, the J-M distance between each pair of sample categories was calculated for every snow cover aging index in the candidate SAP index library (Table 2). As shown in Table 2, the results revealed that only the J–M distances of Snow Grain Index (SGI) and HJ-1A/B CCD band 1 (CCD-1) for every pairwise sample was larger than 1, thus reaching a ‘good’ or ‘excellent’ separability grade. All the other indexes have a separability grade of ‘medium’ or ‘none’ for at least one pair of samples. Snow aging within the MRB is likely to have the maximum differentiation in SGI, followed by reflectance in CCD-1. Therefore, SGI and CCD-1 are the preferred indicators for snow aging.

4.3. Snow-Aging Index (SAI) and Shadow-Adjusted Snow-Aging Index (SASAI)

The weight factors for the optimal snow-aging indexes were calculated using AHP and were determined as 0.583 and −0.417 for SGI and the blue band (CCD-1), respectively. The weight factor of the blue band was set as negative to quantitatively express the continuous snow melting process. From Figure 9a–f, it can be inferred that: (i) The overall values of SAI and SASAI on 6 April 2014 were higher than that on 9 April 2014, which indicates that both the SAI and SASAI have the ability to distinguish the temporal dynamics of snow aging; (ii) the decreasing trend of SAIs with increasing altitude is obvious on 6 and 9 April 2014, especially, in areas adjoining the snowline. This trend was significantly weakened in SASAIs. On 9 April 2014, no change with altitude was observed for SASAI. This is consistent with the characteristics of snow melting rate decreasing with low temperature at high altitudes, and homogeneity of fresh snow during the snowfall period; (iii) SAI has significant stripe-shaped differences in the shadow and bright spot areas near the ridgeline, and the optimized SASAI significantly weakens this difference. SASAI adjusts the high values obtained for hot spot areas in SAI with lower values, and conversely, the lower values obtained for shaded areas are enhanced. This shows that SASAI is superior to SAI in terms of indicating the degree of snow aging.

4.4. Assessment of the Shading Correction and Accuracy of SAI/SASAIs

To evaluate the shade correction effect of SAI and SASAI indexes, Pearson correlation coefficients between the SAIs and the LSIC (see Figure 9g) under different terrain types (see Figure 9h) were compared within a subarea of MRB (shown in the red rectangle in Figure 9c,d). As shown in Figure 10, the overall correlation coefficient R between SAI and LSIC is 0.68. In severely shaded, hot spot, and flat areas, the corresponding R values are 0.921, 0.89, and 0.53, respectively. There is a positive correlation between SAI values and LSIC, especially, in severely shaded and hot spot areas suggesting a directly proportional relationship between them. Comparatively, the overall correlation coefficient R between SASAI and LSIC is 0.19. In the three terrain types, the corresponding R values are 0.15, −0.27, and 0.15, successively. The remarkable decrease in these R values highlights the effectiveness of the shading correction on SAI.

A scatter plot was made to validate the accuracy of the SAI/SASAI with three types of validation samples (Figure 11a). The aim was to demonstrate whether the proposed SAI or SASAI can indicate a significant insight on snow aging. Comparatively, the scatter plot of SASAI shows better separability than SAI. SAI values have a greater degree of overlap when compared with the corresponding SASAI values of the three types of samples. Figure 11b,c show the box plots of SASAI and SAI values for the three types of validation samples, and they indicate that the upper and lower quartile values of SAI and SASAI do not overlap, but the three types of SASAI sample indexes are significantly more concentrated and better distinguished.

To improve the confidence level and explore the adaptability of the proposed algorithm on multiperiod remote sensing images, the average values of SASAI and SCA of the 14 periods during the snowmelt of 2014 were calculated based on HJ data. As shown in Figure 12, SASAI and SCA have a significantly negative correlation. Each decline cycle is characterized by a gradual decrease in snow and an increase in SASAI. After each new spell of snowfall, the SCA suddenly increased, and the SASAI decreased significantly. Therefore, it can be inferred that the snow-aging index proposed in this paper can identify the snow-aging status. The R² between the SASAI and the SCA calculated for 14 periods of data is 0.728, which indirectly illustrates the acceptability of SASAI.

5. Discussion

The snow cover aging process refers to a series of changes in snow cover such as declining area, decreasing depth, and corresponding physical parameters that appear from the end of the snowfall stabilization period to the transition period, and last until the end of the ablation period. The algorithm proposed in this study is based on optical remote sensing and helps in obtaining the status of the snow surface and identifying the changes from fresh snow to severely aging snow. Field work was conducted in two routes within the MRB during winter and spring to collect representative samples. The in-field measurements distinguish the new and old snow based on comprehensive evaluation of particle size, water content, pollution, and time lag from the snowfall day. However, the aging factors and their corresponding effects vary in different seasons and regions. Therefore, reliable aging indexes from previous studies were added to complete the candidate SAI library. The process of snow aging, which exhibits various features, is a compound result of multiple parameters. The expression of SAP in spatiotemporal scale needs to be explored in subsequent studies.

Remote sensing images from two periods were used to extract training and verification sample sets of three types of snow aging. Because of the short time interval between these two periods of images and their spatiotemporal typicality, they are best suited to reflect the abrupt changes in extremely aged snow and fresh snow. Hence, the constructed SAI can be applied to other regions with similar climatic and geographic characteristics. However, the causes for the aging during the snow ablation period can vary from place to place. It is suggested that rescreening of the optimal snow-aging index at different decline periods is necessary for better understanding.

In this study, SGI and blue-band reflectance were optimized to constitute the SAI (see Equation (2)). It can be suggested that aging within the MRB is indicated by the increase in particle size of snow and the decrease in reflectance of blue band. Field work confirmed that the snow is uncontaminated in most of the study area, and hence, the SCI does not discriminate the training samples more than that by SGI and blue-band reflectance. The NDSTI demonstrated good discrimination of training samples because its components (e.g., reflectance in blue band and thermal infrared band) elucidate decreasing shortwave reflectivity and temperature rises in characterizing snow aging. In recent satellites, the applicability of the thermal infrared band is limited, and hence, it was not incorporated in the SAI algorithm. However, the optimal index can be determined for sensors with different satellites. Furthermore, the bands required for calculating SAI (i.e., blue band and SGI) are available in the existing multispectral optical satellites. Samples of this study are representative and can be utilized for other regions as they have the same waveband setting and an approximate bandwidth that can be transferred to multi-sensor platforms such as Landsat-8 and Sentinel-2.

It is worth mentioning that the subsequent adjustment of SAI is based on the premise that the fresh snow is uniformly distributed. Theoretically, the SAI values within every snow-covered area lie within a narrow range. Thus, by enhancing the SAI of the severe shadow area and suppressing the SAI of the hot spot area, the error can be effectively eliminated. The introduction of the SSAI makes the adjustment of the SAI value of the flat area negligible, which can minimize the over-adjustment of the true value of the SAI. However, many factors such as slope and local terrain, precipitation, cloud cover, solar azimuth, etc. may affect the energy distribution within a region and thus, these may cause redistribution of snow aging. The combined effect of these factors on shaping the heterogeneity of snow aging is very complex, and hence, these were not simulated in this study.

In general, based on the proposed method, it is feasible to build a snow cover aging data set. This data set can macrocharacterize the new or old status of snow cover and reflect on the cause and degree of aging by referring to the components and weights of SAI. This data set will complement the existing binary or fractional snow cover analysis modules, and can also provide important input parameters, especially for hydrological (snowmelt–runoff) and climate models.

6. Conclusions

Snow recession is an important component of the global climatic–hydrological processes and it is highly dynamic spatiotemporally. Owing to shortcomings in existing research on the acquisition of snow-aging parameters, a method to grasp the spatiotemporal aging process of snow recession was established. As a comprehensively and quantitatively integrated snow index, the SASAI can describe the snow-aging information induced by various factors and realize the continuous description of the new landing snow, ablation, and exhaustion statuses. This index has a mathematical expression based on the inherent properties of reflectance spectra of different snow statuses that are arrived from the libraries of candidate indexes (i.e., SGI and reflectance of blue band). Using the “time-gap searching method”, the SASAI is parameterized by an AHP calculation using the training data. Using the “extreme value optimization” algorithm, the SASAI removes the effect of shading in the rugged terrain. Assessments of the shaded-area correction, separability of the SAI and SASAI, and adaptability of this algorithm on multiperiod remote sensing images were conducted successively with satisfying results. Considering the rapid snow ablation process, additional evaluations of the SASAI on multiple satellites and multiple sensors will be conducted in the future. The calibration of parameters with component weights of the SAI will be the next step. A snow-aging state field will be established on this basis, thereby facilitating a more refined modeling of the snowmelt–runoff process. The proposed method is reliable, and it can be utilized for identifying the spatiotemporal aging status and distribution of snow ablation.