Droughts over Pakistan: possible cycles, precursors and associated mechanisms

Abstract

In the recent few decades, climate variability had severely affected the socio-economic and environmental conditions worldwide. Frequent shifts in the atmospheric circulation patterns affect large parts of the globe, predominantly the arid and semi-arid regions facing severe to moderate droughts. Therefore, precursors of drought events and their associated mechanisms are important to understand. This study explores the possible cycles and precursor conditions that might be employed for predicting upcoming droughts in Pakistan. Standardize precipitation index and the single Z-index are used to detect and rank the drought years. Moreover, composite analysis is carried out to explore the large-scale circulation anomalies related to extreme drought events. Results demonstrate that extreme drought events are highly correlated with wind patterns and intrinsic weather system in the Pacific and Indian Oceans. This study analyzed air temperature, sea level pressure and geopotential height in the average time-period of January to March, sea surface temperature from October to December, and wind vectors in March to May as precursors that could be employed to predict the occurrence of droughts in Pakistan. This information is of significance for policymakers to plan climate change adaptive measures accordingly.

Keywords:

• Pakistan
• drought severity
• composite analysis
• drought precursors
• prediction

1. Introduction

Droughts are among the most intricate natural hazards (Portela et al. 2015) that generally initiate with a rainfall deficit and further result in enormous losses in water resources (Zhang et al. 2015), agricultural production (Al-Kaisi et al. 2013) and economy of the region (Khan et al. 2019; Xu et al. 2018). In the past few decades, many severe droughts have occurred all over the World. For instance, Texas drought in the USA (2012), drought in California (2012-2015), droughts in East Africa (2010-2011), drought in Balochistan Pakistan (1997–2003), droughts in China (1994, 1997, 2002), the Tibetan Plateau drought in (1997–2000), Indian subcontinent drought (2002), the worst drought in Papua (1997-1998) and the Australian Millennium drought (1997–2010) (McVicar and Bierwirth 2001; Prabhakar and Shaw 2008; Sheffield et al. 2009; van Dijk et al. 2013; Zhang et al. 2013; Ahmed et al. 2016; Hao et al. 2018). The regions most exposed to prolonged and severe droughts during 1951–1970 were the Central United States, the Argentinian Pampas, Russia, and Central Australia; during 1971–1990 they were Southern Chile, the Sahel, and Siberia; during 1991–2010 they were the Amazon Forest, the Congo River Basin, Mongolia, North Eastern China, and Borneo (Spinoni et al. 2014). These large-scale droughts have caused colossal damage to the agriculture sector, devastatingly declining water supply and crop production (Nielsen-Gammon 2011; Peterson et al. 2012; Dutra et al. 2013; Hoerling et al. 2014). These impacts are more intense in a country like Pakistan that is among the most susceptible countries to droughts and famines, considering its massive agricultural dependency.

Climatologically, most parts of Pakistan are arid to semi-arid, with significant spatial and temporal variability in climatic parameters. Monsoon rains are considered as a dominant hydro-meteorological resource for Pakistan as it contributes approximately 59% of annual rainfall. According to the climate risk index (CRI), Pakistan ranked in the 5^th position among the 10 most-affected countries to climate change from 1999 to 2018 (Eckstein et al. 2020). Over the past decade, various studies reported variations in precipitation pattern (Syed et al. 2010; Latif and Syed 2016; Latif et al. 2017) and warming trend of 0.1 °C per decade from 1960 to 2010 in Pakistan has been documented (Farooqi et al., 2005; Afzaal et al. 2009; Wang et al. 2011; Salma et al. 2012; del Rio et al. 2013; Shakoor et al. 2015; Khan et al. 2019). The year 1997-98 was marked by strong El Niño Southern Oscillation (ENSO), and this strong ENSO induced a five-year drought from 1998-2002 in the southwest Asia, including Pakistan, which experienced the most severe drought conditions in the last 50 years (Barlow et al. 2002; Rasul et al. 2008; Ahmad et al. 2010; Ahmad et al. 2014). Given the close relationship between atmospheric circulation patterns and climatic variables, extreme weather phenomena like drought and flood can be attributed to variations of atmospheric circulation patterns (Omidvar et al. 2016).

Future changes in the climate of Pakistan projected by most global climate models would cause a major impact on water resources, and subsequently affect food supply, health and ecosystem sustainability. Problems will most likely arise in southern parts of the country where the water resources are already under stress because that stress would be exacerbated by changes in supply or demand associated with climate change (Shakoor et al. 2015; Eckstein et al. 2020). The extremes in the long-term precipitation variations manifest themselves in the form of large-scale floods and droughts and cause devastating human and economic losses in conjunction with other physical, social, political and economic effects (Solomon et al. 2008).

Recently, many scientists investigated climate change-induced extreme events like droughts in Pakistan. For instance, Khan et al. (2020) provided an understanding and potential of drought prediction models in Pakistan by machine learning techniques (Khan et al. 2020). Jamro et al. (2020) explored the Standardized Precipitation Evapotranspiration Index (SPEI) based Spatio-temporal variability of drought characteristics in Balochistan province as well as in the entire region of Pakistan (Jamro et al. 2019; Jamro et al. 2020). Qaiser et al., (2021) presented a study about the evaluation of seasonal droughts using a composite drought index in Northern Punjab. Also, they described the related atmospheric dynamics (Qaiser et al. 2021). Hina and Saleem (2019) calculated the severity and magnitude of droughts in different provinces of Pakistan and highlighted the drought-prone areas as well (Hina and Saleem 2019). Ahmed et al. (2018) reported the impacts of climate variability on characteristics of seasonal drought in Pakistan and also studied the characterization of droughts in Kharif season (the first sowing season; starts from April to June and harvested from October to December) and Rabi seasons (the second sowing season; began in October to December and harvested from March to May) in the province of Balochistan (Ahmed et al. 2016; Ahmed et al. 2018). Another study conducted by Adnan et al., 2016, highlighted the shifting of precipitation and temperature trends along with drought vulnerability in agro-climatic zones in Pakistan (Adnan et al. 2016).

