Descripción de la cuenca

La cuenca del río Fuerte se localiza en el Noroeste de México, en los estados de Sinaloa, Sonora, Durango y Chihuahua (Figura 1), como parte de la región hidrológica No. 10 Sinaloa, entre 25.68° y 28.24° N, y 106.12° a 109.43° O, con un área de 36 456 km². Nace en la Sierra Madre Occidental y desemboca al norte de Sinaloa en el Golfo de California; en la parte alta de la cuenca está la etnia indígena Rarámuri y es una zona de alta marginación que en la última década tuvo graves problemas debido a la sequía, y en la parte media y baja está el Valle del Fuerte que es una importante zona agrícola de producción de granos y hortalizas para el mercado nacional e internacional. Las altitudes en la cuenca van de 3168 a -9 msnm (^{INEGI, 2014}); su precipitación media anual es de 691 mm y temperatura media de 19.4 °C, de acuerdo con las normales climatológicas de 14 estaciones dentro de la cuenca de 1961 a 2011 (^{SMN, 2014}). La cuenca del río Fuerte es importante por su extensión superficial, el volumen de sus escurrimientos y las presas Miguel Hidalgo y Costilla, Luis Donaldo Colosio y Josefa Ortiz de Domínguez, que se usan para la generar energía eléctrica y son la principal fuente de agua para uso agrícola, urbano e industrial.

Información climatológica

La información climatológica usada para el cálculo de los índices de sequía SPI y SPEI se obtuvo de la base de datos del servicio meteorológico nacional disponible en línea (^{SMN, 2014}), se emplearon los datos mensuales de precipitación (Pt), temperatura mínima (Tmin) y máxima (Tmax) de 14 estaciones meteorológicas durante el periodo 1961-2011, el cual fue el más completo y actualizado durante este estudio. Las estaciones seleccionadas están distribuidas por toda la cuenca; la información más incom-pleta fue la de precipitación y en promedio, se generó el 8 % de datos faltantes de lluvia (Cuadro 1). En uno de los casos más extremos, en Batopilas, Chihuahua, se generó 29 % de los datos de lluvia. Para generar los datos de lluvia se usó el método del inverso del cuadrado de la distancia, y los detalles están en un artículo en el cual se generaron los índices (^{Castillo et al., 2017}).

Metodología

La sequía meteorológica, hidrológica y agrícola están representadas por índice de sequía que es una variable para evaluar el efecto de una sequía e incluye intensidad, duración severidad y extensión espacial (^{Mishra and Singh, 2010}). El SPI (^{McKee et al., 1993}) es uno de los más usados para la detección y monitoreo de sequías, y se bas en la probabilidad de precipitación a cualquier escala de tiempo para cuantificar el déficit de precipitación (^{Velasco et al., 2004}). El SPEI (^{Vicente et al., 2010}), además de la precipitación que usa el SPI, considera la evapotranspiración de referencia (ET₀) y se basa en el balance hídrico (Pt-ET₀) mensual, combinando características de la sensibilidad del Palmer Drought Severity Index (PDSI) y de la simplicidad de cálculo y la naturaleza multitemporal del SPI. Es muy adecuado para detectar, monitorear y explorar las consecuencias del calentamiento global en condiciones de sequía (^{Vicente et al. 2010}).

Los índices SPI y SPEI se desarrollaron para las 14 estaciones seleccionadas, a escalas de 3, 6, 12 y 24 meses. Las series del SPI se generaron con el programa spi_sl_6.exe desarrollado por el National Drought Mitigation Center (^{NDMC, 2014}), y el SPEI se obtuvo con SPEI.R desarrollado por ^{Beguería y Vicente (2014)} para el programa R/RStudio. En el cálculo del SPEI se consideró el método de evapotranspiración de referencia de Hargreaves-Samani (^{Beguería et al., 2014}). Para mayor información de la obtención de SPI y SPEI, y su interpretación como indicadores de sequía, consultar a ^{Castillo et al. (2017)}, quienes detallan el proceso de SPI y SPEI en la cuenca para nuestro estudio.

