A Response of Snow Cover to the Climate in the Northwest Himalaya (NWH) Using Satellite Products

Abstract

The Himalayan region is one of the most crucial mountain systems across the globe, which has significant importance in terms of the largest depository of snow and glaciers for fresh water supply, river runoff, hydropower, rich biodiversity, climate, and many more socioeconomic developments. This region directly or indirectly affects millions of lives and their livelihoods but has been considered one of the most climatically sensitive parts of the world. This study investigates the spatiotemporal variation in maximum extent of snow cover area (SCA) and its response to temperature, precipitation, and elevation over the northwest Himalaya (NWH) during 2000–2019. The analysis uses Moderate Resolution Imaging Spectroradiometer (MODIS)/Terra 8-day composite snow Cover product (MOD10A2), MODIS/Terra/V6 daily land surface temperature product (MOD11A1), Climate Hazards Infrared Precipitation with Station data (CHIRPS) precipitation product, and Shuttle Radar Topography Mission (SRTM) DEM product for the investigation. Modified Mann-Kendall (mMK) test and Spearman’s correlation methods were employed to examine the trends and the interrelationships between SCA and climatic parameters. Results indicate a significant increasing trend in annual mean SCA (663.88 km²/year) between 2000 and 2019. The seasonal and monthly analyses were also carried out for the study region. The Zone-wise analysis showed that the lower Himalaya (184.5 km²/year) and the middle Himalaya (232.1 km²/year) revealed significant increasing mean annual SCA trends. In contrast, the upper Himalaya showed no trend during the study period over the NWH region. Statistically significant negative correlation (−0.81) was observed between annual SCA and temperature, whereas a nonsignificant positive correlation (0.47) existed between annual SCA and precipitation in the past 20 years. It was also noticed that the SCA variability over the past 20 years has mainly been driven by temperature, whereas the influence of precipitation has been limited. A decline in average annual temperature (−0.039 °C/year) and a rise in precipitation (24.56 mm/year) was detected over the region. The results indicate that climate plays a vital role in controlling the SCA over the NWH region. The maximum and minimum snow cover frequency (SCF) was observed during the winter (74.42%) and monsoon (46.01%) season, respectively, while the average SCF was recorded to be 59.11% during the study period. Of the SCA, 54.81% had a SCF above 60% and could be considered as the perennial snow. The elevation-based analysis showed that 84% of the upper Himalaya (UH) experienced perennial snow, while the seasonal snow mostly dominated over the lower Himalaya (LH) and the middle Himalaya (MH).

1. Introduction

Snow cover is one of the most critical land surface parameters in global energy, the water cycle, and the development of the atmospheric boundary layer [1,2]. Thermodynamic properties and radiation characteristics of snow cover can modulate the land-atmosphere interaction, ecosystem functions, and the biological factors at various scales [3,4,5]. Snow cover controls the global and regional climate and releases fresh water in the hydrological cycle [6]. On a regional scale, it becomes more crucial for the river runoff, local water availability, and groundwater recharge that sustain millions of people’s lives and livelihoods [7,8]. Snowmelt variability may affect many sectors, such as hydropower, tourism, and agriculture [9,10,11].

The global warming and climate change scenario highlights the concern about the glaciers and snow cover variability [12]. The quantitative distribution of snow variability, mainly driven by meteorological parameters, depends on the latitude, altitude, and season [13,14]. Snow exhibits a strong negative correlation with atmospheric temperature [3,15]. Thus, a change in snow-covered areas is a direct indicator of climate change at both global and regional scales [16,17,18,19]. Bednorz [20] also reported that the relation between air temperature and precipitation affects snow characteristics. Li et al. [21] discovered that the snow cover is a more dominating factor than the vegetation cover in regulating the Earth’s heat budget on a regional scale. Changes in snow cover alter the land-surface heating and the overlying atmosphere, affecting the land-sea temperature gradient and the monsoonal rainfall intensity [22]. The effects of snow cover indicate that the snow cover in one region can modify the rainfall in the other areas of the world. A significant decrease in the snow cover extent has been observed in the northern hemisphere due to climate warming in recent decades [23]. Recent studies have shown that the snow onset and melt are mainly driven by air temperature, whereas a combined effect of temperature and precipitation influences during the winter [24]. Temperature increment in mountain regions doubles the global average and intensifies with elevation [14]. Brown and Mote [18] analyzed the snow cover response to climate change in the northern hemisphere and reported a complex snow cover response to increasing temperature and precipitation, varying with climate regimes, elevation, and other snow cover variables. They also noted that the snow cover duration has the most potent response to climate warming. Lemke et al. [25] studied the snow cover in the Northern Hemisphere from 1966 to 2005 and found a loss in the snow cover area (SCA) every month except November and December. Notarnicola [26] investigated the global snow cover scenario over mountain regions from 2000 to 2018 and discovered a decreasing trend in both snow cover duration and snow cover area. However, in some regions, they also observed an increased snow cover duration up to 32 days and snow cover area up to 11% in the northern hemisphere, particularly in winter.

