Vegetation expansion in the subnival Hindu Kush Himalaya

Abstract The mountain systems of the Hindu Kush Himalaya (HKH) are changing rapidly due to climatic change, but an overlooked component is the subnival ecosystem (between the treeline and snow line), characterized by short‐stature plants and seasonal snow. Basic information about subnival vegetation distribution and rates of ecosystem change are not known, yet such information is needed to understand relationships between subnival ecology and water/carbon cycles. We show that HKH subnival ecosystems cover five to 15 times the area of permanent glaciers and snow, highlighting their eco‐hydrological importance. Using satellite data from the Landsat 5, 7 and 8 missions, we measured change in the spatial extent of subnival vegetation from 1993 to 2018. The Landsat surface reflectance‐derived Normalized Difference Vegetation Index product was thresholded at 0.1 to indicate the presence/absence of vegetation. Using this product, the strength and direction of time‐series trends in the green pixel fraction were measured within three regions of interest. We controlled for cloud cover, snow cover and evaluated the impact of sensor radiometric differences between Landsat 7 and Landsat 8. Using Google Earth Engine to expedite data processing tasks, we show that there has been a weakly positive increase in the extent of subnival vegetation since 1993. Strongest and most significant trends were found in the height region of 5,000–5,500 m a.s.l. across the HKH extent: R 2 = .302, Kendall's τ = 0.424, p < .05, but this varied regionally, with height, and according to the sensors included in the time series. Positive trends at lower elevations occurred on steeper slopes whilst at higher elevations, flatter areas exhibited stronger trends. We validated our findings using online photographs. Subnival ecological changes have likely impacted HKH carbon and water cycles with impacts on millions of people living downstream, but the strength and direction of impacts of vegetation expansion remain unknown.


| INTRODUC TI ON
Mountain systems are amongst the most dynamic on Earth and are particularly sensitive to climatic change (Dolezal et al., 2016).
The subnival zone, which lies above the treeline and below the permanent snowline, is the most poorly studied vegetation zone in the HKH due to its inaccessibility and high altitude. In this region, ecological responses to climate change are likely to be multidimensional, species-specific and spatially variable (Dolezal et al., 2016) with strong potential for feedbacks between plant cover, carbon, snow and hydrological processes to occur, as evidenced by relatively advanced work from Arctic systems (Myers-Smith et al., 2011Myers-Smith & Hik, 2013). In the HKH, feedbacks have been shown to exist between climate and plant phenology (Chen, Zhu, Wu, Wang, & Peng, 2011) and between shrub cover and snowpack dynamics (Hu et al., 2009). There is, meanwhile, evidence of shrub encroachment into Himalayan grasslands at high altitude (Brandt, Haynes, Kuemmerle, Waller, & Radeloff, 2013;Chophyel, 2009;Qiu, 2016;Wangchuk, Gyaltshen, Yonten, Nirola, & Tshering, 2013), and warming-driven geographical range shifts in endemic plant species (Dolezal et al., 2016;Telwala, Brook, Manish, & Pandit, 2013).
Spatial ecological changes across the full extent of the HKH subnival zone now need scientific attention as we approach, and pass, peak non-renewable water (Gleick & Palaniappan, 2010;Huss & Hock, 2018;Jones, Harrison, Anderson, & Whalley, 2019), and yet across the HKH extent, there is no basic understanding of current subnival vegetation distribution, and no information on rates of change over past decades.
A starting point from which to evaluate the present status, geographical distribution and historical trend in vegetation change across the HKH is to exploit freely available medium resolution (30 m pixel size) satellite data. Whilst this seems relatively straightforward in principle, there are practical challenges caused by the huge data volumes needed to evaluate changes in vegetation cover over thousands of square kilometres and over multi-decadal timescales. Classical remote sensing workflows where satellite data are downloaded from online repositories and then processed locally on computers, are computationally-and time-expensive which prohibits analyses on an area this large. To address this challenge, we exploited the unique freely accessible cloud-processing capabilities of the new Google Earth Engine (GEE) platform, which houses public data catalogues for the whole of the Landsat archive (Gorelick et al., 2017). We focus on answering two questions: 1. What is the extent of the subnival zone? 2. Has the spatial extent of subnival vegetation changed and, if so, at what rate and where?
With regard to answering question 2, our work was concerned with evaluating evidence for subnival system transition from bare ground or sparsely vegetated ground to vegetated ground; not change in species composition.