En nuestro estudio el pronóstico de sequías meteorológicas se realizó mezclando modelos de series de tiempo autorregresivos con el filtro de Kalman discreto (^{Kalman, 1960}), que es un conjunto de ecuaciones matemáticas para estimar el estado de un proceso minimizando la media del error cuadrático. El filtro de Kalman opera por medio de un mecanismo de predicción y corrección, el algoritmo pronostica el nuevo estado a partir de una estimación previa, añadiendo un término de corrección proporcional al error de predicción, minimizándolo estadísticamente (^{Welch y Bishop, 2006}).

Los modelos autorregresivos se usan en la hidrología y meteorología porque tienen dependencia del tiempo y son fáciles de usar (^{World Meteorological Organization, 2011}). En nuestro estudio se propuso la creación de dos modelos autorregresivos, uno de segundo orden (AR2) y otro de segundo orden con entrada exógena (ARX) para pronosticar los estados mensuales futuros de los índices de sequía SPI y SPEI con base en los registros anteriores de la serie; la elección de modelos autorregresivos se demostrará en el Cuadro 2 de Resultados y Discusión. Los modelos AR2 y ARX (ecuaciones 1 y 2, respectivamente) relacionan la entrada del sistema con su salida mediante una ecuación lineal en diferencias con coeficientes constantes:

donde y_t es el valor observado del índice en el tiempo t, el cual representa un mes; r_t es el valor de la variable exógena (Pt, Tmin, Tmax, ET₀) en el tiempo t; e_t+1 es el término de error en la estimación del índice; α_i y β_i son parámetros; los índices na y nb especifican el número de observaciones previas del índice y de las variables exógenas, respectivamente. A partir de la estructura general de ambos modelos autorregresivos, la formulación en espacio de estados permite utilizarlos dentro del algoritmo del filtro de Kalman Discreto así:

donde x_k+1 es el valor del índice (no observado) de tamaño (n x 1); A es la matriz de parámetros α_i de tamaño (n x n); x_k es el índice en el tiempo k de tamaño (n x 1); B es la matriz de parámetros exógenos β_j de tamaño (n x m); v_k es el vector que contiene la variable exógena registrada para el tiempo k; z_k es el índice en el tiempo k de tamaño (m x 1); H es la matriz de transformación que mapea el vector de estados al dominio de la medición con dimensiones (m x n); W_k y v_k son vectores que representan el ruido gaussiano en el proceso y el ruido en la medición para cada observación con tamaños (m x 1), y lo esperado es que tales ruidos gaussianos se distribuyan de manera normal con media 0 y varianzas Q y R, respectivamente:

De acuerdo con ^{Simon (2001)}, se asume que no existe correlación entre Q y R, es decir, son variables aleatorias independientes y pueden variar en el tiempo, pero se suponen constantes por simplicidad, y se pueden definir como:

donde Q es la matriz de covarianza de la perturbación del sistema; R es la matriz de covarianza de la perturbación de la medición; wkT y vkT indican la transpuesta de w_k y v_k , respectivamente, y E[.]significa el valor esperado. La matriz Q es de dimensión (na x na) y contiene el valor ajustado del MSE (Media Cuadrada del Error) estimado por el modelo autorregresivo en la etapa de calibración.

La matriz R es de (1 x 1) y contiene el error esperado de las mediciones, definido como una proporción (α) de la medición anterior, entonces R=α*(x_k -1), en este caso específico α=0.001, valores mayores generan datos extremos que se ubican fuera del rango en que normalmente se presentan los valores de los índices de sequía, reduciendo el poder de predicción del modelo. Con respecto al tamaño de las matrices, el valor de n es igual a na, y el valor de m igual a nb; donde na y nb indican el número de observaciones previas del índice y de las variables exógenas, respectivamente. El algoritmo del filtro de Kalman Discreto se describe en la Figura 2.

Para crear los modelos AR y ARX se programaron rutinas en Matlab^® (^{Math Works 2015}), para obtener los valores de las matrices de parámetros en función del orden del proceso autorregresivo de acuerdo a las ecuaciones 1 y 2. El ARX también quedó como modelo autorregresivo de segundo orden pero con variable exógena.