Although literature reports many studies in the region of Pakistan for drought characterization at national and regional scales (Xie et al. 2013; Ahmed et al. 2016; Haroon et al. 2016; Adnan et al., 2018; Ullah et al. 2021), but none of the studies mentioned above investigated the precursor conditions associated with the occurrence of droughts in the region. Therefore, it is essential to conduct new research analyzing the climatic parameters, some months prior to the summer season with the aim of using their status as a prediction tool for droughts. This study is designed to (Portela et al. 2015) detect the severity and recurrence of drought events in Pakistan (Zhang et al. 2015) identify and ascertain the statistical relationship among summer rainfall of extreme drought years and meteorological parameters like circulation patterns, geopotential height, air temperature, sea surface temperature (SST), and sea level pressure (SLP) prior to rainfall season and (Al-Kaisi et al. 2013) to figure out the drought precursors in Pakistan to anticipate the occurrence of droughts ahead of time.

2. Data and methods

This study is planned to infer Pakistan's rainfall behaviour using ground station precipitation data and deduce the relationship between local precipitation and regional climatic parameters to determine the drought precursors.

The monthly rainfall (mm) data covering the period of 37 years (1981-2017) was obtained from the Pakistan Meteorological Department (PMD). The names of weather stations together with their latitude, longitude and elevation are presented in Table 1. The rain gauge data was used for time series and return level analysis to figure out the general rainfall and drought patterns.

Table 1. List of weather stations in Pakistan.

Station no.	Station name	Elevation (MSL)	Latitude (^0N)	Longitude (^0E)
1	Nokkundi	682	28.81	62.76
2	Dalbadin	848	28.98	64.46
3	Panjgur	968	26.75	64.00
4	Lasbella	87	26.23	66.62
5	Khuzdar	1231	27.50	66.38
6	Kalat	2015	29.17	66.44
7	Sibbi	133	29.55	67.88
8	Quetta	1626	30.2	66.58
9	Zhob	1405	31.21	69.28
10	Barkhan	1097	29.92	69.74
11	Jacobabad	55	28.59	68.47
12	Badin	09	24.61	68.62
13	Hyderabad	28	25.23	68.25
14	Chhor	05	25.67	69.30
15	Padidan	46	26.91	67.98
16	Khanpur	88	28.67	70.66
17	Bahawalpur	110	29.35	71.93
18	Bahawalnagar	161	29.96	73.21
19	Multan	122	30.49	71.33
20	Lahore	214	31.52	74.40
21	Sargodha	187	32.14	72.76
22	Sialkot	255	32.52	74.65
23	Jhelum	287	32.96	73.76
24	Mianwali	210	32.71	71.55
25	D.I.Khan	171	31.93	70.85
26	Peshawar	327	34.11	71.61
27	Balakot	995	34.79	73.45
28	Chitral	1498	36.06	71.74
29	Sakardu	2317	35.39	75.66
30	Gilgit	1460	36.02	74.49

Moreover, to isolate and categorize the drought years based on severity, SPI and single Z-Index were used as themost appropriate method for subsequent analysis.

After the isolation and classification of drought years, the rainfall received during drought years was correlated with regional meteorological parameters to anticipate Pakistan's droughts a few months earlier.

Regional meteorological parameters including air temperature, SST, SLP, wind variables and geopotential height were sourced from the National Center for Environmental Prediction and National Center for Atmospheric Research NCEP-NCAR reanalysis data (with a horizontal resolution of 2.5° × 2.5°), during 1981-2017. A domain covering the region between latitudes 20-140° E and longitude 0-80° N was selected, considering that it is sufficient to cover the influence of circulation patterns responsible for occurrence of droughts in Pakistan (Khan et al. 2019). The detailed methods are presented in the following sections (Figure 1).

Figure 1. (a) Spatial rainfall distribution map of Pakistan (1981-2017), (b) Anomalous precipitation patterns in Pakistan (1981–2017).

2.1. General rainfall behaviour

Pakistan receives less rainfall in comparison with other regions of South Asia and hence predominantly sensitive to rainfall variability (Lu et al. 2020). The average precipitation of Pakistan ranges between 30 and 400 mm/year from south to north with a maximum gradient of 900-1800 mm in the sub-mountainous area (Adnan et al. 2016). The months from June to September are considered as monsoon season in Pakistan which contributes nearly 57% of the total annual rainfall (Hussain and Lee 2016). During summer monsoon season, rainfall is much greater in the northern parts of Pakistan than the southern parts (Hanif et al. 2013). For the purpose of showing the rainfall behavior and to deduce the dry and wet years, the time series of summer rainfall from 1981 to 2017 was plotted (Figure 2).

Figure 2. Pakistan's standardized total rainfall for June-August from 1981 to 2017. The solid blue line represents 11-year running mean.

2.2. Drought indices

In order to isolate and categorize the drought years in the region, two drought indices have been calculated in this study. A brief description of the indices is presented as follows.

2.2.1. Standardized precipitation index (SPI)

SPI index has been widely used to calculate the rainfall deficit and subsequent droughts during a certain time-period in the area of interest (McKee et al. 1993). It can be estimated for different timescales, including monthly, quarterly, biannually and annually (McKee 1995). The SPI well separates the wet and dry years (Soro et al. 2014). According to Mishra & Singh, SPI can be employed for short-term (soil moisture) as well as long-term (streamflow) water resources management (Mishra and Singh 2010).

For the SPI calculation, the long-term precipitation record for a desired period is fitted to a probability distribution, which is then transformed into a normal distribution so that the mean SPI for the location and desired period is zero (McKee et al. 1993; Edwards 1997). Positive SPI values indicate greater than median precipitation, and negative values indicate less than median precipitation. Since SPI is normalized, wetter and drier climates can be represented in the same way.

Thom (1958) found the gamma distribution to fit well to the climatological precipitation time series. The gamma distribution is defined by its frequency or probability density function: (1) $$g(x)=\frac{1}{{\beta }^a\Gamma (a)}{x}^{a-1}{e}^{-x/\beta },\hspace{1em}\text{for\hspace{0.5em}}x>0$$(1) in which α and β are the shape and scale parameters respectively, x is the precipitation amount and Γ(α) is the gamma function. Parameters α and β of the gamma pdf are estimated for each station and for each time scale of interest (1,3, 6, 9, 12 months, etc.). Maximum likelihood estimations of α and β are: (2) $$a=\mathrm{}\frac{1}{4A}\left(1+\mathrm{}\sqrt{1+\mathrm{}\frac{4A}{3}}\right),\hspace{0.5em}\beta =\mathrm{}\frac{\overline{x}}{a},\hspace{0.5em}\mathrm{where}\mathrm{}A=\ln \left(\stackrel{=}{x}\right)-\mathrm{}\frac{\sum \ln \left(x\right)}{n}$$(2) and n is the number of observations.

