Hydrological Drought Index Based on Discharge

Abstract

Drought is a natural phenomenon causing disasters and its period of occurrence can be predicted in recent times based on several methods using the same or different variables. The prediction is usually associated with the climate interactions in the form of rainfall or discharge patterns which can be analyzed using the return period. Therefore, this research was conducted in four different stages of data acquisition and validation, drought analysis method based on the data, drought prediction method based on hydrology, and sample applications to determine the debit availability in other watersheds. Historical rainfall data converted to dependable rainfall at 80% probability were used as input for the rainfall-discharge analysis while the hydrological drought analysis was conducted using the drought threshold value. Moreover, the drought was predicted using an artificial neural network model while historical data were used to verify the hydrological character of the prediction model. The results of the analysis conducted were further used to predict the water balance in different river areas due to the fact that each area has a different hydrological character. Meanwhile, the watersheds used as case research showed that the model has reliability of up to 80%.

1. Introduction

Drought is a condition when there is very little extreme rainfall or no rainfall for a relatively long period outside the dry season [1, 2]. Its occurrence in a watershed is recommended to be analyzed using discharge data in addition to the rainfall data [3]. Moreover, the hydrological drought index (HDI) is normally used to describe the drought severity based on discharge data which is indicated as HDI₁, and simulation which is denoted as HDI₂. The discharge prediction model, which is known as the Artificial Neural Network model, is based on a flexible mathematical structure and has the ability to identify a complex nonlinear relationship between input and output [4]. The characters of the rainfall and discharge data were also observed to be adaptable to the method. The purpose of this research is to answer the following questions:

1. What is the estimated discharge calculation model based on rainfall data?

2. How can the dry duration and deficit be determined?

3. What is the hydrological drought index (HDI) in the watershed of the research location?

4. What are the HDI criteria in the watershed of the research location?

2. Literature review

2.1 Data statistical analysis

The rainfall data consistency was tested using multiple mass curves [5, 6], while the discharge data normality test was conducted using probability density function (pdf) analysis with the criterion being that normally distributed data will form a linear or straight line. Moreover, the discharge probability was analyzed using the cumulative distribution function (CDF).

The statistical analysis is usually conducted to determine:

1. Threshold (X0) is the limit determined based on the analysis requirements [7] and according to the selected distribution.

2. X0 as Q50 or Q80 such that Q50 is a normal Q with a probability of 0.5 or data median, while Q80 is a dependable Q with a probability of 0.2.

The criteria are divided into two, which are as follows:

• The value is dry (D) and wet (W) when Q80 < Q [7].

• The value is extremely dry (ED) and very dry (VD) when Q < Q80 [8].

2.2 Calculation of hydrological drought index (HDI) in data statistical analysis

The drought index is a comparison of the deficit to the watershed area as indicated in the following Equation [8]:

$$\mathit{HDI}=\frac{\mathit{\text{deficit}}}{\mathit{\text{Areal}}}\frac{m^3/\mathit{sec}}{{\mathit{km}}^2}$$E1

HDI is the hydrological drought index, the deficit is the difference between X₀ and X_t (ten daily), X₀ is the dry threshold, and X_t is the ten daily periods of discharge.

Drought severity involves the analysis of the duration and deficit in dry conditions.

2.3 Hydrological modeling with artificial neural network (ANN) model

The rainfall and runoff data were simulated through the artificial neural network method using several inputs and outputs [9, 10]. This method imitates the function of the human nervous system and works in line with human learning patterns [11]. The backpropagation using a binary sigmoid activation function has been discovered to be a good ANN model for hydrology. This is due to the fact that the activation function is the net (network) of the linear combination of the inputs and their weights. It is important to note that this model fits the pattern characteristics of the input and output data which are required to follow a normal distribution in the range of 0–1 (0,10,9) [12, 13].

