Grassland Growth in Response to Climate Variability in the Upper Indus Basin, Pakistan

Grasslands in the upper Indus basin provide a resource base for nomadic livestock grazing which is one of the major traditional livelihood practices in the area. The study presents climate patterns, grassland phenology, productivity and spatio-temporal climate controls on grassland growth using satellite data over the upper Indus basin of the Himalayan region, Pakistan. Phenology and productivity metrics of the grasses were estimated using a combination of derivative and threshold methods applied on fitted seasonal vegetation indices data over the period of 2001–2011. Satellite based rainfall and land surface temperature data are considered as representative explanatory variables to climate variability. The results showed distinct phenology and productivity patterns across four bioclimatic regions: (i) humid subtropical region (HSR)—late start and early end of season with short length of season and low productivity (ii) temperate region (TR)—early start and late end of season with higher length of season and moderate productivity (iii) sub alpine region (SAR)—late start and late end of season with very high length of season and the most 698 productive grasses, and (iv) alpine region (AR)—late start and early end of season with small length of season and least productive grasses. Grassland productivity is constrained by temperature in the alpine region and by rainfall in the humid subtropical region. Spring temperature, winter and summer rainfall has shown significant and varied impact on phenology across different altitudes. The productivity is being influenced by summer and annual rainfall in humid subtropical regions, spring temperature in alpine and sub-alpine regions and both temperature and rainfall are contributing in temperate regions. The results revealing a strong relationship between grassland dynamics and climate variability put forth strong signals for drawing more scientific management of rangelands in the area.


