Agricultural drought periods analysis by using nonhomogeneous poisson models and regionalization of appropriate model parameters

Abstract

Precipitation has a dominant role in agriculture, and a regular rain pattern is usually vital to agriculture; excessive or inadequate rainfall can be harmful. In this paper, an agricultural drought index is utilized to study the agricultural drought periods and analyze them with their intensities at various locations. Some nonhomogeneous Poisson models are also used to calculate the probability of agricultural droughts (number of times occurred) in a time interval of interest. It is to be assumed that the number of agricultural drought event occurrences is a Nonhomogeneous Poisson Process (NHPP) has a rate function, which depends on some parameters that must be estimated. Two cases of these functions are considered: the Weibull and linear intensity function. The Bayesian approach with Gibbs sampling under the Markov Chain Monte Carlo (MCMC) algorithm is used to estimate the parameters of these functions. The most appropriate fitted model is selected by using Deviance Information Criteria (DIC) and use that appropriately fitted model to calculate the accumulated events of agricultural drought in a time interval of interest at each location. Ordinary Kriging (OK) is used to regionalize the parameters and present its spatial behavior. The results based on the DIC indicate that the Power Law Process (PLP) performs better than the linear intensity function, NHPP model. The interpolated parameter values for the appropriate model, their patterns, and fluctuations for the study area are efficiently presented using contour maps. It is a novel and straightforward approach to assess the selected model parameter values used to predict the accumulated drought events at un-sampled locations. The proposed framework might also help to analyze other spatial variables of interest and can be used for climate-change study, ecosystem modeling, etc. The findings can also help to make decisions for sustainable environmental management in Pakistan.

Keywords:

• standard precipitation index
• agricultural drought
• nonhomogeneous Poisson process
• power law process
• linear intensity function
• regionalization
• variogram
• ordinary kriging

1. Introduction

In recent years, it has been observed that the climate changes on the planet, causing the irregularity in the behavior of rainfall or precipitation. Often with little rain in some areas and excessive rain in other regions, which adversely affects agriculture, forests, water reservoirs, and biodiversity. Over the past several decades, extensive and frequent drought events had disastrous effects on the environment, economy, and society (Mishra and Singh, 2010). In Pakistan, this abnormal behavior has been accounted for; some parts of Pakistan face water supply shrinkage. Therefore, rarely observed a shortage of water in water reservoirs that supply water to the agricultural areas.

Low precipitation and unpredictable weather patterns are the leading causes of drought in an area. According to Bayissa et al. (2015), drought is a natural hazard indicated by a remarkable reduction in water availability during an extended period over a large area. Drought frequently occurs because it is a part of the climate, but timely action (i.e., before, during, and after it) could reduce the risk. Precipitation plays a vital role in managing water resources because it is the primary parameter to identify some natural hazards, like drought and floods (Zambrano et al., 2017). In each discipline, drought has its meaning and classifications, which leads to inventing some particular and detailed methodologies or frameworks to analyze the drought events. Some of these methodologies are the Precipitation Data Analysis (PDA) method, Water Balance (WB) method, Stream Flow Analysis (SFA) methods, and methods of using geographical and historical data. The frequently used drought analysis methods are based on the drought indices (Ellahi et al., 2020). Several drought indices have been recommended and proposed by various researchers to analyze drought events. A comprehensive discussion on few drought indices is given in (Mishra and Singh, 2010).

In drought monitoring and analysis, $\mathit{SPI}$ is an essential tool developed by (McKee et al., 1993).

Zhao et al. (2018) utilized $\mathit{SPI}$ for the Spatio-temporal evaluation of drought in Loess Plateau. Nasrollahi et al. (2018) used SPI-3 and SPI-12 at all stations to study the Spatio-temporal pattern of drought hazards in Semnan province (Iran) and constructed the Drought Hazard Index map. SPI-3 was also used by Gidey et al. (2018) to assess the meteorological drought duration, onset, offset, magnitude, severity, frequency, and spatial extent in Raya, including its environs (Northern Ethiopia), and evaluate seasonal rainfall deficit. Mishra and Singh (2010) also computed SPI-3 to analyze and discuss the agricultural droughts. The 3 and 12 months SPI was also calculated and utilized by Merabti et al. (2018) to analyze the agricultural and hydrological drought.

According to the McKee et al. (1993) classification system, the value of SPI less than or equal to −1.0 indicates the occurrence of drought events and ends when the value becomes positive. So, this limit (of SPI less than or equal to −1.0) was considered a threshold by Achcar et al. (2016) to enumerate and analyze the drought periods using the NHPP model. NHPP models are commonly used to analyze the exceedances events like noise pollutions (Guarnaccia et al., 2015), the earthquake ground-motion intensities (Iervolino et al., 2014), and the air quality standards (Achcar et al., 2016; Achcar et al., 2008; Rodrigues et al., 2015). It provides us a way towards the point pattern of exceedance events (Cressie, 1992; Diggle, 2013). The major problem in the NHPP models is to determine and estimate a suitable intensity function for appropriately representing the expected number of failures experienced up to a particular time. The NHPP has various functional forms of intensity functions under different assumptions. So, to obtain the results from NHPP models, the existing parameters of the models should be determined. For the NHPP models used for the count data analysis, the Bayesian approach is considered with various prior distributions to estimate the model parameters, commonly illustrated in previous literature (Guarnaccia et al., 2015; Leininger and Gelfand, 2017). The researcher commonly used the Bayesian approach for the estimation of parameters in various applications of different areas, i.e., in the analysis of air pollution data, in the software reliability analysis, in the drought analysis, and air quality standards (Bar-Lev et al., 1992; Kuo and Yang, 1996; Cid and Achcar, 1999; Achcar et al., 2008, 2011, 2016; Rodrigues et al., 2015) .

Geostatistical techniques are often used to illustrate spatial patterns of the random field (Entin et al., 2000; Mohanty et al., 2000). In geostatistics, a key concept is a variogram, which is used to calculate the local variability of the random field and interpolate the input parameters used for spatial interpolation through kriging (Webster and Oliver, 2007; Chen and Guo, 2017). In the modeling of variogram, the fundamental structural parameters are the sill, range, and nugget. In addition, the variogram models might contain simple models, i.e., Linear, Nugget, Spherical, Exponential, Power, and Gaussian, or the sum of one or more variogram models (Pebesma, 2004).

