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

: Detecting the variations in snow cover aging over undulating alpine regions is challenging owing to the complex snow-aging process and shadow e ﬀ ect 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-ﬁeld measurements were used during the snow ablation and accumulation periods. The in-ﬁeld 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 ﬁnally corrected for the e ﬀ ects of shading and assessed. This study provides both a ﬂexible 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.


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 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.

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.
Remote Sens. 2019, 11, x FOR PEER REVIEW 3 of 18 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.

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.

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.

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. 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.  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.

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)(16)(17)(18)(19)(20)(21)(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

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)(16)(17)(18)(19)(20)(21)(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 Remote Sens. 2020, 12, 1249 5 of 18 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.
Remote Sens. 2019, 11, x FOR PEER REVIEW 5 of 18 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.

"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. 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.
Remote Sens. 2019, 11, x FOR PEER REVIEW 5 of 18 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.

"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

"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.
Remote Sens. 2019, 11, x FOR PEER REVIEW 6 of 18 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.

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 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 Remote Sens. 2020, 12, 1249 7 of 18 40,710 pixels within 306 regions of interest in three categories were collected and equally separated into two groups for training and validation, respectively.

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.

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 shadowadjusted 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, (∆) , and constructing a nonlinear combination of SAI, the shaded and hot spot effects can be controlled. SSAI is defined as: where, for SAI calculated from every period of multispectral images, represents the reflectance data of the red band within the snow-covered area and the maximum reflectance value of the red-band data. To obtain SASAI, the terrain condition adjusted factor, (∆), is used to control the effect of shadows.
SASAI is defined as: The iterative algorithm as an "extreme value optimization" method is introduced to determine (∆) 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

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 (∆), and constructing a nonlinear combination of SAI, the shaded and hot spot effects can be controlled. SSAI is defined as: 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 (∆), is used to control the effect of shadows. SASAI is defined as: Remote Sens. 2020, 12, 1249 8 of 18 The iterative algorithm as an "extreme value optimization" method is introduced to determine f (∆) 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 (∆) was determined using the iterative algorithm as follows: Let f (∆) 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 (∆). If no optimal value is found when the f (∆) 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 (∆) 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 (∆) 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): where i represents the solar angle of incidence of each pixel, σ represents the slope angle, β represents the aspect angle, and θ and ω represent the solar zenith angle and solar azimuth angle.

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. 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. particles revealed that the colors of the new and uncontaminated aging snows are bright and clear, while th 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.  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 (

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). 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.

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.

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 2 between the SASAI and the SCA calculated for 14 periods of data is 0.728, which indirectly illustrates the acceptability of SASAI.   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. Figure 10. Correlation coefficients between the SAI/SASAI and solar incidence cosine in different terrain types.
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 2 between the SASAI and the SCA calculated for 14 periods of data is 0.728, which indirectly illustrates the acceptability of SASAI.

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

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

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.

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. Funding: This work is a part of the project titled, "Key technology research on the monitoring of the snow cover degradation process in the Alpine mountains using optical remote sensing data", which is supported by the National Natural Science Foundation of China (grant number 41801269) for a period of 3 years (1 January 2019-31 December 2021).