Introduction
Grasslands occupy around half of the Hindu Kush Himalayan region, which provides a resource base for nomadic livestock grazing as well as associated ecosystem services such as maintaining biodiversity and avoiding the erosion that might result from cropping systems.The quality of grasslands is declining, primarily due to human-induced modifications such as agriculture, extensive livestock grazing, and invasive species, and they are further threatened by the potential effects of a changing climate.The high altitude grasslands in mountainous areas including the Himalayas are thought to be especially critical in terms of susceptibility to climate change [1,2].The upper Indus basin (UIB) is a particularly complex area of the Hindu Kush Himalayas because of the high climatic variability on a relatively small scale [3].The pastoralists must adapt to the changes in climatic conditions to ensure sustainable grazing plans, but they need support including providing information on vegetation and climate patterns in the area.Unless the local dynamics of the grassland ecosystem are understood and used as guidelines for development policies, interventions are likely to be random activities that comprise development by trial and error [4].There are few studies on rangeland conditions and their relationship with climate factors and change in the Himalayas and the studies that have been done were mostly qualitative and based on focus group discussions and interviews [5].A systematic assessment is needed to understand the grasslands dynamics in the region as a basis for devising appropriate policies for management.
Globally, a significant amount of work has been done on the use of satellite-based products for assessment of vegetation conditions and their related dynamics.Considering its vastness and the inaccessibility of the Hindu Kush Himalayan region, satellite-based measurements offer the only viable methodology for making assessments over anything but very small areas.With reference to vegetation assessment studies, NDVI (normalized difference vegetation index) derived from remote sensing instruments is used as a proxy to measure vegetation biomass and associated dynamics [6,7] Annual production and net carbon dioxide flux are generally assumed to relate directly to NDVI integrated over the growing season [8] and NDVI has shown consistent correlation with vegetation biomass and dynamics in various ecosystems worldwide (e.g., [9][10][11][12][13][14][15][16] Several studies have used remote sensing along with process models to analyze patterns and timing of anomalies in net primary productivity [17,18]. In the nomadic livestock grazing (transhumance) system, the migration cycle is primarily dictated by the local climatic conditions and seasonal forage availability.Phenological parameters such as onset of greenness and length of growing season can be affected by short-term climate fluctuations such as temperature and rainfall, as well as anthropogenic activities, such as groundwater extraction and urbanization.In the longer term, annual phenology may change as a result of climate change and large scale anthropogenic disturbance [18][19][20].Thus, differentiation of annual, inter-annual, and long-term phenological patterns is an important component in monitoring and modelling of global ecosystems [21] and may lead to better understanding of how vegetation cover changes over time [16].
In this context, our research aims to provide an understanding on ongoing climatic conditions and their impact on the grassland growth dynamics in the UIB.The specific questions addressed are: (1) what are the significant spatial patterns and temporal variability of grassland phenological calendars and productivity across the basin?(2) how are the spatial and temporal patterns of grassland calendars and production controlled as a function of climate?

Study Area
Geographically, the study area lies between 31°N to 37°N and 69°E to 78°E with a total area of ~17 million•ha, and is comprised of two provinces of Pakistan: Khyber Pakhtunkhwa (KP) and Gilgit-Balitistan (GB) (Figure 1).The area has diverse topography with alluvial plains, low hills, and high mountains of HKH.The vertical profile of the area ranges from 125 m in the south to about 8000 m in the north.The northern part of the UIB has a typically continental steppe climate, and precipitation is concentrated in winter and spring [22].At high elevations in the Hindu Kush, snowfall is much heavier and consequently large glaciers are a prominent feature of the landscape.Temperature in the valley areas varies between 30 °C in July to as low as 0 °C in January.In south, the climate becomes more typical of the Indian subcontinent, although a considerable proportion of the annual precipitation still comes from frontal cloud zones during the winter months.The combination of a short but powerful summer monsoon with frequent winter cloud zones gives a bimodal rainfall regime.In accordance with ecological regions in Pakistan [23], major grassland regimes in the study area are characterized with four bioclimatic regions (Table S1).

Methods Summary
A brief flow chart of methods are illustrated in Figure 2 (for details please see supplementary material).Normalised Difference Vegetation Index (NDVI) data, from moderate resolution sensor (MODIS) images taken at 16-day intervals and 250 m spatial resolution provided in a standard product (MOD13Q1, Level 3 Product), were used to extract phenology and productivity metrics (PPMs) of the grasslands: start of growing season (SGS), end of growing season (EGS) and length of growing season (LGS), seasonally integrated NDVI (SIN), and maximum seasonal NDVI (MSN) [24] (Table S2).The selected metrics were chosen to reveal spatio-temporal response patterns of grasslands to climatic variables [11,[25][26][27][28]. Grassland was masked to avoid the contribution of numerous vegetation types present in the area.In the region, altitude is considered as a proxy to micro-climatic zones-therefore, altitudinal gradient was stratified into 10 classes of 500 m elevation intervals (ranges from <500 m to >4500 m) to minimize the effect of spatial variance.Land Surface Temperature (LST), from the MODIS MOD11C3 product, which provides monthly estimates of LST with a spatial resolution of 0.05° × 0.05° (~5.6 km at the equator), and Tropical Rainfall Measurement Mission (TRMM) data (3B43 product) were used to derive ten the climate variables (CVs)-seasonal and annual temperature, and rainfall.The PPMs and CVs were aggregated to give a spatial average per elevation zone [24,29].The year-to-year response of PPMs to CVs was evaluated through the regression models [29].For summarizing the results and discussion, the elevation zones were broadly categorized into four bioclimatic regions-that are: Humid subtropical region (HSR), Temperate region (TR), Subalpine region (SAR), Alpine region (AR).

Relationship between Phenology Metrics and Temperature
Predominantly, the topography of the area is characterized with latitudinal and elevation gradient.The rainfall and temperature had large variability among season and across the elevation gradient.Valleys in mountainous areas were warmer ranging between 15 °C to 30 °C as compared to 0 °C to −15 °C on the mountains.Rainfall was relatively low in all elevation zones in autumn; in the humid subtropical region, it was highest in summer, whereas in the upper part of the temperate region it was lower in summer than in winter or spring; and in the sub alpine and alpine regions, it was highest in spring.Generally, temperature exhibited high spatial variability and low temporal variability while rainfall exhibited low spatial variability and high temporal variability (Tables 1 and 2, Figure 3).The Humid Sub tropical zone was found with high rainfall and high temperature and high rainfall variability consisting of three seasons.On a relative scale, the temperate zone has medium rainfall and temperature with less variability of rainfall and temperature between peak and lean seasons.The alpine zone was found contrastingly with high rainfall and low temperature but with high degree of variation in Min-Max rainfall and temperature.This characteristic climatic variability delineated using satellite data across the elevation regimes broadly agrees with reported bioclimatic thresholds [23].Functional relationship between NDVI-based phenology and productivity metrics and the ten climate variables has been examined to elaborate inter-annual variation in vegetation dynamics across the ten elevation zones.In understanding physical meanings of these linear correlations, it is important to consider consequent classes, coefficient of correlation (R), and the corresponding CVs.In addition, plant responses to climate variables may also be specific.Statistically significant correlations (p < 0.05) derived from linear regression analyses of phenology and productivity metrics with climate variables are shown in (Tables S3-S7).
The start of growing season showed a significant positive correlation with spring temperature in the humid subtropical region and winter temperature in the temperate region, and a negative correlation with spring temperature in the sub alpine and alpine regions (Table S3).The end of growing season only showed one significant correlation-with summer temperature in the lower zone of the humid subtropical region (Table S4).The lack of significant relationships especially over EGS could be due to the seasonal averages used for the relationship.We expect that the seasonal minimas of temperature and rainfall, which explain temperature and moisture stress, could have resulted in better explanation of EGS.The length of growing season showed a significant negative correlation with summer temperature in the humid subtropical region, winter temperature in the temperate region-while a positive correlation is observed with annual temperature in the humid subtropical region, temperate region, subalpine, and alpine regions and also with spring temperature in the subalpine and alpine regions.(Table S5).
The maximum seasonal NDVI showed a significant negative correlation with spring and summer temperature in the temperate region and annual temperature in the humid subtropical region and temperate regions (Table S6).The seasonally integrated NDVI showed a significant negative correlation with winter, spring, summer, and annual temperatures in the temperate zone, and a significant positive

SON
correlation with spring and annual temperature in the subalpine and alpine regions, and summer temperature in the higher alpine zone (Table S7).

Relationship between Phenology Metrics and Rainfall
There was only a small amount of correlation between growing season and rainfall.The start of growing season showed a significant positive correlation with winter rainfall and negative correlation with spring rainfall in the humid subtropical while in the subalpine region a positive correlation with winter rainfall (Table S7).The end of growing season only showed significant positive correlation with summer rainfall in the subalpine region (Table S4).The length of growing season showed significant positive correlation with summer rainfall in the subalpine region and negative correlation with winter rainfall in the alpine region (Table S5).
The maximum seasonal NDVI showed a significant positive correlation with winter rainfall in the lower part of the temperate region, with spring rainfall in the temperate region and upper part of the humid subtropical region, and with summer and annual rainfall humid subtropical region.(Table S6).The seasonally integrated NDVI showed a similar pattern to the maximum seasonal NDVI with the addition of a significant positive correlation with summer rainfall in the subalpine region (Table S7.).There was no significant correlation between productivity variables and rainfall in the alpine region.
Collectively, it can be deduced from above results, that seasonal productivity of the most productive grassland communities, zone 7 and zone 8, was influenced by EGS dependent upon JJAR, due to which lengthening or shortening of LGS tended to increase or decrease in Seasonal Integrated NDVI.. Similarly, the negative relation of SIN and LGS with DJFR in high alpine zones (Zone 9 and Zone 10) was unexplained as both EGS and SGS, and MSN do not portray any relation with DJFR.It was interesting to note that the seasonal productivity (SIN) of low (<1000 m) and mid-altitude (1000 m to 3000 m) grassland communities is more dependent upon increase in MSN due to rainfall event rather the variation in LGS.

Phenology and Productivity Patterns
The 11-year means of the phenology metrics (including SGS, EGS, LGS) and productivity metrics (including MSN, SIN) were analysed across the elevation zones (Figure 4).
The earliest growing season found to start in the early days of MAM at TR (zone-3 to 6), in late days of MAM in HSR (zone-1 & 2) and then early days of JJA in high elevations (zone-7 to 10).The HSR was largely influenced by monsoon controlled rainfall and spring temperature regimes.The SGS was found retreating (late start) with rise in spring temperature (MAM) and advancing (early start) with increase in winter and spring rainfall.General late start of season observed in the HSR can be attributed to much less rainfall during the spring season in this region (Tables 1 and 2) while vegetation growth is dependent on the moisture availability.In the rest of the regions, the SGS was found steadily shifting to later parts of the season with increase in elevation at higher elevation regions.The shift ranged from 144th day to 178th day with change in elevation from zone-7 to zone-8.The results show that SGS has different patterns across the bioclimatic regions.
The EGS ranges from August (240 DOY) to November (320 DOY).The difference in EGS by 80 days across elevation ranges clearly shows high variability compared to 64 days of SGS.The earliest and the latest end of the season occur at TR (zone-3) and highest elevations of AR (zone-10).The mid elevations, TR, which have not exhibited variation in SGS were found with increasing trend of EGST with rise of elevation.The range of variation was found by 27 days across the elevation zone-3 to zone-6.The higher elevations, SAR and AR, did not exhibit variation in EGST with increase in elevation as was found in case of SGS.The analysis shows that low and mid elevations have more EGS variability than high elevations.The LGS was found with a range of 123-173 days across the bioclimatic regions.The increase in LGS was observed with rise in elevation till 3500 m (SAR) and then decline at further elevation in Alpine Regions.The maximum range of 150-173 days of LGS was found between elevation zone-6 to zone-7 (SAR) and lowest range of 120-123 days was found in lower elevations (HSR) where late SGS was observed due to limited moisture and early EGS was found due to rise in in-season temperature (JJA).The LGS was positively related with annual temperature (ANNT), except temperate region, consequently contributing in overall grass productivity across the regions.Overall, in terms of bioclimatic regions, the humid sub-tropical region has a late start, early end, and short growing season with low productivity; the temperate region has an early start, early end with moderate productivity that increases with elevation; the sub-alpine region has a late start, late end, and long growing season, with highest productivity; and the alpine region has a late start, early end, and short growing season with low productivity.The grass calendars as observed over the study areas broadly matches with ground-based reports [30,31].
In terms of variability across time (11 years) and space (within the bioclimatic regions), the SGS exhibits high spatial and temporal variability in the humid subtropical and alpine regions and low variability in the temperate and subalpine regions.The end of growing season had greater variability both in terms of days and in differences between bioclimatic regions, with a high temporal variability in the humid subtropical region and high spatial variability in all except the alpine region.The length of growing season also showed different spatial and temporal variability in the different bioclimatic regions, with a higher temporal variability in the subalpine and alpine regions, higher spatial variability in the temperate region, and low temporal and spatial variability in the humid subtropical region (Figure 5).The increase of SGS and decrease of EGS observed over different elevation regions of alpine zones and temperate zones respectively and also increase of LGS with rise in elevation reflects the possible control of climate across different elevations and growth stages of the grasses.Such characteristic variations of across elevation regions affecting grass calendars were also reported in the neighboring areas [28,32].Seasonal maximum productivity, MSN, ranges from 0.15 (in zone-10) to 0.43 (in zone-7), whereas, integrated productivity ranges from 1.30 (in zone-1) to 3.25 (in zone-7).It was notable that LGS in zone-10 was higher than zone-1 but MSN in zone-10 was lower than that of zone-1, therefore, as a by product of MSN and LGS, SIN in elevation zone-1 was relatively higher than elevation zone-10.Analysis indicates that, HSR and AR were characterised with low productivity, TR were moderately productive, and SAR were highly productive.The maximum seasonal NDVI and seasonally integrated NDVI showed a similar pattern of change with elevation zone.Productivity increased with increase in elevation up to a maximum of 0.43 in sub-alpine region and decreased rapidly with elevation above 4000 masl (Alpine region) to a minimum of 0.15 at elevations above 4500 masl.
The maximum NDVI showed the highest temporal variability and lowest spatial variability in the humid subtropical region; temporal variability tended to decrease with elevation and spatial variability to increase, except at the highest elevation.Seasonally integrated NDVI showed low spatial and temporal variability in the humid subtropical region and in the highest zone of the alpine region; temporal variability was only high in the subalpine region, whereas spatial variability was high in the temperate, sub-alpine, and lower zone of the alpine region.

Discussion
The study was undertaken to support understanding of the climate controls on grassland phenology and productivity in the UIB basin using satellite-based information in absence of ground-based climate and vegetation data.The use of satellite based climatology data was helpful in understanding spatial patterns of seasonal and annual differences in rainfall (Rf) and temperature (T) in four bioclimatic regions-humid subtropical, temperate, subalpine, and alpine regions.Each zone showed diverse patterns of rainfall and temperature, the humid subtropical region had high rainfall, high temperature, and high rainfall variability, in three of the four seasons.The temperate region had medium rainfall and temperature with less variability between peak and lean seasons.The subalpine region had high rainfall and low temperature, with high variation in both rainfall and temperature.The alpine region had low temperature and rainfall, with lower variability.These spatio-temporal patterns of rainfall and temperature closely match with the observation made by [33] using ground based climate data across the elevation gradient of upper Indus basin.The role of climate on grassland calendars and production was assessed using these characteristics as a base.The analysis did not show any significant long-term climate trends except for the rainfall in the humid subtropical region over the 11-year time frame.The relationship between climatic patterns and grassland phenology was assessed in terms of the influence at different stages of grass growth and in the different bioclimatic regions.Some significant relationships were identified based on 11 years of seasonal and annual climate phenology relations.In the low altitude humid subtropical region, the growing season was strongly related to preseason rainfall and temperature.The start of the growing season delayed (late) with higher preseason (spring) temperature and advanced (early) with higher preseason (winter, spring) rainfall; and the end of growing season advanced (early) with higher in-season (summer) temperature.This is in line with the fact that this region is largely influenced by monsoon rainfall and temperature regimes.The growing season started earliest overall in the temperate region instead of sub humid regions, which can be attributed to limited rainfall at the start of season and also the frequent moderate drought conditions in the sub-humid region [34].A positive slope for the regression of phenology on autumn temperatures has been reported by several correlative analyses of plant time series [35,36].The slopes identified may reflect chilling requirements, as the UniChill model assumes [37], and there is ample experimental evidence that a period of winter chilling brings budburst forward [38].
The temperate region start of season got delayed with higher preseason (winter) temperature even though grass productivity was not affected by seasonal or annual temperature or rainfall conditions in this region.This can be partly attributed to the good forest cover in the temperate zone, which reduces thermal effects and maintains soil moisture at the local scale.In the humid subtropical region, the productivity and seasonality do not have significant associations with annual temperature and annual rainfall, while it has strong association with seasonal climate conditions.The grass productivity in humid subtropical regions is moisture constrained where productivity increases with increase in summer rainfall.Similarly, seasonality of the grasses have a strong association with seasonal climate conditions where high spring temperature delays the start of season and increase in summer temperature retreats at the end of the season.The growing season in the subalpine and alpine regions starts early in case of higher preseason temperature in spring (MAMT), which is probably due to snow melt triggering the grass growth.However, contribution of snow cover dynamics has not been investigated in the current study.The end of the growing season in the subalpine region was delayed with increased summer rainfall (JJAR), and the growing season became longer consequently.The length of growing season was associated with temperature, especially annual temperature, in most zones, but in different ways in the different bioclimatic regions, with higher annual temperature associated with a longer growing season in the temperate, subalpine, and alpine regions, and with a shorter growing season in the humid subtropical region.The temperate region [39] showed that mean air temperature during late winter and spring, annual mean air temperature, and annual Growing Degree Days (GDD) totals have positive correlations with growing season duration.
No clear relationships were identified that were true for all parameters in all bioclimatic zones.This could be because the growth regime of grass is different in the different zones, and thus affected differently by changes.Equally, the failure to identify significant relationships, especially for the end of the growing season, could also be related to the use of seasonal averages to test for relationships.Seasonal minima for temperature and rainfall might be better indicators of temperature and moisture stress and show clearer relationships with the end of the growing season [40].Overall, several significant relationships identified across bioclimatic regions show that the seasonality of grass has a strong association with seasonal climate conditions, with high spring temperatures advancing the start of the growing season in the sub alpine and alpine region and high summer temperatures leading to an earlier end to the growing season in the humid subtropical region.Further, annual mean temperature has a greater influence than annual precipitation on the grassland dynamics in the study area.
Other studies on long-term trends of climate change in the study area indicate possible impacts based on the grassland phenology and productivity behaviour observed in this study.Researchers observed in [41] a reduction in mean temperature during spring and summer and, as a result, a reduction in mean annual temperature in their analysis of ground data from 1966 to 2004 in the upper Indus basin.Our analysis suggests that start of growing season delays with lower spring temperatures and also productivity in the alpine and sub alpine areas is constrained by temperature.Thus the observed changes in the climate conditions in the recent past suggests that the alpine and sub-alpine grasslands were passing through a reduction in grassland productivity and a delay in the start of growing season.However, these historic trends may not persist in the long term future scenarios [42].The climate-phenology relations observed in the current study are similar to other global and regional studies and [38] reported a link between length of growing season and temperature in Europe; they found that warming in spring by 1 °C advanced the growing season by seven days.In addition, authors in [39] reported an apparent increase in the mean annual growing season in Europe by 10.8 days from an analysis of satellite data from 1981 to 1991.The results of [43] also are consistent with data on the advance in the seasonal cycle (about seven days) which has been inferred from CO2 data taken between the 1960s and early 1990s [44].A study on Tibetan Plateau also observed that the start of growth season (SOG) was delayed gradually, the end of growth season (EOG) advanced slowly, and the length of growth season (LOG) shortened gradually.The study also observed that elevation played an important role in the regional differentiation of phenology.
Grass productivity in the upper Indus basin was found to be rising with elevation to a maximum in the subalpine region.The high productivity in alpine zone is attributed to both seasonal amplitude (max plant growth/photosynthesis) and the extended length of the growing season in the zone.The higher biomass in the subalpine region in spring and summer as compared to the dry temperate or humid subtropical regions, is thought to reflect the increased soil moisture made available from snow melt during spring and early summer, which supports increased plant growth [45,46].The study provides a basis for meaningful discussion on the development of local climate-vegetation models that can be used in grasslands management in the region [47].More detailed analysis of production patterns would add further value for analysing the grazing potential.The information derived from this study can be integrated to assess grassland responses to climate variability and change and to support informed decision making for intra-annual and long term grazing management.

Conclusions
The observed correlations between climate factors, phenology, and productivity metrics should be seen as a first approximation using the best available data.The seasonal climatic parameters could potentially play a role as triggering, enhancing, or prohibitive factors that affect the start and/or end of the growing season.Spring temperature is the major determinant of the start of growing season in the alpine and subalpine region.The length of growing season is mostly being governed by the annual temperature all across the bioclimatic regions.The productivity (season integrated NDVI) is being influenced by summer and annual rainfall in the humid subtropical region, spring temperature in the alpine and sub-alpine region, and both temperature and rainfall are contributing in the temperate region.
In terms of implications for the transhumance grazing system, the large area of sub-alpine and alpine grasslands with higher LGS indicates the greater grazing potential.However, a high degree of temporal and spatial variability of LGS as observed in the study cautions for better ground-based understanding of the causes of variability and associated judicious planning needs for potential grazing.On the other hand, the humid subtropical zones occupying the second largest area with shorter LGS needs to be properly planned for grazing needs.The temperate zones covering the lowest proportion of grassland area but with a greater portion of the study area, supports grasses with early EGS and high variation in LGS across different elevations.The spatial information thus generated on crop calendars could be potentially used for spatial planning of grazing within the zone and regulating the grazing across the zones so as to sustainably meet the herd movements and grazing needs.Stochastically generated future projections can be potential future work to investigate long-term sustainability of grasslands in the region.

Figure 1 .
Figure 1.Map of the Upper Indus basin (study area) showing land cover and location within the Hindu Kush Himalayan region.

Figure 2 .
Figure 2. Summary of methodology (see Supplementary 1 for detailed methods).

Figure 3 .
Figure 3. Distribution of bioclimatic zones against seasonal temperature and precipitation conditions in the study area.

Figure 4 .
Figure 4. Grassland phenology across elevation zones averaged over 11 years.Bars indicate start, end, and length of growing season-Grassland productivity patterns across the elevation zones as indicated by Maximum seasonal NDVI and seasonally integrated NDVI averaged over 11 years.

Figure 5 .
Figure 5. (A,B) Temporal and spatial variability of phenology calendar and (C,D) productivity patterns across the elevation zones.Temporal variability is the average temporal variance of the spatial mean of all pixels in each zone; spatial variability is the average spatial variance of the temporal mean of all pixels in each zone.

Table 1 .
Temperature across the elevation zones-average temperature in each zone with the spatial and temporal variability.

Table 2 .
Rainfall across the elevation zones-average temperature in each zone with the spatial and temporal variability.