Trends of Rainfall Variability and Drought Monitoring Using Standardized Precipitation Index in a Scarcely Gauged Basin of Northern Pakistan

Abstract

This study focused on the trends of rainfall variability and drought monitoring in the northern region of Pakistan (Gilgit-Baltistan). Climate Hazards Group Infrared Precipitation with Stations (CHIRPS) model data were used for the period of 1981 to 2020. The Standardized Precipitation Index (SPI) was applied to assess the dry and wet conditions during the study period. The Mann–Kendall (MK) and Spearman’s rho (SR) trend tests were applied to calculate the trend of drought. A coupled model intercomparison project–global circulation model (CMIP5–GCMs) was used to project the future precipitation in Gilgit-Baltistan (GB) for the 21st century using a multimodel ensemble (MME) technique for representative concentration pathway (RCP) 4.5 and RCP 8.5. From the results, the extreme drought situations were observed in the 12-month SPI series in 1982 in the Diamir, Ghizer, and Gilgit districts, while severe drought in 1982–1983 was observed in Astore, Ghizer, Gilgit, Hunza-Nagar, and Skardu. Similarly, in 2000–2001 severe drought prevailed in Diamir, Ghanche, and Skardu. The results of MK and SR indicate a significant increasing trend of rainfall in the study area, which is showing the conversion of snowfall to rainfall due to climate warming. The future precipitation projections depicted an increase of 4% for the 21st century as compared with the baseline period in the GB region. The results of the midcentury projections depicted an increase in precipitation of about 13%, while future projections for the latter half of the century recorded a decrease in precipitation (about 9%) for both RCPs, which can cause flooding in midcentury and drought in the latter half of the century. The study area is the host of the major glaciers in Pakistan from where the Indus River receives its major tributaries. The area and volume of these glaciers are decreasing due to warming impacts of climate change. Therefore, this study is useful for proper water resource management to cope up with water scarcity in the future.

1. Introduction

Globally, the mean annual temperature is rising [1,2,3] and has resulted in variations in rainfall distribution, trend, and magnitude [3,4,5]. The variations in these climatic elements have put substantial impacts on a hydrological cycle [6,7]. Climate scientists are using varying spans of time series data to quantify changes in rainfall pattern and trend, which has resulted in flood and drought in different regions of the world [8]. Flood brings extensive damages and disruption in a short span of time in the areas of a river basin threatened by floods, well known as flood areas, and gets more attention worldwide as compared with drought. Drought causes significant socioeconomic and environmental impacts in a large area or in a region [9]. Droughts occur as a consequence of rainfall deficiency compared with the average for a season or years [10,11]. Less rainfall and high temperature lead to shortage of water and consequently result in drought, which causes severe impacts on the community and poses serious threats to the environment as well as economy [12,13,14]. It is hard to quantify the drought commencement and cessation because of its diverse nature, the topography of the area, and the spatial variability in rainfall as well as due to its slow-onset nature [15], and its aftermath can even be seen after the end of the drought event [16]. A decline in precipitation results in meteorological drought, and that ultimately leads to agricultural and hydrological drought [17,18,19]. Precipitation is the major factor that causes drought; however, an increase in temperature and a decrease in humidity also exacerbate the drought severity [20,21]. Moreover, the severity of drought can be more intense with less amount of soil moisture content and a longer duration of a dry period [17]. The direct impacts of drought are mortalities and malnutrition as it causes a shortage of water and food in a region [16]. More than 11 million people have lost their lives, and 2 billion people have been affected by drought in the 19th century, which is a substantial loss as compared with other natural hazards [22]. Drought also affects the agriculture sector and results in famine, socioeconomic loss, land degradation, and mass migration [15,20,23].