Several studies have depicted that the Himalayan region is warming faster than the global average over the last century [27,28]. Shekhar et al. [29] revealed a decrease in snowfall with increased seasonal temperatures in the Indian Himalayan region. Numerous snow-related studies [30] have been conducted over the Himalayan region in the past few decades using remote sensing techniques. Different trends in SCA variation have been observed over the different parts of the NWH. Many studies have explored the Himalayan snow cover during the recent past and have reported a decreasing SCA trend in the NWH and some river basins, i.e., the Satluj, Jhelum, Chenab, and Beas river basins [31,32,33]. Singh et al. [16] examined the snow cover variability in the northwest Himalaya and its climatic zones and cited a nonsignificant increasing SCA trend for the NWH and all climatic zones from 2010 to 2016. In contrast to previous studies, they also observed a significant increasing trend over the upper Himalaya. Shafiq et al. [19] investigated the snow cover over the Kashmir Himalaya and concluded a significant increasing trend in annual mean SCA and decreasing trends in minimum and maximum temperature from 2000 to 2016. Several studies have also reported the increasing annual SCA trends over the Western Himalayas [11,34,35].

Monitoring and quantitatively assessing snow cover variability are of utmost importance due to their significant role in modulating socioeconomic sectors and climate [11]. Therefore, we must understand the present snow dynamics and its future trend. The present study examines SCA’s spatiotemporal variability in response to temperature, precipitation, and elevation over the NWH during 2000–2019 using satellite-based observations. Long-term trend analyses were performed using the Modified Mann-Kendall test at monthly, seasonal, and annual levels for the whole study region and the different elevation zones separately. Previous studies [12,19,35,36] analyzing the SCA variation with respect to temperature, precipitation, and DEM have mainly been limited to some river basins or other specific regions. The current study considers the entire NWH region for the most recent 20 years. Such a long-term study could open the path to a new dimension in the field of snow cover variability.

2. Materials and Methods

In the Himalayan region, the conventional field methods are very challenging for data collection due to its vast and mostly inaccessible areas and harsh environmental conditions. Satellite observations provide a more feasible and efficient option due to their synoptic and repetitive coverage and offer several products at different spatial and temporal scales. We employed the Python programming language and Google Earth Engine (GEE) platform to analyze different datasets used in this study.

2.1. Study Area

‘The Himalaya,’ which describes India’s north extent, is the largest source of snow and ice cover outside the polar regions in the world [6,37]. Northwest Himalaya is a significant part of the Indian Himalayan region, covering the Indian states of Jammu and Kashmir, Ladakh, Himachal Pradesh, and Uttarakhand and occupying an area of 331,694.92 km². The Himalaya has direct socioeconomic control over the Indian subcontinent. For example, it contributes 17% and 22% of the Gross Domestic Product (GDP) share in terms of agriculture and hydropower, respectively [33,38]. Several studies have revealed the influence of the Himalayan snow cover over the Indian monsoonal rainfall and its intensity [22,39,40]. A large portion of this region lies above 4000 m altitude (Figure 1b) and receives snow throughout the year [41]. Rivers originating from the higher Himalaya such as Ganga, Satluj, and Chenab, receive about 30–50% of the annual flow from its snowmelt runoff [7,42]. Major rivers and their tributaries originating from the northwest Himalaya (NWH) are key elements of agriculture systems, ecology, hydropower, and numerous socioeconomic and cultural developments in the Indian subcontinent [41,43]. The region receives snow generally from November to March, while snowmelt from April to June provides the predominant runoff source, which forms a significant constituent besides precipitation during the monsoon time [7,28].

2.2. MODIS Snow Product

MODIS is a widely used dataset for snow cover mapping [36,44,45,46]. Accuracy assessment and validation of SCA mapping using MODIS data has shown a good resemblance with ground observations [8,44,47,48,49]. For example, Klein and Barnett [50] compared the in situ SNOTEL (Snowpack Telemetry) measurements with MODIS snow products and described an overall accuracy of 94% in the Upper Rio Grande Basin. Similarly, Parajka and Blöschl [48] reported 95% accuracy of MODIS snow products with daily in situ observations of snow depth from 754 climatic stations over Austria.

