Precipitation changes and their relationships with vegetation responses during 1982–2015 in Kunhar River basin, Pakistan Open Access

Vegetation is a crucial component of the terrestrial ecosystem which plays an essential role in regulating climate systems (Betts 2000; Kleidon et al. 2000; Scholes & Smart 2013; Chen et al. 2015; Zari 2017). Atmospheric water in the form of precipitation is distributed to soil water, runoff, and evapotranspiration via vegetation actions. For instance, vegetation enhances soil water by increasing infiltration opportunities as well as by impeding surface runoff. Furthermore, the accumulated soil water goes back to the atmosphere in the form of transpiration via the actions of roots, stems and leaves (Shaxson & Barber 2003). As a result, vegetation plays the central role of linking atmospheric and land surface processes. In the era of climate change, vegetation has invaluable carbon sequestration and thereby controls global warming (de Vries et al. 2009; Velasco et al. 2016; Atsbha et al. 2019). However, global vegetation diversity is seen to decrease due to farmland expansion, urbanization, and so on (Sax & Gaines 2003; Ceballos et al. 2015).

Consequently, environmental degradation has become a formidable challenge in various places, mainly in developing nations. Despite their roles in climate systems and hydrological processes, vegetation activity, in general, is considered an indicator of environmental health (Florentina & Io 2011; Salmond et al. 2016). Therefore, studying vegetation activity status enables successful assessment of productivity in agricultural and other natural environments, and supports decisions to be made towards ecosystem monitoring and adaptation strategies.

With the growing concern about climate change, researchers and ecosystem monitoring bodies need to understand, model, and predict climate change's impact on vulnerable ecosystems. To better understand the role of climate in ecosystem functioning, it is necessary to assess the relationship, coevolution, and possible feedback between climate components and vegetation. Precipitation is one such crucial climatic component that plays a significant role in vegetation growth. Understanding the linkages of precipitation with vegetation and its impacts on the overall vegetative activity is of paramount importance to prioritize adaptation strategies in favour of vulnerable vegetation to climate variabilities (Hof et al. 2017; Gallagher et al. 2019).

Extensive studies have been conducted to demonstrate the relationships between precipitation changes and vegetation responses in different parts of the world (Nightingale & Phinn 2003; Wang et al. 2003; Hawinkel et al. 2016; Chen et al. 2020). The studies report various degrees of relationship between the two under different conditions. Chen et al. 2020, for instance, reported a strong precipitation–vegetation relationship in China where mean annual precipitation of 150–500 mm and specific precipitation frequency and concentration are observed. Several other studies have shown weak or no relationships between precipitation and vegetation (Jobbágy et al. 2002; Camberlin et al. 2007).

The possibilities unlocked by remote sensing technologies enabled researchers to observe earth processes to a much larger spatial extent. Along with efficient computational ability, remote sensing products also shaped the various approaches to tackle environmental problems (Melesse et al. 2007; Young & Onoda 2017). Remotely sensed normalized difference vegetation indices (NDVI) are widely used proxies for vegetation activities (Tucker et al. 1985; Aguilar et al. 2012). They are derived from red and near infra-red reflectance retrieved from special sensors onboard various satellites (Tucker 1979). The trends in NDVI time series can be used to identify vegetation changes (Tucker et al. 2001; Anyamba & Tucker 2005; Alcaraz-Segura et al. 2010) to inform ecosystem monitoring practices (Stoms & Hargrove 2000).

Kunhar River basin of Pakistan is one of the basins experiencing rapid conversion of forested areas into shrublands and farmlands. The conversion is amplified in recent years due to high population pressure demanding fuel woods and farming spaces (Star Hydropower Limited 2010). Given the sloping and rugged mountainous nature of the basin, the clearance of forested areas and encroachment of slope farming destabilizes the natural ecosystem. As a result, loss of biodiversity, landslides and environmental degradation are becoming common in the basin. The ever-increasing global temperature (IPCC 2018) is also expected to affect the basin's upland cold region biomes and leave them more vulnerable. Therefore, quantifying the possible impacts of climate change on these fragile ecosystems of the basin helps understand the problem's extent and direct management options accordingly. Although the basin is undergoing these alarming situations, limited researches have been conducted to address the issues.

