Zona de estudio

La cuenca del río Yaqui, que abarca una superficie de 65 024 km², pertenece a la región hidrológica número 9 y está situada en la zona centro del estado de Sonora, México. Esta cuenca drena al río Yaqui, considerado como el más grande en el noroeste de México (^{Comisión Nacional del Agua [CONAGUA], 2018}). El río Yaqui se origina en el sur de Arizona, EUA, y desemboca en el Océano Pacífico formando un delta costero que se aprovecha como área agrícola en diferentes municipios del sur de Sonora, lo que se conoce hoy en día como “Valle del Yaqui” (^{Enríquez-Licón & Lorenzana-Durán, 2017}; ^{Padilla-Calderón, 2017}).

A partir de compuestos diarios en malla del sistema CLICOM del Centro de Investigación Científica y de Educación Superior de Ensenada (CICESE), se dividió la cuenca en 12 cuadrantes y se obtuvo un histórico de 32 años de datos meteorológicos (1972-2003). La precipitación acumulada fue de 609 mm, siendo los meses de junio a octubre los que aportaron más en este parámetro (75 %) (Figura 1). Por su parte, la temperatura media anual fue de 16.8 °C, con un máximo de 33 °C (mayo a septiembre) y mínimo de 0 °C (diciembre a febrero) (^{CICESE, 2020}), por lo que el clima se puede considerar de semiárido a árido (^{Mounir et al., 2019}).

Figura 1 Climograma de la cuenca del río Yaqui, México (estadísticas 1974-2003). 

La cuenca del río Yaqui se compone de tres subcuencas hidrológicas (Figura 2), las cuales constituyen las áreas de captación de los embalses Lázaro Cárdenas “La Angostura” (13 956 km²), Plutarco Elías Calles “El Novillo” (39 499 km²) y Álvaro Obregón “El Oviachic” (11 569 km²) (^{Mounir et al., 2019}). En este estudio, se modeló el escurrimiento superficial de cada subcuenca.

Figura 2 Ubicación espacial de subcuencas del río Yaqui. 

Procesamiento de los valores de entrada

La herramienta de evaluación de suelos y agua (SWAT), desarrollada por ^{Arnold et al. (1998)} y ^{Neitsch et al. (2011)}, se utilizó para la modelación hidrológica de la cuenca del río Yaqui, en donde la variable simulada fue el escurrimiento superficial y la escala temporal de información fue mensual. En la Figura 3 se describe a detalle el proceso realizado, el cual aborda aspectos sobre los valores de entrada del modelo, la calibración, la validación y los escenarios de cambio climático.

Figura 3 Diagrama de flujo de la modelación hidrológica del río Yaqui mediante SWAT. 

^{Arnold et al. (2012)} reportan los diferentes datos que alimentan al modelo hidrológico SWAT para realizar la simulación (Cuadro 1), por lo que es necesario tener un orden y referencia en cuanto a su escala, coordenadas y unidades de medida.

Aplicación del modelo y parametrización

En este estudio se utilizó la interfaz QSWAT 2.0 del programa libre QGIS 2.6.1, en la cual se realizó la aplicación y parametrización del modelo hidrológico. Los valores de entrada se procesaron de acuerdo con el tamaño de la cuenca del río Yaqui a partir de un polígono oficial proporcionado por la CONAGUA, al cual se le asignó el sistema de coordenadas Datum WGS84, proyección UTM zona 12. El programa determinó la red de drenaje y las subcuencas que componen a la cuenca hidrográfica. Se estableció un punto de salida en cada presa, ya que cada una es reguladora del escurrimiento aguas arriba. Al mismo tiempo, se crearon unidades de respuesta hidrológica (HRU) para cada subcuenca, las cuales presentaban homogeneidad en cuanto a sus características físicas (tipo de suelo, vegetación y topografía). En este contexto, cada HRU tenía en su matriz variables diferentes que influyeron en la modelación (^{Luo et al., 2011}).

Al momento de parametrizar el modelo, es necesario indicar el periodo de simulación. En este caso, se utilizó un conjunto de datos meteorológicos de 34 años (1970-2003), los cuales se reportaron diariamente y se emplearon para el periodo de calentamiento (Warm up) (1970-1971), calibración (1972-1993) y validación (1994-2003). El periodo de calentamiento es necesario al inicio del modelo, ya que éste tiende a presentar gran incertidumbre debido a que se desconocen las condiciones iniciales de simulación (^{Daggupati et al., 2015}).