This investigation used the most recent version (V6) of the Level-3 MODIS/Terra 8-day composite snow Cover product MOD10A2 [47,51], with a 500 m spatial resolution from March 2000 to February 2020. This product provides the maximum extent of snow cover over an 8-day period in 1200 km × 1200 km tiles. These tiles are generated by compositing MODIS/Terra daily snow cover product (MOD10A1) with 500 m resolution. In the MOD10A2 dataset, a pixel is considered as snow if snow appeared at least once during the 8-day period and cloud if cloud cover is observed on all 8 days. Preferring the 8-day composite dataset would be most useful for analyzing the maximum snow cover, which is a more useful parameter than the minimum or average snow cover. For example, minimum or average products would fail to map storm or some other event deposited new snow cover [52]. The composite dataset also minimizes cloud cover and has been useful in numerous research studies [53]. However, Thapa and Muhammad [54] investigated the original Aqua (MYD10A2), original Terra (MOD10A2), and Improved (cloud cover < 20%) Terra-Aqua (MOYDGL06*) 8-day composite snow products from 2003 to 2018 over the Karakoram Himalaya (part of the NWH). They reported a bias of −1.67% and 1.1% in the mean annual SCA between the improved and Aqua and Terra cloud threshold products, respectively. The MODIS snow cover product is based on a snow mapping algorithm that employs the Normalized Difference Snow Index (NDSI), as shown in the following equation:

$$\mathrm{NDSI}=\frac{\mathrm{VIS}\left(\mathrm{MODIS}\mathrm{band}4\right)-\mathrm{SWIR}\left(\mathrm{MODIS}\mathrm{band}6\right)}{\mathrm{VIS}\left(\mathrm{MODIS}\mathrm{band}4\right)+\mathrm{SWIR}\left(\mathrm{MODIS}\mathrm{band}6\right)}$$

(1)

The NDSI exploits the snow surface’s high reflectivity in the visible (VIS) band (band 4: 0.545–0.565 µm) and very low reflectivity in the shortwave infrared (SWIR) band (band 6: 1.628–1.652 µm). In contrast, the cloud has high reflectance in both the visible and near-infrared bands [55]. The different spectral characteristics of snow and cloud provide an efficient approach to classify the snow, cloud, and other cover types [56]. This algorithm uses the near-infrared (NIR) band 2 (>0.11) and band 4 reflectances (>0.10) along with the NDSI value (>0.40) to identify the snow pixels. Generally, NDSI > 0.4 is considered as snow, but in the forested regions, the NDSI and NDVI are employed together to discriminate between snow and non-snow (if the NDSI is less than the threshold (0.4) and NDVI is about 0.1, then the pixel is classified as snow). The dataset used here uses multiple criteria for snow retrieval, has underwent various quality checks, and has been found to be well suited globally. A brief description of the snow retrieval algorithm is provided by the authors of [51].

The whole study area covers three tiles of MOD10A2 product: h24v03, h24v04, and h24v06. There is a total of 916 images for each tile from March 2000 to February 2020. In this study, each year was considered from the March of a year to February of the next year. For example, the first year was considered from March 2000 to February 2001. The seasonal variations in SCA were analyzed as per the seasons classified by the Indian Meteorological Department (IMD): Winter (December to February), premonsoon (March to May), monsoon (June to September), and postmonsoon (October to November).

All files downloaded from https://search.earthdata.nasa.gov/search (accessed on 10 February 2021) in the GeoTIFF format were mosaiced and then reprojected from sinusoidal to geographic coordinate. The output images were clipped to the study area extent for further analysis.

2.3. Meteorological Data

The MODIS/Terra/MOD11A1 V6 product provides daily land surface temperatures (LST in K) and emissivity values in a 1200 km × 1200 km grid with a spatial resolution of 1 km. Negi et al. [57] analyzed the multidate satellite (MODIS)-derived snow surface temperature and validated with ground data of winter 2004–2005 and 2005–2006, collected at various observation stations located in different zones of the NWH. They stated good correlations for the upper Himalaya (0.98, 0.98), middle Himalaya (0.92, 0.96), and lower Himalaya (0.88, 0.82) for 2004–2005 and 2005–2006 winter, respectively. Climate Hazards Infrared Precipitation with Station data (CHIRPS) is around a 40-year quasi-global daily precipitation dataset [58]. It incorporates rain-gauge data with satellite imagery to generate a 0.05° resolution gridded precipitation (mm/day) dataset. Here, the precipitation dataset includes both rainfall and snow. Banerjee et al. [59] analyzed the rainfall trend and variability in the Uttarakhand Himalaya from 1983 to 2008 and highlighted that CHIRPS has a greater degree of correspondence with observed rainfall than the Precipitation Estimation from Remotely Sensed Information using Artificial Neural Networks-Climate Data Record (PERSIANN-CDR). Several previous studies have analyzed the CHIRPS dataset and reported satisfactory results with other model simulated, satellite, and ground-based observations [60,61,62,63]. Precipitation reflects mainly rainfall in the regions where snowfall is absent. The datasets are available at the GEE platform and have been processed for the study region and duration. Monthly mean surface temperature and precipitation are computed by averaging the daily/daytime temperature and summing the daily precipitation observations. Average monthly datasets are further converted into seasonal and annual data. The datasets are used to evaluate the trends and their interconnection with snow cover variability.

