Hydrological drought forecasting under a changing environment in the Luanhe River basin

Forecasting the occurrence of hydrological drought according to a forecasting system is an important disaster reduction strategy. In this paper, a new drought prediction model adapted to changing environments was constructed. Taking the Luanhe River basin in China as an example, first, nonstationarity analysis of hydrological sequences in the basin was carried out. Then, conditional distribution models with the human activity factor as an exogenous variable were constructed to forecast hydrological drought based on meteorological drought, and the results were compared with the traditional normal distribution model and conditional distribution model. Finally, a scoring mechanism was applied to evaluate the performance of the three drought forecasting models. The results showed that the runoff series of the Luanhe River basin from 1961 to 2010 were nonstationary; moreover, when human activities were not considered, the hydrological drought class tended to be the same as the meteorological drought class. The calculation results of the models involving HI as an exogenous variable were significantly different from the models that did not consider human activities. When the current drought class tended towards less severe or normal, the meteorological drought tended to turn into more severe hydrological drought with the increase in human index values. According to the scores of the three drought forecasting models, the conditional distribution models involving the human index can further improve the forecasting accuracy of drought in the Luanhe River basin.

1 Introduction

Typically, meteorological drought is regarded as the beginning of a drought event; after the occurrence of meteorological drought, other drought phenomena occur, such as hydrological drought (Fendeková et al., 2018; Fuentes et al., 2022; Wang et al., 2021). However, there is a delay period from meteorological drought to hydrological drought (Ding et al., 2021; Xu et al., 2019; Onyutha, 2017; Cammalleri and Vogt, 2018). Therefore, the occurrence of hydrological drought can be forecasted according to meteorological drought monitoring. Accurate hydrological forecast information is beneficial to reduce the losses caused by hydrological drought (Behzad and Hamid, 2019; Oertel et al., 2018 Dixit et al., 2022; Jehanzaib et al., 2020).

To identify the drought characteristics of the region, scholars have developed drought indices. For example, the standardized precipitation index (SPI) is typically used to identify and capture the characteristics of meteorological drought (McKee et al., 1993). Considering the influence of precipitation and temperature, Vicente-Serrano et al. (2010) proposed the standardized precipitation evapotranspiration index (SPEI) to characterize meteorological drought. The standardized runoff index (SRI), which focuses on the surface runoff of catchments, is typically used to indicate hydrological drought (Shukla and Wood, 2008). Aghelpour and Varshavian (2021) proposed the multivariate standardized precipitation index (MSPI) to forecast hydrological drought in Iran.

Statistical technology is an effective prediction method that has been widely used in drought forecasting in recent years (Alquraish et al., 2021; Abbasi et al., 2021; Dehghani et al., 2019; Mallya et al., 2013). For instance, neural network models have been proposed to combine multiple data for drought prediction (Rezaeianzadeh et al., 2016; Mokhtarzad et al., 2017; Ahnadi et al., 2011), and time series models can be used to analyse the variation in time series such as rainfall and runoff to achieve drought prediction (Moghimi et al., 2020; Natsagdorj et al., 2021; Stojković et al., 2020). The conditional probability model was proposed by Cancelliere et al. (2007) and developed for drought forecasting by Bonaccorso et al. (2015). Bonaccorso et al. (2015) showed that the conditional probability model can calculate the transition probabilities from the current drought index values to the future drought classes, and this is a more robust method that can be used to forecast drought than the traditional probability prediction models (such as the multivariate normal distribution model and Markov method).

A change in the environment may lead to the nonstationarity of the relationship between hydrological series (for example, precipitation and runoff series), which also occurs in the Luanhe River basin (Wang et al., 2016, 2018, 2020; Li et al., 2015). Traditional drought prediction methods need to be further improved to adapt to nonstationary conditions (Wang et al., 2022; Zhao et al., 2018; Chen et al., 2021). Ren et al. (2017) found that the conditional distribution model using large-scale climatic indices as covariates can improve the accuracy of meteorological drought forecasting in the Luanhe River basin. Although some progress has been made in the study of drought forecasting, there are relatively few studies considering the impact of the changing environment.

In this paper, to analyse the impact of human activities on hydrological drought, we constructed the human activity index (HI) based on the restoration method. Subsequently, conditional distribution models with the HI as the exogenous variable were developed to forecast hydrological drought based on meteorological drought, and then the results were compared with the traditional normal distribution model and conditional distribution model; as a result, the impact of human activities on transition probabilities was illustrated. A scoring mechanism was applied to the evaluation of the three probability models.

In addition to the introduction, this paper contains the following sections. Section 2 introduces the study area and data. Section 3 briefly describes the methods used in the research. Section 4 introduces the model construction and calculation results and analyses the results. Section 5 presents the prospects.