The resulting parameters are then used to find the cumulative probability of an observed precipitation event for the given month and time scale for the location in question. Since the gamma function is undefined for x = 0 and a precipitation distribution may contain zeros, the cumulative probability becomes: (3) $$H\left(x\right)=q+\left(1-q\right)\mathrm{}G\left(x\right)$$(3) in which q is the probability of zero precipitation and G(x) is the cumulative probability of the incomplete gamma function. If m is the number of zeros in a precipitation time series, then q can be estimated by m/n. The cumulative probability H(x), is then transformed to the standard normal random variable z with mean zero and variance of one (Abramowitz 1964), which is the value of the SPI. (4) $$\mathrm{SPI}=\frac{(X_{ij}-\overline{X})}{\sigma }$$(4)

In the above-mentioned equation, X is seasonal precipitation at the i_th rain gauge and j_th observation, $\overline{X}$ is long-term seasonal mean and $\sigma$ is the standard deviation.

Generally, monthly precipitation is not normally distributed so a transformation is performed such that the derived SPI values follow a normal distribution. The SPI is the number of standard deviations that the observed value would deviate from the long-term mean, for a normally distributed random variable (Tigkas et al. 2015). One interpretation of the resultant values is presented in Table 2.

Table 2. Standard Single Z-Index criteria for categorizing years based on the severity of floods or droughts.

Sr. No	Grades	Single Z-index	Probabilities (%)
1	Extreme flood	Z ≥ 1.645	5
2	Severe flood	1.0367≤ Z < 1.645	10
3	Mild flood	0.5244 < Z < 1.0367	15
4	Normal	−0.5244 ≤ Z ≤ 0.5244	40
5	Mild drought	−1.0367 < Z < −0.5244	15
6	Severe drought	−1.645 < Z < −1.0367	10
7	Extreme drought	Z ≤ −1.645	5

2.2.2. Single Z-Index

China Z-index or single Z-index is a drought index introduced by China Meteorological Administration (CMA) in early 1990s to monitor drought and flood. Z-index was adopted by its used for grading dry/wet (flood) spells as well as their intensity. The single Z-index had been used for drought and flood monitoring in different areas of China (Sun and Gao 2000; Wu et al. 2001). The single Z-index is a statistical measure to identify the frequency of droughts and floods of certain severity that will likely to occur. The seasonal mean rainfall data June to August (JJA) has been used to compute the single Z-index.

The formula for grading dry/wet years is as follows: (5) $$\mathrm{Zij}\mathrm{}=\mathrm{}\frac{6}{C_s}{[\frac{C_{si}}{6}{\mathrm{\varnothing }}_{ij}+1]}^{\frac{1}{3}}-\mathrm{}\frac{6}{C_{si}}+\mathrm{}\frac{C_{si}}{6}$$(5) where Zij is Single Z-Index, i is the timescale, that could be 1, 2, 3,…, 72 months; and j is the current month. The parameter i has not been used in this study because only the Z values for average seasonal rainfall are needed.

Csi denotes the coefficient of skewness and ϕij represents normalized variables and are defined as: (6) $$C_{si}=\mathrm{}\frac{\sum _{i=1}^n(x_{ij}-\mathrm{}\overline{x})}{n{\sigma }_i^3}$$(6) (7) $${\mathrm{\varnothing }}_{ij}=\mathrm{}\frac{x_{ij}-\mathrm{}\overline{x}}{{\sigma }_i}$$(7) where xi represents seasonal rainfall (mm) for j month for i period, stands for mean rainfall (mm), σ is the standard deviation, and n represents the number of years under study.

The severity grades and the corresponding probability of occurrence are presented in Table 2.

2.3. Isolating drought years

For the purpose of isolating the drought years, the Single Z-Index and SPI grading analysis were used. Both the indices were computed to determine drought years and their severity. Afterward, these indices were correlated, and one index had been selected for subsequent analysis. Positive values of indices referred to wet years, whereas negative values stated dry years. Linear regression among the SPI and single Z-index was also carried out. Moreover, Pearson correlation coefficients for Single Z-Index versus SPI were calculated and an excellent relationship (r = 0.96) has been obtained (Figure 3b). Moreover, the Single Z-Index has been selected for further analysis as it has picked larger values in both wet and dry conditions compared to SPI. In addition to this, the SPI index couldn't pick the severity of the year 2002 (a drought year) as considerably captured by the Single Z-Index (Figure 3a). The year 2002 was among the most severe drought years in the region (Hina and Saleem 2019).

Figure 3. (a) Time series of the SPI and Single Z-Index on Pakistan's rainfall. (b) Scatter diagram of the SPI and Z score for Pakistan's rainfall (1981-2017). Positive values indicate wetter than normal, negative values drier than normal.

2.4. Composite analysis

The composite analysis method has been used to discover the spatial patterns of Pakistan's summer rainfall (JJA) throughout the composite drought summers. The single Z-index-derived drought years (severe and extreme) used in this study are 1987, 1991, 2000, 2001, 2002, and 2005. These drought years grouped together and further investigated as a single variable. This technique has been widely used by many scientists all over the World (Adler et al. 2003; Bauer and Del Genio 2006; Thomas and Martin 2007; Nangombe et al. 2018). The composite method has been opted to demonstrate the more pronounced patterns of variables and minimize the number of plots associated with the study (Nangombe et al. 2018). Drought anomalies were calculated and further processed to create composite drought spatial distribution maps.

2.5. Correlation analyses

The statistical relationship among summer rainfall of extreme drought years and meteorological parameters prior to composite drought rainfall season was determined through correlation analysis. The correlation region between 0°S, 80°N and 20°E, 140°W was selected to correlate against Pakistan's JJA rainfall of composite drought years.

Regional meteorological parameters including 2 meter air temperature, SST, SLP, wind variables and geopotential height were sourced from monthly NCEP-NCAR reanalysis data at various levels for the period of 37 years from 1981 to 2017. These meteorological parameters were further subjected to correlation analysis against average JJA rainfall of composite drought years and used to investigate the most correlated parameter and months that could be used as a potential parameter for drought prediction. The geopotential height and wind parameters were investigated at low to middle levels (850,500 and 300 hPa) to assess the depth and strength of convection clouds and circulation systems (Nangombe et al. 2018). The main supposition for this study was to identify the prevailing conditions prior to the rainy season to be used as a drought alert to inform all related stakeholders before heading.

3. Results

3.1. Pakistan rainfall variability

