Influence of large-scale circulation and interannual modes of climate variability on seasonal drought characteristics over Pakistan

The long-term drought monitoring and its assessment are of great importance for meteorological disaster risk management. The recurrent spells of heatwaves and droughts have severely affected the environmental conditions worldwide, including Pakistan. The present work sought to investigate the spatiotemporal changes in drought characteristics over Pakistan during Rabi and Kharif cropping seasons. The role of large-scale circulation and interannual mode of climate variability is further explored to identify the physical mechanisms associated with droughts in the region. Monthly precipitation and temperature data (1983-2019) from 53 meteorological stations were used to study drought characteristics, using the standardized precipitation evapotranspiration index (SPEI). The non-parametric Mann-Kendall (MK), Sen’s Slope (SS), and Sequential Mann-Kendall (SQMK) tests were applied on the drought index to determine the statistical significance and magnitude of the historical trend. The state-of-the-art Bayesian Dynamic Linear (BDL) model was further used to analyze large scale climate drivers of droughts, analysis revealed an increase in drought severity mostly over arid to semi-arid regions for both cropping seasons. Temperature played a significant role in defining droughts over dry and hot seasons, while rainfall is influential over the western disturbances (WD) influenced region. The analysis of atmospheric circulation patterns revealed that large-scale changes in wind speed, air temperature, relative humidity and geopotential height anomalies are the likely drivers of droughts in the region. We found that Niño4, Sea Surface Temperature (SST), and multivariate El Niño-Southern Oscillation (ENSO4.0) Index are the most influential factors for seasonal droughts across Pakistan. Overall, the findings provide a better understanding of drought-prone areas in the region, and this information is of potential use for mitigating and managing drought risks.


