Introducción

Los patrones de las variables hidrometeorológicas están cambiando (^{Gleick, 1989}; ^{Voss, May, & Roeckner, 2002}; ^{Held & Söden, 2006}) y se espera que en el futuro los desastres provocados por eventos hidrometeorológicos sean más frecuentes e intensos debido al calentamiento global (^{IPCC, 2007}). Una forma de reducir el riesgo de tales desastres es por medio de una planeación basada en la probabilidad de los eventos extremos. Por ejemplo, para reducir el riesgo de inundaciones se deben construir estructuras hidráulicas para drenaje y control de escurrimientos, diseñadas a partir de eventos de lluvia asociados con un periodo de retorno determinado.

La metodología clásica para predecir eventos de diseño de variables hidrometeorológicas se basa en la teoría de valores extremos de series de tiempo estacionarias (^{Khaliq, Ouarda, Ondo, Gachon, & Bobée, 2006}; ^{Villarini & Smith, 2010}). Bajo condiciones de estacionaridad, la distribución de la variable de interés es invariante en el tiempo, sin tendencias, cambios ni periodicidades (^{Villarini & Smith, 2010}). No obstante, la estacionaridad de las series de tiempo ha dejado de existir debido al cambio en el clima del planeta, que está alterando la media y los extremos de la precipitación, evapotranspiración y gastos de los ríos (^{Milly et al., 2008}). La no estacionaridad de las series hidrometeorológicas debida al cambio climático por lo usual se manifiesta como una tendencia o un salto repentino en las características estadísticas de los datos (^{Khaliq et al., 2006}). La presencia de una tendencia tiene un efecto considerable en el ajuste de una distribución de probabilidad convencional a una muestra de observaciones no estacionarias; por lo tanto, la comprensión y el diagnóstico de un comportamiento no estacionario es muy importante (^{Khaliq et al., 2006}).

En México existe evidencia del incremento en la frecuencia e intensidad de los eventos extremos (^{Jáuregui, 2000}; ^{Peralta-Hernández, Balling, & Barba-Martínez, 2009}), por lo que las técnicas tradicionales de predicción son incapaces de estimar con precisión su magnitud. Si la predicción de eventos extremos de diseño no es adecuada, las obras hidráulicas tienen mayor riesgo de falla, y sus efectos pudieran ser socialmente inaceptables aguas abajo. Por el contrario, una sobreestimación de los eventos incrementaría la seguridad de la estructura, pero con elevados costos económicos (^{Escalante-Sandoval & Reyes-Chávez, 2010}). Ante esta situación, se requieren emplear técnicas que consideren la falta de estacionaridad de las series. Sin embargo, para México, son pocas las investigaciones (^{Campos-Aranda, 2016}; ^{Álvarez-Olguín & Escalante-Sandoval, 2016}; ^{Álvarez-Olguín & Escalante-Sandoval, 2017}) que incluyen la falta de estacionaridad de las series en los análisis de predicciones hidrometeorológicas. Por lo anterior, el objetivo de este trabajo fue estimar para la república mexicana eventos de lluvia máxima en 24 horas, asociados con periodos de retorno de 2, 5, 10, 50 y 100 años, mediante modelos probabilísticos no estacionarios. Los resultados muestran la influencia de la variabilidad climática en la predicción de la lluvia máxima, lo cual es útil para obtener mejores estimaciones que podrían contribuir a disminuir el riesgo de falla de las obras hidráulicas.

Datos y métodos

Los datos utilizados corresponden a la lluvia diaria medida en las estaciones climatológicas convencionales del Servicio Meteorológico Nacional de la Comisión Nacional del Agua (Conagua). Esta información se extrajo de la base de datos Clicom (Climate Computing Project) y sólo se consideraron estaciones con al menos 64 años de registros en el periodo de 1950 a 2013, y por lo menos 90% de información disponible. La localización de estas estaciones se muestra en la Figura 1. Para deducir datos faltantes se aplicó el método de interpolación de la distancia inversa ponderada (^{Shepard, 1968}), considerando dos estaciones de apoyo y un exponente de distancia igual a dos. Se verificó si existe mejora en la estimación de la media y la varianza de las muestras extendidas, mediante el criterio de información relativa (I), definido con la Ecuación (1), descrita por ^{Escalante-Sandoval y Reyes-Chávez (2002)}.