2 Study area and data

The Luanhe River basin, located in the subtropical monsoon region, covers an area of approximately 33 700 km². Its geographical location is shown in Fig. 1. Due to the influence of geographical location and topography, the annual average north–south temperature difference in the basin is 11.5 ^∘C, and the annual rainfall distribution is uneven. Less rain in spring and winter makes the area prone to meteorological drought and hydrological drought, while there is relatively more rainfall in summer. The average rainfall in summer is approximately 200–560 mm, resulting in highly variable annual runoff in the basin. The concentrated rainfall in summer has also become one of the remarkable features of the climate in this area. In recent years, the precipitation and inflow of the Luanhe River basin have gradually decreased, the water level of the Panjiakou Reservoir in the lower reaches of the basin has decreased, the runoff has decreased, and the frequency of meteorological drought and hydrological drought has significantly increased. Especially after entering the 21st century, the river basin has exhibited continuous drought and even extreme drought. With the change in the global climate and the impact of human activities on the basin environment, drought disasters in the Luanhe River basin occur frequently, causing significant social and economic losses.

Figure 1 The geographical location of the Luanhe River basin.

Influenced by topography, meteorology, hydrology and hydrogeological conditions, the spatial distribution of groundwater resources in the Luanhe River basin is quite different. The recharge and storage conditions of shallow groundwater in plain areas and intermountain basins are relatively superior, and the content of groundwater in mountainous areas is relatively small (the area of mountainous areas in the Luanhe River basin accounts for 98.2 %). Therefore, the total amount of water resources in the Luanhe River basin is mainly considered to be affected by the amount of surface water resources.

In this paper, the monthly rainfall data from 26 stations in the Luanhe River basin from 1961 to 2010 were provided by the Hebei Provincial Hydrology and Water Resources Investigation Bureau. The average monthly rainfall data of the area were obtained by the inverse distance weighting interpolation method. The runoff data from 1961 to 2010 came from the inflow runoff series of the Panjiakou Reservoir. The SPI and SRI can be calculated for 1-month, 3-month, 6-month, and 12-month timescales to characterize meteorological drought and hydrological drought based on these data.

3 Methods

3.1 Nonstationarity test method

In the case of environmental changes, nonstationarity may occur in hydrological series. The Pettitt test, as one of the important methods to test whether there is nonstationarity in time series, can identify whether there are change points in the sample series (Malede et al., 2022). Assuming that the sample sequence is $x=(x_{\mathrm{1}},x_{\mathrm{2}},\mathrm{\dots }{x}_n)$, the formula is as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }{U}_{t,n}=U_{t-\mathrm{1},n}+\sum _{i=\mathrm{1}}^n\text{sgn}(x_t-x_i)(t=\mathrm{2},\mathrm{3},\mathrm{\dots }n)t_{\mathrm{0}}\text{(2)}& {\displaystyle }\text{sgn}(x_t-x_i)=\left\{\begin{array}{rc}\mathrm{1}& x_t-x_i>\mathrm{0}\\ \mathrm{0}& x_t-x_i=\mathrm{0}\\ -\mathrm{1}& x_t-x_i=\mathrm{0}\end{array}\right.,\end{array}$$

where U_t,n is the test statistic, which indicates the cumulative number of the values at time t greater than or less than the values at time i. In addition, if $K_{t_{\mathrm{0}},n}$ satisfies the following,

$$\begin{array}{}\text{(3)}& K_{t_{\mathrm{0}},n}=max\left|U_{t,n}\right|(t=\mathrm{1},\mathrm{2},\mathrm{\dots },n),\end{array}$$

then t₀ is considered to be the change point, and the cumulative probability of possible change is determined by $K_{t_{\mathrm{0}},n}$:

$$\begin{array}{}\text{(4)}& P_{t_{\mathrm{0}},n}=\mathrm{2}\exp \left(-{\displaystyle \frac{\mathrm{6}{K}_{t_{\mathrm{0}},n}^{\mathrm{2}}}{n^{\mathrm{3}}+n^{\mathrm{2}}}}\right).\end{array}$$

Given the significance level α=0.05, if $P_{t_{\mathrm{0}},n}>\mathrm{0.95}$, it means that the point is a significant change point (Li et al., 2022; Koudahe et al., 2018). Furthermore, combined with the Mann–Kendall test, the trend characteristics of the sample series can be obtained (Li et al., 2018).

The sliding T test is a basic method commonly used in statistics. According to the mean and variance of the two sample sequences before and after the change points in the runoff time series, the two sample sequences are tested (Li et al., 2020):

$$\begin{array}{}\text{(5)}& {\displaystyle }t={\displaystyle \frac{{\stackrel{\mathrm{‾}}{x}}_{\mathrm{1}}-{\stackrel{\mathrm{‾}}{x}}_{\mathrm{2}}}{S\sqrt{\frac{\mathrm{1}}{n_{\mathrm{1}}}+\frac{\mathrm{1}}{n_{\mathrm{2}}}}}}\text{(6)}& {\displaystyle }S=\sqrt{{\displaystyle \frac{(n_{\mathrm{1}}-\mathrm{1})S_{\mathrm{1}}^{\mathrm{2}}+(n_{\mathrm{2}}-\mathrm{1})S_{\mathrm{2}}^{\mathrm{2}}}{n_{\mathrm{1}}+n_{\mathrm{2}}-\mathrm{2}}}}\text{(7)}& {\displaystyle }{S}_{\mathrm{1}}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{n_{\mathrm{1}}-\mathrm{1}}}{\sum }_{t=\mathrm{1}}^{n_{\mathrm{1}}}(x_t-{\stackrel{\mathrm{‾}}{x}}_{\mathrm{1}}{)}^{\mathrm{2}}\text{(8)}& {\displaystyle }{S}_{\mathrm{2}}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{n_{\mathrm{2}}-\mathrm{1}}}{\sum }_{t=\mathrm{1}}^{n_{\mathrm{1}}+n_{\mathrm{2}}}(x_t-{\stackrel{\mathrm{‾}}{x}}_{\mathrm{2}}{)}^{\mathrm{2}},\end{array}$$