Here, we provide an efficient, effective, and novel framework for analyzing the drought periods to evaluate the aggregated number of drought events at unsampled locations. In the previous literature for the drought events or period analysis like Ellahi et al. (2020) utilize only the NHPP based on a single intensity function only. But the comparative study of models improves the results significantly in applied statistics. Therefore, our study is based on a comparative assessment of two NHPP models to select the appropriate model for the drought events analysis. Further, it also provides a novel way for the assessment of appropriate model parameters at unsampled locations for the analysis of accumulated drought events at that locations. The current study will support agricultural and irrigation authorities to improve the drought mitigation policies and planning. That will enhance the development of adaptation strategies for long-term drought based on the current limited resources.

The outline of our proposed study is as follows: Section 2 introduced data structure and explanations of methods utilized in our framework. Section 3 presents the descriptive structure of the agricultural drought index at observed locations, the comparison of two different models, and the presentation of the estimated accumulated number of drought events based on the most appropriately selected model. It also consists of the spatial interpolation results for the parameters of the selected model. Section 4 discussed these results, and at the end, conclusions and some future directions are discussed in Section 5.

2. Agricultural drought index formulation

SPI was originally proposed by McKee et al. (1993) and was further modified by Edwards et al. (1997). This index provides detailed information about the relative precipitation variation from normality. The 3-months $\mathit{SPI}$ is commonly used to analyze agricultural droughts; therefore, the current study is mainly concerned with analyzing the $3-$ month time scale $\mathit{SPI}.$ The computation of SPI-3 is a multi-step procedure.

For the calculation of SPI-3 used in the analysis of agricultural drought, in the initial step, summed the monthly precipitation data from month i-2 to month i. The summation value is then allocated to month i (where i = 1,2,3, …. the number of months). Therefore, the values for the first two months will be missed in the case of SPI-3.

The second step involves fitting appropriate distributions to the long-term aggregated precipitation for every considered station of our study region. The R-GUI software is used to determine the appropriately fitted distributions based on multiple goodnesses of fit criteria. The parameters of selected distributions are then estimated by using best-unbiased estimations techniques.

The cumulative probabilities for aggregated monthly average precipitation are then calculated in the next step, using the appropriately fitted distributions and their resulting parameters. The selected distributions with their parameters and estimation methods are given in Table 1.

Table 1. Fitted distributions and methods of parameters estimations.

Distribution	Parameters	Method of Parameter Estimation
Normal	σ(scale), µ(location)	Method of moments
Gamma	α(shape), β(scale)	Method of moments
Gen. Extreme Value	κ(shape), σ(scale), µ(location)	Maximum likelihood method
Gen. Pareto	κ(shape), σ(scale), µ(location)	Method of L-moments
Pearson-V	α(shape), β(scale), γ(location)	Maximum likelihood method

For SPI calculation, the standardization is done using the method proposed by Abramowitz and Stegun (1965), detailed described in Zaho et al. (2018).

The negative values of SPI-3 indicates the occurrence of drought for that time and region. On the other hand, the positive value indicates the wetness dominance for that time and region. Normal, moderate, severe, and extreme, respectively, rank the severity of wetness and drought events. This rank classification is dependent on the negative and positive values of SPI-3, which are presented in Table 2.

Table 2. Classification of droughts based on SPI.

SPI values	SPI Classes
2.0 ≤ SPI	Extremely wet
1.5 ≤ SPI ≤ 1.99	Very wet
1.0 ≤ SPI ≤ 1.49	Moderately wet
−0.99 ≤ SPI ≤ 0.99	Near normal
−1.49 ≤ SPI ≤ −1.0	Moderately dry
−1.99 ≤ SPI ≤ −1.5	Very dry
SPI ≤ −2.0	Extremely dry

2.1. Mathematical description of NHPP models

This section describes the NHPP models mathematically with their intensity and mean value functions. The data of precipitation used in the calculation of SPI-3 are collected from N meteorological stations. These are the areas where we should assess the probability that the severe events of drought occur (i.e., $\mathit{SPI}\le \mathrm{ }-1.0$) a particular number of times, in a specific time interval, along with the total average number of drought events. For the time interval $[0,T_i),$ ${\mathit{T}}_{\mathit{i}}$ > $0,$ let $K_i\ge 0$ denotes the occurrence of the total number of drought events at the location $i$ of interest. Based on K_i, the set of observed months of drought are denoted by $D_{ti}$ = $\left\{d_{1i},d_{2i},\dots ,d_{K_{i}i}\right\}.$

The value of T_i depends on the ith station; on the other hand, in our case, we consider the same value of T_i, i.e., $Ti=T$ for every $i=1,2,\dots ,N.$ The values of $K_i$ and the sets $D_{ti}$ were calculated easily from the considered stations' information where the measurements are taken. So, in our present research, we collected the set of observed data at sampled locations, denoted by $D=\left\{D_{t1},D_{t2},\dots ,D_{tN}\right\}.$ It is to be assumed that for the given time interval $[0,t)\mathrm{ }t\ge 0,$ the number of occurrences of drought events follows NHPP with an intensity (rate) and mean value function $\lambda \left(t\right)\mathrm{ }$>$0$ and $\mathit{m}(\mathit{t})={\int }_0^{\mathit{t}}\mathbf{‍}\mathit{\lambda }(\mathit{s})\mathit{ds}, \mathit{t}\ge 0.$ The cumulative drought event occurrences for a particular time interval $[0,t)\mathrm{ }t\ge 0$ at a location $i$ are denoted by $M_{ti}\ge 0.$ Let, $M_i=\left\{M_{ti}:t\ge 0\right\}$ is a NHPP with ${\lambda }_i(t)$ and $m_i(t)$ as intensity and mean value functions, respectively. So, the Poisson model becomes (2) $$P\left(M_i\right(t)-M_i(u)=k)=\frac{e^{-\left(m_i\right(t)-m_i(u\left)\right)}\left(m_i\right(t)-m_i(u){)}^k}{k!},$$(2) where $k=0,1,2,\dots$ and $t$ > $u\ge 0.$ For the prediction and interpretations of the behavior of $M_i=\left\{M_{ti}:t\ge 0\right\}$ the behavior of the ${\lambda }_i(t)$ should be understood. So our problem reduced to studying ${\lambda }_i(t),$ $t\ge 0,$ to understand and analyze the behavior of $M_i.$ ${\lambda }_i(t)$ that can be assumed in several ways. Here, we used the comparative study of two well-known cases of NHPP for the analysis of agricultural drought periods at 52 locations of Pakistan. These two special cases of NHPP models are usually applied in operating safety policy and reliability. These are the linear intensity function and the power-law process, also known as the Weibull process. The mean value functions are, (3) $${m_i}^L\left(t|{\mathit{\theta }}_i\right)={\alpha }_{i}t+\left(\frac{{\beta }_i{t}^2}{2}\right),{\mathrm{ }\alpha }_i,{\mathrm{ }\beta }_i\in \left(0,\mathrm{\infty }\right),$$(3) (4) $${m_i}^P\left(t\right|{\mathit{\theta }}_i)=(\frac{t}{{\beta }_i}{)}^{{\alpha }_i},{\alpha }_i,{\beta }_i\in \left(0,\mathrm{\infty }\right),$$(4) respectively, where, ${\mathit{\theta }}_i$=(${\alpha }_i,$ ${\beta }_i$), for $i=1,2,\dots ,N,$ is the vector of parameters that rule the behavior of the mean function. the intensity functions related to equations (3)(3) $${m_i}^L\left(t|{\mathit{\theta }}_i\right)={\alpha }_{i}t+\left(\frac{{\beta }_i{t}^2}{2}\right),{\mathrm{ }\alpha }_i,{\mathrm{ }\beta }_i\in \left(0,\mathrm{\infty }\right),$$(3) and (4)(4) $${m_i}^P\left(t\right|{\mathit{\theta }}_i)=(\frac{t}{{\beta }_i}{)}^{{\alpha }_i},{\alpha }_i,{\beta }_i\in \left(0,\mathrm{\infty }\right),$$(4) are given by (5) $${{\lambda }_i}^L\left(t\right|{\mathit{\theta }}_i)={\alpha }_i+{\beta }_{i}t$$(5) (6) $${{\lambda }_i}^P\left(t\right|{\mathit{\theta }}_i)=(\frac{{\alpha }_i}{{\beta }_i})(\frac{t}{{\beta }_i}{)}^{{\alpha }_i-1}$$(6)

