Modified Standardized Precipitation Evapotranspiration Index: spatiotemporal analysis of drought

Abstract

Drought monitoring is a complicated issue as it requires multiple meteorological variables to monitor and anticipate drought accurately. Therefore, developing a method that enables researchers, data scientists, and planners to comprehend drought mitigation policies more accurately is essential. In this research, based on the concepts behind the calculation of the Standardized Precipitation Evapotranspiration Index (SPEI), a new drought index is proposed for regional drought monitoring: the Modified Standardized Precipitation Evapotranspiration Index (MSPEI). The potential of the proposed index is based on the estimation of Reference Evapotranspiration (ET_o). Therefore, the Modified Hargreaves-Samani (MHS) equation based on fuzzy logic calibration is used to estimate ET_o. The proposed index is validated on ten meteorological stations in Pakistan at a one-month time scale. Afterward, based on the Pearson correlation, the performance of the proposed index is compared with the commonly used drought index (SPEI). Results showed a significant correlation (r > 0.7) between the quantitative values of MSPEI and SPEI for all ten stations. Moreover, a modified Tjostheims coefficient is used to estimate and test the spatial correlation between SPEI and MSPEI for different drought classes. According to our findings, the association between the SW, ND, ED, EW, MW, and SD patterns of MSPEI and SPI is 0.74, 0.834, 0.673, 0.592, 0.393, and 0.434, respectively. Meanwhile, considering the significance of future drought trend detection, this research is further extended to detect the future trend of MSPEI by using the Hurst index. In accordance with the results, Bahawalnagar, Sialkot, Lahore, Kotli, and Gilgit all have HI values greater than 0.5 (0.63, 0.58, 0.56, 0.55, and 0.53, respectively). In contrast, Muzaffarabad, Skardu, and Jhelum have HI values 0.47, 0.45 and 0.38, respectively; however, HI values of 0.5 are observed at Dera Ismail Khan (DIK) and Islamabad. Therefore, this research provides a basis for developing and enhancing drought hazard characterization, encouraging researchers and policymakers to monitor and forecast regional droughts using a more accurate drought index.

Keywords:

• Drought monitoring
• reference evapotranspiration
• mitigation policies
• tjostheims coefficient
• hurst index

1. Introduction

Drought is one of the most expensive but least understood natural calamities and causes a substantial impact on natural systems, water supplies, socio-economics, and ecosystems (Jiao et al. (2019); Bayissa et al. (2015)). Although there is no universal definition of drought, in its simplest form, drought has been defined as a scarcity of water, which may occur over a prolonged period in several climatic locations. Unlike floods with a clear start and a sudden end, drought can be characterized by slow progress and long-lasting influences (West et al. (2019)).

Generally, drought is divided into four categories: meteorological, hydrological, agricultural, and socio-economic drought. Meteorological drought is caused by a shortage of rainfall, whereas a lack of stream flow, groundwater, or overall water storage causes hydrological drought. The deficit of moisture in the soil is a critical source of agricultural drought. On the other hand, the demand and supply of water and other societal responses can lead to socio-economic droughts. In addition to the types, drought can also be categorized by severity, intensity, duration, and area extent (Shah and Mishra (2020)).

Since the turn of the century, several droughts have occurred in many regions across the globe with profound effects on agriculture, socio-economic development, and the ecosystem, e.g. Australia (2000-2009), USA (2000-2016), China (2007-2012) and Europe (2007-2010). Moreover, its consequences in famines and crop failures cause each year 6-8 billion dollars (on average) of worldwide economic losses (Chao et al. (2016); Cook et al. (2016)).

Over the past few eras, drought has become more dangerous and fatal in many regions of the world, especially in developing countries like Pakistan, such as; substantial water and food insecurities, which leads to economic losses and financial risks. Furthermore, the drought poses enormous drinking and irrigation water challenges and affects Pakistan’s economy, where agriculture provides a living for more than 60% of the population (Jiao et al. (2019)). Therefore, finding an effective drought monitoring system that integrates many aspects (types and characteristics) of drought at the regional and global scales is essential. Various studies (Wang et al. (2016); Wang et al. (2012); Huang et al. (2016)) have been conducted to monitor and forecast the drought effectively. However, the algorithms of each drought indicator vary by zone.