Donde VarSy12 es la varianza de la varianza de la serie original y VarSy2 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; en consecuencia, la extensión de los registros es adecuada. De esta forma, se obtuvieron datos confiables de lluvia diaria de 769 estaciones, con los que se lograron las series anuales de lluvia máxima en 24 horas, para el periodo comprendido entre 1950 y 2013.

El siguiente paso fue verificar la estacionaridad de las series anuales de lluvia máxima en 24 horas. Para comparar de modo espacial los cambios ocurridos en las series de lluvia, se tomó como base la regionalización de la lluvia obtenida por ^{Álvarez-Olguín y Escalante-Sandoval (2017)}. La homogeneidad de las series se verificó mediante la ausencia de cambios abruptos en la media, a través de la prueba de ^{Pettitt (Pettitt, 1979}). La presencia de tendencias en las series se evaluó con la prueba de Mann- Kendall (^{Kendall, 1975}), aplicada a las series anuales de lluvia máxima en 24 horas. Las tendencias también se evaluaron mediante la gráfica de los residuales de las series obtenidas por el método de descomposición de modos empíricos por conjuntos (EEMD, Ensemble Empirical Mode Decomposition), desarrollado recientemente (^{Huang et al., 1998}; ^{Wu & Huang, 2009}) para descomponer señales en funciones de modos intrínsecos (Intrinsic Mode Function). En este método, el cociente de la desviación estándar de la señal de ruido blanco gaussiano adicionada a los datos fue de 0.2.

