Description of the basin

The Fuerte river watershed is located in the northwest of Mexico, in the states of Sinaloa, Sonora, Durango and Chihuahua (Figure 1), as part of the hydrological region Number 10, Sinaloa, between 25.68° and 28.24 ° N, and 106.12 ° to 109.43 ° W, with an area of 36 456 km². Altitudes in the watershed range from 3168 to -9 masl (^{INEGI, 2014}); its average annual precipitation is 691 mL and average temperature of 19.4 °C, according to the climatological normal values of 14 stations within the watershed, from 1961 to 2011 (^{SMN, 2014}).

The watershed is born in the Sierra Madre Occidental and flows into the north of Sinaloa in the Gulf of California. The Rarámuri indigenous ethnic group inhabits the upper part of the watershed; it is a zone of high marginalization that in the last decade had serious drought problems; the Valle del Fuerte is located in the middle and lower part, which is an important agricultural area of production of grains and vegetables for the national and international market. The Fuerte river watershed is important due to its surface extension, the volume of its runoff and the dams Miguel Hidalgo and Costilla, Luis Donaldo Colosio and Josefa Ortiz de Dominguez, which are used to generate electricity and are the main source of water for agricultural, urban and industrial use.

Climatological information

The climatological information used for the calculation of the SPI and SPEI drought indices was obtained from the national meteorological service database available online (^{SMN, 2014}). Besides, we used the monthly precipitation data (Pt), minimum temperature (Tmin) and maximum temperature (Tmax) of 14 weather stations during the period 1961-2011, which was the most complete and updated during this study.

The stations selected are distributed throughout the watershed; the most incomplete information was that of precipitation and on average 8 % of the rainfall missing data was generated (Table 1). In Batopilas, Chihuahua, one of the most extreme cases, 29 % of the rainfall data were generated. To generate the rain data, we used the inverse square distance method, and the details appear in an article for which the indices were generated (^{Castillo et al., 2017}).

Methodology

The meteorological, hydrological and agricultural droughts are represented by the drought index, which is a variable to assess the effect of a drought and includes intensity, duration, severity and spatial extension (^{Mishra and Singh, 2010}). The SPI (^{McKee et al., 1993}) is one of the most used for the detection and monitoring of droughts and is based on the probability of precipitation at any time scale to quantify the precipitation deficit (^{Velasco et al., 2004}).

The SPEI (^{Vicente et al., 2010}), in addition to the precipitation used by the SPI, considers the reference evapotranspiration (ET₀) based on the monthly water balance (Pt-ET₀), combining characteristics of the Palmer Drought Severity Index (PDSI) sensitivity and the simplicity of calculation, and the multi-temporal nature of SPI. It is very suitable to detect, monitor and explore the consequences of global warming in drought conditions (^{Vicente et al., 2010}).

The SPI and SPEI indices were developed for the 14 selected stations, at scales of 3, 6, 12 and 24 months. The SPI series were generated with the spi_sl_6.exe program developed by the National Drought Mitigation Center (^{NDMC, 2014}), and the SPEI was obtained with SPEI.R developed by ^{Beguería and Vicente (2014)} for the R/RStudio program.

In the SPEI calculation, the reference evapotranspiration method of Hargreaves-Samani was contemplated (^{Beguería et al., 2014}). For more information on obtaining SPI and SPEI, and their interpretation as indicators of drought, consult ^{Castillo et al. (2017)}, who detail the SPI and SPEI process in the watershed for our study.

In our study, the forecast of meteorological droughts was made by mixing autoregressive time series models with the Discrete Kalman filter (^{Kalman, 1960}), which is a set of mathematical equations to estimate the state of a process minimizing the mean of the quadratic error. The Kalman filter operates by means of a prediction and correction mechanism, the algorithm predicts the new state from a previous estimate, adding a correction term proportional to the prediction error, minimizing it statistically (^{Welch and Bishop, 2006}).

Autoregressive models are used in hydrology and meteorology because they are time dependent and easy to use (^{World Meteorological Organization, 2011}). In our study, we proposed the creation of two autoregressive models, one of second order (AR2) and the other of second order with exogenous input (ARX) to forecast the future monthly states of the SPI and SPEI drought indices based on the previous records of the series. The choice of autoregressive models will be shown in Table 2 of Results and Discussion.

The AR2 and ARX models (equations 1 and 2, respectively) relate the input of the system to its output by means of a linear equation in differences with constant coefficients:

where y _t is the observed value of the index at time t, which represents one month; r_t is the value of the exogenous variable (Pt, Tmin, Tmax, ET₀) at time t; e_t+1 is the error term in the index estimate; α_i and β_i are parameters; the indices na and nb specify the number of previous observations of the index and exogenous variables, respectively. From the general structure of both autoregressive models, the formulation in a state space allows to use them within the algorithm of the Discrete Kalman filter as follows:

where x _k+1 is the value of the index (not observed) of size (n x 1); A is the matrix of α_i size parameters (n x n); x _k is the index in time k of size (n x 1); B is the matrix of β_j exogenous size parameters (n x m); v_k is the vector that contains the exogenous variable registered for time k; z_k is the index in time k of size (m x 1); H is the transformation matrix that maps the state vector to the domain of the measurement with dimensions (m x n);W _k and v_k are vectors representing Gaussian noise in the process and noise in the measurement for each observation with sizes (m x 1), and it is expected that such Gaussian noises are distributed in a normal way with mean 0 and Q and R variances, respectively :

According to ^{Simon (2001)}, there is no correlation between Q and R, that is, they are independent random variables and can vary over time, but they are assumed constant for simplicity, and it is defined as:

where Q is the covariance matrix of the system disturbance; R is the covariance matrix of the measurement perturbation; wkT and vkT indicate the transpose of w_k and v_k , respectively, and E[.] means the expected value. The Q matrix is of dimension (na x na) and contains the adjusted value of the MSE (Square Mean Error) estimated by the autoregressive model in the calibration stage.

The matrix R is (1 x 1) and contains the expected error of the measurements, defined as a proportion (α) of the previous measurement, so R=α*(x_k -1), in this specific case α=0.001; higher values generate extreme data that are located outside the range in which values of drought indices normally present, reducing the prediction power of the model. With respect to the size of the matrices, the value of n is equal to na, and the value of m equal to nb; where na and nb indicate the number of previous observations of the index and the exogenous variables, respectively. The algorithm of the Discrete Kalman filter is described in Figure 2.

To create the AR and ARX models, we programmed routines in Matlab^® (^{Math Works 2015}), to obtain the values of the parameter matrices in line with the order of the autoregressive process according to equations 1 and 2. The ARX was also an autoregressive model of second order, but with exogenous variable.

Since the behavior of the drought in the watershed varies depending on the occurrence of the rainy season and phenomena such as El Niño or La Niña, among other global climatic phenomena, the parameters estimated for the AR and ARX models were recalculated after a certain time period (P) to obtain more representative forecasts of the previous conditions of humidity in each station. The forecast began in the period (t₀ + P) until the time (t₀ + 2P), in this point we recalculated the parameters of the models according to period [P: 2P], and so on until ending with the last data group. Thus, the implementation of DKF - AR2 and DKF - ARX became dynamic and incorporated the climatic changes that occurred in the watershed during the study period of each season or certain period.

We did the forecast of the SPI and SPEI indices at four-time scales for each station, and implemented the AR2 and ARX models (Pt, Tmin, Tmax and ET₀ as exogenous variables) with the Discrete Kalman filter. To know the effectiveness of the models in the forecast for L steps forward, that is, with 1, 2, 3 and 4 months in advance, we created a routine in Matlab^® to perform the forecast considering the information in time k and advance L steps in the forecast without updating the information (Figure 3); we conducted this process considering the P period of calibration of the AR2 and ARX models.

To evaluate the results of the forecast, we calculated the main statistics according to ^{Gupta et al. (2009)}; the RMSE (the root of the mean square error) and E (the Nash-Sutcliffe efficiency) were the criteria to be considered as they are the most used for the calibration and evaluation of hydrological models (^{Moriasi et al., 2007}). In addition, a 95% confidence interval was included that follows the route of each state of the series with the objective of establishing the degree of uncertainty associated with the forecast at each time step (L) with respect to the values observed.

The prediction interval (I.P.) was obtained according to the criteria established by ^{Chatfield (2004)}, the I.P. was calculated based on the following general form: 100(1-α)% I.P. for x_t+L is given by:

where etL=xt+L-x^t(L) denotes the forecast error at time t when forecasting with L advance periods; zα/2 denotes the percentage point of a standard normal distribution with a proportion of α/2 over it.

Forecast with DKF - ARX-Pt

