Nota: Este artículo de investigación es parte de Ingeniería–Revista Académica de la Facultad de Ingeniería , Universidad Autónoma de Yucatán, Vol. 22, No. 2, 2018, ISSN: 2448-8364.

INTRODUCCIÓN

La ciudad de Mérida está ubicada al norte de la península de Yucatán, en la región central de una planicie kárstica delimitada por una formación conocida como el anillo de cenotes (Perry et al. 1995, Pope et al. 1996). La altura promedio de la superficie del terreno en metros sobre el nivel medio del mar (msnmm) es de aproximadamente 9 m. La fisiografía que subyace a la ciudad, en la zona insaturada y en la saturada, está formada de horizontes kársticos en los que ocurren flujos de manera preferencial . Estos horizontes kársticos son cavidades de centímetros de espesor que se extienden de manera discontinua y pueden estar comunicados con cavernas (Figura 1). Tres de estos horizontes kársticos se han detectado en los primeros 15 metros de profundidad, debajo del del nivel del terreno. Dos de ellos se localizan en la zona insaturada y otro en la zona saturada entre los 12 y 13 m de profundidad . Subyaciendo a este material calcáreo , se encuentra una caliza semi consolidada formada por foraminíferos y una arenisca calcárea por donde el flujo subterráneo fluye de manera difusa, dando lugar a la característica dual de este sistema acuífero (Sánchez y Pinto et al. 2005).

Tradicionalmente las obras de ingeniería en la ciudad no requerían de excavaciones profundas; sin embargo, en los últimos años, se ha incrementado la construcción de este tipo de obras. Dada la proximidad del nivel freático a la superficie del terreno, de 4.8 a 7 m, en un transecto de norte a sur respectivamente , existe el riesgo que ocurra un ascenso extraordinario del nivel freático que alcance el nivel de piso terminado de estas obras, provocando anegamiento en ellos. La hidrología de los eventos extremos está enfocada al diseño de obras de control del agua para la protección de la vida humana, de la infraestructura en las zonas urbanas, etc. Para ello, se incorpora el análisis estadístico- probabilístico al comportamiento aleatorio de este tipo de eventos, con el fin de evaluar la probabilidad de que ocurra un evento extremo igual o mayor a aquel para el cual fue diseñada la obra. La presencia de un evento mayor al de diseño trae como consecuencia el fallo de la obra o del sistema. Las obras cuyo funcionamiento adecuado depende del análisis durante el cual existe la probabilidad de que se de la aleatoriedad de la variable, son diseñadas para que funcionen debidamente por un período de tiempo (conocido como período de retorno, TR, siglas en inglés) durante el cual existe la probabilidad de presente un evento que iguale o supere al de diseño (Cullen & Frey 1999; Stephenson A. & E. Gilleland 2006, Delignette-Muller &Dutang 2015).

El objetivo de este trabajo fue analizar las elevaciones máximas del nivel freático en el área urbana de la ciudad de Mérida, Yucatán, México, mediante modelos probabilísticos utilizados en el análisis de eventos extremos ; esto con el fin de proporcionar elementos que permitan a los diseñadores evaluar el riesgo a la inundación de obras construidas por debajo del nivel del terreno. Los datos de mediciones de la carga hidráulica (h) fueron realizadas en tres pozos someros. La longitud total de la serie de tiempo de la variación del nivel freático abarcó un período de tiempo entre 1982 y 2014. Debido a que los datos disponibles son limitados , las inferencias seapoyan también en métodos de remuestreo, en este caso: bootstrap paramétrico (Cooley 1997, Cullen & Frey 1999; Sáez 2009, Trichakis et al. 2011, Delignette-Muller & Dutang 2015).

MÉTODO

Los datos para el cálculo de la carga hidráulica se obtuvieron de mediciones utilizando sondas piezométricas (medidor de nivel de agua) y, en los últimos años, sensores de presión. La serie de tiempo consiste de un total de 28 datos de valores máximos anuales del nivel freático, los cuales fueron extraído de un registro semanal realizado en un periodo de 28 años.

El cálculo de la carga hidráulica se hizo tomando como base el nivel medio del mar reportado por el Instituto Nacional de Estadística e Informática (INEGI). Los pozos donde se realizaron las mediciones se localizan en la Ex-Facultad de Ingeniería de la Universidad Autónoma de Yucatán (EX FIUADY), el Observatorio meteorológico de la ciudad de Mérida (Observatorio), y el Campus de Ciencias Exactas e Ingeniería (FIUADY), (Figura 2).