Debido a que el comportamiento de la sequía en la cuenca varía en función de la ocurrencia de la temporada de lluvias y de fenómenos como El Niño o la Niña, entre otros fenómenos climáticos globales, los parámetros estimados para los modelos AR y ARX se recalcularon cada cierto periodo (P) para obtener pronósticos más representativos de las condiciones previas de humedad en cada estación, el pronóstico se inició en el periodo (t₀ + P) hasta el tiempo (t₀ + 2P), en este punto se recalcularon los parámetros de los modelos en función de periodo [P : 2P], así sucesivamente hasta terminar con el último grupo de datos. De esta forma la implementación del DKF - AR2 y del DKF - ARX se hizo dinámica e incorporó los cambios climáticos que se produjeron en la cuenca durante el periodo de estudio cada temporada o cierta época.

El pronóstico de los índices SPI y SPEI se realizó a cuatro escalas temporales para cada estación, y se implementaron los modelos AR2 y ARX (Pt, Tmin, Tmax y ET₀ como variables exógenas) con el filtro de Kalman Discreto. Para conocer la efectividad de los modelos en el pronóstico para L pasos hacia adelante, es decir, con 1, 2, 3 y 4 meses de adelanto, se creó una rutina en Matlab^® para realizar el pronóstico considerando la información en el tiempo k y avanzar L pasos en el pronóstico sin actualizar la información (Figura 3), este proceso se realizó considerando el periodo P de calibración de los modelos AR2 y ARX.

Para evaluar los resultados del pronóstico se calcularon las principales estadísticas de acuerdo con ^{Gupta et al. (2009)}, el RMSE (la raíz del error cuadrático medio) y E (la eficiencia de Nash-Sutcliffe) fueron los criterios a considerar por ser los más usados para la calibración y evaluación de modelos hidrológicos (^{Moriasi et al., 2007}). Además, se incluyó un intervalo de confianza al 95 % que sigue el recorrido de cada estado de la serie con el objetivo de establecer el grado de incertidumbre asociado al pronóstico en cada paso del tiempo (L) con respecto a los valores observados.

El intervalo de predicción (I.P.) se obtuvo de acuerdo con los criterios establecidos por ^{Chatfield (2004)}, el I.P. se calculó con base en la siguiente forma general: Un 100(1-α)% I.P. para x_t+L está dado por:

donde etL=xt+L-x^t(L) denota el error del pronóstico al tiempo t cuando se pronostica con L periodos de adelanto; z_a/2 denota el punto porcentual de una distribución normal estándar con una proporción de α/2 por sobre de él.

Pronóstico con DKF - ARX-Pt

En las estaciones 8167, 25009, 25019, 25025, 25042, 25044, 25100 y 26053 (6 de ellas están en Sinaloa), se obtuvieron mejores resultados en el pronóstico del SPEI con el modelo DKF - ARX-Pt. Los resultados mostraron que el empleo de la precipitación como variable exógena mejora el pronóstico en comparación con el realizado con el modelo DKF-AR2 para estas estaciones. Con el modelo ARX se probaron como variables exógenas la precipitación, la evapotranspiración de referencia, la temperatura máxima y la mínima; sin embargo, en la mayoría de los casos la precipitación aportó mayor información al algoritmo y generó mejores resultados en términos de E y R que las variables restantes (Cuadro 6).

La asociación de los resultados con las unidades climáticas de la cuenca ubica estas estaciones en climas cálidos y secos. Como ya se mostró, las estaciones de clima cálido presentan series con mayor oscilación entre valores negativos y positivos del SPEI que las series de estaciones en climas templados y semifríos. Al tener series menos suavizadas el pronóstico se hace menos preciso, por lo que la inclusión de una variable exógena (precipitación) aporta la información necesaria al modelo para mejorar su desempeño.

Las Figuras 8 y 9 presentan los diagramas de dispersión del SPEI de 12 y 24 meses, respectivamente, para las estaciones donde el modelo ARX-Pt resultó mejor. Estos muestran los valores de MSE, RMSE, E, R y PBE.

Pronóstico con DKF - AR2

