Open AccessArticle

Snow Cover Variability and Trends over Karakoram, Western Himalaya and Kunlun Mountains During the MODIS Era (2001–2024)

Cecilia Delia Almagioni

^1,2,

Veronica Manara

^1,*

Guglielmina Adele Diolaiuti

Maurizio Maugeri

Alessia Spezza

^1,3

and

Davide Fugazza

Environmental Science and Policy Department, University of Milan, via Celoria 10, 20133 Milan, Italy

Department of Physics “Aldo Pontremoli”, University of Milano, via Celoria 16, 20133 Milan, Italy

Environmental Sciences, Informatics and Statistics Department, Ca’ Foscari University, via Torino 155, 30172 Mestre, Italy

Author to whom correspondence should be addressed.

Remote Sens. 2025, 17(5), 914; https://doi.org/10.3390/rs17050914

Submission received: 14 January 2025 / Revised: 28 February 2025 / Accepted: 28 February 2025 / Published: 5 March 2025

(This article belongs to the Section Environmental Remote Sensing)

Download

Browse Figures

Versions Notes

Abstract

Monitoring the snow cover variability and trends is crucial due to its significant contribution to river formation and sustenance. Using gap-filled MODIS data over the 2001–2024 period, the spatial distribution and temporal evolution of three snow cover metrics were studied: number of days, onset and end of the snow cover season across fourteen regions covering the Karakoram, Western Himalayas and Kunlun Mountains. The obtained signals exhibit considerable complexity, making it difficult to find a unique factor explaining their variability, even if elevation emerged as the most important one. The mean values of snow-covered days span from about 14 days in desert regions to about 184 days in the Karakoram region. Given the high interannual variability, the metrics show no significant trend across the study area, even if significant trends were identified in specific regions. The obtained results correlate well with the ERA5 and ERA5-Land values: the Taklamakan Desert and the Kunlun Mountains experienced a significant decrease in the snow cover extent possibly associated with an increase in temperature and a decline in precipitation. Similarly, the Karakoram and Western Himalayas region show a positive snow cover trend possibly associated with a stable temperature and a positive precipitation trend.

Keywords:

snow duration; spatial distribution; trends; MODIS data; ERA5 data; Karakoram; Himalayas; Kunluns

1. Introduction

Snow cover is a crucial component of the Earth’s climate system, with its extent and duration affected by climate-change-induced temperature rise [1]. Numerous studies have indeed evaluated the reduction in snow-covered areas in the Northern Hemisphere [2,3] and at a global scale, reporting that snow-cover extent decreased by 11.8% between 1967 and 2012 [4]. Besides the reduction in snow cover, glaciers at mid-latitudes have been retreating, with over 600 glaciers disappearing since 1894 [4] and with global cumulative ice loss exceeding 20 m water equivalent between 1945 and 2005 [5].

The decrease in snow cover, together with the reduction of glaciers and sea ice, alters the Earth’s albedo, affecting the Earth’s ability to reflect incoming solar radiation and triggering a positive feedback mechanism [6,7]. Moreover, it can have a strong impact on water resources, as in many areas snowmelt feeds a relevant fraction of the water reservoirs, supporting various essential sectors, including drinking water supply, agricultural irrigation, industrial processes, and energy production [7,8,9].

In remote and mountainous regions with limited observational records, satellite data are essential for monitoring snow cover [9,10,11,12]. Since 1967, NOAA satellites have provided long-term snow cover records [13], though with coarse resolution (200 km) and potential underestimation, particularly in complex terrain. Since the mid-1980s, AVHRR [14] sensor data with 1 km resolution have been used to study snow variability, though early datasets faced challenges in distinguishing snow from clouds [15]. Since the launch of NASA’s Earth Observing System mission satellites Terra and Aqua in December 1999 and May 2002, respectively, equipped with the MODIS (Moderate Resolution Imaging Spectroradiometer) sensor, snow cover analyses have made wide use of MODIS products, such as MOD10A1 [16,17,18]. The advantages of MODIS products are several. Hall and Riggs [19] have presented a comprehensive review of studies assessing the accuracy of MODIS snow products, which highlight a rather good agreement (94.2% on average) when compared with corresponding ground data [20]. This accuracy varies seasonally on mountain terrains, with higher errors in winter and lower ones in spring and summer [21,22]. Moreover, the spatial resolution of 500 m and the passage of two satellites per day make the MODIS data an excellent product for snow cover monitoring.

Monitoring the Himalayas and Karakoram region is crucial due to its significant contribution to river formation and sustenance, supporting nearly two billion people through snow and ice melt [23]. Indeed, approximately 50% of the Indus Irrigation Scheme relies on runoff from snowmelt and glacier melt in the Eastern Hindukush, Karakoram, and Western Himalayas [24]. Moreover, the economy of these regions relies upon agriculture, and thus, it is highly dependent upon water availability and irrigation systems [25,26,27]. Additionally, in the context of global warming, the Karakoram region presents a unique anomaly as it has not experienced an overall negative glacier mass balance in the last decades [28,29]. Glaciers in this region seem, in fact, to exhibit stability and even gain mass and thickness; this phenomenon is known in the literature as the “Karakoram anomaly” [23,29,30]. This anomaly is mirrored by a similar anomaly for snow cover [31].

Several studies utilizing MODIS data have investigated snow cover variability in Gilgit Baltistan, Northern Pakistan. They highlight a complex picture: Atif et al. observed an increasing trend for the period 2001–2015 [32], whereas other studies focusing on specific subregions have reported an opposite trend. These include the Upper Indus Basin for the period 2001–2015 [33], the Hunza district for the period 2003–2018 [34], and the Hunza, Shigar, and Shyok districts for the period 2003–2018 [35]. Also, in the Tarim Basin, a decline in snow-covered areas is observed, with varied trends across different subregions [36]. Conversely, in the Northwestern Himalayas, an increase in snow cover area has been reported for the period 2000–2019, primarily in winter [37]. Moreover, Yi et al. [38] investigated the length of the snow season in the Western Upper Indus and Yarkand River basins from 2002 to 2018, reporting an overall lack of significant trends. Situated on the southern and northern sides of the Karakoram range, respectively, these basins exhibit distinct climatic conditions, as emphasized by the authors. Besides studies focusing on specific areas, more comprehensive studies of snow variability patterns in recent decades over High Mountain Asia (HMA) were conducted. Ackroyd et al. [39] a broad snow cover decline (2001–2017). Dharpure et al. [40] found an increasing trend across the Karakoram–Himalayan range, except for the Eastern Himalayas in 2001–2019. Tang et al. [41] observed a slight decreasing trend in monthly snow cover extent across most months in 2002–2020. Xu et al. [42] reported a general decline in the snow season over the Tibetan Plateau in 2002–2022. The first study indicated a broad snow cover decline, the second one an increasing trend for the entire Karakoram–Himalayan range except for the Eastern Himalayas, and the third study indicated that the monthly snow cover extent over HMA exhibited a slight decreasing trend in most months, the latter study reported a general decrease in snow season in the Tibetan Plateau. Other studies analyzed the variation of the Snow Line Altitude at the End of the Melting Season (SLA-EMS) in HMA, with upward trends of it both in the periods 1991–2022 [43] and 2001–2016 [44].

The complexity of the trend signals reported in these papers can indeed be attributed to the diverse factors influencing snow cover in the considered areas; however, the limited length of available time series is also a relevant factor as with such short series, the impact of a single or a few years with low or high snow coverage at the beginning or at the end of the investigated period may be crucial for the results of the trend analysis. Because of these problems, it is important both to investigate the snow cover signals dividing the study area into defined geographic regions and to verify if the observed trends are confirmed by an updated version of the MODIS snow cover records.

This study aims to investigate the spatial distribution and temporal evolution of the snow cover season over the Karakoram, the Western Himalayas, and the Kunlun Mountains. Unlike previous studies, which have often considered either large regions or small sub-regions, a basin-based approach is used, dividing the study area into 14 sub-regions to provide a more detailed spatial characterization. Additionally, an extended MODIS dataset, updated to September 2024, is used, allowing for a more comprehensive assessment of recent snow cover trends. Furthermore, a multi-step gap-filling procedure is applied to mitigate the impact of cloud cover, enhancing data reliability. The snow cover season is characterized by the snow-covered days (SCD), the snow onset date (SOD), and the snow end date (SED), and the inter-and intra-annual variability of these metrics are investigated. Moreover, temperature, snow cover, and precipitation records from ERA5 [45] and ERA5-Land [46] reanalysis data are employed to better understand the results obtained from these snow cover metrics.

Following the Introduction (Section 1), the study area is described (Section 2.1), and the data (Section 2.2) and the methods used for the analysis are presented (Section 2.3 and Section 2.4). Then, the spatial patterns observed in snow cover onset, duration, and melt-out are presented (Section 3.1), the impact of topographic and geographic factors as drivers of snow cover variability are analyzed (Section 3.2), and the interannual variability and trends in snow cover metrics in the 2001–2024 period are analyzed (Section 3.3). Finally, the results are discussed, and the link between the snow cover metrics and the driving meteorological variables is investigated (Section 4), and some final remarks are presented (Section 5).

2. Materials and Methods

2.1. Study Area

The study region is defined by the MODIS tile h24v05: it is an oblique area within the rectangle defined by coordinates 69–91°E in longitude and 29–40°N in latitude (Figure 1). This area encompasses three countries: Pakistan, India, and China. The included regions are part of Gilgit-Baltistan, Kashmir, Himachal Pradesh, Uttaranchal, Punjab, Aksai Chin, Tibet, and Xinjiang. The orography in the area is complex and varied, featuring three distinct mountain ranges: the Karakoram, the westernmost part of the Himalayas, and a substantial portion of the Kunlun Mountains. To the north of the Kunlun range lies the Tarim Basin, including the Taklamakan Desert. To the south of the Karakoram and Himalayas lies the Indo-Gangetic Plain, specifically the Upper Indus Basin (Figure 1). To the north of the Himalayas, south of the Kunlun, and east of the Karakoram lies a portion of the Tibetan Plateau. Given the vastness of the area, it was divided into 14 subregions following a hydrological-based approach (Table 1). These regions cover rather different altitudinal ranges (Figure 2). Considering them all together, 47% of the grid points are situated at elevations between 0 and 2000 m above sea level (a.s.l.), 12% fall within the range of 2000 m.s.l. to 4000 m a.s.l., 40% are located between 4000 m a.s.l. and 6000 m a.s.l., and 1% extend beyond 6000 m a.s.l.