where the change point is x_t, n₁ and n₂ represent the sample size before and after the change point, and $S_{\mathrm{1}}^{\mathrm{2}}$ and $S_{\mathrm{2}}^{\mathrm{2}}$ represent the variance of the samples before and after the change point, respectively. If the t statistic satisfies t>t_α at the significance level of α=0.05, the point can be considered a change point.

The Spearman correlation test can be applied to test the trend of time series, and the specific description refers to the article of Bishara and Hittner (2012).

3.2 Human activity index

The rainfall and runoff series of the watershed are usually strongly correlated. However, under the interference of human activities, the relationship between rainfall and runoff changes.

The double cumulative curve method can test the nonstationarity of the bivariate correlation between rainfall series and runoff series, and the point where the underlying surface is significantly altered by human activities can be determined according to the position of the slope change of the curve. Due to the short data series before and after the change point (20 years before the change point and 30 years after the change point), a linear equation was used to fit the relationship between precipitation and runoff.

The linear regression relationship of the cumulative rainfall and runoff series can be calculated according to the following formula:

$$\begin{array}{}\text{(9)}& \sum x=k\sum y+b,\end{array}$$

where x is the runoff series; y is the rainfall series; k is the correlation coefficient of the regression equation; and b is the intercept of the regression equation.

Human activities are the main reason for the nonstationarity of the runoff series in the watershed, so the HI can be constructed to quantify the impact of human activities on runoff. Based on the linear regression relationship established between the accumulated precipitation and the accumulated runoff before the change point, the theoretical runoff sequence during the human activity period can be calculated from the measured precipitation sequence. SRI^′ represents the standardized runoff index value without human activity interference, and SRI represents the normalized runoff index value calculated based on the measured runoff sequence under the disturbance of human activities. The HI is obtained by subtracting the theoretical SRI^′ and the actual SRI, and the calculation formula is as follows:

$$\begin{array}{}\text{(10)}& \text{HI}={\text{SRI}}^{\prime }-\text{SRI}.\end{array}$$

When HI>0, it can be assumed that human activities exacerbate hydrological drought; HI<0 means that the actual SRI is greater than the theoretical SRI without human activities, and when HI=0, the watershed is considered undisturbed by human activities.

3.3 Multivariate normal distribution model

The SPI is one of the important indicators for evaluating meteorological drought in the basin, and the SRI is an important indicator for evaluating hydrological drought in the basin. According to the rainfall data and runoff data in the basin, the SPI and SRI can be calculated at different timescales. Table 1 provides the drought class classification and corresponding SPI values and SRI values (Kolachian and Saghafian, 2021).

Table 1 Drought class classification and corresponding SPI values and SRI values.

As a traditional drought class forecasting model, the multivariate normal distribution model (Model 1) can forecast the future SRI class according to the current SPI class. Assuming that the current SPI and SRI series both satisfy a multivariable normal distribution, the joint probability density can be expressed as follows (Chang et al., 2022):

$$\begin{array}{}\text{(11)}& f_{Z_{v,\mathit{\lambda }}^{\left(k\right)},W_{v,\mathit{\lambda }+M}^{\left(k\right)}}(t,s)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }\left|\mathrm{\Sigma }\right|}}\cdot \exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{X}^T{\mathrm{\Sigma }}^{-\mathrm{1}}X\right),\end{array}$$