In both models, the parameters vary or depend on each station individually, but the rate functions themselves remain constant. The ${\mathit{\theta }}_i$=(${\alpha }_i,$ ${\beta }_i$), $i=1,2,\dots ,N,$ represents the rate functions parameters that are related to N meteorological stations.

In the current study, we also want to investigate the number of drought events and their behavior in Pakistan for the fixed time interval of interest. Therefore, in the temporal framework for studying the behavior of drought in Pakistan, the rate functions are considered here for a suitable comparative study. So, $\mathbf{\Theta }=({\mathit{\theta }}_i=({\alpha }_i,{\beta }_i);i=1,2,\dots ,N)$ are the vectors of Poisson rate function parameters of each ith sampled location. The statisticians and mathematicians have introduced several parametric and nonparametric estimation methods to estimate the parameters. Still, here we use the Bayesian point of view for the estimation of rate function parameters. For $i^{th}$ the location, ${\mathit{\theta }}_i$=(${\alpha }_i,$ ${\beta }_i$) is an unknown vector of parameters (in both cases of rate functions), estimated under some uncertainties, described by assigning the prior distributions in the Bayesian approach. So our problem is further reduced to estimate ${\mathit{\theta }}_i$=(${\alpha }_i,$ ${\beta }_i$), for every $\mathit{i}=1,2,\dots ,\mathit{N}.$ These estimated parameters represent the most appropriate behavior to present each station's drought measurements.

2.2. Bayesian formulation

This section presents the general mathematical framework of the Bayesian approach to analyze the drought periods. The joint posterior distribution of the vector ${\mathit{\theta }}_i$ is the distribution of interest for each station. This joint probability distribution is obtained by utilizing the likelihood function of the NHPP model and the prior distribution of ${\mathit{\theta }}_i$ (i.e., the prior joint distribution of the parameters ${\alpha }_i$ and ${\beta }_i$). The posterior distribution is $$P\left({\mathit{\theta }}_i|D_{ti}\right)\propto L\left({\mathit{\theta }}_i|D_{ti}\right)P\left({\mathit{\theta }}_i\right),$$

Where $L({\mathit{\theta }}_i|D_{ti})$ is the likelihood function and $P({\mathit{\theta }}_i)$ is the prior distribution of ${\mathit{\theta }}_i.$

So, according to our hypothesis, the NHPP models are as follow; (7) $$L\left({\mathit{\theta }}_i|D_{ti}\right)=\left(\prod _{k=1}^{K_i}\mathrm{‍}{\lambda }_i\left(d_{ki}|{\mathit{\theta }}_i\right)\right)e^{-\left[m_i\left(T|{\mathit{\theta }}_i\right)\right]},$$(7) here, we used ${\lambda }_i(t|{\mathit{\theta }}_i)$ and $m_i(t|{\mathit{\theta }}_i)$ to make it unequivocal the reliance of rate and mean function on the corresponding parameters ${\mathit{\theta }}_i$ (Cox and Lewis, 1966).

The MCMC methods are used to generate the posterior summaries of interest. In the Bayesian approach for the above-discussed models, it is considered that the priors for parameters ${\alpha }_i$ and ${\beta }_i$ are uniform distributions, i.e., ${\alpha }_i\sim U[a_1,b_1]$ and ${\beta }_i\sim U[a_2,b_2].$ The hyperparameter $(a_1,b_1,a_2,b_2)$ values for prior distributions can be considered by the personal judgment of the analyst. Those values of hyperpriors should be considered for which the prior distributions have minimum variance. Here we assume the prior independence between the parameters. The simulated samples of the joint posterior distribution ${\mathit{\theta }}_i$ are obtained by using the MCMC methods. For this purpose, we used the library "R2jags" under R software originally introduced by (Su and Yajima, 2012).

In this work, the prior distributions of hyperparameters are supposed to be known that will be described later. The most appropriate model is selected using the existing Bayesian adequacy measures called DIC; the smaller the DIC, the better the models will be. The performance of the model is observed by comparing the graphs of estimated and actual mean values of the cumulated number of drought events versus the month of occurrence. After selecting the most appropriate model, each station's estimated parameters are regionalized to observe the spatial behavior of the parameters. For this purpose, we used Ordinary Kriging (OK) (see, Malvic, 2008, which is considered as the Best Linear Unbiased Predictor (BLUP).

3. Estimation results

Here, the results of the above-discussed models are presented in this section. The methodology provided in the previous sections is implemented on the average monthly precipitation data of 52 observed locations of Pakistan. Since the objective of our paper is to estimate and investigate agricultural drought events therefore precipitation data is used to calculate $3-$ month $\mathit{SPI}$ values at each sampled location. Figure 1 presents the time series plots of SPI-3 from 1968 to 2016.