The determination of the hydrological data requires that the parameters recorded for January be influenced by the hydrological conditions in January of the previous year and the same trend is expected to continue from February to December. The data were also sorted from January to December at the end of the time series and the estimation was continued for the next few years. Moreover, the runoff model was designed in line with the rainfall data input based on the analysis of variables in the HDI analysis. Meanwhile, the architecture and equations of the backpropagation method used in this research are presented in Figure 1 [14, 15].

Where P₁ is the 1st data input, P_n is the nth data input, Z_1.1 is the 1st auxiliary variable in hidden layer 1, Z_1.2 is the 2nd auxiliary variable in hidden layer 1, Z_2.1 is the 1st auxiliary variable in hidden layer 2, Z_2.2 is the 2nd auxiliary variable in hidden layer 2, b (=1) is the specified bias value which is equal to one, and Qn is the nth data output.

Figure 1.

ANN for discharge calculation.

$$z\_\mathit{ne}{t}_j=v_{j0}+\sum _{i=1}^n{P}_i{v}_{\mathit{ji}}$$E2

$$z_j=f\left(z\_\mathit{ne}{t}_j\right)=\frac{1}{1+e^{-z\_\mathit{ne}{t}_j}}$$E3

$$Q\_\mathit{ne}{t}_k=w_{\mathit{ko}}+\sum _{j=1}^p{z}_j{w}_{\mathit{kj}}$$E4

$$Q_k=f\left(Q\_\mathit{ne}{t}_k\right)=\frac{1}{1+e^{-Q\_\mathit{ne}{t}_k}}$$E5

2.4 Discharge prediction model analysis

The results from the model analysis were examined based on the theory of runs which involves the analysis of two processes of sequential and opposite events at a boundary which is within a certain period [16].

The two kinds of “runs” available are run-length and run-sum with the run-length usually used for hydrological analysis related to the prediction of time which is also known as the frequency analysis while the run-sum is for those related to the prediction of rainfall duration and intensity within the analysis season. The overview of these runs is presented in the following Figure 2.

Figure 2.

Overview of positive run-length, m, positive run-sum, S, negative run-length, n, and negative run-sum, D, on a discrete series [16].

The run test was conducted on the simulated data, and since m is a positive run-length while n is a negative run-length as indicated in Figure 2, the total run-length (r) is m + n. Moreover, the estimated value of r, E(r), was expressed using the probability (q) for the estimated positive run-length, m, and (p) as the estimated probability value for the negative run-length, n. This relationship was expressed as follows [16]:

$$E\left(r_q\right)=\frac{1}{q\left(1-p\right)}=\frac{1}{\mathit{pq}}$$E6

with boundary conditions 0 < q < 1

$$p=\frac{m}{m+n};=\frac{n}{m+n};{\overline{r}}_q=\frac{1}{k_r}{\sum }_{j=1}^{\mathit{kr}}{r}_{q,j}$$E7

where ${\overline{r}}_q$ is the total run-length by q (probability), j = 1, 2, 3, …, kr; and kr is the total number of run-length.

It is also important to note that the estimated output in the model analysis is expected to be within the tolerance limit formulated as follows:

$$\frac{1}{\mathit{pq}}-\frac{t_{\alpha /2}}{\mathit{pq}}{\left(\frac{p^3+q^3}{k_r}\right)}^{1/2}\le {\overline{r}}_q\le \frac{1}{\mathit{pq}}+\frac{t_{\alpha /2}}{\mathit{pq}}{\left(\frac{p^3+q^3}{k_r}\right)}^{1/2}$$E8

where α is the tolerance limit (5%) and t is the normal distribution value from the “t-table.”

2.5 Drought index parameter reliability analysis

The parameter reliability analysis was conducted because the model describes the sample to be generalized to the population and this led to the application of a sample-based model suitability test. In mathematical analysis, reliability is the ratio of the total items to the total variance and was applied to the dry duration and deficit variables of the applied watershed. The general formula used is the Cronbach’s alpha equation as follows [17]:

$${\alpha }_r=\frac{k}{k-1}\left[1-\frac{\Sigma {{\sigma }_i}^2}{{{\sigma }_i}^2+2\left(\Sigma {\sigma }_{\mathit{ij}}\right)}\right]$$E9

Where $=\frac{1}{10}{i}^2$ is the number of variances i (the number of diagonals), $\frac{1}{10}\mathit{ij}$ is the covariance of items i and j, ${{\sigma }_i}^2+2\left(\Sigma {\sigma }_{\mathit{ij}}\right)$ is the total variance, and α_r is the reliability of the model.

3. Research method

The rationale of the research concept was used with the drought processes observed to have occurred sequentially starting from meteorological drought, agronomical drought, hydrological, to socioeconomical droughts, as indicated in Figure 3 [18, 19, 20]. The concept was further used to create an operational research framework that starts from the preparation of the model through the collection and processing of rainfall, discharge, and climate data of the research location followed by the modeling and verification as well as the calculation of HDI. The model was later applied and tested on the selected watershed as presented in Figure 4. The model is built based on historical data in the observation Catchment Area. Then applied to the simulation Catchment Area.

Figure 3.

Concept of research.

Figure 4.

Research framework.

4. Results and discussion

4.1 Data consistency and normality

The data from the rainfall stations in the selected watershed area, including Kedungsangku, Senduro, and Ranupakis stations, were tested for consistency and the results showed that they are consistent. The same process was conducted on the rainfall data from the applied watershed, such as Dampit, Poncokusumo, Sengguruh, Tangkil, and Wagir stations, and the results also showed that the data are consistent.

The normality test was conducted on the discharge data using pdf analysis with the indicator for the normal distribution being the formation of a straight line when the data are plotted visually. The results showed that the discharge data at Automatic Water Level Record (AWLR) and Umbul Dam (observation CA) are not normally distributed and the same trend was also observed for the data from Sutami Dam (Simulation CA).

4.2 Data statistical analysis and threshold

The distribution analysis showed that the discharge data from AWLR and Umbul Dams follow the gamma distribution and the same was also reported for the rainfall data from Kali Asem Catchment Area and Umbul sub-watersheds as well as the discharge data from Sutami Dam.

The next analysis conducted based on the gamma distribution showed the following results for the Kali Asem and Umbul sub-watersheds:

• Qaw50 = 16.27(m³/sec) which is interpreted as normal discharge in AWLR.

• Qd50 = 21.31(m³/sec) which is interpreted as normal discharge in Umbul Dam.

• Qaw80 = 11.31(m³/sec) which is interpreted as dependable discharge in AWLR.

• Qd80 = 14.32(m³/sec) which is interpreted as dependable discharge in Umbul Dam.

• P50 = 36.84(mm) which means normal rainfall.

The same method was used to obtain the threshold values for Sutami Dam as follows:

• P50 = 30.68(mm) which is interpreted as normal rainfall.

• Q50 = 64.34(m³/sec) which is interpreted as normal discharge in Sutami Dam.

• Q80 = 40.03(m³/sec) which is interpreted as normal discharge in Sutami Dam.

• WL50 = +267.91(m) which is interpreted as the dependable elevation of the water level in Sutami Dam.

• WL80 = +264.29(m) which is interpreted as the dependable elevation of the water level in Sutami Dam.

4.3 Deficit and dry duration observational data

The deficit is the difference between the volume of water shortage and the threshold, while duration is defined as the total time the deficit occurred. The results of the analysis showed that the duration and deficit do not have the same visual pattern and this means the deficit value does not describe the duration value and vice versa.

