Zona de estudio

La subregión hidrológica del río Fuerte está localizada al noroeste de México, en la región hidrológica número 10 de Sinaloa. La Figura 1 muestra como su superficie está distribuida en Sinaloa, Chihuahua, Sonora y Durango.

Figura 1 Ubicación de la cuenca del río Fuerte y estaciones hidrometeorológicas. 

En la parte alta de la cuenca se encuentran los asentamientos Tarahumaras, cuyos niveles de pobreza son altos, y en la parte baja se encuentran tres presas (Miguel Hidalgo y Costilla, Luis Donaldo Colosio y Josefa Ortiz de Domínguez) con fines de generación de electricidad e irrigación (^{Castillo-Castillo et al., 2017}). El sistema de presas suministra agua a dos distritos de riego (DR): DR075 El Fuerte y DR076 El Carrizo, con una superficie total de 231 699 y 77 657 ha, respectivamente (^{CONAGUA, 2020}). Además, de acuerdo con la ^{CONAGUA (2015)}, esta subregión está expuesta a la incidencia de huracanes en verano.

Una vez revisada la información disponible sobre la cuenca, ésta se delimitó con la extensión HEC Geo-HMS 10.2 de ArcGis y con base en el modelo de elevación digital del Instituto Nacional de Estadística, Geografía e Informática (^{INEGI, 2019a}), hasta la estación 25009 Bocatoma Sufragio del Servicio Meteorológico Nacional (^{SMN, 2019a}). Con la misma extensión se obtuvieron las características de las subcuencas: pendiente, longitud de cauce, tiempo de retraso y esquema hidrológico HMS (esquema o dibujo con iconos). La zona de estudio se acotó hasta la estación Tubares, la cual comprende una superficie de 17 145.02 km². Las presas Huites, Miguel Hidalgo y Josefa Ortiz de Domínguez están fuera de la zona de estudio, inmediatamente aguas abajo de la estación hidrométrica Tubares. Tubares está en Chihuahua, aproximadamente a 39 km de la presa Huites, Sinaloa.

Información hidrometeorológica

La información de estaciones meteorológicas automáticas (EMAs) fue proporcionada por la Comisión Federal de Electricidad (^{CFE, 2019}) y el ^{SMN (2019b)}. Los meses adquiridos fueron de junio a octubre de 2009, 2011, 2015, 2016 y 2017. También se solicitó información pluviométrica de 24 h del SMN (Cuadro 1). La Figura 1 muestra la ubicación y distribución espacial de las estaciones meteorológicas e hidrométricas.

Una vez ordenada la información, y sabiendo que la serie de caudales horarios es continua en el tiempo, se utilizó el método de interpolación lineal para estimar los datos faltantes, siempre y cuando no fueran más de seis datos continuos. ^{Gordon et al. (2004)} mencionan que para periodos cortos se puede utilizar la interpolación lineal. Para las series de precipitación, al no ser continuas (ni temporal ni espacialmente), se utilizó el método del inverso del cuadrado de la distancia.

Distribución espacial y temporal de la precipitación para el modelo hidrológico

El modelo se trabajó a una escala horaria, ya que los datos de caudales medidos por la CFE así lo permitieron, considerando que gran parte de la cuenca es operada por dicha institución, al menos hasta las presas donde se genera energía eléctrica. Es importante aclarar que las presas quedaron aguas abajo del punto de calibración (Tubares, Chihuahua) del modelo hidrológico.

Se usaron datos horarios de lluvia de ocho estaciones meteorológicas automáticas: cuatro pertenecientes a la CFE y cuatro al SMN. Asimismo, se usó información pluviométrica diaria de cinco estaciones del SMN para realizar una mejor representación espacial de la lluvia en aquellos espacios que quedaban sin cubrir por las ocho EMAs. En las estaciones con valores acumulados de lluvia cada 24 h se empleó la distribución proporcional de precipitación en el mismo lapso que la EMA más cercana. Este manejo de datos fue el utilizado por ^{Haberlandt et al. (2008)}.

Para finalizar esta parte del proceso, la distribución espacial de lluvia se determinó con el método de polígonos de Thiessen. La información sobre los eventos modelados se muestra en el Cuadro 2.

La Organización Meteorológica Mundial (^{WMO, 2018}) considera que la medición de lluvia en un punto dado sólo representa una zona limitada, la cual se acota en función del periodo de acumulación, la homogeneidad fisiográfica de la región, la topografía y los diferentes procesos que propician la precipitación. ^{Espinosa-López et al. (2020)} modelaron la cuenca del río Huaynamota con los datos horarios de las EMAs de la CFE y con los de imágenes satelitales de lluvia horaria. El modelo con datos de las EMAs dio mejores resultados en términos de coeficiente de ajuste al momento de calibrar y validar; por ello, en el presente estudio se utilizaron datos de las EMAs.