Fig. 1. The 3-Month standardized precipitation index (SPI −3) versus time (date) (a) Kotli, (b) Dir, (c) Kalat and (d) Lasbella.

3.1. Agricultural drought

The atypical behavior of SPI-3 (i.e., values of SPI-3 is less than or equal to −2.0) in Fig. 1 indicates the severely dry periods for Dir and Kotli. In Dir, the longest seasonal or agricultural drought period was from August 1969 to July 1971. That was 23 months drought period with −2.3, the peak value. Other severe droughts of 1, 2, 3, and 4 months with peak values ranging from $-2.02$ to $-3.31$ are also observed in Dir. Figure 1a shows the values of SPI-3 of Kotli with temporal variations. It indicates 2, 3, 5, and 7 months of severe agricultural drought periods with different peak values ranging from −2.04 to −3.56. Figure 1c and d indicates that no single month of severe drought exists from 1968 to 2016 for Kalat and Lasbella.

The cumulative number of drought events in each time versus month of occurrence is presented in Fig. 2 for the SPI-3 series($\mathit{SPI}-3\le -1.0$) of Kotli, Dir, Kalat, and Lasbella. The moderate drought periods, other than the severe drought periods, at each station can also be identified from Fig. 1 with various duration and peak values. The moderate and severe drought events classification is provided in Table 2. The summaries of these drought periods are given in Table 3.

Fig. 2. Accumulated numbers of SPI − 3 ≤ −1.0 in each time versus month of © (a) Kotli, (b) Dir, (c) Kalat and (d) Lasbella.

Table 3. Drought periods of each station for 1968-2016 by SPI −3(≤ −1.0).

Kotli					Dir
Start date	End date	Duration	Peak	Average	Start date	End date	Duration	Peak	Average
01/09/1968	01/12/1968	3	−1.13	−0.63	01/07/1969	01/07/1971	24	−2.3	−1.18
01/04/1970	01/06/1970	2	−1.41	−1.38	01/11/1971	01/02/1972	3	−1.58	−1.01
01/01/1971	01/06/1971	5	−1.99	−1.12	01/12/1973	01/04/1975	16	−1.71	−0.74
01/11/1971	01/01/1972	2	−1.42	−0.89	01/12/1976	01/01/1977	1	−1.27	−1.27
01/05/1974	01/02/1975	9	−1.91	−0.72	01/04/1977	01/06/1977	2	−1.87	−1.49
__	__	__	__	__	__	__	__	__	__
__	__	__	__	__	__	__	__	__	__
__	__	__	__	__	__	__	__	__	__
01/07/2015	01/10/2015	3	−2.23	−0.88	01/01/2014	01/03/2014	2	−1.45	−1.23
Lasbella					Kalat
Start date	End date	Duration	Peak	Average	Start date	End date	Duration	Peak	Average
01/08/1968	01/12/1968	4	−1.5	−0.96	01/03/1968	01/06/1969	15	−1.54	−0.82
01/08/1971	01/06/1972	10	−1.08	−0.57	01/01/1971	01/12/1971	11	−1.61	−0.93
01/09/1972	01/12/1972	3	−1.27	−0.78	01/04/1973	01/01/1975	21	−1.73	−0.9
01/07/1974	01/12/1974	5	−1.53	−0.96	01/02/1977	01/04/1978	14	−1.27	−0.56
01/07/1975	01/01/1976	6	−1.02	−0.49	01/01/1979	01/10/1980	21	−1.53	−0.94
__	__	__	__	__	__	__	__	__	__
__	__	__	__	__	__	__	__	__	__
__	__	__	__	__	__	__	__	__	__
01/02/2016	01/05/2016	3	−1.15	−0.46	01/02/2016	01/04/2016	2	−1.34	−0.67

We could consider the Stochastic Volatility(SV) models (see, Ghysels et al., 1996; Meyer and Yu, 2000) under the Bayesian hierarchical approach to model time series SPI-3 data. The drought events can occur any time the SPI-3 continuously negative and can reach an intensity of −1.0 or less. In the present paper, we considered the time(months) of droughts identified by the threshold $L=-1.0$ for all stations from January 1968 to December 2016.

The models discussed in the previous sections (2.2 and 2.3) are applied to the cumulative number of drought events concerning months at each station. For this model approximation, the Bayesian framework is used to estimate parameters for studied models at each station for our study region. We know, according to the classification of SPI values, the drought event occurs if $\mathit{SPI}\le -1.0.$ Therefore we use the threshold $L=<-1.0$ in our study. Many researchers in their research studies recently considered similar threshold values (e.g., see, Achcar et al., 2016). The aggregated drought events at each station for the period 1968–2016 are given in Table 4. For those stations, these aggregated values correspond to the values of $K,$ respectively. The study period corresponds to a total number of months, $T_i=T=588,$ for ($i=1,2,\dots ,N=52).$ In the case of the power-law process, the uniform non-informative priors are assumed for ${\alpha }_i$ and ${\beta }_i,$ i.e., ${\alpha }_i\sim U[0,100]$ and ${\beta }_i\sim U[1,100].$ The uniform non-informative priors are also assumed for ${\alpha }_i$ and ${\beta }_i,$ i.e., ${\alpha }_i\sim U[0,100]$ and ${\beta }_i\sim U[0,100]$ when the linear intensity model is considered. We used those hyperparameter values for which the prior distributions have minimum variance. There is not any hard and fast rule to define hyperparameter's values; these values can be considered based on personal judgment. For the joint posterior distribution of interest, we produced 31,000 Gibbs samples from which the first 1000 were discarded to eliminate or reduce the effect of initial values used in the iteration process in the simulation sample.

Table 4. List of the total number of agricultural droughts at each station (1968-2016).

Stations	Total Events	Stations	Total Events	Stations	Total Events	Stations	Total Events
Cherat	80	Lasbella	69	Padidan	66	Sargodha	98
Peshawar	94	Zohub	92	Rohri	24	Murree	94
Kakul	108	Nokkundi	24	Badin	33	Mianwali	88
DIK	82	Khuzdar	81	Chhor	28	Bahawalnagar	80
Drosh	97	Barkhan	70	Karachi	69	Gilgit	91
Parachinar	81	Faisalabad	93	Nawabshah	50	Sakardu	86
Kohat	102	Jhelum	101	Dalbadin	39	Bunji	98
Risalpur	90	Sialkot	87	Panjgur	42	Chilas	84
Balakot	95	BahawalPur	79	Pasni	33	Astore	95
Chitral	100	Multan	83	Kalat	39	Gupis	73
Dir	90	Khanpur	53	Quetta	56	Kotli	87
Hyderabad	44	Lahore	89	Sibbi	63	Muzafarabad	88
Jacobabad	63	Rawalpindi	98	Jiwani	22	Garhi Dupatta	92