El registro en la EX FIUADY (Antigua Facultad de Ingeniería de la Universidad Autónoma de Yucatán) se realizó de 1982 a 1992; en el OBSERVATORIO, desde el 2002 al 2014 y en el Campus de Ciencias Exactas e Ingeniería (FIUADY), de 1995 a 1999. Por esta razón, los datos fueron agrupados y se decidió una estrategia de análisis que incluye remuestreo paramétrico (bootstrap ). La estrategia general de estimación fue la siguiente:

1.- Preprocesamiento de datos de carga hidráulica y ajuste de las distribuciones de probabilidad utilizadas para el análisis de eventos extremos (Gumbel, Log-Gamma y Weibull) y estimación de parámetros mediante funciones de máxima verosimilitud y la función óptima (fitdistr del paquete MASS: R- project); la distribución Gumbel se ajustó por medio de iteraciones de mínimos cuadrados ponderados usando la función vgam del paquete VGAM de R-project (Yee 2010; Delignette - Muller & Dutang 2015).

Las funciones de distribución se eligieron como resultado de una revisión de literatura (Dupuis 1997, Stephenson & Gilleland 2006, Sáez 2009 , Trichakis et al.2011) cuyas fórmulas son las siguientes:

a) Distribución Gumbell

Tomada de Dupuis (1997) y sustituyendo ξ por α,y α por β de la ecuación original (ecuación 2).

b) Distribución Gamma

Donde xi= Máximo nivel del acuífero en el iésimo año

α = Parámetro de forma

λ = Parámetro de escala

Esta distribución es tratada por Sáez (2009) como sinónimo de la Distribución Pearson Tipo III. La función de verosimilitud para la distribución Gamma es:

xi= Máximo nivel del acuífero en el iésimo año.

Cuando es Log-Gamma:

c) Distribución Weibull

Ricci 2005.

2.-Pruebas de bondad de ajuste (Kolmogorov- Smirnof , Cramer -von Mises y Anderson Darling y los criterios de información de Akaike y el Bayesiano ) utilizando los parámetros previamente estimados. Las pruebas de bondad de ajuste se hicieron con la función gofstat e incluyeron las pruebas de Kolmogorov-Smirnov, Cramer-von Mises y el estadístico de Anderson -Darling así como el uso de los criterios de información de Aikaike y Bayesiano (D’Agostino & Stephens 1986, Ricci 2005, Yee 2010).

Se aplicaron las pruebas de bondad de ajuste de D’Agostino & Stephens (1986), sugeridas en Delignette & Dutang (2015 ), cuyas fórmulas generales se presentan en la Tabla 1.

3.- Estimación de crecidas para diferentes periodos de retorno de 25, 50 y 75 años utilizando bootstrap paramétrico y la mediana de los parámetros de la distribución teórica de probabilidad con 1000 repeticiones (muestreo aleatorio con reemplazamiento de las observaciones de carga hidráulica ; Efron & Tibshirani 1994 ).

RESULTADOS

Las pruebas de bondad de ajuste indican que los tres modelos tienen valores de significancia (P) que no permiten rechazar la hipótesis nula (Ho: Fn (x) = F(x)); sin embargo, el mejor ajuste lo presenta el modelo Log-Gamma con estadístico Kolmogorov- Smirnov = 0.085; estadístico Cramer-Von Mises = 0.029 y estadístico Anderson-Darling = 0.220; Los criterios de información de Aikake = 8.410 y el Bayesiano = 11.074. En la Tabla 2 se presentan los estadísticos de Kolmogorov -Smirnov para los tres modelos , el criterio de selección fue el valor de P estimado en la prueba de (KS).

Para períodos de retorno de 25, 50 y 75 años las crecidas máximas predichas por el modelo Log- Gamma (α=3.5547995 y β= 6.5128622) son de 3.13 m, 3.63 m y 3.94 m respectivamente.

Las estimaciones de crecidas del nivel freático con la distribución Log-Gamma y Tiempo de Retorno usando las medianas de alfa (3.765259) y beta (6.883555) obtenidas con Bootstrap son las siguientes 25 años = 3.07 msnmm; 50 años= 3.55 msnmm; 75 años= 3.85 msnmm y usando una Probabilidad 0. 99 las crecidas del nivel freático no excederían los 4.07 msnm. En la Figura 3 se representa la distribución de densidad, la línea de ajuste y el área de probabilidad de una crecida de la carga hidráulica sobre el nivel medio del mar con tiempo de retorno de 25 años.

