Evaluation of wet and dry event’s trend and instability based on the meteorological drought index

A temporal imbalance in the water availability, which is consistently below average or more than average rainfall, can lead to extremely dry or wet conditions. This impacts on agricultural yields, water resources and human activities. Weather instabilities and trends of wet/dry events have not yet been explored in Pakistan. In this study, we have two-fold objectives: (1) evaluate the weather instabilities, and (2) the trend of dry/wet events of selected stations of Pakistan. To observe weather instabilities, we used Mean Marginal Hilbert Spectrum (MMHS) and Continuous Wavelet Power Spectrum (CWPS) as meteorological series are mostly non-linear and non-stationary. We used Ensemble Empirical Mode Decomposition (EEMD) for the analysis of temporal characteristics of dry/wet events. We found that all stations are facing severe weather instabilities during the short period of 5 and 10 months using MMHS method and CWPS has shown the weather instabilities during 4 to 32 months of periodicity for all stations. Ultimately, the achieved short-term weather instabilities indicated by MMHS is consistent with CWPS. In summary, these findings might be useful for water resource management and policymakers.


INTRODUCTION
In recent decades, the climate changes tend to increase the occurrence and intensity of natural perils such as heavy rainfall, flood, forest fires and droughts that indicate an increasing trend (Estrela & Vargas, 2012;Kreibich et al., 2017;Masud, Qian & Faramarzi, Modified Soil Water Deficit Index (MSWDI) of the Songnen Plain in China. Their study aimed to investigate the agricultural drought frequency and trend by HHT. Several studies suggested that HHT approach found to be more effective for decomposing the series into several components through EEMD as compared to Wavelet Transform (WT) (Massei & Fournier, 2012;Yang et al., 2019). The Mann-Kendall (MK) trend test and Spearman's Rho (SR) test were used for the trend analysis of agricultural droughts. As a result, residuals from HHT showed the increasing trend for the year 1981 and a decreasing trend for the year 1961year to 1980year . Shahid (2010 utilized the rainfall data  of 17 meteorological stations in Bangladesh. Two non-parametric tests, i.e., MK and Sen.'s slop (SS) tests were used on the SPI to detect a significance and magnitude in wet and dry events in Bangladesh. He found a significant decrease of dry months in monsoon and pre-monsoon.
According to a global vulnerability index report, Pakistan is in the list of the top 10 countries that are highly affected by climate change (Ullah, 2016) as its economy is heavily dependent on the agriculture sector which can severely be affected by climate changes (Ali et al., 2017). In this study, we aimed to evaluate the weather instabilities and the trend of dry/wet events of selected stations of Pakistan.

Study area
The study areas include four meteorological stations (Multan, Bahawalpur, Barkhan and Khanpur) of Pakistan. These selected stations are located in the Southern part of Pakistan, mostly very hot and mildly cold areas. Since our objective is to analyze the trend and weather instabilities with respect to time, therefore, selecting four (stations) time-series data is enough, increasing the number of stations may cause a problem in time series analysis. Monthly quantitative data of precipitation from January 1990 to December 2018 has been obtained from the Karachi Data Processing Center through the Pakistan Meteorological Department, Karachi. The geographical presentation of selected stations is displayed in Fig. 1.

Standardized Precipitation Index at 3-months (SPI-3)
Mckee, Doesken & Kleist (1993) proposed a Standardized Precipitation Index (SPI) for defining and monitoring wet and dry events i.e., beginning, ending and intensity. The SPI is used to measure the precipitation shortage from the long-term historical record of precipitation. In this study, the cumulative precipitation series such as 3 months' time scale is used to calculate SPI. Because, SPI-3 with a short time scale describes droughts, which is important for the agriculture sector (Caloiero, 2018;Zhang et al., 2012;Ali et al., 2017).

Empirical Ensemble Mode Decomposition (EEMD)
Ensemble Empirical Mode Decomposition (EEMD) is an extension of Empirical Mode Decomposition (EMD) (Huang et al., 1998). EEMD was proposed by Wu & Huang (2009) and which is used to handle the mode mixing problem. Extension in EMD procedure is done through adding white noise in the time series. EEMD is a very effective procedure that minimizes the presence of mode mixing problems in the decomposition phase and is very helpful in separating frequency scales. EEMD uses the time series to decompose into several stationary Intrinsic Mode Function (IMF) and a residual component. According to Wu & Huang (2009), the EEMD procedure is as follows: • Set the ensemble number E for each case and select the amplitude of white noise.
• Construction of white noise series η i (t ) based on the amplitude of white noise and add the white noise series η i (t ) into the time series y(t ) to create a new time series y i (t ), i.e., After this, use the traditional EMD method to decompose the new time series the IMFs. The detailed methodology of decomposing series into the IMFs by EMD is described in the Huang et al. (1998).

Continuous Wavelet Power Spectrum (CWPS)
Continuous Wavelet Transformation (CWT) is used to represent the frequency in the time domain and determine the temporal variability of the series. Following (Labat, 2005), the mathematical form of CWT is as follows: where y(t ) represents the time series, t is time and ψ * represent the complex conjugate of the mother wavelet function. The scale factor y is also known as the frequency that is related to the location of wavelets in the frequency domain and it is used for controlling the width of the wavelet function. Whereas, τ the parameter is related to the location of the wavelet in the time domain and it is used for the adjustment of wavelet location. Initially, Morlet wavelet is presented by Grossmann & Morlet, (2009) and its mathematical form is defined as: where γ 0 is a non-dimensional frequency and i is the complex values. (Farge, 1992) suggested that the term γ 0 is the non-dimensional frequency (angular frequency) is set 6 to satisfy the admissible condition. Spectrum construction in the wavelet method is like the Fourier method. But it analyzes spectrum time, frequency and decomposition rather than the time and scale (Lim & Lye, 2002). According to (Torrence & Compo, 1998), the Continuous Wavelet Power Spectrum (CWPS) with scale factor is as follows: (4)

Mean Marginal Hilbert Spectrum (MMHS)
Hilbert Hung Transformation (HHT) was offered by Huang, Long & Shen (1996). HHT is a very powerful tool and used for dealing with stationary and non-linear series. IMFs are obtained through a decomposition procedure and there is no difficulty in applying the Hilbert Transformation (HT) to obtain instantaneous frequency (IF) and instantaneous amplitude (IA). HT is only applicable for IMFs and not applicable for residual which is a monotonic function. Following (Huang et al., 2003), we can compute the Hilbert transform z(t ) of time series y(t ) as: where, P is the Cauchy value and HT be present for all the functions of LP class (Bailey & Titchmarsh, 1938). By definition, the analytical signal R(t ) consists of the complex conjugate of and y(t ) and z(t ) where, where, a(t ) represents the superlative local fit of amplitude ϕ(t ) and phase varying trigonometric function to y(t ) (Huang et al., 2003). IF is obtained by taking the first order derivative of phase t to time t . i.e., The IF is a mono-component function in which only one component represents only one frequency. After completion of HT for each IMF, the following structure of the time series is shown as real part: HT is not applicable for the residual term (Non-IMF), because residual is a constant term.
With an explanation of the HS, it is also very important to define the Mean Marginal Hilbert Spectrum (MMHS). i.e., The MMHS presents the contribution of total amplitude which is corresponding to each frequency value.

Non-parametric trend tests
It is very essential to extract the underlying pattern of hydrological and meteorological time series. There are several techniques for identifying patterns and analyzing hydrological meteorological time series trends. The assumptions of parametric test i.e., normality, independence and linearity are not met in the hydrological or meteorological time series. Some non-parametric techniques for dealing with non-normal and non-linear series are available in the literature, such as Mann-Kendall (MK) and Spearman's Rho (SR) trend test. Mann (1945) proposed a non-parametric test which is used to assess the randomness against the trend in hydrology and climatology. According to MK trend test, the null hypothesis H 0 state that the (y 1, y 2 ,...,y n ), are independent and identically distributed sample (there is no trend) and the alternative hypothesis H 1 states that the distribution of y i and y j are not identical (there is a trend) for all i = j. As a resultant, the positive value of the test statistics is an indication of an increasing trend which is increasing with respect to time. Whereas, the negative value of the test statistics is an indication of a decreasing trend which is decreasing with respect to time. Therefore, the test statistic and their theoretical calculation of MK trend test are described in Hirsch, Slack & Smith (1982).

Spearman's Rho (SR) trend test
Spearman's Rho (SR) trend test is based on the rank of observations. SR trend test is used to identify the absence of a linear and non-linear trend. Following (Yang et al., 2019), standardized test statistics are as follows: Where, Ry yi is the rank of ith observation in the time series. As a resultantly negative value z d is an indication of a decreasing trend and the positive value of ith is an indication of an increasing trend. At the 5% of the level of significance (a), the null hypothesis is rejected (trend exist in the series) if |z d | > t (n−2,1−a/2) .

RESULTS
In this paper, the precipitation time series of four meteorological stations (Multan, Bahawalpur, Barkhan and Khanpur) have been used for explaining the wet and dry events periodicity, variability and transient trend. Initially, three months of aggregated precipitation (P-3) series are used for the computation of SPI-3. In the estimation stage, the suitable probability distribution is selected for all series of four meteorological stations. Hence, the detail of fitting suitable probability distributions on P-3 is described in Table 1.

Trend analysis of wet and dry events
Later, the SPI-3 is decomposed into IMFs by using the EEMD procedure. The decomposition level is based on the total number of observations. There are 348 total observations as the SPI-3 series consists of monthly observations from the year 1990 to 2018. According to Sang et al. (2016), log 2 (348) provides eight levels of decomposition. The ensemble number is set to 100 and the amplitude of white noise is about 0.2 standard deviation has been added in SPI-3 for construction of IMFs by EEMD (Wu & Huang, 2009). Finally, by applying the EEMD method, SPI-3 of all stations is decomposed into eight IMFs and a residual term, which is depicted in Fig. 2.  Fig. 2A. Figure 2B shows the residual stability of a normal wet trend before the year 1998. Figure 2C shows the trend of normal wet is increasing from the year 2007. Figure 2D shows the trend of normal wet is decreasing from the year 1990 to 2007, but the trend is gradually started increasing from 2015 to onwards.
The  The significance of IMFs is tested by correlation approach (Peng, Tse & Chu, 2005). Initially, the correlation is computed between the IMFs and SPI-3 for all stations, then the threshold is applied. The correlation approach indicates that all IMFs by EEMD are statistically significant. The reconstruction of SPI-3 is based on the aggregate of all IMFs and a residual. All IMFs are utilized for reconstruction of SPI-3. Therefore, the MSE is computed between SPI-3 and the reconstructed SPI-3. It is found that the MSE of reconstructed SPI-3 has been found in Multan, Bahawalpur, Barkhan and Khanpur is 0.0024, 0.0052, 0.0034 and 0.0020, respectively.
The MK and SR trend tests are applied on SPI-3 of all stations to verify the reliability of residual trends. Initially, SPI-3 is divided into eight periods. So that the detailed information is obtained about the increasing and decreasing trends. Then the most commonly used 5% significance level is used to identify the trend. Hence, the statistical trend results of MK and SR is described in Table 2.
It is shown in Table 2 that the MK and SR indicate similar results of trend identification in all periods of Barkhan station. Whereas, MK and SR do not indicate similarity of trend identification in 2014 for Multan, in 1994 for Bahawalpur and 1998 period for Khanpur at the 5% level of significance. The positive sign is an indication of the increasing trend, while the negative sign is an indication of the decreasing trend. Moreover, the trend identification of different periods by MK and SR test is compared with the residuals of EEMD (Fig. 2) for all stations. It turns out that the SR test indicates the consistency with residuals of EEMD in mostly periods. For instance, the SR test verifies that the residuals of Multan indicate an increasing trend in 2014 while the trend does not verify by the MK test ( Fig. 2A). Likewise, the residuals of Khanpur indicates a decreasing trend in 1998 which is verified by the SR test, while not verified by the MK test (Fig. 2D). SR and MK provide similar trend identification of periods for all stations except the above-mentioned stations i.e., Multan and Khanpur.

Evaluation of weather instabilities
After this, there is no difficulty to apply an HT to obtain the IF and IA. We aim to describe the detailed information of wet and dry events in the time-frequency domain corresponding to each IMF. Therefore, HT is performed on decomposed IMFs of EEMD for all stations. Generally, HT uses the IMFs to present IF and IA corresponding to each IMF. Finally, the minimum, maximum, mean and standard deviation of IA and IF are related to each IMF is presented in Table 3. Table 3 showed that the highest variation (based on the standard deviation of IF) was observed in IMF1 of all stations. It is an indication of non-stationarity of the SPI-3 for all stations.
The MMHS provides an aggregated amount of IA in the time-frequency domain. So MMHS has been applied on the IA and IF of IMFs for all stations. The IA and IF are related to EEMD because this was achieved after HT is applied on the decomposed IMFs of SPI-3. The MMHS plots of four stations are shown in Fig. 3.
In Figs. 3A, 3B and 3D indicate the low IF of wet and dry events at 0.1 cycles per months (about log2(348) months) that is corresponding to mean IF of IMF-2 (see Table 3), while the IFs away from 0.1 are indicating comparatively low IA. Figure 3C shows the low IF of wet and dry events at 0.2 cycles per month (about < 1 months) that is corresponding to the mean IF of IMF-1 (see Table 3). Besides, the IFs away from 0.2 indicate comparatively low IA. It is described that Barkhan station underwent several seasonal instabilities over a short period of 5 months. Whereas, Multan, Bahawalpur and Khanpur stations underwent weather instabilities over a short period of 10 months.
The CWPS is performed on SPI-3 of all stations for the analysis of spectral features at the different time scales. The spectral features of SPI-3 can be seen in Fig. 4. Figure 4 shows the high energy wet and dry events with a periodicity of 32 months are identified during the years (1991-2017) for all stations. Figures 4A and 4B show a spectrum of similar features of wet and dry events with a period of 16 months in Multan and Bahawalpur, respectively. During the years (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010), the high energy wet and dry events are identified with a period of 16 months in Khanpur (Fig. 4D). While the high energy wet and dry events with the periodicity of 16 months are identified during the years (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011) in Barkhan station (Fig. 4C). Besides, some small wavelet spectrums are indicating high energy wet and dry events with 4 to 10 months of periodicity for all stations. But the high strength of wet and dry events is found in Multan, Bahawalpur, Barkhan and Khanpur are 1.49, 1.03, 1.42 and 1.18, respectively. All these high energy wetness and dryness has been identified at the 5% level of significance.

DISCUSSION
In recent decades, the analysis of spatial and temporal fluctuation of wet and dry events is gain more attention (Huang et al., 2014). Pakistan is an agricultural country where 70% of the people's business depends on agriculture (Ahmad et al., 2004). In this paper, SPI-3 was selected to describe the Spatio-temporal fluctuation of wet and dry events.
Mostly meteorological and hydrological series have non-linear and non-stationary characteristics (Di, Yang & Wang, 2014). The SR and MK non-parametric tests were Table 3 The HT is applied on eight IMFs of SPI-3 of all stations. The achieved IMFs are related to EEMD. Therefore, the minimum, maximum, mean, standard deviation of IA are represented by the Min IA, Max IA, MIA and SIA respectively. Whereas, the minimum, maximum, mean and standard deviation of IF are represented by the Min IF, Max IF, MIF and SIF, respectively.  applied on different periods of SPI-3 for the trend analysis. It is found that the SR test is indicating a trend for most of the periods as compared to the MK test (Table 2). In recent years, inter-decadal time scales have been carefully considered about the risk of wet and dry events trend (Li et al., 2015). Besides, EEMD was applied on SPI-3 to decompose into the eight IMFs and a residual term for the trend analysis. The achieved trend of wet and dry events has been compared with SR and MK identification results. It is observed that the SR test indicates consistent results with EEMD residuals for all stations ( Fig. 2 and Table 2). Moreover, the achieved IMFs have been used for reconstruction of SPI-3. As a result, the MSE of reconstructed SPI-3 is found to be less. It means that the decompose components contain the all characteristics of SPI-3. Therefore, the smaller MSE value is an indication of efficient decomposition process (Gaci, 2016). EEMD is a purely empirical procedure and gives complete information of the temporal attributes of agricultural drought. Any non-linear and non-stationary series can be adaptively decomposed into IMF components via EEMD method, without the need to a priori basis as are Fourier and wavelet-based methods. The EEMD basis functions are known as IMF (Gaci, 2016;Ayenu-Prah & Attoh-Okine, 2009).

Figure 4 CWPS is applied on SPI-3 for Multan (A), Bahawalpur (B), Barkhan (C) and Khanpur (D).
The colors of the wavelet power level reveal different levels of energy, and the black line indicates the 5% level of significance.
Full-size DOI: 10.7717/peerj.9729/ fig-4 The standard deviation of IF shows high variability of wet and dry events in IMF-1 (Table 3). This variability is indicating non-stationarity of SPI-3 for all stations. HHT gives detailed information about IA and IF in dissimilar time scales, which is beneficial for discovering temporal variability of the wet and dry events in multiple time scales, but the wavelet analysis cannot give details. Intra wave frequency is a traditional feature of the nonlinear process which represents IF changes in one oscillation cycle. Fourier analysis cannot show the intra wave frequency, because Fourier analysis is a linear structure that applies to the non-linear process. HHT is a very effective way of describing IF in which intra wave frequency is revealed (Huang, 2014). MMHS indicates that during the short period of 5 months, the Barkhan station has experienced severe weather instabilities, while Multan, Khanpur and Bahawalpur stations have suffered severe weather instability during the 10 months (Fig. 3). The different spectrum's of CWPS specifies 4 to 32 months of the periodicity of wetness and dryness during the years (1991-2017) for all stations. CWPS is indicating high energy weather instabilities in Multan, Barkhan, Bahawalpur and Khanpur respectively (Fig. 4). It is observed that the achieved short-term weather instabilities by MMHS provides consistent results with CWPS for all stations. By combining the results of both procedures i.e., MMHS and CWPS, it is concluded that the weather instabilities found in Multan station (1.48) is higher than Khanpur (1.18) station and in Bahawalpur (1.03) it shows least weather instabilities during the short period of 10 months. The intensity of weather instabilities for Barkhan station (1.42) during the short period of 5 months is higher as compared to other stations. It is very important to note that the MMHS provides information about weather instabilities during the short period of a month, but it does not provide information about the years in which the instabilities occurred. Whereas CWPS provides information about the monthly periodicity of weather instabilities for the time domain. The above-mentioned points illustrate the difference between MMHS and CWPS.

CONCLUSION
In this study, we have presented a thorough analysis of four meteorological stations (Multan, Bahawalpur, Barkhan and Khanpur) of Pakistan. We found that all stations are facing severe weather instabilities during the short period of 5 and 10 months using MMHS method.Moreover, the CWPS has shown the weather instabilities during 4 to 32 months of periodicity for all stations. Ultimately, the achieved short-term weather instabilities indicated by MMHS is consistent with CWPS. We also found that the performance of the SR test is better than the MK test for EEMD residuals trend analysis. In summary, these findings might be useful for water resource management and policymakers.