3.2. Model comparison results

It might be at some individual stations that the linear intensity function performs better than the power-law process, but, collectively the power law process model performs better than the linear intensity function. This is confirmed by Monte Carlo estimates for the DIC of each model. The DIC values for the power-law process and linear intensity function models are given in Table 5. As the DIC value of the linear intensity function is greater than the value of the power-law process so, in our case, the power-law process is better than the linear intensity function. Therefore, the estimated parameters of the power-law process are considered to estimate the accumulated number of agricultural droughts at each observed location and time point. The selected model parameters are also estimated at unobserved locations to observe their spatial behavior by using spatial interpolation tools. Figure 3 presents the observed accumulated number $M_t,t\in [0,T]$ of drought events and the estimated mean values of accumulated droughts versus four representative stations using the power law process model.

Fig. 3. Plots of Accumulated numbers of SPI −3 ≤ −1.0 in each time with the estimated mean value function versus month of ©. (a) Kotli, (b) Jiwani, (c) Mianwali, and (d) Barkhan.

Table 5. Values of DIC for power-law process and linear intensity function models.

Models	DIC
Linear Intensity Function	189044
Power Law Process	173266

3.3. Estimation results

The estimates of the parameters ${\alpha }_i$ and ${\beta }_i$ for each location and their 95% credible intervals (the posterior summaries) provide some useful information about the behavior of each parameter and its effect on the power-law model. It has been observed that for the power-law process, the intensity function ${\lambda }_i(t)$ has flexible behavior due to the value of ${\alpha }_i.$ For ${\alpha }_i>1,$ the intensity function is increasing, and for ${\alpha }_i<1$ the intensity function is decreasing and constant for the value of ${\alpha }_i$=1, i.e., the NHPP transformed to HPP. In our study area, It has been observed that ${\alpha }_i$ have similar behavior for the stations of the same geographic region. Moreover, it is to be noted that the stations in the different regions may also have similar behavior. The station of Jewani ${\alpha }_i$ has a maximum value that is 1.704$(>1)$ with a standard deviation of 0.01163, and the 95% credible interval is [1.68, 1.726]. On the other hand, Kalat has a minimum value of ${\alpha }_i$ which is 0.5853 $(<1)$ have respective standard deviation is 0.004609 and the 95% credible interval is [0.5792, 0.597]. It is observed that the larger the value of the parameter $\alpha ,$ the larger will be the standard deviation, and the wider will be the credible interval. Likewise, the smaller value of $\alpha$ has a smaller standard deviation and a little credible interval. It is more important to be noted that the posterior mean of the parameter ${\alpha }_i<1$ for any station implies that the rate at which the drought event occurs is decreasing with respect to time. Therefore, the drought occurrence at the threshold less than or equal to −1.0 is becoming less frequent.

For some stations, the values of the parameter ${\alpha }_i$ may vary significantly even though the stations or locations exist in the same geographical region. The rate function of some stations has similar behavior, which has a similar value of ${\alpha }_i.$ Jewani station also has a maximum value for parameter ${\beta }_i$ with 0.4641 standard deviations and a 95% credible interval is $[98.31,\mathrm{ }99.99].$ However, Mianwali and Bahawalpur both have a minimum (and same) value of ${\beta }_i$ equals to 1.004 and have the same 95% credible interval $[1,\mathrm{ }1.014]$ with their respective standard deviations 0.003669 and 0.003884. The ${\beta }_i$ values indicate that if the drought occurrence is increasing or decreasing, they increase or decrease in different ways. If larger (smaller) the value of ${\beta }_i,$ smaller(larger) will be the value of ${\lambda }_i(t)$ for the fixed value of ${\alpha }_i$ and $t\ge 0.$ So, the considered appropriate model captured the different behavior of the drought occurrence rate at different locations.

The estimated and observed accumulated means of drought occurrence for four stations are presented in Fig. 3. Because for all the stations, these graphs are categorized in three types, i.e., good, medium, and worse estimated, based on estimated accumulated means of drought events. Those estimated means are calculated at each point when the drought events occurred for each station. The categories of those fitted graphs are decided based on their graphical representation and RMSEs. To display all these three categories of stations for our fitted models, we selected four stations, i.e., Kotli, Dir, Lasbella, and Kalat. The observed and estimated accumulated means of these selected stations are presented in Fig. 3. Here, the accumulated mean stand for the mean value function $m_i(t)$ either observed or estimated, calculated at each time where the drought event occurred.

The accumulated mean of drought events at Kotli shown in Fig. 3a is a good fit to estimate the observed mean with RMSE 2.13576. Figure 3b shows the poor fit of the estimated accumulated mean to the empirical average for Dir, with RMSE 6.56612. The accumulated means of drought events for Kalat and Lasbella are shown in Figs. 3c and 3d, which is the average fit of estimated to the observed means with RMSEs 3.470529 and 3.450207. Like Kalat, It may be observed that a model taking into consideration can capture the changepoint, or any alternate rate function could be considered. Similarly, for stations Dir and Lasbella, maybe a model either with an alternate rate function or with the presence of more than one changepoints could give a superior fit. Table 6 presented the RMSEs at each study region.

Table 6. RMSEs of each station for observed vs estimated accumulated drought events.

Stations	RMSE	Stations	RMSE	Stations	RMSE
Cherat	2.8630	Nawabshah	3.1398	Khanpur	1.7831
Peshawar	2.8028	Dalbadin	3.8555	Lahore	4.8889
Kakul	5.8521	Panjgur	3.5737	Rawalpindi	2.4249
DIK	3.2788	Pasni	2.0577	Sargodha	2.7688
Drosh	2.4722	Kalat	3.4705	Murree	6.8652
Parachinar	4.1849	Quetta	4.2837	Mianwali	5.6827
Kohat	4.5213	Sibbi	4.6927	Bahawalnagar	3.5966
Risalpur	2.9361	Jiwani	1.5600	Gilgit	4.4050
Balakot	4.1720	Lasbella	3.4502	Sakardu	2.1707
Chitral	3.2353	Zohub	3.1586	Bunji	3.2027
Dir	6.5661	Nokkundi	1.2927	Chilas	6.3755
Hyderabad	2.6348	Khuzdar	3.4265	Astore	4.0773
Jacobabad	3.3215	Barkhan	4.9927	Gupis	8.4248
Padidan	3.7135	Faisalabad	4.2647	Kotli	2.1358
Rohri	0.9458	Jhelum	2.6521	Muzaffarabad	3.8791
Badin	1.8745	Sialkot	3.7204	GarhiDupatta	3.7050
Chhor	1.9242	BahawalPur	5.4079
Karachi	2.5290	Multan	2.1890

