SciELO - Scientific Electronic Library Online

vol.13 número1Comportamiento espacial de Hemiberlesia lataniae (Signoret) en aguacate ‘Hass’ en el Estado de MéxicoModelación hidrológica para medir el impacto del cambio climático en la cuenca del río Yaqui índice de autoresíndice de materiabúsqueda de artículos
Home Pagelista alfabética de revistas  

Servicios Personalizados




Links relacionados

  • No hay artículos similaresSimilares en SciELO


Ingeniería agrícola y biosistemas

versión On-line ISSN 2007-4026versión impresa ISSN 2007-3925

Ing. agric. biosist. vol.13 no.1 Chapingo ene./jun. 2021  Epub 13-Jun-2022 

Scientific article

Hourly hydrologic modeling in the upper basin of the Fuerte River, Sinaloa, Mexico

Edgar Merino-Jiménez1

Laura Alicia Ibáñez-Castillo1  *

Ramón Arteaga-Ramírez1

Mario Alberto Vázquez-Peña1

1Universidad Autónoma Chapingo. Carretera México-Texcoco km 38.5, Chapingo, Estado de México, C. P. 56230, MÉXICO.



Mexico uses hydrological models to determine floods, evaluate land use change scenarios, evaluate climate change scenarios, and define federal zones, among other applications. However, the models are rarely calibrated beforehand, which increases uncertainty in the design of structures and hydraulic standards.


To build a hydrological model for the watershed of the Fuerte River, Mexico, of extreme rainfall events occurred in 2009, 2011, 2015, 2016 and 2017.


Five extreme rainfall events were considered for this study. The hydrologic model was design using the HEC-HMS program, and calibrated at the Tubares hydrometric station. The runoff curve number methodology and the Clark unit hydrograph were used.


The results collected in four of the five events were positive; the Nash-Sutcliffe efficiency (NSE) ranged between 0.22 and 0.52. The temporal behavior of the river, at times when it moved away from the flow value, preserved the variation trend.

Limitations of the study:

The study reaches the Tubares hydrometric station, Chihuahua, without including the downstream dams in Sinaloa.


There are few hydrological studies that generate a calibrated and, therefore, reliable hourly model.


The hourly hydrologic model had an acceptable performance in four of the five predicted events in terms of NSE and root mean square error (RMSE).

Keywords Clark unit hydrograph; runoff curve number; time of concentration with the California Culvert Practice equation; calibration and validation of surface hydrologic models



En México se utilizan los modelos hidrológicos para determinar avenidas, evaluar escenarios de cambio de uso de suelo, evaluar escenarios de cambio climático, definir zonas federales, entre otras aplicaciones. Sin embargo, rara vez los modelos se calibran previamente, lo que incrementa la incertidumbre en el diseño de estructuras y normas hidráulicas.


Construir un modelo hidrológico para la cuenca del río Fuerte, México, de eventos de lluvia extremos ocurridos en 2009, 2011, 2015, 2016 y 2017.


Se consideraron cinco eventos de lluvia extremos. El modelo hidrológico se construyó con el programa HEC-HMS, y se calibró en la estación hidrométrica Tubares. Se utilizó la metodología de número de curva de escurrimiento y el hidrograma unitario de Clark.


Los resultados obtenidos en cuatro de los cinco eventos fueron buenos; la eficiencia de Nash-Sutcliffe (NSE) osciló entre 0.22 y 0.52. El comportamiento temporal del cauce, en momentos en que se alejaba del valor del caudal, conservaba la tendencia de variación.

Limitaciones del estudio:

El estudio llega hasta la estación hidrométrica Tubares, Chihuahua, sin incluir las presas que están aguas abajo en Sinaloa.


Existen pocos estudios hidrológicos que generen un modelo horario calibrado y, por lo tanto, confiable.


El modelo hidrológico horario construido tuvo un desempeño aceptable en cuatro de los cinco eventos modelados en términos del NSE y la raíz del cuadrado medio del error (RMSE).

Palabras clave hidrograma unitario de Clark; número de curva de escurrimiento; tiempo de concentración con la ecuación California Culvert Practice; calibración y validación de modelos hidrológicos superficiales


Hydrological phenomena are complex and difficult to understand, and in the absence of exact knowledge they can be represented in a simplified form by means of the system concept (Weber et al., 2012). A system can be understood as a set of elements or components related to each other by different processes (Rodríguez et al., 2006).

Due to the need to determine the effects of climate variation, climate change, land use changes and water resource management, several hydrological simulation models have been developed (Van Liew et al., 2005). The United States Department of Agriculture (USDA, 1986) mentions that hydrological studies should be based on historical flow records, but these are rarely available and statistical analysis becomes inaccurate due to land use changes. Therefore, Aparicio-Mijares (2008) indicates that it is necessary to have methods that estimate runoff from the characteristics of the watershed and precipitation, which are known as rainfall-runoff methods.

Although there are different programs for hydrologic modeling, Estrada-Sifontes and Pacheco-Moya (2012) mention that the Hydrologic Modeling System of the Hydrologic Engineering Center of the US Army Corps of Engineers (HEC-HMS) is a flexible program that allows selecting different methods for the calculation of losses, hydrographs, and river traffic; in addition, it allows simulating at the event level and on a continuous basis. Magaña-Hernández et al. (2013) predicted with HEC-HMS, the Escondido river watershed in Coahuila, with an area of 3 242 km2, for three hydrometeorological events, which are among the few with radar data (NEXRAD from NOAA of the USA) for hourly rainfall. These authors had almost perfect Nash-Sutcliffe coefficients (NSE), which was due to the fact that, at the time of calibrating the model, they left the range in which the hydrological parameters could move very flexible. The NSE measures the deviation with respect to the unit of the ratio between the mean square error and the variance of the observations (Singh, 2017).