Moreover, drought indicators based on a single meteorological variable are not sufficient to monitor drought precisely because many factors (rainfall, humidity, wind speed, and temperature) are involved in drought (Shah and Mishra (2020)). Initially, Palmer (1965) introduced a drought index, namely, Palmer Drought Severity Index (PDSI), based on rainfall, temperature, latitude, and available water to characterize and quantify the drought. However, later studies showed that PDSI is unable to compare areas due to the lack of a multi-scaling quality. In contrast to PDSI, the Standardized Precipitation Index (SPI), was proposed by McKee et al. (1993). The SPI is amongst the most extensively used drought indicators and can represent the drought classification states over various time intervals. However, an issue was observed in SPI, i-e. SPI considered only precipitation as an input variable, while many other relevant variables such as temperature, wind speed and Reference Evapotranspiration (ET_o), etc., are not considered in manipulation for drought monitoring (Montazeri et al. (2015)); (Arsenault and Brissette (2014)). Parallel to SPI, Tsakiris et al. (2007) proposed Reconnaissance Drought Index (RDI), incorporating ET_o to characterize drought. Similarly, Vicente-Serrano et al. (2010) proposed the Standardized Precipitation Evapotranspiration Index (SPEI), and the multivariate characteristics of SPEI make it more adequate and superior to SPI, which considers only one climatic variable (i.e. only rain). Thus, SPEI has been extensively used in various research fields, such as investigations of drought variability, drought reconstruction, global warming, and the drought’s impact on hydrology, agriculture, and the environment (Beguería et al. 2014). But, in the original proposal of SPEI, Vicente-Serrano et al. (2010) used Thornthwaite (TH) equation to estimate ET_o (i.e. Penman-Monteith equation (Task Committee on Revision of Manual 70 70 (2016)). However, one of the major flaws of the TH equation is that it produces undefined ET_o values in low-temperature areas (Van Wijk and De Vries 1954; Papadopoulou et al. (2003)). Furthermore, the TH equation underestimates ET_o in dry and semi-arid locations but overestimates ET_o in semi-humid and wet areas (Van der Schrier et al. (2011)).

Therefore, to overcome the abovementioned issues, this research intends to develop and evaluate a new index: the modified Standardized Precipitation Evapotranspiration Index (MSPEI). In MSPEI, the Modified Hargreaves-Samani (MHS) equation, proposed by Habeeb et al. (2021), is used to estimate ET_o. First, the parameters of the calibrated Hargreaves-Samani equation (ah, bh, ch, dh) are empirically estimated using the fuzzy logic interference system. However, it is important to note that, in this manuscript, we used averages instead of the difference between two-time steps (i.e. 0000 UTC and 1200 UTC) to calibrate the parameters. Because the differences might be helpful to estimate ET_o for each of the instantaneous hours, but due to the extreme variations in metrological parameters, these instant hours do not accurately represent daily ET_o.Therefore, daily ET_o is recommended in most studies instead of particular hours, and similarly for other variables on which ET_o depends. Thus, considering the importance of daily ET_o, we used the two-time steps; average wind speed, humidity, and temperature in this research.

Thus, the proposed index’s capability and effectiveness are evaluated using Pearson correlation. Further, the modified Tjostheims coefficient is employed to estimate the spatial correlation between different patterns of SPEI and MSPEI. This research used the Lambert projection method in the Tjostheims coefficient to convert geographic coordinates into rectangular coordinates.

Although MSPEI is an excellent measure of meteorological drought monitoring, it cannot predict future drought patterns. In this context, using the Hurst exponent to forecast drought trends could be a useful technique for better analyzing and interpreting MSPEI time series (Tong et al. 2018; Vega et al. 2019; Zhang et al. 2021). Therefore, the Hurst index is used to predict the future spatiotemporal pattern of MSPEI.

Thus, from the above discussion, it can be concluded that several types of drought indicators have been suggested to monitor distinct types of drought (i.e. agriculture, meteorology, and hydrology) for diverse climatic zones. However, due to the instinctive selection of the probability distribution, geological attributes, and index parameters, there is always uncertainty regarding the precise classification of drought under different techniques (Stagge et al. (2016)). Therefore, it is essential to find a technique that assists researchers, data scientists, and planners in understanding drought mitigation policies more precisely, especially at the provincial level. Thus, considering the importance of accurately estimating drought indices, this study focused on developing and evaluating MHS-based drought indices for drought monitoring under different climate conditions. Meanwhile, considering the importance of the spatial correlation between SPEI and MSPEI, the Modified Tjostheims coefficient is used to estimate and test the spatial correlation between SPEI and MSPEI. In addition, Hurst index based on Rescaled Range (R/S) analysis is used to predict the future trend of MSPEI.

2. Materials and methods

2.1. Standardized precipitation evapotranspiration index (SPEI)

Vicente-Serrano et al. (2010) proposed a new multi-scale drought index SPEI, by considering more than one climatic variable. Both the SPEI and the SPI are calculated using the same mathematical expression. SPEI, on the other hand, has one major advantage over SPI: it considers ET_o in the domain. The SPEI is also known as a water balance model. The mathematical expression of SPEI is: (1) $$D_i=P_i-E{T}_{oi}$$(1)

Where D_i, which represents the moisture deficit for month i, indicates the difference between the P_i and ET_oi. P_i represents total monthly precipitation, whereas ET_oi is the estimated amount of monthly evapotranspiration.

2.2. Modified standardized precipitation evapotranspiration index (MSPEI)

2.2.1. Phase 1: Estimation of reference evapotranspiration

2.2.1.1. Modified Hargreaves-Samani equation

The HS equation is simplified by not considering the effects of humidity and wind speed on the ET_o rate. As a result, under extreme weather conditions, the original HS equation cannot generate accurate location estimates (Patel et al. (2015)). Therefore, the Fuzzy Inference System (FIS) is used in this study to compute the updated values of HS parameters (ah, bh, ch) based on the type of location, such as dry, semi-arid, humid, semi-humid, cold, and so on. The input variables are: WS and Hum represent the average of wind speed and relative humidity observed at 1200 UTC and 0000 UTC, respectively, whereas Temp represents the mean of maximum and minimum temperature and sunshine hours. Figure 1 illustrates how FIS is conceptually organized. For example, the MHS equation is shown as: (2) $$E{T}_{o(\mathit{MHS})}=ah\times R_s{\left(T_{\max }-T_{\min }\right)}^{bh}\left(T_{\mathit{eff}}+ch\right)$$(2)

Figure 1. Architecture of Fuzzy Logic Interference.

Where, T_max and T_min are maximum and minimum temperatures, T_eff = dh × ((3 * T_max) – T_min)/2 (dh is a parameter whose value is also determined by FIS.), R_s denotes extra-terrestrial solar radiation, estimated according to Allen et al. (1998). Figure 2 depicts the block diagram for the fuzzy logic-based calibration of the MHS parameters. For a complete description of FIS, see Habeeb et al. (2021).

Figure 2. Schematic diagram of fuzzy based calibration.

2.2.2. Phase 2: Standardization of index

This phase is concerned with the standardization of the index. The step-by-step procedure is described as follows: The modified water balance model is computed for each selected station in the first step, depending on the difference between precipitation and ET_o(MHS). (3) $$M{D}_i=P_i-E{T}_{o(\mathit{MHS})}$$(3)

In the second step, various probability distributions (more specifically, 32 of the most commonly used) are fitted on MD_i, e.g. Gamma, Logistic, Normal, Laplace, and Gumbel distributions. The propagate R package (Spiess et al. (2018)) is used to apply the Kolmogorov-Smirnov test (Hassani and Silva (2015)) and Anderson Darling tests (Ghosh et al. (2016)) for the goodness of measure. The distributions with the lowest values of Akaike Information Criteria (AIC) and Bayesian Information Criteria (BIC) are then chosen in the third stage. Next, the parameters of each well-fitted distribution are estimated using the appropriate parameter estimation method with the lmom R package. The fourth step of the process involves determining the Empirical Cumulative Distribution Function (ECDF) using the estimated parameters of each well-fitted distribution. (4) $$K\left(x\right)=p+(1-p)F\left(x\right)$$(4)

In order to eliminate the impact of zero values in the MD_i, a little modification in the CDF is done in equation (4)(4) $$K\left(x\right)=p+(1-p)F\left(x\right)$$(4) , e.g. in the case of a Gamma distribution. p is the probability of zero value in the time series data MD_i for each station. If h is the total number of zeros in the MD_i sequence, then p can be calculated using h/T. Where, T denoted the total number of observations in the MD_i time series.

Finally, the ECDF of the probability distribution is converted into a standardized normal distribution with parameters 0 and 1. The approximate transformation method which is proposed by Abramowitz and Stegun (Abramowitz et al.(2012)) is used in this research. Here Z_it is expressed as: (5) $${ Z}_{it}=-W+\frac{u_o+u_{1}W+u_2{W}^2}{1+v_{1}W+v_2{W}^2+v_3{W}^3}$$(5)

Here, $$W=\sqrt{\mathrm{l}\mathrm{n}\left(\frac{1}{\left(K\right(r){)}^2}\right)}$$

When, $$0<K\left(r\right)\le 0.5$$

Similarly, (6) $${ Z}_{it}=+W-\frac{u_o+u_{1}W+u_2{W}^2}{1+v_{1}W+v_2{W}^2+v_3{W}^3}$$(6)

And for, $$W=\sqrt{\ln \hspace{0.17em}\left(\frac{1}{(1-K(r){)}^2}\right)}$$

When, $$0.5<K\left(r\right)\le 1$$

Where, $u_o$ = 2.515517, $u_1$ = 0.802853, $u_2$ = 0.010328, $u_3$ = 1.432788, $v_1$ = 1.432788, $v_2$ = 0.985269, $v_3$ = 0.001308.

2.3. Modified Tjostheims coefficient

Tjostheims coefficient is a non-parametric coefficient, initially introduced by Tjøstheim (1978) and then generalized by Hubert et al. (2010). In a spatial statistic term, Tjostheims coefficient is an index of correlation between two spatial sequences. A standard form of the observed value of the Tjostheims index is given by: (7) $${\text{ Λ}}_{\mathit{OBS}}={\sum }_{i=1}^n d\left(l_F\left(i\right),l_G\left(i\right)\right)$$(7)

Where; $d\left(l_F\left(i\right), l_G(i)\right)$ denotes some appropriate measure of spatial distance, ${ \Lambda }_{\mathit{OBS}}\mathrm{ }$ designates the degree to which homogeneous rank values on F and G occupy spatially close geographical locations. In addition, the untied ranks from $1$ to $n$ serve as the definitions of $F$ and $G;$ the position of rank $i$ on $F$ and $i$ on $G$ are indicated by $l_F(i)$ and $l_G(i)$ respectively, $1 \le i \le n.$ In equation 7(7) $${\text{ Λ}}_{\mathit{OBS}}={\sum }_{i=1}^n d\left(l_F\left(i\right),l_G\left(i\right)\right)$$(7) , the location $l_F(i)$ is represented by two coordinates ${(\mathrm{x}}_{\mathrm{Fi}},\mathrm{ }{\mathrm{y}}_{\mathrm{Fi}}$) (latitude and longitude of the location $l_F(i)$), and $l_G(i)$ by ($x_{Gi}, y_{Gi}$) (latitude and longitude of the location $l_G(i)).$ If $d\left(l_F\left(i\right), l_G\left(i\right)\right)$ is the distance between two points and is defined as the square of the Euclidean distance, then; (8) $$d\left(l_F\left(i\right),l_G\left(i\right)\right)={\left(x_{Fi}-x_{Gi}\right)}^2+{\left(y_{Fi}-y_{Gi}\right)}^2$$(8) (9) $$d\left(l_F\left(i\right),l_G\left(i\right)\right)=x_{Fi}^2+x_{Gi}^2+y_{Fi}^2+y_{Gi}^2-2\left(x_{Fi}{x}_i+y_{Fi}{y}_{Gi}\right)$$(9)

Thus, (10) $${\text{ Λ}}_{\mathit{OBS}}=\mathrm{constant}-2{\sum }_{i=1}^n \left(x_{Fi}{x}_{Gi}+y_{Fi}{y}_{Gi}\right)$$(10)

For convenience, Tjostheims confines his discussion to the distance measure that is the active component of the squared Euclidean distance measure, thus; (11) $$d\left(l_F\left(i\right),l_G\left(i\right)\right)=\left(x_{Fi}{x}_{Gi}+y_{Fi}{y}_{Gi}\right)$$(11)

This index establishes a relationship between two variables using measures that resemble the appearance of ${ \Lambda }_{\mathit{OBS}}.$ However, the earth is spherical, and the geographical coordinates are also spherical. Therefore, it is essential to transform the spherical coordinates into rectangular ones pertaining to spherical distance.

In this research, the modified Tjostheims coefficient is employed to estimate the spatial correlation between different patterns of SPEI and MSPEI. In the modified Tjostheims coefficient, the Lambert projection method converts geographic coordinates into rectangular coordinates. One is referred to see (Hussain et al. 2011) for a comprehensive description and computational procedure involved in Lambert projection. SpatialPack of R package (Osorio et al. (2014)) is used to check the spatial correlation between two indices.

2.4. Trend analysis of MSPEI

The Mann-Kendall (M-K) test is commonly used to analyze climatological and hydrological time series trends. It has two benefits: first, as it is a non-parametric test, the data need not to be normally distributed. The other is that it is less sensitive to sudden interruptions due to non-uniform time series. In this test, the null hypothesis assumes that there is no trend and is compared to the alternative, which infers that there is a trend (Rahmat et al. (2012); Zhang et al. (2021); (Rahmat et al. (2012); Zhang et al. (2021)). Thus, the M-K method was selected to test the drought trend in this study. The M-K statistic is denoted as $S\_MK$ calculated in the following manner: (12) $$S\_MK={\sum }_{i=1}^{n-1} {\sum }_{j=i+1}^n \mathit{sign}\left({\mathit{MSPEI}}_j-{\mathit{MSPEI}}_i\right)$$(12)

Where, sign () is a sign function. S_MK has a mean zero, and the variance can be calculated as follows: (13) $$\mathit{Var}\left(S\_MK\right)=\mathrm{ }\sqrt{\frac{n\left(n-1\right)\left(2n+5\right)}{18}}$$(13)

The test statistic Z _MK is therefore written as: (14) $$Z\_MK=\{\begin{array}{c}\frac{S\_MK-1}{\mathrm{Var}(\mathrm{S}\_\mathrm{MK})}\mathrm{ }\mathrm{ }\mathrm{ }S\_MK<0\\ \mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }0\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }S\_MK=0\\ \frac{S\_MK+1}{\mathrm{Var}(S\_MK)}\mathrm{ }\mathrm{ }S\_MK<0\end{array}$$(14)

A positive Z_MK value represents an upward trend, while a negative one indicates a downward trend.

2.5. Future trend of MSPE based on a Hurst index

Hurst Index (HI) is primarily used to assess a long-term memory of a time series. In other words, it determines how much the series deviates from the random walk. In this research, MSPEI's future trend is examined using HI. If HI = 0.5, the time series is a random walk, which means that if a location is presently experiencing drought, it is impossible to predict whether the drought will persist in the future or not. When 0 ≤ HI < 0.5, the time series is inconsistent, meaning that if a region is dry, it will progressively ease and become wet in the future. When 0.5 < HI ≤ 1, it implies that the time series is continuing in a positive direction, meaning that the drought conditions will continue in the same direction as before (Tong et al. 2018); Vega et al. (2019); Zhang et al. (2021).

Algorithm

The calculation of $HI$ includes the following steps:

Step 1: Consider the MSPEI time series and calculate the mean sequence of the MSPEI ($\mathrm{\kappa }$), where $\mathrm{\kappa }$ = 1, 2, 3, … n. (15) $${\overline{\mathrm{MSPEI}}}_{\left(\mathrm{\kappa }\right)}\mathrm{ }=\frac{1}{\mathrm{\kappa }}{\sum }_{\mathrm{k}=1}^{\mathrm{\kappa }} \mathrm{MSPE}{\mathrm{I}}_{\left(\mathrm{k}\right)}\mathrm{ }\mathrm{\kappa }=1,\mathrm{ }2,3,\mathrm{ }\cdots ,\mathrm{n}$$(15)

Step 2: Calculate the cumulative dispersion using the following formula. (16) $${\text{ D}}_{(\mathrm{k},\mathrm{\kappa })}={\sum }_{\mathrm{k}=1}^{\mathrm{\kappa }} \left(\mathrm{MSPE}{\mathrm{I}}_{\left(\mathrm{k}\right)}-{\overline{\mathrm{MSPEI}}}_{\left(\mathrm{\kappa }\right)}\right)\mathrm{ }1⩽\mathrm{k}⩽\mathrm{\kappa }$$(16)

Step 3: The sequence of the Range is computed as follows: (17) $${ R}_{\left(\kappa \right)}=\max \mathrm{ }\left({\mathrm{D}}_{(\mathrm{k},\mathrm{\kappa })}\right)\text{-}\min \mathrm{ }\left({\mathrm{D}}_{(\mathrm{k},\mathrm{\kappa })}\right)$$(17)

Step 4: The following is the standard deviation sequence. (18) $${ S}_{\left(\kappa \right)}={\left[\frac{1}{\kappa }{\sum }_{\mathrm{k}=1}^{\kappa } {\left(\mathit{MSPE}{I}_{\left(\mathrm{k}\right)}-{\overline{\mathrm{MSPEI}}}_{\left(\mathrm{\kappa }\right)}\right)}^2\right]}^{\frac{1}{2}}\kappa =1,2,\dots ,n$$(18)

Step 5: Finally, the Hurst index is calculated as follows. (19) $$\frac{ R_{\left(\kappa \right)}}{S_{\left(\kappa \right)}}=(\beta \kappa {)}^H$$(19)

In the formula, H is the Hurst exponent (0 < H < 1), and $\beta$ is a constant representing the rescaled range at a unified time scale. (20) $${\mathrm{ }\log [R/S]}_{\left(\kappa \right)}=\log \hspace{0.17em}\left(\beta \right)+\mathit{Hlog}\left(\kappa \right)$$(20)

The logarithm of both sides of the equation (20)(20) $${\mathrm{ }\log [R/S]}_{\left(\kappa \right)}=\log \hspace{0.17em}\left(\beta \right)+\mathit{Hlog}\left(\kappa \right)$$(20) is used to generate the empirical Hurst formula. The Modified Hurst Index is a straight-line slope that exposes the fractal properties of a time series, calculated using time series and Hurst’s practical technique to create a cluster of H values for least-squares fitting.

3. Application

Drought is one of the most perplexing and little-understood natural disasters, affecting more people than any other disaster. Like the rest of the globe, Pakistan experiences several water scarcity and contamination issues. In addition, almost every region of Pakistan has experienced severe droughts, causing serious damage to the industrial, agricultural, and economic sectors Hussain and Mumtaz (2014). Furthermore, several human fatalities have been recorded over the last three decades in Tharpakar (Sindh Province, Pakistan). Therefore, time is needed to fortify drought-tracking modules and alleviation strategies by developing robust drought-monitoring appliances and frameworks. Thus, the proposed index is initially applied to ten meteorological stations in Pakistan to assess the proposed framework’s efficiency. For that purpose, monthly precipitation, wind speed (min and max), temperature (min and max), humidity (min and max), and sunlight hours were all needed. Therefore, secondary data on these variables is collected by the Karachi Data Processing Center (KDPC) through the Pakistan Meteorological Department (PMD) from January 1980 to December 2019. This dataset complies with World Meteorological Organization (WMO) criteria, whereas KDPC does error deduction and quality control.

4. Results and discussion

4.1. Assessment of temporal behavior of precipitation

Droughts are caused by a lack of precipitation and an unpredictable precipitation pattern. The monthly precipitation occurrence in various months of the three meteorological stations of Pakistan, Skardu, Muzaffarabad, and Lahore, is presented in Figures 3–5, respectively. The fluctuations in the precipitation patterns of Skardu, Muzaffarabad, and Lahore are caused by the different altitudes of these meteorological stations. To avoid the multiple figures, we presented the precipitation behavior at three stations; however, monthly precipitation occurrence at other stations can be observed similarly.

Figure 3. The pattern of monthly precipitation at Skardu Statio. Different colours in the figures depict different monthly precipitation levels that occurred in the chosen years. The selected years are indicated on the horizontal axis, while the varying months are shown on the vertical axis.

Figure 4. Muzaffarabad station’s monthly precipitation pattern. Different monthly precipitation levels that occurred in the chosen years are depicted in different colours in the figures.

Figure 5. The monthly precipitation behavior for Lahore.

4.2. Estimation and comparative analysis of modified reference evapotranspiration

4.2.1. Empirical estimation of modified Hargreaves-Samani equation parameters

The fuzzy logic inference system is used to empirically estimate the parameters of MHS equation at ten meteorological stations in Pakistan. The fuzzy inference engine generates the fuzzy output variable for the parameters of MHS. Meanwhile, a suitable defuzzification center of gravity technique converts fuzzy variables into distinct outputs. According to Table 1, the mean value of ah varies between 0.0012 and 0.0022, which is lower than the initial value 0.0023), revealing that the ET_o for this site is likely to be overestimated by the initial equation. Furthermore, Bahawalnagar and Sialkot have higher ah values, which is lower in Gilgit.

