1. Introducción
Un lahar es una mezcla de detritos de diversos tamaños y agua lodosa, que se desplazan desde las laderas de un volcán y puede originarse de forma directa o indirecta por una erupción volcánica (Smith y Fritz, 1989). Los lahares se caracterizan por presentar gran movilidad, por lo que pueden recorrer grandes distancias, siendo altamente destructivos por los daños que pueden causar a los seres vivos y a la infraestructura, sea por el impacto directo, el transporte de sedimentos, el bloqueo de arroyos o el recubrimiento por transporte de escombros (Blong, 1984; Ortiz, 1996).
Se puede distinguir entre lahares primarios, aquellos relacionados directamente con un evento eruptivo, y secundarios, que no tienen relación con una erupción en particular (Vallance, 2000). Los lahares también pueden diferenciarse como: flujo de escombros (debris flow), flujo hiperconcentrado y flujo de arroyada de acuerdo con la proporción de agua y materiales sólidos que contengan, como se menciona más adelante (Smith y Lowe, 1991; Andrés et al., 2014).
Los volcanes activos y sus lahares causan repetidamente pérdida de vidas y daños. Por ejemplo, los eventos laháricos ocurridos en el Monte Santa Helena (E.E.U.U.) en 1980 (Janda et al., 1981; Pierson y Scott, 1985), los asociados a la erupción del Nevado del Ruiz (Colombia) en 1985 (Janda et al.,1981; Sigurdsson y Carey, 1986; Parra y Cepeda, 1990; Pierson et al., 1990) y los generados en el Volcán Pinatubo (Filipinas) en 1991 (Pierson et al., 1992; Rodolfo, 1989; Arboleda y Martínez, 1996) tuvieron efectos de alto impacto en la sociedad e infraestructura. El evento del Nevado del Ruiz (Colombia) fue el segundo desastre volcánico más importante del siglo XX. Para la prevención de este tipo de desastres asociados con volcanes y su actividad, muchas herramientas se han desarrollado para delimitar de forma aproximada las zonas de inundación por lahares alrededor de los volcanes y minimizar la vulnerabilidad de la población e infraestructura ante este fenómeno (Sheridan et al., 2001; Caballero y Capra, 2014; Córdoba et al., 2015). Algunas herramientas computacionales han sido ampliamente utilizadas (i.e. LAHARZ; Sheridan et al., 2001; Schilling, 2014) para mostrar de manera cualitativa y determinista el área de inundación de un lahar. Sin embargo, conforme la sociedad progresa, sus autoridades requieren de información que les permita tomar decisiones cada vez más complejas. Por ello, es fundamental utilizar herramientas que permitan, además de mostrar las áreas de inundación, evaluar probabilidades de parámetros que facilitan la toma de decisiones como son el área de inundación, espesor resultante de depósito y la presión dinámica de los flujos que revelan la vulnerabilidad y amenazan la vida e infraestructura.
Los mapas de peligro (también conocidos como de amenaza) probabilísticos tienen el objetivo de facilitar la cuantificación del riesgo, ya que combinados con la vulnerabilidad de los elementos expuestos ante una amenaza particular permiten este tipo de análisis. En este trabajo se realiza un análisis probabilístico del peligro por lahares en el flanco NE del volcán Popocatépetl, mediante el uso de TITAN2F, herramienta para elaborar modelos computacionales, diseñada para su aplicación en flujos laháricos (ver más en Córdoba et al., 2010; 2015). Aquí, se verifica el pronóstico de TITAN2F con los resultados de investigaciones previas relacionadas con la velocidad del lahar de 2001 realizado por Muñoz-Salinas et al. (2007b). Posteriormente se realiza un muestreo representativo de todos los casos posibles, es decir se generan escenarios hipotéticos que están dentro de las magnitudes posibles registradas en eventos históricos del volcán Popocatépetl, usando la técnica de Muestreo Hipercúbico Latino (LHS por sus siglas en inglés) (Núñez et al., 1999). Finalmente se llega al análisis probabilístico bayesiano (Coles y Tawn, 2005; Bulteau et al., 2015) de esta amenaza en la población de Santiago Xalitzintla.
2. Zona de estudio
El volcán Popocatépetl (5452 msnm) es uno de los estratovolcanes más activos en México, localizado en el sector central del Cinturón Volcánico Transmexicano (CVT) (Capra et al., 2004). En 1994 inició una nueva etapa eruptiva que continúa hasta el momento de escribir este trabajo (Delgado-Granados et al., 2001; Julio-Miranda et al., 2008; Caballero y Capra, 2014), el volcán es considerado como uno de los más peligrosos de México y de los de más alto riesgo en el mundo, puesto que se encuentran grandes poblaciones dentro de su área de influencia, como la Ciudad de México a 65 km con ~ 9 millones de habitantes y la ciudad de Puebla con 1.5 millones de habitantes en 2010, localizada a 45 km (INEGI, 2000; Andrés et al., 2014).
Hasta el año 2001, el volcán tenía un glaciar por debajo de la cumbre en la cara norte, con una superficie de 1.3 x 105 m2 (Delgado-Granados, 1997; Delgado-Granados et al., 2007). Hoy en día permanecen bloques de hielo (séracs) como remanentes de ese antiguo glaciar. Topográficamente abajo de los remanentes del glaciar hay tres gargantas a 4800 msnm, sus nombres de este a oeste son: Tenenepanco, La Espinera y Tepeteloncocone (Figura 1a), según la cartografía del INEGI (2015). Los cauces asociados a estas gargantas drenan la principal barranca de este sector conocida como Huiloac (3500 - 4000 msnm). La barranca Alseseca se localiza en el norte del volcán, inicia prácticamente desde el cráter y se desplaza en sentido N-NE (Capra et al., 2004). Las barrancas Huiloac y Alseseca unen su curso fluvial, tras cruzar el pueblo de Santiago Xalitzintla (Figura 1a) con una población de 2327 habitantes (Marcos et al., 2006; Andrés et al., 2007).
3. Ocurrencia de lahares
Los lahares son eventos frecuentes en el volcán Popocatépetl como lo indican los registros geológicos (Siebe et al., 1995; González-Huesca et al., 1997; Andrés et al., 2014; Castillo et al., 2015; Tabla 1) y los comentarios de los pobladores (conocidos por ellos como “avenidas”). Las magnitudes de estos eventos han variado desde pequeños lahares generados por alta intensidad de lluvias, hasta grandes eventos que por ejemplo, han recorrido 50 km de distancia, como los asociados a la actividad pliniana que dio origen a la pómez rosa (González-Huesca et al., 1997; Caballero y Capra, 2014). Éstos se movilizan a través de las gargantas de Tenenepanco, Huiloac y Alseseca y en ocasiones han sido un peligro para los habitantes de Santiago Xalitzintla (Palacios et al., 2001; Capra et al., 2004; Muñoz-Salinas, 2010; Tabla 1).
En diciembre de 1994 se inició el presente periodo eruptivo (Delgado-Granados et al., 2001), los principales eventos laháricos han ocurrido en 1995, 1997 y 2001 (Capra et al., 2004; Julio-Miranda et al., 2008; Andrés et al., 2014). En abril de 1995 flujos piroclásticos cubrieron aproximadamente el 50 % del glaciar, ocasionando la fusión del mismo canalizando el agua por las gargantas proglaciares (Tepeteloncocone, La Espinera y Tenenepanco) (Andrés et al., 2014).
Entre junio y julio de 1997 se generó una de las mayores erupciones de la reciente actividad (Capra et al., 2004; Julio-Miranda et al., 2008; Andrés et al., 2014). En la noche del 14 de junio se produjo un flujo de escombros que inundó el canal Alseseca (Sheridan et al., 2001). Las condiciones de alta pluviosidad facilitaron el transporte de sedimentos de las gargantas proglaciares generando también un flujo de escombros. Este lahar transportó 1.85 x 105 m3 de material sólido y agua que inundó Santiago Xalitzintla (Capra et al., 2004; Muñoz-Salinas et al., 2009a).
En la erupción del 22 de enero de 2001, se formó un flujo piroclástico que descendió 6 km a lo largo de una de las barrancas del flanco Norte desde el cráter y en un momento dado de su trayectoria, debido al emplazamiento del mismo flujo y la inclusión de agua se facilitó la formación del lahar (Julio-Miranda et al., 2005; Julio-Miranda et al., 2008; Caballero y Capra, 2014). No obstante esta interpretación, dado que el lahar ocurrió 4 horas después del evento explosivo, su origen podría estar asociado a la liberación de agua subglacial acumulada en forma de “bolsas” de acuerdo con observaciones de campo posteriores al evento. Estas bolsas de agua son comunes en el ambiente subglacial de los glaciares mexicanos tal y como se han documentado para el volcán Iztaccíhuatl (Alvarez y Delgado-Granados, 2002). El lahar resultante transportó el agua fusionada del glaciar y los materiales del flujo piroclástico por la garganta Huiloac y alcanzó los 7 m de ancho por 0.5 m de alto. Finalmente, se detuvo 2 km antes de llegar a la población de Santiago Xalitzintla (Capra et al., 2004; Julio-Miranda et al., 2005; Caballero y Capra, 2014). El volumen transportado se estima en 1.6 x 105 m3 y desarrolló velocidades entre 1.3 y 13.8 m/s (Muñoz-Salinas et al., 2007b).
El retroceso de los glaciares en el Popocatépetl ha alcanzado tasas de 40 m a-1 (Delgado-Granados, 1997; Huggel y Delgado-Granados, 2000; Julio-Miranda et al., 2008), lo cual disminuye la probabilidad de ocurrencia de lahares asociados a interacciones entre el glaciar y la actividad eruptiva (Delgado-Granados et al., 2015).
4. Análisis del peligro (Amenaza)
4.1. TITAN Dos Fases (TITAN2F)
TITAN2F es un modelo computacional basado en las leyes de balance de masa y momentum, tanto para la fase líquida, como para la fase sólida. Asume que el material granular obedece a la relación constitutiva de Coulomb, mientras que la fase fluida se asume que es no viscosa. A fin de tener en cuenta la fricción basal en la parte fluida, se usa la formulación de Darcy-Weisbach, junto con un coeficiente de arrastre semi-empírico que transfiera las fuerzas entre la fase sólida y líquida (Córdoba et al., 2010; 2015).
El programa requiere como datos de ingreso la localización de la pila o fuente del flujo, el volumen total, el porcentaje en volumen de la concentración inicial de sólidos, la velocidad inicial. El programa facilita el ahorro del consumo computacional porque éste se adapta a la malla 3D sobre la cual trabaja, además permite conocer campos de velocidad, concentración y presión dinámica del flujo en cada punto a lo largo de su trayectoria, siendo ésta una información muy útil para los tomadores de decisiones y para los procesos de gestión del riesgo. Ha sido puesto a prueba en distintos escenarios como el lahar del 2005 en el Valle de Vazcún, del volcán Tungurahua, Ecuador y los lahares de Villa La Angostura, Argentina (Córdoba et al., 2015).
4.2. Mapa probabilístico
La probabilidad es una herramienta que permite evaluar la confiabilidad de un estudio basado en información de muestras colectadas de una población (Mendenhall, 2008).
Para la evaluación del peligro, generalmente se propone que la probabilidad (H) dependa de la probabilidad de ocurrencia del evento (P) y la probabilidad de que un punto de interés, sea alcanzado por un nivel determinado de intensidad (I), expresándose como el producto de estas condiciones H = P × I (Maskrey, 1993).
En este estudio, se asume una probabilidad de ocurrencia P = 1, haciendo que el análisis probabilístico se enfoque únicamente en la determinación de los niveles de intensidad de los lahares (I), por lo cual se trata de un nivel de probabilidad condicional, es decir en caso de que ocurra el evento (Dalbey et al., 2008).
Debido a la complejidad del estudio de los peligros, existen muchas variables que condicionan su proceso, requiriendo como en este caso escenarios hipotéticos con diferentes combinaciones de las condiciones que puedan presentarse en la formación de lahares. De hecho, estas combinaciones, pueden ser potencialmente infinitas, pero para hacer algo práctico es necesario un muestreo que sea representativo de todos los casos dentro de los rangos de las condiciones de inicio, esto es lo que se conoce como el “Problema del Muestreo” (Goodman, 1960). En el caso de TITAN2F se deben estimar los rangos para las condiciones de: volumen, velocidad y concentración.
4.3. TITAN2F y DEMS
TITAN2F realiza la simulación del lahar sobre un modelo de elevación digital (DEM, por sus siglas en inglés) y la calidad del modelo tiene implicaciones en los resultados obtenidos por el programa. En consecuencia, como primer paso, es importante comprobar si TITAN2F reproduce las condiciones reales de la amenaza sobre el DEM.
Inicialmente se trabajó sobre un DEM del volcán Popocatépetl de 3 m de resolución, que aunque proporciona información a gran detalle, por otra parte, contiene diferentes elementos artificiales o ruido que generan incongruencias en el modelo del proceso lahárico y situaciones complejas que afectan los pronósticos de TITAN2F. La plataforma del Sistema de Información Geográfica (SIG) GRASS (Neteler et al., 2012) ofrece diferentes herramientas para tratar de corregir las imperfecciones de un DEM, sin embargo, la depuración del DEM de 3 m implica un trabajo más allá de los propósitos de esta investigación, por lo que se decidió usar un DEM alternativo, el DEM de 15 m de resolución, con menores elementos de ruido pero de menor resolución. Al ejecutar TITAN2F sobre este modelo, se comprobó que el flujo seguía trayectorias observadas en eventos pasados. En GRASS se procesó el modelo con el fin de disminuir imperfecciones bastante sutiles, utilizando la herramienta “r.denoise” que suaviza la superficie del DEM, eliminando el ruido y mejorando la claridad y calidad del modelo en cuanto a pendiente y comportamiento hidráulico (Neteler et al., 2012).
Para verificar los pronósticos que hace TITAN2F se compararon los resultados de velocidad del flujo modelado con los valores estimados por Muñoz-Salinas et al.(2007b) quien con el método de sobre-elevación (Chow, 1959) estima la velocidad a partir de las mediciones de campo realizadas para el lahar de 2001, como se verá más adelante.
El volumen del lahar de 2001 fue estimado por Julio-Miranda et al. (2005) como 4 x 105 m3 a partir de los depósitos de 2.4 x 105 m3, asumiendo una concentración de sólidos del 60 %. Sin embargo, con una concentración del 50 % (que resulta en un volumen de 4.7 x 105 m3) se obtienen los mejores resultados del modelo.
Los resultados de velocidad obtenidos para el modelado del lahar de 2001, fueron comparados con los reportados por Muñoz-Salinas et al.(2007b), para 11 sitios de depósitos evaluados, como se indica en la Tabla 2.
El análisis de confiabilidad se realiza aplicando la prueba Z no paramétrica de Kolmogorov-Smirnov (K-S), que consiste en determinar la diferencia mayor en valor absoluto entre la distribución acumulativa de una muestra Sn (xi) (observada en campo) y la distribución acumulativa teórica Fo(xi) que realiza el modelo numérico (Pérez-López, 2001; Arriaza-Balmón, 2006).
La diferencia es estadísticamente representativa si Sn (xi)-Fo(xi) supera el límite de errores aleatorios, este límite depende del número de muestras (N) y el nivel de significancia Ɛ. Este valor, conocido como K-S crítico, es tomado de las tablas de valores críticos sobre bondad de ajuste para el método K-S. En nuestro caso, para un número de muestras N = 11 y un nivel de significancia Ɛ = 0.01 corresponde a K-S = 0.47 (Kolmogorov, 1933; Smirnov, 1939). Se puede ver en la Tabla 2 que el valor de mayor diferencia es |-0.26|, al ser este valor menor al K-S = 0.47 tenemos que los resultados con TITAN2F son confiables.
4.4. Condiciones de inicio
4.4.1. Localización de las pilas
En la Figura 1b se observa que las barrancas Huiloac y Alseseca, cruzan y confluyen después de la población de Santiago Xalitzintla. Debido a que por cualquiera de ellas pudieran descender lahares, se decide localizar las pilas de material inicial en cada una de las gargantas.
4.4.2. Rangos de condiciones de inicio
En la Tabla 1 se presentan algunos de los registros geológicos de depósitos asociados a la ocurrencia de lahares en el volcán Popocatépetl, observándose que en él, se han generado flujos de lahares hasta del orden de 108 m3 (Sheridan et al., 2001) de acuerdo al glaciar existente entre 1919 - 1927. Sheridan et al. (2001) modelaron lahares entre 107 m3 y 108 m3. Por su parte, Huggel y Delgado-Granados (2000) consideran que volúmenes de este orden tienen una baja probabilidad de ocurrencia. El volumen de los lahares asociados a una erupción más grande (pliniana) es de 7 x 107 m3 (González-Huesca et al., 1997; Julio-Miranda y Delgado-Granados, 2003).
En investigaciones anteriores se han elaborado modelos con el software LAHARZ para diferentes escenarios de acuerdo al porcentaje de fusión del glaciar en la formación de lahares (Huggel et al., 2000) y también de acuerdo a las mediciones de los depósitos de lahares desde 1997 y 2001 los cuales han generado volúmenes que están en el orden de 105 m3 (González-Huesca et al., 1997; Huggel y Delgado-Granados 2000; Julio-Miranda y Delgado-Granados, 2003; Julio-Miranda et al., 2008; Muñoz-Salinas et al., 2009a, 2009b, 2010).
Huggel et al. (2000) indican que el volumen de sedimentos para el lahar de 1997 está entre 3.3 x 105 a 4 x 105 m3 y Julio-Miranda et al. (2005) establecen que el lahar de 2001 tuvo un volumen de 2.4 x 105 m3. Estos volúmenes se refieren al depósito, por lo cual asumiendo en el caso de mayor volumen de depósito (4 x 105 m3; Huggel et al., 2000) y una concentración del 50 % de sólidos, el rango de volúmenes de los flujos que pudieran presentarse estaría entre 2 x 105 m3 y 8 x 105 m3.
Con el fin de disminuir el costo computacional que representa modelar el flujo desde la cima con una velocidad inicial igual a cero (0), se decide ubicar la fuente del material en una parte intermedia, con un volumen distribuido en 3 pilas continuas (Figura 1b). Al no saber la velocidad que lleva el flujo en ese punto, se genera como incertidumbre adicional el valor de este parámetro, para determinar los rangos de velocidades mínima y máxima se hace necesario realizar simulaciones de prueba para los volúmenes mínimo y máximo. Resultando así que para los puntos de localización finalmente establecidos en este estudio, el rango de velocidades corresponde entre 1.8 - 10.2 m/s (Figura 2), las velocidades estimadas se encuentran dentro del rango del lahar de 2001 de 1.3 a 13.8 m/s (Muñoz-Salinas et al., 2007b).
El rango de concentración de sólidos en los lahares gobierna su naturaleza. En el caso de flujos hiperconcentrados, se tiene una concentración que varía de 35 % a 40 % en volumen (Fei, 1983; Major y Pierson, 1992) y en los flujos de escombros se tiene entre un 55 % y 60 % de sólidos (Pierson y Scott, 1985; Pierson, 1986). La concentración de sólidos en los flujos reportados por Julio-Miranda et al. (2005) y Muñoz-Salinas (2007a) en sus estudios de los lahares de 1997, 2001 y 2002 permiten considerar una concentración de sólidos entre 30 % y 50 % en volumen, cubriendo un espectro de flujos que van de hiperconcentrados a flujos de escombros (Pierson et al., 1985), los cuales reproducen los eventos mencionados. Concentraciones mayores no se consideran porque el modelado con concentraciones mayores no permiten alcanzar las condiciones observadas.
4.5. Muestreo probabilístico
Existen diferentes métodos de muestreo para generar los parámetros de entrada y éstos reflejan diferentes incertidumbres en los resultados del modelo computacional de los flujos de peligros naturales (Dalbey et al., 2008). Uno de ellos es el método de Monte Carlo (MC), que consiste en un muestreo aleatorio de los parámetros de entrada, lo cual hace fácil su implementación en cualquier problema sin importar cuán complejo sea. Sin embargo, MC necesita muchas combinaciones para tener validez estadística, lo cual implica que para la solución de problemas complejos se necesitan ordenadores de alta gama, por su costo computacional elevado (Dalbey et al., 2008) que haría intratable el problema. Por lo anterior, buscamos otro modelo que sea más eficiente computacionalmente, pero que tenga la misma precisión que MC. LHS es un muestreo estratificado desarrollado por McKay (1992). LHS implementa mejores propiedades de llenado de espacio mejorando la eficiencia en comparación con MC, sobretodo en la aplicación de modelos que requieren gran capacidad de almacenamiento, donde se necesita un número pequeño de muestras con un diseño óptimo (Núñez et al., 1999; Dalbey et al., 2008; Baalousha, 2009).
En la Figura 3 se indica una comparación de la distribución que hacen los dos métodos de muestreo anteriormente citados, para 16 muestras. Se identifica en la Figura 3a que MC hace una distribución no uniforme de las muestras, dejando ciertos vacíos en el espacio de muestreo, que son evidentes en los segmentos 0.2 - 0.4 y 0.6 - 0.8. También puede sobrevalorar ciertas zonas como en la última sección (0.8 - 1) en donde está cerca del 30 % del tamaño total de la muestra, razón por la cual serían necesarias más combinaciones para llenar estos espacios y tener una muestra que sea representativa.
El método LHS desarrolla una mejor distribución (Figura 3b) que MC, se observa que en cada división del espacio de muestreo se tienen igual número de combinaciones. Por lo anterior en este estudio usamos LHS y se realiza el muestreo mediante un script del software de libre distribución GNU OCTAVE que está diseñado para cálculos numéricos (Eaton, 1997).
4.6. Definición del tamaño de la muestra
El método geocronológico consiste en el reconocimiento y datación de los lahares a partir de sus depósitos, proporcionando conocimiento de algunos escenarios de eventos pasados (Muñoz-Salinas et al., 2010). Sin embargo, la cantidad de escenarios contemplados siguiendo este método (Tabla 2) no son una muestra representativa que permita llegar al cálculo de probabilidades, por lo cual se recurre a la opción de modelar escenarios hipotéticos a partir de distintas combinaciones de las condiciones de inicio, que estén dentro de los rangos factibles reconocidos a partir de la información geológica y observaciones sobre el volcán, pero asegurando su validez estadística.
En este trabajo se simularon 100 escenarios con las distintas combinaciones de condiciones iniciales para cada una de las fuentes de flujo, equivalentes a 10000 combinaciones de MC (Giner-Rubio, 1997). El error esperado se puede estimar con varias aproximaciones de diferente complejidad (Matala, 2008; Tang, 1993), siendo la más simple y más útil en nuestro caso, la expresada de la siguiente forma (Stein, 1987)
Dónde: εLHS = Error de LHS; N= Tamaño de la muestra.
Para este número de simulaciones y usando la Ecuación 1, el error en el muestreo LHS es de 10-2. Bajo el concepto de nivel de peligro aceptable, Connor (2011) estima que se deberían tomar medidas de mitigación para las estructuras vulnerables a partir de un nivel de probabilidad de ocurrencia de 10-1, bajo este parámetro el error asumido en este estudio es aceptable.
Por medio de un script en OCTAVE se generaron las distintas combinaciones de las condiciones de inicio (volumen, velocidad y concentración) para cada una de las pilas. La sistematización del ingreso de datos en cada una de las simulaciones se realizó con ayuda de rutinas adicionales de programación en OCTAVE, en las cuales el software se encarga de llamar a TITAN2F, proveer los datos de ingreso para cada simulación y a su vez organizar los resultados en diferentes directorios para facilitar su procesamiento posterior.
TITAN2F genera resultados agrupados en diferentes archivos de texto simple (.txt), en los cuales está contenida la información de los campos de velocidades, presión dinámica, concentración, así como altura de inundación en cada uno de los puntos de la malla computacional. La visualización de estos datos se hace utilizando los archivos de salida dentro de la plataforma GRASS GIS.
El procesamiento posterior de los resultados para el cálculo de probabilidad se hace mediante rutinas de OCTAVE, que tienen como base los principios del análisis de la probabilidad bayesiana descritos a continuación.
4.7. Análisis probabilístico del modelado
La determinación de la probabilidad depende del enfoque que se le dé al estudio, generalmente se hace una evaluación clásica que consiste en el cálculo de la frecuencia relativa que depende únicamente del número de ensayos y sus resultados repetidos. Sin embargo el hecho de tener dos fuentes del flujo con localizaciones distintas antes de confluir en la población, representa una partición del espacio de muestreo en dos subconjuntos independientes, generando las condiciones para que su análisis probabilístico apropiado sea el de tipo bayesiano (Walpole et al., 2007).
La probabilidad Bayesiana consiste en utilizar datos históricos como información previa para estimar una distribución a priori para los parámetros de distribución de probabilidad (Bulteau et al., 2015). Esta técnica proporciona un marco adecuado debido a su capacidad natural para el manejo de las incertidumbres, conduciendo a la representación probabilística de las estimaciones dadas por los modelos; además, permite la construcción modular en modelos complejos facilitando la inferencia en estos sistemas (Coles y Tawn, 2005; Bulteau et al., 2015).
El cálculo de la probabilidad se hace a partir de la inferencia bayesiana la cual se basa en la Regla de Bayes y permite actualizar o inferir la probabilidad de que una hipótesis pueda ser cierta según las evidencias u observaciones, que en este caso son los resultados de las simulaciones realizadas en TITAN2F (Bretthorst, 2013). Adaptando la Regla de Bayes (Mendenhall, 2008) a este estudio, el análisis probabilístico realizado se explica a continuación.
Se asume inicialmente que todas las fuentes tienen igual probabilidad de llegada a Santiago Xalitzintla puesto que el análisis se hace para cualquiera de las barrancas Huiloac o Alseseca y que en un primer momento no se cuenta con evidencias en datos simulados, expresamos su probabilidad como:
Dónde: P (Fi) = Probabilidad a priori (hipótesis) de que un flujo provenga de la fuente Fi; j = Número de fuentes de simulación.
Una vez que se realiza el modelado de las simulaciones requeridas para cada fuente, se calcula la probabilidad frecuencial:
Dónde: P(L|Fi) = Probabilidad frecuencial de la llegada del flujo al lugar L desde la fuente Fi; nf = Número de veces que llega el flujo al lugar L; ne = Número de ensayos de la fuente Fi.
Con lo anterior es posible tener una primera aproximación de la probabilidad total, que se calcula como la sumatoria de las probabilidades de las distintas fuentes desde i a j.
Dónde: P(Li)= Probabilidad total inicial de la llegada del flujo al lugar de estudio L.
Por inferencia bayesiana tenemos la probabilidad a posteriori que permite conocer la probabilidad de la hipótesis o primera distribución P (Fi) con la evidencia dada por P (L|Fi), así:
Dónde: P (Fi|L) = Probabilidad a posterioride llegada desde la fuente Fi dado el flujo en el lugar L.
El valor de P(Fi|L) indica verdaderamente la probabilidad de la fuente Fi, por lo cual es necesario realizar ese ajuste y volver a calcular la probabilidad total P(L) con la probabilidad de cada fuente P(Fi|L):
Dónde: P (L) es la probabilidad total ajustada por llegada de flujos laháricos para la población de Santiago Xalitzintla.
5. Resultados
5.1. Modelado del peligro
Las rutinas OCTAVE utilizadas para el análisis probabilístico generan un resumen de resultados en un archivo de tipo .txt que puede procesarse visualmente en la plataforma GRASS GIS facilitando la interpretación de esta información.
El script permite extraer los resultados de probabilidad para diferentes niveles o valores necesarios para las diferentes características del flujo (velocidad, presión dinámica, concentración, altura de inundación). Por ejemplo, para observar la probabilidad bayesiana de llegada de un flujo con una altura mayor a un valor particular o con un límite de presión dinámica que se considere importante y así, para cualquiera de las características que se requiera analizar.
Una vez procesados los resultados en el software GRASS GIS, se exportan estos datos para ser mostrados en un archivo de tipo Keyhole Markup Language (KML) que puede observarse sobre la plataforma Google Earth.
La Figura 4 presenta la distribución de probabilidad bayesiana para la llegada de un flujo con altura de inundación mayor a 0.2 m. Se observa que hay una probabilidad mayor del 80 % de que gran parte de la población de Santiago Xalitzintla sea inundada con esta profundidad.
Una herramienta muy útil en la gestión del riesgo es la información sobre la presión dinámica, que es la relación entre densidad del flujo y su velocidad al cuadrado. La presión dinámica o sobrepresión es la que produce daños en humanos y estructuras (Valentine, 1998; Spence et al., 2004; Zanchetta et al., 2004; Jenkins et al., 2015). En este trabajo, extraemos los datos de presión dinámica que arroja TITAN2F a fin de que puedan ser usados en conjunto con estudios de vulnerabilidad estructural como el realizado por Rodríguez et al. (2015).
En la Figura 5 se indican los resultados de la distribución de la probabilidad bayesiana para diferentes niveles de presión dinámica. Gráficamente se representa una gama de colores equivalentes a valores de probabilidad desde 0 a 1, en donde 0 corresponde a una probabilidad nula y 1 es el 100 % de probabilidad de llegada de un lahar con el valor de presión dinámica indicado en cada imagen.
En la Figura 5a se muestra la distribución de probabilidad para una presión dinámica mayor a 10 kPa, se observa que las barrancas que atraviesan tienen una probabilidad alrededor del 80 %. Sin embargo, esta área representa aproximadamente el 50 % del área de afectación por inundación mostrada en la Figura 4, lo cual indica que un área de inundación grande no significa que el flujo lleve necesariamente una presión dinámica considerable.
La Figura 5b muestra la distribución de probabilidades para un lahar con una presión dinámica mayor a 20 kPa, mostrando predominancia de una probabilidad con valores cercanos al 50 %. En las Figuras 5c y 5d, se muestran probabilidades para presiones dinámicas mayores a 30 kPa y 35 kPa respectivamente, para estas condiciones las probabilidades son menores estando entre 10 % a 20 %.
A medida que se incrementan los niveles de presión dinámica, el área afectada con altos valores de probabilidad disminuye, sin embargo la zona correspondiente a la barranca Huiloac presenta una alta probabilidad de llegada de lahares con presiones dinámicas superiores a los 35 kPa (Figura 5d). Con este nivel de presión dinámica las estructuras se verían altamente afectadas (Valentine, 1998; Spence et al., 2004; Jenkins et al., 2015) haciendo que en esta zona se deba tener especial consideración ante este peligro.
5.2. Consideraciones peligro-riesgo
Los resultados de probabilidad bayesiana para niveles de presión dinámica bajos como 10 kPa y 20 kPa, indican que las zonas comunes localizadas en el centro de la población se encuentran altamente amenazadas (Figura 5). Pese a que son niveles bajos, con esta presión dinámica hay mayor susceptibilidad de arrastre de seres humanos y una leve afectación en las edificaciones.
Para una presión dinámica que sobrepase los 35 kPa, gran parte de la población presenta una probabilidad por debajo del 25 %, sin embargo encontramos una zona crítica localizada en la desembocadura de la barranca Huiloac, en ella, la probabilidad de sobrepasar una presión dinámica de 35 kPa es cercana al 80 %. En esta zona se encuentra un campo deportivo y algunas viviendas, teniendo en cuenta que este nivel de presión dinámica genera daños severos en edificaciones (Valentine, 1998; Spence et al., 2004; Jenkins et al., 2015), este estudio otorga información de utilidad para que los tomadores de decisiones puedan considerar medidas al respecto.
6. Conclusiones
En este trabajo desarrollamos una metodología para el análisis probabilístico del peligro por lahares en la población de Santiago Xalitzintla. Mediante la utilización de las condiciones iniciales (físicas y geológicas) presentadas en diferentes estudios acerca de eventos anteriores, se logró hacer la verificación del software TITAN2F y determinar los rangos de condiciones iniciales. El software TITAN2F generó resultados aceptables que permitieron hacer un uso confiable del programa en un ambiente más amplio como lo es el análisis probabilístico del peligro.
Gracias al muestreo de tipo estratificado con propiedades adecuadas de llenado se generaron escenarios hipotéticos que abarcaron todos los casos posibles con un error en el muestreo del orden de 10-2 a partir del cual se deberían tomar medidas de mitigación para las estructuras. Cabe mencionar que una muestra diseñada para un nivel en el que se consideren las comunidades, implicaba un costo computacional mayor debido al incremento en la muestra requerida. Sin embargo, para las características mencionadas los escenarios modelados (100 por cada Barranca) fueron representativos de las condiciones de inicio y estaban dentro de los rangos físicos y condiciones geológicas posibles del volcán Popocatépetl.
El uso de herramientas computacionales permite realizar cálculos de gran complejidad y a su vez conocer características físicas de los lahares con mayor nivel de precisión, puesto que además de mostrar la zona de afectación, se puede conocer la intensidad del fenómeno en cada punto de ella. En este estudio se presenta una primera aproximación, teniendo en cuenta que el DEM usado no tenía implementada la configuración urbanística de Santiago Xalitzintla, sin embargo la información resultante es un primer insumo de gran utilidad para el posterior análisis de la vulnerabilidad y gestión del riesgo de estos fenómenos.
TITAN2F respondió adecuadamente sobre el Modelo de Elevación Digital (DEM) utilizado, representando situaciones complejas como la que ocurre en la Barranca Huiloac cuando comienza a cruzar la población de Santiago Xalitzintla, en donde se observa que debido a la intervención humana, la topografía ha cambiado de tal forma que facilita el direccionamiento de los flujos hacia zonas comunes de la población.
El presente trabajo se constituye como una herramienta en el conocimiento de la amenaza por lahares en la población de Santiago Xalitzintla, y proyecta una metodología que puede ser aplicada en otros escenarios y poblaciones que se vean afectados por amenazas geológicas de este tipo.