Given the extent and diverse nature of the study area, various climate influences are relevant within it. The most important include the mid-latitude westerlies, the Indian Summer Monsoon, and the Inner Asian high-pressure system [47].

The Karakoram climate is mainly influenced by the westerlies, which bring winter precipitation and persist into the summer months when the monsoon also plays a role. Indeed, monsoonal air and storms occasionally penetrate the Himalayan Front Ranges with sufficient force to bring significant precipitation to the Central Karakoram region [48,49]. The Indian Summer Monsoon does, however, weaken going northwest, making the intrusions rare on the northern sides of the mountains in the westernmost part of the range [48]. Most rainfall in the Karakoram occurs above 4000 m a.s.l., with snowfall ranging from 1000 mm to 2000 mm in water equivalent in glacier source areas, contrasting with 150–300 mm in the valleys [50]. Around 90% of annual precipitation falls as snow at altitudes near 5000 m a.s.l., while snow contributes less than 10% in lower valley regions [24].

The climate of the Himalayan and Northern Indus Plain is mainly influenced by the monsoon, which brings significant summer precipitation, whereas winter precipitation, brought by westerly winds, is more prominent in the Western Himalayas, particularly at higher elevations. Overall, total annual precipitation increases from west to east, with average annual rainfall on south slopes ranging from about 1500 mm in the west to about 3000 mm in the east [51]. As observed in the Karakoram, the lower elevation areas in the Western Himalayas exhibit a drier climate [52], with total precipitation differing by almost one order of magnitude between low altitudes and high altitudes [24].

The Himalayan range blocks the northward transport of water vapor, resulting in arid conditions on the Tibetan Plateau. The average annual precipitation there is around 460 mm, falling mostly during the summer months [53]. The western part of the Tibetan region has a drier climate with annual precipitation averaging between 100 and 250 mm, whereas the eastern part receives higher precipitation, ranging from 400 to 1000 mm annually [54].

The Kunlun Mountains are predominantly influenced by the continental air masses, experiencing significant temperature fluctuations annually and daily. The average annual precipitation of the area is below 300 mm, and it occurs mostly during summer months [55].

Most of the rainfall in the Tarim Basin occurs during June, July, and August. Mountainous areas receive precipitation influenced by the westerlies, with some areas receiving over 300 mm annually, primarily as snowfall. However, the basin’s average annual precipitation remains below 50 mm [56]. Situated in the middle of the Tarim Basin, the Taklamakan Desert exhibits an arid climate, with minimal precipitation averaging 12 mm annually, ranging from annual means of 38 mm in the west to 10 mm in the east [57].

To achieve a synoptical coverage of precipitation distribution over the study area, without being limited by the lack of observational data, the ERA5 reanalysis dataset was used [45], which has a spatial resolution of 0.25°. The total annual precipitation for each pixel was calculated and averaged over the period 1991–2020, and the precipitation distribution throughout the year for individual subregions was analyzed. The High Asia Refined Analysis dataset (HAR) was also analyzed [58], which is a dynamic downscaling of ERA5 using the Weather Research and Forecasting Model (WRF), and a similar distribution was found. The only notable difference was slightly more pronounced precipitation on mountain reliefs. Figure 3 and Table 1 highlight the variations in annual rainfall contributions using ERA5 data across the different areas and the coherence between the adjacent ones.

2.2. Snow Cover and Ancillary Data

