The Impacts of Climate Change on Vegetation Productivity and the Trends of the Impacts in Pakistan over the Last Two Decades

1. Introduction

Global terrestrial vegetation has exhibited a pronounced greening tendency in a rapidly changing climate, characterized by an overall increase in net primary productivity (NPP) [1]. However, this tendency displays strong spatial heterogeneity, with some regions experiencing declines in vegetation productivity [2]. Consequently, identifying the regional disparities in the responses of vegetation to climate change requires further in-depth investigation. Pakistan stands out as one of the few regions globally that possesses an enormous elevational gradient within a relatively compact geographic area [3]. From the southern coastal zone to K2 (Mount Godwin-Austen), the world’s second-highest peak on the China–Pakistan border, the elevation increases sharply from sea level to 8,611 m. This extensive vertical gradient causes the development of distinct climatic zones [4,5,6], serving as an ideal “natural laboratory” for investigating the effects of climate change on vegetation. Elucidating the spatial heterogeneity of the relationship between climate change and vegetation productivity is pivotal for formulating targeted strategies for vegetation conservation and sustainable management [7,8].

Previous studies have explored the relationships between climate and vegetation in Pakistan, yet the findings remain inconsistent. At the national scale, Chen et al. [9] reported that the NPP in Pakistan was mainly negatively correlated with temperature and precipitation. In contrast, Ahmad et al. [10] reported that the normalized difference vegetation index (NDVI) was mainly positively correlated with both factors and was strongly associated with precipitation. Additionally, Zheng et al. [11] and Anees et al. [12] reported that the NDVI generally exhibited a negative correlation with temperature but a positive correlation with precipitation across the country. Given such the dramatic elevational gradient within a compact geographic area, identical climatic changes can trigger markedly different ecological responses across regions, suggesting that the impacts of climate change on vegetation cannot be generalized nationwide. To address this, several studies have examined climate‒vegetation relationships within specific geomorphological zones or have conducted comparative analyses across subregions [13,14,15,16]. For instance, Mehmood et al. [17] reported that in the northern mountainous region, the NDVI was negatively correlated with precipitation and positively correlated with temperature, whereas the opposite relationships, exhibiting positive correlations with precipitation and negative correlations with temperature, were observed elsewhere. However, correlation analyses reveal only statistical associations rather than causality. While they identify which climatic factors are related to vegetation dynamics, they fail to quantify the magnitude or actual direction of climatic forcing on the NPP. Crucially, a positive correlation does not guarantee a greater NPP if the driving climatic factor itself is declining. Consequently, correlation-based analyses alone are insufficient to clarify the actual contributions of climate change to vegetation dynamics, and whether its overall impact on vegetation is beneficial or detrimental remains to be determined.

Moreover, vegetation growth is an ever-changing dynamic process [18,19], the impacts of climate change on vegetation are not static but are dynamically evolving [20,21,22,23]. Existing global-scale studies examining the effects of single environmental factors have shown that the correlation between solar-induced chlorophyll fluorescence (SIF) and precipitation has generally strengthened, with substantial regional heterogeneity in Pakistan [24]; moreover, the NDVI in dryland ecosystems has been shown to increase with precipitation [25], and the positive correlation between gross primary productivity (GPP) and summer temperature in the Northern Hemisphere may undergo a future reversal [26]. However, current research on the changing relationships between climate and vegetation generally focuses on single climatic variables and typically covers broad spatial scales, overlooking the complex regional variations within Pakistan. Consequently, the dynamic responses of Pakistan’s vegetation to the combined impacts of temperature and precipitation have yet to be fully explored.

In this study, the spatiotemporal changes in climate and vegetation, their relationships, and the evolving dynamics of these relationships were investigated. Specifically, we aim to address the following core questions: Across Pakistan’s pronounced elevational gradient, in which regions does climate change enhance vegetation productivity, and where does it inhibit it? How are these impacts expected to evolve in the future? By answering these questions, this study provides theoretical support for predicting future vegetation dynamics and mitigating climate-related risks, offering a scientific foundation for the development of region-specific vegetation management strategies.

2. Materials and Methods

2.1. Study Area

Pakistan is located in the northwestern part of the South Asian subcontinent (23.5°N–37°N, 61°E–78°E) (Figure 1). In the southern coastal regions, the mean annual temperature reaches approximately 30°C, whereas the temperatures decrease progressively with increasing elevation and decrease to below 0°C in the northern mountainous areas. The central and southern parts of the country are mainly arid, with mean annual precipitation amounts of less than 500 mm [27]. In contrast, owing to the abrupt topographic relief, the Indian Ocean monsoon brings abundant rainfall to the northern region. The annual precipitation in these areas exceeds 1,000 mm and can locally surpass 2,000 mm.