In the hydrological region of the Fuerte River, Mexico, there is an important infrastructure for storing and supplying water for different uses, and for generating electricity and controlling floods. However, the incidence of extreme phenomena (such as hurricanes and droughts) is present, being unpredictable and with effects that can cause very serious economic damages (Comisión Nacional del Agua [ CONAGUA ], 2015).

Therefore, the objective of this paper was to build a calibrated and validated hydrological model, at hourly scale in subbasins belonging to the hydrological region of the Fuerte River, for extreme rainfall events of 2009, 2011, 2015, 2016 and 2017.

Materials and methods

Study area

The hydrological subregion of the Fuerte River is located in northwestern Mexico, in the hydrological region number 10 of Sinaloa. Figure 1 shows how its surface area is distributed in Sinaloa, Chihuahua, Sonora, and Durango.

Figure 1 Location of the Fuerte river watershed and hydrometeorological stations. 

In the upper part of the watershed there are Tarahumara settlements, whose poverty levels are high, and in the lower part there are three dams (Miguel Hidalgo y Costilla, Luis Donaldo Colosio and Josefa Ortiz de Domínguez) for electricity generation and irrigation purposes (Castillo-Castillo et al., 2017). The dam system supplies water to two irrigation districts (RD): DR075 El Fuerte and DR076 El Carrizo, with a total area of 231 699 and 77 657 ha, respectively (CONAGUA, 2020). In addition, according to CONAGUA (2015), this subregion is exposed to the incidence of hurricanes in summer.

Once the available information on the watershed was reviewed, it was delimited with the HEC Geo-HMS 10.2 extension of ArcGis and based on the digital elevation model of the Instituto Nacional de Estadística, Geografía e Informática (INEGI, 2019a), up to station 25009 Bocatoma Sufragio of the Servicio Meteorológico Nacional (SMN, 2019a). With the same extension, the characteristics of the subbasins were obtained: slope, river length, lag time and HMS hydrological scheme (scheme or drawing with icons). The study area was delimited up to the Tubares station, which covers an area of 17 145.02 km2. The Huites, Miguel Hidalgo and Josefa Ortiz de Dominguez dams are outside the study area, immediately downstream of the Tubares hydrometric station. Tubares is in Chihuahua, approximately 39 km from the Huites Dam, Sinaloa.

Hydrometeorological information

Information from automatic weather stations (AWS) was provided by the Comisión Federal de Electricidad (CFE, 2019) and the SMN (2019b). The months reported were from June to October 2009, 2011, 2015, 2016 and 2017. Twenty-four h rainfall information was also requested from the SMN (Table 1). Figure 1 shows the location and spatial distribution of meteorological and hydrometric stations.

Table 1 Hydrometeorological stations in the upper watershed of the Fuerte River up to the Tubares station. 

Station Administration/Type State Coordinate X, UTM (m) Coordinate Y, UTM (m) Information Scale
Tubares CFE/AWS Chihuahua 800062.98 2983759.21 Flow Precipitation 1 h
Urique CFE/AWS Chihuahua 805824.31 3014204.43 Flow Precipitation 1 h
Batopilas CFE/AWS Chihuahua 823729.33 2993087.05 Flow Precipitation 1 h
Guerachic CFE/AWS Chihuahua 873443.76 2950383.58 Flow Precipitation 1 h
Urique SMN/AWS Chihuahua 805386.86 3014070.40 Precipitation 10 min
Chinatu SMN/AWS Chihuahua 922682.80 2907997.05 Precipitation 10 min
El Vergel SMN/AWS Chihuahua 959766.68 2936360.79 Precipitation 10 min
Guachochi SMN/AWS Chihuahua 890436.31 2971802.14 Precipitation 10 min
Maguarichi SMN/AWS Chihuahua 795976.38 3085139.09 Precipitation 10 min
8038 Creel SMN/ Pluviometric Chihuahua 831477.88 3074041.35 Precipitation 24 h
8106 Norogachi SMN/ Pluviometric Chihuahua 883587.54 3024230.93 Precipitation 24 h
8172 Guadalupe y Calvo SMN/ Pluviometric Chihuahua 902646.84 2893908.09 Precipitation 24 h
8266 Batovira SMN/ Pluviometric Chihuahua 804700.54 3091986.56 Precipitation 24 h
10125 La Rosilla SMN/ Pluviometric Durango 966840.86 2898310.47 Precipitation 24 h

CFE = Comisión Federal de Electricidad; AWS = automatic weather stations; SMN = Servicio Meteorológico Nacional; UTM = Transverse Mercator Coordinates.

Once the information was sorted and knowing that the series of hourly flow rates is continuous over time, the linear interpolation method was used to estimate the missing data, as long as there were no more than six continuous data. Gordon et al. (2004) mention that for short periods, linear interpolation can be used. For precipitation series, since they are not continuous (neither temporally nor spatially), the inverse of the square of the distance method was used.

Spatial and temporal distribution of precipitation for the hydrological model

The model was developed on an hourly scale, since the flow data measured by the CFE allowed it, taking into account that a large part of the watershed is operated by this institution, at least up to the dams where electric power is generated. It is important to clarify that the dams were located downstream of the calibration point (Tubares, Chihuahua) of the hydrological model.

Hourly rainfall data from eight automatic weather stations were used: four belonging to the CFE and four to the SMN. In addition, daily rainfall data from five SMN stations were used for a better spatial representation of rainfall in those areas that were not covered by the eight AWS. In the stations with cumulative rainfall values every 24 h, the proportional distribution of precipitation in the same period as the nearest AWS was used. This data management was used by Haberlandt et al. (2008).

To complete this part of the process, the spatial distribution of rainfall was determined using the Thiessen polygon method. The information on the predicted events is shown in Table 2.