where Σ is the covariance matrix, and $X=[t,s{]}^T$. The form of the covariance matrix is as follows:

$$\begin{array}{}\text{(12)}& \mathrm{\Sigma }=\left[\begin{array}{cc}\mathrm{1}& \text{cov}\left[Z_{v,\mathit{\lambda }}^{\left(k\right)},W_{v,\mathit{\lambda }+M}^{\left(k\right)}\right]\\ \text{cov}\left[Z_{v,\mathit{\lambda }}^{\left(k\right)},W_{v,\mathit{\lambda }+M}^{\left(k\right)}\right]& \mathrm{1}\end{array}\right].\end{array}$$

Furthermore, according to the joint probability density function of the SPI value $Z_{v,\mathit{\lambda }}^{\left(k\right)}$ at year v and month λ and the future M month's SRI value $W_{v,\mathit{\lambda }+M}^{\left(k\right)}$, the analytical formula of the transition probability of the future SRI drought class can be obtained (Zhang et al., 2017):

$$\begin{array}{}\text{(13)}& P\left[W_{v,\mathit{\lambda }+M}^{\left(k\right)}\in C_M\right]={\displaystyle \frac{\int {\int }_{C_N{C}_M}{f}_{Z_{v,\mathit{\lambda }}^{\left(k\right)},W_{v,\mathit{\lambda }+M}^{\left(k\right)}}(t,s)\cdot \mathrm{d}t\cdot \mathrm{d}s}{{\int }_{C_N}{f}_{Z_{v,\mathit{\lambda }}^{\left(k\right)}}\left(t\right)\cdot \mathrm{d}t}},\end{array}$$

where C_M represents the drought class, and $f_{Z_{v,\mathit{\lambda }}^{\left(k\right)}}\left(t\right)$ represents the marginal density function of $Z_{v,\mathit{\lambda }}^{\left(k\right)}$ in the current λ month.

3.4 Conditional distribution model

The conditional distribution model (Model 2) proposed by Bonaccorso et al. (2015) is described as follows: when one group of sample data X obeys a normal distribution and satisfies $X\sim N({\mathit{\mu }}_{\mathrm{1}},{\mathrm{\Sigma }}_{\mathrm{1}})$, and another group of sample data Y also obeys a normal distribution, namely, $Y\sim N({\mathit{\mu }}_{\mathrm{2}},{\mathrm{\Sigma }}_{\mathrm{2}})$, then the total sequence can be written as follows:

$$\begin{array}{}\text{(14)}& B=\left[\begin{array}{c}X\\ Y\end{array}\right]\begin{array}{c}r\\ p-r\end{array}\sim N_p\left(\left[\begin{array}{c}{\mathit{\mu }}_{\mathrm{1}}\\ {\mathit{\mu }}_{\mathrm{2}}\end{array}\right],\left[\begin{array}{cc}{\mathrm{\Sigma }}_{\mathrm{11}}& {\mathrm{\Sigma }}_{\mathrm{12}}\\ {\mathrm{\Sigma }}_{\mathrm{21}}& {\mathrm{\Sigma }}_{\mathrm{22}}\end{array}\right]\right).\end{array}$$

When sequence Y obeys a normal distribution, the distribution of sequence X under the Y condition still satisfies a normal distribution; namely, the distribution of (X|Y) is as follows (Gong et al., 2021):

$$\begin{array}{}\text{(15)}& \left(X\mathrm{|}Y\right)\sim N({\mathit{\mu }}_{\mathrm{3}},{\mathrm{\Sigma }}_{\mathrm{3}}),\end{array}$$

where μ₃ represents the expected value under the conditional distribution, and Σ₃ is the conditional covariance matrix:

$$\begin{array}{}\text{(16)}& {\displaystyle }{\mathit{\mu }}_{\mathrm{3}}={\mathit{\mu }}_{\mathrm{1}}+{\mathrm{\Sigma }}_{\mathrm{12}}{\mathrm{\Sigma }}_{\mathrm{22}}^{-\mathrm{1}}(y-{\mathit{\mu }}_{\mathrm{2}})\text{(17)}& {\displaystyle }{\mathrm{\Sigma }}_{\mathrm{3}}={\mathrm{\Sigma }}_{\mathrm{11}}-{\mathrm{\Sigma }}_{\mathrm{12}}{\mathrm{\Sigma }}_{\mathrm{22}}^{-\mathrm{1}}{\mathrm{\Sigma }}_{\mathrm{21}}.\end{array}$$

