Our study covers the region from 25 to 45°N and from 65 to 105°E, running across some of the most densely populated regions of the world. Several key watersheds, such as the Indus, Syr Darya, Amu Darya, Yangtze, Salween, and Ganges/Brahmaputra drain from HMA (Figure 1A, blue outlines).

Precipitation in HMA is driven by three main weather systems—the Indian Summer Monsoon, the East Asian Summer Monsoon, and the Winter Westerly Disturbances (Bookhagen and Burbank, 2010; Cannon et al., 2016b). These major weather systems interact to bring a heterogeneous mix of snow and rain to different regions of HMA. Recent changes in the timing and intensity of precipitation from these major weather systems have been observed (e.g., Kitoh et al., 2013; Menon et al., 2013; Vaughan et al., 2013; Cannon et al., 2015; Singh et al., 2014; Cannon et al., 2016b; Malik et al., 2016; Norris et al., 2020), and have been shown to impact the timing and volume of snow-water storage and snowmelt (Kapnick et al., 2014; Smith et al., 2017; Smith and Bookhagen, 2018).

Snow has been extensively studied with passive microwave data—albeit at low spatial resolutions (e.g., 0.25° × 0.25°) (Chang et al., 1987; Kelly et al., 2003; Smith and Bookhagen, 2016; Smith and Bookhagen, 2018). Recent image processing advances have allowed researchers to take advantage of the elliptical nature of passive microwave footprints to re-process the data onto a much finer spatial grid than previous approaches had allowed. In this study, we use the EASE-grid 2.0 high-resolution passive microwave product (1987–2016) (Early and Long, 2001; Brodzik et al., 2012; Brodzik et al., 2016; Long and Brodzik, 2016), which provides the 19 and 37 GHz passive microwave frequencies at spatial resolutions of 6.25 and 3.125 km, respectively. This dataset has been carefully cross-calibrated between the various SSMI and SSMI/S satellite platforms to provide consistent and homogenized data through the entire time series (Brodzik et al., 2016).

To produce consistent snow-water equivalent (SWE) estimates over the entire study region, we further re-grid the 19 GHz passive microwave data to a 3.125 km spatial resolution. We then remove areas near lakes and areas with shallow or infrequent snow-cover, as these areas are not suitable for long-term SWE trend analysis. Finally, following the methodology of Smith and Bookhagen (2018), we use the computationally efficient algorithm proposed by Chang et al. (1987) to convert the passive microwave data to snow depth: S D = 1.59 × ( T b 19 V − T b 36 V ) (1)

We then convert those snow depth estimates to SWE using a constant snow density of 0.24 g/cm³, which has been shown to be a reasonable global average (Sturm et al., 2010; Takala et al., 2011). In short, the Chang et al. (1987) algorithm uses the difference between the 19 and 37 GHz passive microwave channels to estimate snow depth based on a comparison with extensive snow survey data collected throughout the Canadian and Russian Arctic (Chang et al., 1982; Chang et al., 1987). This algorithm is widely used to estimate SWE over diverse terrain, and has served as the basis for further updates to passive microwave SWE retrieval algorithms which take advantage of additional passive microwave channels not carried on SSMI/S to better constrain the impacts of vegetation cover on SWE estimates (Chang et al., 1991; Chang et al., 1996; Foster et al., 2005; Derksen, 2008; Kelly, 2009; Langlois et al., 2011; Smith and Bookhagen, 2016). In our low-vegetation study area, we rely on the Chang et al. (1987) algorithm to take advantage of the full SSMI/S time series.

For parts of the passive microwave time series, there are multiple overlapping satellite overpasses. For consistency, we aggregate all night-time overpasses (October 1987–September 2016) into an average daily SWE estimate over HMA (number of grid cells = 1,027,847) using the median of all available night-time measurements per day. As the various SSMI satellite platforms have been carefully cross-calibrated (Brodzik et al., 2016), this step serves simply to homogenize the temporal sampling of the dataset over the entire time period. For computational efficiency, we then further resample each individual daily SWE time series to a temporal frequency of three days before computing trends; in our tests this does not significantly modify computed long-term trends.

Before trend analysis, we first remove the seasonal component of each individual time series via Seasonal Trend Decomposition by Loess (Cleveland et al., 1990), using a decomposition window of 365 days (Figure 2). This method yields a seasonal signal, long-term signal, and residual short-term signal from a given time series by removing oscillations at the chosen decomposition time frequency. We then test the resulting de-seasoned time series for significant increasing or decreasing trends using the Mann-Kendall test (Mann, 1945; Kendall, 1948). If there exists a significant trend, we use Sen’s slope method to capture the overall trend at that grid cell (Sen, 1968). We thus use a conservative approach by testing for significance both with the Mann-Kendall test and via Sen’s slope method. We only present results from statistically significant (p < 0.05) trends in this study.

Aggregate trends over the entirety of HMA are slightly negative (sum: −55.5 mm/yr, average: −0.01%) over 3,618 km²× 10³ km², including only trends with p < 0.05 and areas at least 500 m above sea level). While the aggregate trends appear to be small, we emphasize that trends are measured over 9.75 km² grid cells, and represent a snow-water storage loss of 5.41 m³ × 10⁵ m³ of water per year (Figure 3).

It is clear that SWE trends are spatially diverse—positive SWE trends are concentrated in the Karakoram, Pamir, Kunlun Shan, and the high Himalaya (Figures 3 and 4A–D). The most negative trends are concentrated in the south-eastern Tibetan Plateau (Figures 4E,F), at the headwaters of the Yangtze, Salween, and Mekong rivers. There also exist many small-scale features; for example, there are clear alternating positive-negative SWE trend patterns along the front of the Himalaya. While there are multiple possible causes for such small-scale variability, the extreme topography and the microclimates it creates can drastically alter snow-loading and snow-water storage on nearby slopes—particularly when there are large differences in the sun exposure and overall aspect of neighboring slopes (Supplemental Figures S1–S5).

When the SWE data are further divided into seasonal components, clear differences in trend direction and magnitude appear (Figure 5). Strong positive winter (December-January-February) trends are visible in most of the highest peaks of HMA—running along the Himalaya, into the Pamir-Karakoram-Kunlun Shan region, as well as through the Tien Shan. These positive trends are visible through the spring (March-April-May) and fall (September-October-November) seasons as well; the only region, however, to maintain positive trends through the full year and in each seasonal slice is the Karakoram-Kunlun Shan region, which has been noted for glacier stability and growth in recent years (Hewitt, 2005; Kääb et al., 2012; Kapnick et al., 2014; Treichler et al., 2019; Shean et al., 2020). These positive trends are offset by large negative SWE trends in lower-elevation regions of HMA and along the eastern edge of the Tibetan Plateau which has seen rapidly decreasing SWE—particularly in the December-January-February and September-October-November periods (Figure 5).

To further explore the dynamics of SWE trends in HMA, we have performed a second set of regressions using monthly SWE percentiles (Figure 6). In short, we calculate the 10th, 25th, 50th, 75th, and 90th percentile SWE value at each pixel over each month using daily-averaged SWE data (October 1987–September 2016), remove the long-term monthly mean value for each given month to reduce the impacts of seasonality, and perform regressions through each SWE percentile separately. This yields a set of SWE trend results based on only the lowest (e.g., 10th percentile) or highest (e.g., 90th percentile) SWE value for each month.

When the SWE trend magnitudes at each percentile are compared, differences between high- and low-percentile trends are apparent (Figure 6). In almost all cases, the trends in high-percentile SWE are steeper than those in low-percentile SWE. In positive SWE-trend regions, this indicates that high SWE amounts are becoming relatively more frequent. For example, along the border of India and Pakistan, positive SWE-trend regions (see Figure 3) have a more than six-fold higher trends in 90th percentile monthly SWE than in 10th percentile SWE, indicating a drastic increase in monthly high-SWE day frequency or magnitude. This could be due to the increased strength of the Winter Westerlies, which has been previously reported (Cannon et al., 2016b). In contrast, positive SWE-trend areas in the Pamir have higher 10th percentile magnitudes than 90th, indicating a that positive SWE trends are driven by increases in low-magnitude SWE rather than in high-magnitude SWE.

In the majority of HMA, however, SWE trends are negative (Figure 3). Thus, the high 90th/10th (75th/25th) percentile trend ratios indicate that declines in SWE have been steeper in the higher end of the monthly SWE distribution, and that high-SWE values are becoming less common overall. This agrees well with previous research, which points to overall increasing temperatures in HMA, particularly on the Tibetan Plateau (e.g., Wang et al., 2018).

Previous work by Smith and Bookhagen (2018) used data at 0.25° × 0.25° spatial resolution from only the SSMI-series of satellites (1987–2009) to establish trends in SWE over HMA. The higher spatial resolution data used in this study yields only slight differences in SWE trend when the same period (e.g., 1987–2009) is considered. However, there are clear differences in the trends presented by Smith and Bookhagen (2018) and those shown in Figure 3, which are due to the difference in analysis time window. To test the sensitivity of SWE trends to the analysis window, we first break our dataset into three decade-long slices, as seen in Figure 7.

It is clear that SWE trends are highly variable in time. The long-term reversal from negative to positive SWE trends seen in Figure 7, however, is supported by analysis of other related climate variables. Previous work has noted changes in regional precipitation and temperature patterns (e.g., Archer and Fowler, 2004; Yao et al., 2012; Palazzi et al., 2013; Lutz et al., 2014; Cannon et al., 2016a, Cannon et al., 2016b; Zhang et al., 2017; Wang et al., 2018; Treichler et al., 2019; Norris et al., 2020) and increases in high-elevation snowcover (Kapnick et al., 2014; Tahir et al., 2015) in recent years. Furthermore, Treichler et al. (2019) showed that increasing lake levels on the Tibetan Plateau are strongly correlated with regions of increased precipitation; modeled precipitation data suggest stepwise increases in mean annual precipitation on the Tibetan Plateau between the ∼1980s–1990s and 2000s–onwards (e.g., Kääb et al., 2018).

Recent analysis also indicates that the timing of the snowmelt season has changed over the past decades (Smith et al., 2017). While the long-term trends (1987–2016) were found to be generally negative (e.g., earlier and more rapid snowmelt), recent trends (e.g., 2004–2016) were found to be much more positive (later and slower onset of snowmelt) (Smith et al., 2017). It is therefore possible that there has been a reversal of the long-term losses in SWE storage in HMA; however, it is not clear if this is a temporary or long-term shift in the snow dynamics of HMA.

The timing and magnitude of SWE trend variability can be further explored by performing the same trend analysis on a set of time windows and start dates. We use time windows of 5, 10, 15, 20, 25, and 29 years, along with each possible combination of start years (e.g., 1987–2011) to examine changes in SWE trends through time (Figure 8).

Long-term (>25 years) SWE trends are generally negative; these trends also correspond to a particularly negative period of short-term trends starting in the late 1980s (Figure 8). More recent trends (e.g., 15–20 years) are generally positive starting in the early 1990s. One possible explanation for this phenomena is the previously proposed large-scale changes in regional precipitation over the past decades (e.g., Kääb et al., 2018). However, the impacts of changes in temperature cannot be ruled out—increasing regional temperatures can have highly variable positive and negative impacts on snow-water storage, for example, by enhancing atmospheric water content, snow density, and snowmelt rates.

There also exist strong regional variations in windowed trends (Figures 8B–D), driven by differences in climatic conditions, major weather systems, snow accumulation and ablation regimes, and dust and aerosol melt forcing between regions (Fujita, 2008; Kaspari et al., 2014; Sarangi et al., 2019). Generally positive SWE trends in the Kunlun Shan region are contrasted by mixed trends in the Karakoram, and majority negative trends in Eastern Tibet (Figures 3 and 4).

Many recent studies have investigated changes in HMA’s glaciers using a range of satellite (Bolch et al., 2012; Kääb et al., 2012; Loomis et al., 2019; Treichler et al., 2019; Shean et al., 2020) and modeling (Kapnick et al., 2014; Rounce et al., 2019) approaches to derive spatial patterns in glacier gains and losses. Using the Randolph Glacier Inventory (Arendt et al., 2015), we can measure the areal extent of glaciers within each passive microwave pixel, and—where glaciers are large enough—derive SWE trends over only glaciated areas, defined here as areas with at least 10% glacier coverage (Figure 9).

In general, SWE trends over glaciated terrain are negative, outside of parts of the Tien Shan, Karakoram, and Kunlun Shan. Areas with positive SWE trends agree well with regions of positive glacier mass balance, as presented by Shean et al. (2020) and Treichler et al. (2019). While there are many factors that influence glacier dynamics, it is likely that changes in snowfall are one of the key drivers of glacier mass gain and loss over HMA (Fujita, 2008; Fujita and Nuimura, 2011; Kapnick et al., 2014).

It is important to mention caveats to the trend analysis presented in this study. The largest caveat is that passive microwave SWE estimates are often uncertain—especially over large and complex regions such as HMA (Kelly, 2009; Takala et al., 2011; Smith and Bookhagen, 2016). We also cannot rule out the impacts of both natural seasonality and regional temperature changes on snow densities, which could also modify passive microwave SWE estimates over the course of our time series, and thus are part of the trends that we present as changes in SWE in this study (Judson and Doesken, 2000; Chen et al., 2011; Dai et al., 2012). The impact of seasonal oscillations in snow density is somewhat mitigated by removing the seasonal cycle from our data before trend fitting, as some of the seasonality in SWE estimates will be driven by changes in snow density. However, without a more in-depth understanding of snow-density evolution in HMA, we cannot fully constrain this part of our analysis.

Passive microwave signal saturation could bias the presented SWE trends in deep-snow areas, as previous work has suggested that passive microwave SWE estimates saturate around 200 mm of SWE (Takala et al., 2011; Vander Jagt et al., 2013; Dozier et al., 2016). In examining our dataset, we find that the vast majority of HMA is not severely impacted by passive microwave signal saturation (Figure 10); however, it is likely that passive microwave signal saturation still biases some of our trend results. In our percentile regressions, we find that SWE trends generally maintain a consistent direction between high- and low-percentiles, indicating that saturation biases don’t drastically influence SWE trend direction (Figure 6).

The third caveat is that trends are somewhat biased by the first and last values of the time series—this could also play a role in the trend reversals seen between previous studies of SWE trends in HMA (e.g., Smith and Bookhagen, 2018; Wang et al., 2018) and this study (Figures 3 and 7). We attempt to minimize this effect by using Sen’s slope estimator, which is less sensitive to the first and last values of the time series (Sen, 1968). We further attempt to mitigate the impact of the time window over which the trend is calculated by using multiple overlapping time windows and window lengths (Figure 8). Finally, we compare our results to a percentile-regression approach and find that while the magnitudes of the trends vary between percentiles and between the de-seasoned trend and percentile approaches, the trend directions found are consistent. However, snowfall can have high inter-annual variability, and we do not preclude the possibility that variations in the timing of large snowfall events between years, or shifts in the timing of snowfall and snowmelt (Smith et al., 2017) could impact estimated annual and seasonal SWE trends.

Lastly, it is important to emphasize that the trends presented here are relative to the SWE time series as estimated, and are not calibrated by in-situ measurements. While the SWE algorithms used here have been extensively validated throughout the world and have been shown to be generally reliable in low-vegetation areas (Chang et al., 1982; Chang et al., 1987; Chang et al., 1991; Chang et al., 1996; Armstrong and Brodzik, 2002; Foster et al., 2005; Derksen, 2008; Kelly, 2009; Langlois et al., 2011; Dai et al., 2015; Smith and Bookhagen, 2016), the complex topography and inaccessibility of HMA poses unique challenges for in-situ data collection. Unfortunately, calibration data of sufficient spatial and temporal resolution to properly assess our SWE estimates and trend results is not currently available in our study region. Future work could consider other approaches to constraining the estimated SWE trends, for example, by using watershed-level snowmelt runoff measurements across HMA.