Characterization and Prediction of Water Stress Using Time Series and Artificial Intelligence Models

Abstract

In agroecosystems, drought is a critical climatic phenomenon that affects evapotranspiration and induces water stress in plants. The objective in this study was to characterize and forecast water stress in the Hyderabad region of India using artificial intelligence models. The monthly precipitation data for the period 1982–2021 was characterized by the standardized precipitation index (SPI) and modeled using the classical autoregressive integrated moving average (ARIMA) model and artificial intelligence (AI), i.e., artificial neural network (ANN) and support vector regression (SVR) model. The results show that on the short-term SPI3 time scale the studied region experienced extreme water deficit in 1983, 1992, 1993, 2007, 2015, and 2018, while on the mid-term SPI6 time scale, 1983, 1991, 2011, and 2016 were extremely dry. In addition, the prediction of drought at both SPI3 and SPI6 time scales by AI models outperformed the classical ARIMA models in both, training and validation data sets. Among applied models, the SVR model performed better than other models in modeling and predicting drought (confirmed by root mean square error—RMSE), while the Diebold–Mariano test confirmed that SVR output was significantly superior. A reduction in the prediction error of SVR by 48% and 32% (vs. ARIMA), and by 21% and 26% (vs. ANN) was observed in the test data sets for both SPI3 and SPI6 time scales. These results may be due to the ability of the SVR model to account for the nonlinear and complex patterns in the input data sets against the classical linear ARIMA model. These results may contribute to more sustainable and efficient management of water resources/stress in cropping systems.

1. Introduction

Meteorological and hydrological droughts commonly cause agricultural drought, i.e., water stress in agroecosystems. The three longest droughts globally occurred between July 1928 and May 1942 (Dust Bowl drought of the 1930s), July 1949 and September 1957 (1950s drought), and between June 1998 and December 2014 (early 21st-century drought). However, the frequency of droughts has increased in recent decades, reducing the productivity of the agri-food sector [1,2,3]. For example, in 2017 at least 3 percent of the world’s land area was affected by extreme drought [4]. Drought is a complex natural hazard that affects global food security through water stress, making the agri-food sector one of the most affected by drought. Namely, water stress negatively affects many physiological functions of plants, including photosynthesis, transpiration, and nutrient uptake, and reduces vegetative growth and crop yield. Water stress is projected to worsen under global climate change and thus further threaten global food security [5,6,7]. Climate change is affecting the production of most agricultural commodities and it is projected to decline further [8].

Drought detection, characterization, and prediction via different approaches are critical for decision support systems and more successful drought, i.e., water management, and thus food security. Moreover, drought modeling and prediction can be a very useful tool for more sustainable water planning to avoid future uncertainties that may arise from water scarcity and stress [9], not only in agroecosystems but also in other sectors.

As a result, much attention has been dedicated to drought characterization and modeling, and continuous progress has been recorded in this field. However, a major research challenge still remains in the development of an appropriate approach to predicting the onset and termination of droughts, given the critical stage in mitigating the effects of drought is the inability to accurately predict drought conditions for the upcoming short-/medium-term period [9]. In addition, despite significant improvements in drought modeling in recent decades, not many studies have used advanced approaches, such as artificial intelligence (AI) modeling, to characterize and predict droughts. For example, many studies have used the standardized precipitation index (SPI) model developed by [10] to characterize drought. The author of [11] modeled climatological precipitation using the gamma distribution, ref. [12] studied hydrological drought using the SPI for medium and long time periods, while ref. [13] used the SPI to calculate long-term spatial and temporal trends in drought. Precipitation patterns were analyzed by [14] who studied drought severity using SPI, while ref. [15] characterized drought in Italy using SPI at different time scales. Furthermore, ref. [16] used SPI to study meteorological droughts for different climate models at different time scales. SPI for drought characterization using multi-satellite precipitation data from the Tropical Precipitation Measurement Mission were computed by [17]. SPI was developed by [18] for characterizing short and intermediate droughts, while ref. [19] characterized drought in the Lower Mekong Basin in Southeast Asia using SPI at different time scales. Drought conditions for maize in different growing seasons using SPI values in the Adana, Mersin, and Osmaniye regions of Çukurova were determined by [20], whereas ref. [21] applied SPI to calculate drought and heatwaves in China.