| Study system: The subnival zone
Ecosystems in the HKH region are diverse due to the variety of complex climatic conditions, alongside stark variations in altitude (Dorji, Olesen, Bøcher, & Seidenkrantz, 2016) and slope aspect (Anthwal, Bhatt, Nautiyal, & Anthwal, 2011). The treeline in the Himalaya is typically found at approximately 4,000 m above sea level (m a.s.l.; Gaire, Koirala, Bhuju, & Borgaonkar, 2014) although certain tree species have been shown to grow at altitudes of up to 4,900 m a.s.l.
A useful description of the major Himalayan plant zones is provided by Ives and Messerli (1989), who describe the upper timberline as occurring between 4,000 and 4,500 m a.s.l., above which a rhododendron-shrub belt grades into alpine meadows, with a 'subnival belt of extensive bare ground and scattered dwarf plants, mosses and lichens'. They describe the 5,000-5,500 m a.s.l. elevation zone as being characterized by permanent ice and snow with steep rocky outcrops. In this work, we use the upper elevational limit of B. utilis

| Satellite data analysis
We utilized data from the NASA Landsat satellite mission, which provides long time-series (1972-present) medium resolution (~30 m pixels) imaging data over the region (Cohen & Goward, 2004;Roy et al., 2014). Such data have a long history of use for monitoring land cover change (Aplin, 2004) including in mountainous regions (Cingolani, Renison, Zak, & Cabido, 2004;Helmer, Brown, & Cohen, 2000;Levin, Shmida, Levanoni, Tamari, & Kark, 2007;White, Running, Nemani, Keane, & Ryan, 1997), in areas of short-sward vegetation (Chopping et al., 2008;Goslee, Havstad, Peters, Rango, & Schlesinger, 2003;Laliberte et al., 2004;Skowno et al., 2017;Tian, Brandt, Liu, Rasmussen, & Fensholt, 2017), and at sites within the HKH (Brandt et al., 2013). Landsat missions also offer long time-series data with a similar sensor, a significant benefit because some work has shown the lack of consistency when comparing vegetation metrics between different sensors (Abuzar, Sheffield, Whitfield, O'connell, & Mcallister, 2014). Indeed, a recent study suggested that freely available data from satellites such as NASA's Landsat were 'remarkably congruent' for broadscale monitoring of shrub encroachment in African savannahs (Marston, Aplin, Wilkinson, Field, & O'regan, 2017), which have short-stature, patchy vegetation structurally similar to that within the HKH subnival zone. At 30 m pixel resolution, Landsat data offer considerable benefits for vegetation time-series analysis in mountainous areas, particularly as compared to other satellite-derived vegetation products. For example, the GIMMS 3g fAPAR product used by Keenan and Riley (2018) for evaluation of global cold climate ecosystem change has a spatial resolution of 8 km, such that, across the HKH, one would find significant variations in elevation and ecology over the space of a single pixel. Corroborating this, recent work by Fassnacht, Schiller, Kattenborn, Zhao, and Qu (2019) evidences the lack of suitability of coarse-grained 500 m resolution MODIS vegetation products for mapping changes in high altitude vegetation on the Tibetan Plateau. In contrast, they showed that Landsat data could allow for 'detailed spatial patterns of land cover changes' to be identified (Fassnacht et al., 2019). Importantly, Fassnacht et al. (2019) corrected for radiometric differences between Landsat 7 Enhanced Thematic Mapper (TM) and Landsat 8 Operational Land Imager (OLI) detectors using published coefficients from Roy et al. (2016).
Our study analysed data from three Landsat missions (5, 7 and 8) to evaluate time-series dynamics in the presence/absence of subnival vegetation cover over the 26 year period from 1993 to 2018.
Data from the period prior to 1993 were not included due to patchy time-series coverage and because of the complexity of issues associated with integrating data from earlier missions (i.e. Landsat 1-4).
We subsequently split the analysis to study vegetation changes across three regions encompassing different spatial extents: (a) Sagarmatha (Everest) National Park, Nepal; (b) the national extent of Nepal; and (c) the entire HKH region. For each of these regions, we measured vegetation trends in four elevation bands: • H1 (4,150-4,500 m a.s.l.); • H2 (4,500-5,000 m a.s.l.); • H3 (5,000-5,500 m a.s.l.); • H4 (5,500-6,000 m a.s.l.). Subsequent sections explain how the GEE methodology was implemented for each of these in turn.

| Sagarmatha National Park area, Nepal ('P140-R40/41 region')
The first part of the analysis focused on analysing time-series change in two tiles of Landsat data (path 140, rows 40 and 41; Figure 2), covering a region on the Nepal/Tibet border centred on Mount Everest. This was chosen for initial exploration for several reasons. This region contains the greatest elevation range in the Himalayas, and includes areas of glacierized (i.e. ice-covered) terrain and recently deglaciated terrain (including the Khumbu region, south of Everest; Byers, 2007). In the northernmost part of the image extent are areas on the Tibetan plateau exceeding 6,000 m a.s.l. where previous reports have suggested vegetation change is occurring rapidly (e.g. Qiu, 2016). Furthermore, these data also cover a broad range of climatic conditions including precipitation gradients caused by winds and orographic barriers (Bookhagen & Burbank, 2006). The area thus provides a good starting point for investigating vegetation change because it is a microcosm of the range of conditions (particularly elevation and aridity) experienced across the wider HKH. Figure 2 shows the location of the two tiles of Landsat data analysed in this initial phase, and throughout the rest of this paper, we refer to this region as P140-R40/41.

| National extent of Nepal
The second part of the analysis focused on time-series change across the areal extent of the country of Nepal. For this, a KMZ shapefile delineating the boundary of Nepal was uploaded into GEE and used to constrain the analysis.

| HKH-wide analysis
Google Earth Engine allocates users a fixed processing capacity, so to be able to measure change over the entire HKH, a random sampling method using regions of interest (ROIs) was necessary. We defined 100 circular ROIs with a 5 km radius and randomly deployed these within each of the four height bands previously described. The total area covered by the ROIs that were used to sample the satellite data record equalled 31,416 km 2 (overlap of ROIs and height bands not taken into account). Figure 3 shows the spatial distribution of the different height bands across the HKH sampled using the circular ROIs.

| Permanent snow to seasonal snow cover calculations
To answer question 1 (What is the extent of the HKH subnival zone?), the freely available moderate resolution imaging spectrometer (MODIS) F I G U R E 2 Location of the P140-R40/41 region of interest, which utilized data from Landsat path 140, rows 41 and 42 on the Nepal/Tibet border [Colour figure can be viewed at wileyonlinelibrary.com] fractional snow product (Hall, Riggs, Salomonson, Digirolamo, & Bayr, 2002) was used and analyses were run at two scales: first for the P140/R40-41 region, and second, for the country of Nepal.
We selected recent years (2013-2017) to generate a product describing the median snow-covered area in late summer (August and September), when the snow cover is at a minimum (determined using: http://geoap ps.icimod.org/HKHSn owCov er/) averaged over a 1 km pixel scale. This output was used to represent permanent snow-covered areas for the ROI. The permanent snow-covered area was combined with the Randolph Glacier Inventory (RGI ;Pfeffer et al., 2014; which was at a finer spatial resolution) to estimate the spatial extent of permanent snow and ice cover. This was compared to the extent of the entire subnival zone, represented by the total area above 4,150 m a.s.l. using the Shuttle Radar Topography Mission (SRTM) 30 m gridded data set as a measurement of height above mean sea level (Farr et al., 2007). To estimate the total area of permanent snow and ice, the snow area on glaciers was subtracted from permanent snow cover + permanent ice cover. The final ratio between the subnival habitat area which can contain seasonal snow and ice, and permanent snow and ice were calculated by first, subtracting the total area of permanent snow and ice from the total land area above 4,150 m a.s.l., and then dividing the two values.

| Landsat surface reflectance product procedure for vegetation fractional cover determination
To answer question 2 (Has the spatial extent of subnival vegetation changed; and if so, at what rate and where?), Landsat Tier 1 surface reflectance (SR) data processed by USGS were used. The SR products contain pixel quality attribute information resulting from the atmospheric correction process. Likely cloud, snow and shadow pixels are flagged and can be used to mask undesirable data, whilst the remaining SR values are expected to be more consistent than standard top of atmosphere radiance products, because differences in the atmospheric composition are taken into account. However, USGS caution that SR products could be adversely affected by hyper-arid or snow-covered conditions, which exist in certain regions of the area studied here (USGS, 2018). We needed a simple vegetation product to allow vegetated pixels to be discriminated from unvegetated pixels, but a complex land cover analysis or retrieval of biophysical vegetation properties was not required. For this reason, a simple and widely used index, the Normalized Difference Vegetation Index (NDVI), was chosen. Other studies have shown that NDVI is a relatively robust indicator of green biomass even in sparsely vegetated systems (Gamon et al., 1995), and that specifically, Landsat TM and enhanced TM (ETM+) derived NDVI delivered consistent data to surfacemeasured NDVI in Arctic systems (Pattison, Jorgenson, Raynolds, & Welker, 2015), which have similar low-stature plants. Recent work by Fassnacht et al. (2019) on the Tibetan plateau has also shown the suitability of Landsat NDVI products for measuring high altitude vegetation change. For these reasons, we argue that NDVI derived from SR data will produce a robust presence/absence indicator for vegetation in subnival HKH systems. The description of the applied workflow is given in the below paragraph with key factors marked in italic.
First, data from Landsat 5, 7 and 8 Tier 1 SR image collections were filtered to extract data from specific post-monsoon months (October and November). This period in the post-monsoon season is when the HKH does not receive a significant volume of rainfall and when cloud and snow cover are expected to be at an annual minimum (Bookhagen & Burbank, 2006). By constraining the analysis to this short post-monsoon period, phenological differences between individual acquisitions were also limited. Pixels flagged as cloud, snow or shadow in the perpixel quality assessment (QA) layer were then masked in the individual images. To verify whether there were any significant changes in snow fractional cover over time over the years analysed, a snow cover product based on pixels flagged as snow within the Landsat QA flags was also generated to derive a fraction of snow-covered pixels per year for the three spatial extents considered throughout this paper.
Landsat 5, 7 and 8 collections were then merged and NDVI was computed for each image in the collection. To composite multiple images available for the observed time period and fill regions of missing data due to cloud cover, an annual vegetation index average was generated using the median NDVI value per pixel from image sets captured within a single season. Finally, an aspect mask was built F I G U R E 3 Defining the Hindu Kush Himalaya-wide sample region, highlighting elevation bands explored in the rest of the paper [Colour figure can be viewed at wileyonlinelibrary.com] by computing aspect based on the SRTM global digital elevation model (DEM) at 30 m resolution, but the analysis was rescaled to 90 m resolution during the pixel counting phase (to calculate green fraction, described below). We masked all values <45° and >315° to remove north facing slopes that were likely to be shaded by neighbouring topography (i.e. causing the most significant bidirectional reflectance distribution function [BRDF] effects) and these were excluded from all subsequent analyses. The aspect and elevations outside the height band analysed (using SRTM DEM data) were masked before summing the total number of unmasked pixels as well as the number of unmasked pixels with NDVI values above a 0.1 threshold and deriving the 'green fraction' per year by dividing the number of pixels with an NDVI > 0.1 by the total number. An NDVI > 0.1 threshold was chosen to indicate the presence of vegetation-set at a low level due to the potential spatial sparsity of high altitude vegetation ( Figure 1).
The same procedure was used for generating green fractions per year for both regional analyses (P140-R40/41, Nepal national extent) and to ROIs over the entire HKH.

| Considerations surrounding Landsat 8 OLI to Landsat 7 ETM+ correction
Landsat 5 and Landsat 7 both use a 'thematic mapper' type sensor with very similar band positions and widths. Claverie, Vermote, Franch and Masek (2015) showed that Landsat SR products as used here were comparable between mission 5 and 7, with the exception of the blue band (not used in NDVI derivation). In contrast, the OLI sensor on board Landsat 8 differs from these previous missions: it still offers a 30 m spatial resolution, but the bands are spectrally narrower and cover different spectral ranges (Roy et al., 2016). To be more specific, blue, green and red Landsat 8 OLI band spectral response functions intersect with 82.76%, 71.08% and 60.63% of the corresponding Landsat 7 ETM+ band spectral response functions (Roy et al., 2016). Critically, for derivation of NDVI, which uses a near infra-red (NIR) and red band, the L8 OLI's NIR band (centred on 850 nm) misses a water absorption feature that occurs in the corresponding ETM+ NIR band (Irons, Dwyer, & Barsi, 2012;Roy et al., 2016). These differences can impact the quality of comparisons between L7 and L8: Holden and Woodcock (2016) report that Landsat 8 is 'consistently darker in red, green and blue bands' compared to Landsat 7, and Roy et al. (2016) report that OLI NDVI was positively biased relative to ETM+ NDVI. Biases reported by Roy et al. (2016) were lowest with NDVI derived from the Landsat SR product, as used in our study. They report an average difference of +0.0165 (4.86%) between OLI and ETM+ derived NDVI, with minimal biases at NDVI values around the 0.1 threshold used here.
In GEE, a further consideration is that the Landsat 7 and 8 SR products available use slightly different algorithms for atmospheric correction and calibration to SR: Landsat 8 uses a scheme called 'land surface reflectance code (LaSRC)' (described in Li, Roy, Zhang, Vermote, & Huang, 2019;Vermote, Justice, Claverie, & Franch, 2016) whilst Landsat 7 data are calibrated using the 'landsat ecosystem disturbance adaptive processing system (LEDAPS)' scheme (Masek et al., 2006). Both LaSRC and LEDAPS use the same 6SV radiative transfer code (Vermote, Tanre, Deuze, Herman, & Morcette, 1997) with the primary difference being in the way that the aerosol content in the atmosphere is retrieved.
For robustness, we applied the L7/L8 correction method described by Roy et al. (2016) to our data prior to statistical analysis.
This calibration is robust since it was developed using 59 million corresponding sensor observations extracted from 6,317 Landsat-7 ETM+ and Landsat-8 OLI images over the United States in winter and summer seasons (Roy et al., 2016). Alpine and high altitude vegetation with similar traits to the plants found in the HKH will be represented in this calibration model owing to the presence of high mountain systems across the conterminous United States. Fassnacht et al. (2019) have also applied the Roy et al. (2016)  Remaining uncertainties which will not be accounted for by this approach would be BRDF and topographic effects not removed by the aspect masking (see Section 2.3.2), and aerosol differences not removed by the SR quality masking. Also, if pixels contained nonvegetated dry terrain, the aerosol retrieval for LEDAPS and LaSRC could be quite different, which would manifest in blue, green and red band LEDAPS (Landsat 7) versus LaSRC (Landsat 8) differences (D. Roy, personal communication, 2019). To test the extent to which the latter effect impacted on our data, we used a glacial lake outburst flood scar near Panboche Nepal, which has remained unvegetated since the flood occurred in 1977 (Buchroithner, Jentsch, & Wanivenhaus, 1982) and is thus suitable for testing bare ground-  Table S1) so from hereon, we assumed bare ground-related biases to be of minimal concern.

| Statistical analysis
Time-series trend analysis was undertaken using OLS regression this manuscript since they carry different information. The Kendall's τ parameter conveys information about the time-series pattern, independently of the magnitude of change, whilst OLS coefficients allow trends to be modelled.

| Validation
The lack of ecological studies across the subnival HKH means that there were no existing spatially distributed ecological data sets available for validation. However, there is a relatively new, rich re-

| What is the extent of the subnival zone?
Results showed that an average area of 27,388.4 km 2 of Nepal lies above 4,150 m a.s.l. (Figure 4; Table 1). According to the RGI (Pfeffer et al., 2014) and MODIS analysis, 4,520 km 2 of Nepal is covered by permanent ice and snow, resulting in a ratio of seasonal snow to permanent snow/ice cover of 5.1:1. In the P140-R40/41 region, the ratio is higher at 15.1:1 (Table 1).  Table S3. The higher fractional vegetation cover in the H2 height band compared to the H1 height band ( Figure 5) is somewhat counterintuitive but closer investigation revealed that this was caused by a local feature within H1 which contained a large meandering riverbed with little vegetation. Table 2 gives the OLS and Mann-Kendall trend test results for these time-series trends. At all height bands for the P140-R40/41 region, we found positive, significant (p = .05) trends using data corrected according to Roy et al. (2016), with strongest slopes in H1 (OLS slope = 0.0066) and H3 (OLS slope = 0.0067). Trends at the highest elevations (H4) were very weak (OLS slope = 0.002) compared to stronger trends in lower height bands. The non-parametric

| Nepal
We repeated the analysis across the entire spatial extent of Nepal ( Figure 6). The OLS slope coefficients were smaller than those found for the P140-R40/41 region but trends in all four height bands were weakly positive. The only significant, positive trend was found in H1 and H2 using OLS (Table 3). Contrastingly, the Mann-Kendall trend analysis found significant positive trends in H1, H2 and H4. The value for Kendall's τ was quite similar for H1, H2 and H4 (0.355, 0.357, 0.344 respectively)-suggesting proportional rates of change for these height regions (Table 3). The outlier shown in 1997 ( Figure 6) probably affected OLS analyses for H3 and H4. This outlier was caused by missing data for roughly one-third of the total area due to missing scenes and cloud/snow cover, and it likely influenced the computation of the green pixel fraction. It only affects data in the higher height bands and  is not noticeable in H2 and H1, explained by more snow and cloud at higher elevations.

| HKH region
For the broader HKH region, the overall pattern of mean green fraction over time show weakly positive time-series trends ( Figure 7; Table 4). As with the more localized investigations, these trends exhibit differences with height with the highest OLS slope coefficients measured in H2 and H3 (0.0061 and 0.008 respectively; H1 shows the lowest percentage of significant positive trends (OLS = 56%; Mann-Kendall = 58%). The Mann-Kendall trend analysis found a higher proportion of significant trends compared to OLS in H1 and H2, but less in H4.
To provide an independent test of time-series trends where sensor configuration was constant, for the HKH-wide data set, we also calculated the time-series trend using SR data from Landsat 5 and 7 TM and ETM+ sensors only (i.e. excluding L8 OLI data).
Results are shared in the Supplementary Information ( Figure S5; Table S4). We found similar positive time-series trends to the full analysis shown in Figure 7, although these trends were weaker when L8 data were excluded ( Figure S5) as compared to when they were included (Figure 7). OLS slope coefficients were positive, and Kendall's τ similarly showed the highest value in H3 (Table S4) ; Table S4).

| Geographic effects: Slope and aspect
A final analysis determined whether there were differences in the slope or direction of the time-series vegetation response given different underlying geographic conditions-specifically slope and aspect. This analysis was performed with the HKH-wide data set. First, slope angle was queried within the four height bands (Figure 8). Whilst there is a lot of scatter, the slope of green fraction trends is more positive for steeper terrain slopes in H1 (Figure 8a) whilst for H2-H4, the patterns with slope trend slightly more towards negative behaviour as elevation increases (Figure 8a through d). Analysis with aspect showed no visually discernible or statistically meaningful patterns, so we have not included a plot.

| Snow fractional cover change over time
The time series of snow fractional cover showed no significant trend when analysed with the Mann-Kendall rank test for the P140/R40-41 region and Nepal extent (see Figure S3). For the 100 ROIs distributed HKH-wide, we found an increasing number of ROIs that showed weak, yet significant, negative trends in snow fraction, with increasing height.
In H1, 7% of ROIs showed a significant negative snow fraction trend over the time series, increasing to 23% in H4 (see Figure S4; Table S2).

F I G U R E 8
Green fraction  plotted against the median terrain slope value per ROI for the L5 + L7 + L8 collection over the HKH-wide extent within four height bands (

| Validation
We evaluated the utility of photographic data for providing presence/absence information about vegetation in the subnival HKH.
We report that photographs available through Google's 'Streetview' and 'Photosphere' products, and our own field photographs, contained useful ecological information about the presence or absence of broad vegetation groups (grasses, mosses, shrubs and bare ground). Figure 9 shows

TA B L E 5
Confusion matrix for the reference and classified vegetated and unvegetated pixels commission for the two classes are provided in Table 6. In the pixels falsely classified as unvegetated (n = 13), 53.8% contained shrubs, 30.8% grasses and 38.5% mosses.

| D ISCUSS I ON
A wealth of scientific work has highlighted that the HKH region is changing in response to climate (Keenan & Riley, 2018;Pachauri et al., 2014;Schickhoff, 2011), and whilst great research efforts have been invested in high altitude cryospheric and climatological work (Bajracharya et al., 2015;Gardner et al., 2013;Gautam, Hsu, Lau, Tsay, & Kafatos, 2009;Sano, Ramesh, Sheshshayee, & Sukumar, 2012;Shrestha, Wake, Mayewski, & Dibb, 1999), studies that address high altitude ecological responses are comparatively lacking (Gaire et al., 2014). Our work has advanced ecological understanding of the poorly studied HKH subnival system by being the first to measure the extent of changes in the fractional cover of subnival vegetation in the HKH (4,150-6,000 m a.s.l.). Specifically, we have shown that significant changes have occurred in vegetation distribution measured through 'green fraction' analyses of Landsat 5, 7 and 8 data since 1993. Our analyses have incorporated a robust cross-sensor calibration to correct for sensor differences between the Landsat 7 ETM+ sensor and the Landsat 8 OLI, which have different radiometric sensitivities. Taking a fully precautionary approach and excluding data from L8 OLI, across the HKH extent, we still found significant positive trends in vegetation fractional cover over time, although those trends were weaker without, than with L8 data included in the time series ( Figure S5; Table S4). We have shown that the satellite-derived green fraction product can be validated using freely available photographic data with good accuracy (79.6%; Tables 5 and 6).

| What is the extent of the subnival zone?
Subnival habitat in the Nepalese and P140/R40-41 regions, which are microcosms of the broader HKH system, covers between five and 15 times the area of permanent glaciers and snow (Table 1). These ratios indicate the extent of seasonal snow cover that could be impacted by future changes in vegetation cover. There is a much larger proportion of the high altitude Qinghai-Tibetan plateau (QTP) within the region of the P140-R40/41 scene (15.1:1) than in Nepal, which exhibits a smaller ratio (5.1:1); and the QTP exhibits greater seasonal snow storage than the steep mountainous areas of Nepal. Across the HKH region generally, it is reasonable to assume that subnival snow to permanent snow/ice spatial coverage will likely vary between the Nepalese and P140-R40/41 regional ratios: areas on the plateau exhibiting higher ratios whilst areas with steeper topography will have lower ratios. Of importance, beyond areal estimates are snow depth and snow water equivalence which impact hydrological regimes: we could not evaluate those with the satellite data used here, but further studies should seek to advance such understanding. Regardless of geographical location, these results indicate that the subnival zone should be recognized for its likely role in influencing non-base flow water supplies, with impacts on downstream communities. The lack of science questioning the role of subnival vegetation in HKH ecohydrology, with impacts scaling from headwater catchments to the spatial extent of the 'third pole' (http://www.icimod.org/?q=3487) is a gap that requires urgent scientific attention.
Whilst rapidly changing Arctic ecosystems have received considerable investment in ecological research, the same cannot be said for Himalayan ecosystems (Gurung & Bajracharya, 2012), despite the much larger number of people reliant on mountain water supplies from the HKH than in the Arctic (Immerzeel et al., 2010;Pritchard, 2019).
Yet, there is a potential nexus of impacts (particularly hydrological impacts) that could arise from vegetation expansion at high elevation. snow melt rates Pomeroy et al., 2006;Wookey et al., 2009), since Arctic work has shown that vegetation changes can impact these processes profoundly. Readers should note that Arctic work focuses strongly on changes in species/functional composition rather than the conversion of bare ground to vegetated as explored here, but we argue that the latter might be expected to deliver more profound hydrological impacts. Importantly, work emerging from Tibet shows a contrasting picture to Arctic plant/snow feedbacks, with evapotranspiration-driven cooling outweighing albedo-driven warming in shrub-dominated systems (Shen et al., 2015). Future in situ work that uncovers the functional relationships between Himalayan subnival plants and hydrology would provide timely answers to these questions, whilst modelling approaches which test the eco-hydrological impacts of vegetation expansion will deliver new insights in the HKH (e.g. following work such as Rasouli, Pomeroy, and Whitfield (2019)). We also suggest that future remote sensing work that delivers spatial information about plant traits (e.g. grass/shrub land cover mapping) over large areas could deliver new understanding about the wider mechanisms and impacts of vegetation expansion.

| Has the spatial extent of subnival vegetation changed and, if so, at what rate and where?
Using NDVI derived from Landsat data records since 1993 to distinguish between vegetated and non-vegetated pixels, we showed an increase in green fractional cover over time. Exploring these trends TA B L E 6 Omission and commission error percentages for the vegetated and unvegetated pixels, including the number of reference observations for the two classes

Vegetated Unvegetated
Reference n 78 20 Error of omission 16.7% 35.0% Error of commission 9.7% 50.0% over three spatial extents and within four elevation bands, we found that at all scales and elevations, positive trends in vegetation were found between 1993 and 2018, but these were not always statistically significant. In summary: 1. An area centred on Mount Everest (P140/R40-41) showed significant positive trends in all four height bands (OLS and Kendall's τ; Table 2); 2. At the national scale of Nepal vegetation trends in H1 and H2 (OLS and Kendall's τ) and H4 (Kendall's τ only) showed significant positive trends (Table 3); 3. Across the HKH extent, using 100 randomly distributed ROIs to sample time series all height bands between 56% (H1; OLS) and 76% (H3; OLS and Kendall's τ) of ROIs showed positive trends significantly different from 0 (p < .05; Table 4).
From this work, we conclude that subnival vegetation expansion is occurring HKH-wide, although there are clearly geographic differences in the strength and significance of that trend. This is to be expected given the topographic and climatic diversity of the HKH (Bookhagen & Burbank, 2006). The observed trends are consistent with modelling work showing a decline in temperature-limited areas resulting from warming of the Earth's cold regions (Keenan & Riley, 2018), and these results also corroborate the work of Gonzalez, Neilson, Lenihan, and Drapek (2010)  and H3 height bands cover most of the QTP area for which models predict an increase in shrub vegetation at the cost of alpine meadow, steppe and desert (Zhao, Wu, Yin, & Yin, 2011). It is possible that the observed positive trends are linked to this regime shift, also reported by Brandt et al. (2013); however, it is unclear whether the NDVI threshold-based analysis is sensitive to a change to a different vegetation type.
Trends were weakest in the highest elevation band, H4 (5,500-6,000 m a.s.l.) across all three scales of our analysis. We hypothesize that this was due to temperature limitation and overall steeper slopes which are less suitable for vegetation colonization due to regular disturbance from peri-or paraglacial activity (e.g. rockfalls and landslides; Ives & Messerli, 1989). Where increases were found to occur in H4, these were associated with gentler slopes (Figure 8d), as these may be less susceptible to erosion and plant colonization.
The opposite was found to be true for areas at lower elevations (e.g. H1; Figure 8a)-with higher green fraction model slopes on steeper terrain slopes. We suggest that this is a function of lower lying areas having high vegetation cover at the start of the time series, and thus, vegetation expansion is more likely to occur on steeper terrain at lower elevations.

| Critical evaluation of Landsat time-series data
The linear model-based Roy et al. (2016) correction is a simple approach to resolve the L7-L8 sensor transition and has been used by others in similar settings (Fassnacht et al., 2019). As shown in

| Validation
We to obtain a large sample size with which to evaluate the accuracy of Landsat products for mapping bare ground. Also, the positional uncertainty of these geotagged photos is not known, which is why we used a 100 m buffer to eliminate potential replicate observations. Validating the 2017 NDVI thresholded product (which contained observations from both L7 ETM+ and L8 OLI), we found that 79.6% of pixels were accurately classified, highlighting at least for the recent data, some of the uncertainties in the method. We assert that online repositories of photographs contain a rich source of archival data for answering basic ecological questions or performing validation in hard-to-reach areas, and yet they remain almost completely untapped as a source of data for poorly instrumented regions of the world.

| Google Earth Engine
This work evidences the power of using cloud-based computation of remote sensing analyses. To undertake the analyses described in this manuscript using traditional remote sensing workflows would have necessitated the use of a very powerful computer and required the downloading of thousands of individual Landsat scenes. For ecologists interested in studying large-scale ecosystem processes over extensive spatial areas and over multi-decadal timescales, and for those working in places without access to supercomputing facilities, the capabilities offered by GEE are unparalleled. We used Earth Engine v0.1.201 and the Javascript Application Programme Interface for all processing. We share our GEE code openly in the spirit of science as a transparent endeavour. A link to the code is provided in the Supplementary Information.

| SUMMARY
The HKH are Asia's water towers (Immerzeel et al., 2010). It is accepted that water supplies in the HKH are threatened as a result of climatic shifts (Bolch et al., 2019;Pritchard, 2019;Shannon et al., 2019), and that plants and the hydrological cycle are coupled (Fatichi et al., 2016).
Despite this, science has not so far questioned the extent to which high altitude ecological shifts or transitional processes across the HKH will impact water or carbon cycle processes. A first step towards understanding the role of vegetation in HKH ecohydrology is to measure the extent of the ecosystem relative to snow and ice cover, and to determine whether the ecosystem has changed over time. Our results show vegetation expansion is occurring at high altitudes (>4,150 m a.s.l.) across the HKH and that subnival systems cover between five and 15 times the area of permanent ice and snow. Supported by this evidence, we argue that subnival ecological systems play an important role in HKH hydrology and their role will increase as snowlines ascend and glaciers melt. There is an urgent need for new science to uncover the status, role and fate of high altitude ecosystems in the unique setting of the HKH for modulating seasonal non-base flow water supplies and cycling carbon. Further scientific work in this inaccessible region will undoubtedly rely heavily on the Earth observation data analysis, as we have performed here, but there is an accompanying, pressing need for new in situ studies. Overall, we need an improved functional understanding of the state and fate of subnival ecosystems and their ecohydrology and such data would also prove valuable for validating satellite-based analyses.

ACK N OWLED G EM ENTS
We are grateful to the following people who contributed data to our crowdsourcing methodology to identify vegetation types