For the snow cover analysis, MODIS satellite image products from October 2000 to September 2024 downloaded from the NASA Earth Data portal (https://earthdata.nasa.gov/, last accessed on 10 December 2024) were used. The products used are MODIS Terra Snow Cover Daily L3 Global 500 m SIN Grid V061 (MOD10A1 V61) and MODIS Aqua Snow Cover Daily L3 Global 500 m SIN Grid V061 (MYD10A1 V61), distributed by the US National Snow and Ice Data Center (NSIDC), which are daily products with a spatial resolution of 500 m [18]. The satellites operate in a polar orbit, crossing the equator at 10:30 local time descending southward for Terra and at 13:30 local time ascending northward for Aqua.

The algorithm that identifies snow on a pixel relies on the difference in reflectivity of different types of surfaces (e.g., snow, vegetation, clouds, or water). Snow has a higher reflectivity in the visible band and a high absorbance in the shortwave infrared. The Normalized Difference Snow Index (

N D S I

) utilizes these characteristics to measure reflectivity differences between the visible band (

B_{4}

, 0.545–0.565 µm) and the infrared band (

B_{6}

, 1.628–1.672 µm for MODIS Terra or B7, 2.105–2.155 µm for MODIS Aqua [19]) as follows:

N D S I = \frac{B_{4} - B_{6}}{B_{4} + B_{6}},

(1)

where

B_{7}

is substituted in place of

B_{6}

for the MODIS Aqua product. Pixels with

N D S I

≤ 0.0 are considered snow-free, whereas those with 0.0 <

N D S I

≤ 1.0 are considered at least partially snow-covered. However, in this study, which relies on a binary snow classification, only pixels with

N D S I

index between 0.4 and 1 were classified as snow-covered, as the 0.4 threshold is reported to indicate pixels that are snow-covered for at least 50% of their area [59]. The pixels with 0.0 <

N D S I

< 0.4 were therefore classified as snow-free.

Tests based on visible, shortwave, and longwave infrared channels reduce false detection errors [60,61]. Additional filters are applied to all pixels to highlight specific conditions that may impact the overall data quality or require further investigations.

In this work, the

N D S I

Snow Cover and the

N D S I

Snow Cover Basic QA layers were used. Only the data obtained in optimal conditions were retained, while poor quality data were categorized as missing. The

N D S I

Snow Cover layer classifies the pixel with its corresponding

N D S I

index in percentage if

N D S I

> 0, or as cloudy (250) or water covered (237 for internal waters and 239 for oceans), while other values represent error conditions.

A water mask combining three different masks was developed because for inland waters the algorithm often produced unreliable data, and no single mask was sufficient for complete water categorization. The first mask was created by identifying all points marked as water in at least one year by the

N D S I

Snow Cover product as water covered over the whole period. To be classified as water-covered in one year, a pixel should have at least 10% of the days with available data throughout the year as water-covered. The second mask was derived from the MODIS land cover yearly product MCD12Q1, which includes several classification schemes to distinguish between land cover types. If a pixel was classified as water in any scheme, it was considered to be water-covered for that year. If a pixel was water-covered in any year, it was considered water-covered for all years. The third water mask was downloaded from the Diva-Gis data store (https://diva-gis.org/ last accessed on 10 September 2023) in shapefile format and converted to raster with the same grid size as the MODIS products. A grid point that was considered water-covered by any of the three masks was removed from the analysis.

Lastly, ERA5-Land [46] records for monthly means of snow cover area and air temperature with a spatial resolution of 0.1° and ERA5 [45] records for monthly precipitation with a spatial resolution of 0.25° (Figure 3) were downloaded from the Copernicus data store (https://cds.climate.copernicus.eu, last accessed on 10 December 2024).

The elevation of the grid points from the ALOS (Advanced Land Observing Satellite) Global Digital Surface Model “ALOS World 3D-30 m (AW3D30)” Digital Elevation Model [62] was further obtained for the entire study area with a resolution of 30 m and resampled the Digital Elevation Model (DEM) to 1/240° for consistency with the MODIS products. From the resampled DEM, the slope and aspect to analyze the contribution of the topography to the variability in snow cover were derived. Latitude and longitude of each pixel were extracted using the grid centroid.

2.3. Cloud Cover and Gap Filtering

MODIS snow products are obtained from data in the visible and near-infrared spectral range, making it impossible to observe surface conditions beneath clouds. The analysis of this work revealed high percentages of cloud cover over the study area, especially in mountainous regions, where the percentage of missing data could exceed 80% after combining the MOD10A1 and MYD10A1 datasets. Several methods have been proposed for cloud removal and assigning values to the obscured grid points [63,64,65,66]. The approach proposed by Gafurov and Bárdossy [67], and adapted by Fugazza et al. [68] was adopted. To filter cloud cover, each MOD10A1 and MYD10A1 scene was reclassified into three categories: snow-covered, snow-free, and cloud-covered. Grid values from 40 to 100, where the number represents the NDSI value multiplied by 100, were categorized as snow-covered (assigned a value of 1) [69]; land, ocean, and water areas were marked as snow-free (value of 0); cloud-covered areas and error values were altogether given a value of 3. Any data gaps were treated as days with 100% cloud cover. Gafurov and Bárdossy’s method involves six steps with progressively broader spatial and temporal interpolations to ultimately produce a binary grid where each point is classified as either snow-covered or snow-free. The detailed steps followed are described in the work of Fugazza et al. [68]. While similar approaches have been widely used and validated in previous studies [67,69,70,71], it is important to acknowledge their potential limitations. The accuracy of the gap-filling method is inherently dependent on the assumptions made during interpolation, and in regions with persistently high cloud cover, errors may propagate due to the lack of reliable reference data. Additionally, complex terrain effects, such as topographic shading and heterogeneous snow distribution, may introduce further uncertainties that are difficult to correct without independent ground-based validation.

2.4. Calculation of Snow Cover Metrics

To describe the spatial and temporal extension of annual snow cover, we considered 24 hydrological years in the period October 2000–September 2024. A hydrological year starts on 1 October ends on 30 September, taking the year of its end as its reference. Throughout this study, all years mentioned refer to hydrological years.

Three metrics were used: snow-covered days (SCD), snow onset date (SOD), and snow end date (SED). The snow-covered days SCD is defined as the number of days in which a certain grid point is covered by snow in a hydrological year (starting on the 1st of October and ending on the 30th of September). It is calculated with the following equation:

{S C D}_{i, j} = \sum_{k = 1}^{n} S_{i, j, k},

(2)

where k is the day of the hydrological year, n = 365 (or 366 for leap years) and

S_{i, j, k}

assumes values of 0 if the grid point is snow-free and 1 if it is snow covered. The start and end of the snow cover season were defined based on the fixed date method that was described in the work by Wang and Xie [1]. The snow onset date SOD is defined with respect to the number of snow cover days between the 1st of October and the 25th of January. The SOD is calculated as follows:

{S O D}_{i, j} = D - \sum_{k = 1}^{D} S_{i, j, k},

(3)

where D = 117 days which is the number of days between the 1st of October and the 25th of January. The snow end date SED is defined with respect to the number of snow days between the 26th of January and the 30th of September as:

{S E D}_{i, j} = D + \sum_{k = D}^{n} S_{i, j, k}

(4)

This method for calculating the SOD and the SED does not aim to find a specific date for the start and end of the snow cover season. Indeed, it does not consider transient snowfall and melting: thus, for low-elevation regions, the SOD and SED reflect the number of snow cover days in the early part of the hydrological year (SOD, before the 25th of January) and late part of the hydrological year (SED, between the 25th of January and 30th of September). For high-elevation regions, where snow cover is reasonably continuous, the SOD and SED tend to correspond to the effective date of the start and the end of the snow cover season. For these regions, the correspondence between the day of the hydrological year and the date is reported in Table S1.

The trends over the period 2001–2024 for the three metrics were further analyzed. To calculate the slope and intercept of these trends, the Sen–Theil estimator was employed [72], known for its robustness against outliers compared to standard ordinary least-squares linear regression. The significance of the trend was evaluated with the nonparametric Mann–Kendall test [73] that identifies monotonic trends (not necessarily linear) in a time series.

3. Results

3.1. Snow Cover Metrics over the Study Area

The average SCD in the study area during the investigated period is 46 ± 1 days. This value was obtained by calculating the mean SCD value over all the grid points for each year and then averaging these values over the 24-year period. The uncertainty is expressed as the error of the 24-year mean and gives an indication of the robustness of the mean SCD. Figure 4a shows the spatial distribution of the average SCD across the entire area. At small spatial scales, the SCD gradients accurately reflect the orography of the region, with higher values at higher elevations and lower values at lower elevations, highlighting the three major mountain ranges (Karakoram, Himalaya, and Kunlun) and their valleys. However, on a larger spatial scale, there is also a strong influence of the different climatic regions present in the study area.

Starting from the south, the Upper Indus Basin is discernible with SCD values equal to about 9 ± 4 days for points situated below 2000 m a.s.l., terminating on the southern side of the Himalayas. In the Himalayan subregions, SCD reaches 97 ± 3 days at 3000–3500 m a.s.l., 129 ± 3 days at 4000–4500 m a.s.l., and 189 ± 4 days at 5000–5500 m a.s.l. At elevations above 6000 m a.s.l., the average SCD is 287 ± 3 days. Moving northwest, the Karakoram range is characterized by similar SCD values even though with greater changes with elevation: in these subregions, the SCD reaches 63 ± 2 days at 3000–3500 m a.s.l., 139 ± 3 days at 4000–4500 m a.s.l., 194 ± 3 days at 5000–5500 m a.s.l., and 348 ± 1 days above 6000 m a.s.l. On the northern side, the Kunlun Mountains are distinguishable with lower SCD values for the same elevation bands compared to the other two main chains of the considered area, reaching 19 ± 2 days at 3000–3500 m a.s.l., 37 ± 2 days at 4000–4500 m a.s.l., and 70 ± 2 days at 5000–5500 m a.s.l. Higher SCD values in these subregions are found at elevations above 6000 m a.s.l., reaching 340 ± 1 days. In the northern subregions (AT, ET, CK, and WT), the SCD drops to 10 ± 2 days for elevations lower than 1500 m a.s.l. Lastly, on the eastern side of the area, the Tibetan Plateau has lower SCD values despite altitudes exceeding 4000 m a.s.l.; for example, in SWTP at 5000–5500 m a.s.l., the SCD is 30 ± 4 days.

In agreement with this description, when the mean values of the individual subregions are examined (Table 2), those with the highest SCD mean values (EK, NK, WK, and NH) can be distinguished from those with the lowest values (AT, ET, CK, and WT). EK has the highest SCD value of 184 ± 3 days, while ET and WT have the lowest SCD value, standing at 14 ± 1 and 14 ± 2 days, respectively. Indeed, most part of these two subregions is covered by the Taklamakan Desert. However, the SCD is not zero due to the presence of the Kunlun Mountains in their southern part.

The maps of SOD (Figure 4b) and SED (Figure 4c) correspond well with the map of SCD. The average values of SOD and SED for each subregion are displayed in Table 2. The average SOD for all grid points is day 100 ± 1 of the hydrological year, corresponding to the 7th of January. ET and WT exhibit the latest average SOD on day 110 ± 1 (18th of January), while EK has the earliest average SOD on day 53 ± 2 (21st of November). The average SED for all grid points is day 146 ± 1 of the hydrological year, corresponding to the 23rd of February. ET has the earliest average SED on day 125 ± 1 (3rd of February), and EK the latest average SED on day 236 ± 2 (25th of May).

3.2. Spatial Variability in the Snow Cover

The study area exhibits a remarkable diversity of climates, which complicates the identification of the factors influencing snow cover metrics. Figure 5 shows a complex relationship between elevation and SCD; many different SCD values exist at the same elevation. Specifically, for elevations between approximately 4000 and 5000 m a.s.l., SCD values vary significantly, ranging from 0 to 365 days. This variability underscores the importance of considering additional factors beyond elevation, such as regional climatic differences, when analyzing snow cover patterns across the area. Figure 6 illustrates the pixel-by-pixel difference between SCD values and the 5 m elevation band average. Classifying the points in 5 m elevation intervals allows for an accurate description of elevation-driven variations in the snow metrics across the region. At the same time, this approach ensures that each elevation interval is populated with a sufficient number of grid cells, with more than 1000 for any of the considered bins up to over 6000 m a.s.l. As observed from Figure 5 and Figure 6, the analysis reveals significant variability in SCD values at the same elevation. The three main chains of the area present higher SCD values than the overall average. Conversely, the Tibetan Plateau shows lower SCD values. Indeed, about 70% of the grid points above 2000 m a.s.l. in the WK, IH, GH, NH, and EK subregions have higher SCD values (+100 days on average). In contrast, the remaining 30% have lower SCD values (−36 days on average).

At the same time, 90% of the grid points above 2000 m a.s.l. in ST, NT, SWTP, NWTP, WT, CK, ET, AT, and NK subregions have lower SCD values (about −50 days), while only 10% have higher SCD values, most of them in the Kunlun chain.

The spread of the SCD values below 2000 m a.s.l. is more limited, with about 86% of the grid points showing a difference between about −10 and +10 days. This is due to the fact that the snow cover period is generally rather short throughout this area.

To better highlight these findings, Figure 7 illustrates the quantiles of the SCD distribution across the 14 subregions. Each figure depicts seven points per 5 m elevation interval, representing the 1st, 5th, 25th, 50th, 75th, 95th, and 99th percentiles of the SCD distributions of the considered elevation interval. Comparing the behavior of the 14 subregions is easier when focusing on the medians, i.e., the 50th percentiles of the distributions, which are displayed together in Figure 8. These two figures reveal distinct behaviors among the regions. However, they also show that some of the individual subregions exhibit complex patterns. This highlights that the influence of different climates can lead to a non-monotonic relationship between snow cover and elevation, even at the scale of the selected subregions. This is, e.g., evident in IH, where SCD values increase with elevation up to approximately 4000 m a.s.l., above which there is a sharp decline. GH exhibits a similar elevation-dependent increase in SCD, even with a lower slope and no abrupt decrease. This suggests that SCD variability in the Indus and Ganges basins is influenced by factors beyond elevation alone. Above 4500 m a.s.l., both GH and IH SCD values continue to rise at a similar rate; however, GH values exceed those of IH, further confirming the role of other influencing variables.

As one moves from the southern to the northern face of the Himalayan chain, NH initially exhibits a rate of SCD increases similar to GH for elevations below 4000 m a.s.l. Beyond this point, NH exhibits a sharp decline in SCD before increasing again at the same rate as EK, its neighboring region. At approximately 5000 m a.s.l., the SCD in EK and NH increases at a similar rate. It is noteworthy that the SCD curves for elevations above 5500 m a.s.l. become noisier due to the reduced number of available grid points. At lower elevations, EK displays a different pattern compared to the previously discussed regions. SCD increases gradually until about 4000 m a.s.l., followed by a steeper rise. Moving further north, EK’s behavior seems to be the typical one: it is characterized by a gradual increase up to an elevation between 4000 and 5000 m a.s.l., followed by a steep rise to 365 days at around 6000 m a.s.l. Other regions, such as NK, AT, ET, CK, SWTP, NWTP, NT, and ST, exhibit this pattern with minor peculiarities. WK and WT also follow this trend but at different rates of increase. WK SCD starts increasing at around 2000 m a.s.l., while WT begins its rise at 4000 m a.s.l. Both reach an SCD of 365 days at approximately 5500 m a.s.l.

A similar complex behavior is observed for the other snow cover metrics, SOD and SED. These metrics also exhibit significant spatial variability and do not always follow a linear relationship with elevation across the different subregions.

3.3. Interannual Variability in Snow Cover Metrics and Trend Analysis

Figure 9 illustrates the records of the deviation of annual SCD, SOD, and SED values from their mean over the whole considered period for the entire study area. The overall variability across the entire area is relatively low, with standard deviations of 5 days for SCD, 3 days for SOD, and 4 days for SED. This limited variability is likely due to the large size of the study area. The Mann–Kendall test does not indicate any significant trends in the time series of these metrics, while the Sen–Theil slope estimator returns values close to zero for all records. The figure also highlights that the SCD is more in agreement with the SED than with the SOD. This is also evident considering the correlation coefficients of SCD and SOD and of SCD and SED, which are −0.77 and 0.81, respectively.

In order to give a clearer picture of the temporal series of the SCD in the different subregions, they were clustered by means of a Principal Component Analysis (Figure 1). The results highlight that four principal components explain up to 92% of the variance of the SCD records over the fourteen areas. The clustered records are represented in Figure 10, together with the average SCD records of the areas they belong to. None of the clustered series nor the single region series display significant trends. It is, however, worth noticing that the most relevant variation concerns the cluster of the northern subregions (Group 2 in Figure 10). Compared to the average over the entire area, the year-to-year variability is more pronounced when examining individual regions. NWTP exhibits the highest variability at 17 days, while ET shows the lowest at 4 days. Despite these regional differences, all temporal series are characterized by substantial year-to-year variability, with no discernible trends observed over the 24-year period.

As observed before, the regional SCD, SOD, and SED temporal series do not show significant trends. Considering separate 1000 m elevation bands, the SCD variability is highest at 3000–6000 m a.s.l. (Figure 11) and only a few significant trends were identified. Specifically, in the 1000–2000 m a.s.l. elevation band, significant decreasing SCD trends were observed in AT and ET, reaching up to −3.8 days per decade (p-value < 0.05). In the 6000–7000 m a.s.l. band, a decreasing trend of −4.0 days per decade is also observed in AT (p-value < 0.1), while SWTP shows a decrease of −7.6 days per decade in the same range (p-value < 0.1). NK reveals significant negative trends, ranging from −0.4 to −2.5 days per decade, for the 1000–3000 m a.s.l. and 6000–8000 m a.s.l. elevation bands. Significant positive trends (p-value < 0.05) for the SOD were observed in AT at 1000–5000 m a.s.l., ranging between +2.7 and +5.1 days per decade, and in the 2000–3000 m a.s.l. band in NK, with a trend of +1.0 days per decade. For SED, only a few trends reached a p-value lower than 0.1: at 6000–7000 m a.s.l. in AT of −4.0 days per decade and NK of −1.9 days per decade, both with p-0.05.

SCD trends and their statistical significance were computed also for each grid point (Figure 12). Half (53%) of the grid points exhibit no trend, mostly where the SCD remains constant at zero throughout the study period. The remaining half is divided between positive (22%) and negative (26%) trends. Among these, significant positive trends are found for 1% of grid points, significant negative trends in 3%, and the remaining 96% show no significant trend. In the Western Himalayas, many grid points exhibit significant positive trends, between +5 and +25 days per decade. The Karakoram range shows a mix of significant positive and negative trends. Specifically, the northern part of the Karakoram range generally exhibits increasing trends, between +5 and +30 days per decade, while the southern side displays an overall decreasing trend of lower magnitude. Consistent significant negative trends are observed along the borders of the Taklamakan Desert, in the Kunlun Mountains, and in Altyn Tagh, particularly on its eastern side, ranging from −2 days to −40 days per decade. The Northern Tibetan Plateau has predominantly negative trends in reliefs and positive trends in valleys. Around some lakes, highly positive trends (more than +25 days per decade) suggest that, even though the water-covered pixels were carefully masked, possible false snow detection due to the presence of water can be present for a few pixels.

4. Discussion

4.1. Spatial Variability

The analysis of the spatial distributions of SCD, SOD, and SED across the entire study area (see Section 3.1 and Section 3.2) reveals the complexity of the snow cover signal. This makes it challenging to identify the factors influencing the snow cover metrics. In particular, it is difficult to explain the variability of the metrics when only elevation is taken into account, as only a small portion of the variability can be explained by this factor alone.

To analyze this aspect in depth, the possible role of elevation and other geographic and topographic potential drivers (e.g., aspect, slope, latitude, and longitude) on the spatial pattern of the SCD was investigated by means of linear regressions (see Table S2). When all variables are considered, the explained variance reaches 48%. For individual variables, the explained variance was 34% for elevation and 27% for slope. The explained variance generally increases when the subregions are considered individually. Indeed, the multilinear regression in single subregions shows generally an explained variance higher than 50%, reaching more than 60% in nine of them. When considering the individual variables, elevation typically explains the highest fraction of the variance of the SCD spatial signal, with values over 50% in GH, IH, WK, WT, and NT. Conversely, in CK and NH, where there is no clear increase in SCD quantiles with elevation (Figure 7), the variance explained by this variable is much lower (less than 33%). Besides elevation, in regions like IH, AT, ET, WT, and NWTP, slope emerges as the second most important factor, explaining between 13% and 37% of the variance. This highlights the role of topographic factors in driving SCD variability. Longitude and/or latitude emerge as a significant factor in GH, NH, WK, AT, ET, WT, CK, SWTP, and NT, explaining 11–23% of the variability. This is likely due to the presence of areas with different climates within these sub-regions. Longitude and latitude can act as proxies for meteorological variables. This suggests that a division of regions that also considers climatic factors could improve the description of SCD variability within individual regions.

This is particularly evident in a region such as IH, which includes both the southern slopes of the Himalayas on its western side and the northern slopes on its eastern side. The Himalayas act as a barrier for air masses, causing distinct climates in these two areas. The southern slopes are influenced by moist air masses from the south, resulting in higher precipitation levels. For this reason, the mean elevation of the northern slopes of the Himalayas included in IH is 4500 m a.s.l., Figure 7 and Figure 8 show a sharp decline in SCD after 4000 m a.s.l.

SOD and SED show similar pictures when their spatial signal is compared to elevation and to the other geographic and topographic potential drivers.

4.2. Temporal Variability and Meteorological Variables

To gain deeper insights into the causes underlying the SCD trends (Figure 12), snow cover extent (SCE) and temperature and precipitation records from the ERA5 and ERA5-Land reanalyses over the study area were also considered. Reanalyses are the only possible source for meteorological data when ground stations are not available. In other areas of High Mountain Asia with available station observations, temperature has been identified as the dominant factor for changes in snow cover phenology, though its influence varies significantly across different locations [74].

Previous studies highlight that ERA5 and ERA5-Land snow cover records generally overestimate snow on the ground [75]. Moreover, they give evidence of possible inhomogeneities in the reanalysis records. This problem seems to be particularly evident for ERA5, which assimilates snow cover data and is therefore influenced by the temporal changes in the availability of the assimilated observations. ERA5-Land, on the contrary, does not directly assimilate snow cover data, and it is, therefore, more stable over time [76]. However, this reanalysis can also be affected by this problem due to ERA5 atmospheric forcing. In any case, it is important to highlight that the IMS (interactive multi-sensor snow and ice mapping system) fields, which produce a sharp inhomogeneity in the ERA5 SCE record in 2004, are assimilated only for points with elevation lower than 1500 m a.s.l. and thus only for a rather small fraction of the mountain ranges considered in this study. Therefore, the effect of such data on the ERA5-Land SCE records in the study area may not be so relevant. A good performance of the ERA5 and ERA5-Land snow cover records over a wide region, including almost all the areas considered in this paper, is also reported by Yan et al. [77] for the 2001–2017 period.

The temporal series of ERA5-Land SCE were compared with the MODIS SCD temporal series (Figure 13) for the four PCA groups defined in Section 3.3 (see Figure 1). The comparison shows a good agreement, especially in areas with higher snow coverage. The correlation coefficients between the SCD and SCE temporal series range from 0.73 to 0.94. This strong agreement, along with previous observations about the minimal effect of well-known issues with the reanalysis of snow cover data, makes the ERA5-Land SCE records valuable for studying snow cover duration trends over a longer period than the MODIS era alone. The trends in the ERA5-Land subregional SCE records from 1951 were therefore estimated. The results show that, in 10 out of 14 subregions, the changes are not significant. However, four subregions (GH, IH, NH, and WK) exhibit significant negative trends at the 90% confidence level. These results should be interpreted with caution, as there may be significant inhomogeneities in the SCE reanalysis records. Nevertheless, considering that the literature often reports an overestimation of the ERA5-Land SCE negative trends, these findings provide additional evidence that the investigated area is experiencing an anomalous trend in snow cover. This trend likely extends beyond the MODIS era and could span a much longer period.

To strengthen this evidence and to better highlight the anomaly of the investigated area relative to global trends, the ERA5-Land SCE trends of the 25–45°N latitudinal band were also compared with that of the study region (Figure 14). This approach allows for a comparison between series potentially affected by similar biases and instabilities, as they originate from the same dataset. For the period 1951–2024, the negative trend in the 25–45°N band is significant, with a decrease of −0.34% per decade (p < 0.05). In contrast, the trend for the study region is much less evident (−0.23% per decade), and it is not significant. During the MODIS era (2001–2024), the Sen’s slope for the 25–45°N band is −0.70% per decade (p < 0.05). In the same period, the trend for the study area is −0.48% per decade, which is more pronounced than in the 1951–2024 period but still remains non-significant.

The spatial distribution of the trend of the ERA5-Land SCE data in the 2001–2024 period (Figure 15) was then analyzed. The Western Kunlun Mountains exhibited the most substantial loss of SCE, with some slopes showing a decrease of 10–20% in SCE per decade. This loss of SCE is also present in the surrounding areas of the range in the Taklamakan Desert, with lower magnitudes. In the Karakoram range, some points exhibit a decrease in SCE of less than −10% per decade. Despite the difference in resolution between MODIS (approximately 500 m) and ERA5-Land (about 9 km), the overall patterns in the SCE trend closely match those observed in the SCD trend map, particularly in the northern regions.

The same analysis for ERA5-Land temperature data (Figure 16) and for ERA5 precipitation data (Figure 17) was then performed.

As far as temperature is concerned, in the north of the Kunlun range, it has increased by approximately 0.1–1 °C per decade in the 2001–2024 period. This rise is statistically significant across a broad region that includes the Western Taklamakan Desert, the northern slopes of the entire Kunlun Range, and the Altyn Tagh. The Westernmost Taklamakan Desert and Northern Altyn Tagh also experienced, for some grid points, a reduction in total annual precipitation ranging between 15 and 30 mm per decade. In particular, when dividing the hydrological year into two parts, it was observed that the precipitation decrease occurred primarily between October and March.

In agreement with the rise in temperature and decrease in precipitation, 5% of the grid points in WT, CK, ET, and AT show a significant negative SCD trend, with values ranging from −3 to −25 days per decade. In contrast, only 0.7% of the grid points show a significant positive trend. However, the number of grid points showing a statistically significant trend is not sufficient to establish a regional trend. These observations reflect the findings on SCE.

In the southern part of the study area, very few grid points show a statistically significant temperature trend, either positive or negative. Regarding the precipitation trends, the Karakoram region does not show statistically significant trends. However, the Himalayan range gives evidence of a statistically significant increase of approximately +200 mm per decade, occurring mainly between April and September. In agreement with the increase in precipitation and a non-significant temperature variation, SCD trends show significant positive values in 1.3% of grid points in the IH region, with an increase of about 12 days per decade. In contrast, 0.4% of the points show a significant negative trend. This is despite the precipitation increase that is happening during the warmer months. Indeed, this tendency mainly concerns points in the Himalayan range where precipitation falls as snow above 5000 m a.s.l. even during summer [24].

4.3. Comparison with Previous Studies

The trend results are in good agreement with the existing literature. On a larger scale, Ackroyd et al. [39] and Dharpure et al. [40] found no significant trend in snow metrics in the study area. When comparing the temporal series of the SCD in the Himalayan regions (First Principal Component: IH, GH, NH, ST, NT) with the WH basin analyzed by Dharpure et al., the analysis of extreme years also shows similar findings. Specifically, there was a shorter than average SCD during the winters of 2015/2016 and 2017/2018, and a longer than average SCD in 2004/2005 and 2014/2015. At a grid point level, the results on SCD trends fit also well with the results of Tang et al. [41].

The agreement of the results of this work with the existing literature is rather good also on a smaller spatial scale, where the trends are either not statistically significant or show contrasting situations. For example, Yi et al. [38] found that snow metrics in the Western Upper Indus and Yarkand River basins showed no significant trends from 2002 to 2018. Regarding snow cover extent in the Karakoram and Western Himalayas, the literature reveals similar findings. Choudhury et al. [37] and Notarnicola [31] found trends that are either not clear and robust or slightly positive. For parts of these areas, Atif et al. [32], Dharpure et al. [40], Hussain et al. [33], and Moazzam et al. [34] showed that the results are often not significant and strongly depend on the specific basin and timeframe considered. Conversely, in the Tarim Basin, Zhang et al. [36] showed that the trends in snow cover extent exhibit a clearer sign of decline, which agrees with what was observed in this work. Indeed, in the Kunlun Mountains, a significant number of grid points exhibit a decreasing SCD.

The analysis of snow cover area percentage from monthly mean time series over 1985–2011 indicated no significant linear trend for the entire alpine region, too. The same is true for regional analysis, even though some tendencies toward increasing/decreasing were found by Husler et al. [78]. For the Alpine region, however, Fugazza et al. found significant trends in snow cover metrics for points above 3000 m a.s.l that reported significant negative trends [68]. Differently, snow cover extent decreased in the Central Andes of Chile and Argentina by about 13% between 2000 and 2016. This decline occurred throughout the year, with the largest decrease seen at the onset and end of the snow season [79].

5. Conclusions

In this study, the MODIS snow cover product was used to describe the snow cover across Karakoram, Western Himalayas, and Kunlun Mountains, dividing the area into fourteen subregions with a basin-based approach. To delineate the spatial and temporal extent of annual snow cover, three metrics (SCD, SOD, SED) were computed for each grid point and for the subregions.

Average SCD in the investigated period varied significantly over the study region: elevation was the key factor, but strong variability at fixed elevation levels highlighted a rather complex pattern. This gives evidence that other geographic and climatic factors have to be considered to capture the spatial distribution of snow cover.

The SCD subregional average temporal records were clustered using Principal Component Analysis, revealing that four components explain up to 92% of their variance. In general, trend analysis using the Sen–Theil and Mann–Kendall tests did not indicate statistically significant changes in the SCD across the study area. However, significant trends were identified in specific regions and elevation bands, as in the Eastern Taklamakan Desert (1000–2000 m a.s.l. and 6000–7000 m a.s.l. elevation bands) in the Southwest Tibetan Plateau (6000–7000 m a.s.l. elevation band) and in Northern Karakoram (1000–3000 m a.s.l. and 6000–8000 m a.s.l. elevation bands). The SCD regional temporal series remained, however, generally stable across the study area from 2001 to 2024. This distinguishes the investigated region from others in the world, where the impact of climate change has more frequently led to decreasing SCD trends and affected other snow parameters.

The average SCD records over the areas identified by the Principal Component Analysis were also compared with corresponding ERA5-Land snow cover extent records, identifying correlation coefficients from 0.78 to 0.95, with the highest values for the areas with longer snow seasons. The good agreement between MODIS-derived SCD and ERA5-Land snow cover extent temporal series allows for the use of reanalysis data to better interpret the drivers behind the observed trends in snow cover season metrics. The reanalysis data were therefore used to better understand the changes in the snow cover metrics, highlighting a good coherence between the trends of the snow cover metrics and the corresponding trends of temperature and precipitation.

Globally, the results give evidence that the MODIS data updated to 2024 confirm that the investigated area generally highlights low trends in the snow cover metrics, providing another piece of evidence on its anomalous behavior with respect to most areas of the Planet, which shows a clear decrease in snow cover over the last decades. Reanalysis of snow cover extend data seems to confirm that this anomalous behavior concerns a much longer period than the MODIS era only.

The relatively short timeframe of this study and the prevalence of year-to-year variability suggest a complex interplay of factors influencing climate dynamics in these regions. Exploring longer satellite datasets, further investigating the problems of reanalysis snow cover data, and validating both satellite and reanalysis data with observational records could enhance the understanding of the snow cover season in the area. Future research could also use climate models to project climate change’s impact on snow extent and regional climate dynamics, improving future water management strategies.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/rs17050914/s1, Table S1: Correspondence between the calendar date and the day of the hydrological year (starting 1st of October and ending 30th of September). Table S2: Variance explained by multilinear regression and linear regression of individual variables against the three metrics in the fourteen subregions and over the entire study area. Each cell reports the R2. The acronyms of the subregions follow the definition of Table 1.

Author Contributions

Conceptualization, C.D.A., V.M., M.M. and D.F.; methodology, C.D.A., V.M., M.M. and D.F.; software, C.D.A. and A.S.; validation, C.D.A.; formal analysis, C.D.A. and A.S.; investigation, C.D.A.; writing—original draft preparation, C.D.A., V.M., M.M. and D.F.; writing—review and editing, C.D.A., V.M., M.M., G.A.D., A.S. and D.F.; visualization, C.D.A., V.M., M.M. and D.F.; supervision, V.M., M.M. and D.F.; project administration, M.M. and G.A.D.; funding acquisition, M.M. and G.A.D. All authors have read and agreed to the published version of the manuscript.

Funding

This study was funded within “Glaciers and Students” project (project number: 00144462) and “Water for Development (W4D)” project (Award ID: 1274884), funded by the Ministry of Foreign Affairs and International Cooperation and the Italian Agency for Development Cooperation (AICS), executed by the United Nations Development Programme (UNDP) and implemented by EvK2CNR.

Data Availability Statement

The results of the analysis are available upon request to the authors of the paper. The MOD10A1 V61 and the MYD10A1 V61 daily snow cover products are freely available from https://earthdata.nasa.gov/ (accessed on 10 December 2024). ERA5 and ERA5 land data are openly available from Copernicus Climate Change Service (C3S) in the Climate Data Store (CDS), respectively: “ERA5 monthly averaged data on single levels from 1940 to present” at https://cds.climate.copernicus.eu/datasets/reanalysis-era5-single-levels-monthly-means (last accessed on 10 December 2024) and “ERA5-Land monthly averaged data from 1950 to present” at https://cds.climate.copernicus.eu/datasets/reanalysis-era5-land-monthly-means (accessed on 10 December 2024). Diva-Gis water masks are available from https://diva-gis.org/ (last accessed on 10 November 2023). ALOS AW3D30 Digital Elevation Models are available after registration from https://www.eorc.jaxa.jp/ALOS/en/aw3d30/data/index.htm (accessed on 10 September 2023).

Acknowledgments

This study was performed in the framework of the “Glaciers and Students” and “Water for Development (W4D)” projects, funded by the Ministry of Foreign Affairs and International Cooperation and the Italian Agency for Development Cooperation (AICS), executed by the United Nations Development Programme (UNDP) and implemented by EvK2CNR. Veronica Manara was supported by the “Ministero dell’Università e della Ricerca” of Italy [grant FSE—REACT EU, DM 10/08/2021 n. 1062].

Conflicts of Interest

The authors declare no conflicts of interest.

References

Intergovernmental Panel on Climate Change (IPCC). Climate Change 2021—The Physical Science Basis; Cambridge University Press: Cambridge, UK, 2023; ISBN 9781009157896.
Kunkel, K.E.; Robinson, D.A.; Champion, S.; Yin, X.; Estilow, T.; Frankson, R.M. Trends and Extremes in Northern Hemisphere Snow Characteristics. Curr. Clim. Change Rep. 2016, 2, 65–73. [Google Scholar] [CrossRef]
Lemke, P.; Ren, J.; Alley, R.B.; Allison, I.; Carrasco, J.; Flato, G.; Fujii, Y.; Kaser, G.; Mote, P.; Thomas, R.H.; et al. Observations: Changes in Snow, Ice and Frozen Ground. In Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change; Solomon, S.D., Qin, M., Manning, Z., Chen, M., Marquis, K.B., Averyt, M., Tignor, H.L., Miller, H.L., Eds.; Cambridge University Press: Cambridge, UK, 2007. [Google Scholar]
Vaughan, D.G.; Comiso, J.; Allison, I.; Carrasco, J.; Kaser, G.; Kwok, R.; Mote, P.; Murray, T.; Paul, F.; Ren, J.; et al. Observations: Cryosphere. In Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change; Stocker, T.F., Qin, D., Plattner, G.-K., Tignor, M., Allen, S.K., Boschung, J., Nauels, A., Xia, Y., Bex, V., Midgley, P.M., Eds.; Cambridge University Press: Cambridge, UK, 2013. [Google Scholar]
Zemp, M.; Roer, I.; Kääb, A.; Hoelzle, M.; Paul, F.; Haeberli, W. Global Glacier Changes: Facts and Figures. In World Glacier Monitoring Service; UNEP—UN Environment Programme: Nairobi, Kenya, 2008. [Google Scholar]
Flanner, M.G.; Shell, K.M.; Barlage, M.; Perovich, D.K.; Tschudi, M.A. Radiative Forcing and Albedo Feedback from the Northern Hemisphere Cryosphere between 1979 and 2008. Nat. Geosci. 2011, 4, 151–155. [Google Scholar] [CrossRef]
Barnett, T.P.; Dümenil, L.; Schlese, U.; Roeckner, E.; Latif, M. The Effect of Eurasian Snow Cover on Regional and Global Climate Variations. J. Atmos. Sci. 1989, 46, 661–686. [Google Scholar] [CrossRef]
Huss, M.; Bookhagen, B.; Huggel, C.; Jacobsen, D.; Bradley, R.S.; Clague, J.J.; Vuille, M.; Buytaert, W.; Cayan, D.R.; Greenwood, G.; et al. Toward Mountains without Permanent Snow and Ice. Earths Future 2017, 5, 418–435. [Google Scholar] [CrossRef]
Bormann, K.J.; Brown, R.D.; Derksen, C.; Painter, T.H. Estimating Snow-Cover Trends from Space. Nat. Clim. Change 2018, 8, 924–928. [Google Scholar] [CrossRef]
Notarnicola, C. Hotspots of Snow Cover Changes in Global Mountain Regions over 2000–2018. Remote Sens Environ. 2020, 243, 111781. [Google Scholar] [CrossRef]
Matiu, M.; Crespi, A.; Bertoldi, G.; Carmagnola, C.M.; Marty, C.; Morin, S.; Schöner, W.; Cat Berro, D.; Chiogna, G.; De Gregorio, L.; et al. Observed Snow Depth Trends in the European Alps: 1971 to 2019. Cryosphere 2021, 15, 1343–1382. [Google Scholar] [CrossRef]
Salzmann, N.; Huggel, C.; Rohrer, M.; Stoffel, M. Data and Knowledge Gaps in Glacier, Snow and Related Runoff Research—A Climate Change Adaptation Perspective. J. Hydrol. 2014, 518, 225–234. [Google Scholar] [CrossRef]
National Research Council. Division on Earth and Life Studies; Board on Atmospheric Sciences and Climate; Committee on Climate Data Records from NOAA Operational Satellites. In Climate Data Records from Environmental Satellites; National Academies Press: Washington, DC, USA, 2004.
Harrison, A.R.; Lucas, R.M. Multi-Spectral Classification of Snow Using NOAA AVHRR Imagery. Int. J. Remote Sens. 1989, 10, 907–916. [Google Scholar] [CrossRef]
Dietz, A.; Conrad, C.; Kuenzer, C.; Gesell, G.; Dech, S. Identifying Changing Snow Cover Characteristics in Central Asia between 1986 and 2014 from Remote Sensing Data. Remote Sens. 2014, 6, 12752–12775. [Google Scholar] [CrossRef]
Immerzeel, W.W.; Droogers, P.; de Jong, S.M.; Bierkens, M.F.P. Large-Scale Monitoring of Snow Cover and Runoff Simulation in Himalayan River Basins Using Remote Sensing. Remote Sens. Environ. 2009, 113, 40–49. [Google Scholar] [CrossRef]
Tekeli, A.E.; Akyürek, Z.; Arda Şorman, A.; Şensoy, A.; Ünal Şorman, A. Using MODIS Snow Cover Maps in Modeling Snowmelt Runoff Process in the Eastern Part of Turkey. Remote Sens. Environ. 2005, 97, 216–230. [Google Scholar] [CrossRef]
Hall, D.K.; Riggs, G.A.; Solomonson, V. NASA MODAPS SIPS MODIS/Terra Snow Cover Daily L3 Global 500m SIN Grid. Available online: https://data.nasa.gov/dataset/MODIS-Terra-Snow-Cover-Daily-L3-Global-500m-SIN-Gr/ih6z-662z/about_data (accessed on 27 February 2025).
Hall, D.K.; Riggs, G.A. Accuracy Assessment of the MODIS Snow Products. Hydrol. Process. 2007, 21, 1534–1547. [Google Scholar] [CrossRef]
Coll, J.; Li, X. Comprehensive Accuracy Assessment of MODIS Daily Snow Cover Products and Gap Filling Methods. ISPRS J. Photogramm. Remote Sens. 2018, 144, 435–452. [Google Scholar] [CrossRef]
Parajka, J.; Holko, L.; Kostka, Z.; Blöschl, G. MODIS Snow Cover Mapping Accuracy in a Small Mountain Catchment—Comparison between Open and Forest Sites. Hydrol. Earth Syst. Sci. 2012, 16, 2365–2377. [Google Scholar] [CrossRef]
Notarnicola, C.; Duguay, M.; Moelg, N.; Schellenberger, T.; Tetzlaff, A.; Monsorno, R.; Costa, A.; Steurer, C.; Zebisch, M. Snow Cover Maps from MODIS Images at 250 m Resolution, Part 2: Validation. Remote Sens. 2013, 5, 1568–1587. [Google Scholar] [CrossRef]
Hewitt, K. The Karakoram Anomaly? Glacier Expansion and the ‘Elevation Effect,’ Karakoram Himalaya. Mt. Res. Dev. 2005, 25, 332–340. [Google Scholar] [CrossRef]
Winiger, M.; Gumpert, M.; Yamout, H. Karakorum–Hindukush–Western Himalaya: Assessing High-altitude Water Resources. Hydrol. Process. 2005, 19, 2329–2338. [Google Scholar] [CrossRef]
Akhtar, M.; Ahmad, N.; Booij, M.J. The Impact of Climate Change on the Water Resources of Hindukush–Karakorum–Himalaya Region under Different Glacier Coverage Scenarios. J. Hydrol. 2008, 355, 148–163. [Google Scholar] [CrossRef]
Kahlown, M.A.; Raoof, A.; Zubair, M.; Kemper, W.D. Water Use Efficiency and Economic Feasibility of Growing Rice and Wheat with Sprinkler Irrigation in the Indus Basin of Pakistan. Agric. Water Manag. 2007, 87, 292–298. [Google Scholar] [CrossRef]
Aggarwal, P.K.; Joshi, P.K.; Ingram, J.S.I.; Gupta, R.K. Adapting Food Systems of the Indo-Gangetic Plains to Global Environmental Change: Key Information Needs to Improve Policy Formulation. Environ. Sci. Policy 2004, 7, 487–498. [Google Scholar] [CrossRef]
Shean, D.E.; Bhushan, S.; Montesano, P.; Rounce, D.R.; Arendt, A.; Osmanoglu, B. A Systematic, Regional Assessment of High Mountain Asia Glacier Mass Balance. Front. Earth Sci. 2020, 7, 363. [Google Scholar] [CrossRef]
Minora, U.; Bocchiola, D.; D’Agata, C.; Maragno, D.; Mayer, C.; Lambrecht, A.; Vuillermoz, E.; Senese, A.; Compostella, C.; Smiraglia, C.; et al. Glacier Area Stability in the Central Karakoram National Park (Pakistan) in 2001–2010. Prog. Phys. Geogr. Earth Environ. 2016, 40, 629–660. [Google Scholar] [CrossRef]
Li, J.; Sun, M.; Yao, X.; Duan, H.; Zhang, C.; Wang, S.; Niu, S.; Yan, X. A Review of Karakoram Glacier Anomalies in High Mountains Asia. Water 2023, 15, 3215. [Google Scholar] [CrossRef]
Notarnicola, C. Observing Snow Cover and Water Resource Changes in the High Mountain Asia Region in Comparison with Global Mountain Trends over 2000–2018. Remote Sens. 2020, 12, 3913. [Google Scholar] [CrossRef]
Atif, I.; Iqbal, J.; Mahboob, M. Investigating Snow Cover and Hydrometeorological Trends in Contrasting Hydrological Regimes of the Upper Indus Basin. Atmosphere 2018, 9, 162. [Google Scholar] [CrossRef]
Hussain, D.; Kuo, C.-Y.; Hameed, A.; Tseng, K.-H.; Jan, B.; Abbas, N.; Kao, H.-C.; Lan, W.-H.; Imani, M. Spaceborne Satellite for Snow Cover and Hydrological Characteristic of the Gilgit River Basin, Hindukush–Karakoram Mountains, Pakistan. Sensors 2019, 19, 531. [Google Scholar] [CrossRef]
Moazzam, M.F.U.; Rahman, G.; Ali, K.S.; Lee, B.G. Spatio-Temporal Snow Cover Change in the Early Twenty-First Century Using Improved MODIS Dataset: A Case Study of District Hunza, Pakistan. Clim. Dyn. 2023, 60, 3417–3433. [Google Scholar] [CrossRef]
Thapa, A.; Muhammad, S. Contemporary Snow Changes in the Karakoram Region Attributed to Improved MODIS Data between 2003 and 2018. Water 2020, 12, 2681. [Google Scholar] [CrossRef]
Zhang, J.; Jia, L.; Menenti, M.; Zhou, J.; Ren, S. Glacier Area and Snow Cover Changes in the Range System Surrounding Tarim from 2000 to 2020 Using Google Earth Engine. Remote Sens. 2021, 13, 5117. [Google Scholar] [CrossRef]
Choudhury, A.; Yadav, A.C.; Bonafoni, S. A Response of Snow Cover to the Climate in the Northwest Himalaya (NWH) Using Satellite Products. Remote Sens. 2021, 13, 655. [Google Scholar] [CrossRef]
Yi, Y.; Liu, S.; Zhu, Y.; Wu, K.; Xie, F.; Saifullah, M. Spatiotemporal Heterogeneity of Snow Cover in the Central and Western Karakoram Mountains Based on a Refined MODIS Product during 2002–2018. Atmos. Res. 2021, 250, 105402. [Google Scholar] [CrossRef]
Ackroyd, C.; Skiles, S.M.; Rittger, K.; Meyer, J. Trends in Snow Cover Duration Across River Basins in High Mountain Asia From Daily Gap-Filled MODIS Fractional Snow Covered Area. Front. Earth Sci. 2021, 9, 713145. [Google Scholar] [CrossRef]
Dharpure, J.K.; Goswami, A.; Patel, A.; Kulkarni, A.V. Snehmani Assessment of Snow Cover Variability and Its Sensitivity to Hydrometeorological Factors in the Karakoram and Himalayan Region. Hydrol. Sci. J. 2021, 66, 2198–2215. [Google Scholar] [CrossRef]
Tang, Z.; Deng, G.; Hu, G.; Zhang, H.; Pan, H.; Sang, G. Satellite Observed Spatiotemporal Variability of Snow Cover and Snow Phenology over High Mountain Asia from 2002 to 2021. J. Hydrol. 2022, 613, 128438. [Google Scholar] [CrossRef]
Xu, J.; Tang, Y.; Dong, L.; Wang, S.; Yu, B.; Wu, J.; Zheng, Z.; Huang, Y. Temperature-Dominated Spatiotemporal Variability in Snow Phenology on the Tibetan Plateau from 2002 to 2022. Cryosphere 2024, 18, 1817–1834. [Google Scholar] [CrossRef]
Wang, J.; Tang, Z.; Deng, G.; Hu, G.; You, Y.; Zhao, Y. Landsat Satellites Observed Dynamics of Snowline Altitude at the End of the Melting Season, Himalayas, 1991–2022. Remote Sens. 2023, 15, 2534. [Google Scholar] [CrossRef]
Tang, Z.; Wang, X.; Deng, G.; Wang, X.; Jiang, Z.; Sang, G. Spatiotemporal Variation of Snowline Altitude at the End of Melting Season across High Mountain Asia, Using MODIS Snow Cover Product. Adv. Space Res. 2020, 66, 2629–2645. [Google Scholar] [CrossRef]
Hersbach, H.; Bell, B.; Berrisford, P.; Hirahara, S.; Horányi, A.; Muñoz-Sabater, J.; Nicolas, J.; Peubey, C.; Radu, R.; Schepers, D.; et al. The ERA5 Global Reanalysis. Q. J. R. Meteorol. Soc. 2020, 146, 1999–2049. [Google Scholar] [CrossRef]
Muñoz-Sabater, J.; ERA5-Land Hourly Data from 1950 to Present. Copernicus Climate Change Service (C3S) Climate Data Store (CDS). 2019. Available online: https://cds.climate.copernicus.eu/datasets/reanalysis-era5-land?tab=overview (accessed on 27 February 2025).
Hewitt, K. Glaciers of the Karakoram Himalaya. In Encyclopedia of Snow, Ice and Glaciers; Singh, V.P., Singh, P., Haritashya, U.K., Eds.; Springer: Berlin, Germany, 2011; pp. 429–436. [Google Scholar]
Young, G.; Hewitt, K. Hydrology Research in the Upper Indus Basin, Karakora, Himalaya, Pakistan. Hydrol. Mt. Area 1990, 190, 139–152. [Google Scholar]
Mayewski, P.A.; Jeschke, P.A. Himalayan and Trans-Himalayan Glacier Fluctuations Since AD 1812. Arct. Alp. Res. 1979, 11, 267–287. [Google Scholar] [CrossRef]
Hewitt, K. Glacier Change, Concentration, and Elevation Effects in the Karakoram Himalaya, Upper Indus Basin. Mt. Res. Dev. 2011, 31, 188–200. [Google Scholar] [CrossRef]
Chatterjee, S.P.; Bishop, B.C. Himalayas; Encyclopedia Britannica: Chicago, IL, USA, 2021. [Google Scholar]
Rees, H.G.; Collins, D.N. Regional Differences in Response of Flow in Glacier-fed Himalayan Rivers to Climatic Warming. Hydrol. Process. 2006, 20, 2157–2169. [Google Scholar] [CrossRef]
Wylie, T.V.; Shakabpa, T.W.D.; Falkenheim, V.C.; Richardson, H.E. Tibet; Encyclopedia Britannica: Chicago, IL, USA, 2025. [Google Scholar]
Wang, X.; Pang, G.; Yang, M. Precipitation over the Tibetan Plateau during Recent Decades: A Review Based on Observations and Simulations. Int. J. Climatol. 2018, 38, 1116–1131. [Google Scholar] [CrossRef]
Xia, Z.; Chen, Y.; Zhang, X.; Li, Z.; Fang, G.; Zhu, C.; Li, Y.; Li, J.; Xia, Q.; Liang, Q. Precipitation Changes on the Northern Slope of the Kunlun Mountains in the Past 42 Years. Water 2024, 16, 1203. [Google Scholar] [CrossRef]
Tao, H.; Gemmer, M.; Bai, Y.; Su, B.; Mao, W. Trends of Streamflow in the Tarim River Basin during the Past 50years: Human Impact or Climate Change? J. Hydrol. 2011, 400, 1–9. [Google Scholar] [CrossRef]
Xu, Z.; Liu, Z.; Fu, G.; Chen, Y. Trends of Major Hydroclimatic Variables in the Tarim River Basin during the Past 50 Years. J. Arid Environ. 2010, 74, 256–267. [Google Scholar] [CrossRef]
Wang, X.; Tolksdorf, V.; Otto, M.; Scherer, D. WRF-based Dynamical Downscaling of ERA5 Reanalysis Data for High Mountain Asia: Towards a New Version of the High Asia Refined Analysis. Int. J. Climatol. 2021, 41, 743–762. [Google Scholar] [CrossRef]
Salomonson, V.V.; Appel, I. Development of the Aqua MODIS NDSI Fractional Snow Cover Algorithm and Validation Results. IEEE Trans. Geosci. Remote Sens. 2006, 44, 1747–1756. [Google Scholar] [CrossRef]
Platnick, S.; King, M.D.; Ackerman, S.A.; Menzel, W.P.; Baum, B.A.; Riedi, J.C.; Frey, R.A. The MODIS Cloud Products: Algorithms and Examples from Terra. IEEE Trans. Geosci. Remote Sens. 2003, 41, 459–473. [Google Scholar] [CrossRef]
Platnick, S.; Meyer, K.G.; King, M.D.; Wind, G.; Amarasinghe, N.; Marchant, B.; Arnold, G.T.; Zhang, Z.; Hubanks, P.A.; Holz, R.E.; et al. The MODIS Cloud Optical and Microphysical Products: Collection 6 Updates and Examples from Terra and Aqua. IEEE Trans. Geosci. Remote Sens. 2017, 55, 502–525. [Google Scholar] [CrossRef] [PubMed]
Tadono, T.; Ishida, H.; Oda, F.; Naito, S.; Minakawa, K.; Iwamoto, H. Precise Global DEM Generation by ALOS PRISM. ISPRS Ann. Photogramm. Remote Sens. Spat. Inf. Sci. 2014, II–4, 71–76. [Google Scholar] [CrossRef]
Tran, H.; Nguyen, P.; Ombadi, M.; Hsu, K.; Sorooshian, S.; Qing, X. A Cloud-Free MODIS Snow Cover Dataset for the Contiguous United States from 2000 to 2017. Sci. Data 2019, 6, 180300. [Google Scholar] [CrossRef]
Maskey, S.; Uhlenbrook, S.; Ojha, S. An Analysis of Snow Cover Changes in the Himalayan Region Using MODIS Snow Products and In-Situ Temperature Data. Clim. Change 2011, 108, 391–400. [Google Scholar] [CrossRef]
Li, X.; Jing, Y.; Shen, H.; Zhang, L. The Recent Developments in Cloud Removal Approaches of MODIS Snow Cover Product. Hydrol. Earth Syst. Sci. 2019, 23, 2401–2416. [Google Scholar] [CrossRef]
Deng, G.; Tang, Z.; Dong, C.; Shao, D.; Wang, X. Development and Evaluation of a Cloud-Gap-Filled MODIS Normalized Difference Snow Index Product over High Mountain Asia. Remote Sens. 2024, 16, 192. [Google Scholar] [CrossRef]
Gafurov, A.; Bárdossy, A. Cloud Removal Methodology from MODIS Snow Cover Product. Hydrol. Earth Syst. Sci. 2009, 13, 1361–1373. [Google Scholar] [CrossRef]
Fugazza, D.; Manara, V.; Senese, A.; Diolaiuti, G.; Maugeri, M. Snow Cover Variability in the Greater Alpine Region in the MODIS Era (2000–2019). Remote Sens. 2021, 13, 2945. [Google Scholar] [CrossRef]
Dietz, A.J.; Wohner, C.; Kuenzer, C. European Snow Cover Characteristics between 2000 and 2011 Derived from Improved MODIS Daily Snow Cover Products. Remote Sens. 2012, 4, 2432–2454. [Google Scholar] [CrossRef]
Parajka, J.; Blöschl, G. Spatio-temporal Combination of MODIS Images—Potential for Snow Cover Mapping. Water Resour. Res. 2008, 44, W03406. [Google Scholar] [CrossRef]
Wang, X.; Xie, H. New Methods for Studying the Spatiotemporal Variation of Snow Cover Based on Combination Products of MODIS Terra and Aqua. J. Hydrol. 2009, 371, 192–200. [Google Scholar] [CrossRef]
Theil, H. A Rank-Invariant Method of Linear and Polynomial Regression Analysis. In Henri Theil’s Contributions to Economics and Econometrics: Econometric Theory and Methodology; Springer: Berlin, Germany, 1992; pp. 345–381. [Google Scholar]
Sneyers, R. Über Den Einsatz von Statistischen Methoden Zum Objektiven Nachweis von Klimaschwankungen. Meteorol. Z. 1992, 1, 247–256. [Google Scholar] [CrossRef]
Ma, X.; Lin, K.; Sun, X.; Luo, L.; Ma, N.; Zha, H.; Zhang, L.; Tang, S.; Tang, Z.; Zhang, H. Investigating Climatic Drivers of Snow Phenology by Considering Key-Substage Heterogeneity. J. Hydrol. 2024, 645, 132215. [Google Scholar] [CrossRef]
Kouki, K.; Luojus, K.; Riihelä, A. Evaluation of Snow Cover Properties in ERA5 and ERA5-Land with Several Satellite-Based Datasets in the Northern Hemisphere in Spring 1982–2018. Cryosphere 2023, 17, 5007–5026. [Google Scholar] [CrossRef]
Urraca, R.; Gobron, N. Temporal Stability of Long-Term Satellite and Reanalysis Products to Monitor Snow Cover Trends. Cryosphere 2023, 17, 1023–1052. [Google Scholar] [CrossRef]
Yan, S.; Chen, Y.; Hou, Y.; Liu, K.; Li, X.; Xing, Y.; Wu, D.; Cui, J.; Zhou, Y.; Pu, W.; et al. Which Global Reanalysis Dataset Has Better Representativeness in Snow Cover on the Tibetan Plateau? Cryosphere 2024, 18, 4089–4109. [Google Scholar] [CrossRef]
Hüsler, F.; Jonas, T.; Riffler, M.; Musial, J.P.; Wunderle, S. A Satellite-Based Snow Cover Climatology (1985–2011) for the European Alps Derived from AVHRR Data. Cryosphere 2014, 8, 73–90. [Google Scholar] [CrossRef]
Malmros, J.K.; Mernild, S.H.; Wilson, R.; Tagesson, T.; Fensholt, R. Snow Cover and Snow Albedo Changes in the Central Andes of Chile and Argentina from Daily MODIS Observations (2000–2016). Remote Sens. Environ. 2018, 209, 240–252. [Google Scholar] [CrossRef]

Figure 1. Study area. The border of the MODIS tile is represented with a bold black line and the fourteen subregions with a thin black line. The subregions are named following the acronyms defined in Table 1, and the four different colors cluster them into four different groups by means of the PCA discussed in Section 3.3. The color scheme of the labels represents these four groups: Group 1 is yellow, Group 2 is green, Group 3 is red, and Group 4 is blue.

Figure 2. Orography of the study area, alongside elevation distribution of the fourteen subregions represented by means of the percentage of pixels for each 500 m elevation band. The acronyms of the subregions are defined in Table 1.

Figure 3. Total precipitation averaged over the 1991–2020 period using ERA5 data for the study area, alongside monthly precipitation distributions for the fourteen subregions represented by means of the percentage of precipitation with respect to the annual total. The acronyms of the subregions are defined in Table 1.

Figure 4. (a) Average snow-covered days (SCD); (b) average snow onset date (SOD) and (c) average snow end date (SED). The three metrics refer to the 2001–2024 period and are expressed in days. The acronyms of the subregions follow the definition of Table 1.

Figure 5. Snow-covered days (SCD) distribution with respect to the elevation of all the grid points (black points) and the mean value for each subregion (colored points). The acronyms of the subregions follow the definition of Table 1.

Figure 6. Difference between the snow-covered days (SCD) of each grid point and the corresponding average SCD value over the whole considered area of points in the same 5 m elevation band. The acronyms of the subregions follow the definition of Table 1.

Figure 7. Quantiles (colored points) of the snow-covered days (SCD) distributions for every 5 m elevation band in the different subregions. The acronyms of the subregions follow the definition of Table 1.

Figure 8. Median (colored points) of the snow-covered days (SCD) distributions for every 5 m elevation band in the different regions. The acronyms of the subregions follow the definition of Table 1.

Figure 9. Additive anomaly series with respect to the whole period of the snow-covered days (SCD—blue line), snow onset date (SOD—orange line), and snow end date (SED—yellow line) over the whole area in the 24-year study period.

Figure 10. Additive anomaly series with respect to the whole period of the snow-covered days (SCD) in the 24-year study period for the fourteen subregions (thin lines) clustered into 4 groups using PCA. The bold red lines represent the average series within each cluster. The acronyms of the subregions follow the definition of Table 1.

Figure 11. Snow-covered days (SCD) series for each subregion and 1000 m elevation band. The acronyms of the subregions follow the definition of Table 1.

Figure 12. Sen’s slope of the significant SCD trends (p-value < 0.1) for each grid point, expressed in days per decade. The acronyms of the subregions follow the definition of Table 1.

Figure 13. Comparison between SCD (MODIS) and SCE (ERA5-Land) 24-year temporal series for the four groups of regions defined with the PCA (see Figure 1).

Figure 14. Comparison between SCE (ERA5-Land) average anomalies in the study area (blue) and in the latitudinal band 25–45°N (orange) over 1951–2024 and 2001–2024.

Figure 15. Sen’s slope of the significant SCE (ERA5-Land) trends (p-value < 0.1) for each grid point, expressed in % per decade. The acronyms of the subregions follow the definition of Table 1.

Figure 16. Sen’s slope of the significant temperature (ERA5-Land) trends (p-value < 0.1) for each grid point, expressed in °C per decade. The acronyms of the subregions follow the definition of Table 1.

Figure 17. Sen’s slope of the significant precipitation (ERA5) trends (p-value < 0.1) for each grid point, expressed in mm per decade. The acronyms of the subregions follow the definition of Table 1.

Table 1. Details on the name, location, and elevation of the fourteen subregions. Information (yearly cumulated value, minimum and maximum yearly value averaged over the 1991–2020 period) about precipitation values obtained with ERA5 is also given.

Name and Acronym	Extension [km²]	Average Elevation [m. a.s.l.]	Elevation Range [m. a.s.l.]	Averaged Yearly Cumulated Precipitation Value [mm]	Month with Min. Precipitation [mm]	Month with Max. Precipitation [mm]
Ganges-Himalayas (GH)	62,143	1974	170–7641	1260	November (16)	July (308)
Indus-Himalayas (IH)	313,434	1550	111–7037	905	October (23)	July (179)
Northern Himalayas (NH)	47,608	4643	2285–6916	573	October (20)	August (68)
Northern Transhimalayas (NT)	35,207	5030	3605–6679	305	November (6)	August (78)
Southern Transhimalayas (ST)	25,087	5068	4220–6445	334	November (4)	August (97)
South Western Tibetan Plateau (SWTP)	50,231	4976	4344–6577	321	November (2)	August (99)
North Western Tibetan Plateau (NWTP)	88,002	5130	4625–6880	371	December (4)	August (85)
Western Karakoram (WK)	40,327	3612	278–8564	1021	October (50)	June (111)
Eastern Karakoram (EK)	31,244	5078	2287–7593	508	October (21)	August (62)
Northern Karakoram (NK)	32,099	4788	1860–8260	469	November (15)	August (85)
Altyn Tagh (AT)	126,556	2196	770–6280	163	November (4)	June (31)
Eastern Taklamakan (ET)	177,078	2196	906–6891	131	November (3)	June (24)
Central Kunluns (CK)	103,580	3623	1051–7051	265	December (4)	August (53)
Western Taklamakan (WT)	102,761	1644	1049–6316	174	December (3)	August (32)
All	1,235,357	2818	111–8564

Table 2. The 2001–2024 average snow-covered days (SCD), snow onset date (SOD), and snow end date (SED) for each subregion and for the whole area. The acronyms of the subregions follow the definition of Table 1. The mean elevation of each region is reported, too.

Region	Mean Elevation [m a.s.l.]	SCD (Days)	SOD (Day)	SED (Day)
GH	1974	48 ± 1	101 ± 1	148 ± 1
IH	1550	41 ± 1	102 ± 1	143 ± 1
NH	4643	136 ± 3	71 ± 2	206 ± 2
NT	5030	58 ± 3	97 ± 2	154 ± 2
ST	5068	35 ± 3	106 ± 1	140 ± 2
SWTP	4976	27 ± 3	106 ± 2	132 ± 1
NWTP	5130	39 ± 4	100 ± 3	139 ± 1
WK	3612	162 ± 2	57 ± 2	219 ± 2
EK	5078	184 ± 3	53 ± 2	236 ± 2
NK	4788	125 ± 2	75 ± 1	199 ± 2
AT	2196	19 ± 1	107 ± 1	126 ± 1
ET	2196	14 ± 1	110 ± 1	125 ± 1
CK	3623	53 ± 1	99 ± 1	151 ± 1
WT	1644	14 ± 2	110 ± 1	124 ± 1
All		46 ±1	100 ± 1	146 ± 1

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2025 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Almagioni, C.D.; Manara, V.; Diolaiuti, G.A.; Maugeri, M.; Spezza, A.; Fugazza, D. Snow Cover Variability and Trends over Karakoram, Western Himalaya and Kunlun Mountains During the MODIS Era (2001–2024). Remote Sens. 2025, 17, 914. https://doi.org/10.3390/rs17050914

AMA Style

Almagioni CD, Manara V, Diolaiuti GA, Maugeri M, Spezza A, Fugazza D. Snow Cover Variability and Trends over Karakoram, Western Himalaya and Kunlun Mountains During the MODIS Era (2001–2024). Remote Sensing. 2025; 17(5):914. https://doi.org/10.3390/rs17050914

Chicago/Turabian Style

Almagioni, Cecilia Delia, Veronica Manara, Guglielmina Adele Diolaiuti, Maurizio Maugeri, Alessia Spezza, and Davide Fugazza. 2025. "Snow Cover Variability and Trends over Karakoram, Western Himalaya and Kunlun Mountains During the MODIS Era (2001–2024)" Remote Sensing 17, no. 5: 914. https://doi.org/10.3390/rs17050914

APA Style

Almagioni, C. D., Manara, V., Diolaiuti, G. A., Maugeri, M., Spezza, A., & Fugazza, D. (2025). Snow Cover Variability and Trends over Karakoram, Western Himalaya and Kunlun Mountains During the MODIS Era (2001–2024). Remote Sensing, 17(5), 914. https://doi.org/10.3390/rs17050914

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu