Introducción

La predicción de eventos de lluvia mediante análisis de frecuencias es de suma importancia en diversos estudios de ingeniería, como el diseño de obras hidráulicas para el control de escurrimientos, la planeación adecuada de uso de suelo y la evaluación de la disponibilidad de agua. No obstante, debido a las implicaciones económicas y sociales, el análisis de frecuencias debe realizarse con alta precisión, pues la subestimación de la disponibilidad de agua podría implicar fallas en las estrategias de asignación de agua a los usuarios, en particular durante periodos prolongados de sequías.

La metodología clásica de análisis y predicción de eventos de variables hidrológicas se basa en la teoría de valores extremos de series estacionarias. Bajo condiciones de estacionarei-dad, la distribución de la variable de interés es invariante en el tiempo, sin tendencias, cambios y periodicidades (^{Villarini & Smith, 2010}). No obstante, la validez de la estacionareidad de las series de tiempo debido al cambio del clima está siendo cuestionada por la comunidad científica (^{Milly et al., 2008}; ^{Gleick, 1989}; ^{Voss, May & Roeckner, 2002}; ^{Held & Soden, 2006}; ^{Webster, Holland, Curry, & Chang, 2005}). Según el Grupo Intergubernamental de Expertos sobre Cambio Climático (^{IPCC, 2007}), las actividades humanas están provocando el calentamiento del planeta, por lo que se han registrado numerosos cambios de largo plazo en el clima, que incluyen posibles aumentos en la intensidad de los ciclones tropicales y las ondas de calor, así como en la intensidad y frecuencia de eventos extremos, como sequías y lluvias torrenciales. Se espera que el cambio climático intensifique el estrés actual sobre los recursos hídricos debido al crecimiento de la población, a los cambios económico y de uso de la tierra, y en particular a la urbanización (^{IPCC, 2007}). Por lo anterior, conocer la disponibilidad del agua es un factor clave en la creación de estrategias para afrontar periodos de escasez, no obstante, el proceso de asignación de agua a los usuarios podría ser ineficiente si la disponibilidad se evalúa a través de valores medios de escurrimientos o precipitación, sin considerar el cambio en el comportamiento temporal de estas variables. Por el contrario, si la disponibilidad se estima a partir de eventos de lluvia asociados con una determinada probabilidad y escenarios que consideren la disminución de la lluvia, se podrán crear mejores estrategias para la distribución del agua.

En México, los escenarios climáticos construidos hasta ahora sugieren que la precipitación en muchas regiones del país podría disminuir; por lo tanto, bajo un escenario de cambio climático con condiciones de sequía prolongada, el desarrollo sostenido de la sociedad estaría en peligro (^{Magaña-Rueda & Gay-García, 2002}); principalmente en zonas del noroeste del país, donde existe una alta presión sobre los recursos hídricos (^{Conagua, 2012}). Ante este panorama, el objetivo de este trabajo fue estimar eventos de lluvia acumulada anual del noroeste de la república mexicana, asociados con diferentes periodos de retorno y escenarios futuros, a través de modelos probabilísticos válidos para series de tiempo no estacionarias. Los resultados obtenidos podrán tomarse como base para una adecuada evaluación de la disponibilidad de agua en la región.

Materiales y métodos

La zona de estudio abarcó los estados de Baja California Norte, Baja California Sur, Sonora, Sinaloa, Chihuahua y Durango, ubicados en el noroeste de la república mexicana (Figura 1). Estas entidades fueron seleccionadas por tener altas presiones sobre los recursos hídricos (^{Conagua, 2012}) y zonas con alta vulnerabilidad ante las sequías.

Los datos pluviométricos utilizados en este trabajo corresponden a los registros diarios de las estaciones climatológicas convencionales ubicadas en la zona de estudio, contenidas en la base de datos CLICOM (Climate Computing Project) del Servicio Meteorológico Nacional de la Comisión Nacional del Agua (Congua). Se consideraron estaciones con al menos 90% de registros completos entre el periodo de 1950 a 2013. Mediante la prueba de ^{Grubbs y Beck (1972)} se identificaron valores de lluvia anual que se desvían de modo considerable de la distribución de los datos, y fueron comparados con los registros de estaciones aledañas para corroborar su magnitud. Los valores que resultaron claramente erróneos fueron eliminados. Para deducir datos faltantes se aplicó el método de interpolación de la distancia inversa ponderada (^{Shepard, 1968}), en el cual se consideraron dos estaciones de apoyo y un exponente de distancia igual a dos. Para verificar si existe mejora en la estimación de la media y la varianza de las muestras extendidas, se utilizó el criterio de información relativa (I), definido por ^{Escalante-Sandoval y Reyes-Chávez (2002)} como:

Donde VarSy12 es la varianza de la varianza de la serie original, y VarSy12 es la varianza de la varianza de la serie extendida.

Si I > 1, entonces la varianza de la varianza de la serie extendida no excede a la original, por lo que la extensión de los registros es adecuada. Con base en lo anterior, se obtuvieron datos confiables de lluvia diaria de 221 estaciones.

Para poder comparar espacialmente los cambios ocurridos en las series de lluvia, se llevó a cabo una regionalización de la lluvia mensual basada en un análisis de componentes principales (ACP) en modo S (^{Richman, 1986}). Se seleccionaron las primeras cinco componentes que explican 73% de la varianza de la precipitación. Las componentes retenidas se sometieron a una rotación oblicua directa, con el parámetro de oblicuidad δ = 0. En la delimitación de las regiones homogéneas se asignó a cada estación la componente principal de mayor carga en valor absoluto, considerando sólo aquellas con valores mayores a 0.4. Las regiones se definieron de acuerdo con el área de influencia de cada estación obtenida con el método de polígonos de Thiessen.

El siguiente paso fue verificar la presencia de cambios graduales en las series (tendencias) o puntos de cambio abruptos en la media y/o la varianza de la distribución de la variable de interés (^{Villarini et al., 2009}). Se aplicó la prueba de Pettitt (^{Pettitt, 1979}), que permite detectar cambios significativos en la media (mediana) de las series cuando el punto de cambio es desconocido. Las tendencias en las series se analizaron con el método de descomposición de modos empíricos (DME) desarrollado por ^{Huang et al. (1998)}, el cual es un método intuitivo, directo y adaptativo a procesos no lineales y no estacionarios, cuya efectividad para el análisis de tendencias y señales de periodo largo se ha demostrado en estudios recientes (^{Qian, Wu, Fu, & Zhou, 2010}; ^{Carmona, 2010}; ^{Wu, Huang, Long, & Peng, 2007}). Para determinar la significancia de las tendencias, los residuales obtenidos por DME se analizaron con la prueba de Mann-Kendall (^{Kendall, 1975}).

De los fenómenos de gran escala que influyen en la variabilidad climática en México, se seleccionó la oscilación decadal del Pacífico (PDO, por sus siglas en inglés) para ser incluida en el análisis de frecuencias de la lluvia, esto debido a la asociación existente entre episodios persistentes de valores PDO negativos y condiciones de sequías severas y prolongadas a través del occidente y centro de América del Norte (^{Cook, Woodhouse, Eakin, Meko, & Stahle, 2004}). Por tal razón, se analizó la correlación entre las series anuales regionales de lluvia acumulada y el índice PDO, a través del coeficiente de correlación ρ (rho), propuesto por ^{Spearman (1904)}, el cual es una prueba no paramétrica que mide la asociación entre dos variables discretas y está definido por la ecuación (2):

Donde los valores de ρ varían en el intervalo [-1, 1]; d es la diferencia entre los valores correspondientes de x-y; N, el número de parejas, para N > 20, sigue una aproximación a la distribución t de Student (3):

En el análisis de frecuencias, para considerar la falta de estacionareidad de la series se pueden asociar los parámetros de las distribuciones con más de una covariable. La inclusión de covariables ha permitido incorporar tendencias o índices de circulación atmosférica de gran escala como ENSO (El Niño-Southern Oscillation), PDO (Pacific Decadal Oscillation o NAO (North Atlantic Oscillation), en la modelación de datos de precipitación extrema (^{Katz, Parlange, & Naveau, 2002}; ^{El Adlouni, Ouarda, Zhang, Roy, & Bobée, 2007}; ^{Villarini, Smith, & Napolitano, 2010}; ^{Vovoras & Tsokos, 2009}; ^{Ouarda & El Adlouni, 2011}), temperaturas máximas (^{Villarini et al., 2010}), gastos máximos (^{Villarini et al., 2009}; ^{Towler et al., 2010}; ^{Vogel, Yaindl, & Walter, 2011}), nivel extremo del mar (^{Coles, 2001}) e intensidades de huracanes (^{Mestre & Hallegatte, 2009}).

No existe una teoría general establecida para procesos no estacionarios, por lo que es usual utilizar como base modelos de valores extremos estándares y extenderlos (^{Coles, 2001}). Por ejemplo, de la distribución log normal (LN3) con parámetros x₀, μ y σ (ubicación, escala y forma, respectivamente), se puede obtener un modelo adecuado para estimar una variable x_t en función del tiempo, de la siguiente forma:

Donde cada parámetro x₀(t), μ(t), σ(t) tienen una expresión en términos del tiempo. El parámetro μ(t) se puede expresar con la ecuación (5):

Donde t corresponde a la tasa anual de cambio en el tiempo. De esta forma, variaciones a través del tiempo en los procesos observados se modelan como una tendencia lineal en el parámetro de ubicación de la distribución.

La estimación de los parámetros de un modelo en función del tiempo, o de otra covariable, se puede realizar mediante el método de máxima verosimilitud (^{Katz et al., 2002}), considerado como el más eficiente, pues proporciona la menor varianza muestral de los parámetros estimados y, por ende, de los eventos estimados, en comparación con otros métodos; además, una ventaja es su adaptabilidad a cambios en la estructura del modelo.

Para un modelo no estacionario propuesto de la función LN3, la función de verosimilitud es la siguiente:

Donde x₀(t), μ(t) y σ(t) son los parámetros de ubicación, escala y forma, respectivamente, en función del tiempo o de otra covariable.

En el análisis de lluvia acumulada anual del noroeste de México es de interés conocer cuáles serán los escenarios más adversos de escasez de agua en el futuro. Por lo tanto, la predicción de eventos de lluvia se debe realizar con funciones de probabilidad para valores mínimos. Sin embargo, no existen modelos probabilísticos no estacionarios para funciones de distribución de valores mínimos. Por tal razón, se propusieron los modelos descritos en el Cuadro 1. Los modelos M₀ corresponden al análisis de frecuencias convencional, mientras que M₁, M₂ y M₃ corresponden al caso no estacionario y sus parámetros se expresaron en términos del tiempo y del índice PDO. La serie promedio anual del índice PDO se obtuvo a partir de los valores mensuales del periodo de 1950 a 2013, descargados de la página del National Climatic Data Center de los Estados Unidos de América (^{NCDC, 2013}).

Para cada estación, los estimadores de los parámetros de las funciones se obtuvieron tras maximizar las funciones logarítmicas de verosimilitud respectivas. Los modelos que mejor describen la variabilidad de los datos se seleccionaron mediante el Criterio de Información de ^{Akaike (1974}), en el cual se utilizó la expresión (7). El mejor modelo es el que tuvo el menor valor AIC. En el caso de tener valores similares de AIC, se seleccionó el modelo con menos parámetros:

Donde l es el máximo valor de la función de verosimilitud y K es el número de parámetros estimados.

Se utilizó el método descrito por ^{Coles (2001}) para comparar la validez de un modelo M₁ contra otro M₀, tal que M₀ ⊂ M₁, en el cual se utiliza la medida de discordancia definida por la expresión (8).

Donde l(M_i) es el máximo valor de la función logarítmica de verosimilitud del modelo M_i .

El estadístico D se distribuye de acuerdo con la distribución chi-cuadrada Xν2. El parámetro ν es la diferencia entre el número de parámetros de los modelos M₁ y M₀. Valores de D más grandes que las cantidades de la distribución Xν2 para un nivel de confianza particular se consideran significantes; entonces se rechaza M₀ en favor del modelo M₁.

Una vez seleccionado el mejor modelo fue necesario verificar que se ajustara bien a los datos observados. Debido a que en una serie no estacionaria los datos no se distribuyen de manera idéntica, éstos deben ser transformados (^{Coles, 2001}; ^{Katz et al., 2002}). La transformación de los datos consistió en la normalización de los mismos, proceso al que ^{Krzysztofowicz (1997)} denominó como transformación normal de cuantiles. La normalización se realizó mediante la expresión (9), en la que una variable aleatoria X, con una determinada distribución de probabilidad F, se transforma en una variable Y, con una distribución normal estándar G:

Donde G^-1 es la función inversa de la distribución normal estándar definida por la ecuación

De esta forma, tras ordenar los valores estandarizados y asociarlos con los correspondientes de la distribución empírica de Weibull, se obtuvieron los siguientes pares de datos:

Con los pares de datos anteriores se elaboraron las gráficas de cuantil-cuantil (Q-Q) y "worm plot" (^{Buuren & Fredriks, 2001}). En una gráfica "worm plot", el eje vertical es la diferencia entre los valores empíricos y teóricos, y contiene el intervalo de confianza a 95% por la ecuación (12):

Donde f(z) es la función de densidad normal; z es un evento asociado con una probabilidad, y n es el tamaño de la muestra.

Finalmente, se estimaron eventos de diseño X^, asociados con diferentes periodos de retorno. Para la función LN3, con el modelo M₃, los eventos de diseño se estimaron al resolver la ecuación (13):

Donde x^₀ es el estimador del parámetro de ubicación; σ^, el estimador del parámetro de forma; μ^(t), el parámetro de escala en función de los estimadores ^{μ^1}₀, ^{μ^1}₀ y ^{μ^1}₂; t₀, el valor del índice del tiempo para un escenario determinado (t₀ = 65 para el escenario del año 2014 y t₀ = 95 para el escenario del año 2044); PDO, un valor promedio del índice PDO, y X es el evento de diseño asociado con un periodo de retorno T y un escenario futuro t₀.

El análisis de datos abarcó el periodo comprendido entre los años 1950 y 2013, por lo cual, al primer año (1950) se le asignó el valor t₀ = 1, mientras que al último (2013), le correspondió un valor t₀ = 64. De esta forma, para años posteriores, los valores del índice de tiempo aumentaron una unidad de forma consecutiva.

Debido a la asociación de valores PDO negativos y periodos con déficit de lluvia, el valor PDO^ se estimó como la media aritmética de los valores PDO correspondientes al lapso de anomalías con mayor déficit de lluvia acumulada anual analizada.

Resultados y discusión

La distribución espacial de la lluvia acumulada anual en la zona de estudio se encuentra en la Figura 2. Los valores más bajos se presentan en la península de California, con menos de 200 mm al año, siendo la zona con menor precipitación de todo el país; en contraste, entre los límites de Sonora, Sinaloa, Chihuahua y Durango se tienen valores de más 1 000 mm.

El análisis de componentes principales arrojó cinco componentes que explican 73% de la varianza de la precipitación. Tras asignar la componente de mayor carga al polígono de Thiessen de cada estación, se pudieron identificar espacialmente grupos de polígonos asociados con la misma componente; cada grupo se consideró como una región homogénea; de esta forma, se generaron las cinco regiones mostradas en la Figura 2.

De acuerdo con el Cuadro 2, 19% de las series analizadas es no estacionaria, de las cuales, 47% se ubica en la región 1, que abarca principalmente los estados de Sonora y Sinaloa. Del total de series que presentan puntos de cambio descendentes significativos, 73% corresponde a la región 1; en esta región, además, se localizan las tres series que presentan puntos de cambio descendentes y tendencias decrecientes, así como dos de las tres series con tendencia decreciente.

Según el análisis de correlación de Spearman, existe una correlación de 0.43, estadísticamente significativa a cualquier nivel, entre el índice PDO y la serie de lluvia acumulada de la región 3; mientras que con las regiones 4 y 5, la correlación es de 0.2, significativa a un nivel de 0.1. Para las tres regiones anteriores, la correlación es positiva; por lo tanto, al estar el PDO en fase negativa de manera persistente, se tendrán periodos con déficit de lluvia, tal como ocurrió en la región 3 entre 1950 y 1975 (Figura 3a). La asociación entre episodios prolongados de valores PDO negativos y condiciones de sequías severas y prolongadas a través del occidente y centro de Norteamérica fue mencionada por ^{Cook et al. (2004)}. Entre 1977 y 1997 se presentó un régimen PDO positivo y aparentemente a partir de 1998 se tiende hacia una fase negativa, por lo que es muy probable que en las próximas décadas cambie de polaridad (^{Mantua, Hare, Zhang, Wallace & Francis, 1997}; ^{MacDonald & Case, 2005}). Debido a lo anterior, es posible que la pérdida de la estacionareidad de las series por el incremento de la lluvia se deba en gran medida a la influencia de la fase positiva del PDO; por lo tanto, la persistencia de valores PDO negativos que pudieran causar sequías prolongadas se consideró en el análisis de predicción de la lluvia.

Del análisis de frecuencias realizado a las 221 series de lluvia acumulada anual, las pruebas de bondad de ajuste indicaron que para 84 de éstas, los modelos no estacionarios estadísticamente representan mejor la variabilidad de los datos (Cuadro 3). En el caso de las series no estacionarias, en 83% de los casos los modelos no estacionarios se ajustó mejor a los datos; para las series no estacionarias restantes, aunque el estadístico AIC indicó que los modelos no estacionarios fueron mejores, la medida de discordancia (D) fue inferior a 3.84 (valor de Xν2 a un nivel ɑ = 0.05), por lo que no se justificó el uso de un modelo M₁, M₂ o M₃ en lugar de un M₀.

La mayoría (60%) de los modelos no estacionarios seleccionados es M₂ o M₃ de la función LN3. Los resultados obtenidos son consistentes con las pruebas de estacionareidad y de correlación realizadas; por ejemplo, para la región 3, que tiene la correlación más alta con PDO, se seleccionaron modelos M₂ en 20 casos. Además, para la región 1, donde se localiza 47% de las estaciones con series no estacionarias, se seleccionaron 22 modelos M₃. Es importante señalar que el modelo M₀ de la distribución GW se seleccionó para 61 de las series estacionarias; por lo tanto, las distribuciones mezcladas representan mejor la variabilidad de los datos cuando éstos son estacionarios.

Como ejemplo del procedimiento para la selección del mejor modelo, en los Cuadros 4 y 5 se muestran los resultados de las pruebas de bondad de ajuste y los estimadores de los parámetros de los modelos ajustados a las series de las estaciones 2060 ("Santa Cruz", Ensenada, Baja California) y 25082 ("San Blas", El Fuerte, Sinaloa), respectivamente. Cabe mencionar que la serie de la estación 2060, según la prueba de Pettitt, presentó en el año 2002 un punto de quiebre significativo a un nivel de 0.05, con un cambio descendente en la media de la serie (Figura 3b). Por otro lado, para la estación 25082, el residual obtenido por el método DME (Figura 3c) indicó que existe una tendencia monótona decreciente, estadísticamente significativa (a un nivel de 0.05), según la prueba de Mann-Kendall.

Para la estación 2060, el modelo con el menor valor AIC fue el M₂ de la función LN3, en el cual el parámetro de escala está en función del índice PDO; además, el valor del estadístico D es superior a 3.84, lo cual indica que es más adecuado y explica mejor la variabilidad de los datos que el modelo M₀ . Para la serie de lluvia de la estación 25082, el modelo M ₁ de la función GA3, que incluyó una tendencia en el parámetro de ubicación, representa mejor la variabilidad de los datos que los modelos restantes; el valor de AIC del modelo M₁ fue el menor, con un valor D igual a 5.0, estadísticamente significativo.

La Figura 4 corresponde a las gráficas "worm plot" y de cuantil-cuantil (Q-Q) de las estaciones 2060 y 25082, para los modelos mencionados. Se observa para ambos casos que los datos de las gráficas caen dentro de los límites de confianza y están cerca de la diagonal unitaria; por tal razón, se determina que los modelos seleccionados se ajustan de forma adecuada a los datos.

Para la estación 2060, el valor del índice PDO utilizado para estimar los eventos de diseño fue de -0.86, promedio del periodo comprendido entre 2007 y 2013, caracterizado por la persistencia de valores negativos y que abarcó el periodo con mayor déficit de lluvia (Figura 3b). Por lo anterior, para el modelo M₂ de la función LN3 (Cuadro 4), los estimadores de los parámetros fueron los siguientes:

Los valores de los eventos de diseño asociados con diferentes probabilidades de no excedencia, p, se estimaron al resolver la siguiente expresión:

Donde LN3^-1 es la función de distribución de probabilidad log normal inversa.

Por otro lado, para la estación 25082 (Cuadro 5), los parámetros estimados para el modelo M₁ de la función GA3 fueron los siguientes:

Para el escenario del año 2014:

Para el escenario del año 2044:

Las expresiones utilizadas para estimar los eventos de diseño fueron las siguientes:

Para el escenario del año 2014:

Para el escenario del año 2044:

Donde GA3^-1 es la función de distribución de probabilidad gamma inversa.

La representación gráfica de los eventos de diseño estimados para las estaciones 2060 y 25082, a partir de las expresiones anteriores, se muestran en la Figura 5.

Para la estación 2060, los eventos estimados G con una probabilidad de ocurrencia de 50% (pegriodo de retorno de dos años) fueron de 271.8 mm, con el modelo estacionario, y 266.8 mm con el modelo no estacionario. Con el enfoque tradicional, los eventos están sobreestimados, mientras que con el análisis no estacionario se predijeron valores de lluvia para las condiciones históricas más adversas de las sequías originadas por la persistencia de valores PDO negativos.

Al considerar la presencia de una tendencia decreciente de la lluvia en la estación 25082 fue posible simular el decremento en la magnitud de los eventos de diseño para escenarios futuros. Por ejemplo, para un periodo de retorno de dos años se predice una disminución de la lluvia de 77 mm en el año 2044 con respecto a la lluvia del escenario de 2014. Debido a que las pruebas estadísticas aplicadas a la serie de esta estación mostraron una tendencia decreciente significativa, la predicción del decremento en la lluvia para el escenario 2044 es consistente con el comportamiento de la serie.

De las demás estaciones cuyos modelos seleccionados (M₁ o M₃) están en función del tiempo, se estimó el incremento (decremento) de la lluvia que se tendrá para el año 2044 con respecto al año 2014, para un periodo de retorno de dos años. Este periodo de retorno corresponde al de la lluvia mediana (lluvia asociada con 50% de probabilidad de ocurrencia).

En la Figura 6 se muestra la variación espacial del incremento (decremento) de la lluvia, como porcentaje de la lluvia media anual; para las estaciones que no presentan cambios, los eventos de diseño se obtuvieron con modelos M₀ o modelos M₂, que no dependen del tiempo. En 29 estaciones de las regiones 1, 2, 3 y 5, ubicadas en la península de California, Sonora y Sinaloa, se predice una disminución menor a 20% de lluvia media anual (Cuadro 6). Cabe mencionar que en esta zona se tienen los valores de lluvia media anual más bajos del país, por lo cual estos resultados sugieren que será más susceptible a sequías en el futuro. En contraste, en ocho estaciones de la región 4 se esperan aumentos en la lluvia media de hasta 30%.

Conclusiones

En 83% de las series no estacionarias analizadas los modelos seleccionados para la estimación de eventos de diseño de lluvia acumulada en el noroeste de México incluyó una tendencia lineal y al índice PDO como covariable; por lo tanto, se determina que ante la falta de estacionareidad en las series de lluvia, los modelos no estacionarios propuestos explican mejor la variabilidad de los datos que los modelos convencionales. Para las series no estacionarias restantes, los modelos no estacionarios fueron mejores; sin embargo, la medida de discordancia indicó que para un nivel de significancia de 0.05, no se justifica el uso de un modelo no estacionario en lugar de un convencional. Por tal razón, el cambio en el comportamiento temporal de la lluvia no es una condición para todas las estaciones.

Los modelos no estacionarios permitieron la estimación de eventos de diseño asociados con las condiciones más adversas del índice PDO que han influido en los periodos históricos con déficit de lluvia más severos. Es posible que en las próximas décadas se presente de nuevo un periodo persistente de sequías, pues aparentemente a partir de 1998 se tiende hacia una fase negativa del PDO; por lo tanto, la persistencia de valores PDO negativos deben considerarse en los análisis de predicción de la lluvia.

Con la inclusión de una tendencia en el análisis de frecuencias se obtuvieron eventos de diseño asociados con diferentes periodos de retorno y escenarios futuros, con los cuales, en 29 estaciones ubicadas en la península de California, Sonora y Sinaloa, se predice una disminución para al año 2044 de menos de 20% de la lluvia media anual. En esta zona se tienen los valores de lluvia media anual más bajos del país, por lo cual tales resultados sugieren que será más susceptible a sequías en el futuro.