3.1.1. Inter-annual rainfall structure

Pakistan receives less precipitation compared to the other South Asian countries; hence it is predominantly sensitive to rainfall variability. For that reason, several studies have been conducted to understand the variability and spatiotemporal features of Pakistan's rainfall (Adnan et al. 2016; Hussain and Lee 2016; Latif et al. 2017).

The seasonal climate of Pakistan is based on four symmetric divisions of a year: winter (December–February or DJF), spring (March–May or MAM), summer (June–August or JJA), and fall (September–November or SON). The annual average precipitation ranges between 30 and 400 mm from south to north with a maximum gradient of 900-1800 mm in sub-mountainous area (Adnan et al. 2016; Arshad et al. 2021). Almost 57% of the annual rainfall is due to summer monsoon rains; a dominant hydro-meteorological resource for Pakistan. Summer monsoon rain fall contributes 57%, winter precipitation contributes 30%, and remaining 13% is received by tropical cyclones (Hussain and Lee 2009; Saeed et al. 2011). The Southern Himalayan Mountains receive 200 mm/month precipitation during the summer monsoon and winter precipitation in the form of snow. Melting of this accumulated snow and ice during summer keeps the rivers perennial throughout the year. Precipitation caused by western disturbances in winter provides critical mass to glaciers and snowpacks (Moors and Stoffel 2013). The country's remaining extents receive very little precipitation (20–30 mm/month) during the year (Farooqi et al., 2005).

To give insight into Pakistan's summer precipitation patterns, rainfall analysis was carried outfrom the monthly rainfall data of 30 weather observatories (1981 to 2017) and rainfall spatial distribution map was produced. Figure 1(a) displays the spatial rainfall distribution map of Pakistan and from this map, it is witnessed that northern highlands of Pakistan receive more rainfall than southern areas. Moreover, the precipitation anomalies were also calculated for the same time-period and presented in Figure 1(b). The positive anomalies can be observed in core monsoon regions of Pakistan (central and some northern areas), whereas in the southern areas, the strong negative anomaly has been witnessed. Most of the southern and some of the country's extreme northern parts were most hit by droughts with precipitation anomaly of −30 mm.

On interannual timescale, the upper tropospheric cyclonic or anticyclonic anomalies over West Asian region are closely associated to less or more rainfall over Pakistan and India (Saeed et al. 2011; Hong and Lu 2016; Latif et al. 2017; Yadav 2017). These tropospheric cyclonic or anticyclonic anomalies are closely linked with the two teleconnection patterns over the Eurasian continent (Zhou et al. 2019) including circumglobal teleconnection (CGT) or Silk Road Pattern (SRP) (Lu et al. 2002; Ding and Wang 2005). The CGT and summer North Atlantic Oscillation (SNAO) affect the summer rainfall in North of Pakistan and the Northwest of India (Ding and Wang 2005; Syed et al. 2012). The positive phases of CGT and SNAO are related to excessive rainfall during the summer season over Pakistan. The water vapour transportation (Arabian Sea to South Asia) associated with the monsoon trough, is critical for dry and wet conditions over Pakistan during summer season (Bollasina and Nigam 2011).

Some studies have revealed that these tropical and extratropical interactions have a key role in causing extreme rainfall events in Pakistan and South Asia (Hong et al. 2011; Vellore et al. 2016).

3.1.2. Summer rainfall variability

Pakistan, situated at the northernmost edge of the Asian monsoonal region, experiences intricate climate variability in summer rainfall (Lu et al. 2020). The time-period from June to September is categorized as the summer monsoon season and yields about 60% of total rainfall in Pakistan (Hussain and Lee 2016). Summer precipitation in Pakistan demonstrates strong interannual variability, mostly related to the seasonal rainfall. According to Hussain and Lee (2016), the coefficient of variance for summer rainfall is 26% for Pakistan.

Lu et al. (2020), investigated the interannual variation of monthly summer precipitation (June to September) in Pakistan and concluded that July and August are the months of pronounced rainfall compared to June and September. Moreover, the monthly rainfall patterns (JJA) are similar in the upper and lower troposphere (Lu et al. 2020).

The time series of standardized total summer rainfall over Pakistan for JJA is presented in Figure 2(a) for the study period of 37 years (1981- 2017). The standardized rainfall value of 1.5 is considered as the threshold value for extreme above-normal and extreme below-normal rainfall.

According to this threshold value, the years 1994 and 2014 are categorized as the extreme above-normal rainfall years, whereas 1987, 1991, and 2002 fall into the extreme below-normal rainfall rank.

The periodic variability of Pakistan's summer rainfall can be better characterized as quasi-periodic oscillation with non-random rainfall characteristics on decadal basis, as shown in Figure 2(a). This compares well with the past regional study of dry and wet years in Pakistan (1981-2017), which had also shown that dry years were dominated in the first two decades, whereas the last two decades displayed the dominance of wet years (Hina and Saleem 2019). This oscillation could be understood by closely inspecting Figure 2(a), which presented the time series of Pakistan's standardized rainfall for JJA from 1981-2017. On this time series, 11 years running mean has also been calculated. According to this running mean, positive peak nodes have been observed in 1999, 2014 and 2017, whereas negative peak nodes can be seen in 2002 and 2007.

In the first decade of the study period (1981-1990), Pakistan receives below and much below rainfall (1982, 1985, and 1987) whereas the succeeding decade had 5 years of above-normal rainfall among which only the year (1994) has much above-normal rainfall. The third decade has the dominance of below and much below-rainfall, especially in the first three consecutive years (2001, 2002, 2003), which were ENSO years in the region. Many studies had also confirmed this time-period as severe drought (Barlow et al. 2002; Barlow et al. 2016). Finally, in the last decade, the above-normal rainfall dominance (2011, 2012, 2013) can be observed, whereas 2014 had much above-normal rainfall. However, this analysis revealed that how Pakistan's rainfall pattern evolved over nearly 40 years, it is essential to opt for more robust and verified indices that will provide better understanding of moisture deficit. Henceforth, in this study, the Single Z-Index has been used to isolate dry and wet years along with their severity.

3.2. Drought indices for isolating drought years