Table 1. The calibrated Hargreaves-Samani equation’s empirical parameters.

Station	ah	bh	ch	dh
Jhelum	0.0017	0.49	16.8	0.55
Lahore	0.0019	0.45	16.3	0.58
Sialkot	0.0022	0.48	13.7	0.56
Muzaffarabad	0.0015	0.52	21.6	0.44
Kotli	0.0019	0.48	13.4	0.53
Skardu	0.0015	0.58	22.4	0.22
Islamabad	0.0013	0.53	21.1	0.57
Gilgit	0.0012	0.6	13.6	0.51
DIK	0.0014	0.57	16.2	0.59
Bahawalnagar	0.0022	0.45	12.7	0.61

In the same way, the value of bh fluctuates between 0.45 and 0.6 and is higher in five stations (Gilgit, Skardu, DIK, and Islamabad), indicating that the initial HS equation tends to overestimate the ET_o for these stations while being lower in the remaining stations than the initial value 0.5, indicating that the initial HS equation will underestimate the ET_o for the leftover stations. Meanwhile, ch is observed to be varied from 12.7 to 22.4, with higher values in Skardu, Muzaffarabad and Islamabad, and lower values in the remaining stations in contrast to the actual value 17.8, and can be interpreted in the same way as the preceding. The dh is a calibration coefficient of effective temperature whose initial value varies from research to research. Camargo et al. (1999) used dh = 0.72, Pereira and Pruitt (2004) recommended dh = 0.69, while Ahmadi and Fooladmand (2008) showed that the dh value will be changed for different locations. Similarly, in this research, the value of dh is changed according to the station, dh is observed to be varied from 0.22-0.61. Bahawalnagar has the greatest dh value, while Muzaffarabad has the lowest one.