2.2. Data Acquisition and Processing

2.2.1. Datasets

Monthly mean temperature and monthly total precipitation data from 2001 to 2022 were obtained from the ERA5-Land climate reanalysis dataset provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) [28] at a spatial resolution of 0.1° (approximately 9 km). The NPP data were derived from NASA’s MODIS MOD17A3HGF (Version 6.1) annual net primary productivity product for the period 2001–2022 [29,30], with a spatial resolution of 500 m. Land cover information was obtained from the MODIS MCD12Q1 (Version 6.1) annual land cover product [31,32], with a temporal coverage of 2001–2022 and a spatial resolution of 500 m.

2.2.2. Data Processing

To ensure spatial consistency, the NPP and land cover datasets were resampled to 0.1°. To minimize the confounding effects of anthropogenic activities on climate‒vegetation relationships, pixels that experienced land cover transitions between 2001 and 2022, as well as those classified as urban or built-up lands [33], were excluded from the analysis of the relationships between climate and NPP.

2.3. Statistical Analysis

2.3.1. Analysis of Temporal Trends and Dynamic Relationships

Pixel-based linear regression was applied to estimate the temporal trends in mean annual temperature, total annual precipitation, and NPP. A partial correlation analysis was conducted for each grid cell to quantify the partial correlation between the NPP and each climatic driver (mean annual temperature and total annual precipitation). To capture the temporal evolution of these climate‒vegetation relationships, a moving partial correlation analysis was performed using a 5-year moving window [34].

2.3.2. Classification of Climatic Impacts on Vegetation

Multiple linear regression was used to assess the contributions of temperature (temp) and precipitation (prec) to NPP variability. The decomposition of NPP changes is expressed as follows:

$${\text{∆}\mathrm{NPP}}_{total}={\text{∆}\mathrm{NPP}}_{temp}+{\text{∆}\mathrm{NPP}}_{prec}+{\text{∆}\mathrm{NPP}}_{res}$$

(1)

∆NPP_total represents the total linear trend of the actual NPP during the study period; ∆NPP_temp and ∆NPP_prec denote the individual contributions of temperature and precipitation variations, respectively, to NPP dynamics; and ∆NPP_res refers to the residual contribution. We identified the main climatic driver of NPP variability by comparing the absolute magnitudes of NPP changes that were attributable to temperature (temp) and precipitation (prec). On the basis of the temporal trends of the main climatic driver in conjunction with its moving partial correlation with the NPP (Table 1), we classified the specific modes of climatic impacts on vegetation productivity.

3. Results

3.1. Regional Differentiation of Temperature and Precipitation Variations from 2001 to 2022

In Pakistan, 42.57% of the land area exhibited a continuous increase in mean annual temperature, with an average warming rate of 0.120°C per decade. These areas were located mainly across the northern mountainous regions approximately north of 35°N and east of 72°E, the western parts of the Gandhara Plain and the Potwar Plateau, and regions south of 26°N. Areas with rapid warming were concentrated in the mountainous areas north of 35°N, the southwestern part of the Potwar Plateau, and the western parts of the Makran Cosat Range, where the maximum warming rate reached 0.505°C per decade. Conversely, 57.43% of the land area showed a continuous decrease in mean annual temperature, with an average cooling rate of 0.117°C per decade. These areas were located mainly in the regions between approximately 26°N and 32°N. Rapid cooling occurred east of the Sulaiman Mountains and west of the Indus River, where the maximum cooling rate was 0.450°C per decade (Figure 2).

In parallel, 76.20% of the land area exhibited a continuous increase in total annual precipitation, with an average rate of increase of 36.82 mm per decade. These areas were distributed mainly south of 32°N. Areas with rapid precipitation increases were concentrated northeast of Lahore, east of the Sulaiman Mountains, west of the Indus River, southeast of the Central Brahui Range, and in the Indus River Delta, where the maximum rate of increase reached 140.08 mm per decade. Conversely, 23.80% of the land area showed a continuous decrease in total annual precipitation, with an average rate of decrease of 25.74 mm per decade. These areas were located mainly in the Hindu Kush Mountains, the Gandhara Plain, the western part of the Potwar Plateau, and south of the Chagai Desert and Chagai Hills. Rapid decreases occurred in the Hindu Kush Mountains, the Gandhara Plain and the western Potwar Plateau, where the maximum rate of decrease was 129.35 mm per decade (Figure 3).