To verify the reliability of the meteorological parameters, we performed further comparisons. For instance, we compared the MODIS LST data with ERA5 skin temperature (see Appendix A, Figure A1a) and found a strong correlation (0.97) between them. Like MODIS (−0.039 °C/year), ERA5 (−0.044 °C/year) also shows a decreasing temperature trend over the NWH during the study period. Besides, both datasets are in close agreement with SCA variability (Appendix A,

Table A1. Spearman’s correlation coefficients between SCA (maximum extent) and ERA5 (skin temperature) and MODIS (LST).

	Annual	Pre-Monsoon	Monsoon	Post-Monsoon	Winter
ERA5	−0.708 * p < 0.001	−0.608 * p = 0.004	−0.681 * p < 0.001	−0.559 * p = 0.01	−0.632 * p = 0.002
MODIS	−0.81 * p < 0.001	−0.80 * p < 0.001	−0.45 * p = 0.043	−0.73 * p < 0.001	−0.56 * p = 0.01

* Statistically significant.

), even if their spatial resolutions are different (1 km pixel size for MODIS, 10 km grid for ERA5). Concerning precipitation, we also compared the CHIRPS observations with the PERSIANN-CDR (Appendix A, Figure A1b) over the study period and found that a strong correlation (0.866) exists between them. Both the datasets are in close agreement with SCA results (Appendix A,

Table A2. Spearman’s correlation coefficients between SCA (maximum extent) and CHIRPS and PERSIANN-CDR precipitation.

	Annual	Pre-Monsoon	Monsoon	Post-Monsoon	Winter
CHIRPS	0.47 * p = 0.03	0.41 p = 0.07	0.31 p = 0.17	0.74 * p < 0.001	0.22 p = 0.35
PERSIANN-CDR	0.42 p = 0.06	0.42 p = 0.066	0.19 p = 0.23	0.65 * p < 0.001	0.45 * p = 0.042

* Statistically significant.

).

2.4. Digital Elevation Model (DEM) Analysis

The SRTM [64] digital elevation V3 product is provided by NASA/USGS/JPL at a resolution of 30 m. This dataset is also available at the GEE platform and processed for the study region. The entire region is categorized into three elevation zones (Figure 1b): The lower Himalaya or LH (<2000 m), the middle Himalaya or MH (2000–4000 m), and the upper Himalaya or UH (>4000 m), contributing 23.4%, 23.85%, and 52.74% of the total area of the NWH, respectively.

2.5. Snow Cover Area (SCA)

In the MOD10A2 product, pixels having labeled values of 200 are considered as snow. The total count of snow pixels in the study area multiplied by pixel resolution computes the SCA for the 8-day period. A pixel that appears as snow in any of the four 8-day composite datasets present in a month is considered as snow for that month. The annual and seasonal mean SCA are computed from monthly SCA values. In the present study, SCA represents the maximum snow cover extent rather than the average SCA.

2.6. Snow Cover Frequency (SCF)

SCF is a quantitative assessment of SCA’s spatial distribution and variation over time, defined as the percentage of time an area remains snow-covered during a specific period, as shown in the following equation:

F = (F_snow/F_total) × 100

(2)

where F is the snow cover frequency (computed for each pixel) for a period, F_snow is the total number of months with snow, and F_total is the total number of months (in our study F_total = 240). Seasonal SCF is computed using monthly snow variation for a particular season.

2.7. Statistical Analysis Methods

2.7.1. Modified Mann-Kendall Test

The Mann-Kendall (MK) test is performed to analyze the statistical significance of snow cover, temperature, and precipitation trends, which is one of the robust tools to analyze the data over time for consistent trends. It is a nonparametric test, i.e., it works for all distributions. It is based on a null hypothesis, H₀, i.e., there is no monotonic trend, and is tested against 3 possible alternative hypotheses: (i) There is an increasing monotonic trend, (ii) there is a decreasing monotonic trend, or (iii) there is neither an increasing nor a decreasing monotonic trend. The null hypothesis is rejected when the p-value is less than the predefined significance level (0.05). If a series contains a serial correlation or seasonality effect [65], the original MK-test fails to consider it, which affects the trend test results. However, seasonality and autocorrelation effects exist in many real situations, such as climatic, hydrologic, air quality, and other natural timeseries. Hamed and Ramachandra Rao [66] proposed a modified MK-test (mMK) to address the serial correlation problem and suggested a variance correction approach to improve the trend analysis. For this study, the ‘pyMannkendall,’ a python package, was used, which has 11 different MK-tests and 2 Sen’s slope estimator functions [67]. Here, a combination of Hamed and Rao modified MK-test and Theil–Sen’s Slope Estimator was used to identify the trend and its magnitude, respectively, at different timescales (annual, seasonal, and monthly). Several studies have suggested that the modified MK-test provides consistent results with Sens’s slope [68,69,70].

2.7.2. Spearman’s Correlation Method