4.2.2. Comparative analysis of modified reference evapotranspiration

This subsection presents a comparative evaluation of the non-calibrated and calibrated HS equation for ten meteorological stations in Pakistan. The original HS and MHS equations were evaluated on the basis of statistical indicators Root Mean Square Error (RMSE), Mean Absolute Error (MAE), and Performance Index (PI) classification. For a complete description of these statistical indicators, see Habeeb et al. (2021).

Table 2 demonstrates the comparison of HS and MHS with PM equation for estimating ET_o. In Jhelum, with a PI value of 0.17, the non-calibrated HS equation has the highest RMSE and MAE values. Therefore, according to the classification criteria stated in Habeeb et al. (2021), the performance of the non-calibrated HS equation is categorized as not good. By using the single parameter (dh) calibrated HS equation, the results were significantly improved, MAE decreased from 257.49 to 106.58, RMSE decreased from 270.4 to 118.12, and PI increased from 0.17 to 0.56. Therefore, the performance of a single parameter calibrated HS equation is improved. Meanwhile, we obtained better results after calibrating two, three, and four parameters of HS equation, demonstrating the effectiveness of fuzzy logic calibration for this study site.

Table 2. Comparison of HS and calibirated HS equations with PM equation for estimating ET_o.

Stations	Equation	MAE	RMSE	PI	Classification
Jhelum	HS	257.49	270.4	0.17	Not Good
	MHS (dh)	106.58	118.12	0.564	Good
	MHS(ah + dh)	48.71	57.86	0.79	Very Good
	MHS(ah + bh + dh)	49.51	57.18	0.8	Very Good
	MHS(ah + bh + ch + dh)	50.32	57	0.81	Very Good