Introduction
Drought is recognized as one of the most severe hydrometeorological disasters, caused by prolonged dry spells due to the below-average rainfall amounts (Dai 2013;Yang et al. 2020).
Several studies have recently reported an increase in drought occurrence because of shifting climatic patterns (i.e., rainfall and temperature) (Duffy et al. 2015;Ahmed et al. 2018;Hui-Mean et al. 2018;Himayoun and Roshni 2019). The frequency and intensity of drought is projected to increase in many regions by the end of the 21 st century (Dai 2013). It is essential to study the climate-drought relationship for risk reduction and climate change adaptation (Ahmed et al. 2018). A recent study has reported the global warming-induced soil moisture deficit during summer and spring seasons in semi-arid regions, leading to greater drought risk vulnerability (Wu et al. 2018). The rising temperature patterns with below-average rainfall could decrease soil moisture contents over arid to semi-arid regions, making it more vulnerable to drought risks (Himayoun and Roshni 2019;Ain et al. 2020).
Recently, Pakistan has been severely impacted by extreme events such as heatwaves, floods, and droughts (Adnan et al. 2016;Ahmed et al. 2018;Ain et al. 2020). However, the country is situated in a predominantly arid region where the temperature is extremely high and precipitation is low (Haroon et al. 2016). The economy is highly dependent on agriculture, which plays a significant role in economic growth (Ashraf and Routray 2015). In winter and summer seasons, Pakistan receives enough amount of rainfall, while early spring season the moisture is transported from the Mediterranean and Caspian Sea developing western weather system (Ahmed et al. 2018;Anjum et al. 2018). Drought is a well-known natural hazard in Pakistan, which occurs at least every four out of ten years (Ahmed et al. 2018). The agriculture sector contributes about 25% to Gross Domestic Product (GDP) and employs approximately 45% of the total labor-workers (Ullah et al. 2017;Ain et al. 2020).
Most of the country population lives in rural areas and depend on agriculture production as their primary income source.
Numerous studies have been conducted in Pakistan to identify the spatiotemporal characteristics of droughts on annual and seasonal timescales as well as long term trends in these characteristics (Ahmed et al. 2018;Hina and Saleem 2019;Adnan et al. 2020;Naz et al. 2020). A recent study has reported that the seasonal drought exhibits a cumulative trend in south-eastern Pakistan, while a declining trend is evident in the northwest coastal belt (Naz et al. 2020). Moreover, the drought events in Rabi and Kharif seasons have shown an upward trend in the south-eastern Pakistan, while a downward trend in the central region (Ahmed et al. 2018). The interannual variability of seasonal droughts over Pakistan is projected to increase in the near future, which may lead to adverse effects on the environment and society (Chatterjee et al. 2016;You et al. 2017). The studies mentioned above have provided detailed information about the spatiotemporal changes in annual and seasonal droughts characteristics over Pakistan, however, none of the studies have reported the large-scale climate drivers of the seasonal drought variability.
The El Niño Southern Oscillation (ENSO) is the primary mode of interannual variability in the tropical Pacific with significant impacts on global weather and climate via atmospheric circulations. Given this, efforts have been made to develop different forecasting methods to predict ENSO evolution several seasons in advance (Endris et al. 2019;Kang and Lee 2019).
It is well documented in the literature that interannual rainfall variability over Pakistan is influenced by several dominant large-scale modes of climate variability (Ahmed et al. 2018;Sajjad and Ghaffar 2019), including ENSO new Version-4 (ENSO4.0) defined by the sea surface temperature (SST) anomalies in the central and eastern equatorial Pacific Ocean (McBride and Nicholls 1983), the Indian Ocean Dipole (IOD) defined by SST gradients in the west-eastern tropical Indian Ocean (Saji and Yamagata 2003), and the Southern Annular Mode (SAM), which refers to the atmospheric circulation in the mid to high latitudes of the southern hemisphere on inter-annual timescales (Marshall 2003).
On longer timescales, the temperature of Pakistan has been increasing due to anthropogenic climate change, yet there have been no significant changes in precipitation.
For instance, the mean surface temperature has increased by 0.1℃/decade during 1960-2010 over Pakistan (IPCC 2014;Khan et al. 2018). ENSOs is one of the dominant factors of climate variability in the tropical Indian Ocean, showing a possible influence on Pakistan's rainfall variability (Ain et al. 2020). In contrast, most of the tropical and semi-arid regions suffer from droughts during El Niño events. The phenomena of drought are widespread in the arid and semi-arid regions during La Niña (Bates et al. 2008;Ahmed et al. 2018). In addition, the SST based (Niño3.4-SST, and Niño4-SST) in the south-eastern Indian Ocean may also cause heavy rainfall (Ali et al. 2019). Most recently, Das et al. (2020), and Mishra (2020) found that Niño4-SSTs are the most influential factor in triggering heavy precipitation over the southeasterm Indian Ocean in near future. This is, in part, because the negative phase of IOD is associated with La Niña events showing a dipole pattern of SST with warming (dry) and cooling (wet) in the east-western Indian Ocean (Mishra et al. , 2021Ehsan et al. 2020). It is anticipated that the complex interactions between internal modes of variability and climate change has, and will continoue, to change Pakistan's drought characteristics in space and time .
Though the aforementioned studies have investigated the seasonal drought characteristics and its trend over Pakistan; however, several aspects of the seasonal droughts and its associated mechanisms are not sufficiently understood. The current study explores the underlying physical mechanisms of seasonal droughts and different modes of climate variability in the region (Pakistan) to fill these gaps. We used state-of-the-art Bayesian Dynamic Linear (BDL) model to study the influences of large-scale climate indices on seasonal droughts over Pakistan. The finding of the study would provide better understanding of the seasonal drought characteristics over Pakistan, essential for improving the drought risk reduction and management strategies. The rest of the paper is organized as follows: Section 2 provides the details about the study area while Section 3 presents datasets and methodologies.
The results of the investigation are shown in Section 4. Section 5 elucidates discussions and conclusions based on significant outcomes.