To isolate the wet and dry years, the Single Z-Index and SPI were statistically computed and compared to select which one is more appropriate for the subsequent analysis. Figure 3(a) presenting a time series of both indices capturing the same drought years (1987, 1991, 2000, 2001, 2005) with almost the same severity. But Single Z-Index well-captured the whole drought cycle in contrast to SPI with more robust results (1987, 1991, 2000, 2001, 2002 and 2005). Figure 3(b) presents a scatter plot analysis between the Single Z-Index and SPI. From this plot, an excellent relationship (r = 0.96) can be investigated. Moreover, it is also observed that the Single Z-Index has picked larger values in both wet and dry conditions. In previous studies, many authors have also compared the robustness of SPI and Single Z-Index and chosen the Single Z-Index (Naresh Kumar et al. 2009; Nangombe et al. 2018). The main reason is, although the SPI picked out the same drought and wet years but it couldn't pick the severity of wet and dry years as robustly as Single Z-Index. SPI's failure to investigate the considerable severity was consistent with the study conducted by Naresh Kumar et al. (2009), who stipulated that SPI underestimated the intensity of wetness and dryness when the rainfall is too high or too low, respectively. In addition to this, the SPI index couldn't pick the severity of the year 2002 (a drought year) as considerably captured by the Single Z-Index (Naresh Kumar et al. 2009). The year 2002 was among the most severe drought years in the region (Hina and Saleem 2019). According to another study conducted by Jain et al. (2015), Single Z-Index can capture the dry conditions that SPI cannot spot (Jain et al. 2015).

3.3. Pakistan's drought years

In this section, the entire study period (1981-2017) has been graded into normal, flood and drought years, using the Single Z-index. The results presented in Table 3 have revealed that 2.7% and 5.4% fall into an extreme flood and extreme drought categories, respectively. Overall, 27.1% of the years are categorized as normal years, 40.5% years as flood years and 32.4% as drought years out of 37 years of the study period. Furthermore, the severe and extreme drought and flood years have equal probability, i.e., 16.2%.

Table 3. Grading of 1981 – 2017 period into flood, normal and drought years based on Pakistan's average summer rainfall and their probabilities.

Grade	Extreme drought	Severe drought	Mild drought	Normal	Mild flood	Severe flood	Extreme flood
Flood/Drought years	1987 2002	1991 2000 2001 2005	1982,1985, 1993,1999, 2004, 2009	1981, 1984, 1986, 1990, 1996, 1997, 1998, 2015, 2016, 2017	1983, 1989, 1992, 2003, 2006, 2008, 2011, 2012, 2013	1988 1995 2007 2010 2014	1994
Total	2	4	6	10	9	5	1
Probability (%)	5.4	10.8	16.2	27.1	24.3	13.5	2.7

From the Table 3, the drought years (years in the severe and extreme drought categories) identified are 1987, 1991, 2000, 2001, 2002, and 2005.

3.4. Precursor conditions to drought years

Investigation of drought precursor conditions is imperative for the prediction of drought events prior to their inception and hence can be used for the early warning system. The changes in climatic parameters occurring few months before the summer drought seasons behave differently and can be used as forecast. Henceforth, it can subsequently be used to prepare against droughts. The precursor period includes the winter months preceding the summer season. For this study, JJA has been chosen as summer season as these are the months receiving maximum rainfall throughout the year. Any fluctuation in rainfall received within these 3 months can subsequently affect the amount of summer season rainfall. The deficit in the summer monsoon rainfall caused moisture stress and impacted the crop yield in the region (Qaiser et al. 2021).

In Pakistan, major crops are grown in two cropping seasons, namely "Kharif" and "Rabi". Kharif, the first sowing season, starts from April to June and is harvested during the months from October to December. The major Kharif crops include rice, maize, soybean, cotton, sugarcane, jowar and bajra. The second sowing season Rabi, began in October to December and is harvested from March to May. Wheat, oat, barley, potato, tomato and onion are major Rabi crops (Ain et al. 2020). The summer monsoon improves the moisture stress during the Kharif season and meets the water demand for Rabi crops (Adnan et al., 2018). Hence any variability in the rainfall received during summer can negatively impact the GDP of an agriculture-based economy like Pakistan (Latif and Syed 2016).

The correlation coefficients derived from correlation analyses between Pakistan composite drought JJA rainfall and meteorological parameters precursor period are presented in Table 4.

Table 4. Correlation coefficients of Pakistan JJA composite droughts and selected meteorological parameters. The correlations are significant at 95%confidence level. All the values with asterisks are considered to be significant at 95% confidence interval. Shaded cells are chosen for further analysis.

	OND	NDJ	DJF	JFM	FMA	MAM
Air Temperature (2 m)	0.60	0.01	−0.01	−0.32*	0.12	0.50
SLP	−0.37	−0.81	−0.85	−0.71*	−0.53*	−0.58*
SST	0.83*	0.60	0.53	0.70	0.74	0.72
Geopotential Height (500 hPa)	0.06	0.06	−0.09	−0.49*	−0.25	0.33
U wind (500 hPa)	0.01	0.55	0.71	0.52	0.52	−0.68*
V wind (500 hPa)	−0.51*	0.26	0.01	−0.34	−0.45*	−0.56*

The precursor period is defined as 3-month averages of OND, NDJ, DJF, JFM, FMA, MAM representing Oct to November, November to January, December to February, January to March, February to April and March to May, respectively. Based on extreme and severe drought grading, the years 1987, 1991, 2000, 2001, 2002, and 2005 were identified as the composite drought years.

The regional meteorological parameters, including air temperature, SST, SLP, wind variables and geopotential height, were averaged over the domain 0°S, 80°N and 20°E, 140°W to correlate against Pakistan's JJA rainfall of composite drought years. The specified domain (0°S, 80°N; 20°E, 140°W) for this study covers both land oceans that have a direct influence on Pakistan's rainfall.

From Table 4, it is depicted that air temperature (2 m), SLP, and geopotential height (500 hPa) have the highest correlations (-0.32, −0.71, and −0.49, respectively) in the JFM period. However, the SLP also had a good correlation (>0.5) during FMA and MAM. SST also revealed a good correlation in all of the periods under study and showed a significant value (0.83) in OND. The correlation coefficients for zonal wind are good in the periods of NDJ, DJF, JFM, FMA, but the highest correlation was found in MAM. Whereas for meridional wind, the period of OND and MAM has shown a significant correlation.

3.4.1. Air temperature

Air temperature is among the various significant physical drivers used to deduce the extreme climatic conditions and is generally negatively related to the rainfall in the South Asian region. Hence, air temperature can be used as a proxy indicator for below-normal rainfall events. Trenberth and Shea (2005) suggested that both the temperature and rainfall records must be interpreted by considering the robust co-variability among them (Trenberth and Shea 2005). Chang (1997) has used temperature to investigate the teleconnections among droughts and spatial patterns that might be responsible for drought event mechanism during summer (Chang 1997).