DISCUSIÓN

El acuífero que subyace a la ciudad de Mérida es somero del tipo libre y está formado por una caliza del terciario. Su naturaleza kárstica acentúa la heterogeneidad de este, lo cual también es un factor que determina la porosidad dual que lo tipifica. La karsticidad también juega un papel importante en la elevada transmisividad del acuífero que lo hace ser considerado como altamente productivo (Sánchez y Pinto et al. 2005). Las crecientes del agua subterránea en el karst de la planicie norte de Yucatán se explican mediante un fenómeno similar a los movimientos de “marea” debidas a diferencias de carga hidráulica entre localidades situadas en un medio complejo que incluye flujos subterráneos en medios granulares y flujos turbulentos en canales cavernosos; esto se debe a que la composición del karst a diferentes escalas, desde metros a decenas, cientos de metros o kilómetros, tiene una estructura definida por procesos biogeoquímicos que ocurrieron en escalas de tiempo de cientos , miles y millones de años de tal modo que la complejidad de la estructura rocosa del karst es muy alta (Schmitter-Soto et al. 2002). En acuíferos extensos (escalas mayores a cientos de metros) las diferencias en la profundidad del nivel del agua en pozos, con respecto a un mismo plano de referencia, puede explicarse por varias causas: el gradiente hidráulico, zonas de mayor recarga o extracción y áreas con distinta permeabilidad. La teoría del movimiento del agua subterránea, basada en la ecuación de Darcy, establece que para que fluya un mismo caudal en acuíferos con la misma litología, pero diferente conductividad hidráulica debe incrementarse el gradiente hidráulico de aquel con menor conductividad hidráulica (Fetter , 2001 ). Para que esto suceda , es preciso aumentar la diferencia de la carga hidráulica entre un punto y otro. Es decir, cuando un mismo volumen de agua se recarga en dos zonas con diferente permeabilidad, la elevación del nivel del agua será mayor en aquella menos permeable. Lo anterior explica como la heterogeneidad del medio juega un papel importante en la profundidad del nivel del agua registrada en pozos ubicados en zonas con distinta permeabilidad. Por otra parte, los fenómenos hidrológicos como la precipitación, infiltración, evapotranspiración, movimientos horizontales y verticales de las aguas subterráneas del karst de la planicie norte de Yucatán no son homogéneos y se han trastocado por la intervención humana que incluye crecimiento urbano, inyecciones industriales y extracción minera, entre otros ( Hernández Terrones et al. 2011).

Toda esta complejidad estructural del karst y funcional de la cinética del agua subterránea se refleja en los cambios en la carga hidráulica en puntos discretos (pozos) y aunque las variaciones son continuas tanto en el tiempo como en él espacio, los métodos de observación de la carga hidráulica permiten solamente mediciones discretas tanto en el tiempo como en el espacio, de tal manera que existe la posibilidad de que los niveles máximos de la carga hidráulica no hayan sido observados debido a que las mediciones fueron realizados, posiblemente, antes o después del evento máximo.

Los análisis de series de valores extremos tradicional requieren al menos de un registro de 20 años a fin de reducir la incertidumbre en las predicciones que pueden hacerse a través de las funciones de distribución de probabilidad que mejor se ajusten. En este estudio no se contó con series de datos de 20 años para alguno de los pozos; en consecuencia, el análisis se hizo juntando las series de valores máximos registrados del nivel del agua de los pozos estudiados, y de esa manera alcanzar una serie de 28 observaciones de máximos anuales del nivel del agua, de manera adicional se utilizaron métodos numéricos de remuestreo (bootstrap) para hacer más robustas las estimaciones y dar mayor consistencia a los cálculos de los niveles del acuífero ante diferentes períodos de retorno. Los parámetros del modelo Log- Gamma estimados en este trabajo engloban el efecto que la heterogeneidad del medio tiene sobre el nivel del agua subterránea debido a la recarga pluvial del acuífero de la ciudad de Mérida (Cooley 1997, Trichakis et al. 2011).