Table 2 Hydrometeorological events selected for modeling the Fuerte river watershed. 

Storms Start Eind Duration (days) Maximum flow measured in Tubares (m3·s-1) Total rainfall over the watershed (mm)
2009 06/08/2009 08:00 30/08/2009 08:00 24 1 284 3 039
2011 18/08/2011 08:00 24/08/2011 08:00 6 296 504
2015 03/09/2015 08:00 18/09/2015 08:00 15 507 992
2016 16/07/2016 08:00 08/08/2016 08:00 23 935 2 074
2017 01/08/2017 09:00 09/08/2017 08:00 8 1 500 534

The World Meteorological Organization (WMO, 2018) considers that rainfall measurement at a given point only represents a limited area, which is bounded according to the accumulation period, physiographic homogeneity of the region, topography, and the different processes that favor precipitation. Espinosa-López et al. (2020) predicted the Huaynamota river watershed with hourly data from the CFE AWS and hourly satellite rainfall images. The model with AWS data gave better results in terms of fitting coefficient at the time of calibration and validation; therefore, data from the AWS were used in the present study.

Hydrologic modeling

Hydrologic modeling was performed using the HEC-HMS 4.2.1 program, with the U.S. Soil Conservation Service (SCS) method to transform rainfall into runoff, and the Clark unit hydrograph (U.S. Army Corps of Engineers [ USACE ], 2000) was used to calculate the hydrographs at the outflow of the subbasins. The modeling time interval was 1 h. The importance of building an hourly model lies in the relevance of using this model as a forecasting tool for civil protection purposes and to be more certain of the time when the highest flow will occur. It is also important to have an hourly model at those points in the watershed where the runoff concentration time is less than 24 h, since otherwise the model would not be useful for forecasting the time of the maximum flow. The calibrated model, in hourly terms, would be useful for decision making.

According to McCuen (2016), the equation for estimating the runoff is as follows:

Q=P-Ia2P+0.8S (1)

where Q is the runoff or excess rainfall (mm), P is the precipitated rainfall (mm), S is the maximum retention potential of the watershed (mm) and I a is the initial abstraction (mm). Empirical evidence leads to the assumption that I a = 0.2S; considering this, Equation 1 is conditioned to P ≥ I a , otherwise Q = 0. Furthermore, the maximum retention potential is related to the curve number (CN) by the following equation:

S=25400CN-254 (2)

According to USDA (1986), the runoff curve number is mainly a function of the hydrological group of the soil (A, B, C or D), type of cover and antecedent runoff condition. Curve number values were obtained from the same source for the different watershed conditions. CNA (1987) explains how to choose the hydrologic group according to the soil unit reported in the edaphological map of INEGI (2019b).

To calculate the runoff curve number, a curve number raster layer was first created using the vector data of soil type series II, and land use and vegetation series VI; both were downloaded from the INEGI digital platform (2019c). The curve number raster was performed in the HEC-GeoHMS extension, as indicated by USACE (2013) in the user manual of that extension. Once the curve numbers were calculated, the correction for slope was performed according to Neitsch et al. (2011):

CN2S=CN3-CN23×1-2×exp-13.86slp+CN2 (3)

CN3=CN2×exp0.00673100-CN2 (4)

where CN 2S is the curve number for condition II (CN 2 ) adjusted for slope, i.e., the condition when moisture content is considered average, and slp is the average slope of the watershed (%). USACE (2013) defines three antecedent moisture conditions: I dry (wilting point), II average moisture and III field capacity.

For each storm, the correction was made for antecedent rainfall accumulated five days before (Aparicio-Mijares, 2008) (Table 3). The conditions used for this purpose were: 1) correction A: when the antecedent rainfall is lower than 25 mm, 2) no correction: when the antecedent rainfall is greater than 25 and lower than 50 mm, and 3) correction B: when the antecedent rainfall is greater than 50 mm. A practical way to apply Equations 3 and 4 is shown in Table 3.

Table 3 Corrections of the curve number for antecedent moisture. 

No correction Correction A Correction B
0 0 0
10 4 22
20 9 37
30 15 50
40 22 60
50 31 70
60 40 78
70 51 85
80 63 91
90 78 96
100 100 100

Source: Aparicio-Mijares (2008).

Synthetic unit hydrographs (SUH) can be used to calculate the maximum flow and runoff volume in a hydrograph, especially in ungauged watersheds. This study used the Clark unit hydrograph, which represents two processes of transformation of excess rainfall to runoff: 1) translation of excess rainfall from its origin through the drainage network to the outflow point and 2) attenuation of the amount of discharge once the excess is stored in the watershed (USACE, 2000). Straub et al. (2000) mention that flow attenuation can be represented as a linear reservoir, where storage is related to the outflow as follows:

S=RO (5)

where S is the storage in the watershed (m3), R is the storage coefficient (h) and O is the outflow (m3·s-1). From this relationship, the basic Clark unit hydrograph equation can be derived, where I is added as the inflow (m3·s-1) and C as the coefficient that is a function of the modeling step (h):

O2=CI+1-CO (6)

C=2t2R+t (7)

Hydrological models commonly use a time-related parameter, and the most used parameter is time of concentration (Campos-Aranda, 2010). Smith-Quintero and Velásquez-Henao (1995) indicate that the time of concentration is the time required for water to reach from the farthest point of the watershed to the exit point. Grimaldi et al. (2012) mention that, in practi/ce, empirical formulas are used to determine this parameter; one is that of Giandotti, which is highly used in Italy, another is that of Kirpich and that of NRCS, which are used in the USA. Table 4 shows different formulas to calculate this parameter.

Table 4 An overview of empirical and semi-empirical time-of-concentration estimation methods in hours. 