The anomalous patterns of 2-meter air temperature during JFM indicate the lowest negative anomalies over most regions over Pakistan and India, as depicted in Figure 4. The JFM period has been chosen for composite analysis as it has a significant correlation as compared to other periods. However, the OND and MAM periods also have a good correlation. The results also have shown that during the JFM period, there is remarkably cooler than normal weather in the south Asian regions, including Pakistan. These low temperatures can inhibit the convection activity and henceforth disturb normal rainfall in the region. In addition to this, another major driver (ENSO) is responsible for affecting the summer monsoon rainfall (Azad and Rajeevan 2016; Gadgil et al. 2004; Saleem et al. 2021; Sein et al. 2021). Still, local temperature variations have more significant effect on frequency and intensity of extreme rainfall spells over Indian region. Therefore, the local atmospheric warming can be considered another vital process that influences extremes' characteristics (Sharma and Mujumdar 2019).

Figure 4. January to March (JFM) 2 m air temperature anomaly for composite drought years (units: °C).

3.4.2. Sea level pressure

The interannual variability of summer monsoon in south Asia is associated with large-scale atmospheric and oceanic variability, particularly ENSO and Indian Ocean Dipole (IOD) (Mishra, 2020a; Fan et al. 2017; Mishra, 2020a). ENSO influences rainfall deficits over South Asia; an earlier signal of highly positive SLP anomalies has shown a connection to droughts in the sub-region.

The spatial patterns of JFM composite mean sea level pressure anomalies for drought years are presented in Figure 5. The spatial pattern portrays a substantial high-pressure system (positive anomalies) over the Himalayan region in the north of Pakistan. In contrast to northern region, there is a negative anomaly trend over some parts of Pakistan's southern region.

Figure 5. January to March (JFM) mean sea level pressure anomaly for composite drought years (units: Pa).

The monsoon trough (extended low-pressure zone over the Indo-Gangetic Plains from the northwest Indian region up to the Bay of Bengal) controls the region's precipitation distribution. The substantial increasing trend in SLP from northwest India up to the Bay of Bengal suggests a weakening of monsoon trough resulting in subdued precipitation activity over northern parts of the India region (Preethi et al. 2017).

Moreover, a very high negative SLP anomaly can be seen over the East Asian region (55-65 N, 40-70 E). North Pacific Subtropical High (NPSH) also plays a significant role in the year-to-year variation of the interface between East and South Asian monsoons (SASM) (Cao et al. 2012). The anomalous eastward (westward) location of this interface is associated with a strong East Asian monsoon and a weak South Asian Monsoon (vice -versa).

3.4.3. Sea surface temperature (SST)

SST is one of the most important oceanic climate variables that can be employed to interpret the interaction among various oceanic and atmospheric phenomena (Wmo, n.d; Bojinski et al. 2014; Reynolds et al. 2007). Various studies have reported that variations in SST have triggered drought events in different parts of the World (Fan et al. 2017; Rajagopalan et al. 2000; Lu and Delworth 2005; Feng et al. 2008; Annamalai et al. 2013; Tefera et al. 2020).

ENSO is the predominant phenomenon that can influence the global climate, especially South Asian summer monsoon. Several previous studies have addressed the relationship among ENSO and SASM (Grove 1998; Krishnamurthy and Goswami 2000; Fan et al. 2010). Given the fact that ENSO influences SASM through the Walker circulation (east–west displacement), the variation in spatial distribution of SST anomalies results in the transformation of the position and strength of Walker circulation anomalies, that subsequently alters the spatial configurations of precipitation anomalies over the South and Southeast Asia (Roxy, Ritika, Terray, Masson 2015; Roxy, Ritika, Terray, Murtugudde, et al. 2015). This might also elucidate why climate scientists nowadays rely on the status of ENSO for seasonal rainfall prediction (Nangombe et al. 2018).

In this section, the SST anomalies have been presented for the period of OND, derived from the Table 4. The trends in SST anomalies during composite dry years over the Indo-Pacific region have been displayed in Figure 6. The red contours display positive SST anomalies designating warm pools whereas blue contours show negative anomalies (cool pools). The OND SST anomaly shows an increasing trend and a warm pool in the Arabian Sea and Bay of Bengal. The warming trend over Arabian Sea has a high magnitude (exceeding 0.3 °C) compared to the Bay of Bengal. This suggests the development of convection in the Arabian Sea, contrary to the Bay of Bengal that might be a reason for the precipitation deficit and consequently cause severe droughts in the Indian subcontinent (Preethi et al. 2017). The trends of SST anomalies reported here agree with the findings of Barlow et al. (2016). According to Barlow, SST anomalies (cold) in the central and eastern Pacific and pattern of SST anomalies (warm) in the western Pacific and eastern Indian Ocean are relative to severe droughts. Moreover, the large-scale climate variability, predominantly La Niña, seems to play a vital role in regional droughts, including the most severe drought episode of the last 50 years (1999–2001) (Barlow et al. 2016).

Figure 6. October to December (OND) mean sea level pressure anomaly for composite drought years (units: °C).

Annamalai et al. (2013) also suggested that the increasing SST trend over the West Pacific has altered atmospheric circulation over the Bay of Bengal that led to additional cool and dry air from the northeast region, resulted in higher rainfall over the west Pacific and reduced over South Asia (Annamalai et al. 2013).

3.4.4. Geopotential height

To investigate the role of geopotential height for droughts' prediction, the geopotential height anomalies of the extreme and severe dry years (1987, 1991, 2000, 2001, 2002, and 2005) have been calculated from Jan to Mar (JFM). The spatial trends of geopotential height anomalies at the middle and low levels (850, 500 & 300 hPa) are investigated, but the geopotential height anomaly trend at 850 hPa tends to be weak and hence not presented in the paper.

Geopotential height anomalies for composite dry years at 500 hPa and 300 hPa are presented in Figure 7a and b. A strong positive anomaly with a magnitude of 5 gpm, has been seen over the north of Pakistan with an eastward extension at both 500 and 300 hPa levels (Figure 7). Anticyclones or positive height anomalies are associated with fine, dry, and settled weather, consequently producing a descending air motion (Vellore et al. 2016). This anticyclone over northern Pakistan fetched cold, dry air from the northeast that resulted in dry conditions. Moreover, there is a strong negative high anomaly (< −5 gpm) with the cyclonic flow that can be observed over most parts of Pakistan, which in turn transport moisture out of the region and cause dry weather. Henceforth, the presence of dipole-like structure having positive/negative anomalies over the north/east of Pakistan might be associated with drought events.