No obstante estas limitaciones en las observaciones de la carga hidráulica máxima, el uso de modelos estocásticos permite aproximaciones que sirven como criterio de diseño de obras de ingeniería, y en este caso particular, las estimaciones de los parámetros de la distribución Log- Gamma (para calcular los máximos de carga hidráulica sobre el nivel medio del mar) se hicieron siguiendo una estrategia que incorpora toda la incertidumbre y variabilidad de las mediciones de tal manera que aunque los valores de los parámetros de la distribución Log- Gamma permiten el cálculo de las cargas hidráulicas sobre el nivel medio del mar en el área de la ciudad de Mérida se debe tener en cuenta que se trata de aproximaciones probabilísticas basadas en los datos disponibles hasta este momento.

CONCLUSIONES

El modelo más apropiado para calcular los valores máximos de la carga hidráulica, para diferentes tiempos de retorno en años en el área de la ciudad de Mérida es el Log-Gamma y los parámetros recomendables para estos cálculos son alfa (α = 3.765259) y beta (β = 6.883555). Las simulaciones realizadas para diferentes periodos de retorno ponen de manifiesto el riesgo que corre, por las posibles inundaciones, la infraestructura urbana subterránea en la región. Dada la heterogeneidad del medio estudiado y las características de los datos analizados se sugiere la aplicación de este modelo con las reservas del caso.

Referencias

Cooley RL. (1997). Confidence intervals for ground-water models using linearization, likelihood, and bootstrap methods. Ground Water 35: 869–880.

Cullen A.C. and Frey H.C. (1999). Probabilistic techniques in exposure assessment. Plenum Press, USA, pp. 181-241.

D’Agostino RB and Stephens MA (1986). Goodness-of-Fit Techniques. 1st edition. Dekker.

Delignette-Muller M.L. and Dutang C. (2015), fitdistrplus: An R Package for Fitting Distributions. Journal of Statistical Software, 64(4), 1-34.

Dupuis D.J. (1997). Extreme value theory based on the r largest annual events: a robust approach. Journal of Hydrology 200 (1997) 295-306

Efron B and Tibshirani RJ (1994). An Introduction to the Bootstrap. 1st edition. Chapman & Hall.

Fetter C.W. (2001). Applied hydrogeology. Fourth edition. Prentice Hall.

Hernández-Terrones L., Rebolledo-Vieyra M., Merino-Ibarra M., Soto M, Le-Cossec A. y E. Monroy-Ríos (2011). Groundwater Pollution in a Karstic Region (NE Yucatan): Baseline Nutrient Content and Flux to Coastal Ecosystems. Water Air Soil Pollut (2011) 218:517–528.

Perry, E.C., Marin, L.E., McClain, J., Velazquez, G., (1995). The ring of cenotes (sinkholes) northwest Yucatan, Mexico: its hydrogeologic characteristics and association with the Chicxulub impact crater. Geology 23, 17– 20.

Pope, K.O., Ocampo, A.C., Kinsland, G.L., Smith, R., (1996). Surface expression of the Chicxulub crater. Geology 24, 527–530.

Ricci V (2005). “Fitting Distributions with R.” Contributed documentation available on CRAN, URL http://CRAN.R-project.org/doc/contrib/Ricci-distributions-en. pdf.

Sáez A.J. 2009. Modelización estocástica de precipitaciones máximas para el cálculo de eventos extremos a partir de los periodos de retorno mediante R. http://www4.ujaen.es/~ajsaez/informe.html.LyXconv/informe.html

Sánchez y Pinto I., González-Herrera R. y E. Perry. (2005). Hydrodynamic behavior of the Yucatán aquifer. A perspective on the hydraulic conductivity estimation. Espeleo@digital 2 Ciudad de la Habana, Cuba. Pág. 8-20.

Schmitter-Soto J.J., Comín F.A., Escobar-Briones E., Herrera-Silveira J., Alcocer J., Suárez- Morales E., Elías-Gutiérrez M., Díaz-Arce V., Maín L.E. and B. Steinch. (2002). Hydrogeochemical and biological characteristics of cenotes in the Yucatan Peninsula (SE Mexico). Hydrobiologia 2002;467:215–28.

Stephenson A. & E. Gilleland (2006). Software for the analysis of extreme events: The current state and future directions. Extremes (2006) 8: 87–109. DOI 10.1007/s10687-006-7962-0

Trichakis, I., Nikolos, I. & Karatzas, G.P. 2011, "Comparison of bootstrap confidence intervals for an ANN model of a karstic aquifer response", Hydrological Processes, vol. 25, no. 18, pp. 2827-2836.

Yee T.W. (2010). The VGAM Package for Categorical Data Analysis. Journal of Statistical Software, 32(10), 1-34. URL http://www.jstatsoft.org/v32/i10/.