Assessing the Performance of Two Hydrologic Models for Forecasting Daily Streamflows in the Cazones River Basin ( Mexico )

Floods have caused significant human and economic losses in the Cazones River Basin, located on the Gulf of Mexico. Despite this knowledge, steps towards the design and implementation of an early warning system for the Cazones are still a pending task. In this study we contributed by establishing a hydrological scheme for forecasting mean daily discharges in the Cazones Basin. For these purposes, we calibrated, validated and compared the HyMod model (HM) which is physics-based, and an autoregressive-based model coupled with the Discrete Kalman Filter (ARX-DKF). The ability of both models to accurately predict discharges proved satisfactory results during the validation period with RMSEHYMOD = 2.77 [mm/day]; and RMSEARX-DKF = [2.38 mm/day]. Further analysis based on a Streamflow Assimilation Ratio (SAR) revealed that both models underestimate the discharges in a similar proportion. This evaluation also showed that, under the most common conditions, the simpler stochastic model (ARX-DKF) performs better; however, under extreme hydrological conditions the deterministic HM model reveals a better performance. These results are discussed under the context of future applications and additional requirements needed to implement an early warning hydrologic system for the Cazones Basin.


Introduction
The interaction between global warming and natural variability affects the magnitude and frequency of extreme hydrological events, i.e. droughts and floods [1]- [3].Future hydroclimatic scenarios associated with these impacts are influenced by numerous sources of uncertainty, motivating the scientific community to work on better forecasting of hydrological events [4].Hydrological modeling techniques are an essential tool for the development of efficient water governance schemes of various types i.e., hydraulic works management, mitigation of natural disasters, and risk management, among others [4]- [11].
In Mexico, communities around river basins within humid-subtropical climates could especially benefit from accurate hydrological modeling.These humid-subtropical climates result from the interaction between several ocean-atmosphere teleconnections, i.e., the El Niño Southern Oscillation (ENSO), the Pacific Decadal Oscillation (PDO), the Arctic Oscillation (AO), the North Atlantic Oscillation (NAO), and the Atlantic Multidecadal Oscillation (NAO).These ocean-atmosphere interactions generate flooding and landslides that can result in important socio-economic losses, particularly in the floodplains of the basins.
To mitigate the impacts associated to these extreme events, during 2002 the Mexican government established an Early Warning System for Tropical Cyclones (EWSTC), which is the main tool used to monitor the intensity, trajectory and distance of the cyclones.The information generated by EWSTC is then used by federal and state agencies who take supportive measures.Despite this centralized system has demonstrated to be effective for early warning of meteorological disasters, the role of local agencies is still pending in terms of expanding the "hydrological knowledge" about it e.g.calibration and validation of hydrologic models that deal with specific physical characteristics of the basins or those sustained in new algorithms of data assimilation.
In order to advance towards the development of an early warning hydrological system for the Cazones Basin, in this study we have established a first approximation for forecasting mean daily discharges in the basin.For these purposes, we calibrated, validated, and compared two lumped hydrological models The first model is HyMod (HM), a physics-based model that accounts for all components of the water balance [12]- [16].The second is an autoregressive-based model coupled with the Discrete Kalman Filter (ARX-DKF); the latter is considered as the optimal sequential data assimilation method for linear dynamics and measurement processes with Gaussian error statistics [17]- [21].In the following sections of this paper we describe the study area and dataset, the calibration and validation steps for both models, and the discussion and conclusions concerning the results of our study.

Study Area
The Cazones River Basin is located between the coordinates 20˚18' and 21˚15'N and 97˚17' and 98˚32'W, in the Hydrological Region N˚27 North, called Veracruz Tuxpan-Nautla [22] [23].The mean annual flow of the Cazones River is 54 m 3 /s, but for a single day this value can reach a maximum of about 1250 m 3 /s [24].The river has a total drainage area of approximately 2688 km 2 distributed in three states of Mexico (Figure 1(a)): Puebla (56%), Hidalgo (13.5%), and Veracruz (30.5%) [25].For this study we considered the Poza Rica stream gauge as the outlet of the basin (Figure 1(b)), accounting for a total area of 1,605 km 2 .Grasslands and agricultural fields represent more than 70% of the total land area of the basin (Figure 1(c)).The average slope of the basin is about 27% (Figure 1(d)), and the highest elevation is about 2880 m.a.s.l.(Figure 1(e)).The maximum height of the main channel is about 2,100 m.a.s.l.These geomorphological characteristics translate into an estimated time of concentration of about 16 hours [25] [26].