Then, the probability of the current SPI value transitioning to the future SRI drought class can be deduced as follows (Ren et al., 2017):

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& P\left[W_{v,\mathit{\lambda }+M}\in C_M\mathrm{|}{Z}_{v,\mathit{\lambda }}=z_{\mathrm{0}}\right]={\int }_{C_{M_i}}^{C_{M_s}}\\ & \cdot {\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{\sigma }}_Z}}{e}^{-\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{x-\mathit{\rho }{z}_{\mathrm{0}}}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}\right)}^{\mathrm{2}}}\mathrm{d}x,\end{array}\end{array}$$

where Z_v,λ represents the SPI value of the current month λ, $W_{v,\mathit{\lambda }+M}$ represents the SRI value of the λ+M month, $C_{M_s}$ and $C_{M_i}$ are the upper and lower limits of drought class C_M, and the correlation coefficient between the current SPI value and the future SRI value is ρ. Furthermore, the current SPI and future SRI can be expressed as the standard normal cumulative distribution function Φ:

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& P\left[W_{v,\mathit{\lambda }+M}\in C_M\mathrm{|}{Z}_{v,\mathit{\lambda }}=z_{\mathrm{0}}\right]=\mathrm{\Phi }\left[{\displaystyle \frac{C_{M_s}-\mathit{\rho }\cdot z_{\mathrm{0}}}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}\right]\\ & -\mathrm{\Phi }\left[{\displaystyle \frac{C_{M_i}-\mathit{\rho }\cdot z_{\mathrm{0}}}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}\right].\end{array}\end{array}$$

The calculation of the correlation coefficient ρ is as follows:

$$\begin{array}{}\text{(20)}& \mathit{\rho }={\displaystyle \frac{\text{cov}\left[Z_{v,\mathit{\lambda }}^{\left(k\right)},W_{v,\mathit{\lambda }+M}^{\left(k\right)}\right]}{\sqrt{\text{var}\left(Z_{v,\mathit{\lambda }}^{\left(k\right)}\right)\text{var}\left(W_{v,\mathit{\lambda }+M}^{\left(k\right)}\right)}}},\end{array}$$

where K represents the timescale of the drought index. Assuming that the cumulative rainfall Y and runoff X satisfy a normal distribution, after the standardization process, the SPI value $Z_{v,\mathit{\lambda }}^{\left(k\right)}$ corresponding to cumulative rainfall Y and SRI value $W_{v,\mathit{\lambda }+M}$ corresponding to runoff X obey the standard normal distribution:

$$\begin{array}{}\text{(21)}& \text{var}\left(Z_{v,\mathit{\lambda }}^{\left(k\right)}\right)=\text{var}\left(W_{v,\mathit{\lambda }+M}^{\left(k\right)}\right)=\mathrm{1}.\end{array}$$

$\text{cov}\left[Z_{v,\mathit{\lambda }}^{\left(k\right)},W_{v,\mathit{\lambda }+M}^{\left(k\right)}\right]$ represents the covariance between the current SPI and the Sri value with a forecast period of M months. The calculation is as follows:

$$\begin{array}{}\text{(22)}& \begin{array}{rl}& \text{cov}\left[Z_{v,\mathit{\lambda }}^{\left(k\right)},W_{v,\mathit{\lambda }+M}^{\left(k\right)}\right]={\displaystyle \frac{\mathrm{1}}{\sqrt{\sum _{i=\mathrm{0}}^{k-\mathrm{1}}{\mathit{\sigma }}_{\mathit{\lambda }+M-i}^{\mathrm{2}}\sum _{j=\mathrm{0}}^{k-\mathrm{1}}{\mathit{\sigma }}_{\mathit{\lambda }-j}^{\mathrm{2}}}}}\\ & \cdot \sum _{i=\mathrm{0}}^{k-\mathrm{1}}\sum _{j=\mathrm{0}}^{k-\mathrm{1}}\text{cov}\left[X_{v,\mathit{\lambda }+M-j},Y_{v,\mathit{\lambda }-i}\right].\end{array}\end{array}$$

3.5 Conditional distribution model involving the HI as an exogenous variable

According to the above conditional probability model, when considering the HI as an exogenous variable, the model (Model 3) can be extended as follows:

$$\begin{array}{}\text{(23)}& {\displaystyle }\begin{array}{rl}& P\left[W_{v,\mathit{\lambda }+M}\in C_M\mathrm{|}{Z}_{v,\mathit{\lambda }}=z_{\mathrm{0}},H_{v,\mathit{\lambda }}=h_{\mathrm{0}}\right]\\ & ={\int }_{C_{M_i}}^{C_{M_s}}{\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{\sigma }}_z}}{e}^{-\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{x-{\mathit{\mu }}_z}{{\mathit{\sigma }}_z}\right)}^{\mathrm{2}}}\mathrm{d}x\end{array}\text{(24)}& {\displaystyle }{\mathit{\mu }}_z=E\left[W_{v,\mathit{\lambda }+M}\mathrm{|}{Z}_{v,\mathit{\lambda }},H_{v,\mathit{\lambda }}\right]={\mathrm{\Sigma }}_{\mathrm{12}}^{\prime }({\mathrm{\Sigma }}_{\mathrm{22}}^{\prime }{)}^{-\mathrm{1}}\left[\begin{array}{c}{z}_{\mathrm{0}}\\ h_{\mathrm{0}}\end{array}\right]\text{(25)}& {\displaystyle }{\mathit{\sigma }}_z^{\mathrm{2}}=\text{var}\left[W_{v,\mathit{\lambda }+M}\mathrm{|}{Z}_{v,\mathit{\lambda }},H_{v,\mathit{\lambda }}\right]=\mathrm{1}-{\mathrm{\Sigma }}_{\mathrm{12}}^{\prime }({\mathrm{\Sigma }}_{\mathrm{22}}^{\prime }{)}^{-\mathrm{1}}{\mathrm{\Sigma }}_{\mathrm{21}}^{\prime },\end{array}$$

where

$$\begin{array}{}\text{(26)}& {\displaystyle }{\mathrm{\Sigma }}_{\mathrm{12}}^{\prime }=\left[\text{cov}\left(W_{\mathit{\nu },\mathit{\lambda }+M},Z_{\mathit{\nu },\mathit{\lambda }}\right)\text{cov}\left(W_{\mathit{\nu },\mathit{\lambda }+M},H_{\mathit{\nu },\mathit{\lambda }}\right)\right]\text{(27)}& {\displaystyle }{\mathrm{\Sigma }}_{\mathrm{22}}^{\prime }=\left[\begin{array}{cc}\mathrm{1}& \text{cov}(Z_{\mathit{\nu },\mathit{\lambda }},H_{\mathit{\nu },\mathit{\lambda }})\\ \text{cov}(H_{\mathit{\nu },\mathit{\lambda }},Z_{\mathit{\nu },\mathit{\lambda }})& \mathrm{1}\end{array}\right]\text{(28)}& {\displaystyle }{\mathrm{\Sigma }}_{\mathrm{21}}^{\prime }=({\mathrm{\Sigma }}_{\mathrm{12}}^{\prime }{)}^T.\end{array}$$

3.6 Scoring mechanism

A scoring mechanism was applied to evaluate the performance of the drought forecasting models. In this method, the monthly drought transition probability is summed to evaluate the model (Chen et al., 2013), where p_s,t characterizes the transition probability in month t of year s, and n is the length of the validation period.

$$\begin{array}{}\text{(29)}& \text{Score}={\displaystyle \frac{\mathrm{1}}{\mathrm{12}n}}{\sum }_{t-\mathrm{1}}^{\mathrm{12}}{\sum }_{s=\mathrm{1}}^n{p}_{s,t}\end{array}$$

4 Results and discussion

4.1 Nonstationarity analysis

In this paper, the area average monthly rainfall data of the Luanhe River basin from 1961 to 2010 were obtained by spatial interpolation. The runoff data came from the inflow runoff series of the Panjiakou Reservoir. Given the significance level α=0.05, the nonstationarity test results are shown in Fig. 2.

Figure 2 The change points of the runoff series.

Figure 2a shows that the years of possible runoff change were 1979, 1996, 1997, 1998 and 1999. The P values in 1979 and 1998 were infinitely close to 1, which were considered to be extremely significant runoff change points. Among all the possible points satisfying t>t_α, there were two maximum points (Fig. 2b), namely, 1979 and 1998, which were considered to be possible runoff change points. The final change point needs to be judged based on the actual situation of the watershed.

The results of the Spearman correlation test (Table 2) indicate that the runoff series showed an upwards trend before 1979, but the trend was not significant. However, there was a significant downwards trend in the series after 1979. In general, the runoff series showed a significant downwards trend.

Table 2 Spearman correlation test results of runoff series trend.

In addition, according to historical records, local human activities (such as land use change and reservoir construction) are regarded as the main factors influencing runoff (Yan et al., 2018; Chen et al., 2021). Synthesizing the above analysis, 1979 was determined as the change point for the runoff sequence in the basin, and this conclusion was consistent with Li et al. (2015) and Wang et al. (2015).

4.2 Transition probabilities from current SPI values to future SRI classes

According to the normality test results of rainfall and runoff series, it was reasonable to apply the conditional distribution model. To analyse the influence of different timescales of the SPI on the transition probabilities, using the forecast period as 1 month and the SPI timescales at 1, 3, 6 and 12 months as examples, the probabilities of converting SPI values to SRI classes were calculated (Fig. 3).