Figure 7. January to March (JFM) geopotential height anomalies at 500 hPa (a) and 300 hPa (b) for composite drought years (units: gpm).

The results of this section are well consistent with the studies conducted by Qaiser et al. (2021), Ahmed et al. (2020) and Ain et al. (2020).

3.4.5. Wind circulation

In this section, the wind anomaly has been derived from zonal and meridional (U & V-wind) wind components in the lower and middle levels (300 hPa, 500 hPa & 850 hPa), but the result of low-level analysis (850 hPa) is not shown here (weak patterns). For both of the wind components (U-wind & V-wind), correlation analysis has been done for different periods, as presented in Table 4. The U-wind component has a good correlation with Pakistan's JJA rainfall in all periods (NDJ, DJF, JFM, FMA & MAM) except OND. In contrast, the V-wind component has a good correlation in OND, FMA & MAM periods. The MAM period is selected to acquire wind vector for this study, as this is the period where both of the wind components have shown a significant correlation.

Figure 8 shows March to May (MAM) wind vector anomalies at 500 hPa (a) and 300 hPa (b) for composite drought years. The outflow of wind vectors from Pakistan is observed at both levels (500 hPa & 300 hPa) that has been caused by strong wind anomalies (clockwise circulation) over the region. Such an outward flow of wind vectors from the Indian subcontinent, inhibit moisture transport and cause a dry condition in the region. Such a circulation system indicates enhanced sinking motions associated with precipitation suppuration and, consequently, led to drought conditions. The above-mentioned results are consistent with some previous studies (Ain et al. 2020; Qaiser et al. 2021).

Figure 8. March to May (MAM) wind vector anomalies at 500 hPa (a) and 300 hPa (b) for composite drought years (units: m² s^-2).

4. Discussion

Understanding drought risk in an agricultural country like Pakistan is vital to combat water scarcity, food security and economic challenges. Although literature reports many studies of drought characterization in the region of Pakistan (Xie et al. 2013; Ahmed et al. 2016; Haroon et al. 2016; Adnan et al., 2018) yet it is still unclear how the occurrence of extreme droughts in the observed climate will change under the warming climate. The extreme droughts in Pakistan had posed serious threats to all spectrum of lives, whereas little is known about their characteristics and drivers under changing climate.

Pakistan is situated at the northernmost edge of the Asian monsoonal region and experiences intricate climate variability in summer rainfall (Lu et al. 2020). The time-period from June to September is categorized as the summer monsoon season and yields about 57% of total rainfall in Pakistan (Hussain and Lee 2016). The spatial rainfall distribution map of Pakistan has shown more rainfall over northern highlands as compared to southern areas (Figure 1). The positive Rainfall anomalies (Figure 2) can be observed in core monsoon regions of Pakistan (central and some northern area) whereas a strong negative anomaly has been witnessed in southern regions.

Most of the southern and some of the country's extreme northern parts were most hit by droughts with precipitation anomaly of −30 mm. In view of global warming, changes in precipitation patterns due to large-scale atmospheric phenomena including South Asian summer monsoon and western disturbances might have caused severe droughts in Pakistan. The monsoon precipitation has been significantly increased in the north of Pakistan (Hanif et al. 2013; Ahmed et al. 2018), whereas a substantial decline has been reported over some western and southern parts of the country that are already declared as drought-prone areas (Salma et al. 2012; Hina and Saleem 2019).

In this paper, we tried to explore the possible cycles and precursor conditions of upcoming droughts in Pakistan. Therefore, in this study, first, we detected the severity and recurrence of drought events in Pakistan, then statistical relationship among summer rainfall of extreme drought years and meteorological parameters (air temperature, SST, SLP, wind variables and geopotential height) prior to rainfall season had been calculated, and finally, drought precursors in Pakistan were figured out to anticipate the occurrence of droughts ahead of time. In order to isolate and categorize the drought years based on severity, SPI and single Z-Index were calculated, and the single Z-index was selected for subsequent analysis. The key reason is that the single Z-Index can capture the dry conditions that cannot be spotted by SPI (Jain et al. 2015). The severe and extreme drought years identified in this study are 1987, 1991, 2000, 2001, 2002, and 2005 (Table 3). Similar drought years were also found in various other Pakistan studies (Hina and Saleem 2019; Ain et al. 2020).

The statistical relationship among summer rainfall of extreme drought years and meteorological parameters prior to rainfall season was determined through correlation analysis (Table 4) and used to investigate the most correlated parameter and months that could be used as a potential parameter for drought prediction. The drought years obtained from single Z-Index have been grouped together and further investigated as a single variable (Adler et al. 2003; Bauer and Del Genio 2006; Thomas and Martin 2007). The region between 0°S, 80°N and 20°E, 140°W was selected to correlate regional meteorological parameters against Pakistan's JJA rainfall of composite drought years. Drought anomalies were calculated and further processed to create composite drought spatial distribution maps.

It was supposed that the atmospheric circulation patterns and associated weather systems few months before the drought summer seasons behave differently from those before the normal or wet seasons. Hence, climatic fluctuations occurring before the dry summer season can be taken into account to forecast droughts well in advance. The results demonstrated that 2 m air temperature sea level pressure, and geopotential height had the highest correlations in JFM period. In contrast, sea surface temperature and the wind had significant correlations in OND and MAM, respectively.

Air temperature is considered as a significant physical driver for extreme climatic conditions and generally negatively correlated with South Asian rainfall; it hence might be used as proxy indicator for the onset of drought in some cases. This study established that temperature status in JFM period in Pakistan can be employed to predict the nature of the upcoming season. The anomalous air temperature patterns in JFM (Figure 5) indicated remarkably cooler than normal weather in the south Asian regions, including Pakistan. These low temperatures can inhibit the convection activity and henceforth disturb normal rainfall in the region. In addition to this, the local atmospheric warming can be considered as another vital process that influences the character of weather extremes (Sharma and Mujumdar 2019).