Hydroclimatology of the Cazones Basin
Much of Mexico has a monsoon climate, with a rainy season during the summer months (JJA) and a relatively dry winter season (DJF) [27].Climate anomalies in Mexico depend mostly on the interaction between ENSO, PDO, AO, NAO, and AMO [28]- [33].In the positive phase of the El Niño phenomenon (i.e., a warmer Equatorial Pacific Ocean), the winter is humid throughout Mexico, whereas during the summer, conditions are dry in the north and humid in the south [31].A humid-subtropical climate dominates the Cazones River Basin, with an average annual temperature of 24˚C and an average annual rainfall of ~1500 mm [34].This precipitation accumulation is associated with frontal systems during winter, and extreme meteorological events such as tropical storms or hurricanes that form in the Atlantic over the Gulf of Mexico during the summer [29] [30] [35] [36].

Socio-Economic Issues in the Cazones Basin
The Cazones Basin drains over its west-facing slopes.The drainage enters the lowlands around the town of Poza Rica in the state of Veracruz, an important petroleum production area.Extreme meteorological events affecting the area cause high runoff over the saturated soils of the basin, resulting in significant damage [25] [26] [34] [37].A review of the history of the state of Veracruz reveals extensive human losses due to flooding and landslides.For instance, in October 1999 a tropical depression associated with a polar front produced a 50-year flooding event that affected more than 200,000 inhabitants and caused more than $150 million (US) in economic losses [38].In 2005, four tropical storms affected about 1.5 million people, and 130,000 houses were damaged by flooding [33].During 2010, two hurricanes in the same month, Karl (September 17-18) and Mathew (September 26-27,) caused flooding in different parts of the State and economic losses of over $5 billion (US) [33]- [39].During the last 41 years, flooding in Veracruz has demonstrated a significant positive exponential trend (r = 0.81) (see Figure 2), which can be attributed to climate change (due to global warming), other anthropogenic factors, and natural variability [39] [40].

Hydrologic Modeling Using HM and ARX-DKF
Daily streamflow forecasting in the Cazones Basin was conducted using two lumped models (one physics-based and one systems-based).In the definition of HM model we considered the formal steps of model building as stated by [41] i.e., the general model formulation (conceptual model) was performed defining the state variables; the model fluxes; and the adjustable and derived parameters.For the ARX model, we defined the exogenous inputs and the rainfall-runoff transfer function component (see i.e. [42]- [44]).This model was then coupled with a Discrete Kalman Filter, following the steps outlined by [17] [45]- [48].

Forcing Data
The input forcing data used in this study include daily, spatially aggregated rainfall amounts (2008-2011) obtained from a total of nine raingauges, mean daily flows obtained from Poza Rica stream gauge (Table 1), and   potential evapotranspiration (PET) calculated using Hargreaves Method [49].These records were obtained from the database CLICOM (Climate Computerized) administered by the National Meteorological Service of Mexico [50].Mean daily streamflows from the Cazones Basin were obtained from the BANDAS database [51] and then used in the calibration and validation schemes.