Using satellite-derived NDVI and globally available interpolated precipitation time series, this study aims to (1) characterize the spatiotemporal patterns of precipitation–vegetation trends over the Kunhar River basin, and (2) identify the temporal extent and spatial patterns of vegetation response to precipitation by considering the time-lag effect.

The Kunhar River is a 171 km long river originating from northwestern Himalaya and flowing southwards to join the Jhelum river (Akbar & Gheewala 2020). Its entire watershed area is located in Pakistan and is an important source of water to Mangla Dam. It contributes up to 11% of its flow to the Mangla reservoir (Khan et al. 2019). It is one of the most important rivers administered under the Pakistan Water and Power Development Authority in hydrological monitoring. The Kunhar River basin covers an area of around 2,600 km² with an elevation ranging from 631 m.a.s.l to 5,214 m.a.s.l (Figure 1). The topography is characterized by highly rugged mountainous features with steep slopes. The area has mild summers and cold winters. It also receives substantial rainfall during the summer and monsoon seasons (April–October). The average annual rainfall at Balakot, for instance, is 1,538.5 mm and at Muzaffarabad is 1,351.9 mm (Star Hydropower Limited 2010). In winter (December–March), the basin receives significant precipitation in the form of snow.

The dominant vegetation in the Kunhar River basin is alpine meadows and coniferous forests covering the basin's northern and southern portions (Figure 2). Snow cover and agricultural lands also cover significant portions of the basin. The extent of natural forests is significantly shrinking in recent years to satisfy fuelwood demands due to the expansion of the population towards the high hills (Schickhoff 1995). They are being lost at a substantially greater rate than regeneration. As a result, the forested areas are now being progressively replaced by a mixed grassland and shrubland community (Star Hydropower Limited 2010).

In recent years, slope agriculture is also becoming a common agricultural practice in which sloping uplands are being used for farming. The rapid changes in land cover coupled with climatic variabilities and changes will ultimately affect the basin's flora and fauna's well-being. Due to the steep slopes, landslides, erosion and other environmental disturbances are waiting to happen unless a well-informed decision has been made and wise ecosystem management has been put in place to preserve the natural environment.

Spatiotemporal studies of climatic variation usually require accurate precipitation datasets. In most cases, these precipitation observations come from in situ rain gauges. However, the data obtained from rain gauges lack the required quality and spatial density (Rana et al. 2015; Kidd et al. 2017). The problem is more severe, especially for developing countries such as Pakistan. For this reason, various alternative sources of precipitation datasets are used for environmental studies. Such datasets usually come from interpolation from better quality rain gauges, remote sensing sources or model simulations (Xie & Arkin 1996, 1997; Yang et al. 2014). In many cases, these alternative data sources have shown a better representation of precipitation phenomena in various parts of the globe at different spatiotemporal scales (Stampoulis & Anagnostou 2012; Peña-Arancibia et al. 2013; Huang et al. 2016; Abebe et al. 2020). Globally available interpolated precipitation datasets are frequently used to study environmental processes and validate precipitation products obtained from remote sensing and model simulations (Chen et al. 2002; Xie et al. 2007; Harris et al. 2014).