3.2. Regional Differentiation of NPP Variations from 2001 to 2022

From 2001 to 2022, the spatial distribution of the NPP variations in Pakistan exhibited considerable heterogeneity. Specifically, 62.61% of the land area exhibited a continuous increase in NPP, with an average increase of 20.28 gCm^-2 per decade. These areas were located mainly in the northern mountainous regions, the Gandhara Plain, the Potwar Plateau, most parts of the Indus Plain, the Sulaiman Mountains, the Central Brahui Range, and the Kirthar Mountains. The areas where the NPP increased rapidly were concentrated in the regions between 33°N and 35°N, where the maximum rate of increase reached 261.78 gCm^-2 per decade. In contrast, 2.39% of the land area showed a continuous decrease in NPP, with an average decrease of 4.37 gCm^-2 per decade. These areas were distributed sporadically across the northern mountainous regions, the section of the Indus Plain between 31°N and 32°N, south of the Kirthar Mountains, the Indus River Delta, and the southern part of the Thar Desert, where the maximum rate of decrease reached 55.67 gCm^-2 per decade. Areas with persistently near-zero NPP accounted for 35.14% of the land area and were located mainly in the Karakoram Mountains, west of the Kirthar Mountains, and the north part of the Thar Desert (Figure 4).

3.3. Partial Correlations and Their Temporal Trends Between Climatic Drivers and NPP from 2001 to 2022

A positive partial correlation between mean annual temperature and NPP was observed in 16.43% of the land area (Figure 5a). Specifically, 11.99% of the land area exhibited an increasing trend in this positive partial correlation (Figure 5b), which is distributed primarily in the mountainous regions north of 35°N, west of the Sulaiman Mountains, and the Indus Plain between 26°N and 31°N. Conversely, 4.44% of the land area showed a decreasing trend in this positive partial correlation (Figure 5c), which is distributed mainly in the mountainous regions between 34°N and 36°N. In contrast, a negative partial correlation between mean annual temperature and NPP was observed in 21.25% of the land area (Figure 5a). Specifically, 6.83% of the land area exhibited an increasing trend for this negative partial correlation (Figure 5d), which is distributed primarily on the Gandhara Plain, the northern Potwar Plateau, the northeastern Upper Indus Plain, and the western flank of the Kirthar Mountains. Conversely, 14.42% of the land area showed a decreasing trend for this negative partial correlation (Figure 5e), which is distributed mainly in the Upper Indus Plain between 30°N and 33°N, the Central Brahui Range, the Kirthar Mountains, the Indus Basin and the southern part of the Thar Desert.

A positive partial correlation between total annual precipitation and NPP was observed in 27.39% of the land area (Figure 6a). Specifically, 16.66% of the land area exhibited a strengthening trend in this positive partial correlation (Figure 6b), which is distributed primarily in the Sulaiman Mountains, the Central Brahui Range, the northern Kirthar Mountains, and the Indus Plain between 25°N and 31°N. Conversely, 10.73% of the land area showed a weakening trend in this positive partial correlation (Figure 6c), which is located mainly in the Gandhara Plain and the northern Potwar Plateau, the southern Kirthar Mountains, the Indus Plain in southern Sindh, and the southern part of the Thar Desert. In contrast, a negative partial correlation between total annual precipitation and NPP was detected in 10.29% of the land area (Figure 6a). Specifically, 4.71% of the land area exhibited a strengthening trend in this negative partial correlation (Figure 6d), which is distributed primarily in the northern mountainous regions and the northeastern Upper Indus Plain. Conversely, 5.58% of the land area showed a weakening trend in this negative partial correlation (Figure 6e), which is distributed mainly in the northern mountainous regions and the southern part of the Upper Indus Plain.

3.4. Regional Differentiation of the Responses of Vegetation Productivity to Temperature and Precipitation from 2001 to 2022

As shown in Table 1, temperature and precipitation affected the vegetation productivity in 16 distinct ways, six of which exhibited relatively clustered distributions, with each covering more than 1% of the land area. The regions in which the mean annual temperature was the main driver accounted for 12.22% of the land area. The areas in which warming enhanced the vegetation productivity, with this positive impact intensifying, accounted for 2.70% of the land area; these areas were distributed mainly in the northern mountainous regions (elevation > 3,500 m a.s.l.). The areas in which warming enhanced the vegetation productivity, with this positive impact diminishing, accounted for 1.92% of the land area; these areas were distributed mainly in the northern mountainous regions (elevation < 3,500 m a.s.l.). The areas in which cooling enhanced the vegetation productivity, with this positive impact diminishing, accounted for 4.38% of the land area; these areas were distributed mainly in the northern Upper Indus Plain (Figure 7).