Conceptual Representation (Process Model)
The physical or mathematical principles governing the entire behavior of the hydrologic system were summarized for each model through the use of simple conceptual Directed Graphs (DG), which are an explicit representation of the major processes occurring within the hydrological system, their structural organization, and their results [41].Because the HM is a physics-based model (see e.g., [12] [14] [16] [52], the major physical processes leading to the behavior of the hydrologic system were defined in order to clearly describe the dynam-ics by which the state variables change over time [41].The HM conceptualizes the basin as a simple, nonlinear soil moisture reservoir connected with two series of linear reservoirs for channeling excess rainfall, two or three identical quick-flow tanks (depending on the calibration scheme applied) representing short-term surface detention, and one, slow-flow tank representing groundwater storage (Figure 3(a)).
On the other hand, physical conditions are not considered in the formulation of the ARX-DKF model, which is intended to find the best mathematical or causal relation(s) between the input and the output of the hydrological system [46] [53]- [55].This coupled model is based on an autoregressive component (ARX) for the streamflow, in which rainfall is considered the exogenous input.This component is then merged with the Discrete Kalman Filter (DKF) algorithm, following the steps described by [17] [48] [56]- [58].The resulting model is then used to forecast the mean daily stream flow records (Figure 3(b)).

System Parametrization
This stage consisted of developing a group of hypotheses regarding the mathematical forms of the process equations that are believed to describe the physical processes linking the subsystem components (see i.e. [41] [59]).We used the HM as a conceptual non-distributed (lumped) hydrological model with six parameters (see Table 2) accounting for all subcomponents of the water balance ( [12] [16] [60]).The ARX-DKF model contains timevarying parameters (i.e., they evolve over time) obtained from the autoregressive component, which are then used in the Kalman implementation (Table 2).The autoregressive characteristics of the streamflow records allow for the development of a matrix of parameters that are used to forecast streamflow [55].

Quick Routing
Reservoir 1 (Sq1) e + is the error for the streamflow forecast in time t + 1  α is a parameter calculated from the autoregressive component applied to streamflow. β is a parameter calculated from the autoregressive component applied to rainfall. na is the lag number for streamflow. nb is the lag number for rainfall.

Discrete Kalman
Filter Component x + is a vector containing the present streamflow (not observed). k x is a vector containing the streamflow in the time k. k u is a vector containing the rainfall in the time k. k w and k v are vectors containing the Gaussian noise for the process and measurements.

( )
 k z is a vector containing the streamflow measurement.
 A and B are matrices containing α and β parameters from the series of streamflow and rainfall data in the ARX model. H is a transformed matrix that contains the states of the measurements. Q and R are matrices containing the covariance for the noise in the process and the measurements.

Computational Model, Calibration, and Validation
The computational implementation of both models was coded using a Matlab TM routine, which was packaged in a toolbox1 format to facilitate usage.From the historical records of mean daily streamflows, the year 2008 was selected to calibrate the initial parameters of both models, because it contained the largest number of flood events (see Figure 2).The Mean Squared Error (MSE), Root Mean Squared Error (RMSE), and Nash-Sutcliffe Efficiency (NSE) [61] [62] were used as error measures for both the calibration and validation schemes (Table 3).The parameters of HM were calibrated from 2008 records using the Shuffle Complex Evolution Algorithm (SCE-UA) [63]- [65].Since the parameters of the ARX-DKF model are time-varying, they were initially calibrated using 2008 records, and then updated using 2009-2010 records (see Figure 3(b)).For these purposes we used an autoregressive component (ARX) taking into account the principles of parsimony [66]- [68].To ensure consistent analysis, both models were using the 2009-2011 period.During this validation, the performance of each model was further evaluated by calculating a Streamflow Assimilation Ratio Index (SAR) which is proposed in this study.This index shows the relationship between modeled (Q Mod ) and observed (Q Obs ) streamflows for each day i under analysis . In this study, SAR was used for measuring the predictive capacity of each model; however, it can also be used to apply bias correction schemes in a similar way of that performed over Global or Regional Climate Models (GCMs or RCMs); that is, by calculating an averaged version of SAR (averaged Bias Correction Factor) over all n days under analysis: and then using this factor to correct raw data resulting from hydrological models.

Sensitivity Analysis and Calibration of Parameters
For the HM model the selection of the final set of parameters from the 2008 calibration was used to define the physical processes associated with hydrological responses in the Cazones Basin.For instance, the MSE HM obtained fell as low as 8.79 mm, and all the parameters converged after completing the optimization scheme (Figure 4).For example, the calibrated maximum storage capacity in the basin was found to be about 39.4 mm, which may be considered a small capacity given the humid-subtropical characteristics of the Basin and the wide range of capacity (5 -600 mm) used to calibrate the max m S parameter.Additionally, as denoted by the α-parameter, the surface runoff in the basin is partitioned in slightly dissimilar proportions between the surface (53%) and the subsurface component (47%) revealing a higher runoff component compared to infiltration rates.In fact, Table 3. Forecasting performance measures used in this study.

Test Equation Description
(1) Mean Squared Error (MSE) ( ) (2) Root Mean Squared Error (RMSE) ( ) The evolution scheme was applied using five complexes and 25 loops at each complex (for details see: [63] [64] [73]. the land use in the basin which is dominated by grasslands and agricultural fields (>70%) must play a significant role in modulating the hydraulic properties and consequently the infiltration capacity of soils.For instance, some studies have found that grasslands have lower (higher) infiltration (runoff) rates than forests [69].
For the ARX-DKF model the best set of parameters α and β calculated during the calibration period and then used to predict daily streamflows was identified when we combined the streamflow of the previous day (Q t-1 ) and the rainfall that occurred in the prior two days (P t-2 ).This means that the largest autocorrelation of the observed streamflows is achieved at lag t − 1 (days), and the largest cross-correlation between rainfall and observed streamflows is achieved at lag t − 2 (days).Similar findings have been reported by [70]- [72].

Classical Validation of Models
During the validation period, the performance of the chosen parameters was tested, revealing that both models are able to satisfactorily reproduce the mean daily discharges in the Basin.In fact, the forecasting performance of both models was improved during the validation period as evaluated through error measures.The MSE HM was about 35% larger than MSE ARX-DKF and these results were also consistent when looking at RMSE and NSE measures (Table 4).The relationship between modeled and observed discharges shows significant correlations: 0.88 for HM and 0.90 for ARX-DKF (Figure 5); however, as indicated by the slight deviation of the least squares line (LSL) from the 1:1 line, the predicted discharges from both models tend to slightly underestimate the observed discharges.In general, all previous findings suggest, at least at this stage, that ARX-DKF timevariable parametrization performs better than HM parametrization.

Streamflow Assimilation Ratio (SAR)
A significant relationship between SAR HM and SAR ARX-DKF (r = 0.68) revealed interesting additional details about the performance of the models (in a validation period of 982 days).For instance, the joint distribution of SAR revealed four different regions (Figure 6(a)) that allow further evaluation of their performance.The Overestimation (Underestimation) Region is the one in which both models overestimate (underestimate) the observed discharges.In the Mixed Region, one model overestimates while the other underestimates or vice versa.The analysis of the Overestimation Region denoted similar performance by both models (Figure 6(b)).However, when the Underestimation Region was analyzed, ARX-DKF performed better than HM (Figure 6(c)): in about 76% of the days in the validation period, the latter model underestimated the observed discharges, with an average relationship represented by HM 0.42 Obs Q Q = × .In contrast, the ARX-DKF model underestimated the discharges in only about 66% of the days in the validation period, with an average relationship represented by . Focusing on the ability of the models to predict extreme conditions, we found that HM performs slightly better than ARX-DKF when predicting large discharge values: the average relationship between predicted and observed discharges (>20 mm/day) for each model is represented by HM 0.88 . The implementation and analysis of SAR seems to be a practical and more detailed approach to evaluate and compare the performance of hydrological models; however, more data and new applications will be required to determine its real significance in a context of models comparisons.

Hydrological Uncertainty
Judging from the hydrographs of both models (Figure 7), the predicted discharges can realistically forecast the mean daily discharges.The uncertainty bands calculated through a Monte Carlo approach (i.e., bootstrapping of errors) can be satisfactorily used to minimize both models' risk of inaccuracy.However, during extreme flooding events-i.e. an event in 2011 that resulted in an historic peak of ~1250 m 3 /s (~67.2 mm/day) on July 17 th (Figure 7)-the uncertainty bands can be not able to explain the observed discharge; this can be partially attributed to the fact that point-based rainfall cannot account for the spatial distribution of rainfall in the entire basin.Therefore, as a way of improving the confidence of the predictions and reduce the hydrological uncertainty, the design of an early warning hydrological system will additionally require the validation of combined new satellite-based precipitation products (gridded datasets) and lumped or distributed models.

Discussion
The use of predictions from the centralized EWSTC allows for the generation of pre-alert scenarios with larger lead time; however, this scheme does not contain information about the states of the hydrological system i.e. it cannot be used to determine possible damages of hydraulic infrastructure.This is the most important pending task of local agencies as a way to expand the hydrological knowledge of EWSTC.This knowledge will have to develop new early warning hydrological systems dedicated to represent hydrological response based on the intrinsic characteristics of the Basins or dependent on new algorithms of data assimilation.In this study we contributed to establish a hydrological scheme for forecasting mean daily discharges in the Cazones Basin.For instance, based on classical validation schemes and the application of SAR, we determined that under normal con- ditions ARX-DKF model performs slightly better than HM; however, this latter model performs better under extreme conditions.These hydrological statements and the information generated from both models can be used to warn the emergency coordinators in a different way than that currently offered by EWSTC, and thus are of great value to pre-manage a probable hydraulic emergency.An additional and important component of this validation process will be based on the coordination between public agencies, decision-makers, and stakeholders, and how these members make use of these tools and the information generated by the models to define appropriate policies towards the development of early warning hydrological systems at catchment scale i.e. temporal resolution of interest, selection of models, and ranges of applicability.

Conclusion
This study evaluated two lumped hydrological models to contribute in advancing towards the development of an early warning hydrologic system in the Cazones River Basin.HM and ARX-DKF were calibrated, validated, and compared to determine the best modeling system for forecasting mean daily discharges in the basin.The results from error measures and SAR revealed that the performance of the autoregressive-based ARX-DKF model (less underestimation) was slightly better than the performance of the physics-based HM model (greater underestimation).Despite these differences, under extreme hydrological conditions HM performs slightly better than ARX-DKF.The implementation of SAR proposed for this study showed to be a practical alternative to evaluate and compare the performance of hydrological models, and it seems to be promising for applying bias correction schemes over raw predictions obtained from hydrological models.The application of these methods is aimed to improve the predictions and reduce the uncertainty of hydrological models of course future studies will have to deal with the validation RCMs for the basin i.e. precipitation estimates provided by real-time satellite-based products as CMORPH, PERSSIANN or TMPA, combined with different lumped and spatially distributed hydrological models i.e.SACRAMENTO, HYMOD, HBV, TOMODEL, among others.The success or failure of these tools and the use of the information they provide will continue being a challenge since hydrologists have to also deal with climate models outputs coming with an additional source of uncertainty; therefore, the level of coordination attained between national agencies, decision-makers and stakeholders, in developing policies aimed to adequately implement new early warning systems including maybe not all but several scenarios of uncertainty, will play a central role for improving water governance schemes at catchment scale.

Figure 1 .
Figure 1.(a) Cazones River Basin; (b) rain gauges and Poza Rica stream gauge; (c) land use in the basin; (d) slope map of the basin; (e) digital elevation map (30 meters) of the basin.

Figure 2 .
Figure 2. Annual precipitation obtained from the Poza Rica rain gauge and the number of flood events per year in the Cazones Basin for the period 1970-2011.Source: Flood occurrence data was obtained from DesInventar (2015).

Table 1 .
Rain and stream gauges used in this study.The code in the first column indicates the gauges as numbered in Figure 1(b).Source: Clicom database [50].

Figure 3 .
Figure 3. Simple conceptual directed graphs for the (a) HM and (b) ARX-DKF models.

Figure 4 .
Figure 4. Evolution of parameters and objective function (MSE) by applying the SCE-UA optimization scheme during 2008.The evolution scheme was applied using five complexes and 25 loops at each complex (for details see:[63] [64][73].

Figure 5 .
Figure 5. Linear relationships between predicted and observed streamflows in the Cazones Basin.

Figure 6 .
Figure 6.(a) Scatter plot between the SAR calculated for both HM and ARX-DKF models; (b) distribution of the Overestimation Region of both models; (c) Distribution of the Underestimation Region of both models.

Figure 7 .
Figure 7. Observed and simulated mean daily river discharge [mm] for the 2009-2011 validation period at the Poza Rica stream gauge in the Cazones River Basin.The shaded areas represent 95% uncertainty levels (confidence intervals) obtained from Monte Carlo resampling (1000 samples) of the absolute errors of prediction under the assumption that they are normally distributed.When the lower boundaries reached negative values, they were adjusted to be 20% of the modeled discharge.
is the total number of days.