Modelación hidrológica

La modelación hidrológica se realizó en el programa HEC-HMS 4.2.1, con el método del Servicio de Conservación de Suelo (SCS) de EUA para transformar la lluvia en lámina escurrida, y para calcular los hidrogramas a la salida de las subcuencas se utilizó el hidrograma unitario de Clark (^{U.S. Army Corps of Engineers [ USACE ], 2000}). El intervalo de tiempo de la modelación fue de 1 h. La importancia de construir un modelo horario radica en la relevancia que se tiene al momento de utilizar este modelo como una herramienta de pronóstico para fines de protección civil y para tener más certeza del momento en que ocurrirá el mayor caudal. También resulta importante tener un modelo horario en aquellos puntos de la cuenca donde el tiempo de concentración del escurrimiento es menor a 24 h, ya que, de lo contrario, el modelo no sería útil para pronosticar el tiempo del caudal máximo. El modelo calibrado, en términos horarios, serviría para tomar decisiones.

De acuerdo con ^{McCuen (2016)}, la ecuación para estimar la lámina escurrida es la siguiente:

donde Q es la lámina escurrida o lluvia exceso (mm), P es la lámina precipitada (mm), S es el potencial máximo de retención de la cuenca (mm) e I _a es la abstracción inicial (mm). La evidencia empírica lleva al supuesto de que I _a = 0.2S; considerando esto, la Ecuación 1 está condicionada a P ≥ I _a , de otra manera Q = 0. Además, el potencial máximo de retención se relaciona con el número de curva (CN) mediante la siguiente ecuación:

De acuerdo con ^{USDA (1986)}, el número de curva de escurrimiento está, principalmente, en función del grupo hidrológico del suelo (A, B, C o D), tipo de cobertura y condición de escurrimiento antecedente. De esta misma fuente se obtuvieron los valores de número de curva para las diferentes condiciones de la cuenca. La ^{CNA (1987)} explica la manera de elegir el grupo hidrológico acorde con la unidad de suelo reportada en la carta edafológica del ^{INEGI (2019b)}.

Para calcular el número de curva de escurrimiento, primero se creó una capa raster de número de curva utilizando los datos vectoriales de tipo de suelo serie II, y uso de suelo y vegetación serie VI; ambos se descargaron de la plataforma digital del ^{INEGI (2019c)}. El raster de número de curva se realizó en la extensión HEC-GeoHMS, como lo indican ^{USACE (2013)} en el manual de usuario de dicha extensión. Una vez calculados los números de curva, se realizó la corrección por pendiente de acuerdo con ^{Neitsch et al. (2011)}:

donde CN _2S es el número de curva para la condición II (CN ₂ ) ajustado por pendiente, es decir, la condición cuando el contenido de humedad se considera promedio, y slp es la pendiente media de la cuenca (%). ^{USACE (2013)} define tres condiciones antecedentes de humedad: I seco (punto de marchitamiento), II humedad promedio y III capacidad de campo.

Para cada tormenta se realizó la corrección por lluvia antecedente acumulada cinco días antes (^{Aparicio-Mijares, 2008}) (Cuadro 3). Las condiciones usadas para tal fin fueron: 1) corrección A: cuando la lluvia antecedente es menor a 25 mm, 2) sin corrección: cuando la lluvia antecedente es mayor a 25 y menor a 50 mm, y 3) corrección B: cuando la lluvia antecedente es mayor a 50 mm. Una manera práctica de aplicar las Ecuaciones 3 y 4 es la expuesta en el Cuadro 3.

Fuente: ^{Aparicio-Mijares (2008)}.

En el cálculo del caudal máximo y el volumen escurrido en un hidrograma se pueden utilizar los hidrogramas unitarios sintéticos (HUS), especialmente en cuencas no aforadas. En este estudio se utilizó el hidrograma unitario de Clark, el cual representa dos procesos de transformación de lluvia en exceso a escurrimiento: 1) traslación de la lluvia en exceso desde su origen a través de la red drenaje hasta la salida y 2) atenuación de la cantidad de la descarga una vez que el exceso se va almacenando en la cuenca (^{USACE, 2000}). ^{Straub et al. (2000)} mencionan que la atenuación del flujo se puede representar como un depósito lineal, cuyo almacenamiento se relaciona con la salida de la siguiente manera:

donde S es el almacenamiento en la cuenca (m³), R es el coeficiente de almacenamiento (h) y O es la salida del flujo (m³·s^-1). A partir de esta relación se puede llegar a la ecuación básica de hidrograma unitario de Clark, donde se agrega I como el flujo de entrada (m³·s^-1) y C como el coeficiente que está en función del paso de modelación (h):

Comúnmente los modelos hidrológicos utilizan un parámetro relacionado al tiempo, y el más utilizado es el tiempo de concentración (^{Campos-Aranda, 2010}). ^{Smith-Quintero y Velásquez-Henao (1995)} señalan que el tiempo de concentración es el tiempo que se requiere para que el agua llegue desde el punto más alejado de la cuenca hasta su punto de salida. ^{Grimaldi et al. (2012)} mencionan que, en la práctica, se usan fórmulas/ empíricas para determinar este parámetro; una es la de Giandotti, que es altamente usada en Italia, otra la de Kirpich y la del NRCS, que son usadas en EUA. En el Cuadro 4 se muestran diferentes fórmulas para el cálculo de dicho parámetro.

Cuadro 4 Resumen de métodos de estimación, en horas, del tiempo de concentración empíricos y semi-empíricos. 

Método	Ecuación	Observaciones
Carter	Tc=0.0015476 L0.6S-0.3	Desarrollada para cuencas urbanas.Área menor a 20.72 km2Longitud de canales menor a 11.27 kmImportante flujo de tuberías en la cuencaDesarrollada a partir de datos de drenaje de aeropuertos
FAA	Tc=0.01656261.1-CL0.5100S-0.333	Válida para pequeñas cuencas donde dominan el flujo de lámina y el terrestre
Fórmula de onda cinemática	Tc=0.0319639 L0.6n0.6i-0.4S-0.3	Desarrollada a partir del análisis de ondas cinemáticas de la escorrentía superficial
Kerby-Hathaway	Tc=0.024068Lc0.47n0.47Sc-0.235	Área < 0.0404686 km2 Pendiente < 1 % n < 0.8
Kirpich (Tenesee)	Tc=0.0003245 L0.77S-0.385	Área de 0.00405 a 0.4532 km2 Pendiente de 3 a 10 %
US SCS	Tc=0.0002269 L0.81000CN-90.7S-0.5	Para cuencas rurales donde predomine el flujo terrestreÁrea de cuenca < 8.0937 km2
Arizona DOT	Tc=0.0097956 A0.1L0.25Lca0.25S-0.2	Forma modificada de FAACuencas agrícolas
Giandotti-Fang	Tc=4A0.5+1.5Lc/0.8H-E00.5	Desarrollada para pequeñas cuencas agrícolas en Italia
California Culvert Practice	Tc=0.871(L')3∆H-1.00.385	Desarrollada para pequeñas cuencas montañosas en California
Témez	Tc=0.3LcS0.250.76

A = área (km²); C = coeficiente de escurrimiento; CN = número de curva de escurrimiento; E _o = elevación de la salida de la cuenca (m); H = elevación media de la cuenca (m); ΔH = diferencia de elevación del punto inicial y final de un canal (m); i = intensidad de lluvia (mm·h^-1); L = longitud de la cuenca a lo largo del canal principal desde el punto más distante hidráulicamente hasta la salida (m); L _c = longitud del canal principal (m); L _ca = longitud medida desde el punto de concentración a lo largo de L hasta el punto en L que es perpendicular al centroide de la cuenca (m); L’= Longitud del curso del agua (km); n = coeficiente de rugosidad de Manning; S = pendiente media de la cuenca (m·m^-1); Sc = pendiente media del canal (m·m^-1). Fuente: modificado de ^{Sharifi y Housenni (2011)}.

Mediante la extensión HEC-GeoHMS se puede calcular el tiempo de retraso con la siguiente expresión (^{USACE, 2013}):

donde Lag es el tiempo de retraso (h), CN es el número de curva, L es la longitud hidráulica de la cuenca (m) y y es la pendiente media de la cuenca (%). A este método se le llamará SCS de aquí en adelante. Es importante recordar que Lag = 0.6T _c , donde T _c es el tiempo de concentración. L es la longitud que corresponde al tiempo que tarda el agua en viajar desde el punto hidráulicamente más lejano, esto en términos de distancia y pendiente, ya que la hacen tomar más tiempo en llegar a la salida.

Existen diferentes estudios en los que tratan de encontrar la mejor fórmula para calcular el tiempo de concentración, como el de ^{Sharifi y Hosseini (2011)} y el ^{Vélez-Upegui y Botero-Gutiérrez (2011)}. En este trabajo, se modelaron los eventos con diferentes tiempos de concentración, y se eligió el que permitió obtener mejores coeficientes de ajuste en una modelación inicial, específicamente el que presentó mayores coeficientes NSE y menores valores de raíz del cuadrado medio del error (RMSE, por sus siglas en inglés) para el mismo evento.

En lo que respecta al hidrograma unitario de Clark (Ecuación 5), el ^{USACE (1967)} recomienda usar un coeficiente de almacenamiento (R) con valor de 0.8 veces el tiempo de concentración. ^{Zimmermann (2003)} menciona que, ante la ausencia de datos de caudales, se deben utilizar fórmulas como R = cT _c , donde c varía de 0.5 a 0.8. ^{Magaña-Hernández et al. (2013)} usaron un c de 0.75 al momento de modelar una cuenca en el norte de México. En este trabajo, se utilizó un valor de 0.75 veces el tiempo de concentración.

El tránsito en cauces se determinó con el método de Muskingum (Ecuación 9). De acuerdo con ^{Karamouz et al. (2013)}, en dicho método el hidrograma de salida en la parte baja de un río se calcula para un hidrograma determinado en la parte alta del mismo.

donde K es la constante de tiempo del viaje y X es el factor de ponderación que varía de 0 a 1. De acuerdo con ^{Aparicio-Mijares (2008)}, K es igual a:

donde L es la longitud del tramo (m) y ꙍ es la velocidad de la onda (m·s^-1). ^{Aparicio-Mijares (2008)} menciona que K se puede relacionar como 1.5 veces la velocidad media del agua en el cauce. Para este trabajo, se consideró una velocidad de referencia de 2.8 m·s^-1, que es mayor a la velocidad media de canales naturales no bien definidos en pendientes de 12 a 15 % (2.4 m·s^-1) de acuerdo con ^{Aparicio-Mijares (2008)}. El mismo autor considera que, a falta de información, se recomienda usar un valor de 0.2 para el factor de ponderación. Según ^{Heggen (1984)}, el valor de X comúnmente varía de 0.2 a 0.3; por ello, el valor utilizado fue de 0.2.

Finalmente, se realizó una auto calibración en el evento de 2009, lo que permitió un ajuste máximo de 20 % en los números de curva, tiempos de concentración y, por lo tanto, en el coeficiente de almacenamiento R. Esto se llevó a cabo con el mecanismo optimization trial en HEC-HMS. Posteriormente, se realizó una validación en los eventos de 2011, 2015, 2016 y 2017. La calibración y validación se llevó a cabo en tres estaciones hidrometeorológicas de la CFE: Urique, Guerachic y Tubares; sin embargo, no se tuvo éxito en las dos primeras (correspondientes a escurrimientos de subcuencas de la parte alta), por lo que se consideró como referencia el ajuste en la estación Tubares. Cabe mencionar que el mecanismo optimization trial, de manera general, sigue el procedimiento que se muestra en la Figura 2, en donde se seleccionó el método de gradiente univariado, y la función objetivo se evaluó con el coeficiente de NSE.

Figura 2 Esquema del procedimiento de calibración del modelo hidrológico (^{USACE, 2000}). 

De acuerdo con ^{USACE (2000)}, el algoritmo del gradiente univariado realiza correcciones consecutivas a los parámetros estimados. Si x ^k es el parámetro que se estima con la función objetivo f(x ^k ) en la iteración k, la búsqueda define un nuevo x ^k+1 en la iteración k + 1 como:

donde Δx ^k es la corrección del parámetro y debe ser tal que las estimaciones tiendan hacia el valor mínimo de la función objetivo. La función objetivo puede ser el ajuste de bondad medido con el coeficiente de NSE, RMSE u otros.

Ajuste de bondad del modelo hidrológico

Los indicadores más utilizados para valorar los ajustes de modelos son: el coeficiente NSE, la RMSE, el coeficiente de determinación (R ² ) y el error cuadrático medio (MSE, por sus siglas en inglés) (^{Vargas-Castañeda et al., 2015}). En el presente trabajo se utilizaron el NSE y el RMSE. Cabe mencionar que no se encontró un referente, estadísticamente hablando, de un buen modelo hidrológico horario. ^{Moriasi et al. (2007)} presentan una clasificación general de estadísticos recomendados para evaluar simulaciones hidrológicas con tamaño de paso mensual, la cual se divide de la siguiente manera: muy buena 0.75 < NSE ≤ 1, buena 0.65 < NSE ≤ 0.75, satisfactoria 0.50 < NSE ≤ 0.65 y no satisfactoria NSE ≤ 0.50. Estos autores también mencionan que valores entre 0 y 1 son considerados como niveles aceptables de ajuste, ya que un NSE > 0 significa que, por lo menos, el modelo es mejor predictor que usar simplemente la media aritmética.