Method Equation Comments
Carter Tc=0.0015476 L0.6S-0.3 Developed for urban watersheds.Area lower than 20.72 km2Length of channels lower than 11.27 kmSignificant pipe flow in the watershedDeveloped from airport drainage data.
FAA Tc=0.01656261.1-CL0.5100S-0.333 Valid for small watersheds where laminar and overland flow dominate
Fórmula de onda cinemática Tc=0.0319639 L0.6n0.6i-0.4S-0.3 Developed from the analysis of surface runoff kinematic waves
Kerby-Hathaway Tc=0.024068Lc0.47n0.47Sc-0.235 Area < 0.0404686 km2 Slope < 1 % n < 0.8
Kirpich (Tenesee) Tc=0.0003245 L0.77S-0.385 Area of 0.00405 a 0.4532 km2 Slope from 3 to 10 %
US SCS Tc=0.0002269 L0.81000CN-90.7S-0.5 For rural watersheds where overland flow predominates Watershed area < 8.0937 km2
Arizona DOT Tc=0.0097956 A0.1L0.25Lca0.25S-0.2 Modified form of FAAAgricultural watersheds
Giandotti-Fang Tc=4A0.5+1.5Lc/0.8H-E00.5 Developed for small agricultural watersheds in Italy
California Culvert Practice Tc=0.871(L')3H-1.00.385 Developed for small mountainous watersheds in California
Témez Tc=0.3LcS0.250.76

A = area (km2); C = runoff coefficient; CN = runoff curve number; E o = watershed outlet elevation (m); H = average watershed elevation (m); ΔH = difference in elevation between the start and end point of a channel (m); i = rainfall intensity (mm·h-1); L = length of the watershed along the main channel from the most hydraulically distant point to the outflow point (m); L c = length of the main channel (m); L ca = length measured from the concentration point along L to the point in L that is perpendicular to the watershed centroid (m); L’= Watercourse length (km); n = Manning's roughness coefficient; S = average watershed slope (m·m-1); Sc = average channel slope (m·m-1). Source: modified by Sharifi and Housenni (2011).

Using the HEC-GeoHMS extension, delay time can be calculated with the following expression (USACE, 2013):

Lag=0.00136L0.81000CN-90.7y-0.5 (8)

Where Lag is the lag time (h), CN is the curve number, L is the hydraulic length of the watershed (m) and y is the average slope of the watershed (%). Hereafter, this method will be called SCS. It is important to remember that Lag = 0.6T c , where T c is the time of concentration. L is the length that corresponds to the time it takes for the water to travel from the hydraulically farthest point, this in terms of distance and slope, since water takes more time to reach the exit.

There are different studies that try to find the best formula to calculate the time of concentration, such as Sharifi and Hosseini (2011) and Vélez-Upegui and Botero-Gutiérrez (2011). The events with different times of concentration were predicted in this study, and the one with the best fitting coefficients in an initial modeling was chosen, specifically the one with higher NSE coefficients and lower root mean square error (RMSE) values for the same event.

Regarding the Clark unit hydrograph (Equation 5), the USACE (1967) recommends using a storage coefficient (R) with a value of 0.8 times the time of concentration. Zimmermann (2003) mentions that, in the absence of flow data, formulas such as R = cT c should be used, where c varies from 0.5 to 0.8. Magaña-Hernández et al. (2013) used a c of 0.75 when modeling a watershed in northern Mexico. This study used a value of 0.75 times the time of concentration.

Streamflow was determined using the Muskingum method (Equation 9). According to Karamouz et al. (2013), in that method the outflow hydrograph in the lower part of a river is calculated for a hydrograph determined in the upper part of the river.

S=KXI+1-XQ (9)

where K is the travel time constant, and X is the weighting factor that varies from 0 to 1. According to Aparicio-Mijares (2008), K is equal to:

K=Lω (10)

where L is the length of the section (m) and is the wave velocity (m·s-1). Aparicio-Mijares (2008) mentions that K can be related as 1.5 times the mean water velocity in the channel. This study considered a reference velocity of 2.8 m·s-1, which is higher than the average velocity of natural channels not well defined on slopes of 12 to 15 % (2.4 m·s-1) according to Aparicio-Mijares (2008). The same author considers that, in the absence of information, it is recommended to use a value of 0.2 for the weighting factor. According to Heggen (1984), the value of X commonly varies from 0.2 to 0.3; therefore, the value used was 0.2.

Finally, a self-calibration was performed in the 2009 event, which allowed a maximum fit of 20 % in curve numbers, concentration times and, therefore, in the storage coefficient R. This was carried out with the optimization trial mechanism in HEC-HMS. Subsequently, validation was performed on the 2011, 2015, 2016 and 2017 events. Calibration and validation were carried out in three hydrometeorological stations of the CFE: Urique, Guerachic and Tubares; however, it was not successful in the first two (corresponding to runoff from subbasins of the upper part), so fitting in the Tubares station was considered as a reference. It is worth mentioning that the optimization trial mechanism follows the procedure shown in Figure 2, where the univariate gradient method was selected, and the objective function was evaluated with the NSE coefficient.

Figure 2 Diagram of the hydrologic model calibration procedure (USACE, 2000). 

According to USACE (2000), the univariate gradient algorithm performs consecutive corrections to the estimated parameters. If x k is the parameter estimated by the objective function f(x k ) at iteration k, the search defines a new x k+1 at iteration k + 1 as:

xk+1=xk+xk (11)

where Δx k is the parameter correction and should be such that the estimates tend toward the minimum value of the objective function. The objective function can be the goodness-of-fit measured with the NSE coefficient, RMSE or others.

Goodness-of-fit of the hydrological model

The most used indicators to assess model fitting are NSE coefficient, RMSE, coefficient of determination (R 2 ) and mean squared error (MSE) (Vargas-Castañeda et al., 2015). The present study used NSE and RMSE. No reference of a good hourly hydrological model was found. Moriasi et al. (2007) present a general classification of recommended statistics to evaluate hydrological simulations with monthly step size, which is divided as follows: very good 0.75 < NSE ≤ 1, good 0.65 < NSE ≤ 0.75, satisfactory 0.50 < NSE ≤ 0.65 and unsatisfactory NSE ≤ 0.50. These authors also mention that values between 0 and 1 are considered acceptable levels of fit, since an NSE > 0 means that, at least, the model is a better predictor than simply using the arithmetic mean. Moriasi et al. (2007) show a general classification of recommended statistics for evaluating hydrologic simulations with monthly step size, which is divided as follows: very good 0.75 < NSE ≤ 1, good 0.65 < NSE ≤ 0.75, satisfactory 0.50 < NSE ≤ 0.65 and unsatisfactory NSE ≤ 0.50. These authors also mention that values between 0 and 1 are considered acceptable levels of fit, since an NSE > 0 means that, at least, the model is a better predictor than using the arithmetic mean.

Results and discussion

Figure 3 shows the subbasins corresponding to the modeling carried out, and their characteristics are shown in Table 5. The subbasins, up to the Tubares hydrometric station, have an average slope of 40.47 % and an average curve number of 79.04.

Figure 3 Fuerte river watershed and location of hydrometric stations. 

Table 5 Characteristics of the subbasins of the study area. 

Subbasin Direct contribution to Area (km2) Slope (%) CN2 CN2S Channel lenght (km) Channel slope (m·m-1)
W1050 Urique 4 263.48 33.75 75.09 79.58 222.99 0.007
W510 Batopilas 2 001.54 37.81 75.02 79.55 90.35 0.016
W680 Guerachic 1 588.99 42.39 74.65 79.26 68.41 0.005
W710 Guerachic 1 632.40 42.02 74.62 79.22 80.90 0.014
W720 Guerachic 1 661.52 36.10 74.76 79.31 86.16 0.012
W770 Guerachic 1 395.40 32.06 77.51 81.63 75.01 0.017
W1250 Guerachic 42.57 62.11 77.54 81.76 2.93 0.009
W1040 Tubares 423.17 58.32 72.11 77.07 53.57 0.005
W840 Tubares 360.12 55.54 73.75 78.50 48.83 0.005
W580 Tubares 356.40 47.58 72.18 77.12 28.95 0.002
W610 Tubares 462.07 44.76 73.33 78.11 10.08 0.032
W620 Tubares 588.85 50.76 67.98 73.42 32.68 0.001
W640 Tubares 95.94 55.97 65.53 71.23 14.75 0.007
W650 Tubares 1 195.16 48.14 72.62 77.51 57.77 0.008
W660 Tubares 506.87 56.11 71.06 76.15 31.15 0.004
W690 Tubares 396.43 46.47 78.68 82.71 58.64 0.000
W1240 Tubares 161.21 53.98 74.02 78.72 19.62 0.008
W900 Tubares 12.92 30.91 72.50 77.28 5.43 0.001

CN2 = curve number for condition II; CN2S = curve number for condition II fitted per slope

Concentration times calculated with different formulas are shown in Table 6. The maximum values are those calculated with the Giandotti and California Culvert formula, and the minimum values are those obtained with the SCS formula. Regarding the above, the events were predicted using these three times of concentration.

Table 6 Time of concentration, in hours, derived from different formulas. 

Subbasin Giandotti California Culverts Kirpich Témez SCS
W1040 13.81 13.67 11.02 7.89 5.50
W510 12.30 16.12 10.39 14.77 11.24
W840 13.28 12.60 10.34 7.34 4.95
W580 28.11 21.50 10.33 7.22 5.32
W610 11.42 9.80 1.48 7.01 5.09
W620 33.46 21.98 14.18 7.04 5.65
W640 9.27 5.91 3.66 3.55 2.84
W650 16.43 17.78 9.57 11.30 8.40
W660 19.21 16.39 7.51 7.64 5.52
W680 19.87 18.73 12.99 11.15 8.15
W690 39.35 29.31 28.60 7.79 4.87
W1240 9.21 6.08 4.27 4.09 2.68
W710 12.28 14.88 10.09 12.87 9.52
W720 14.44 19.41 11.32 15.35 11.96
W770 11.29 14.58 8.84 13.77 10.28
W900 18.16 7.75 3.92 2.12 1.66
W1050 21.52 37.18 28.93 26.49 21.50
W1250 12.11 6.43 0.94 2.64 1.48

SCS = Soil Conservation Service.

Table 7 shows the NSE coefficients derived from different concentration times. In 14 predicted events a better fit was obtained when the California Culvert method was used, three were better with SCS and two with Giandotti. The above coincides with the RMSE estimated, so the time of concentration calculated with the California Culvert Practice formula was used in this study.

Table 7 Nash-Sutcliffe coefficients calculated with the concentration times of Soil Conservation Service (SCS), California Culvert Practice and Giandotti. 

Event Station Nash-Sutcliffe coefficients Best fit method
SCS California Giandotti Maximum
2009 Tubares -7.34 -4.36 -4.89 -4.36 California
2011 Tubares -44.14 -17.37 -18.68 -17.37 California
2015 Tubares -92.51 -74.89 -85.20 -74.89 California
2016 Tubares -22.32 -20.07 -22.31 -20.07 California
2017 Tubares -2.73 -2.31 -2.31 -2.31 California

The values calculated in the 2009 event calibration are shown in Table 8, which were applied to each event and resulted in the curve numbers shown in Table 9.

Table 8 Calibration results for the 2009 event. 

Subbasin Curve number Time of concentration Storage coefficient R (h)
Modification (%) Value Modification (%) Value (h)
W1040 -1.71 * -11.49 12.10 9.07
W1050 0.00 * -8.11 34.15 25.61
W1240 4.86 * -12.23 5.34 4.00
W1250 1.70 * 0.74 6.48 4.86
W510 14.87 * 6.25 17.12 12.84
W580 8.58 * -11.41 19.04 14.28
W610 -1.28 * -8.50 8.96 6.72
W620 0.00 * -11.98 19.34 14.51
W640 0.00 * -10.79 5.28 3.96
W650 -4.23 * -8.38 16.29 12.22
W660 -1.76 * -9.98 14.75 11.06
W680 0.80 * -6.09 17.58 13.19
W690 -2.31 * -9.63 26.48 19.86
W710 -19.11 * -10.42 13.33 10.00
W720 0.00 * -8.78 17.70 13.27
W770 -1.00 * -7.63 13.46 10.10
W840 -0.01 * -8.01 11.59 8.69
W900 14.99 * 8.14 8.37 6.28

*Variable value depending on antecedent rainfall (see Table 9).

Table 9 Curve number fitted for each event. 

Subbasin 2009 2011 2015 2016 2017
W1040 58.47 58.47 75.75 58.47 87.71
W1050 62.49 62.49 62.49 62.49 79.58
W1240 64.46 94.62 64.46 64.46 64.46
W1250 66.75 93.44 66.75 66.75 66.75
W510 71.74 91.37 71.74 71.74 91.37
W580 64.65 64.65 64.65 64.65 83.73
W610 59.96 77.11 59.96 59.96 77.11
W620 55.11 55.11 73.42 55.11 73.42
W640 52.47 71.22 52.47 52.47 71.22
W650 57.47 85.72 57.47 57.47 74.23
W660 57.35 87.13 57.35 57.35 57.35
W680 62.61 79.89 62.61 62.61 62.61
W690 65.52 90.22 65.52 65.52 65.52
W710 50.21 64.08 50.21 50.21 64.08
W720 79.31 62.17 62.17 62.17 90.58
W770 64.79 64.79 64.79 64.79 64.79
W840 61.20 78.50 61.20 61.20 78.50
W900 68.69 68.69 88.86 88.86 88.86

Fittings on the curve numbers showed three cases: 1) decrease in eight subbasins (two of them runoff directly into the Guerachic station and the others directly into Tubares), 2) four subbasins had no changes and 3) two subbasins had direct runoff into Tubares, and there was increase in the Batopilas subbasin, in two subbasins of Guerachic and in three subbasins of Tubares. On average, the decrease was 1.76 % and the increase was 7.63 %.