In stations 8167, 25009, 25019, 25025, 25042, 25044, 25100 and 26053 (6 of them are in Sinaloa), we obtained better results in the SPEI forecast with the DKF-ARX-Pt model. The results showed that the use of precipitation as an exogenous variable improves the forecast compared to the one made with the DKF-AR2 model for these stations.

With the ARX model, precipitation, reference evapotranspiration, maximum and minimum temperature were tested as exogenous variables; however, in most cases precipitation contributed more information to the algorithm and generated better results in terms of E and R than the remaining variables (Table 6).

The association of the results with the climatic units of the watershed locates these stations in hot and dry climates. As already shown, the hot weather stations show series with greater oscillation between positive and negative values of the SPEI than the series of stations in temperate and semi-cold climates. By having less smoothed series the forecast becomes less precise, so the inclusion of an exogenous variable (precipitation) provides the necessary information to the model to improve its performance.

Figures 8 and 9 show the dispersion diagrams of SPEI of 12 and 24 months, respectively, for the stations where the ARX-Pt model was better. These Figures show the values of MSE, RMSE, E, R and PBE.

Forecast with DKF - AR2

The stations 8038, 8106, 8161, 8172, 8182 and 8267, all located in the state of Chihuahua on the Sierra Madre Oriental, obtained better results in the SPEI forecast with the DKF-AR2 model (Table 6). In these stations the inclusion of an exogenous variable such as precipitation, reference evapotranspiration, maximum and minimum temperature did not provide enough information to the algorithm to improve the results obtained with the AR2 model; if anything, it matched them.

The association of the results with the type of climate of the station locates these stations in semi-cold and temperate semi-humid climates.

The value of the PBE (Percentage Bias Error) statistic is negative in most stations, which indicates that the model DKF-AR2 underestimates the value of the SPEI index (Table 6).

In both forecast models of the SPEI, DKF-AR and DKF-ARX, the 12-month series show higher error and lower value of the Nash coefficient than the 24-month series; this is due to the fact that on a larger scale the time series assimilates more slowly the changes given in the balance Pt-ET₀ over time, and therefore presents more smoothness in its route, improving the forecast.

However, both forecasts of droughts of 12 and 24-month duration are useful. But given the characteristics of the drought indices (SPI and SPEI) at a certain scale (duration of drought), in terms of operability, the 12-month time scale is considered the most appropriate for the forecast of droughts. The latter is particularly suitable in the middle and lower parts of the watershed due to its agricultural and energy values, where annual policies of operation are analyzed when operating the dams, although the analysis of drought periods of 24 months is very useful as well.

In addition to the forecast of the drought index for the time t+1|t, we predicted the indices with 2, 3 and 4 months in advance; that is, given an observed value of the index in time t, we predicted the index values for the time t+1, t+2, t+3 and t+4 without updating the information, with the purpose of knowing the power of prediction of the models for different steps of advancement in time. However, as a sample, we present here only the forecast results for the 12-month SPEI. Figure 10 shows the forecast at the Miguel Hidalgo Dam station in Sinaloa, with the 12-month SPEI using the DKF-ARX-Pt model with 1, 2, 3 and 4 months in advance.

The forecast focuses on the point estimate of a 12-month SPEI value for some future month, but according to ^{Chatfield (2004)}: “the point forecast is adequate for many purposes, but a prediction interval is often very helpful to give a better indicator of future uncertainty”. Therefore, we calculated the confidence interval of the forecast at 95 % for the advance steps (L), as indicated in equation (9).

In Figure 10, in addition to the observed and predicted values, there is a stripe that follows the forecast path (the more negative the index, the greater the drought), and this is the confidence interval of the forecast at 95 %. As the forecast progresses, the interval becomes wider. For example, at time t there is less uncertainty from where the observed value of the SPEI can be located at time t+1 than from the value at time t+4 in which uncertainty is much greater.

It is also important to add that in the end we observed how the total errors of the prediction (of the process + of measurement) were adjusted to a t-Student distribution of mean 0, which approaches a normal distribution.

It is worth emphasizing that the Kalman filter is important for the forecast of hydrometeorological variables because daily mean flows can be predicted (^{Gonzalez et al, 2015}), as well as sub-hourly flows (^{Morales et al., 2014}), and predict indicators of droughts with some months of anticipation. The prediction of droughts allows preparing contingency plans to reduce their negative impact.