Regions where the total annual precipitation was the main driver accounted for 25.46% of the land area. The areas in which decreasing total annual precipitation amounts enhanced the vegetation productivity, with this positive impact diminishing, accounted for 1.88% of the land area; these areas were distributed mainly in the Gandhara Plain and the northern Potwar Plateau. The areas in which increasing total annual precipitation amounts enhanced the vegetation productivity, with this positive impact intensifying, accounted for 13.92% of the land area; these areas were distributed mainly in the Sulaiman Mountains, the Central Brahui Range, the northern Kirthar Mountains, and the Indus Plain between 26°N and 31°N. The areas in which increasing total annual precipitation amounts enhanced the vegetation productivity, with this positive impact diminishing, accounted for 8.35% of the land area; these areas were distributed mainly in the southern Kirthar Mountains, the Indus River Delta, and southeastern Sindh Province (Figure 7).

4. Discussion

Owing to the pronounced elevational gradient in Pakistan, the ecological impacts of climate change on vegetation vary spatially, even when the magnitude of climate change is the same. In the northern mountains, the altitudinal belt at approximately 3,500 m a.s.l. acts as a transitional zone for vegetation responses to temperature. Above this zone, warming enhances the vegetation productivity, with this positive impact intensifying. If climate change continues at its current pace, it is unlikely to limit the vegetation productivity. These regions sustain several IUCN-listed threatened species, including the snow leopard (Panthera uncia) [35] and the Kashmir musk deer (Moschus cupreus) [36]. Enhanced NPP likely improves forage availability and habitat quality for these species, facilitating their survival and reproduction. Notably, the population of Pakistan’s national animal, the markhor (Capra falconeri), has steadily recovered. As a result, the markhor was downlisted to "Near Threatened" on the IUCN Red List in 2014 [37]. In contrast, below that zone, warming enhances the vegetation productivity, but this positive impact is diminishing, coinciding with an observed decline in total annual precipitation. Zheng et al. [38] and He et al. [39] reported that in regions with limited water availability, warming-induced drought can lead to a decline in tree productivity. Thus, precipitation may become the limiting driver for the vegetation productivity in this region, where vegetation may undergo risks of degradation. This region remains a vital forestry base [40], hosting abundant resources and plant species of significant medicinal and economic value, such as the Himalayan yew (Taxus wallichiana) and Himalayan Mayapple (Podophyllum hexandrum) [41], both of which are protected under Pakistan’s National Biodiversity Strategy and Action Plan. Therefore, the implementation of adaptive vegetation management strategies is urgently needed to address these climatic shifts.

Eastern Pakistan is an important agricultural region, particularly within the Indus River Basin. In regions such as Potwar, croplands are entirely rainfed and lack irrigation infrastructure, rendering crop production systems highly sensitive to climate change, as even slight climatic variations can lead to severe crop losses [42]. In this region, the positive impact of total annual precipitation on the vegetation productivity is weakening, while the region is experiencing a warming trend. Continued warming coupled with limited water availability may lead to drought, thereby imposing additional effects on vegetation productivity. In the northern Upper Indus Plain, cooling enhances the vegetation productivity; however, this promoting impact is weakening. Sustained cooling can rapidly alter or even reverse the relationship between temperature and vegetation productivity when temperatures deviate from the optimal growth window of crops [43]. Therefore, if climate change continues at its current pace, it may limit the vegetation productivity in this region. Conversely, on the Indus Plain between 26°N and 31°N, increasing total annual precipitation amounts enhance the vegetation productivity, and this beneficial impact is intensifying. If climate change continues at its current pace, it is unlikely to limit the vegetation productivity.

5. Conclusions

The NPP in Pakistan exhibited an overall increasing trend from 2001 to 2022. The spatial extent in which NPP is primarily affected by precipitation is approximately twice as large as the area primarily affected by temperature. The impacts of climate change on vegetation are mainly positive. If climate change continues at its current pace, the vegetation distributed in the northern mountainous regions below 3,500 m a.s.l., the Gandhara Plain, the northern Potwar Plateau, the northern Upper Indus Plain and regions south of approximately 26°N may undergo risks of degradation. Our study provides empirical evidence highlighting the importance of spatial differentiation when assessing the impacts of climate change on vegetation productivity. Furthermore, even if the current climate has a positive impact on vegetation growth, monitoring the temporal dynamics of the impact is essential to accurately predict whether vegetation might face the risk of degradation in the future.