The analysis of sea level pressure in JFM period is presented in Figure 6. The spatial pattern portrays a substantial high-pressure system (positive anomalies) over the Himalayan region in the north of Pakistan while low-pressure system (negative anomalies) over some parts of Pakistan's southern region. ENSO influences rainfall deficits over South Asia, and an earlier signal of positive SLP anomalies has shown a connection to droughts in the sub-region (Mishra, 2020a). The substantial increasing trend in SLP from northwest India up till the Bay of Bengal suggests a weakening of monsoon trough compares well with the study conducted by Preethi et al. (2017) and could result in subdued precipitation activity over northern parts of the India region (Preethi et al. 2017).

SST is a key climate variable that defines the interaction among oceanic and atmospheric phenomena (Wmo, n.d; Bojinski et al. 2014), and SST variations had triggered many drought events worldwide (Fan et al. 2017; Tefera et al. 2020). The OND SST anomalies (Figure 7) showed increasing trend and a warm pool in the Arabian Sea and Bay of Bengal which proposes development of convection in Arabian Sea that might be a reason of precipitation deficit which consequently cause severe droughts in the Indian subcontinent (Preethi et al. 2017). The trends of SST anomalies in this study agree with the study conducted by Barlow et al. (2016). Literature reported that ENSO influences SASM through the Walker circulation (east–west displacement), the changes in SST anomalies can result in the change in the strength and position of Walker circulation anomalies, that subsequently alter the spatial configurations of precipitation anomalies over the South and Southeast Asia (Roxy, Ritika, Terray, Masson 2015; Roxy, Ritika, Terray, Murtugudde, et al. 2015). The composite dry years in this study were under the influence of ENSO, confirming the upshot of ENSO in causing droughts in the study area. From the six drought years (1987, 1991, 2000, 2001, 2002, and 2005), five of them were under the influence of ENSO conditions and only one year (1991) was categorized as neutral year (see Figure 9).

Figure 9. Hot and dry extremes over Pakistan during the core monsoon season (JJA) and associated links with the ENSO during 1981-2017. (a) Three-month (June to August) precipitation anomalies, (b) the single Z-index anomalies, (c) Nino 3.4 SST anomalies during the same time period. Grey vertical bars denote the dry monsoon anomalies associated with the 6 most severe and extreme drought years.

For geopotential height anomaly in JFM (Figure 8a and b), a strong Anticyclone has been seen over the north of Pakistan with an eastward extension. In contrast, a strong negative high anomaly with the cyclonic flow can be observed over most parts of Pakistan, which in turn transport moisture out of the region and cause dry weather. These results are consistent with the findings of previous studies (Qaiser et al. 2021; Ain et al. 2020).

The wind circulation patterns in and around Pakistan in the middle levels can also be used as a precursor of drought. Notably, during the MAM period of composite dry years, the outflow of wind vectors from Pakistan is observed at both levels (500 hPa & 300 hPa) that has been caused by strong wind anomalies (clockwise circulation) over the region. Such an outward flow of wind vectors from the Indian subcontinent inhibits moisture transport and cause a dry condition in the region.

In summary, this study identified the periods at which the patterns of various meteorological parameters could be employed to forecast seasonal changes (weather patterns) a few months in advance (see Figure 10). Hence, the inception of drought events can be assessed.

Figure 10. Schematic diagram of drought severity assessment and drought prediction.

5. Conclusion

Drought is a complex phenomenon that must be investigated more comprehensively under climate change. There are many studies regarding the characterization of droughts in Pakistan and their relation to large-scale atmospheric dynamics. Yet, there is a lack of studies concerning the precursor conditions responsible for the dry (wet) spells in the region. The severity and recurrence of drought events have been investigated, along with its precursors for the region of Pakistan, and the following conclusions can be made:

1. The spatial patterns of Pakistan's annual summer rainfall anomalies have shown positive anomalies over core monsoon regions of Pakistan (central and northern areas), whereas in the southern areas, a strong negative anomaly has been witnessed.

2. The years 1994 and 2014 are categorized as the extreme wet years, whereas 1987, 1991, and 2002 fall into the extreme dry years, ranked on the basis of standardized total summer rainfall threshold value. The Single Z-index also well captures the dry episodes during 1987, 1991, 2000, 2001, 2002 and 2005, and is a good indicator for the identification of dry (wet) spells in the region.

3. The JFM air temperature anomalies of composite dry years have been remarkably cooler than normal weather in the study region that may inhibit the convection activity and consequently disturb the normal rainfall.

4. The spatial pattern of SLP for JFM period portrays a substantial high-pressure system (positive anomalies) over the Himalayan region and a negative anomaly over southern parts of Pakistan. The substantial increasing trend in SLP from northwest India up to the Bay of Bengal suggests a weakening of monsoon trough resulting in subdued precipitation.

5. The OND SST anomaly shows an increasing trend and a warm pool in the Arabian Sea and Bay of Bengal. The warming trend over the Arabian Sea has a high magnitude (exceeding 0.3 °C) compared to the Bay of Bengal. This suggests the development of convection in the Arabian Sea, contrary to the Bay of Bengal that might be a reason for the precipitation deficit and consequently cause severe droughts in the Indian subcontinent.

6. Geopotential height anomalies of JFM for composite dry years show a strong positive anomaly (anticyclone) over the north of Pakistan that can fetch cold, dry air from the northeast and resulted in dry conditions.

7. The outflow of wind vectors (during the MAM period of composite dry years) from Pakistan is observed at both levels (500 hPa & 300 hPa). Such an outward flow of wind vectors from the Indian subcontinent inhibits moisture transport and causes a precipitation suppuration that consequently led to drought conditions.

The present study results, including drought characteristics, trends, periodicity, and associations with large-scale climatic variables, are useful for recognizing the drought mechanism at a regional scale and will offer the scientific basis to develop an early warning system to mitigate and adapt against the drought impacts.

Acknowledgements

The authors gratefully acknowledge the CAS-TWAS President's PhD Fellowship Program. We thank the Pakistan Meteorological Department for providing the datasets for the study. We are also grateful for the selfless guaidance of Professor Dr. Riyu Lu and financial support by The Second Tibetan Plateau Scientific Expedition and Research (STEP) program.

Data availability statement

The monthly reanalysis data sets that support the findings of this study are publically availabl at: https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5

The station data was obtained from Pakista meterological department https://www.pmd.gov.pk/en/

Correction Statement

This article has been republished with minor changes. These changes do not impact the academic content of the article.

Additional information

Funding

This study was supported by the Second Tibetan Plateau Scientific Expedition and Research (STEP) program (grant no. 2019QZKK0102).