Calibración y validación

La calibración consiste en estimar los parámetros de un modelo hidrológico y comparar los datos simulados contra los datos observados (^{Moriasi et al., 2007}). Para tal efecto, se utilizó el programa SWAT-CUP 2019 (Calibration and Uncertainty Program) y se siguió el protocolo de calibración automatizada recomendado por ^{Abbaspour et al. (2015)} y ^{Arnold et al. (2012)}. En dicho protocolo se menciona que la parametrización es un paso importante en el modelo de calibración, y se debe basar en el conocimiento de los procesos hidrológicos y la variabilidad del tipo de suelo, uso del suelo, pendiente y ubicación de la subcuenca. Se seleccionó el algoritmo Sequential Uncertainty Fitting (SUFI2) con el fin de ajustar los parámetros de la modelación a partir de la disponibilidad de valores relacionados con el modelado del balance hídrico (^{Gholami et al., 2016}). Para cada subcuenca, se ejecutaron 1 500 simulaciones, las cuales se dividieron en tres periodos de iteraciones (500 simulaciones cada una) y se reajustaron manualmente los parámetros después de cada iteración (^{Montecelos-Zamora et al., 2018}).

Los parámetros de sensibilidad que se utilizaron en el proceso de calibración (Cuadro 2) representan los procesos de escurrimiento superficial en la subcuenca hidrológica (^{Niraula et al., 2015}; ^{Niraula et al., 2012}; ^{Norman et al., 2019}; ^{Norman & Niraula, 2016}; ^{Veith et al., 2010}).

Evaluación del modelo

Para evaluar el rendimiento del modelo hidrológico se emplean métodos estadísticos utilizados en las aplicaciones hidrológicas; esto con la finalidad de obtener el grado de semejanza entre los datos observados y los generados. En el presente estudio, se empleó el coeficiente de determinación (R²), el cual describe la proporción de la varianza de los datos observados y simulados. R² varía de 0 a 1, y se consideran aceptables valores mayores a 0.5 (^{Arnold et al., 2012}). Por su parte, el coeficiente de correlación de Pearson (r) presenta variaciones de 0 a 1 (Ecuación 1), donde 0 indica que no hay correlación y 1 representa una correlación perfecta (^{Abbaspour et al., 2015}; ^{Arnold et al., 2012}; ^{Moriasi et al., 2007}).

donde Q _obs es el flujo observado (m³·s^-1), Q _sim es el flujo simulado (m³·s^-1), Q _obs,m es la media del flujo observado (m³·s^-1) y Q _sim,m es la media del flujo simulado (m³·s^-1).

De igual manera, se evaluó el coeficiente de eficiencia de Nash-Suttcliffe (NSE, por sus siglas en inglés), el cual varía de -∞ a 1. Los valores entre 0 y 1 indican una buena predicción por parte del modelo sobre los valores observados (^{Nash & Sutcliffe, 1970}). ^{Moriasi et al. (2007)} clasifican el NSE de la siguiente manera: < 0.2 insuficiente, 0.2 - 0.4 satisfactorio, 0.4 - 0.6 bueno, 0.6 - 0.8 muy bueno y > 0.8 excelente (Ecuación 2).

Como último indicador estadístico se consideró el sesgo porcentual (PBIAS, por sus siglas en inglés) como medida cuantitativa (^{Gupta et al., 1999}), el cual evalúa si los valores simulados por el modelo fueron mayores o menores que los valores observados en el punto de aforo (Ecuación 3). El valor óptimo es el 0, y los valores cercanos a éste representan simulaciones con errores mínimos en el balance de masa; un valor positivo de PBIAS indica que el modelo está subestimando los valores medidos, mientras que un valor negativo indica que el modelo sobreestima lo medido en el punto de aforo (^{Abbaspour et al., 2015}; ^{Moriasi et al., 2007}).

El rendimiento del modelo también se verificó a partir de un gráfico de regresión lineal y una función de identidad (1:1). Para ello, se comparó el flujo simulado contra el observado en los periodos de calibración y validación, ubicando la línea de identidad como referencia para observar la tendencia de los puntos de dispersión.