It is easy to calculate the correlation between 2 variables when both have a normal distribution that can be quantified using standard techniques such as Pearson’s correlation. For a nonparametric or non-Gaussian distribution, the rank-correlation method is preferred. This method first converts the values into rank data, assigns an integer rank value for each variable, and then computes the correlation coefficient between the 2 ranked variables. In this study, it was assumed that a monotonic relationship exists between the variables. For further analysis, Spearman’s correlation method [71] was used to evaluate the significance of the interrelationship between SCA and climatic parameters at the 0.05 significance level.

3. Results

3.1. SCA Analysis

Figure 2 summarizes the interannual variability of SCA between 2000 and 2019. The analysis shows that the average annual mean SCA (maximum extent) from 2000 to 2019 was 19.5 × 10⁴ km², accounting for 58.8% of the entire study region. During the same period, the mean surface temperature and precipitation were 11.06 °C and 600 mm, respectively. The years 2015, 2019, and 2009 showed relatively more snow cover, with 20.6 × 10⁴ km², 20.5 × 10⁴ km², and 20.3 × 10⁴ km² in an area on average, respectively, of which SCA in 2015 showed the maximum snow cover, with 1 × 10⁴ km² more area than the overall average (2000–2019). Whereas 2001, 2016, and 2000 indicated relatively low snow cover, SCA in 2000 was the lowest, with 17.9 × 10⁴ km² area and 1.6 × 10⁴ km² less than average. The years depicting high SCA values are in good agreement with low-temperature observations and vice versa. A sudden increase in precipitation was experienced after 2008–2009, which may have been a reason for the recently increased snow cover. However, during 2016–2017, a sudden peak in temperature was detected with a sudden dip in snow cover despite high precipitation. The sudden peak indicates that temperature has a stronger correlation with snow cover compared to precipitation. Dharpure et al. [72] also demonstrated a sudden peak in temperature in 2016 and highlighted that the SCA is more sensitive to temperature over the Chenab river basin in the western Himalaya. In a recent study, Krishnan et al. [73] inferred an increasing variability in western disturbances and expected a higher probability of snowfall in the Karakoram and western Himalayan regions, which will likely contribute to an increase in glacier mass.

Figure 3 summarizes the average maximum snow cover extent during each season in the NWH region from March 2000 to February 2020. Among seasons, the average SCA (maximum extent) in winter was the largest, reaching 23.9 × 10⁴ km² (72% of total area), followed by premonsoon (65%), postmonsoon (56.4%), and monsoon season (45.5%). The lowest average land surface temperature (−2.24 °C) was observed during the winter season, and the highest average surface temperature (20.86) during the monsoon period, which are in close resemblance with the highest SCA during winter and lowest during the monsoon.

In the premonsoon season, 2005 showed the maximum positive SCA anomaly (difference with respect to the overall average), with 9827 km² area, which is in strong agreement with temperature as 2005 indicated the lowest temperature during the premonsoon. Maximum precipitation and relatively low temperature in 2015 caused the second-highest positive anomaly during the premonsoon. The highest negative SCA anomaly in premonsoon season occurred in 2000, with 15,745.58 km², followed by 2001 and 2018, with 11,913.73 km² and 11,850.61 km² less area than average, showing a relatively high temperature and low precipitation (Figure 3a).

In the monsoon season, 2015 showed the largest SCA, with 15,002.4 km² more area than average, followed by 2019, with a positive anomaly of 13,098.2 km². The SCA in 2000 was least with 18,334.24 km² less area than average, followed by 2001 with a negative anomaly of 17,308.55 km² (Figure 3b).

In the postmonsoon season, 2018 showed the highest seasonal positive SCA anomaly, with a 26,870.37 km² area, followed by 2009, with 22,658.86 km² more area than average. In, 2016 lowest seasonal negative anomaly was observed, with a 36,918 km² area. In other years, postmonsoon absolute bias was less than 20,000 km² (Figure 3c).

In winter, the SCA was largest in 2001 with 12,411 km² more area than average, followed by 2019, with 10,111 km² area above the average. The winter’s lowest SCA was in 2000, having 12,828.58 km² below the average. However, the average SCA differences in winter were low compared to other seasons (Figure 3d).

SCA varies during the different months of the year due to a change in snow accumulation and ablation patterns. Monthly analysis over the past 20 years shows that January received the maximum SCA and minimum temperature, whereas August occupied the minimum SCA and second-highest temperature.

3.2. Trend Analysis

3.2.1. SCA Trend Analysis

SCA trend analyses at different scales were carried out for the period 2000–2019. The mMK-test shows a statistically significant increasing trend in annual SCA, with a rate of 663.88 km²/year (

Table 1. Summary of trend analysis of mean SCA (maximum extent) over the NWH during 2000–2019.

	Trend	p-Value	Slope (km2/year)
Annual *	Increasing	0.040	663.88
Premonsoon	No trend	0.205	339.73
Monsoon	No trend	0.183	720.12
Postmonsoon	No trend	0.205	728.36
Winter *	Increasing	0.006	453.39
March	No trend	0.314	320.46
April	No trend	0.417	355.14
May	No trend	0.581	290.67
June	No trend	0.085	732.55
July *	Increasing	0.029	1191.79
August	No trend	1	−25.20
September	No trend	0.922	−109.49
October	No trend	0.673	352.19
November	No trend	0.495	834.37
December	No trend	0.205	570.59
January	No trend	0.314	692.65
February	No trend	0.770	−192.15

* Statistically significant.

).