With respect to time of concentration, an average decrease of 9.56 % was recorded in 15 of the 18 subbasins, and an increase of 5.04 % in three subbasins (Batopilas, one in Tubares and one in Guerachic), which represent 12 % of the total study area.

The NSE coefficients and RMSE for calibration and validations are shown in Table 10. Figures 4, 5, and 6 show three of the five predicted events. The results show an adequate fit at the Tubares station, with coefficients above zero. A negative value was obtained in the 2016 storm; however, the predicted hydrograph is temporally similar to that measured (Figure 6).

Table 10 Hydrologic model fit coefficients for the Fuerte River at the Tubares station. 

Coefficient 2009 2011 2015 2016 2017
NSE 0.42 0.26 0.52 -0.28 0.22
RMSE (m3·s-1) 248.13 58.28 74.03 243.84 337.08

NSE = Nash-Sutcliffe efficiency; RMSE = root mean square error.

Figure 4 Hydrograph of Tubares station, Chihuahua (August 2009 event). 

Figure 5 Hydrograph of Tubares station, Chihuahua (September 2015 event). 

Figure 6 Hydrograph of Tubares station, Chihuahua (July-August 2016 event). 

The maximum flow rates measured and calculated are shown in Table 11. To contrast the visual interpretation (Figures 4, 5 and 6) and the model performance efficiencies (Table 10), Figures 7, 8, 9 and 10 show the behavior of that measured versus that predicted flow rates of the 2009, 2015, 2016 and 2017 events. Although the 2016 event had an NSE of -0.28, Figure 9 shows that the R2 is 0.59. This is important because if there were no error between that measured and predicted the trend of this relationship would be a straight line; therefore, the modeling of that 2016 event is not as negative as the NSE would suggest.