Las estaciones 8038, 8106, 8161, 8172, 8182 y 8267, toda ubicadas en el estado de Chihuahua sobre la Sierra Madre Oriental, obtuvieron mejores resultados en el pronóstico del SPEI con el modelo DKF - AR2 (Cuadro 6). Los resultados mostraron que en estas estaciones la inclusión de una variable exógena como la precipitación, la evapotranspiración de referencia, la temperatura máxima y la mínima no aporta la suficiente información al algoritmo para mejorar los resultados obtenidos con el modelo AR2, en el mejor caso los iguala. La asociación de los resultados con el tipo de clima de la estación ubica estas estaciones en climas semifríos y templados semihúmedos.

El valor del estadístico PBE (Percentage Bias Error, Porcentaje de Error de Sesgo) presenta valor negativo en la mayoría de las estaciones, lo que indica que el modelo DKF-AR2 subestima el valor del índice SPEI (Cuadro 6).

En ambos modelos de pronóstico del SPEI, DKF-AR y DKF-ARX, las series de 12 meses presentan mayor error y menor valor del coeficiente de Nash que las de 24 meses, esto se debe a que a una escala mayor la serie temporal asimila más lentamente los cambios dados en el balance Pt-ET₀ a lo largo del tiempo, y, por lo tanto, presenta mayor suavidad en su recorrido mejorando el pronóstico. Sin embargo, ambos pronósticos de sequías de duraciones de 12 y 24 meses son útiles. Pero dadas las características de los índices de sequía (SPI y SPEI) a una determinada escala (duración de sequía), en términos de operatividad, la escala temporal de 12 meses se considera la más adecuada para el pronóstico de las sequías. Esto último, es especialmente apropiado en la parte media y baja de la cuenca, por su vocación agrícola y energética, donde al operar las presas se revisan políticas anuales de operación. Aunque el análisis de periodos de sequía de 24 meses es muy útil también.

Además del pronóstico del índice de sequía para el tiempo t+1|t, se pronosticaron los índices con 2, 3 y 4 meses de adelanto, es decir, dado un valor observado del índice en el tiempo t se pronosticaron los valores del índice para el tiempo t+1, t+2, t+3 y t+4, sin actualizar la información, esto con el propósito de conocer el poder de predicción de los modelos para diferentes pasos de adelanto en el tiempo.

Sin embargo, a manera de muestra, aquí se presentan solo los resultados del pronóstico para el SPEI de 12 meses. La Figura 10 muestra el pronóstico en la estación de la Presa Miguel Hidalgo en Sinaloa, con el SPEI de 12 meses con el modelo DKF - ARX-Pt con 1, 2, 3 y 4 meses de anticipación.

El pronóstico se concentra en la estimación puntual de un valor del SPEI de 12 meses para algún mes futuro, pero de acuerdo con ^{Chatfield (2004)}: “el pronóstico puntual es adecuado para muchos propósitos, pero un intervalo de predicción es frecuentemente de gran ayuda para dar un mejor indicador de la incertidumbre futura”. Por lo tanto, se calculó el intervalo de confianza del pronóstico al 95 % para los pasos de adelanto (L) como se indica en la ecuación (9).

En la Figura 10, además de los valores observados y pronosticados, hay una franja que sigue el recorrido del pronóstico (entre más negativo el índice, mayor es la sequía), y este es el intervalo de confianza del pronóstico al 95 %. Conforme se avanza en el pronóstico el intervalo se hace más ancho, esto es, por ejemplo, que en el tiempo t hay menor incertidumbre de donde puede localizarse el valor observado del SPEI en el tiempo t+1 que del valor en el tiempo t+4 en el cual la incertidumbre es mucho mayor. También es importante agregar que al final se observó como los errores totales de la predicción (del proceso + de medición) se ajustaron a una distribución t-Student de media 0, la cual se acerca a una distribución Normal.

Es importante enfatizar que el filtro de Kalman es importante para el pronóstico de variables hidrometeorológicas porque se pueden predecir caudales medios diarios (^{Gonzalez et al, 2015}), caudales subhorarios (^{Morales et al., 2014}), y predecir indicadores de sequías con algunos meses de anticipación. La predicción de sequías permite preparar planes de contingencia para reducir su impacto negativo.