A positive SCA trend and a decreasing temperature trend over the region have also been observed in previous studies [11,16,35,72].

At the seasonal scale, only winter showed a significant increasing trend in SCA over the past 20 years, with a rate of 453.38 km²/year. Other seasons, excluding winter, did not show a clear significant increasing or decreasing trend (i.e., no trend) for the same period (Table 1).

On a monthly scale, July was the only month that showed a significant increasing trend in monthly SCA, with a 1191.97 km²/year rate over the study period. This may be due to the significant increasing precipitation trend in the July months compared to other months over the region. No significant trend was observed for the rest of the months (Table 1).

In this study, snow cover variability with elevation was also explored. The analysis shows statistically significant increasing SCA trends at all seasonal and annual scales in the lower Himalaya. Whereas in the middle Himalaya, only annual and monsoon SCAs showed significant increasing trends. However, the upper Himalaya did not show any clear trend either at an annual or seasonal scale (

Table 2. Zonewise SCA (maximum extent) trend analysis (slope in km²/year).

	Lower Himalaya	Middle Himalaya	Upper Himalaya
Annual	Increasing * p < 0.001 slope = 184.5	Increasing * p = 0.014 slope = 232.1	No trend p = 0.581 slope = 119.3
Premonsoon	Increasing * p = 0.021 slope = 135.1	No trend p = 0.267 slope = 203.6	No trend p = 1.0 slope = −0.474
Monsoon	Increasing * p = 0.002 slope = 147.4	Increasing * p = 0.034 slope = 251.6	No trend p = 0.314 slope = 188.9
Postmonsoon	Increasing * p = 0.012 slope = 182.9	No trend p = 0.205 slope = 337.4	No trend p = 0.581 slope = 291.6
Winter	Increasing * p = 0.001 slope = 241.5	No trend p = 0.205 slope = 163.8	No trend p = 0.217 slope = −52.158

* Statistically significant.

).

3.2.2. Temperature Trend Analysis

The mMK-test shows a nonsignificant decreasing trend in annual average temperature, with a rate of −0.04 °C/year for the whole NWH region during 2000–2019. At the seasonal scale, premonsoon (−0.06 °C/year) showed the highest rate of temperature change, followed by postmonsoon (−0.05 °C/year), winter (−0.02 °C/year), and monsoon (−0.01 °C/year). Among all the seasons, only the monsoon season indicated a significant decreasing temperature trend (

Table 3. Summary of temperature trend analysis over the NWH during 2000–2019.

	Trend	p-Value	Slope (°C/year)
Annual	No trend	0.085	−0.039
Premonsoon	No trend	0.346	−0.066
Monsoon *	Decreasing	0.001	−0.013
Postmonsoon	No trend	0.537	−0.053
Winter	No trend	0.673	−0.023

* Statistically significant.

), which may explain the significant increasing SCA trend during this season in the lower and middle Himalayan regions. A similar result was cited by the authors of [19] over the Kashmir Himalaya from 2000 to 2016 at a 0.05 significance level with −0.05 °C/year.

3.2.3. Precipitation Trend Analysis

The statistical analysis of precipitation over the NWH region reveals a significant increasing trend in the annual precipitation with a rate of 24.56 mm/year during the study period. The seasonal analysis shows that the monsoon had the highest increasing rate, followed by premonsoon, winter, and postmonsoon, with 16.39 mm/year, 5.41 mm/year, 2.75 mm/year, and 1.28 mm/year, respectively. In contrast to temperature, precipitation showed statistically significant increasing trends for all the seasons (

Table 4. Summary of precipitation trend analysis over the NWH during 2000–2019.

	Trend	p-Value	Slope (mm/year)
Annual *	Increasing	0.003	24.56
Premonsoon *	Increasing	<0.001	5.41
Monsoon *	Increasing	0.003	16.39
Postmonsoon *	Increasing	0.003	1.28
Winter *	Increasing	0.029	2.75

* Statistically significant.

).

3.3. Association of SCA with Climate Variables

The strength of the SCA’s association with temperature and precipitation was evaluated using Spearman’s correlation method. The analysis shows a significant, strong negative correlation between SCA and annual mean temperature, with a correlation coefficient (r) of −0.81 over the past 20 years (

Table 5. Spearman’s correlation coefficients between SCA (maximum extent) and climate variables over the NWH during 2000–2019.

	Annual	Pre-Monsoon	Monsoon	Post-Monsoon	Winter
Temperature	−0.81 * p < 0.001	−0.80 * p < 0.001	−0.45 * p = 0.043	−0.73 * p < 0.001	−0.56 * p = 0.01
Precipitation	0.47 * p = 0.03	0.41 p = 0.07	0.31 p = 0.17	0.74 * p < 0.001	0.22 p = 0.35

* Statistically significant.

). It implies that an increase or decrease in temperature indicates a decrease or increase in SCA. Figure 2 illustrates that the enhanced SCA coincides with decreased temperature and increased precipitation over the NWH region.