Box–Jenkins autoregressive integrated moving average (ARIMA) is the most commonly used time series model for predicting linear time series data [22]. There are many cases where ARIMA is used for modeling and forecasting drought and other agricultural phenomena [23,24,25]. However, the classical time series ARIMA model is not able to model the nonlinear, heterogeneous, and complex patterns present in the time series and does not provide reliable forecasts. Many nonparametric machine learning models have been developed and are available in the literature to solve the problem of modeling nonlinear components in time series data. However, the only approach to model and predict such time series is still the use of AI techniques [26,27,28,29]. Artificial neural networks (ANN) and support vector regression (SVR) is the most commonly used techniques for modeling and prediction of drought and other agricultural variables [23,30,31,32,33,34].

The Hyderabad region in the Indian state of Telangana is characterized by intensive agroecosystems affected by frequent drought and water stress [35]. According to recent reports, 7/10 districts were affected by drought (Association for India’s development, 2016 survey report), threatening the water supply for food production and other industrial and municipal sectors. Therefore, the development of drought characterization and its prediction would greatly improve drought management at the multi-sectoral scale. For the selected Hyderabad region, this is the first attempt to model and predict drought, i.e., water stress, using AI techniques such as ANN and SVR. Moreover, this work was the first attempt to significantly compare the modeling performance of the studied techniques using the test of Diebold and Mariano [36]. The methodological framework begins with the characterization of drought by SPI. The AI-based ANN and SVR models presented here were developed to model and predict drought, not only in the Hyderabad region, but also in similar water-scarce and drought-affected agroecosystems. Thus, the results could serve as a valuable tool and extension in drought, i.e., water management, and food security.

2. Materials and Methods

2.1. Description of the Data and Study Area

The study pertains to Hyderabad, the capital city of Telangana state of India. Telangana has an area of 112,077 sq. km and has a population of 35,003,674. The state receives an average annual rainfall of about 1000 mm. The Hyderabad (17.3850° N, 78.4867° E) region of the Hyderabad region of the Telangana state is selected as the study area, as the average rainfall received by the region is usually less than the state average. The annual rainfall of the Hyderabad region was observed to be around 750 mm, the region receives most of its rainfall from the southwest monsoon. The monthly rainfall data (millimeter (mm)) in Hyderabad from January 1988 to May 2021, was obtained from the Directorate of Economics and Statistics, Government of Telangana. Further, the rainfall data is used to calculate the SPI for three and six-month time scales.

2.2. Standardized Precipitation Index (SPI)

The rainfall data were used to calculate the SPI indices for time scales of 3 (SPI3) and 6 (SPI6) months. The SPI was calculated based on the probability of precipitation for a specific time scale. The gamma distribution is defined by its frequency or probability density function.

$$f\left(x,\alpha ,\beta \right)=\frac{1}{{\beta }^{\alpha }\mathsf{\Gamma }\left(\alpha \right)}{x}^{\alpha -1}{e}^{-\left(\frac{x}{\beta }\right)}$$

(1)

where, α and β are the shape and scale parameters respectively; x is the rainfall amount (millimeter), and Γ(α) is the gamma function, α > 0,

$$\mathsf{\Gamma }\left(\alpha \right)={{\displaystyle \int }}_0^{\infty }{x}^{\alpha -1}{\mathrm{e}}^{-x}$$

(2)

where, Γ(α) is the gamma function.

Computation of the SPI was mainly involved in fitting a gamma probability density function to a given frequency distribution of precipitation. The maximum likelihood solutions used for optimal estimation of α and $\beta$ are given as follows [9]:

$$\alpha =\frac{1}{4A}\left(1+\sqrt{1+\frac{4A}{3}}\right)·\alpha >0,$$

(3)

$$\beta =\frac{\overline{x}}{\alpha },\beta >0$$

(4)

$$A=\mathit{ln}\left(\overline{x}\right)-\frac{{{\displaystyle \sum }}_{i=1}^n\mathit{ln}\left(x\right)}{n},$$

(5)

where, n = number of precipitation observations.

Integrating the probability density function with respect to x and attaching α and β parameters yields the cumulative probability distribution function F(x):

$$F\left(x,\alpha ,\beta \right)={{\displaystyle \int }}_0^{x}f\left(x,\alpha ,\beta \right)=\frac{1}{{\beta }^{\alpha }\mathsf{\Gamma }\left(\alpha \right)}{{\displaystyle \int }}_0^x{x}^{\alpha -1}{e}^{-\left(\frac{x}{\beta }\right)}dx$$

(6)

Equation (6) can be expressed as follows:

$$F\left(x,\alpha ,\beta \right)=\frac{1}{\mathsf{\Gamma }\left(\alpha \right)}{{\displaystyle \int }}_0^x{t}^{\alpha -1}{e}^{-\left(t\right)}dx$$

(7)

where, t = x/β.