Figure 3 Influence of the SPI timescale on transition probabilities (z₀: initial value of SPI).

As shown in Fig. 3, when meteorological drought is categorized as extreme drought, the probabilities of maintaining the SRI class in extreme drought increased with the increasing SPI timescale. While the SPI had a short timescale, the response of the future SRI class to rainfall was fast, so the hydrological drought was more likely to tend to a normal status. This situation also occurred when the current meteorological drought was in another status.

In addition, the transition probabilities of drought were distinct for different forecast periods. As seen in Fig. 4, when the forecast periods were short (M=1 or 2), the hydrological drought classes obtained from the transition of meteorological drought tended to be the same as those of meteorological drought. With the extension of the forecast period (M=2 or 3), the hydrological drought classes obtained from the transition tended to be lower than the meteorological drought or the normal status.

Figure 4 Influence of forecast period on transition probabilities (z₀: initial value of SPI).

4.3 Transition probabilities involving the HI as the covariate

The effects of human activities are complex. To quantify the impact of human activities, the change point was identified, and then it was believed that the difference in the relationship between precipitation and runoff before and after the change point was caused by human activities. Moreover, the HI is easy to calculate and can approximately replace the influence of human activities. According to historical records, local human activities (such as land use change and reservoir construction) were regarded as the main factors influencing runoff (Yan et al., 2018; Chen et al., 2021). According to the above nonstationarity test results, 1979 was the change point, and the linear regression relationship of the cumulative rainfall and runoff series before and after the change point was established. The calculation results are shown in Table 3.

Table 3 Linear regression relationship between cumulative precipitation (x per mm) and cumulative runoff (y per 10⁶ m³).

The HI results for different timescales are shown in Fig. 5.

Figure 5 Different average periods of the HI (HI-1: HI with a 1-month timescale; HI-3: HI with a 3-month timescale; HI-6: HI with a 6-month timescale; HI-12: HI with a 12-month timescale).

As shown in Fig. 5, the HI at all monthly scales generally ranged upwards, which means that human activities have intensified the occurrence of hydrological drought. According to historical statistics, many water conservancy projects were built in the basin from 1980 to 2000, and the construction and operation of large reservoirs in the mid-1990s may be the main reason for the serious negative values of the HI.

The HIs of different monthly scales were standardized, taking the 12-month timescale as an example, and the results were calculated as shown in Table 4.

Table 4 HI-12 monthly mean and standard deviation (SD).

Furthermore, the drought transition probabilities involving the HI can be calculated from Eq. (23). Using the forecast period of 1 month from December and the SPI timescale of 12 months as an example, the drought transition probabilities from the current SPI values to the future SRI classes were calculated (Fig. 6). To analyse the effect of human activities on the drought transition probability more clearly, the calculation results of the three models are compared here separately. The horizontal coordinate indicates the drought classes corresponding to the SRI for the coming month, and the vertical coordinate is the drought transition probability.

Figure 6 Drought transition probability under the influence of human activities (C₀ denotes the initial drought class of the SPI in the multivariate normal model; z₀ represents the initial value of the SPI in the conditional distribution model; Model 1: the normal distribution model; Model 2: the conditional distribution model; Model 3: the conditional distribution model involving the HI).

In Fig. 6a, when the initial z₀=0.75 and C₀=N, the results shown in Model 1 and Model 2 were similar, and the probabilities of the SPI values transitioning to the SRI classes in the future month in the normal class were close to 1. However, the results of Model 3 indicated that the probabilities of maintaining the SRI in the normal class in the future decreased as the HI increased. When HI =2, the future hydrological drought classes were more likely to transition to severe drought or extreme drought.

From the initial $z_{\mathrm{0}}=-\mathrm{1.25}$ and C₀= Mo (Fig. 6b), the results of Model 3 showed that the transition probabilities of the SPI values to a normal SRI class in the coming month were higher when the HI was less than 1. As the HI increased, the transition probabilities of the SPI values to a moderate drought or even a more severe drought in the future increased. In addition, the probabilities of maintaining moderate drought were the highest when human activities were not considered, and Model 2 showed a higher probability than Model 1.

While the initial meteorological drought class was severe drought (Fig. 6c), the probabilities of the future SRI drought class being in the normal class became larger as the HI decreased. When the effect of human activities was not considered, the probability that the current SPI value transitioned to the SRI class under severe drought in the future month was the highest, and the probability of being in the normal class was the lowest. For Model 2, the probability of the SRI classes transitioning to severe drought was higher than the result of Model 1.