In a similar study, Islam and Rao [74] also inferred that periods of high rainfall coincide with periods of low temperature and vice versa over the Kashmir Himalaya (part of the NWH). The correlation coefficient value (r = 0.47) indicates a weak positive correlation between SCA and annual precipitation for the same period (Table 5). Shafiq et al. [19] inspected the snow cover variability with climate and cited a negative correlation (r = −0.49) between annual average temperature and mean SCA, while they found a positive correlation (r = 0.44) between annual rainfall and mean SCA over the Kashmir Himalaya during 2000–2016. In the same study, they also reported a decline in temperature (−0.05 °C/year) and a rise in precipitation (19.1 mm/year).

In the premonsoon season, SCA and temperature are strongly correlated; during the season, the temperature remains relatively high, and snow accumulated in winter starts melting. In contrast, a low and non-significant positive correlation exists between SCA and precipitation. During the monsoon, the correlation between SCA and temperature is lower than the other seasons, which is statistically significant at the 95% confidence level. Simultaneously, a nonsignificant positive correlation exists between SCA and precipitation, which suggests that temperature dominates over precipitation even in monsoon time over the NWH region. In postmonsoon time, SCA’s association with both temperature and precipitation is relatively high and significant. During this time, the temperature starts lowering, and the accumulated water starts converting into ice. Similarly, there is a negative correlation between SCA and temperature in winter, which is significant at a 95% confidence level. During the same period, the correlation between SCA and precipitation is lowest and non-significant.

3.4. Snow Cover Frequency (SCF) Analysis

The spatial distribution of monthly (over the entire 240 months, see Section 2.6) and seasonal SCF over the NWH was investigated for each pixel, and results are shown in Figure 4 and Figure 5.

From the results, it is evident that SCF was unevenly distributed over the NWH region. The mean monthly SCF was around 60%, and around one-fourth of the total study area had more than 90% SCF. Approximately half of the NWH had SCF of more than 90% for the premonsoon and winter seasons, while the average SCF during the same period remained 64.37% and 71.42%, respectively. The monsoon season had the lowest average SCF, and during this season, around 23% of the entire study region had SCF equal to or less than 10%. On a seasonal scale, the average SCF in premonsoon season was 64.37% of the whole study region, second to winter (71.42%), followed by postmonsoon (56.43%) and monsoon season with 46.01%. If we consider the SCF below 60% and above 60% to classify seasonal and perennial snow cover, the perennial snow cover was estimated to be 54.81%, and seasonal snow cover accounted for 45.19% of the total NWH region.

The whole study region is classified into three elevation zones: Upper (UH), middle (MH), and lower Himalaya (LH), where permanent, seasonal, and occasional snow is generally observed. Our focus was to find the effects of climate change on permanent, seasonal, and occasional snow in this study. The SCF distribution at different elevation zones is shown in Figure 6a. The average SCFs of LH, MH, and UH were 11.60%, 56.07%, and 81.22%, respectively. Around 65% of the total LH region had SCF less than 10%, accounting for approximately 15% of the total NWH. The majority of the LH (99.91%) and MH (56.49%) experienced seasonal snow, while 84% of the UH experienced perennial snow. Of the total UH, 46% had SCF more than 90%. One-fourth (25.16%) of the area of the entire NWH region had SCF of more than 90%, and most of that was in the UH (97.31%). The results show that higher SCF was associated with higher elevation (Figure 6b). SCF below 20% was observed mainly in the lower Himalaya (LH), SCF from 21–60% was experienced in the middle Himalaya (MH), and SFC > 60% was mostly present in the Upper Himalaya (UH).

4. Discussion

The annual and seasonal snow cover variability and its response to main climate variables were examined in depth over the NWH region from 2000–2019. The MODIS snow product was used for the analysis due to its higher accuracy and spatiotemporal coverage. The study revealed several fluctuations in annual average SCA (maximum extent), with a significant increasing magnitude from 2000 to 2019. The average annual temperature showed a decreasing trend (0.04 °C/year), and precipitation increased (24.56 mm/year) during the study period. Decreased temperature and increased precipitation are the main factors for the recent increasing SCA trends over the region. Most of the scientific reports have highlighted the NWH as the most vulnerable region due to the rapidly receding snow cover under the propensity of rising temperatures at a significant rate [27,28,29]. The snow cover fluctuations are in more substantial alignment with temperature fluctuations than the precipitation over the region. The influence of elevation on snow cover was also examined over the region. We perceived that the UH showed nonsignificant results, whereas the MH showed a positive SCA trend annually and during monsoon season. Statistically significant increasing SCA trends was observed annually and for all the seasons over the LH. The overall snow cover frequency at the higher altitudes was high compared to the lower altitudes, signifying the presence of perennial snow at the higher altitudes. Shekhar et al. [75] depicted a significant increasing trend in winter precipitation in the middle altitude (1500–3500 m) over the western Himalaya from 1971–2013, which may explain the increasing snow cover in the lower Himalaya.