Droughts are divided into four types: meteorological, agricultural, hydrological, and socioeconomic [24,25]. Drought can be determined by its drastic effects; however, various methods are used to calculate the drought severity, occurrence, and frequency [26]. Information obtained about drought using various indices is valuable for drought monitoring, preparedness, and mitigation. Therefore, many methods/indices are developed by meteorologists and drought scientists (i.e., Standardized Precipitation Index (SPI) [27], Palmer drought severity index (PDSI) [28], Standardized Precipitation Evapotranspiration Index (SPEI) [29], Vegetation Cover Index (VCI) [30], Temperature Condition Index [31], and Precipitation Condition Index [32]. However, SPI is the simple and efficient method to calculate the drought and its intensity using precipitation data. Moreover, the significance of the SPI method is to estimate the possibility of drought occurrence regardless of the location, season, and climate of the area. For this reason, in this research, SPI is used to assess the drought occurrence, severity, and duration.

In 1900−2011, the precipitation in the eastern parts of North and South America, Northern Europe, and Northern and Central Asia increased, whereas in the Mediterranean, Southern Africa, and parts of Southern Asia, it decreased [33]. Various studies have analyzed rainfall trend with different techniques and historical climatic data. Khan et al., (2021) [14] assessed the precipitation and drought in Southwest Xinjiang, China, using SPI and MK tests. They observed that the study area faced a dry period from 1997–2015. The years 2001, 2011, and 2014 were reported as extremely dry periods. The MK test revealed the wet condition for a few met stations in north of the Tianshan Mountains. Rahman et al. (2021) [11] also assessed the rainfall pattern and drought condition in Khyber Pakhtunkhwa (KP) using the SR test and SPEI. The authors observed that 8 out of 15 met stations fall in a humid region, while 7 met stations fall in subhumid, semiarid, and arid regions. They observed dry conditions in all met stations during the study period, particularly in the northern and northeastern parts of the region, recording an increasing trend of drought. Rao [34] outlined substantial warming trends in the minimum, maximum, and mean temperatures in the Mahanadi River Basin (India). Salma et al. [35] estimated rainfall trends in the different climatic zones of Pakistan from 1976 to 2005 using the analysis of variance (ANOVA) method and Dennett T3 tests. They observed a decreasing trend of rainfall of −1.18 mm/decade throughout the country and revealed the adverse impact of decrease in rainfall in the form of drought. Raza et al. [36] used temperature data for the period of 1955–2010 for GB and concluded that temperature has increased in the recent decade (2001–2010) and become 0.33 °C warmer than in the period of 1955–2000. Khattak and Ali [37] carried out a detailed study to find trends in maximum monthly temperature and rainfall over Punjab Province using the MK test. They found increasing trends in both maximum monthly temperature (0.002 °C/year) and total annual rainfall (3.23 mm/year) over the region. Khan et al. (2020) [38] studied drought hazard assessment for the period of 2000–2015 using remote sensing and drought indices over a part of Punjab. The authors noted a meteorological drought spell in 2000–2002, 2004, 2009–2010, and 2012. They noticed multiple factors (temperature, precipitation, and soil moisture) involved in drought severity.

GCMs are the mathematical models widely used to project the future climate at the global scale [39]. Multiple GCMs have been developed to simulate coarse resolution CMIP5, which is a set of GCMs from the Intergovernmental Panel on Climate Change’s Fifth Assessment Report (IPCC AR5) [40]. CMIP5 has more than 50 GCM projections available for the global climate and has individual uncertainties for different regions. To overcome the issue of GCM uncertainties, multiple GCMs are coupled using MMEs to enhance the reliability of climate projections using different GCMs [41].

From the abovementioned literature, it is apparent that rainfall and temperature are the main driving factors of drought; thus, there is a dire need to assess the temperature and precipitation patterns of GB [36]. Many studies have been conducted on drought throughout Pakistan [11,14,20,38,42,43,44,45]; however, there is lack of drought and future prediction of rainfall studies in the mountainous region, especially in Gilgit-Baltistan. The GB region is part of the Upper Indus Basin (UIB), which is the source of water for the downstream region due to a nest of snow cover and glaciers. Therefore, this study has significance particularly in terms of proper water resource management. The shortage of water for the agriculture sector can lead to agriculture drought. Pakistan is an agroeconomic country and largely depends on snow and glacier melt water in the summer season. The aim of this study is to assess the temporal distribution of rainfall and meteorological drought, its severity, intensity, and frequency in Gilgit-Baltistan (GB) using Standardized Precipitation Index (SPI), and its trend was assessed using Spearman’s Rho and the Mann–Kendall trend method for the period of 1981–2020, and CMIP5 data were used for rainfall projection.

2. Materials and Methods

2.1. Study Area

Gilgit-Baltistan (GB) was formerly known as the Northern Areas of Pakistan. Geographically, GB lies between the latitudes of 34° N and 37° N and longitudes of 72° E to 77° E. GB is divided into seven administrative districts (Hunza-Nagar, Skardu, Astore, Gilgit, Diamir, Ghanche, and Ghizer). GB shares borders with China, Afghanistan, and India (Figure 1). GB covers the mountainous region of extreme North Pakistan, which includes the high-altitude areas of the Himalayas, Karakorum, and Hindu Kush ranges covering an area of 72,971 km² [36]. The GB climate is hot in summer and cold in winter. GB has huge masses of snow, glaciers, and more than 50 peaks over 6500 m high, including 5 of the world’s 17 highest peaks. K2 is the highest peak (8611 m) in Pakistan and the second highest in the world, located in Gilgit-Baltistan. The region has more glacial ice than anywhere on earth outside the polar region, and thus, it is considered a third pole [46]. GB hosts more than 5000 glaciers, which include three longest glaciers (Biafo, Baltoro, and Batura) outside the polar regions [47]. The water flow in the Indus River is characterized by snow and glacier melt within GB, which is substantial in summers. Snowmelts in GB contribute 50% of the total water flow in the Indus River [48]. Indus is the largest river of Pakistan, and the country’s agriculture and hydropower production depend on it. GB also contains the nation’s most important natural forests, extensive mineral reserves, a wealth of biodiversity, and a rich cultural and archaeological heritage. In the mountains of GB, the occurrence of natural hazards, such as flash floods, avalanches, landslides, debris flows, and glacial lake outburst floods (GLOFs), are common [48]. Like other mountainous regions of the world, the climate of GB Province is also changing.

2.2. Methodology

The CHIRPS precipitation data were acquired from the CHIRPS website [49] for the period of 1981–2020 on a monthly basis. The gridded netcdf file (commonly known as network common data form, it has self-describing software libraries and machine independent data formats to create, access, and share the scientific data) was opened in a QGIS environment [50] and was extracted for the region of interest with an extract by a mask tool. The data were saved in a tiff format to obtain the monthly precipitation data. Eventually, the raster data were converted from raster to point data and further used for analysis. We calculated the descriptive statistics (standard deviation (SD), skewness, and kurtosis) in XLSTAT (XLSTAT is an add-on for Excel use for data analysis in Excel) Addinsoft in Excel [51]. In the next step, we evaluated the drought on a 1- and 12-month time scale using the SPI method. The drought trend was calculated using the MK and SR trend test methods. Drought characteristics are also calculated using run theory. Eventually, the rainfall future projection data were used to observe the pattern of rainfall in mid- and late century.

2.3. Source of Data

The gridded precipitation data of CHIRPS were used for the baseline period (1981–2020) with a spatial resolution of 0.05° × 0.05° [49]. The CHIRPS data were used due for the scarcity of meteorological stations in the study area. CHIRPS is a satellite-based precipitation product in combination with a ground-based gauge network that provides good correlation with the observed precipitation data [52]. The district wise annual and monthly precipitation data were extracted for GB. The monthly precipitation data for six GCMs were extracted from the IPCC data distribution center for the period of 2021–2100 (

Table 1. Selected GCMs from CMIP5.

Modelling Center	GCM	Resolution
National Center for Atmospheric Research, USA [54]	CCSM4	0.9° × 1.25°
Meteorological Office, UK [55]	HadCM3	2.5° × 3.75°
Geophysical Fluid Dynamics Laboratory, USA [56]	GFDL-CM3	2° × 2.5°
National Institute for Environmental Studies, Japan [57]	MRI-CGCM3	1.12° × 1.12°
Canadian Centre for Climate Modelling and Analysis, Canada [58]	CanESM2	2.79° × 2.8°
Beijing Climate Center, China [59]	BCC-CSM1–1	2.81° × 2.81°

). The GCMs’ precipitation data were rescaled using a bilinear interpolation method to have similar resolution for the datasets [53].

2.4. Multimodel Ensembles (MMEs) for GCMs

MMEs are used for incorporating six GCMs’ precipitation data for the 21st century to reduce the uncertainties in the projections [40]. There are two approaches for MME (i.e., simple ensemble mean (SEM) and weighted ensemble mean (WEM) [60]. The SEM approach was applied, which is relatively simple, and its performance is better than the individual GCM on the rescaled GCMs [61]. SEM is calculated using simple arithmetic ensemble means for the precipitation of six GCMs for the study area (Equation (1)):

$$\mathit{S}\mathit{E}\mathit{M}\mathit{=}\frac{\mathit{1}}{\mathit{n}}{{\displaystyle \mathit{\sum }}}_{\mathit{p}\mathit{=}\mathit{1}}^{\mathit{n}}\mathit{G}\mathit{C}{\mathit{M}}_{\mathit{p}}$$

(1)

where n is the number of GCMs used for the SEM, and $\mathit{G}\mathit{C}{\mathit{M}}_{\mathit{p}}$ refers to the precipitation simulations of GCMs for the future. The future precipitation projections are reclassified into midcentury (2021–2060) and late century (2061–2100) for each district of the study area.

2.5. Standardized Precipitation Index (SPI)

SPI solely depends on precipitation data for evaluating drought irrespective of any other hydrometeorological factors. The CHIRPS precipitation dataset covering the 1981–2020 temporal period was used for the calculation of SPI. The classification system of SPI defines the intensity of drought (

Table 2. SPI index classification system [27].

SPI Value	Category	Abbreviation
≥2	Extreme wet condition	EXW
1.5 to 1.99	Severe wet condition	VEW
1.0 to 1.49	Moderate wet condition	MOW
−0.99 to 0.99	Near normal condition	N
−1.0 to −1.49	Moderate drought condition	MOD
−1.5 to −1.99	Severe drought condition	SED
≤−2.0	Extreme drought condition	EXD

). McKee et al. [27] created SPI in 1993. SPI is the number of standard deviations that can observe the deviation of rainfall from its mean and distribute it normally [27]. The SPI values range from −2 to 2. A value below zero indicates a drought condition, and a value above zero specifies a wet condition. SPI was designed to estimate the deficiency of rainfall at multiple time scales. It is the most efficient and reliable method for diverse topographical areas [62]. The mathematical formulation of SPI [27] is shown below.

$$\mathit{S}\mathit{P}\mathit{I}=\frac{\mathit{X}-\overline{\mathit{X}}}{\mathit{\sigma }}$$

(2)

Here, X is the monthly rainfall, $\overline{\mathit{X}}$ is the mean of the monthly rainfall, and σ represents the standard deviation. SPI was estimated using a drought indices calculator (DrinC) open-source software (it is developed by a center for the assessment of natural hazards and proactive planning and laboratory of reclamation works and water resource management, National University of Athens. It has a user-friendly interface to calculate meteorological, hydrological, and agricultural drought [63]. After calculation of SPI, the dry and wet conditions were categorized. The calculation can be made to evaluate drought events for short or long periods [64]. Short-term SPI is significant for agricultural purposes, whereas long-term SPI is for the study of water resources, such as groundwater supplies, surface water flow, and reservoir level. In this paper, the SPI results for 1 and 12 months were calculated and discussed.

2.6. Characterization of Drought

Drought was characterized using run theory to identify the drought events. Run theory is used on the basis of threshold values to identify the runs [65]. When the values are below or above the defined threshold level, the drought event will be defined accordingly. We defined a threshold value of −1 in XLSTAT Addinsoft. When the value goes below the defined threshold, it will be considered a drought, and a drought will be finished when the values surpass the threshold. Drought was characterized into drought duration, drought intensity, drought severity, and drought peak. The drought period starts when the SPI value is below −1, and the drought period is finished when the SPI value is above −1, which is defined as a drought duration and measure in months. Drought severity (DS) is a sum of negative SPI values for a drought duration. The average drought values in a drought period are drought intensity (DI), and the peak of the drought (DP) is the lowest value reached in a drought duration.

$$\mathit{D}\mathit{S}={{\displaystyle \sum }}_{\mathit{K}=\mathit{1}}^{\mathit{D}\mathit{D}} \mathit{i}\mathit{n}\mathit{d}\mathit{e}{\mathit{x}}_{\mathit{K}}$$

(3)

$$\mathit{D}\mathit{I}=\frac{\mathit{D}\mathit{S}}{\mathit{D}\mathit{D}}$$

(4)

Drought severity can be calculated using Equation (3) and the intensity of a drought can be derived using Equation (4), where, K represents month, and SPI values represent index_K in K month.

2.7. The Mann–Kendall Trend Test

The World Meteorological Organization (WMO) endorses the MK trend test, which is a nonparametric trend test used to assess the presence of significant trends in variables using time-series data. According to the MK test, the null hypothesis infers that the variable is independent and identically distributed (no trend), with the alternative hypothesis being that a monotonic trend occurs over time. Thus, the statistical X variables tested are calculated as follows:

$$\mathit{X}={{\displaystyle \sum }}_{\mathit{j}=\mathit{k}}^{\mathit{m}-\mathit{1}} {{\displaystyle \sum }}_{\mathit{k}=\mathit{j}+\mathit{1}}^{\mathit{m}} \mathit{s}\mathit{g}\mathit{n}\left(\mathit{y}\mathit{k}-\mathit{y}\mathit{j}\right),$$

(5)

where

$$\mathit{s}\mathit{g}\mathit{n}\left(\mathit{x}\mathit{k}-\mathit{x}\mathit{j}\right),=\left\{\begin{array}{cc}+\mathit{1},& \mathit{y}\mathit{k}>\mathit{y}\mathit{j}\\ \mathit{0},& \mathit{x}\mathit{k}=\mathit{y}\mathit{j}\\ -\mathit{1},& \mathit{x}\mathit{k}<\mathit{x}\mathit{j}\end{array}\right\}$$

(6)

X is normally distributed. If the mean is 0, then the variance is:

$$\mathit{V}\mathit{a}\mathit{r}\left(\mathit{S}\right)=\frac{\mathit{1}}{\mathit{18}}[\mathit{n}\left(\mathit{n}-\mathit{1}\right)\left(\mathit{2}\mathit{n}+\mathit{5}\right)-{{\displaystyle \sum }}_{\mathit{p}=\mathit{1}}^{\mathit{q}}{\mathit{t}}_{\mathit{p}}\left({\mathit{t}}_{\mathit{p}}-\mathit{1}\right)\left(\mathit{2}{\mathit{t}}_{\mathit{p}}+\mathit{5}\right)]$$

(7)

Here, ‘t’ is the range of any node. When n > 10, Z_MK joins a normal distribution and is computed as follows:

$${\mathit{Z}}_{\mathit{M}\mathit{K}}=\{\begin{array}{cc}\frac{\mathit{X}-\mathit{1}}{\sqrt{\mathit{V}\mathit{a}\mathit{r}\left(\mathit{X}\right)}}& \mathit{i}\mathit{f} \mathit{X}>\mathit{0}\\ \mathit{0}& \mathit{i}\mathit{f} \mathit{X}=\mathit{0}\\ \frac{\mathit{X}+\mathit{1}}{\sqrt{\mathit{V}\mathit{a}\mathit{r}\left(\mathit{X}\right)}}& \mathit{i}\mathit{f} \mathit{X}<\mathit{0}\end{array}$$

(8)

The positive values of the X statistic specify a decreasing drought trend, while the negative ones specify an increasing drought trend. The null hypothesis of no trend is excluded if |Z_MK| > [Z_1-a/2] [20]. A comprehensive depiction of the test can be found elsewhere [66,67]. The MK test has two benefits: (1) it does not need the data to be normally distributed and (2) sensitivity is very low to break the data due to inhomogeneity.

2.8. Spearman’s Rho Test

Spearman’s Rho is a nonparametric method used for trend analysis because it is very easy and simple method. The ‘B’ statistic and standardized test statistics Y_B are shown below in Equation (9). This method was developed by Spearman in 1904 [68].

$$\mathit{B}=\mathit{1}-\frac{\mathit{6}{{\displaystyle \sum }}_{\mathit{i}=\mathit{0}}^{\mathit{n}}{[\mathit{J}\left({\mathit{L}}_{\mathit{i}}\right)-\mathit{M}]}^{\mathit{2}}}{\mathit{n}\left({\mathit{n}}^{\mathit{2}}-\mathit{1}\right)}$$

(9)

$${\mathit{Y}}_{\mathit{B}}=\mathit{B}\sqrt{\frac{\mathit{n}-\mathit{2}}{\mathit{1}-{\mathit{B}}^{\mathit{2}}}}$$

(10)

Here, J(L_i) is used to rank the observations, and n indicates the length of the observations. The decreasing drought trend can be observed with positive values of Y_B and vice versa. The 95% of significance level was used in this study. When a positive or negative trend detects it, the null hypothesis is rejected.

3. Results and Discussion

3.1. Spatiotemporal Variability of Rainfall

The evaluation of rainfall variability in the seven districts of Gilgit-Baltistan was based on precipitation dataset of CHIRPS. According to the classification of rainfall regions, the Hunza-Nagar and Ghanche districts fall under cold arid region, while the rest of the districts fall under cold semiarid regions [69] (

Table 3. Statistics of annual rainfall at the district level.

Station	Min Rainfall (mm)	Max Rainfall (mm)	Mean Rainfall (mm)	SD	Skewness	Kurtosis
Astore	329.47	674.92	497.36	86.66	0.28	−0.47
Diamir	376.78	730.41	551.54	85.26	0.08	−0.52
Ghanche	157.25	326.67	227.19	42.40	0.33	−0.55
Ghizer	182.63	395.46	279.43	51.42	0.27	−0.55
Gilgit	190.43	391.20	299.49	46.50	−0.19	−0.32
Hunza-Nagar	172.48	343.13	244.37	44.02	0.48	−0.49
Skardu	214	477.35	332.55	61.67	0.34	−0.26

). The data were statistically evaluated to observe the mean, standard deviation, skewness, and kurtosis in the total annual precipitation for each district. The mean annual precipitation varies considerably over the study area, from 227.19 mm (Ghanche) to 551.54 mm in the Diamir district (Table 3). The standard deviation represents the variability in annual rainfall. The highest rainfall variability was observed in Astore (86.66 mm), followed by the Diamir (85.26 mm) and Skardu (61.67 mm) districts, while lowest rainfall variability was observed in the Ghanche district (42.40 mm), followed by Hunza-Nagar (44.02 mm), which represent the high rainfall variability in the region (Table 3). The data are considered normally distributed if the skewness and kurtosis have the value of zero or near zero, respectively. The results revealed that the data are positively skewed for all the districts except Gilgit (Table 3), which indicates an increasing trend of precipitation from what is normal. The results revealed that all districts have negative kurtosis, which implies a flat peak near the mean.

3.2. Temporal Analysis of SPI

The long-term SPI period is an integrated effect of short-term SPI; therefore, both 1-month and 12-month SPIs were calculated. The results of 12-month SPI is considered for the month of December, as presented in Figure 2. The results of 1-month and 12-month SPIs indicated drought and wet conditions for all districts during the study period. The 1-month SPI results revealed moderate to extreme drought in the period of 1981–2020 at different locations in Gilgit-Baltistan. In 1-month SPI, the Astore district witnessed 14 months of moderate drought and 6 months of severe drought from 1981 to 1990. It was observed that after the year of 1990 (1991–2020), the Astore district experienced 5 months of extreme drought, 25 months of moderate drought, and 3 months of severe drought (

Table 4. Frequency of drought for 1-month SPI.

District	Extreme Dryness (<−2.0)	Severe Dryness (−1.5 to −1.99)	Moderate Dryness (−1 to −1.49)
Astore	5	9	39
Diamir	3	17	46
Ghanche	3	7	37
Ghizer	4	5	47
Gilgit	5	14	45
Hunza-Nagar	2	20	45
Skardu	5	9	50

). The Diamir district went through 2 months of extreme drought, 11 months of severe drought, and 33 months of moderate drought in the period of 1981–2000; however, in the 21st century the number of drought months decreased, and the district witnessed 13 months of moderate drought, 6 months of severe drought, and 1 month of extreme drought (Table 4). The Ghanche district experienced 37 months of moderate drought, and 7 and 3 months of severe and extreme droughts, respectively, in the period of 1981–2020. The Ghizer district suffered 47 months of moderate drought, 5 months of severe drought, and 4 months of extreme drought conditions during the period of 1981–2020.

It was examined that the Gilgit district faced 30 months of moderate drought, 11 months of severe drought, and 2 months of extreme drought in 1981–1995. The dry condition prevailed in Gilgit in 1996–2020, and the district experienced 15 months of moderate drought, 3 months of severe drought, and 2 months of extreme drought (Table 4). Moreover, the Skardu district faced 64% of the moderate drought events in the first half of the study period (1981–2000), while 36% of the moderate drought events occurred in the period of 2001–2020. Similarly, the district went through 77% of the severe drought events in the period of 1981–2000, and only 23% (2) of the drought events happened in 2001–2020. However, it was determined that extreme drought events were high (60%) in 2001–2020 as compared with the period of 1981–2000 (40%) (Table 4). As compared with the Skardu district, Hunza-Nagar witnessed 93% of the moderate drought events in 1981–2000, while only 7% of the moderate drought events occurred in the period of 2001–2020. Additionally, 85% of the severe and 50% of the extreme dry spells occurred in 1981–2000, while the rest happened in the period of 2001–2020 (Table 4). The highest number of extreme drought events in 1-month SPI was experienced by Astore, Gilgit, and Skardu. The lowest number of severe drought events was observed in the Ghizer district, followed by Skardu and Astore. On the other hand, the highest number of moderate drought events was examined in Skardu, followed by Ghizer. The highest number of total drought months was observed in the Hunza-Nagar district (67), followed by Diamir (65) and Gilgit (64), and the lowest drought months were observed in the Ghanche district (47), followed by Astore (53) and Ghizer (56) in 1-month SPI.

From Figure 2, it can be clearly observed that GB faced various dry spells in the study period. The 12-month SPI results revealed severe to extreme drought conditions from 1982 to 1983 in almost all districts of GB (Figure 2). Another extreme drought spell started in 1989 and was finished in 1991 in all districts of GB except Ghanche. The dry spell was followed by a wet condition in 1996, while in 1998, a moderate drought event occurred in the Astore, Ghanche, and Skardu districts. The Astore, Diamir, Ghanche, and Skardu districts suffered from severe to extreme drought conditions in the year of 2000 and 2001; similarly, Gilgit and Hunza-Nagar went through a moderate drought condition (2000–2001). The results revealed three dry spells: first from 1982 to 1983, second from 1989 to 1991, and third from 1997 to 2002 (

Table 5. Frequency of drought for 12-Month SPI.

Districts	Extreme Drought	Severe Drought	Moderate Drought
Astore	2000	1982	1985, 1991, 1998. 2001
Diamir	1982	2000	1983, 1985, 1990, 1991, 2001
Ghanche		2000	1982, 1995, 1997, 1998, 2001, 2002, 2004
Ghizer	1982	1983, 1985	1988, 1989, 1990, 2001
Gilgit	1982	1983, 1985, 1990, 1991	2000
Hunza-Nagar		1982, 1983	1985, 1987, 1988, 1989, 1990, 1991, 1995
Skardu	2000	1982, 2001	1997, 1998, 2004

). The highest frequency of severe drought was observed in Gilgit, followed by Ghizer, Hunza-Nagar, and Skardu. On the other hand, highest number of moderate drought events was examined in the Hunza-Nagar and Ghanche districts, followed by Diamir, Astore, and Ghizer. The overall highest drought events were found in the Hunza-Nagar district, followed by Ghanche, Diamir, and Ghizer (Figure 3).

3.3. Characteristics of Drought

The longest duration of drought was observed in the Astore district in the year of 2000 with a drought severity (intensity) of −7.01 (−1.40). The Astore district witnessed another 4 months of drought period in the year of 1982 and 1985 with drought severities of −5.07 and −5.73, respectively (Figure 4, Figure 5 and Figure 6). The Diamir district went through 5 months of drought in 1982 and 1985 with drought severities (intensity) of −7.71 (−1.54) and −7.73 (−2.1), respectively (Figure 4, Figure 5 and Figure 6). The peak drought values were −1.7 and −1.76 in 1982 and 1985, respectively (Figure 7). The results revealed that the Ghanche district witnessed the longest drought duration in 1982 and 1983 with severities (intensity) of −4.81 (−1.2) and −4.65 (−1.6), respectively. However, the drought peak was observed in the year of 1992 with a value of −2.54 (Figure 7), and its intensity was −1.27 (Figure 6). From Figure 4, it can be examined that the Gilgit district witnessed the longest duration of drought in the years of 1982 and 1985 for 5 months each with a severity (intensity) of −7.79 (−1.55) and −8.04 (−1.63), respectively (Figure 5, Figure 6 and Figure 7).

The drought peak in Gilgit (−2.99) was observed in 2004 with an intensity of −2.99 (Figure 6 and Figure 7). Similarly, the highest drought severity was −7.94 in the Ghizer district in 1985 (Figure 5) when the district went through 5 months of dry spell (Figure 4). Four and 6 months of dry spell were observed in the Hunza-Nagar district in 1983 and 1989 (Figure 4) with severities (intensity) of −8.03 (−1.33) and −6.37 (−1.22), respectively (Figure 5 and Figure 6). The peak droughts in the Hunza-Nagar district were −1.69 and −1.88 in 1983 and 1989, respectively (Figure 7). The highest drought severity in the Skardu district (−6.06) was observed in the year of 1982 (Figure 5) with an intensity of −1.2 (Figure 6) with a dry spell of 5 months (Figure 4). The results of this study are in line with those of previous studies [45,70,71].

3.4. Mann–Kendall and Spearman’s Rho Trend Test of SPI

The Mann–Kendall and Spearman’s Rho trend tests were applied to 1- and 12-month SPI results, and the trend was evaluated with a 95% significance level. The results of the MK and SR tests are presented in

Table 6. SPI trend for 1 month and 12 months in Gilgit-Baltistan using MK and SR. Note: ZMK, Mann–Kendall; ZSR, Spearman’s Rho.

	Astore		Diamir		Ghanche		Ghizer		Gilgit		Hunza-Nagar		Skardu
	ZMK	ZSR	ZMK	ZSR	ZMK	ZSR	ZMK	ZSR	ZMK	ZSR	ZMK	ZSR	ZMK	ZSR
January	1.03	1.51	1.03	1.01	1.68	1.82	1.4	1.57	1.69	1.77	2.31	2.68	1.26	1.39
February	1.36	1.46	1.01	1.04	1.06	1.2	2.03	2.09	2.2	2.22	2.71	2.85	1.73	1.77
March	−0.6	−0.6	−0.8	−0.86	0.1	0.002	0.32	0.36	−0.02	0.002	1.71	1.92	−0.24	−0.27
April	2.15	2.25	1.75	1.84	1.59	1.64	2.08	2.21	2.08	2.23	2.29	2.94	1.89	2.02
May	2.52	2.67	2.2	2.39	0.52	0.61	2.14	2.41	2.69	2.82	2.48	2.66	2.31	2.36
June	1.43	1.33	1.92	2.01	2.1	1.95	1.95	2	1.96	1.87	1.09	1.17	1.26	1.12
July	−0.8	−0.8	0.01	0.005	−1.21	−1.01	1.63	1.5	1.82	2.1	1.17	1.43	−0.68	−0.51
August	0.8	0.88	0.89	0.94	0.4	0.49	0.09	−0.01	2.38	2.25	1.8	1.89	1.06	1.1
September	2.24	2.38	3.55	4.35	3.2	3.63	3.76	4.46	3.58	4.48	3.48	3.98	2.99	3.43
October	0.96	0.87	1.45	1.32	0.94	0.77	2.66	2.75	1.42	1.27	1.96	1.89	1.1	0.88
November	1.68	1.83	1.82	1.77	2.9	3.41	2.54	2.76	2.41	2.53	3.08	3.39	2.05	2.19
December	−0.7	−0.6	−0.99	−1.05	0.4	0.3	−0.9	−0.94	−0.22	−0.22	0.57	0.61	0.08	0.06
Annual	2.83	3.44	3.34	3.78	1.99	2.13	4.29	5.42	4.76	6.5	4.46	5.87	2.27	2.64

Note: The bold values indicate the existence of a significant trend in SPI results with a 95% significance level.

. The MK and SR values for 1-month SPI revealed a significant positive trend in the months of September and November for all districts. Similarly, except Skardu and Ghanche, other districts had a significantly decreasing drought trend in the month of April. In January only, Ghanche, Gilgit, and Hunza-Nagar experienced a significant positive trend. The months of March, July, and December did not have any significant positive or negative trend. The results revealed that most of the districts experienced an increase in the amount of rainfall. The increasing amount of rainfall is due to a decreasing snowfall in the winter season, which is creating an alarming situation throughout Pakistan. Similarly, the 12-month SPI results indicate a significant positive trend for all districts. The highest significant trend was observed for the Gilgit, Hunza-Nagar, and Ghizer districts. The increasing rainfall in the area is not a positive sign as it is just because of the warming effect of the climate, which transforms the form of precipitation from snowfall to rainfall. Not only will the decreasing snowfall and increasing rainfall create issues of water shortage in Gilgit-Baltistan, but also it will affect the agriculture sector in the whole country.

3.5. Future Precipitation Projections for 21st Century

Future projections of rainfall for each district of GB are classified and plotted into the midcentury (2021–2060) and late-century (2061–2100) projections (Figure 8). Rainfall is a complex climatic process, and its projections are also complex on spatiotemporal basis especially for mountainous regions such as GB. The temperature projections mostly endorse the storyline of RCPs predicting a smaller temperature rise under RCP 4.5 as compared with severe changes under RCP 8.5 [39]. Overall, the GB region rainfall projections depicted an increase of 4% for the 21st century as compared with the baseline period of 1981–2020. To analyze the spatial and temporal variability of rainfall, the projections are divided into midcentury and late century for each district (Figure 8). The scenario of RCP 4.5 projected a slight increase of 1.9%, while the scenario of RCP 8.5 depicted a robust increase of 11.8% for the midcentury period, which may cause excessive flow of surface water in the region for this time (Figure 8).

The future projection model depicted a severe decrease in rainfall for the latter half of the century under RCP 4.5 (7.2%) and RCP 8.5 (2.39%), which can cause scarcity of water in the region in the latter half of the century (Figure 8). GB has a diverse and unique topography that has impacted the rainfall projections of each district under both RCPs (Figure 8). The future projection data of rainfall were divided into two timelines (i.e., midcentury (2020–2060) and late-century (2061–2099) projections). The correlation coefficient of precipitation projections is weak for both RCPs because of its complex phenomena. It was observed that in some decades, the precipitation has increasing trend, while in others, a decreasing trend was examined. The Astore, Ghizer, Gilgit, and Skardu districts have a severe decreasing trend under both RCPs, particularly in the latter half of the century. On the other hand, the Diamir, Ghanche, and Hunza-Nagar districts projected a slight increase in precipitation under both RCPs especially in the midcentury regime (Figure 8a–n). Overall, GB depicted an increasing trend for the midcentury projections and a decreasing trend for the latter half of the century, which can cause floods in midcentury and a drought in the latter half of the century. The GB region is a major and a significant source of water supply for the whole country; therefore, the variation in rainfall projections will have severe impacts on the hydrological dynamics of the country. To efficiently deal with the rainfall variation projection in the upcoming decades, it is necessary to make better water management policies for a sustainable future of the country.

3.6. Discussion

The CHIRPS precipitation data for the period of 1981–2020 were used for the scarcely gauged area (GB) to quantify the drought conditions based on 1- and 12-month time scales using SPI. CHIRPS data have good agreement with station data over Pakistan and for the northern mountainous region [72,73]. The results of this study revealed that all districts of GB experienced short- and long-term drought ranging from moderate to extreme intensity. Extreme drought events were found in the year of 1982 and 2000, while severe droughts were mostly found in 1982, 1983–1985, 1990–1991, and 2000 for the 12-month time scale; moreover, moderate drought events were also found and classified into the groups of 1982–1983, 1985–1989, 1990–1991, 1997–1998, 2000–2002, and 2004. Some of the drought events found in this study are in line with previous studies [43,71,74,75,76]. Jamro et al. (2019) [71] also revealed three drought periods (i.e., 1946–1953, 1984–1989, and 1998–2002) in GB. Severe drought conditions were found in Asia during the period of 1999–2002 [70], which was also witnessed in this study for GB. Interseasonal rainfall variability due to atmospheric circulation is one of the possible reasons of drought and floods in Pakistan [77,78,79,80]. Rainfall variability defines droughts in Pakistan particularly in the cropping season; however, changes in monsoon and westerlies influence rainfall due to global warming, which can intensify the drought severity in Pakistan [78]. A shifting pattern of rainfall in the monsoon season due to global warming in South Asian countries will have severe impacts on the ecosystem in the 21st century [79]. It is expected that rainfall-dependent regions will face famine due to prolonged drought episodes as a result of climate variability [43,81]. The drought episodes of 1997–1998 and 2000–2002 in the region were influenced by a strong ENSO event for 5 years in South Asia. The possible reason for the lengthened drought episode (1998–2002) could be the warm ocean water, although in 1999 a cold wet episode occurred due to the influence of La Nina. This complex mechanism is defined by scientists as regional and large-scale climatic variability [82,83]. Due to longer-period drought episodes in 1997–1998 and 2000–2002 in the region, the agriculture sector of the country had a −2.6% growth rate. Adnan and Ullah (2020) [42] revealed that the GB Province receives 78% of its annual precipitation in the monsoon season. Ahmed et al. (2017) [44] observed a significant increasing trend of rainfall in the northern region of Pakistan; however, Bhatti et al. (2020) [84] found that the mean rainfall in the northern mountainous region is decreasing with rates of 0.17 and 0.8 mm with 95% and 99% percentile, respectively. However, glacier ice melt acts as a drought-resilient source of water in the northern mountainous region and reduces drought vulnerability for a larger population [85]. The results of the present research are consistent with the findings of some of the previous studies conducted in Pakistan.

The southern and central parts of Pakistan are considered to be more vulnerable to drought as compared with other regions of the country [78,86,87]. However, Ahmed et al. (2017) [44] and Ullah et al. (2018) [88] observed significant drought conditions in the mountainous region of Pakistan due to rainfall variability. Meteorological and agricultural drought conditions were also observed in the period of 2000–2018 for the Chitral Kabul River Basin [45]. It was observed by various researchers that the northern mountainous region of Pakistan was affected by moderate drought conditions during period of 1951–2010 [73,88,89]. A recent study also found drought conditions in the northern part of Pakistan due to the impact of less-than-normal rainfall [38]. The Khyber Pakhtunkhwa region experienced the worst drought event in 2000–2001 due to rainfall deficiency as compared to its normal rainfall [42].

We incorporated GCMs’ precipitation projections for the GB region under RCP 4.5 and RCP 8.5 to predict the precipitation condition in the future, which may help to adopt mitigation strategies to cope with future water dynamics in the region. The results revealed that during midcentury (2021–2060), the precipitation will increase by 1.9% and 11.8% under RCP 4.5 and RCP 8.5, respectively. However, in the latter half of the century (2061–2100) the climatic models showed a decreasing trend of precipitation with 7.2% and 2.39% under RCP 4.5 and RCP 8.5 respectively. Previously, there were no studies that specifically dealt with precipitation projections of the GB region using RCPs for the 21st century. Ali et al. (2017) [90] incorporated a precipitation projection of GB while studying the precipitation projections of Pakistan for the 21st century using old emission-based scenarios. The findings for the GB region depicted an overall increasing trend in precipitation for the 21st century. However, this recent study depicted the precipitation projection trends for each district of the region for the 21st century using multiple GCMs to avoid biases for a particular GCM under recent concentration pathways. This study is useful to national and provincial disaster management authorities, policymakers, water resource management departments, and farmers. It is necessary to monitor and assess long-term drought-affected areas for better water resource management. This study gives information about the future precipitation and provides a basis to properly manage water resources. To study the climatic variability, it is necessary to use a longer period of data; however, lack of long-period data and adequately spatially distributed meteorological station data is a big hindrance in the mountainous region. We planned to study the extreme climate indices of the region with respect to cryosphere changes because the downstream region is highly dependent on snow cover and glaciers nested in the Upper Indus Basin.

4. Conclusions

This study investigated the precipitation variability and assessed drought patterns in the northern mountainous region of Pakistan using a CHIRPS gridded dataset. The results of SPI indicated extreme to severe droughts in the years of 1982–1983, 1985, 1990–1991, 2000–2001, while the moderately dry years were 1977, 1984, 2007, 2014–2015, and 2017–2018 in the GB region. Overall, the results of the MK and SR trend tests deduced the increasing rainfall trend in historical data (1981–2020). The MMEs’ precipitation projections depicted an overall increase of 4% for the 21st century under both RCPs of the GB region. The increasing rainfall trend is dominant for the midcentury projections as compared with the latter half of the century. The district wise precipitation projections depicted an increasing trend for the northern districts of the region, while the southern districts of GB depicted a declining trend in precipitation for the 21st century. The increase in the amount of rainfall in the study area, especially in January, February, and March, indicates a shift of precipitation from snowfall to rainfall, which will increase the threat of dryness in the future, especially in the summer seasons in the Indus River Basin.