It was noteworthy that when the initial $z_{\mathrm{0}}=-\mathrm{2.5}$ and C₀= Es (Fig. 6d), the probabilities of transition of the SPI values to the future SRI classes at the normal class were close to 1 as HI<0. However, hydrological drought was more likely to be moderate drought or severe drought, as the HIs were greater than 0, and the transition probabilities exceeded 0.25. For Model 1 and Model 2, the probabilities of transition of the current SPI values or classes to the future month SRI classes in extreme drought were both higher than 0.75, and Model 2 showed a higher probability than Model 1.

In general, for the evaluation of drought transition probabilities in the future month, hydrological drought classes tended to be the same as meteorological drought when human activities were not considered, and this situation was more significant in Model 2 than in Model 1. The calculation results of the model involving the HI as an exogenous variable were significantly different from those of the models that did not consider human activities. The calculation results of Model 1 and Model 2 showed that the future hydrological drought classes were more likely to be the same as the meteorological drought classes in the current period, and they were more significant in Model 2. In addition, it was obvious that the drought transition probabilities of Model 3 were significantly different from those of Model 1 and Model 2. Taking Fig. 6b as an example, when $z_{\mathrm{0}}=-\mathrm{1.25}$ and C₀= Mo, the result of Model 1 showed that the probability of the SPI values transitioning to the SRI classes in the future month in the normal class was close to 0.15, the result of Model 2 was close to 0, and the result of Model 3 (HI =0) was close to 0.95. The results of Model 3 (HI =0) indicated that hydrological drought was likely to remain at the normal class in the future month. Moreover, the value of the HI had a great impact on the results of Model 3; for example, when HI $=-\mathrm{2}$ or −1, the probabilities of the SPI values transitioning to the SRI classes in the future month in the normal class were both close to 1, but the probability was close to 0.65 and 0.17 when HI =1 and 2, respectively.The results further indicated that meteorological drought tended to turn into more severe hydrological drought with increasing HI values.

4.4 Model evaluation and analysis

To quantitatively evaluate the prediction accuracy of Model 1, Model 2 and Model 3, the study period was divided into a correction period (1961–2003) and a verification period (2004–2010), and then the drought transition probability from the SPI value or class to the SRI class in the future M-month was calculated. The calculation results are shown in Table 5.

Table 5 Model evaluation (Model 1: multivariate normal distribution model; Model 2: conditional distribution model; Model 3: conditional distribution model with the HI).

With the same timescale of the SPI, the model scores of Model 1 and Model 2 decreased as the forecast period M lengthened, while the model scores of Model 3 were not significantly affected by the forecast period M. Model 1 had the highest rating of 0.36 at an SPI of a 1-month timescale and a forecast period of 1 month; Model 2 reached the highest model rating of 0.74 at a 12-month timescale and a forecast period of 1 month; and Model 3 performed well at an SPI of 1-month timescale and a 12-month timescale. Overall, Model 3 had the highest rating, and Model 1 had the lowest rating for the same SPI timescale and the same forecast period, which also indicated that the forecast accuracy of the conditional distribution model considering the HI was higher for short-term forecasts with a forecast period of 3 months or less, and including the HI could further improve the forecast accuracy of the model.

5 Conclusions

Many studies have noted that human activities have a significant impact on watershed runoff in the Luanhe River basin. In this paper, three probability models were constructed to calculate the transition probabilities from the current SPI classes or values to the future SRI classes; then, a scoring mechanism was applied to evaluate the performance of the models.

The calculation results of Model 1 and Model 2 showed that the future hydrological drought classes were more likely to be the same as the meteorological drought classes in the current period, and they were more significant in Model 2. In addition, it was obvious that the drought transition probabilities of Model 3 were significantly different from those of Model 1 and Model 2. Under the condition of considering the HI, the results of the drought transition probability showed that when HI <0, the future hydrological drought classes tended to normal status, and this situation was more obvious with the decrease in the HI values, which indicates that human activities mitigate the degree of hydrological drought when HI <0. However, when HI >0, the future hydrological drought classes generally transitioned to more severe drought with increasing HI values. Thus, it was indicated that human activities exacerbate the degree of hydrological drought as HI >0.

Finally, a scoring mechanism was applied to the evaluation of the models, and the forecast results of the three models were evaluated. The results demonstrate that when the SPI timescale was the same, the scores of Model 1 and Model 2 decreased as the forecast period lengthened. In most cases, Model 2 performed better than Model 1, and the performance of Model 3 was the most stable of the three models and had the highest score. The conditional probability model considering the HI was more suitable for the Luanhe River basin, where human activities have a high influence.

Although this study has made some progress in the forecasting of hydrological drought in a changing environment, only the HI was considered as the exogenous variable in this paper, and human activities were generalized. In future studies, the HI can be analysed specifically. For example, the impact of land use and socioeconomics on drought prediction can be specifically analysed. In addition, climate factors can be further considered in future research.