The precipitation dataset used in this study is acquired from the WorldClim database (https://www.worldclim.org/data/monthlywth.html#). The dataset consists of historical monthly precipitation data for 1960–2018 at a spatial resolution of 2.5 min (approximately 4.6 km).

The data are downscaled from Climatic Research Unit gridded Time Series version 4.03 (CRU-TS-4.03) by the Climatic Research Unit, University of East Anglia, using WorldClim 2.1 for bias correction (Harris et al. 2014; Fick & Hijmans 2017). The original CRU-TS-4.03 dataset is derived from in situ observations distributed all over the world. An angular distance weighting interpolation has been applied to these gauge observations to produce global gridded precipitation information. In this study, we used the data of 1982–2015 at its native spatial resolution. Figure 3(a) shows the spatial pattern of mean growing season WorldClim precipitation values for 1982–2015 over the Kunhar River basin.

Canopy reflectances and NDVI have been widely used as an indicator of vegetation activities (Sellers 1985; Myneni et al. 1995). One of the most used NDVI products in many studies is the third generation Global Inventory Monitoring and Modeling Systems (GIMMS3 g.v1). The dataset has received wide acclaim due to its high-quality and long-term temporal coverage. GIMMS3 g.v1 is derived from the advanced very high-resolution radiometer instruments onboard National Oceanic and Atmospheric Administration satellites (Tucker et al. 2005). It has a spatial resolution of around 8 km and a bi-monthly temporal resolution. Various noise corrections on sensor view angle, volcanic aerosols, satellite drift and so on have been applied to maintain its quality (Pinzón et al. 2005; Forzieri et al. 2011). Consequently, the GIMMS3 g.v1 NDVI comes with a global extent and a temporal coverage from July 1981 to December 2015.

In this study, we used GIMMS3 g.v1 (https://ecocast.arc.nasa.gov/data/pub/gimms/3g.v1/) NDVI dataset from January 1982 to December 2015. We have generated monthly NDVI datasets by combining each 15-day raster file into a monthly scale using the Maximum Value Composite technique or MVC (Holben 1986). Then we resampled each monthly data to nearly 4.6 km to match it to the resolution of the WorldClim precipitation dataset used in this study. Next, we retained growing season NDVI values (April–October) to avoid the impact of winter snow and better reflect the vegetation growth in summer and monsoon seasons. Figure 3(b) shows the spatial pattern of mean growing season NDVI values for 1982–2015 over the Kunhar River basin.

Land cover maps are crucial datasets to assess the impacts of atmospheric fluctuations over a given area. They allow the existing vegetation types to be located and the identifcation of which vegetation types are more vulnerable to climatic changes. In this study, we used the Global Land Cover dataset of the year 2000 (GLC2000). The GLC2000 was developed by an international partnership coordinated by the European Commission's Joint Research Centre (Bartholomé & Belward 2005) and is made freely accessible to users (https://forobs.jrc.ec.europa.eu/products/glc2000/products.php). The dataset is produced globally by using imageries obtained from the VEGETATION sensor onboard the SPOT4 satellite (Mayaux et al. 2006). As a procedure for the generation of GLC2000, an unsupervised classification technique has been applied on the raw imageries to identify land cover types. The dataset has a native spatial resolution of nearly 1 km and consists of six land cover classes over the Kunhar River basin (Figure 2). Although one land cover map cannot represent the land cover changes throughout our study period 1982–2000, we choose to use GLC2000 due to its high-quality nature and the year 2000 is in the middle of our analysis period.

To identify possible trends in both precipitation and NDVI over the Kunhar River basin for 1982–2015, we applied the nonparametric Mann–Kendall trend test (Mann 1945; Kendall 1955) along with the Theil–Sen trend estimation technique (Sen 1968; Theil 1992). We chose the nonparametric tests and estimation methods to eliminate the influence of outliers that might affect datasets’ statistical properties.

The Mann–Kendall trend test is widely used in hydrometeorological studies to identify monotonic trends (Rosmann et al. 2016; Forootan 2019; Yagbasan et al. 2020). It has a null hypothesis (Ho) of ‘no trend’, meaning that the data come from an independent population and are identically distributed (iid). The alternative hypothesis (Ha) supposes the data follow either an increasing or a decreasing trend over time. The test has its strengths and limitations. The test's strengths lie with its conceptual simplicity and no requirement for the underlying data to have a specific statistical distribution. As a result, the Mann–Kendall trend test can handle time series with outliers (Yue et al. 2002). The limitations of the Mann–Kendall trend test originate from its null hypothesis. It assumes the underlying data to come from a population with independent realizations and identically distributed. Therefore, rejecting the null hypothesis can only mean an absence of iid nature in the dataset. However, since there is no trend in iid datasets, we can safely apply Mann-Kendall trend tests to detect trends once the null hypothesis is rejected (Chandler & Scott 2011).

To estimate magnitudes of possible trends, we employed the Theil–Sen trend estimator (Sen 1968; Theil 1992). Theil–Sen is a robust linear regression method that chooses the median slope among all lines through pairs of sample points. The estimator is nonparametric and is capable of dealing with datasets consisting of outliers. Figure 4 shows the geographical spread of trends identified by the Mann–Kendall trend test and estimated by Theil–Sen estimator for both precipitation and NDVI time series over Kunhar River basin for 1982–2015.

Nonparametric correlation analyses also provide information on the distribution-free relationship between variables. In this study, we applied one such nonparametric correlation method called Spearman's rank correlation (Best & Roberts 1975). The Spearman's rank correlation coefficient is a nonparametric measure of correlation strength between two variables. In that sense, an arbitrary monotonic function is set to describe the relationship between the two variables without making assumptions about the underlying distribution they possess (Maritz 1995). The Spearman's rank correlation coefficient is calculated similarly as the usual parametric Pearson correlation coefficient, but the computations are based on the ranks rather than the actual data values.

The null hypothesis (Ho) of Spearman's rank test assumes ‘no correlation’, meaning one variable's rank is not covarying with the other variable's rank. In other words, as the rank of one variable increases, the rank of the other variable is not more likely to increase (or decrease). The alternative hypothesis (Ha) considers a positive or negative correlation between the two variables. The optimum value for Spearman's correlation is unity with a correlation value of zero, indicating the absence of correlation.

This study applied Spearman's rank correlation analysis at pixel scale to examine the relationship between precipitation and NDVI. We chose the Spearman's rank correlation against the other competitive Pearson's correlation mainly because of two reasons. The Spearman's correlation evaluates monotonic relationships in which the changes between the variables are not required to be constant. A constant change between the correlating variables is mandatory to apply Pearson's correlation. On the other hand, Spearman's rank correlation is a nonparametric method that is suitable for variables with possible outliers. We performed four correlation analyses between them to identify simultaneous and lagged relationships. We considered the lagged correlation analysis because the precipitation–vegetation relationship is better inferred from lagged rather than the simultaneous correlations due to strong atmospheric variability (Svoray & Karnieli 2011). Lagged correlations enable us to characterize the delayed response of vegetation to precipitation variations. Consequently, we analyzed the correlation between simultaneous monthly precipitation and NDVI pixels, and the correlation between current NDVI values and previous 1, 2 and 3 months of precipitation values.

Precipitation and vegetation distributions vary spatially depending on the various influencing factors affecting them. Characterizing their spatial variation helps guide ecosystem monitoring activities and direct the necessary attention towards certain sensitive vegetation types. The spatial distribution of mean growing season precipitation and NDVI for 1982–2015 over the Kunhar River basin is shown in Figure 3.

The mean growing season precipitation over the Kunhar River basin follows a decreasing pattern from south to north (Figure 3(a)). The downstream portion receives from 700 mm to 1,050 mm precipitation. As a result, this area is mostly dominated by temperate coniferous forests and agricultural vegetation types. The central portion of the basin receives moderate growing season precipitation between 350 and 700 mm, while the northern upstream part receives relatively low precipitation, amounting to below 350 mm.

Similarly, the growing season average NDVI distribution over the Kunhar River basin follows the general decreasing trend from the southern to the northern portion of the basin (Figure 3(b)). In the downstream part, the forested and agricultural areas show NDVI values ranging from 0.5 to 0.7. In the central portion, however, a more considerable variation in NDVI is observed. In this region, the typical NDVI values range from 0.3 to 0.6. In the basin's northern highlands, growing season NDVI falls as low as 0.2, indicating limited vegetation activity. The lower NDVI values in the north can be due to the dominant alpine meadows and snow covers. Brightly coloured flowers of the alpine meadows and the radiant snows could affect reflectances captured by the satellite sensors and, as a result, lower greenness is observed in the NDVI. The higher NDVI values to the south are also due to the green coniferous forests and agricultural vegetation.

The south to the north decreasing pattern in mean growing season precipitation is in agreement with the findings of (de Scally 1992). Their study reported that the average monsoon precipitation (May–October) in the Kunhar River basin decreases markedly in the northeast from 870 to 200 mm (de Scally 1992). This finding is almost equivalent to ours. The same south to north decreasing precipitation trend seems to guide vegetation distribution in the Kunhar River basin. As a result, as represented by NDVI, the vegetation activity also falls in the same direction.

Vegetation and precipitation conditions are usually subject to interannual variations. Certain factors such as atmospheric conditions, topography and land cover types can influence the changes’ patterns. Interannual changes in precipitation and NDVI patterns of the Kunhar River basin are shown in Figure 4. As shown in Figure 4(a), the whole basin shows an increasing precipitation trend for 1982–2015. However, the precipitation change rate is higher in the southern portions of the basin than in the northern parts. The central parts of the basin experience a moderate interannual precipitation change. Interannual growing season precipitation changes over the Kunhar River basin for the years 1982–2015 exhibit a trend as low as 0.046 mm/year in the upstream portion and as high as 0.521 mm/year in the downstream parts. The interannual precipitation changes seem to increase in the opposite direction of an increase in altitude. As a result, higher altitude regions had a smaller year-to-year precipitation fluctuation than lower altitude regions.

On the other hand, interannual NDVI changes over the Kunhar River basin for 1982–2015 show slightly decreasing and increasing trends in certain portions of the basin (Figure 4(b)). Northern portions of the basin exhibit a slightly decreasing trend, whereas southern and central parts increase. The magnitude of the interannual decreasing NDVI change reaches as low as 0.003 in the upstream portion and the increasing trend extends up to 0.003 in the central and downstream parts of the basin. The vegetation activity is also seeming to be influenced by the topography. On the one hand, high altitude regions show smaller and decreasing year-to-year changes in vegetation activity. On the other hand, low altitude regions show an increasing trend in vegetation activity.

The trends observed in NDVI over the Kunhar River basin also agree with the land cover type distributions depicted in Figure 2. Interestingly, the direction of increase in NDVI trend seems to follow the south to north encroachment of the slope agriculture. Due to population growth and the need for additional farmlands in the basin, slope agriculture is increasing northwards. As a result, vegetation activity tends to increase following this human influence.

The nature of vegetation can partly explain the relatively higher precipitation trend over the basin's southern section. The taller and highly vegetated forests could significantly contribute to the release of water vapour (necessary for cloud formation) and affect solar radiation absorption and emission (Reynolds et al. 1988). As a result, the forests play a regulatory role in the surrounding microclimate. Therefore, the interannual precipitation fluctuation remains relatively high. On the other hand, the northern highlands are usually cold and harbour short herbaceous vegetation. These dwarf meadows are physically constrained to influence the local climate condition as the taller forests do. As a result, these areas show a more or less similar year-to-year interannual precipitation change.

Although other factors such as atmospheric conditions and biogeographic features like topography and soil type play a significant role in determining vegetation responses, precipitation is critical in influencing vegetation activities. The magnitude, intensity and distribution of precipitation highly affect the associated patterns and vegetation activities. In some cases, vegetation may respond quickly to precipitation variation. This quick response especially makes sense in arid and semi-arid regions where water availability is limited. In humid areas, however, vegetation may exhibit a delayed response to precipitation fluctuations. The same is true for the Kunhar River basin. The spatiotemporal correlation between vegetation and precipitation over the Kunhar River basin for 1982–2015 is shown in Figure 5. The figure shows the relationship between NDVI and simultaneous month and the previous 1–3 months’ lagged precipitation values.

In general, it is rare to find a higher correlation (r > 0.75) between precipitation and NDVI in the Kunhar River basin. Except for lag1 (when NDVI is correlated with previous month precipitation values), all other lagged correlations yield correlation values less than 0.5. For lag0 (when NDVI is correlated with simultaneous month precipitation values), upstream portions show the lowest correlation (r < 0). Central portions of the basin exhibit poor correlations (0 < r < 0.25) followed by southern portions with correlation coefficient between 0.25 and 0.5. The generally poor correlation between same month NDVI and precipitation indicates that the response of vegetation in most of the basin is not quick. Since the basin is located in moist, high altitudes, it is apparent that vegetation in these areas tends to utilize precipitation in a delayed manner compared to vegetation in water-scarce arid and semi-arid regions.

For the lag1 scenario, the correlation coefficient values follow the similar south–north decreasing pattern of lag0 condition. However, the magnitude of the correlation coefficients is one level higher than the lag0 situation. Moderately high correlations (0.5 < r < 0.75) are observed in the downstream portions of the basin. Agricultural and forested areas show the highest lag1 correlation while the alpine meadows show the lowest correlation. Therefore, the basin's cultivated and forested lands tend to respond at least a month later after the precipitation fall. The northern meadows still exhibit a low correlation with precipitation, indicating the pronounced influence of other growth factors such as temperature and topography than precipitation alone.

For the rest of the lags (lag2 and lag3), the correlation further diminishes, implying the absence of additional delayed vegetation response to precipitation changes.

An example scatter plot is given in Figure 6 to illustrate the relationship between monthly NDVI and the corresponding precipitation values for lag1 condition. Slope agriculture has shown a better correlation (with r = 0.673). Irrigated agriculture and coniferous forests exhibited a similar correlation (with r ≈ 0.6). Alpine meadows, on the other hand, showed the least correlation (with r = 0.448).

The general absence of a higher correlation coefficient (r > 0.75) over the basin might suggest that precipitation is not the only determining factor for vegetation activities. Other factors such as temperature, solar radiation, topography and soil properties might also play an essential role. The pronounced influence of these additional factors on vegetation is also reported in Granger & Schulze (1977) and De Long et al. (2015).

In this study, spatiotemporal vegetation response patterns to precipitation changes were examined in the Kunhar River basin using satellite-derived NDVI and gauge interpolated precipitation datasets between 1982 and 2015. Meanwhile, The long-term precipitation and vegetation trends were also explored using the nonparametric Mann–Kendall trend test and Theil–Sen slope estimator. Finally, the correlation between precipitation and NDVI was investigated using the nonparametric Spearman's rank correlation.

The mean growing season precipitation and NDVI exhibit a south–north decreasing pattern in the basin. Forested and agricultural areas of the basin received higher precipitation amounts compared to northern alpine meadows. Precipitation has shown a total increasing trend during 1982–2015. NDVI changes, on the other hand, offer an increasing trend in the basin's south and central sections and a decreasing trend in the northern ones.

A relatively better correlation (0.5 < r < 0.75) is obtained for lag1 condition in the lagged correlation analysis. Consequently, vegetation in the basin's southern sections responds at least 1 month after the precipitation event. On the other lags (lag2 and lag3), a weaker precipitation-NDVI correlation is observed, suggesting either vegetation in the Kunhar River basin responds quickly with a maximum of 1-month precipitation delay or environmental factors other than precipitation play a significant role in vegetation activities.

This study is part of the first author's doctoral research work, conducted at Zhengzhou University, China. The author wishes to acknowledge and is thankful to Pakistan Meteorological Department and the Water Development Authority of Pakistan for providing important and valuable data for the research.

The first author conducted this research during his doctoral degree program under the supervision of second author. The third author supervised the writing of this article and also reviewed it.

The authors declare no conflict of interest.

This research received no external funding.

Data cannot be made publicly available; readers should contact the corresponding author for details.