The debit deficit on the Very Dry criteria (based on Qaw80), consecutive (run-length) (d), and the total deficit (run-sum) was recorded from January 3 to October 3, 1991, with a total deficit of 142,928 (m³/sec. Ten daily) and 27 ten daily followed by May 1 to Dec 1, 1997 with 15.88 (m³/sec. Ten daily) and 22 ten daily, and Jan 1 to Aug 2, 2003 with 46,725 (m³/sec. Ten daily) and 25 ten daily, respectively, as indicated in Figures 5 and 6.

Figure 5.

Deficit based on Qaw50 and Qaw80 in the AWLR of the Kali Asem sub-watershed [21].

Figure 6.

Deficit of discharge, Qaw in 1991, 1997, 2003 [21].

4.4 HDI₁ and HDI₂ in Kali Asem and Umbul sub-watersheds

The calculation showed that the HDI₁ for the Kali Asem sub-watershed at Qaw₈₀ was 0.00, the HDI₁ limit at Qaw₅₀ was 0.0180, and at 70% Qaw80 was recorded to be −0.0123. Moreover, the result for the HDI₁ in the Umbul sub-watershed at Qd80 was found to be 0.00, the HDI₁ limits at Qd50 were 0.0200 while at 70% Qd80 was −0.0123.

The discharge prediction simulation model implemented the backpropagation ANN using the following parameters:

1. The model consists of one input, two hidden, and two output layers.

2. The network was formed using Descant Gradient Learning (trainingdm) with logsig used for activation in the hidden layer and purelin in the output layer.

3. The model simulation stops at the specified epoh of 1000 epoh or a mean square error (MSE) of 0.05.

The “scatterplot” of the Q_simulated and Q_observed in Figure 7 shows a straight line and coincides. This means the simulation data has a statistical character that is not significantly different from the observation data.

Figure 7.

Suitability analysis for the simulation results based on the discharge at Umbul dam [21].

The HDI₁ value also ranged from 0.002 to 0.024 while the HDI₂ value ranged from 0.001 to 0.022 (Extremely Dry). Moreover, a shift was observed in the dry time and this means it is difficult for the simulation to accurately predict the start or end time for the drought condition (see Figure 8).

Figure 8.