Escenarios de cambio climático y configuración del modelo

Con la propuesta de forzar el modelo SWAT ante escenarios de cambio climático, se generaron proyecciones mediante la interfaz LARS-WG 6.0, la cual incorpora predicciones de 18 modelos climáticos globales del conjunto CMIP5, utilizados en el quinto informe de evaluación del IPCC (^{Semenov & Stratonovitch, 2010}; ^{Semenov & Stratonovitch, 2015}). Además, se utilizó el modelo HadGEM2-ES (Hadley Center Global Environment Model ver. 2) para proyectar las trayectorias de concentración representativas durante tres periodos de tiempo: 2021-2040, 2041-2060 y 2061-2080 (RCP 2.6, RCP 4.5 y RCP 8.5, respectivamente) (^{IPCC, 2014}). Estos escenarios de cambio climático utilizaron como entrada la temperatura y la precipitación para el modelo SWAT, previamente calibrado y validado, y la simulación del escurrimiento superficial del río Yaqui en períodos de tiempo futuros.

Calibración y validación

Los resultados obtenidos en las tres subcuencas (La Angostura, El Novillo y El Oviachic) muestran que los parámetros estadísticos de rendimiento, para la calibración y la validación de los modelos, están en rangos aceptables (> 0.5); esto de acuerdo con los criterios establecidos por ^{Abbaspour et al. (2015)} y ^{Moriasi et al. (2007)}. Así, los parámetros estadísticos prueban que existe una buena relación entre el flujo observado y simulado después de ajustar los parámetros más sensibles del modelo (^{Arnold et al., 2012}). Con respecto a la evaluación del modelo hidrológico, el Cuadro 3 presenta los resultados obtenidos de acuerdo con el periodo de calibración y validación.

Como se observa en la Figura 4, los hidrogramas de calibración y validación de las tres subcuencas muestran que el flujo observado está ligeramente por encima de las estimaciones simuladas, aunque en algunos eventos es lo contrario. ^{Fontes y Montenegro (2019)} señalan que estas diferencias son resultado de los principales eventos de precipitación que preceden a los impulsos de escurrimiento superficial. Lo anterior debido a que las cuencas con régimen árido son intermitentes en los meses húmedos y suelen presentar eventos de flujos efímeros en los meses secos, cuando no hay escurrimiento.

Figura 4 Escurrimiento observado vs simulado por SWAT para la calibración y validación de tres subcuencas: a) Angostura, b) El Novillo y c) El Oviachic. 

Como se puede observar en la Figura 5a y 5b, la línea de regresión se encuentra por debajo de la línea de identidad, con lo que se verifica que el modelo subestima los datos observados en las subcuencas La Angostura y El Novillo. Por su parte, la línea de regresión en El Oviachic se encuentra por encima a la línea de identidad, lo que evidencia que el modelo sobrestima los datos observados en ambos periodos (Figura 5c). En esta figura también se puede observar que el modelo SWAT tiende a ser más preciso en flujos pequeños que en flujos grandes (mayores a 20 m³·s^-1), de ahí la menor dispersión en flujos menores.

Figura 5 Gráfica escalada del flujo observado y simulado para la calibración y la validación: a) La Angostura, b) El Novillo y c) El Oviachic. 

Es evidente que durante algunos años en la simulación con SWAT ocurrieron algunos fenómenos extremos (sequias e inundaciones) en la subcuenca, por lo que el modelo no obtuvo una buena respuesta de simulación. En este sentido, ^{de Oliveira-Serrão et al. (2019)} mencionan que el modelo no reproduce con precisión algunos periodos de flujo, principalmente los resultantes de fenómenos atípicos.

Escenarios de cambio climático

En La Angostura, el escurrimiento superficial y la precipitación mantuvieron una relación lineal (Figura 6a): R² = 0.80 (RCP 2.6), R² = 0.73 (RCP 4.5) y R² = 0.66 (RCP 8.5), con una disminución en el escenario RCP 8.5. En El Novillo, los valores de R² fueron de 0.78, 0.79 y 0.73, para RCP 2.6, RCP 4.5 y RCP 8.5, respectivamente, mostrándose una reducción del escurrimiento superficial en el último nivel (Figura 6b). Finalmente, en El Oviachic se obtuvieron valores de R² de 0.69, 0.68 y 0.65, para RCP 2.6, RCP 4.5 y RCP 8.5, respectivamente, observándose una diminución en el escenario 8.5 (Figura 6c).