3.4. Regionalization and interpretation results

In this section, we implemented the above-discussed models for the regionalization or spatial distribution of power-law process model parameters ($\alpha$ and $\beta$) by using the estimated posterior mean values of these parameters for 52 locations in Pakistan. These model parameters were used to fit the accumulated number of drought events evaluated by using $\mathit{SPI}-3.$ In the next sub-sections, we deal with the parameters $\alpha$ and $\beta$ results individually.

3.4.1. For parameter $\mathit{\alpha }$

In the case of parameter "λ" for normalization, the Box-Cox transformation technique with an estimated transformation parameter $\lambda =-0.8689718$ is considered for our data. The variogram model can be calculated with either modulus or classical procedure. We calculated spatial variogram of $\alpha$ by the method of the modulus. Figure 4 shows the graphical assessment of the empirical variogram for ordinary kriging calculated for the $\alpha$ parameter. Before applying the spatial interpolation technique, it should be noted that the spatial correlation must exist in the data. Suppose the spatial correlation does not exist in sample point data. In that case, we cannot use any of the spatial interpolation techniques and need to search for another interpolation method for this data.

Fig. 4. Observed variogram, variogram envelope, and fitted variogram for $\mathrm{\alpha }.$

Various techniques have been introduced in the literature to investigate the spatial correlation of sampled point data. Here we used the variogram envelope technique to check the spatial correlation in the data. Figure 4a represents the graphical representation of the variogram envelope. Dashed lines are the 95% pointwise variogram envelope, based on the $2.5^{th}$ and $97.5^{th}$ percentile of semi-variance estimates. The graph clearly shows that most of the variogram points lie in the envelope and four points lying on the envelope lines; to some extent, these points are also considered in the variogram envelope range. So, it indicates that there is a strong spatial correlation in the sampled points.