Study area
Geographically, Pakistan (23°-37° and 60°-78°) is located in central-south Asia, covering a total area of 881,913 km 2 (Fig. 1). The country receives maximum rainfall in the summer season, accounting for more than 65% of the annual total precipitation (Waqas and Athar 2019). The mean annual maximum and minimum temperatures are 15 to 35℃ and 0 to 20℃ across the country. Pakistan has two main seasons: summer and winter monsoon; the summer monsoon precipitation spans four months from June to September, and winter precipitation extents over December to March. The westerly disturbance in the Mediterranean Sea caused by winter monsoon precipitation, while the monsoon winds mostly cause the summer precipitation originated from the Bay of Bengal (Ahmed et al. 2018;Khan et al. 2020). The Rabi (winter) and Kharif (summer) are the two distinct cropping seasons in the region. Kharif season starts in May and ends in October, and is characterized by soil moisture plentiful enough to support the rain-fed crops. The main Kharif season crops include rice, maize, sugarcane, and cotton. Similarly, the Rabi season fluctuates with respect to latitude and is mostly influenced by the monsoon departure. Furthermore, the Rabi season begins in November and ends in April and the Rabi season crops include peas, gram, wheat, and mustard (Miyan 2015;Adnan et al. 2020).

In-situ observations
The monthly rainfall and temperature datasets (1983-2019) of 53 meteorological stations were archived from the Pakistan Meteorological Department (PMD). A brief summary of the selected meteorological stations is presented in Table S1. The monthly values were summed up to make the seasonal and annual time series of temperature and precipitation datasets. The seasonal drought was categorized as follows: Kharif season (May to October) and Rabi season (November to April). Moreover, the seasonal mean precipitation and temperature were calculated by averaging the mean rainfall and temperature of the selected months. The standardized normalized homogeneity test (SNHT) is among the commonly used statistical technique for estimating the uniformity of climatic data records; therefore, the main purpose of SNHT and homogenization assessment is to find out the peaks and outliers in datasets attributed by non-climatic factors (Toreti et al. 2011;Akhtar and Athar 2019). During the climatic signal's preservation, only (< 2%) of the data was missing (a few stations) that do not substantially affect the outcomes of the study to a broader range.

Reanalysis dataset
Recently, the European Centre for Medium-Range Weather Forecasts (ECMWF) developed a global reanalysis dataset, namely ERA-5, which was used in the present study. ERA-5 is the 5th global atmospheric reanalysis product developed by ECMWF following FGGE, ERA-15, ERA-40, and ERA-Interim, with key strengths and advances over ERA-Interim. Large-scale treatment of precipitation and clouds, the inclusion of analytical variables for precipitating and snow, model boundary layer conditions, and radioactive data are among the key aspects of ERA-5 dataset (Hersbach et al. 2019). The output appeared to be better with these changes and improvements, representing the large-scale tropospheric circulations since 1979 (Olauson 2018). In the present study, meridional wind speed, soil moisture, geopotential height 850 hPa, total precipitation, 2m relative humidity 850 hPa, and 2m air temperature were used.