Figura 6 Relación lineal del escurrimiento observado vs precipitación de tres subcuencas: a) La Angostura, b) El Novillo y c) El Oviachic. 

Porcentaje de cambio en el escurrimiento superficial

En las tres subcuencas estudiadas se observaron cambios en el escurrimiento superficial. En La Angostura, considerando el promedio anual, el cambio negativo con mayor impacto fue de -67 % para el periodo 2021-2040 (RCP 2.6), y se espera un aumento de 60 % durante el periodo 2041-2060 (RCP 8.5). En El Novillo se predice una disminución promedio anual de hasta -65 % para el periodo 2061-2080 (RCP 4.5) y un aumento de 42 % durante el periodo 2061-2080 (RCP 2.6), y en El Oviachic se espera una disminución de flujo de -63 % para el periodo 2041-2060 (RCP 2.6 y 4.5) y un aumento de 47 % en el periodo 2041-2060 (RCP 4.5).

Estos resultados contrastan con los observados por ^{Boru et al. (2019)} al identificar un cambio porcentual en el caudal del río Nilo de -24.3 a 49.5 % en el periodo 2021-2040, y de -29 a 20 % en el periodo 2081-2100 para escenarios RCP 4.5. En el mismo estudio, el cambio proyectado en el caudal fue de -18 a 31.9 % en el periodo 2021-2040 y -33 a 22.3 % en el periodo 2081-2100 para RCP 8.5.

Como se observa en el Cuadro 4, en las tres subcuencas la disminución porcentual del flujo fue mayor que el aumento de éste, por lo que se espera un periodo seco más prolongado, principalmente en El Oviachic, seguido de La Angostura y El Novillo (^{Sánchez-Cohen et al., 2020}). ^{Chen et al. (2017)} mencionan que la frecuencia de ocurrencia y la fuerza de las sequías en la estación seca aumentará por el calentamiento global, mientras que las inundaciones en la estación húmeda tienden a aumentar. Específicamente en los embalses de la cuenca del río Yaqui, es importante considerar el volumen total escurrido en las tres subcuencas, ya que el distrito de riego efectúa transvases para proveer de agua al distrito de riego número 41 y a las colonias Yaquis, aguas debajo de las presas. Asimismo, las presas proveen de agua para la generación de energía.

Respecto de la representatividad del modelo al medio físico, el modelo SWAT en México ha sido usado en algunas cuencas como herramienta de decisión para su manejo. ^{Cruz-Arévalo et al. (2021)} usaron este modelo en la cuenca del río Chapingo para evaluar el impacto del uso del suelo en escurrimiento. Los resultados de estos autores muestran la necesidad de contar con datos observados de las variables del medio físico para una adecuada calibración. De igual manera, ^{Torres-Benites et al. (2004)} calibraron el modelo SWAT en una cuenca del Estado de México (El Tejocote) para evaluar su capacidad predictiva y concluyeron que el modelo subestimó en 18 % el escurrimiento y en 38 % los sedimentos. Otros estudios aseguran que el modelo SWAT subestima el factor pendiente en las subcuencas cuando este valor es mayor al 25 % (^{Rivera-Toral et al., 2012}).

Existen diferentes estudios que concuerdan con los resultados del presente análisis. ^{Chang et al. (2017)} reportan cambios similares en el escurrimiento de una cuenca en China en la que evaluaron diferentes escenarios climáticos y sus variaciones en precipitación. En las regiones Alpinas, ^{Braun et al (2020)} estudiaron el efecto del cambio en los patrones de escurrimiento causado por el incremento en temperaturas y observaron diferentes escenarios en la producción de escurrimientos.

El uso del modelo SWAT para predecir los posibles impactos producidos por el cambio climático hacia los recursos hídricos se ha recomendado como apoyo para la toma de decisiones en diferentes partes del mundo (^{Abbaspour et al., 2015}). Por ello, se considera una herramienta importante para la predicción de posibles escenarios futuros en la cuenca del río Yaqui.