For spatial interpolation and structural analysis, variogram modeling and estimation methods used to estimate variogram model parameters are essential (Burrough and McDonnell, 1998). So, for this purpose, various estimation techniques are available in the literature. In this study, the Ordinary Least Square (OLS) estimation, Maximum Likelihood (ML) estimation, Weighted Least Squares (WLS) estimation, and Restricted Maximum Likelihood (REML) estimation are utilized for fitting the variogram model. Various analytical functions (models) are used to fit the variogram model to the experimental variogram, parameters estimation techniques (ML, REML, OLS, and WLS), and their Root Mean Squared Errors (RMSE's) are given in Table 6. The main aim of using various analytical functions (models) and parameter estimation methods is to find a suitable model and appropriate estimate of parameters for that model.

The model estimated results were validated by using the leave-one-out cross-validation technique. In Table 7, the minimum RMSE= 0.19201 is for the wave model, and the estimated method used is REML. The plot of an appropriate fitted model is shown in Fig. 4b. So, we used this model and parameter estimation technique for fitting the empirical variogram and used that fitted model in OK for spatial interpolation of $\alpha$ parameter values.

Table 7. Parameters and RMSE's of variogram models in case of $\mathrm{\alpha }.$

Methods	Models	Nugget	Sill	Range	RMSE	Methods	Models	Nugget	Sill	Range	RMSE
ML	Metern	0.0311	0.0519	5.4865	0.201258	OLS	Metern	0.0327	32.6976	5313.904	0.19324
	Spherical	0.0331	0.0427	10.2613	0.199141		Spherical	0.0328	15.409	3769.407	0.19326
	Cubic	0.0331	0.036	5.3648	0.208773		Cubic	0.0491	157.193	1722.901	0.22421
	Exponential	0.0311	0.0519	5.4865	0.201258		Exponential	0.0327	32.6976	5313.904	0.19324
	Cauchy	0.0397	0.0451	454.329	0.197553		Cauchy	0.0492	268.05	602.6732	0.22515
	Wave	0.0401	0.0277	2.0221	0.20096		Wave	0.0491	64.4585	170.6228	0.22516
REML	Metern	0.0313	2.2563	295.287	0.192497	WLS	Metern	0.0249	0.0768	6.1163	0.19979
	Spherical	0.0319	0.0524	10.5521	0.198267		Spherical	0.024	0.0626	9.2055	0.20142
	Cubic	0.0329	0.0413	5.6293	0.208115		Cubic	0.0331	0.0544	11.1678	0.20452
	Exponential	0.0313	2.2563	295.287	0.192497		Exponential	0.0249	0.0768	6.1163	0.19979
	Cauchy	0.039	0.2022	8.0518	0.199316		Cauchy	0.0355	0.0806	3.7688	0.19907
	Wave	0.0401	0.0368	2.1237	0.192013		Wave	0.0327	0.0472	1.8221	0.21026

Finally, efforts are made to predict the parameter $\alpha$ at un-sampled locations by using the above appropriate selected variogram model and estimation technique in OK that provides the map of Pakistan for the spatial interpolation results with their associated prediction errors. These maps are shown in Fig. 5. The contour map of spatial prediction of parameter $\alpha$ posterior mean values is presented in Fig. 5a, and Fig. 5b shows the spatial interpolation errors map for the parameter $\alpha .$ These maps only show an understanding of the behavior or, to some extent, the confidence intervals of predicted values and their prediction errors in our study area. Using the primary and prominent colors in the construction of contour maps, we can figure out the estimated results for specific points or locations because color fluctuates at every degree. From Fig. 5a and b, the contour intervals for the estimated or predicted $\mathit{\alpha }$ with their prediction or estimation errors are figured out with color keys. Contour intervals for predicting $\alpha$ and prediction errors are approximately [0.7, 1.4] and [0.11, 0.37] i.e., the interval lengths are 0.7 and 0.26. The predicted mean values of the parameter $\alpha$ and their standard deviations have minimum values in the area under latitude [26.5, 28.5] and longitude [69.5, 71.8].

Fig. 5. Spatial interpolation and associated prediction error for parameter $\mathrm{\alpha }.$

Similarly, The predicted mean values of the parameter $\alpha$ and their standard deviations have their maximum values in the area under latitude [25, 26.2] and longitude [61.8, 62.8]. In both cases, i.e., predicted $\alpha$ parameters and their standard error contour maps, the green color represents the minimum values, and the red color shows the maximum values. The prominent green, blue, yellow, and red-colored areas have $\alpha$ values 0.7, 0.95, 1.15, and 1.35 with standard deviations of 0.12, 0.20, 0.28, and 0.36, approximately.

3.4.2. Regionalization and interpretation of parameter $\mathit{\beta }$

Before applying the spatial interpolation technique to examine the spatial correlation of sampled point data, several techniques have been introduced in the literature. The spatial correlation is checked by using the variogram envelope technique. Figure 6a provides the graphical presentation of the variogram envelope. The dashed lines are 95% pointwise variogram envelope, based on the $2.5^{th}$ and $97.5^{th}$ percentile of semi-variance estimates. From the figure, it has been observed that a maximum number of variogram points lying in the envelope and two points lying on the envelope lines; to some extent, these points are also considered in the range of variogram envelope. So, it is an indication of the existence of a strong spatial correlation in sampled points. The distribution of the parameter $\beta$ was non-normal, so we transformed the data parameter $\beta$ to make its behavior normal. For normalization, the Box-Cox transformation technique with an estimated transformation parameter $\lambda =-0.2229$ is considered in our data.

Fig. 6. Observed variogram, variogram envelope, and fitted variogram for $\mathrm{\beta }.$

Here, the modulus method is used for the calculation of empirical variogram for "$\beta$ ". Figure 6 shows the visual assessment of the empirical variogram calculated from the $\beta$ parameter values. After plotting the experimental variogram of the data, different variogram models were fitted to the experimental variogram as the experimental variogram model does not give us better structural analysis and spatial interpolation of the study variable at un-sampled locations. In this study, the parameters of variogram models were estimated by using ML, REML, OLS, and WLS. The fitted variogram model having the smallest RMSE was used to interpolate the parameter $\beta$ at unobserved locations. Different analytical function models fitted to the experimental variogram, their estimated parameters, parameter estimation techniques, and their Root Mean Squared Errors (RMSE's) are given in Table 8. The main aim of using various analytical functions (models) and parameter estimation methods is to find the appropriate variogram model and the appropriate estimate of parameters for that model.

Table 8. Parameters and RMSE's of variogram models in case of $\mathrm{\beta }.$

Methods	Models	Nugget	Sill	Range	RMSE	Methods	Models	Nugget	Sill	Range	RMSE
ML	Metern	0.3249	0.4358	6.6889	13.0437	OLS	Metern	0.1322	11658	130663.6	8.63995
	Spherical	0.3365	0.3411	10.5881	13.1787		Spherical	0.1327	783.9052	13187.99	8.6436
	Cubic	0.3815	0.3386	13.9475	12.5888		Cubic	0.3491	4179.781	2284.168	14.6308
	Exponential	0.3249	0.4358	6.689	13.0437		Exponential	0.1326	5841.627	65506.27	8.6433
	Cauchy	0.3867	0.3574	377.526	12.3193		Cauchy	0.3506	3126.162	529.8673	14.7507
	Wave	0.3806	0.3541	5.8317	14.7532		Wave	0.35	909.1843	164.8734	14.7566
REML	Metern	0.3252	21.6513	400.669	10.4671	WLS	Metern	0.3569	1185.485	26079.79	10.9931
	Spherical	0.3251	8.8044	245.01	10.4208		Spherical	0.3573	118.7585	3925.583	10.9984
	Cubic	0.3806	0.5533	16.8201	10.7795		Cubic	0.4651	1780.037	1934.114	14.7747
	Exponential	0.3252	22.7968	422.036	10.4655		Exponential	0.3569	1303.168	28669.36	10.9931
	Cauchy	0.3804	1.4352	7.915	10.0929		Cauchy	0.4657	1422.229	463.7862	14.8388
	Wave	0.3987	0.3486	2.4321	11.7545		Wave	0.4655	440.6046	149.0217	14.8421

The leave-one-out cross-validation technique was used to validate the estimated model results. The criteria for selecting an optimal model is set by RMSE, i.e., the smaller the RMSE better will be the interpolated results. In Table 7, the minimum $\mathit{RMSE}=8.639952$ is for the Matern model, and the most appropriate estimation method used is OLS. Figure 6b shows the plot of the suitable fitted variogram model. So, we used this model and param eter estimation technique for fitting the experimental variogram and used that fitted model in OK for spatial interpolation of $\beta$ parameter values at unobserved locations.

Finally, the above appropriate selected variogram model and estimation technique were used in OK to predict the parameter $\beta$ at unobserved locations for the distribution maps of the spatial interpolation results and associated prediction error. Figure 7 shows the spatial distribution of the parameter $\beta$ and their prediction errors, i.e., standard deviations. The contour map of spatial prediction of parameter $\beta$ posterior mean values is presented in Fig. 7a, and Fig. 7b shows the spatial interpolation errors map of the parameter $\beta .$ Prominent colors, i.e., green, blue, yellow, and red, were used to construct contour maps, which helps to figure out the estimated results for specific points or locations because color fluctuates at every degree.

Fig. 7. Spatial interpolation and associated prediction error for parameter $\mathrm{\beta }.$

Contour intervals of predicted $\beta$ and prediction errors are approximately [0, 145] and [0, 90], i.e., the interval lengths are 145 and 90. The predicted mean values of the parameter $\beta$ and their standard deviations have minimum values for the maximum study area, i.e., the area under latitude [27, 31] and longitude [68, 74] with their surroundings. Similarly, The predicted mean values of the parameter $\beta$ and their standard deviations have their maximum values in the area under latitude [25, 30] and longitude [61.5, 62.5]. In both cases, i.e., predicted $\beta$ parameters and their standard error contour maps, the green color represents the minimum values, and the red color shows the maximum values. The prominent green, blue, yellow, and red-colored areas have $\beta$ values 5, 50, 95, and 140 with standard deviations 0, 30, 60, and 90, approximately.

4. Discussion

Drought is a natural disaster, and it differs from other natural disasters like floods, earthquakes, and tornadoes, etc. Almost every region of Pakistan has a history of facing severe drought, which substantially damaged the industrial, agricultural, and economic sectors (Hussain and Mumtaz, 2014). The paper provides an assessment of agricultural drought, as categorized by the $\mathit{SPI}-3.$ The advantage in consideration of SPI for drought monitoring is that it can easily estimate the precipitation deficiency at multiple time scales. These time scales measure the consequences of drought events on different water resources. Therefore, in various countries around the world, SPI is now used in multiple operational monitoring systems. It is also recognized by the World Health Organization (WHO). Therefore, this research study may provide a better understanding for policymakers, which helps construct effective drought mitigation and adaptation policies.

In our study, the time series plots of $\mathit{SPI}-3$ (see, Fig. 1) presents the descriptive structure of agricultural drought indices. Table 3 provides the complete and detailed list of start and end dates of agricultural drought periods with their durations (magnitudes), peak and average intensities of agricultural droughts that occurred at Kotli, Dir, Lasbella, and Kalat. Whereas Table 4 has the list of aggregated agricultural drought events at each region or station under study from 1968 to 2016. The results of these tables indicate that every location or region of our study area has a history of severe agricultural drought events. Achcar et al. (2016) used the NHPP for the drought periods analysis at a single station by using the power law process. But here in our study, we compare the two different NHPP approaches, i.e., the Power-law process and the Linear intensity function, and select the appropriate model to analyze the agricultural drought events at 52 meteorological stations. For this purpose, The appropriate model is selected using DIC, which indicates that power-law process performs better than the linear intensity function (as the DIC value of the power-law process model is less than the linear intensity function model) to evaluate the cumulative drought events in the study area(see, Table 5). The parameters were estimated using the joint posterior distribution of ${\mathit{\theta }}_i$ under MCMC simulation. The means of marginal posterior distributions for each parameter of the power-law process at each station are used to estimate and predict the agricultural drought events at that station. Here, we presented and discussed the observed and estimated four stations' drought events at each time point by the power-law process model. The RMSEs and graphical presentation in Fig. 3 provide a clear image of the best, average, and worse estimated results of the accumulated number of drought events models. Table 6 presented the RMSEs of the fitted model for all 52 regions. The presented results showed that approximately for maximum stations, the model performed significantly well. The worse performance of the model is due to the abrupt changes in the behavior of drought occurrence events at some locations. At some regions where the performance of the model was not as satisfactory is due to the abrupt changes in drought occurrences could be improved, which may be considered in our future project. The novelty in this research work is the use of the geostatistical interpolation method to interpolate the appropriate model (i.e., power-law process) parameters for the study area. The estimation of variogram models, including their parameters, is essential for spatial prediction and structural analysis. Different variogram models were estimated and fitted to the experimental variogram by using ML, REML, OLS, and WLS. The variogram model estimated results were validated by using the Leave-one-out cross-validation technique. The RMSE is used as a measure of selecting the optimal model. The wave model with the REML estimation method has minimum RMSE (i.e., 0.19201, see Table 7 and Fig. 4b).

On the other hand, in the case of the parameter $\beta ,$ the Matern model with the OLS estimation method has minimum RMSE (i.e., 8.639952, see Table 8 and Fig. 6b). These estimated variogram models are then used in OK to predict the parameters $\alpha$ and $\beta$ at un-sampled locations that provide the detailed structures or spatial variation mirrors (contour maps) of these parameters by generating an accurate and reliable spatial distribution. From these maps, we can observe a clear view of the shape and scale of the selected model over the study area (see Figs. 5 and 7). The rate of change of parameter $\alpha ,$ which is the shape parameter of selected models, is very high in the Western Balochistan region. On the other hand, the Eastern Balochistan, the Western highlands, the sub-Montane, the greater Himalayas, and the lower Indus areas have lower rates of variation of the shape parameter $\alpha ,$ whereas the remaining study area, which is the central and Southern Punjab, has a nominal rate of spatial variation with minimum values of the parameter $\alpha .$ The contour maps of their standard deviation also present approximately the same spatial behavior.

Similarly, the rate of variation of parameter $\beta ,$ which is the scale parameter of models, has the highest rate of variation in the area of the Western Balochistan plateau. In other regions, the variation rate of scale parameters is very low, which is negligible and has minimum values of the parameter $\beta .$ The contour maps of their prediction errors also have the same spatial behavior.

5. Conclusion

Natural disasters cannot be stopped; every disaster provides information to ensure well-planned management and useful information for preventive measures. This paper mainly analyzed the agricultural drought periods and regionalized the appropriately selected model parameters. It delivers a better research framework or structure in understanding the agricultural drought with various intensities at multiple locations. It also provides a framework for the analysis and estimation of cumulative agricultural drought events at each location by an appropriately selected model and regionalizes the selected model parameters for the study area. For this purpose, our research work is divided into three essential phases; each phase has its concept. In the first phase, the concept of $3-$ month time scale $\mathit{SPI}$ is used to calculate and exploratory analyses the agricultural droughts with their intensities at each location. That performs well in presenting and estimating the agricultural drought events with their intensities. At the second phase, the concept of point processes is utilized, i.e., NHPP with linear intensity function and Weibull process. The most appropriate model, i.e., the Weibull process model, is selected using DIC to analyze and estimate the accumulated drought events at each time point for each station, which is obtained by introducing a threshold in drought intensities. In the end, the geostatistical interpolation technique (OK) is used to regionalize the appropriate model parameters which are very useful to estimate the parameters at unsampled locations. The spatial variations of the estimated parameters and their respective errors are also presented in contour maps.

We have identified the fruitful directions for the associated future research built on the work developed in this paper. There are no criteria for comparing the drought indices because every drought index is based on a different conceptual framework. However, more climatic variables in calculating the drought index can increase the drought index accuracy to perform better in drought prediction and monitoring. A class of geostatistical models can be used to determine the comparative study to analyze the spatial behavior of model parameters, which will improve the efficiency and accuracy of the proposed research.

Competing interests

The authors declare there are no competing interests.

Consent to participate

All authors are voluntarily agreed to participate in this research study.

Ethical statement

All procedures followed were in accordance with the ethical standards with the Helsinki Declaration of 1975, as revised in 2000.

Consent to publish

All authors are agreed for publication; there is no legal constraint in publishing the data used in the manuscript.

Availability of data and materials

The data used for the preparation of the manuscript is available with the corresponding author and can be provided upon request.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

The author extends his appreciation to the Deanship of Scientific Research at King Khalid University for funding this work under grant number (RGP.1/26/42), received by Mohammed M. Almazah (www.kku.edu.sa).