Large scale indices
Climate change is not the only factor affecting droughts severity in the region; therefore, interannual modes of climate variability, defined by common large-scale indices, are further investigated (Gao et al. 2017;Yang et al. 2019Yang et al. , 2020. These large-scale indices include oceanic SST fluctuation, atmospheric circulation, and air pressure fluctuation, etc. In recent decades, many studies have been conducted at regional and global scale associated with the ENSO4.0, IOD, and monsoon variability towards extreme climate events (i.e., floods and droughts) conditions in the monsoon region (Gao et al. 2017;Joshi and Kar 2018;Yang et al. 2020). However, in the region of Pakistan, the westernmost boundary of South Asian monsoon system is not explored in terms of precipitation and temperature variability. This study chooses seven indices as a potential driver, representing the large-scale mode of climate variability. Detailed information on each index is given in Table 1. The monthly data (1983-2019) for each index were archived from (https://psl.noaa.gov/data/climateindices/list/).

Standardized precipitation evapotranspiration index (SPEI)
The SPEI considers potential evapotranspiration (PET) in addition to precipitation for evaluating drought conditions (Vicente-Serrano et al. 2010), and therefore found to be more efficient in constructing the temporal variability of droughts with respect to climate change perspective (Himayoun and Roshni 2019;Adnan et al. 2020). Furthermore, in SPEI the collected monthly precipitation and PET differences at various timescales are fitted with a 3parameter probability distribution function (PDF). The best fitted PDF of those estimated parameters are used to calculate the SPEI values. The detailed information about SPEI is given by Vicente-Serrano et al. (2010). The state of drought is classified by SPEI values as: mild drought (SPEI ≤ -1.0), moderate drought (-1.5 < SPEI ≤ -1.0), severe drought (-2.0 < SPEI ≤ -1.5 ), and extreme drought refer to (SPEI ≤ -2.0), respectively. In this study, we used the 6-month SPEI values for seasonal analysis and 12-month SPEI values were adopted as the annual timescale.

Trend descriptive statistics
In current study, the rank-based non-parametric modified-MK (Kendall 1955), Sequential Mann-kendall (SQMK) (Sneyers 1990), and Sen's Slope (SS) estimator (Sen 1968) were applied to detect secular trends in auto-correlated SPEI data (Ahmed et al. 2018;Yang et al. 2020). The presence of autocorrelation in the data could affect the probability of trend detection. Hamed and Rao (1998) suggested that the trend estimator should be subtracted from the original time series to estimate the autocorrelation. After that, modified-MK test was employed to detect the significance of the monotonic trend in SPEI, temperature and precipitation (Ahmed et al. 2018;Yang and Zhang 2018;Yu et al. 2018). Moreover, the modified-MK test was further used for SPEI values and secular trends with -value < 0.05 were considered as statistically significant. Likewise, the SS estimator was applied to determine the slope of the trend in SPEI time series (Ain et al. 2020;Adnan et al. 2020). In addition, the SQMK test was used to perceive the abrupt changes in SPEI time series at seasonal and annual timescales over time (Himayoun and Roshni 2019). The aforesaid statistical approaches are briefly described in previous research studies (Ahmed et al. 2018;Himayoun and Roshni 2019;Yang et al. 2020). The mathematical equations of modified-MK test statistic (S) is given below as: where shows the length of time series, and are the sequential data values at time observation and , and denotes the symbolic function that accepts the values 1, 0, or -1; In Eq. 3, ( ) is the variance of whereas, the computed standardized ( ) values with variance "1" and average "0" follow normal distribution time series data, while it is applied to measure the significance of the trend(Bayer Altın and Barak 2017).
where, is the data points and Φ −1 indicate the inverse standard normal distribution, while the correlation matrix of the time series was calculated from the Hurst coefficient ( ) (Koutsoyiannis 2003). is the rank of equivalent normal variants of the de-trended series.
( ) = [ | − | ], = 1: ; = 1: In Eq. 5 and 6, , , indicates the length of time series and show autocorrelation function of lag for a given . The significance level of was defined by the mean and standard deviation of = 0.5 . In addition, to remove bias in ( ) ′ estimation we used bias correction factor. Where indicate the function of , while calculating the significance of modified-MK we replace ( ) with ( ) .
For SS estimation (Eq. 9), the changes in time series are calculated as the median of slopes computed from two consecutive points of the series.
For linear regression (Eq. 10), represents the precipitation and temperature trend; shows the input precipitation and temperature over time; while and are indicated the regression coefficients, and intercepts, accordingly.
The SQMK test is calculated through applying ranked values ( ) of the data time series such as 1 , 2 … . . . during analysis. The magnitude of the value, ( = 1, 2 …. ) are calculated with as the value for ( = 1 …. -1). However, for assessment, each case where values greater than are calculated and represented by . A mathematical statistic ( ) can be defined as, (Ahmed et al. 2018;Chen et al. 2019) in Eq. (11).
The calculation of test statistics , has a mean as in Eq. (12), For variance, the statistical test computed as in Eq. (13), For mitigation and standardized variable of the sequential values, statistical test ( ) is computed for each of the test statistic variables as follow in Eq. (14).
In Eq. 14, the null hypothesis values will be accepted at the significant level, if | ( )| ≤ ( ) 1− /2 , where ( ) 1− /2 of the critical value of standard normal distribution whose probability value exceeds /2 and where value set as 0.05. Moreover, the decrease values of UF show as a negative trend in time series and vice versa (Zhong et al. 2018;Zhang et al. 2019).

Bayesian dynamic linear model
Bayesian where indicate response variable with respect to drought index, is the covariate in term of climate pattern, and represent the dynamic intercept and slope coefficient at time , respectively, whereas the error sequences , , , and , are independent, both within them and between them.

Statistical validation approach
Furthermore, the Pearson's correlation coefficient ( ), normalized Root Mean Square Error  Based on the evaluation metrics shown in Table S2, the seven large-scale climate indices were applied as a potential driver for drought analysis over Pakistan. The overall statistical validation for those climate indices is in line with earlier studies (Ain et al. 2020;Yang et al. 2020;Feng et al. 2020;Shah and Mishra 2020).

Jointly seasonal droughts return period
The seasonal drought is determined by estimating the SPEI for a time-ratio equivalent to the last months seasonal time span of the corresponding season. We used 6-month SPEI index for the analysis of Rabi and Kharif seasons drought over Pakistan. Seasonal maxima series of SPEI were employed to calculate the return period in seasonal SPEI values by using frequency analysis. In our case, we considered SPEI values ≤ -0.5 as a drought, and years with no droughts are replaced by zero. We preferred drought frequency analysis based on non-zero values and non-exceedance probability ( ′ ) that is adopted in account of zero values, where indicates non-exceedance probability value attained by applying frequency analysis, and is the probability of zero values computed as the ratio of the number of time intervals without drought occurrences to the total number of recording period (Ahmed et al. 2018;Mohsenipour et al. 2018  The annual SPEI analysis depicts increasing and decreasing patterns of drought severity across Pakistan during 1983-2019 (Fig. 2c). The south-eastern and central parts of the country exhibited the highest drought severity with an amount of -0.18 to -0.12 over arid to semi-arid regions of the country. For instance, significant decrease in annual SPEI is evident over northern and south-western parts of the country, which is dominated by westerlies monsoon regions. It can be inferred from the findings that the drought severity would rise in the region due to global warming, where the larger risk is associated with arid and semi-arid regions (Fig. S1).

Spatial distribution of annual and seasonal drought trends
The spatial distribution of the SPEI trends during Rabi season, Kharif season, and annual timescale is shown in Fig. 3. The light-yellow dots in the figure depict the significance of increasing and decreasing trend at the 95% confidence level. In Rabi season (Fig. 3a), the negative trend is obvious with a significant rise in drought severity over the major portion of the country (increase in negative SPEI), however, it is observed that a significant increase in Rabi SPEI is dominant over arid and semi-arid regions, which are highly vulnerable During the annual SPEI trend analysis (Fig. 3c), the result indicates an increasing  1994, 2003, and 2014 (1984, 1989, and 2010). In the Kharif season (Fig. 4b), a relatively slight decrease in trend as depicted by a regression coefficient of -0.84 with dry and wet years appeared, possibly causing upward (downward) trends. The driest (wettest) years of the Kharif season were 1984, 1994, and 2003 (1986, 1995, and 2011). On an annual timescale (Fig. 4c), the smallest decreasing trend is evident with a regression coefficient of -0.86 for the study region. Moreover, the highest (lowest) annual dry and wet years have been observed during 1984, 1991, 1994, 2003, and 2010 (1986, 1988, 1995, 2004, and 2014) in the target region.
From Fig. 4, it can be elucidating that regional variability in seasonal drought (Rabi and Kharif) over Pakistan is obvious during 1983-2019. The extent of regional drought variability could be associated with the geographical location, monsoon circulation patterns, and/or large-scale climate indices. The temporal trend over the Asian and Indian monsoon system showed a strong variation at different time scale; from monthly, inter-annual, seasonal, and annual time scale.

Mutations in temporal trends at the annual and seasonal scale
To document the abrupt changes in seasonal temporal trends and annual SPEI time series, the secular SQMK trend test was applied (Fig. 5). The progressive and retrograde series were In Kharif season (Fig. 5b), the trend appeared to be showing a decreasing pattern during the initial study period, however, a rapid rise in the middle years has been observed.
The overall trend increased with significant abrupt shifts that appeared during 1991-2001 with two turning points during 1995 and 2016. In the annual trend (Fig. 5c), the abrupt negative changes were detected for the first decade, whereas, a rapid increase in the middle and later years suggesting an abrupt increasing trend. The annual trend is increasing with a significant increase that appeared during 1991-2001 with two turning points during 1997 and 2014, respectively.

Composite analysis
To assess the seasonal drought variability and its relationship with large-scale climate indices a composite analysis was performed. To do so, the regional scale monthly, and annual 6month robust SPEI index was developed at the interannual scale. The years with deviation above/below ±1 were chosen as dry and wet. A total of 8 years each were selected as dry and wet years for seasonal and annual timescales over Pakistan 1983-2019 as shown in Table S3.
Composites anomalies of geopotential height (850 hPa), wind speed, 2m air temperature, and relative humidity (850 hPa) for seasonal and the annual timescale over Pakistan are given in Fig. 6 and S2. During Rabi dry years (Fig. 6a), a relatively weak subtropical jet is observed over the target region, while WD shows a dominant pattern, exhibiting a negative anomaly during dry years. A similar pattern can be seen for 2m air temperature and relative humidity, further clarifying dry conditions with elevated air temperature anomalies due to weak precipitation (Fig. S2). This negative anomaly is associated with strong east-western westerlies, transporting moisture from the Bay of Bengal, and the Arabian Sea trajectory to the northern hemisphere during Rabi dry years. A similar pattern can be seen for Rabi wet years (Fig. 6b), suggesting a negative anomaly during wet years. It is observed that the anticyclone westerly pattern (Tibetan Plateau) plays a significant role in transporting moisture incursion from the Bay of Bengal, which enhances the wet conditions in the region. We find that the development of the significance moisture flux convergence is mainly associated with the Tibetan Plateau's anticyclonic pattern, causing wet conditions over the study region.
During Kharif dry years (Fig. 6c), the subtropical and easterly jets appeared to be pushing the oceanic moisture-laid winds towards the continental landmass. The anomalous dry pattern is evident over South Asia (SA) western parts that engulfs Pakistan with enhanced wind speed, underlying dry conditions. On the other hand, the anomalous geopotential height is obvious over western SA that includes Pakistan, with enhanced relative humidity and reduced air temperature (Fig. S2), indicating increased precipitation. The anticyclonic pattern appeared to be moved to Bay of Bengal region, strengthening sub-tropical jets intensity further into continental land masses. A similar pattern can be seen for Kharif wet years (Fig.   6d), suggesting an increased anomalous pattern over the Tibetan Plateau and foothills of Himalayan, which is mainly responsible for south-western summer monsoon in the region.
Likewise, Fig. 7b represents the variations between wind speed and geopotential height at 850 hPa over SA, indicating less amount of wind speed over northern Pakistan, while increasing wind speed is obvious in the southern parts. The increased geopotential height at 850 hPa is observed over the foothill's mountainous region with an amount of 6-8gpm, while the constrained geopotential height at 850 hPa was detected in the south-western parts of Pakistan, which fluctuating from 8 to 10 gpm, respectively. Such uncertainties of increase/decrease in geopotential height at 850 hPa, and low/high wind speed may be linked with the increasing (decreasing) trends of extremely high temperatures over the target region.
Similarly, Fig. 7c indicated that the soil moisture is decreased over the northwest parts of the study region, while an increase in soil moisture was found in the southern parts.
Notably, a decreasing trend of soil moisture is evident over the target region, suggesting a lack of precipitation and higher temperature intensity in the study area. Moreover, Fig. 7d shows that the relative humidity at 850 hPa in the southeast part of Pakistan is obviously higher than that in the southeast and central parts of the SA. This instability could possibly associate with the increased frequency and intensity of temperature extremes over the study region.
The drought return period further clarifies the results obtained in the previous section for Rabi and Kharif seasons over Pakistan. We used the average moving window for calculating return periods over the analysis period of 1983-2019 (Fig. 8). The return periods of drought during seasonal and annual timescales are found in the range of 2-5 (mild), 6-10 (moderate), 11-20 (severe), and 21-35 (extremely severe) years, respectively. During Rabi season (Fig. 8a), the relative changes of 2-5 years return period of drought severity are obvious over most stations, however, the north-eastern part of the country exhibits severe to an extreme severe return period of drought severity. Few stations in eastern Pakistan detected a severe 11-20 years return period of drought severity, suggesting a gradual decline in SPEI values in recent years. In Kharif season (Fig. 8b), a fluctuation in the severe drought return period (6-10 and 10-20 years) is dominant over semi-arid and arid regions of the study area, which is historically more prone to droughts hazard. Moreover, the northern and southeastern parts of the region observed extreme severe return periods, indicating that the frequency of severe drought is high during Kharif season. Similarly, during the annual timescale an obvious increase in drought severity over arid region is appeared (Fig. 8c), particularly in few stations of the eastern part exhibit extreme severe return periods. Overall, severe to extreme droughts return periods are found to decrease gradually over time in the majority of the stations. It can be inferred from figure results that increasing temperature in most parts of the country could decrease SPEI values and thus increase the frequency of droughts.

Relationship between seasonal SPEI and large-scale climate indices
We calculated the Pearson's correlation coefficient to obtain the preliminary outlook of the relationship between all the meteorological stations of seasonal SPEI and each selected climate indices over Pakistan (Fig. 9). Overall, DMI and Niño4 were relatively negatively Figure 10 revealed the appraised dynamic regression slope coefficients of seven climate drivers on two distinct seasons of Rabi and Kharif, and annual SPEI. Fig. 10(a) showed that the association between DMI and SPEI is obvious during Rabi season and exhibit positive change before 1987 to negative until 2017. Furthermore, DMI was negatively correlated with SPEI during Kharif season before 1980s (Fig. 10h), while a positive increase was found during 1990s and 2010s. The slope showed a decadal variability for the annual timescale and gradually changed from negative to positive (Fig. 10o). In addition, the DMI slope coefficients tend to approach 0 before 2010s, suggesting that the effects of the DMI weakened during both seasons of Pakistan in recent decades. Pakistan. For annual timescale (Fig. 10p), ENSO4.0 slope coefficient exhibit incline above 0 after 1990s, indicating that ENSO4.0 has a strong influence during annual timescale in recent decades over Pakistan.
For IOD, the slope coefficient of Rabi season varies from positive to negative during the 2000s (Fig. 10c), while a similar pattern can be seen for Kharif season where IOD depicts a negative correlation with SPEI but its influence is decreased from 2000 to 2019 (Fig. 10j).
Both seasons are found under the two troughs circa 2000, implying that IOD's effect was stronger at that time. Similar pattern can also be seen for the annual timescale, further clarifying IOD's effect across Pakistan (Fig. 10q). Our analysis shows that IOD has a relative negative association with the SPEI series, suggesting that a positive IOD attribute to the drier conditions in the region.
Moreover, Niño4 depicts a predominantly negative impact on the SPEI during Rabi season over 2000-2019 (Fig. 10d), while persistent negative and positive behavior can be seen for Kharif season (Fig. 10k). Typically the weak influence of Niño4 have been found for both seasons, whereas, the correlation for Rabi season evident a phase change circa the 2000s, and the effect of the Niño4 for Kharif season was mettle after 2005, inferring a decrease in the strength and could possibly phase change in near future. In addition, the predominately positive slope coefficient in annual timescale (Fig. 10r), SPEI is expected to decline, which could result in drier climate in the southern and western Pakistan.
During PDO, the slope coefficient increased from negative to positive circa of 2000s during Rabi season (Fig. 10e), while the opposite association is found for the Kharif season in 2010s can be seen in Fig. 10 . The relatively weak influence of PDO in Rabi season over the years, even though it remains negative. Nevertheless, the annual timescale's slope coefficient changed the circa of 2010s from negative to the positive phase (Fig. 10s). On the other hand, the PDO warmer phase after the 2000s with respect to a negative slope contributing dryer conditions for Rabi season, though the cool phase (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) depict positive correlation and wetting conditions during Kharif season.
The influence of SST fluctuates from negative to positive during both Rabi and Kharif seasons in 1990s and 2000s, respectively. During Rabi season, a circa appeared in the 1990s (Fig. 10f), suggesting a stronger impact of SST at that period, while the effect is decreasing after 2000s to a minimum point. Moreover, the slope coefficients tend to increase, indicating a emphasizing influence of SST over Kharif season (Fig. 10m). The correlation fluctuates during 1990-2019 for annual timescale and depicts a positive slope after 1990s (Fig. 10t).
The SST is a leading and potential climate driver influencing the summer temperature in southern and western Pakistan, where summer tends to be warmer and drier when SST is positive and vice versa.
Distinct effects of sunspots on SPEI appeared relatively stable and predominantly negative to positive for Rabi season, Kharif season, and annual timescale, however, during Rabi season, the impact changed negative during the late 1980s and 1990s, while an enhanced positive slope were detected during 2010s (Fig. 10g). It can be remarked that the slope coefficient tends to incline after the 2000s for both seasons and annual timescale, inferring that the influence of sunspots has reinforced in the last few decades. Furthermore, sunspots relatively show a negative relationship with the SPEI from 1980s to 2010s at all scales. Further clarifying that increase sunspots could lead to a decline in SPEI, depicting that the periods with high sunspots amount are more likely to occur dry events in the region.

Discussion
This study aims to provide a comprehensive assessment of climate change influence on seasonal drought and its potential drivers across Pakistan, during 1983-2019. The Rabi and Kharif are evaluated as drought seasons in Pakistan, using modified-MK, SQMK, SS, Linear regression, drought return period, correlation coefficient, nRMSE, and BDL model. Overall, the results inferred that most of the stations across Pakistan exhibited a rise in drought severity during Kharif season and annual timescale, especially those located in arid southeastern and northern parts of the region. The increase/decrease in these seasons and overall variability is somehow agreed with the findings of regional studies (Adnan et al. 2016;Gao et al. 2017;Ahmed et al. 2018;Ain et al. 2020). The drought severity and frequency appeared to be shifting from the arid/semi-arid regions towards the country's western humid region. This shift and increase in drought severity altogether can impact the food security with the provision of abundant/less water for rain-fed regions during early/late summer, though the risk of floods/droughts can also be increased and thus need to be seen.
During seasonal and annual timescale drought severity appeared to be somehow complex in reporting the shifts in trends. The Kharif season drought severity is similar to what has been observed at annual timescale, however, such a pattern was relatively weaker during Rabi season. The possible reason could be statistical operations that average/sum the seasonal droughts, and thus such indicators are lost due to a few extreme events in the months with increasing trends or vice versa. The second possible reason could be Asian summer monsoon and/ or WD dominancy during the Rabi season (Rehman et al. 2018;Khan et al. 2020;Adnan et al. 2020;Feng et al. 2020). The positive trend in Kharif season and annual timescale may positively affect agriculture production, water management and built environment in Pakistan. The study results support the findings of the earlier studies at regional and global scales (Ahmed et al. 2018;Guo et al. 2018;Wang et al. 2019;Alamgir et al. 2020;Shah and Mishra 2020). According to Wang et al. (2014), the varying characteristics of severe droughts based on SPEI exhibit significantly increase trend in terms of global drought zones during 1902-2008. Moreover, recent severe climatic changes caused several extreme droughts events over many regions of Pakistan (Adnan et al. 2016;Ali et al. 2017;Ahmed et al. 2018;Hina and Saleem 2019;Ain et al. 2020;Naz et al. 2020 During Rabi season, climate patterns usually exerted more negative than positive impact on the drought variability in terms of dynamic regression slope coefficients after the 1990s. For Kharif season, the majority of the climate indices were negatively associated with the SPEI after the 2000s, when DMI, ENSO4.0, and sunspots exhibited a phase change circa 1980s, respectively. On an annual timescale, all climate indices are positively correlated with SPEI except IOD and PDO after the 1990s. In general, a relative decrease in severe and extreme droughts is noticed across Pakistan in the recent few decades, which can be accredited to the regional influence of global warming.