Table 11 Flows measured and predicted from the Tubares gauging station. 

Event Hydrograph Maximum flow (m3·s-1)
2009 Predicted 1 282.50
Measured 1 283.70
2011 Predicted 394.10
Measured 296.00
2015 Predicted 402.50
Measured 506.50
2016 Predicted 1 387.10
Measured 934.80
2017 Predicted 853.90
Measured 1 499.90

Figures 7, 8, 9 and 10 show the goodness-of-fit for three hydrometeorological events from 2009, 2015, 2016 and 2017, respectively.

Figure 7 Goodness of fit between that measured and predicted in the August 2009 event in the upper Rio Fuerte watershed. 

Figure 8 Goodness of fit between that measured and predicted in the September 2015 event in the upper Rio Fuerte watershed. 

Figure 9 Goodness of fit between that measured and predicted in the July/August 2016 event in the upper Rio Fuerte watershed. 

Figure 10 Goodness of fit between that measured and predicted in the August 2017 event in the upper Fuerte River watershed. 

The model fit would be expected to have a relationship with the maximum flow measured; however, the results show no influence. For example, in the calibration of the 2009 event, an NSE coefficient of 0.42 was obtained for a maximum discharge measured of 1 283.70 m3·s-1 (Figure 7), while the 2015 validation had an NSE of 0.52, corresponding to a maximum flow measured of 506.50 m3·s-1 (60.54 % lower) (Figure 8), and the 2017 event, with a flow of 1 499.90 m3·s-1 (16.84 % higher), had an NSE of 0.22 (Figure 10). Moriasi et al. (2007) consider, for a calibration and validation on a daily scale, a value greater than 0.4 as satisfactory.