Lahore	HS	150.29	150.93	0.51	Medium
	MHS (dh)	117.24	124.93	0.52	Medium
	MHS(ah + dh)	54.82	68.22	0.71	Very Good
	MHS(ah + bh + dh)	42.67	51.37	0.8	Very Good
	MHS(ah + bh + ch + dh)	42.89	49.96	0.82	Very Good
Sialkot	HS	142.67	143.3	0.55	Medium
	MHS(dh)	114.96	123.54	0.56	Good
	MHS(ah + dh)	97.96	108.26	0.6	Good
	MHS(ah + bh + dh)	79.51	92.17	0.65	Good
	MHS(ah + bh + ch + dh)	50.15	61.67	0.79	Very Good
Muzaffarabad	HS	236.26	253.22	0.23	Not Good
	MHS(dh)	77.83	94.28	0.61	Good
	MHS(ah + dh)	70.01	84.01	0.7	Good
	MHS(ah + bh + dh)	65.39	76.21	0.73	Very Good
	MHS(ah + bh + ch + dh)	60.83	67.81	0.74	Very Good
Kotli	HS	225.32	237.83	0.18	Not Good
	MHS(dh)	110.38	119.9	0.52	Medium
	MHS(ah + dh)	54.64	67.65	0.71	Very Good
	MHS(ah + bh + dh)	48.54	59.28	0.75	Very Good
	MHS(ah + bh + ch + dh)	45	51.53	0.82	Very Good