Es común que la no estacionaridad de una serie hidrometeorológica se deba a la presencia de una tendencia y/o un salto repentino en las características estadísticas de los datos (^{Yevjevich & Jeng, 1969}; ^{Khaliq et al., 2006}). La presencia de una tendencia tiene un efecto considerable en la interpretación de los resultados cuando se ajusta una distribución de probabilidad convencional a una muestra de observaciones no estacionarias (^{Khaliq et al., 2006}). No se ha desarrollado una teoría general para el análisis de frecuencia de procesos no estacionarios, pero una práctica aceptada consiste en modificar las distribuciones de probabilidad estándar para obtener versiones no estacionarias de las mismas (^{Coles, 2001}). Una distribución no estacionaria combina parámetros de la distribución estacionaria con algunas covariables, como el tiempo o índices de circulación atmosférica a gran escala (ENOS, DOP, OMA y NAO), para generar un nuevo tipo de parámetro para la estimación. Por ejemplo, de la distribución general de valores extremos (GVE), con parámetros (, ( y β (ubicación, escala y forma, respectivamente), se puede obtener un modelo para estimar una variable x _t en función del tiempo, como se indica en la Ecuación (2).

Donde cada parámetro t, αt y βt tiene una expresión en términos del tiempo. El parámetro t se puede expresar como se muestra en la Ecuación (3).

Donde t es un índice que depende del tiempo. Así, las variaciones a través del tiempo en los procesos observados se modelan como una tendencia lineal.

En esta investigación, el análisis de las series se realizó con las funciones para valores máximos descritos en la Tabla 1 y los modelos incluidos en la Tabla 2. Para el caso estacionario, se utilizaron los modelos M ₀; los modelos M ₁, M ₂ y M ₃ corresponden al caso no estacionario, y sus parámetros se expresan en términos del tiempo e índice de oscilación decadal del Pacífico (PDO). Este índice se seleccionó por ser una de las fluctuaciones de periodo largo de mayor correlación con la lluvia en México (^{Álvarez-Olguín & Escalante-Sandoval, 2017}). Las series anuales del índice PDO se obtuvieron del promedio de los valores mensuales para el periodo comprendido entre 1950 y 2013 (^{NOAA, 2014}).

La estimación de los parámetros de los modelos se realizó por medio del método de máxima verosimilitud, que se considera el más eficiente, pues proporciona la menor varianza muestral de los parámetros y de los eventos estimados, en comparación con otros métodos. Además, tiene la ventaja de adaptarse a cambios en la estructura del modelo (^{Katz, Parlange, & Naveau, 2002}).

Para el caso de la función de densidad de probabilidad GVE (general de valores extremos), la función logarítmica de verosimilitud se expresa como:

Donde t,α y β son los parámetros de ubicación, escala y forma, respectivamente.

Para el modelo GVE_M ₃ (modelo M ₃ de la función general de valores extremos), se consideró al parámetro de escala en función del tiempo y del índice PDO, por lo tanto, la función logarítmica de verosimilitud propuesta es:

Donde υ0,υ1 y υ2 son los parámetros de ubicación; α es el parámetro de escala y β es el parámetro de forma; t es un índice que depende del tiempo (t = 1 para el año 1950), y PDO _t es el valor promedio anual del índice PDO correspondiente a t.

Los modelos que mejor describen la variabilidad de los datos se seleccionaron mediante el criterio de información de Akaike (^{Akaike, 1974}), donde se utilizó la Ecuación (6). El mejor modelo es el que tuvo el menor valor de AIC. En el caso de tener valores similares, 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.

Además, 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 Ecuación (7). Valores grandes de D indican que el modelo M ₁ es más adecuado y explica mejor la variación de los datos, en comparación con el modelo M ₀. El estadístico D se distribuye de acuerdo con la distribución chi-cuadrada ((² ₍). 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 (² ₍ para un nivel de confianza particular se consideran significantes, entonces se rechaza M ₀ en favor del modelo M ₁:

Donde liMi es la función logarítmica de verosimilitud maximizada del modelo M _i .

Para verificar que los modelos seleccionados se ajustan de forma adecuada a los datos observados, fue necesario transformarlos mediante una estandarización (^{Krzysztofowicz, 1997}). De esta manera, tras ordenar los valores estandarizados y asociarlos con los correspondientes de la distribución empírica de Weibull, se elaboraron las gráficas de residuales (Q-Q) y worm plots (^{Buuren & Fredriks, 2001}). En una “worm plot”, el eje vertical es la diferencia entre los valores empíricos y los teóricos, y contiene el intervalo de confianza al 95%, estimado como:

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

Por último, se estimaron eventos de diseño relacionados con diferentes periodos de retorno. Para el modelo GVE_M ₃, los eventos de diseño se estimaron por medio de la Ecuación (9):

Donde υ^0, υ^1 y υ^2 son los estimadores del parámetro de ubicación α^ es el estimador del parámetro de escala; β^ es el estimador del parámetro de forma; t ₀ es el valor del índice del tiempo para un año determinado; PDO^ es un valor promedio del índice PDO yX^ es el evento de diseño asociado con una probabilidad de no excedencia p. El valor PDO^ se estimó como la media aritmética de los valores PDO correspondientes al mayor periodo de anomalías, es decir, al mayor periodo con persistencia de eventos superiores a cero.

Los datos utilizados en los análisis son del periodo 1950-2013. Por lo tanto, al primer año (1950) le corresponde el valor de t ₀ = 1; para los años siguientes, el valor de t ₀ aumenta de modo consecutivo una unidad cada año. De esta forma, al escenario del año 2013 se asignó un valor de t ₀ = 64 y al escenario de 2040 un valor de t ₀ = 94.

Resultados y discusión

A través de la prueba de Pettitt se verificó la presencia de cambios abruptos en las series, mientras que con las pruebas de Mann-Kendall y EEMD se evaluaron tendencias. Las series identificadas como no estacionarias a un nivel de significancia de 0.05, se clasificaron en los siguientes tipos:

En total, 22% (167) de las series analizadas se consideró no estacionaria (Figura 2). Corresponde a estaciones ubicadas en las regiones 1 y 2, en el centro y noroeste del país, donde se tienen 67 y 32 estaciones, respectivamente (Tabla 3). De estas series, 89% tiene puntos de cambio ascendentes y/o tendencias monótonas crecientes significativas.

En la Figura 3 se muestra la variación temporal de la lluvia máxima y los residuales de algunas series no estacionarias; para éstas se determinaron tendencias monótonas crecientes significativas. También se identificaron cambios ascendentes significativos en las series de las estaciones Acámbaro, El Novillo, Presa Rodríguez, La Servilleta y Cerro Ortega; en éstas, tanto la prueba de Pettitt como de Mann-Kendall indicaron incrementos abruptos y graduales, respectivamente. En el caso de Acámbaro, se encontró un punto de cambio en 1987 y una tendencia creciente (Tabla 4). Además, con la gráfica del residual obtenido por EEMD se aprecia que la tendencia es monótona.

Las tendencias crecientes indican que en los próximos años las lluvias extremas serán más intensas y habrá mayor probabilidad de que ocurran desastres relacionados con estos eventos, sobre todo en las regiones del centro y noroeste del país. El aumento en las intensidades de lluvia en Norteamérica está referido por ^{Peterson, Zhang, Brunet‐India y Vázquez‐Aguirre (2008)}, y ^{Cavazos et al. (2008)}. Entre las causas probables de dicho cambio están el aumento del contenido de humedad atmosférica derivado del calentamiento global (^{IPCC, 2007}; ^{Meehl, Arblaster, & Tebaldi, 2005}) y la influencia de fenómenos oscilatorios que regulan la variabilidad climática global.

Hay evidencias de la relación entre el comportamiento de la lluvia en México y fenómenos oscilatorios de gran escala como ENSO (El Niño Southern Oscillation), PDO (Pacific Decadal Oscillation) (^{Magaña, Vázquez, Pérez, & Pérez, 2003}; ^{Méndez-González, Návar-Cháidez, González-Rodríguez, & Treviño-Garza, 2007}; ^{Méndez-González, Ramírez-Leyva, Cornejo-Oviedo, Zárate-Lupercio, & Cavazos-Pérez, 2010}; ^{Méndez & Magaña, 2009}). Por lo tanto, es de esperarse que influyan en las características estadísticas de las series. Ejemplo de la influencia de la PDO sobre las lluvias máximas en el noroeste del país se muestra en la serie de la estación San Vicente, Ensenada, Baja California (Figura 4), que se localiza en la región 5, la única del país con clima mediterráneo, donde 46.7% de la aportación a la precipitación anual se presenta durante el invierno (^{Álvarez-Olguín & Escalante-Sandoval, 2017}). Para la serie, se identificó un punto de cambio ascendente significativo en el año 1975 (Tabla 4), el cual coincide con el inicio de la fase positiva de la PDO, que abarcó el periodo 1976-1998. Se observó que durante esta fase se presentaron tres de los eventos más intensos, en los años 1978, 1993 y 1997. La fase positiva de PDO favorece las lluvias invernales en el norte de México (^{Mantua, Hare, Zhang, Wallace, & Francis, 1997}; ^{Méndez et al., 2010}), por lo cual la falta de estacionaridad de la serie se puede atribuir a la influencia de la PDO.

Como resultado del análisis de frecuencias, se determinó que para 83% (139) de las series no estacionarias se seleccionó modelos no estacionarios, como los que mejor representan la variabilidad de los datos. Para las 28 series no estacionarias restantes, si bien con el estadístico AIC los modelos no estacionarios aparentemente fueron mejores, el estadístico D fue inferior a 3.84 (valor de (² ₍ a un nivel ( = 0.05), por lo que no se justificó estadísticamente la utilización de un modelo M ₁, M ₂ o M ₃ en lugar de un M ₀. Por otro lado, predominan los modelos de tipo GVE_M ₁ elegidos para 91 series (Figura 5). La mayoría de éstas pertenece a estaciones localizadas en las regiones 1 y 7 (centro del país), con un total de 42 y 10, respectivamente (Tabla 5).

En la región 5 de la península de California se tiene la mayor cantidad de estaciones (23), con modelos LN3_M ₂ y LN3_M ₃, en los que el parámetro de escala está en función del índice PDO. La selección de modelos del tipo M ₂ y M ₃ ratifica la asociación entre el índice PDO y las lluvias máximas en la península.

Al comparar la Figura 2 y la Figura 5 se aprecia una correspondencia entre los sitios del país con series no estacionarias y los sitios para los que, a su vez, se han seleccionado modelos no estacionarios, por lo cual se puede determinar que los resultados bajo el enfoque no estacionario son consistentes con las pruebas de estacionaridad realizadas.

En la Tabla 6 se muestran los resultados de las pruebas de bondad de ajuste y los estimadores de los parámetros de las estaciones Acámbaro y San Vicente, elegidas para mostrar la metodología seguida. Para San Vicente, los estadísticos de bondad de ajuste indicaron que el valor más bajo del AIC es de 517.53 y corresponde al modelo LN3_M ₂. Además, el estadístico D del modelo M ₂ con respecto al M ₀ es de 12.68, superior al valor (² = 3.84. Por lo tanto, el modelo M ₂ es más adecuado y explica mejor la variación de los datos que el modelo M ₀. En el caso de la estación Acámbaro, tanto el modelo G_M ₁, como el modelo GVE_M ₁ tienen valores AIC bajos, iguales a 496.8; sin embargo, el estadístico D para este último es de 14.0 (superior al de la función Gumbel), por lo cual se consideró que el mejor modelo es el GVE_M ₁.

Los modelos seleccionados para las series de las estaciones San Vicente y Acámbaro (Figura 6) indican que, para las dos estaciones, los datos de las worm plot están dentro de los límites de confianza, y en las gráficas cuantil-cuantil (Q-Q) los datos están cerca de la diagonal unitaria. Por tal razón, los modelos seleccionados se ajustan de forma adecuada con los datos.

Los eventos de diseño de la estación San Vicente asociados con diferentes periodos de retorno se estimaron al resolver la Ecuación (10), Ecuación (11) y Ecuación (12). El valor PDO^ de la Ecuación (11) se estimó como la media aritmética de los valores PDO mensuales de la fase positiva de esta oscilación, comprendida en el periodo 1976-1998. A pesar de que el modelo M ₂ seleccionado no depende del tiempo, se consideró que es un modelo no estacionario. De acuerdo con ^{Villarini y Smith (2010)}, bajo condiciones de estacionaridad, la distribución de la variable de interés es invariante en el tiempo, sin tendencias, cambios ni periodicidades, por lo que se considera válido asumir que un modelo tipo M ₂ es no estacionario debido a la periodicidad intrínseca del PDO incluido como covariable:

Donde F ^-1 es la función de distribución de probabilidad log normal inversa, y p es la probabilidad de no excedencia en función de un periodo de retorno T, p = 1-1/T.

Los eventos de diseño obtenidos con el modelo convencional son inferiores a los estimados con el modelo no estacionario (Figura 7). La diferencia es mayor en periodos de retorno de 2, 5 y 10 años, en comparación con los periodos de 50 y 100 años. Por ejemplo, los eventos estimados para un periodo de retorno de 10 años fueron de 60.02 y 67.03 mm, para los casos estacionario y no estacionario, respectivamente. Con lo anterior, se puede determinar que al no considerar la influencia del PDO en la variabilidad de la lluvia en la estación se tendría una subestimación de los eventos de lluvia máxima, principalmente en los de menor periodo de retorno.

Por otro lado, los eventos de la estación Acámbaro se obtuvieron con la Ecuación (13) y la Ecuación (14). Al considerar una dependencia lineal con respecto al tiempo en el parámetro de ubicación de la función GVE, se pudo predecir el incremento de la lluvia máxima de diseño en el futuro. Así, la lluvia máxima de diseño para un periodo de retorno de 10 años será de 73.1 mm en el año 2020; para el año 2050, para el mismo periodo de retorno, se predice una lluvia máxima de 79.6 mm (Figura 8), es decir, la lluvia se incrementará 6.5 unidades en un lapso de 30 años:

Donde t ₀ es el valor del índice del tiempo para un año determinado, t ₀ = 71 para el año 2020 y t ₀ = 101 para el año 2050.

Conclusiones

Se encontró evidencia de que 22% de las series de lluvia máxima analizadas es no estacionaria. Para 83% de éstas, los modelos no estacionarios seleccionados explican mejor la variabilidad de los datos. Con respecto a las series no estacionarias restantes, el criterio de información de Akaike indicó que los modelos no estacionarios fueron mejores, sin embargo, el estadístico de discordancia no justificó estadísticamente la utilización de un modelo M ₁, M ₂ o M ₃ en lugar de un M ₀.

Los modelos que incluyen al PDO como covariable son mejores que los estacionarios para el noroeste del país. Por otro lado, la incorporación del tiempo como un índice en los parámetros de las funciones genera mejores ajustes en las series que presentan puntos de cambio y tendencias, de las cuales, en 89% fuer ascendente y creciente, respectivamente. Para estas series, los modelos no estacionarios del tipo M ₁ se ajustan mejor a los datos. Para la estación Acámbaro, los eventos de diseño asociados con diferentes periodos de retorno indican que con el modelo no estacionario GVE_M ₁ se predice un incremento de los valores de la lluvia máxima en el futuro.

El análisis de las series de las estaciones San Vicente y Acámbaro evidenció que, bajo el enfoque clásico estacionario, se subestima la lluvia de diseño y, en consecuencia, su uso en el diseño de alguna obra hidráulica implicaría mayor riesgo de falla.