Since the gamma function is undefined for x = 0 and the rainfall time series data may contain zero rainfalls, the cumulative probability of zero and non-zero rainfalls, H(x) can be calculated as:

$$H\left(x\right)=q+\left(1-q\right)F\left(x,\alpha ,\beta \right)$$

(8)

where, q is the probability of zero rainfall. If m is the number of zeros present in a rainfall time series, then q is estimated by m/n. The cumulative probability was then transformed into a standard normal random variable Z, then its mean is 0 and variance is 1, which is the value of the SPI.

Ref. [37] derived the following equations to transform the cumulative probability distribution into a standardized normal distribution to calculate the SPI as follows:

$$\mathrm{SPI}=-\left(K-\frac{c_0+c_{1}K+c_2{K}^2}{1+d_1+d_2{k}^2+d_3{K}^3}\right)\mathrm{when}K=\sqrt{ln\left(\frac{1}{{(H\left(x\right))}^2}\right)}$$

(9)

for 0 < H(x) ≤ 0.5

$$\mathrm{SPI}=+\left(K-\frac{c_0+c_{1}K+c_2{K}^2}{1+d_1+d_2{k}^2+d_3{K}^3}\right)\mathrm{when}K=\sqrt{ln\left(\frac{1}{1-{(H\left(x\right))}^2}\right)}$$

(10)

for 0.5 < H(x) ≤ 1

where, C₀ = 2.515517, C₁ = 0.802853, C₂ = 0.010328, d₁ = 1.432788, d₂ =0.189269, and d₃ = 0.001308. The SPI for time scales 3 (SPI3) and 6 (SPI6) represent short-term and mid-term meteorological droughts, respectively [37]. The SPI is categorized in to different groups as mention in

Table 1. Category of Standardized precipitation index.

SPI Range	Category
+2 to more	Extremely wet
1.5 to1.99	Very wet
1.0 to 1.49	Moderately wet
−0.99 to 0.99	Near normal
−1.0 to −1.49	Moderately dry
−1.5 to −1.99	Severely dry
−2 to less	Extremely dry

SPI: standardized precipitation index.

below:

After the development of SPI3 and SPI6, the indexes from January 1982 to May 2020 were considered as the training dataset (data for model building), while the indexes from June 2020 to May 2021 were considered as the testing dataset (data for model validation) in this study.

2.3. ARIMA Model

The ARIMA model was characterized by three terms: p, d, and q, where, ‘p’ is the order of the ‘auto regressive’ (AR) term, referring to the number of lags of Y_t time series to be used as predictors. ‘q’ is the order of the ‘moving average’ (MA) term, assuming the number of lagged forecast errors that should go into the ARIMA model. ‘d’ is the number of differences required to make the time series stationary; the most common approach of differencing is, subtract the previous value from the current value.

$$\varnothing \left(B\right){\left(1-B\right)}^d{Y}_t=\theta \left(B\right)e_t$$

(11)

$$Y_t={\varnothing }_1{Y}_{t-1}+{\varnothing }_2{Y}_{t-2}+\cdots +{\varnothing }_p{Y}_{t-p}+{\epsilon }_t-{\theta }_1{e}_{t-1}-{\theta }_2{e}_{t-2}-\cdots -{\theta }_q{e}_{t-q}$$

(12)

where, $Y_t$ is the time series, ${\varnothing }_i$ and θj are model parameters for AR and MA terms, respectively, $e_t$ is the random error, p is a number of autoregressive terms, q is number of lagged forecast errors, and B is the backshift operator such that $B{Y}_t=Y_{t-1}$. The ARIMA model building consists of four stages, viz. identification, estimation, diagnostic checking, and forecasting [22].

2.4. Artificial Neural Network (ANN)

The ANN is one of the most widely used artificial intelligence techniques in the last two decades, for both classification and regression tasks. As time-series falls under regression problem, the lagged autoregressive observations were used as input to the network with an implicit functional representation of time. The mathematical expression for the output Y_t of a single-layer or multi-layer feed-forward autoregressive network is expressed as follows:

$$Y_t={\alpha }_0+{\displaystyle \sum }_{j=1}^q{\alpha }_{j}g\left({\beta }_{0j}+{\displaystyle \sum }_{i=1}^p{\beta }_{ij}{Y}_{t-i}\right)+e_t$$

(13)

where, ${\alpha }_j\left(j=0,1,2,\dots ,q\right)$ and ${\beta }_{ij}\left(i=0,1,2,\dots ,p,j=0,1,2,\dots ,q\right)$ are the model parameters, p is the number of input nodes, q is the number of hidden nodes, g is the activation function, and $e_t$ is the random error [38]. The logistic function is the most commonly used activation function in time series modeling.