Islamabad	HS	134.3	135.03	0.56	Good
	MHS (dh)	130.96	140.3	0.49	Medium
	MHS(ah + dh)	66.48	80.36	0.67	Good
	MHS(ah + bh + dh)	57.9	67.56	0.74	Very Good
Gilgit	HS	146.82	163.86	0.51	Medium
	MHS (dh)	121.32	132	0.63	Good
	MHS(ah + dh)	85.51	110.94	0.77	Very Good
	MHS(ah + bh + dh)	65.63	80.81	0.85	Very Good
	MHS(ah + bh + ch + dh)	62.95	71.55	0.87	Excellent
DIK	HS	130.37	132.74	0.61	Good
	MHS (dh)	118.56	133.76	0.51	Medium
	MHS(ah + dh)	79.79	96.38	0.65	Good
	MHS(ah + bh + dh)	71.44	82.31	0.71	Very Good
	MHS(ah + bh + ch + dh)	65.65	72.88	0.75	Very Good
Bahawalnagar	HS	181.4	190.09	0.26	Not Good
	MHS (dh)	146.19	147.46	0.58	Good
	MHS(ah + dh)	129.3	139.57	0.52	Medium
	MHS(ah + bh + dh)	63.03	79.38	0.7	Good
	MHS(ah + bh + ch + dh)	48.71	56.64	0.83	Very Good
Skardu	HS	122.39	123.84	0.7	Good
	MHS (dh)	72.52	84.15	0.61	Good
	MHS(ah + dh)	68.91	84.13	0.68	Good
	MHS(ah + bh + dh)	61.87	69.83	0.72	Very Good
	MHS(ah + bh + ch + dh)	62.28	74.81	0.75	Very Good

The outcomes associated with other stations can also interpret similarly. However, to summarize the results, the non-calibrated HS equation showed good performance in Skardu, DIK, Islamabad, medium in Gilgit, Lahore, and Sialkot, and not good in Bahawalnagar, Jhelum, and Muzaffarabad. The calibration of only dh ranked good for Jhelum, Sialkot, Muzaffarabad, Gilgit, Skardu, and Bahawalnagar and as medium for the leftover stations, making it possible to obtain improved results after the calibration. Similarly, the calibration of four parameters is ranked as good for Gilgit, while other stations are ranked very good. Meanwhile, the results of the study can also be verified by analyzing Figure 6. Figure 6 showed that, in comparison to the non-calibrated HS equation, the calibrated HS equations are closer to PM equation, indicating the certainty that all calibrations, especially four parameters calibrated HS equation promoted better ET_o estimate, and the success of fuzzy logic calibration is the reason behind this.

Figure 6. Comparison of HS and calibrated HS equations with PM equation for estimating ET_o.

4.3. Estimation and comparative analysis of standardized indices

4.3.1. Fitting probability distribution

The appropriate probability distribution for each station is considered when calculating the drought index. Furthermore, the formulation incorporates 32 of the most commonly used probability distributions, and the most appropriate distributions are identified using the R package "propagate", which provides a list of these distributions. Minimum AIC and BIC values are used to determine the suitable probability distribution for each station. In Figures 7 and 8, six meteorological stations of Pakistan, namely Bahawalnagar, Kotli, Gilgit, DIK, Islamabad, and Jhelum, are chosen randomly to illustrate the probability distributions that are most suitable for the SPEI and MSPEI.

Figure 7. Appropriate Probability Distributions of SPEI.

Figure 8. Appropriate Probability Distributions of MSPEI.

4.3.2. Comparative analysis of SPEI and MSPEI

The effectiveness and capability of the proposed index (MSPEI) are measured by including the results of the existing index (SPEI). Therefore, this subsection provides the comparative analysis of SPEI and MSPEI for four randomly selected stations of Pakistan. The temporal deviation between the estimated value of SPEI and MSPEI is illustrated in Figure 9. To avert the multiple figures, the behavior of the four stations is presented here. Thus, a minimal variation has been found between the quantitative values of both methods in all selected stations. Moreover, the comparative analysis of SPEI and MSPEI revealed that although both indices’ patterns are identical, the variation has been slightly reduced in MSPEI. Meanwhile, in Table 3, the significant correlation between the quantitative values of SPEI and MSPEI for the selected stations can be observed.

Figure 9. Comparative analysis of SPEI and MSPEI.

Table 3. Pearson’s product-moment correlation between quantitative values of SPEI and MSPEI.

Stations	Correlation	Test-Statistic	Confidence Interval (95%)
Jhelum	0.885	41.7	0.864-0.903
Lahore	0.887	42.21	0.867-0.905
Sialkot	0.727	23.19	0.682-0.767
Muzaffarabad	0.897	44.41	0.878-0.913
Kotli	0.907	47.16	0.890-0.921
Skardu	0.835	33.3	0.806-0.861
Islamabad	0.890	42.89	0.870-0.908
Gilgit	0.809	30.18	0.776-0.838
DIK	0.682	20.37	0.631-0.727
Bahawalnagar	0.801	29.34	0.767-0.831

4.3.3. Spatial comparative analysis