If we compare the NSE of 0.4 of this study with the 0.95 of the hourly hydrological model of Magaña-Hernández et al. (2013), for the Escondido river watershed, Coahuila, it would seem that this study is at a disadvantage. However, in the Escondido river model, too much correction margin was given to the hydrological parameters at the time of calibration. That is, the hydrological parameters were allowed to vary up to any value and at any level of increase or decrease, regardless of the reality of their values, which gave them an almost perfect model. Another aspect is that in that study they had some radar images of rainfall provided by the USA due to the proximity to Texas. It is convenient to improve the hourly hydrological modeling due to its importance in civil protection and in those points where the time of concentration of the watershed is less than 24 h. However, modeling should be done preserving the principles of the modeler representing the reality it intends to simulate.

The visual analysis of the Tubares station shows temporal similarity between the predicted and measured hydrographs. Considering the coefficients obtained, it can be said that the model found at the outlet of the Tubares subbasin represents, at an acceptable level (Moriassi, et al, 2007), the hydrometeorological phenomenon in the study area, because it shows a satisfactory fit in calibration and in one of the validations.

An attempt was also made to calibrate the model at two upstream hydrometric stations but was unsuccessful. One possible explanation is that rainfall varies more spatially in mountain areas. Moreover, WMO (1972) recommends one station per 100 to 250 km2 for mountainous areas, and in the study watershed there are 14 stations (eight AWS and four pluviometric stations with daily data) distributed over 17 145 km2; that is, one station covers an average of 1 224.64 km2, which leads us to think that rainfall may not be well represented in this area.

On the other hand, good fitting at the outlet is attributed to two aspects: 1) the transit of the flow from the upstream stations until reaching the outlet and 2) the contributions of 11 subbasins (lower part), which are not calibrated individually (as is the case of the upstream subbasins Urique and Batopilas), but up to the Tubares station, which allows compensating the bad fitting caused by the contributions of the upstream subbasins.

Haberlandt et al. (2008) mention that simple hydrological models provide an objective global estimate, but at the cost of smoothing the hydrograph, which leads to underestimation or overestimation; this is observed in the hydrographs measured in the present study. Perhaps this model could be improved by considering a distributed model that uses information from satellite images of hourly rainfall and that is calibrated with data from the CFE. Another improvement to the model would be to generate measured unit hydrographs, instead of the synthetic Clark unit hydrographs used in this study.

Although one might believe that it is not worthwhile to build hourly hydrological models in Mexico because there is not enough measured rainfall and hourly flow data, it is important to prepare for the future and strive to develop a sufficient and reliable hydro-meteorological network.


The hydrologic model developed at the outlet of the Tubares watershed using the curve number method to calculate losses and Clark's synthetic unit hydrograph method for hydrographs, is acceptable based on the Nash-Sutcliffe coefficients.

The time of concentration that yielded the best modeling results is the one calculated with the California Culvert Practice formula, which overestimates the correct value by about 10 %. An alternative to improve modeling could be the use of an aggregated model that uses information from satellite or radar images to estimate rainfall. It is also recommended to model individually the upper subbasins with different methods than those used in the present study, for example, a different unit hydrograph or a different spatial distribution of rainfall. All this to obtain a general model that properly represents hydrology in the entire watershed.


The authors thank the Consejo Nacional de Ciencia y Tecnología (CONACyT) for funding this study, the Comisión Federal de Electricidad and the Servicio Meteorológico Nacional for providing hydrometeorological information from their databases


Aparicio-Mijares, F. J. (2008). Fundamentos de hidrología de la superficie. Limusa. ]

Campos-Aranda, D. F. (2010). Introducción a la hidrología urbana. Universidad Autónoma de San Luis Potosí. [ Links ]

Castillo-Castillo, M., Ibáñez-Castillo, L. A., Valdés, J. B., Arteaga-Ramírez, R., & Vázquez-Peña, M. A. (2017). Análisis de sequías meteorológicas en la cuenca del río Fuerte, México. Tecnología y Ciencias del Agua, 8(1), 35-52. [ Links ]

Comisión Federal de Electricidad (CFE). (2019). Información hidrometeorológica. [ Links ]

Comisión Nacional del Agua (CNA). (1987). Instructivo de Hidrología para determinar la avenida máxima ordinaria asociada a la delimitación de la zona federal. ]

Comisión Nacional del Agua (CONAGUA) (2015). Acuerdo por el que se da a conocer el resultado de los estudios técnicos de aguas nacionales superficiales de la Subregión Hidrológica Río Fuerte de la Región Hidrológica número 10 Sinaloa. ]

Comisión Nacional del Agua (CONAGUA). (2020, October 2). Si stema Nacional de Información del Agua. ]

Espinosa-López, J. A., Ibáñez-Castillo, L. A., Arteaga-Ramírez, R., & Galeana-Pizaña, J. M. (2020). Modelo hidrológico distribuido con imágenes GPM-IMERG en la cuenca río Hueynamota, Nayarit, México. Tecnologia y Ciencias del Agua, 11(5), 276-302. [ Links ]

Estrada-Sifontes, V., & Pacheco-Moya, R. M. (2012). Modelación hidrológica con HEC-HMS en cuencas montañosas de la región oriental de Cuba. Ingeniería Hidráulica y Ambiental, 33(1), 71-80. [ Links ]

Gordon, N. D., McMahon, T. A., Finlayson, B. L., Gippel, C. J., & Nathan, R. J. (2004). Stream hydrology: an introduction for ecologists. Jhon Wiley & Sons Ltd. [ Links ]