Numerous other recent studies [12,33,34,76] in the Himalayan region have attempted to illustrate the snow cover variability and its response to climate for a better understanding of snow dynamics. Generally, these studies have focused on a specific part and period and have reported various aspects of snow variability. Shekhar et al. [29] reported a decline in total seasonal snowfall with increased seasonal mean temperature over the entire western Himalaya. Several other studies monitoring the Himalayan snow cover have reported a decreasing SCA trend over the NWH and some river basins, i.e., the Satluj, Jhelum, Chenab, and Beas river basins [31,32,33]. In contrast to the researches reporting a decreasing SCA trend, many studies have reported a positive SCA trend in different parts of the NWH. Shafiq et al. [19] examined the snow cover status over the Kashmir Himalaya from 2000 to 2016 and concluded a significant increasing annual mean SCA trend, with a positive trend in annual precipitation and decreasing trends in minimum and maximum temperature. Singh et al. [11] cited an increasing SCA trend between 2000 to 2011 over the Indus basin (part of the NWH). Many other studies have also revealed positive annual SCA trends over the western Himalaya [16,34,35]. Our observations support increasing SCA, precipitation, and decreasing surface temperature over the NWH region in the past 20 years.

The CHIRPS dataset provides total precipitation rather than snowfall and rainfall separately. The snowfall might have a direct impact on the SCA compared to the total precipitation, but snowfall is more prominent, especially in the higher altitudes. However, in the lower altitudes, snowfall and rainfall are both significant contributors to SCA. It has also been observed that CHIRPS data have a better performance in summer than winter [77], which should be taken into account in future studies for a better understanding of the correlation between SCA and precipitation.

Although we observed a statistically significant increasing SCA trend during the winter season, the constituent months did not show any significant trend over the NWH region individually. Further analysis of the winter SCA and the constituent months showed a nonsignificant trend from 2001–2018. This shows the influence of extreme SCA values at the beginning and end of our study period in the trend statistics. The component months SCA were analyzed in-depth, and we observed that in most years, December reported low SCA compared to winter. We also observed that mostly low SCA values in December were followed by higher SCA in January, and many times, February SCA compensated the both. Thus, at least one constituent month replenished the others, thereby maintaining an overall increasing trend during the winter. We also found a nonsignificant increasing trend for the annual, winter, and July SCA during 2001–2018, although all months showed a significant increasing trend during 2000–2019. So, the timeframe should be considered as a crucial factor in SCA trend analysis that may influence the robustness of the trend statistics. Therefore, future studies should include long-term SCA observations to examine and establish robust trend statistics.

Zonewise monthly SCA analysis (Appendix A, Figure A2) shows a continuous decline in SCA from May to August in both the UH and MH. During the same period (May to August), an increasing SCA was observed in the LH. The average SCA during May was around 73% of the total NWH region, and it continued to decrease until August (48%). During this period, the snowmelt directly impacted the downstream river flow, hydropower generation, agriculture, flood mitigation, and other environmental aspects.

Several other important aspects should be taken into consideration that may affect the results over the study region. Many studies have shown that snow impurities such as black carbon can significantly reduce the snow albedo [78,79,80], enhancing snow melting, which may alter the SCA and its relationship with temperature and precipitation.

Overall, our results suggest that SCA over the NWH has multimode variability and is influenced by the climatic variables. Further, an in-depth analysis is required for more information about the cause and predictability of the SCA in the past and future scenarios, which could be helpful for agriculture, climatology, and energy resource management.

5. Conclusions

The Himalaya, which houses the most extensive ice cover outside the poles, influences the environment at many scales, from the river runoff to the regional climate. In the above sections, we discussed how it affects millions of daily lives in many ways. The present study aimed to highlight the current scenario of snow cover over the NWH region so that the necessary steps could be taken accordingly and to motivate the research community to conduct further research. A multiscale analysis was carried out to understand the variability and trends in SCA and climatic parameters, especially temperature and precipitation over the region. The analysis revealed the strength of intercorrelation between SCA and climatic parameters. It was observed that SCA increased with a decrease in temperature and increased precipitation, whereas it decreased with an increase in temperature and a decrease in precipitation. We conclude that SCA is quite sensitive to the climate over the region. However, better understanding and accurately estimating the SCA analysis in the NWH region are still challenging. Besides, more in-depth research is needed to understand the consequences and their potential future implications. For instance, snow depth and snow physical parameters are necessary for a better understanding but not yet available as routine remote sensing products.