The results of the comparative spatial analysis are based on Tjostheims coefficient index. Tjostheims coefficient measures the correlation between two stochastic processes observed over space. In this research, the Modified Tjostheims coefficient is employed to estimate the spatial correlation between different patterns of SPEI and MSPEI. For that purpose, the quantitative values of MSPEI are classified as Extremely Wet (EW), Sever Wet (SW), Moderate Wet (MW), Normal Drought (ND), Moderate Drought (MD), Sever Drought (SD) and Extreme Drought (ED) by using Oliver’s classification criteria, given in Table 4. Hence, the total counts for SPEI and MSPEI are given in Tables 5 and 6, respectively. The overall results for this index demonstrated that the frequencies of MD in both indices are high as compared to other patterns for all stations. Thus, outcomes associated with this research show little deviation in the counts of annual variability patterns determined by MSPEI. In addition, this section also describes the results of the comparative spatial analysis of EW, SW, MW, ND, MD, SD, and ED classes under both indices (MSPEI and SPEI). Table 7, demonstrated that there is a perfect spatial correlation between MSPEI and SPEI for the ND (r = 0.834), for SW and the MD spatial correlation is also high (r = 0.740, r = 0.718), for EW and ED spatial correlation is moderate, r = 0.673 and r = 0.592 respectively, whereas, for SD and MW, spatial correlation is slightly low (r = 0.434, r = 0.393).

Table 4. Drought classification criteria of SPEI and MSPEI.

SPEI/MSPEI Values	Classes
≥ 2	EW
1.5 to 1.99	SW
1 to 1.49	MW
0.99 to −0.99	ND
−1 to −1.49	MD
−1.5 to −1.99	SD
≤ −2	ED

Table 5. Total counts for SPEI of EW, SW, MW., ND, MD, SD and ED patterns.

Station	EW	SW	MW	ND	MD	SD	ED
Jhelum	3	9	33	264	70	61	40
Lahore	2	6	14	299	58	45	56
Sialkot	6	10	17	283	70	64	30
Muzaffarabad	5	9	15	283	68	56	44
Kotli	7	11	30	250	87	68	27
Skardu	2	4	18	267	55	82	52
Islamabad	6	11	30	270	68	73	22
Gilgit	0	0	11	281	52	78	58
DIK	4	4	12	286	59	50	65
Bahawalnagar	5	6	10	280	67	59	53

Table 6. Total counts for MSPEI of EW, SW, MW, ND, MD, SD and ED patterns.

Station	EW	SW	MW	ND	MD	SD	ED
Jhelum	1	7	38	256	69	69	40
Lahore	4	9	18	272	73	63	41
Sialkot	4	4	8	288	84	47	45
Muzaffarabad	4	9	12	272	65	74	44
Kotli	1	12	25	262	71	62	47
Skardu	0	6	28	244	78	86	38
Islamabad	1	10	27	258	73	74	37
Gilgit	3	4	17	278	78	60	40
DIK	6	7	14	265	84	65	39
Bahawalnagar	5	7	17	276	67	57	51

Table 7. Summary statistics of modified tjostheims coefficient.

Index	Pattern	Spatial Correlation
SPEI vs MSPEI	EW	0.592
	SW	0.740
	MW	0.393
	ND	0.834
	MD	0.718
	SD	0.434
	ED	0.673

4.4. Trend analysis

This subsection presents the current and future trend analysis of meteorological drought. The M-K test is first applied to the MSPEI time series on a one-month time scale. The Z_MK statistic obtained from M-K test for selected stations in Pakistan is depicted in Table 8. A Positive value of Z_MK indicates an upward trend and whereas a downward trend is indicated by a negative value of Z_ MK. In addition, the null hypothesis is tested at 5% level of significance. From Table 8, it can be demonstrated that only station Muzaffarabad showed a significantly increasing (wet) trend with the Z_MK value 1.69, Islamabad showed a trend that only slightly increased with Z_MK value 0.60, while the remaining stations showed a significantly decreasing (drying) trend. In Muzaffarabad, the minimum value of MSPEI is −1.413 in October 2010, which was a moderate drought, and the maximum value was 2.336 in July 1980, which was extremely wet. We can also interpret the results for other stations in the same way.

Table 8. Summary of the Mann–Kendall (M.K.) test of MSPEI.

Stations	Years	Min MSPEI	Max MSPEI	Z_MK
Muzaffarabad	1980-2019	−1.413 (Oct,2010)	2.336 (July,1980)	1.69
Bahawalnagar	1980-2019	−2.511 (Nov,1995)	3.620 (Sep,2012)	−2.20
Islamabad	1980-2019	−2.422 (May,1985)	2.243 (Jul,1998)	0.60
Sialkot	1980-2019	−2.516 (May,1993)	2.283 (Aug,2011)	−1.77
Lahore	1980-2019	−2.496 (May,1988)	3.221 (Aug,1999)	−1.10
Skardu	1980-2019	−2.474 (Jul,1997)	1.933 (Feb,2013)	−1.57
Jehlum	1980-2019	−2.516 (May,1984)	2.282 (Jul,2000)	−0.87
Gilgit	1980-2019	−2.516 (Jun,1985)	2.682 (Aug,2013)	−1.41
Kotli	1980-2019	−2.516 (May,1994)	2.244 (Jul,1991)	−1.36
DIK	1980-2019	−2.449 (Jun,1985)	5.731 (Aug,2013)	−1.27

In addition, considering the significance of future drought trend detection, the R/S analysis method is used to calculate the HI values of MSPEI for selected stations of Pakistan from 1980-2019 and to predict the trend of MSPEI. The result of the HI is demonstrated in Table 9. According to the statistics, the HI value of MSPEI in Bahawalnagar Sialkot, Lahore, Kotli, and Gilgit is 0.63, 0.58, 0.56, 0.55, and 0.53, respectively. In particular, HI> 0.5 indicates that present and future MSPEI tend to persist, at least in the near future. In other words, the MSPEI will continue in the same direction as before. On the other hand, Muzaffarabad, Skardu, and Jhelum had HI values of 0.47, 0.45, and 0.38, respectively, indicating moderate anti-persistence with the trend from 1980-2019. Meanwhile, two stations DIK and Islamabad, both have HI = 0.5, indicating that the MSPEI of DIK and Islamabad are a random walk. Therefore, in these two stations, it is not certain whether the trend of MSPEI in the future is consistent or anti-consistent with the current state.