$$g=\frac{1}{1+exp\left(-x\right)}$$

(14)

The ANN model building was generally divided into the model building part (training data set) and validation part (testing set). In model building, the algorithm minimizes the error function between actual and predicted values. The error function was expressed as follows:

$$E=\frac{1}{N}{{\displaystyle \sum }}_{t=1}^N{(e_t)}^2=\frac{1}{N}{{\displaystyle \sum }}_{t=1}^N{\left\{X_t-(w_0+({{\displaystyle \sum }}_{J=1}^Q{w}_{J}g\left(w_{oj}+{{\displaystyle \sum }}_{i=1}^P{w}_{ij}{X}_{t-i}))\right)\right\}}^2$$

(15)

where, N is the total number of error terms and $w_{ij}$ are the weight parameters of the neural network which are changed by a number of changes in $\Delta w_{ij}=-\eta \frac{\partial E}{\partial w_{ij}}$ where $\eta$ is the learning rate which is a user-defined hyper-parameter [39,40]. The schematic representation of the autoregressive feed-forward neural network is depicted in Figure 1.

2.5. Support Vector Regression (SVR)

The support vector regression (SVR) is a supervised machine learning technique originally developed for binary classification and later extended to nonlinear regression estimation problems by introducing an ε-insensitive loss function [41]. The basic principle involved in NLSVR (nonlinear SVR) is to transform the original input time series into a high dimensional feature space and then build the regression model in a new feature space.

Consider a vector of data set $Z={\left\{x_{i }{y}_i\right\}}_{i=1}^N$ where $x_i\in R^n$ is an input vector, y_i is the scalar output, and N is the size of data set. The general equation SVR can be written as follows [41]:

$$f\left(x\right)=W^T\varphi \left(x\right)+b$$

(16)

where, $\varphi \left(x\right)$ is function of x, W is weight vector, b is bias term, and superscript T denotes the transpose. The coefficients W and b were estimated by minimizing the following regularized risk function:

$$R\left(\theta \right)=\frac{1}{2}\Vert w{\Vert }^2+C\left[\frac{1}{N}{{\displaystyle \sum }}_{i=1}^N{L}_{\epsilon }\left(y_i,f\left(x_i\right)\right)\right]$$

(17)

The first term $\frac{1}{2}\Vert w{\Vert }^2$ is called ‘regularized term’, which measures the flatness of the function. Minimization of the regularized terms will make a function as flat as possible. The second term $\frac{1}{N}{\displaystyle \sum }_{i=1}^N{L}_{\epsilon }\left(y_i,f\left(x_i\right)\right)$ is called as ‘empirical error’ or loss function $L_{\epsilon }$ which was estimated by Vapnik ε-insensitive loss function as follows:

$$L_{\epsilon }\left(y_i,f\left(x_i\right)\right)=f\left(x\right)=\left\{\begin{array}{c}\left|y_i,f\left(x_i\right)-\epsilon \right|; \left|y_i-f\left(x_i\right)\right|\ge \epsilon ,\\ 0; \left|y_i-f\left(x_i\right)\right|<\epsilon ,\end{array}\right.$$

(18)

where, y_i is the actual value, $f\left(x_i\right)$ is the estimated value, and $\epsilon$ is the margin of tolerance where no penalty is given to errors.

The kernel function in SVR plays a critical role in mapping the input to output, where the information is processed using a suitable kernel function. The most commonly used kernel functions radial basis function (RBF) K is given as follows:

$$K\left(Y_1,Y_2\right)=\mathit{exp}\left(-\frac{{\left(\parallel x_1-x_2\parallel \right)}^2}{2{\sigma }^2}\right)$$

(19)

where, ${\sigma }^2$ is the variance, and $\Vert x_1-x_2\Vert$ is the Euclidean distance between two sample points $x_1$ and $x_2$. The performance of SVR depends upon the optimization of two hyper-parameters; regularization parameter C also called the cost function, which balances the complexity and approximation accuracy of the model and kernel bandwidth parameter $\gamma =\frac{1}{2{\sigma }^2}$, which represents the variance of the kernel function.

Finally, a general flow chart of the methodology followed in this study is delineated in Figure 2. The flow chart depicts the development of drought indices SPI3 and SPI6 using rainfall data, modeling, and forecasting of drought using different models and performance of models in both training and testing data sets in terms of RMSE and Diebold–Mariano (DM) test.

3. Results

3.1. Calculated SPI

The calculated water stress indexes for SPI3 and SPI6 timescales (Table S1) are plotted in Figure 3a,b, respectively. The drought events were developed from 1982 to 2021 for both SPI3 and SPI6 timescales (

Table 2. Calculated drought months by short-term SPI3.

Year	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec
1982										−1.21
1983	−1.56	−2.04	−2.25	−3.15	−2.13				−1.30	−1.48
1984						−1.56	−1.23	−1.89
1985												−1.15
1988												−1.30
1989	−1.19	−1.43		−1.22							−1.08	−1.34
1990	−1.34	−1.64	−1.78	−1.17
1991			−1.51
1992			−2.06			−1.11	−1.49
1993			−2.04			−1.12	−1.71
1994										−1.21
1997								−1.30	−1.48
1998					−1.60	−1.00	−1.92
1999											−1.14
2002										−1.26
2003			−1.37			−1.48	−1.67
2004									−1.15		−1.13
2007			−2.08
2010				−1.32	−1.67
2015							−2.25
2018	−2.08	−2.49							−1.27	−1.97	−1.11	−1.80
2019	−1.12				−1.82	−1.66

Jan: January, Feb: February, Mar: March, Apr: April, Jun: June, Jul: July, Aug: August, Sep: September, Nov: November, Dec: December.

and

Table 3. Calculated drought months by mid-term SPI6.

Year	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec
1983	−1.18	−1.08			−2.66
1985						−1.37	−1.33	−1.36
1986										−1.86	−1.53
1989	−1.27	−1.35	−1.19
1990	−1.11	−1.18
1991					−2.42
1996							−1.37	−1.28	−1.07
2000									−1.08		−1.02	−1.04
2002								−1.21	−1.86	−1.56		−1.98
2003	−1.72	−1.24
2007						−1.53
2011		−2.20	−1.85
2016			−2.66
2018											−1.95	−1.62
2019	−1.66

Jan: January, Feb: February, Mar: March, Apr: April, Jun: June, Jul: July, Aug: August, Sep: September, Nov: November, Dec: December.

). For instance, ten months in 1983 (−2.24 SPI3_Feb, −2.25 SPI3_Mar, −3.15 SPI3_Apr, −2.13 SPI3_May), 1992 (−2.06 SPI3_Mar), 1993 (−2.04 SPI3_ Feb), 2007 (−2.08 SPI3_ Mar), 2015 (−2.25 SPI3_Jul), and 2018 (−2.08 SPI3_Jan, −2.49 SPI3_Feb) experienced extreme water stress on SPI3 time scale (Table 2, Figure 3a). In addition, on the same scale during 1982, 1983, 1984, 1985, 1988, 1989, 1990, 1992, 1993, 1994, 1997, 1998, 1999, 2002, 2003, 2004, 2010, 2018, and 2019 were recorded as moderately dry, while in 1983, 1984, 1990, 1991, 1998, 2003, 2010, 2018, and 2019, severely dry conditions were observed (Table 2).

On the other hand, SPI6 indicates medium-term drought conditions, and on this time scale there were four months with extreme water stress in 1983 (−2.66 SPI6_May), 1991 (−2.42 SPI6_May), 2011 (−2,20 SPI6_June), and 2016 (−2.38 SPI6_Mar), moderately dry conditions were confirmed in 1983, 1985, 1989, 1990, 1996, 2000, 2002, and 2003, while severely dry conditions were recorded in 1986, 2002, 2003, 2007, 2011, 2018, and 2019 (Table 3, Figure 3b).

3.2. ARIMA Model

The time series plots of SPI at 3 and 6 timescales are depicted in Figure 4a,b. The time series modeling ‘Box–Pierce test’ was performed to check the autocorrelation of the series using the following formula [22]:

$$Q=n\left(n+2\right){{\displaystyle \sum }}_{k=1}^h\frac{\widehat{\rho }{2}_k}{n-k}$$

(20)

where n is the sample size, ${\widehat{\rho }}_k$ is the sample autocorrelation at lag k, and h is the number of lags being tested; and Q follows ${\chi }^2$, a distribution with h degrees of freedom. Both SPI3 and SPI6 are autocorrelated as the Box–Pierce test was significant (p < 0.0001) (

Table 4. Parameter estimation of ARIMA by maximum likelihood method for SPI3 and SPI6 time scales of Hyderabad region.

Time Scales	Model	Parameters	Estimation	S. E	Z Value	Prob.	Model Fitting		Box–Pierce Non-Correlation Test
SPI3	ARIMA (1,0,2)	AR1	−0.018	0.130	−0.143	p < 0.0001	Log likelihood	−544.14	Original	Residuals
		MA1	0.674	0.121	5.559	p < 0.0001	BIC	1112.79	Chi-square = 148.39 (p < 0.001)	Chi-square = 0.005 (p = 0.938)
		MA2	0.393	0.066	5.956	p < 0.0001	AIC	1096.27
SPI6	ARIMA (1,0,0)	AR1	0.717	0.032	22.03	p < 0.0001	Log likelihood	−479.53	Chi-square = 233.19 (p < 0.001)	Chi-square = 0.078 (p = 0.779)
							AIC	963.07
							BIC	971.31

SPI: standardized precipitation index, ARIMA: auto-regressive integrated moving average, AR: autoregressive, MA: moving average, S.E: standard error, prob. (p): probability, BIC: Bayesian information criterion, AIC: Akaike information criterion.

). The main step in the ARIMA model is the identification of the different combinations of parameters, such as autoregressive terms (p), differencing terms (d), and moving average terms (q). A combination of the order of the parameters which has maximum log-likelihood and lowest values of Akaike information criteria (AIC) and Bayesian information criteria (BIC) is considered the best model. Along with log-likelihood, AIC and BIC criteria’s ACF and PACF plots were also used (Figures S1–S4) for model order identification. The results are presented in Table 4. For SPI3 and SPI6 the best fitted ARIMA model was (1,0,2) and (1,0,0) with maximum log-likelihood of −544.14 and −479.53, minimum AIC and BIC values are 1096.27 and 1112.79 for SPI3 and 963.07 and 971.31 for SPI6. After ARIMA model fitting, the diagnostic checking of residuals was carried out using the Box–Pierce non-correlation test, and residuals under consideration are non-correlated in nature as the probability of significance are p = 0.938 and p = 0.779 for SPI3 and SPI6, respectively.

For building the AI model, the nonlinearity of the original data set was tested using the BDS test [42]. The results of the BDS test are shown in

Table 5. BDS test for non-linearity.

Sample	Dimension	SPI3		SPI6
		Statistics	Probability	Statistics	Probability
eps(1)	m = 2	15.63	p < 0.0001	15.68	p < 0.0001
	m = 3	15.37	p < 0.0001	18.09	p < 0.0001
eps(2)	m = 2	14.23	p < 0.0001	12.72	p < 0.0001
	m = 3	13.11	p < 0.0001	12.09	p < 0.0001
eps(3)	m = 2	13.90	p < 0.0001	12.51	p < 0.0001
	m = 3	12.34	p < 0.0001	11.18	p < 0.0001
eps(4)	m = 2	15.40	p < 0.0001	13.92	p < 0.0001
	m = 3	13.76	p < 0.0001	12.43	p < 0.0001

eps: epsilon values (standard deviation of input variable as sample), m: embedded dimension, p: probability.

and confirm that the data under consideration is nonlinear in nature as the probability of significance of residuals is p < 0.0001.

3.3. ANN Model

ANN model parameters are given in

Table 6. Parameter specification of ANN model.

Parameters	SPI3	SPI6
Input lag	5	4
Output variable/dependent	1	1
Hidden nodes	1	1
Hidden layers	2	4
Model	(5:2S:1L)	(4:4S:1L)
Total number of parameters	15	25
Network type	Feed forward	Feed forward
Activation function (I:H)	Sigmoidal	Sigmoidal
Activation function (H:O)	Identity	Identity
Box–pierce non-correlation test for residuals	0.018	0.124
p-value	0.873	0.724

I: input node, H: hidden node, S: sigmoidal activation function used in input to the hidden layer, L: linear unit activation function used in hidden to the output layer, p-value: probability value.

for both SPI3 and SPI6 time scales. The drought indices of time series were subjected to ANN model with 5 and 4 tapped time delays and 2 and 4 optimum nodes for SPI3 and SPI6, respectively. Sigmoidal activation function in input to the hidden layer and linear identity function in hidden layers to output layer was used with feed-forward network architecture. The total number of parameters or weights obtained was 15 and 25 for SPI3 and SPI6, respectively. Diagnosis checking of residuals by Box–Pierce non-correlation test, which indicates the non-auto correlated or random nature of the residuals (p = 0.83 and p = 0.72) for SPI3 and SPI6 was done.

3.4. SVR Model

For the SVR model, parameters were considered based on the least RMSE values. The parameter specifications are listed in

Table 7. Parameter specification of SVM model.

Model	SPI3	SPI6
Parameters	Values	Values
Kernel function	RBF	RBF
No. of S. V’s	451	337
Cost	7.9	7.7
Gamma	0.16	0.16
Epsilon	0.01	0.1
Cross validation error	0.043	0.029
Box–Pierce non-correlation test for residuals	0.19	0.31
p-value	0.27	0.57

RBF: radial basis function, p-value: probability value.

for both SPI3 and SPI6 timescales, where the radial basis function (RBF) was used as a kernel function, and the number of support vectors obtained are 451 and 337, respectively. The cost function, gamma, and epsilon are 7.9, 0.16, and 0.01 for SPI3, and 7.7, 0.16, and 0.1 for SPI6, respectively. Box–Pierce tests done for diagnostic checking of residuals are random and white noise in nature as p = 0.27 and 0.53 for both SPI3 and SPI6, respectively.

4. Discussion

Drought characterization by using SPI3 and SPI6 timescales indicated different categories of water stresses, i.e., moderately dry, severely dry, and extremely dry. Ten months showed extreme water stress in SPI3 and four months showed water stress on the SPI6 time scale, which indicates high water scarcity during the reported period (

Table 8. Classification of drought categories in SPI3 and SPI6.

Category	SPI3	SPI6
Extremely wet	-	-
Very wet	-	-
Moderately wet	-	-
Near normal	-	-
Moderately dry	Oct-1982, Sep-1983, Oct-1983, Jul-1984, Dec-1985, Dec-1988, Jan-1989, Feb-1989, Apr-1989, Nov-1989, Dec-1989, Jan-1990, Apr-1990, Jun-1992, Jul-1992, Jun-1993, Oct-1994, Aug-1997, Sep-1997, Jun-1998, Nov-1999, Oct-2002, Mar-2003, Jun-2003, Sep-2004, Nov-2004, Apr-2010, Sep-2018, Nov-2018, Jan-2019	Jan-1983, Feb-1983, Jun-1985, Jul-1985, Aug-1985, Jan-1989, Feb-1989, Mar-1989, Jan-1990, Feb-1990, Jul-1996, Aug-1996, Sep-1996, Fep-2000, Nov-2000, Dec-2000, Aug-2002, Feb-2003
Severely dry	Jan-1983, Jun-1984, Aug-1984, Feb-1990, Mar-1990, Mar-1991, May-1998, Jul-1998, Jul-2003, May-2010, Oct-2018, Dec-2018, May-2019, Jun-2019	Oct-1986, Nov-1986, Sep-2002, Oct-2002, Dec-2002, Jan-2003, Jun-2007, Mar-2011, Nov-2018, Dec-2018, Jan-2019
Extremely dry	Feb-1983, Mar-1983, May-1983, Mar-1992, Mar-1993, Jul-1993, Mar-2007, Jul-2015, Jan-2018, Feb-2018	May-1983, May-1991, Feb-2011, Mar-2016

). Similar results for different classes of droughts and different time scales were also reported in [43,44,45,46,47,48,49,50].

The modelling and forecasting of SPI3 and SPI6 were compared in terms of MSE, RMSE for both training and testing sets for SPI3 and SPI6, respectively. AI models outperformed the classical model, SVR model was found to be superior over ARIMA and ANN models in training and testing sets for both SPEI3 and SPEI6. Performance hierarchy of these models is as follows: SVR > ANN > ARIMA in both training and validation sets. The detailed comparison of modelling and forecasting performance of different models in training and testing sets are depicted in

Table 9. Model comparison for training data for SPI3 and SPI6.

Model	Parameter	ARIMA	ANN	SVR
SPI3	MSE	0.438	0.418	0.378
	RMSE	0.496	0.459	0.413
SPI6	MSE	0.328	0.275	0.237
	RMSE	0.363	0.317	0.289

MSE: mean square error, RMSE: root mean square error.

and

Table 10. Model comparison for testing data for SPI3 and SPI6.

Period	Actual		Forecasted
			ARIMA		ANN		SVR
	SPI3	SPI6	SPI3	SPI6	SPI3	SPI6	SPI3	SPI6
Jun-20	0.917	0.851	0.241	0.083	0.132	0.254	0.197	0.064
Jul-20	1.477	1.418	0.194	0.11	0.07	0.182	0.713	0.669
Aug-20	0.133	0.12	−0.004	0.082	−0.089	0.131	0.554	0.83
Sep-20	0.226	0.262	−0.689	0.085	−0.068	0.094	−0.085	0.352
Oct-20	0.279	0.400	−0.812	0.059	−0.083	0.067	−0.186	−0.501
Nov-20	1.734	0.369	−0.987	0.433	0.357	0.048	0.403	0.338
Dec-20	1.142	0.322	−0.988	0.025	−0.072	0.035	0.87	0.294
Jan-21	−0.342	0.23	−0.567	0.096	−0.069	0.025	−0.549	0.25
Feb-21	−0.636	1.646	−0.288	−0.028	−0.074	0.018	−0.804	0.839
Mar-21	−0.853	0.84	−0.755	−0.013	−0.07	0.013	−0.025	1.416
Apr-21	−0.138	−0.526	1.988	−0.02	−0.072	0.009	−0.994	0.399
May-21	1.709	1.491	1.832	−0.025	−0.07	0.007	−0.111	0.227
MSE			1.317	0.72	0.867	0.672	0.686	0.495
RMSE			1.734	0.848	0.931	0.82	0.828	0.703

MSE: mean square error, RMSE: root mean square error.

. Forecasted values and model comparison criteria of testing data sets are given in Table 10. Further, actual vs. fitted data of both training and testing sets are given in Table S2 and are also plotted in Figure 5a,b for drought in both SPI3 and SPI6 time scales, respectively. Similar results were obtained in [51,52,53,54], where the SVR model outperformed ARIMA and ANN models.

The ANN for water stress modeling was developed with the ‘Levenberg–Marquardt backpropagation algorithm’ in a feed-forward network based on repetitive iterations. The user-defined hyper-parameters for ANN model, such as learning rate and momentum, were fixed based on repetitive experimentation as 0.04 and 0.001, respectively. To obtain the stable model, the network was repeated 30 times with 1250 iterations. The candidate input lags and hidden nodes were determined based on the lowest training errors. Tenfold cross-validation was carried out for each model and the lowest cross-validation error obtained is reported in Table 6 and Table 7. For SPI3 and SPI6 time scales we tuned the model with different permutations and combinations of hyper-parameters and chose the candidate parameters based on the lowest training MSE and RMSE. A similar procedure was followed in [55] for building ANN and SVR models.

The comparison criteria, viz., MSE and RMSE exhibit only the observed differences between the predicted values of the models. The DM test statistics were used to determine the statistically significant differences among the models used. The ARIMA (M1) and ANN (M2) models are significantly different with respect to SVR (M3) model for both SPI3 and SPI6. Inter combinational significances are clearly mentioned in

Table 11. Diebold–Mariano (DM) test for comparing performance of different time series models for training and testing data set for both SPI3 and SPI6.

Model	Data Type	M1, M2	M1, M3	M2, M3
SPI3	Training set	5.62 (p < 0.0001)	5.05 (p < 0.0001)	−0.38 (p < 0.0001)
	Testing set	3.67 (p < 0.0001)	2.57 (p < 0.0001)	−0.12 (p < 0.0001)
SPI6	Training set	7.77 (p < 0.0001)	6.24 (p < 0.0001)	4.07 (p < 0.0001)
	Testing set	0.75 (p < 0.0001)	0.62 (p < 0.0001)	0.52 (p < 0.0001)

M1: ARIMA, M2: ANN, M3: SVR, p: probability value.

for both SPI3 and SPI6. The performance of SVR model is significantly superior to the other two models in the input and testing data sets for both SPI3 and SPI6. These results may be due to the ability of the SVR model to account for the nonlinear and complex patterns in the input data sets, which can be a comprehensive tool for managing and alleviating water stress in plants. Similarly, the DM test was applied in [27,56,57] for significant comparison of modeling performance of the developed model in both training and validation data sets.

5. Conclusions

The present study was carried out to understand the short and medium-term water stress in agroecosystems of the Hyderabad region using SPI3 and SPI6 timescales. The results showed that during the period 1982–2021, 13 months were identified as severe drought, 10 months as extreme drought, and 29 months as moderate drought on the short-term scale (SPI3). Similarly, the medium-term drought index (SPI6) identified 4 months as extreme drought, 11 months as severe drought, and 18 months as moderate water stress during the period 1982–2021. In addition, classical models and AI models were used to predict water stress, confirming that among the AI models, SVR performs better vs. ANN, and the classical ARIMA model in terms of MSE and RMSE. The results can be a very useful tool for management and planning at the regional/national level of water (stress) as one of the most important and critical resources in cropping systems. Due to the lack of availability of continuous historical data on other weather parameters for the study region, only rainfall data was used for estimating the drought index in this study. However, future studies need to be expanded with either weather (e.g., air temperature, solar radiation), hydrologic (e.g., water table), agronomic (e.g., yield reduction), or in combination, parameters to allow detailed characterization of drought using various machine learning-based models.