FDC analysis on the Q observation and simulation (Q_observation and Q_simulation, as well as IK observation and IK simulation [21].

The calculation showed that the equation can be written in the form of a matrix as follows:

• $\text{Input weight},\left(\mathrm{P}\to {\mathrm{Z}}_1\right)=\left[\begin{array}{cc}-\mathrm{2,434}& \mathrm{2,385}\\ -\mathrm{1,422}& -\mathrm{1,312}\\ -\mathrm{2,730}& -\mathrm{1,428}\\ \mathrm{2,198}& -\mathrm{2,700}\\ \mathrm{1,798}& -\mathrm{0,613}\\ -\mathrm{0,823}& \mathrm{3,194}\\ \mathrm{1,206}& -\mathrm{0,166}\\ -\mathrm{4,729}& -\mathrm{3,937}\end{array}\right]$

• $\text{Input bias weight},\mathrm{b}=\left[\mathrm{5,893}\mathrm{5,369}\right]$

• $\text{Hidden layer weight}-1,\left({\mathrm{Z}}_1\to {\mathrm{Z}}_2\right)=\left[\begin{array}{cc}\mathrm{0,755}& \mathrm{7,685}\\ -\mathrm{7,359}& \mathrm{2,927}\end{array}\right]$

• $\text{Hidden layer bias weight}-1,\mathrm{b}=\left[\begin{array}{cc}-\mathrm{8,477}& -\mathrm{1,755}\end{array}\right]$

• $\text{Hidden layer weight}-2\left({\mathrm{Z}}_2\to \mathrm{Q}\right)=\left[\begin{array}{cccccccc}\mathrm{0,021}& -\mathrm{0,033}& \mathrm{0,891}& \mathrm{0,762}& -\mathrm{0,087}& \mathrm{0,068}& \mathrm{0,748}& -\mathrm{0,101}\\ -\mathrm{0,022}& \mathrm{0,562}& -\mathrm{0,739}& \mathrm{0,987}& \mathrm{0,999}& -\mathrm{0,648}& -\mathrm{0,179}& \mathrm{0,332}\end{array}\right]$

• $\text{Hidden layer bias weight}-2\mathrm{b}=\left[\begin{array}{cccccccc}\mathrm{0,713}& \mathrm{0,129}& \mathrm{0,236}& \mathrm{0,267}& \mathrm{0,681}& \mathrm{0,644}& \mathrm{0,389}& \mathrm{0,383}\end{array}\right]$

The model was simulated at Umbul Dam for verification and the results showed that the simulation data is within the tolerance limit (α) of 5% which indicates that it is not significantly different. A run test was also conducted on the drought parameter.

• The statistical tests conducted using analysis of variance showed that the statistical value in Levene’s test for the data that was not necessarily normally distributed was 1.38 with p = 0.123. The p-value α > (0.05) means the null hypothesis was not rejected and this indicates the same variance.

• The Bartlett test also showed that p = 0.995 > α = 0.05 and this means the simulation data are not the same as the observation data while the median analysis conducted using Levene’s test also obtained p = 0.123 > α (0.05) and this means the null hypothesis was not rejected (same median). Moreover, the run analysis showed that E(r) Qdsim4 was 4.049 and the E(r)Qdsim was 4.016, which lies within the tolerance limit of 4.701–3.331 at α of 5%. This means the simulation data are not the same as the observation data.

4.5 HDI₁ in Sutami watershed

The threshold analysis for the Sutami sub-watershed using observational data from 1991 to 2006 showed the value for Q50 was 64.34 m³/sec and Q80 was 40.03 m³/sec. Meanwhile, HDI₁ for Q80 was calculated to be 0.00 while the limit at Q50 was 0.0119 and at 70% Q80 was −0.0059.

The calculation showed that HDI₁ ranged from 0.000 to 0.007 while HDI₂ was from 0.000 to 0.007 in line with the severe drought criteria and this means the simulation was unable to determine the exact time the drought starts and ends. Moreover, the dry deficit duration calculation showed that the simulation data was accepted at α = 5% while the median test indicated that the simulation data have the same characteristics and are not significantly different from the observation data. The “scatterplot” test of the simulation data presented in Figure 9 and the Flow Duration Curve test in Figure 10 also showed that there was a match between the observation and simulation data.

Figure 9.

The suitability analysis of the simulated scatterplot in Sutami dam [21].

Figure 10.

Flow duration curve conformity analysis of the simulation results at Sutami dam [21].

4.6 HDI₂ at Sutami sub-watershed

The Z test analysis showed that p = 0.064 > α = 0.05 and this means the simulation and observation are not the same as the mean test while Mann–Whitney median test also showed the value of p = 0.0547 > α = 0.05 and this also indicates the same trend as indicated in Table 1.

Variable	N	Mean	StDev	SE-Mean	95% CI		Z	P
Q	288	71,676	35,894	2,115	(67,531;	75,821)	0	0,999
Qsim	288	75,603	33,244	2,115	(71,458;	79,748)	1,86	0,064

Table 1.

Results of the Z-test Q_Observation and Q_Simulation in Sutami dam [21].

The Q_simulation mean value was observed to be within the tolerance range of the Q_observation while the p-value > α = 0.05 indicates the simulation is not the same as the observation data but the analysis conducted at Sutami Dam showed that the simulation data can be accepted.

The results showed that the simulation can be used to predict the drought in the Sutami sub-watershed for 2007–2014 with the dry limit observed to be at 0.000 < HDI < 0.0119, the very dry limit at −0.0059 < HDI < 0.000, and the extremely dry limit at HDI < −0.0059. The analysis further indicates that the lowest HDI value that has ever occurred was −0.006 which is included in the extremely dry criteria.

4.7 Reliability of the simulation results in Sutami sub-watershed

The reliability calculated using Eq. (9) showed that the analysis value for Q₅₀ was 79% and Q₈₀ was 81%.

5. Conclusion

The research concludes that:

1. The calculation model for the estimated discharge is based on the ten daily rainfall data from the selected watershed which includes Kali Asem Catchment Area and Umbul sub-watersheds. This is in the form of a mathematical equation which is stated as follows

   Q1 = −0,021 Z₂₁₋0,022 Z₂₂ + 0,713

   where:

   Z₂₁ = 0,755 Z₁₁ + 7685 Z₁₂–8477

   Z₂₂ = − 7359 Z₁₁ + 2927 Z₁₂–1755

   And,

   Z₁₁ = −2434 P1–1422 P2–2730 P3 + 2198 P4 + 1798 P5–0,823 P6 + 1206 P7–4729 P8 + 5893

   Z₁₂ = 2385 P1–1312 P2–1428 P3–2700 P4–0,613 P5 + 3194 P6–0,166 P7–3937 P8 + 5369

   Q2 = −0,033 Z₂₁ + 0,562 Z₂₂ + 0,129

   Q3 = 0,891 Z₂₁₋0,739 Z₂₂ + 0,236

   Q4 = 0,762 Z₂₁ + 0,987 Z₂₂ + 0,267

   Q5 = −0,087 Z₂₁ + 0,999 Z₂₂ + 0,681

   Q6 = 0,068 Z₂₁₋0,648 Z₂₂ + 0,644

   Q7 = 0,748 Z₂₁₋0,179 Z₂₂ + 0,389

   Q8 = −0,101 Z₂₁ + 0,332 Z₂₂ + 0,383

   The simulation test showed that the model was accepted at α = 5% or 95% confidence level and the reliability of its parameters, when applied to Sutami sub-watershed, was found to be only 80%.:

2. The dry duration and deficit were determined based on the threshold value obtained from historical discharge data using pdf and CDF analysis. The criteria used were that the watershed is believed to be dry to wet when the discharge is above the threshold and dry when it is below the threshold. Moreover, the analysis showed that the data were gamma distributed with a threshold value of:

   • P50 = 36.84 which means the normal rainfall in the analysis period was 36.84 mm/ten daily.

   • Qaw50 = 16.27 which means the normal discharge at AWLR during the analysis period was 16.27 m³/sec. Ten daily.

   • Qd50 = 21.31 which means the normal discharge at AWLR in the analysis period was 16.27 m³/sec. Ten daily.

   • Qaw80 = 11.31 which means that the dependable discharge in AWLR in the analysis period is 11.31 m³/sec. Ten daily.

   • Qd80 = 14.32 which means the dependable discharge at Umbul Dam during the analysis period was 14.32 m³/sec. Ten daily.

3. HDI for Kali Asem and Umbul sub-watersheds have the following limitations:

   • Dry when 0.000 < HDI < 0.018,

   • Very Dry when −0.012 < HDI < 0.00,

   • Extremely Dry when HDI < −0.012.

     HDI in Sutami sub-watershed has the following limits:

   • Dry when 0,000 < HDI < 0,012,

   • Very Dry when −0,006 < HDI < 0,00,

   • Extremely Dry when HDI < −0,006.

4. The criteria for IK at Kali Asem and Umbul sub-watersheds based on the research results showed the extremely dry (ED) with duration sharpness in zone 2 which occurred in 1991 and very dry (VD) with a sharpness of zone 2 in 2003. Meanwhile, the lowest HDI value in 1991 was −0.024 and the duration was 23 ten daily. In 2003, the lowest HDI value was −0.016 with a mean value of −0.007 and a duration of 23 ten daily. This means the drought occurred twice during the analysis period, in 1991 and 2003, with a duration sharpness of more than 7 months which exceeds the usual 6 months for the dry season.