Table 9. R/S analysis table for MSPEI 1980–2019.

Stations	Regression Equation	HS
Muzaffarabad	y = −0.0834 + 0.47x	0.47
Bahawalnagar	y = −0.1321 + 0.63x	0.63
Islamabad	y = −0.0653 + 0.50x	0.50
Sialkot	y=-0.1230 + 0.58x	0.58
Lahore	y = −0.1502 + 0.56x	0.56
Skardu	y = −0.0537 + 0.45x	0.45
Jhelum	y = −0.2301 + 0.38x	0.38
Gilgit	y = −0.0265 + 0.53x	0.53
Kotli	y = −0.0667 + 0.55x	0.55
DIK	y = −0.0151 + 0.50x	0.50

4.5. Discussion

This research provides an assessment of meteorological drought by using the modified drought index, MSPEI. The results of this study were compared with the previous studies conducted on meteorological drought monitoring, such as; Vicente-Serrano et al. (2010) proposed SPEI, and the multivariate characteristics of SPEI make it more adequate and superior to SPI. But, in SPEI Vicente-Serrano et al. (2010) used TH equation to estimate ET_o. However, the TH equation underestimates ET_o in dry and semi-arid locations but overestimates ET_o in semi-humid and wet areas (Van der Schrier et al. (2011)). Therefore, this study focused on developing and evaluating a novel drought index for drought monitoring under different climate conditions to overcome the aforementioned issue. In MSPEI, we used a fuzzy logic-based calibrated HS equation to estimate ET_o, see, Habeeb et al. (2021). Thus, the efficiency and capability of the proposed index (MSPEI) are measured by including the results of the existing index (SPEI). Our analysis shows a substantial positive correlation between SPEI and MSPEI at each station. Moreover, the comparative analysis of SPEI and MSPEI (Figure 9) revealed that although both indices’ patterns are identical, the variation has been slightly reduced in MSPEI. Therefore, it can be concluded that the SPEI's inadequacies can be addressed by employing the MSPEI.

Subsequently, we used M-K test to investigate spatial and temporal variations in Pakistan’s drought characteristics. Out of ten stations, only Muzaffarabad showed a significantly increasing trend or wet trend, and Islamabad showed a semi-drying trend. In contrast, the remaining stations indicated a drying trend due to the significant decrease in precipitation and a considerable increase in temperature occurring in the southern part of Pakistan, leading to severe drought. Moreover, according to Slop Difference (SD) classification (Nabeel and Athar 2018), 66% of Pakistan’s land is classified as arid, whereas 30% (4%) as semi-arid (humid). Meanwhile, Hurst exponent predicted the future drought trend of Pakistan. Statistics showed that the present and future MSPEI of 5 out of 10 stations tend to sustain, at least shortly; on the other hand, three stations indicated mild anti-persistence with the trend in 1980-2019, whereas HI values of DIK and Islamabad indicated that their MSPEI's are a random walk. However, this research is unable to determine how long the anticipated drought trend will continue in the future, nor can it forecast future drought patterns in various categories.

In summary, our experimental results indicate that the performance of MSPEI was very well in monitoring the meteorological drought of the study area. Moreover, the outcomes of this research may accurately and timely inform decision-makers to anticipate sound strategies and plans to mitigate the adverse effects of drought. But, despite all these advantages, some limitations could affect the performance of MSPEI. First of all, there is no absolute "true" drought measure. We established a benchmark for comparisons by utilizing SPEI as the reference but cannot provide a true validation of MSPEI. Secondly, MSPEI was used to monitor drought only on a monthly scale. But the severity and spatial distribution of drought is continuously changing. This requires a smaller time scale, such as 8 days or 16 days. Thirdly, although there is a strong positive correlation between each quantitative value of SPI and SPEI, it does not guarantee that each indicator produces the same drought class in a particular month. Fourthly, the trend statistics presented in this research can be different if the trend analysis was conducted over a longer or shorter period. In order to develop a more credible region-specific drought trend, future studies should investigate the turning or change point in the drought trend.

5. Conclusion

In order to avert economic losses and natural disasters, it is essential to monitor the drought effectively. Various studies (Palmer (1965); McKee et al. (1993); Tsakiris et al. (2007); Vicente-Serrano et al. (2010)) have suggested different indices to effectively monitor and forecast the drought. However, each drought index’s methodology differed by zone. Therefore, in this research, we introduced a new drought index based on MHS for regional drought monitoring, MSPEI, and then, spatio-temporal analysis of drought is performed. From the analysis of ten meteorological stations of Pakistan, we came to the following conclusions:

• The effectiveness and capability of the proposed index (MSPEI) are measured by including the existing index (SPEI) results, and significant correlation between the quantitative values of MSPEI and SPEI has been found.

• Therefore, we have concluded that the inadequacies of SPEI can be addressed by using this modified index MSPEI. Consequently, MSPEI can be used to develop effective drought mitigation policies.

• Modified Tjostheims coefficient showed a strong positive spatial correlation between SPEI and MSPEI in SW, ND, and MD for ED and EW spatial correlation is also positive but slightly low, whereas the spatial correlation for SD and MW is low.

• Hurst index for Bahawalnagar, Sialkot, Lahore, Kotli and Gilgit showed a strongly persistence long term trend, whereas Muzaffarabad, Skardu and Jhelum indicated moderate anti-persistence with the trend of MSPEI in 1980-2019. HI values of DIK and Islamabad indicated that their MSPEI's are random walk.

• However, the results of HI are unable to foretell how long the predicted drought trend would last in the future or how future drought patterns will behave in different categories.

Concerning future work, the inferences and mathematical calculations can be generalized for other time scales. Moreover, for a more precise estimate of the ET_o, it is recommended to utilize the actual Rs data instead of relying on the estimates.

Acknowledgments

Their authors extend Their appreciation to the Deanship of Scientific Research at King Khalid University for funding this work through Larg Groups. (Project under grant number (RGP.2/23/44). And this study is supported via funding from Prince Sattam bin Abdulaziz University project number (PSAU/2023/R/1444).

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

Data will be provided upon request.