Grimaldi, S., Petroselli, A., Tauro, F., & Porfiri, M. (2012). Time of concentration: a paradox in modern hydrology. Hydrological Sciences Journal, 57(2), 217-228. [ Links ]

Haberlandt, U., Ebner von Eschenbach, A. D., & Buchwald, I. (2008). A space-time hybrid hourly rainfall model for derived flood frequency analysis. Hydrology and Earth System Sciences, 12(6), 1353-1367. [ Links ]

Heggen, R. J. (1984). Univariate least squares Muskingum Flood Routing. Water Resources Bulletin, 20(1), 103-107. [ Links ]

Instituto Nacional de Estadística y Geografía (INEGI). (2019a, January 12). Continuo de elevaciones mexicano 3.0 (CEM 3.0). ]

Instituto Nacional de Estadística y Geografía (INEGI). (2019b, January 12). Conjunto de datos vectoriales edafológicos serie II. ]

Instituto Nacional de Estadística y Geografía (INEGI). (2019c, January 12). Conjunto de datos vectoriales de uso de suelo y vegetación serie VI. ]

Karamouz, M., Nazif, S., & Falahi, M. (2013). Hydrology and hydroclimatology: Principles and applications. CRC Press. [ Links ]

Magaña-Hernández, F., Bâ, K. M., & Guerra-Cobián, V. H. (2013). Estimación del hidrograma de crecientes con modelación determinística y precipitación derivada de radar. Agrociencia, 47(8), 739-752. ]

McCuen, R. H. (2016). Hydrologic analysis and design. Pearson. [ Links ]

Moriasi, D. N., Arnold, J. G., Van Liew, M. W., Bingner, R. L., Harmel, R. D., & Veith, T. L. (2007). Model evaluation guidelines for systematic quantification of accuracy in watershed simulations. American Society of Agricultural and Biological Engineers, 50(3), 885-900. [ Links ]

Neitsch, S. L., Arnold, J. G., Kiniry, J. R., & Williams, J. R. (2011). Soil & water assessment tool theorical documentatino version 2009. Texas A&M University. ]

Rodríguez, A., Pagot, M. R., Hillman, G. D., Pozzi, C. E., Plencovich, G. E., Caamaño, N. G., … Bucher, E. H. (2006). Modelo de simulación hidrológica. In Bucher, E. H. (Ed), Bañados del río dulce y laguna mar Chiquita (pp. 57-77). Academia Nacional de Ciencias. ]

Servicio Meteorológico Nacional (SMN). (2019a, January 12). Información estadística climatológica. ]

Servicio Meteorológico Nacional (SMN). (2019b, November). Información meteorológica. [ Links ]

Sharifi, S., & Hosseini, S. M. (2011). Methodology for identifying the best equations for estimating the time of concentration of watersheds in a particular region. Journal of Irrigation and Drainage Engineering, 137(11), 712-719. [ Links ]

Singh, V. P. (2017). Handbook of applied hydrology. McGraw-Hill Education. [ Links ]

Smith-Quintero, R. A., & Velásquez-Henao, J. D. (1995). Modelos conceptuales en hidrología para la estimación del hidrograma unitario. Avances en Recursos Hidráulicos, 3, 55-72. ]

Straub, T. D., Melching, C. S., & Kocher, K. E. (2000). Equation for estimating Clark Unit-Hydrograph parameters for small rural watersheds in Illinois. U.S. Department of Interior - U.S. Geological Survey. [ Links ]

U.S. Army Corps of Engineers (USACE). (2000). Hydrologic modeling system HEC-HSM technical reference manual. U.S. Army Corps of Engineering, Hydrologic Engineering Center. ]

U.S. Army Corps of Engineers (USACE). (2013). HEC-GeoHMS geospatial hydrologic modeling extension user’s manual. U.S. Army Corps of Engineers. ]

U.S. Army Corps of Engineers (USACE). (1967). Generalized standard project rain flood criteria southern California Coastal Streams, Sacramento. U.S. Army Corps of Engineers. ]

United States Department of Agriculture (USDA). (1986). Urban hydrology for small wathersheds TR-55. United States Department of Agriculture. ]

Van Liew, M. W., Arnold, J. G., & Bosch, D. D. (2005). Problems and potential of autocalibrating a hydrologic model. American Society of Agricultural Engineers, 48(3), 1025-1040. [ Links ]

Vargas-Castañeda, G., Ibáñez-Castillo, L. A., & Arteaga-Ramírez, R. (2015). Development, classification and trends in rainfall-runoff modeling. Ingeniería Agrícola y Biosistemas, 7(1), 5-21. [ Links ]

Vélez-Upegui, J. J., & Botero-Gutiérrez, A. (2011). Estimation of the time of concentration and the lag time at San Luis creek basin, Manizales. Dyna, 78(165), 58-71. ]

Weber, J. F., Jorquera, E., & Stehli, P. T. (2012). Simulación de la respuesta hidrológica de una cuenca ante el movimiento de una tormenta. Mecánica Computacional, 31(4), 579-603. ]

World Meteorological Organization (WMO). (1972). Casebook on hydrological network design practice. World Meteorological Organization. ]

World Meteorological Organization (WMO). (2018). Guía de instrumentos y métodos de observación meteorológicos. World Meteorological Organization. ]

Zimmermann, E. D. (2003). A generalization of Clark´s IUH for flatland areas with strong human interventions. Journal of Environmental Hydrology, 11(2), 1-14. ]

Received: November 23, 2020; Accepted: April 29, 2021

*Corresponding author:, tel. 595 12 55 581.

Creative Commons License This is an open-access article distributed under the terms of the Creative Commons Attribution License