1. INTRODUCCIÓN
El Cotopaxi es un volcán activo que tiene erupciones devastadoras, con periodo recurrente de alrededor de un siglo. La última erupción que ocurrió en el año de 1877, destruyó parcialmente Latacunga y llegó hasta el valle de los Chillos. Actualmente, las zonas amenazadas son altamente pobladas y de gran interés económico para el Ecuador, por esta razón, es necesario determinar con mayor precisión las zonas susceptibles al flujo de lahares (flujos de lodo, rocas, y escombros mezclados con agua provocados por la fusión casi instantánea de volúmenes de hielo importantes) (Cáceres, y otros, 2004).
En el año 2004, se realizó la modelación numérica unidimensional en el flanco norte del volcán, utilizando un modelo digital de terreno (MDT) con resolución espacial de 30 por 30 metros de tamaño de pixel (EMAAP & EPN, 2004). En la presente investigación, se utiliza un MDT con resolución de 3 por 3 metros de tamaño de celda, lo cual, permitirá representar con mayor especificidad las zonas analizadas, aumentando la fiabilidad de resultados. De igual manera, se utiliza información del área del glaciar, procedente de estudios recientes (Cáceres B. , 2010) y se aplican herramientas computacionales de última generación, tales como: HEC-RAS 5.0, IBER 2.3 y FLOW 3D 11.1.
La modelación contempla una longitud de cauce de alrededor de 65 kilómetros, donde posiblemente el flujo de escombros, producto de una erupción, cruce las conducciones de tres sistemas de suministro de agua potable para la ciudad de Quito (Pita – Tambo, La Mica y Papallacta), lo cual pone en riesgo el abastecimiento del recurso hídrico para el 80% de la población. Además, existen varias zonas con alta densidad poblacional, importantes vías de comunicación, numerosas industrias y fábricas.
2. OBJETIVOS
Actualizar la modelación numérica unidimensional del flujo generado por un evento eruptivo en el flanco norte del volcán Cotopaxi, hasta la zona densamente poblada de Sangolquí.
Calibrar el modelo numérico unidimensional con base en los vestigios de la erupción del 26 de junio de 1877.
Determinar los niveles y zonas de inundación producidas por el flujo de escombros.
Analizar bidimensionalmente el comportamiento del flujo en un tramo de la zona densamente poblada de Sangolquí.
3. METODOLOGÍA
3.1 Base Teórica
El análisis de problemas de Mecánica de Fluidos se realiza mediante: la investigación experimental, el desarrollo teórico y la aplicación de la metodología de Dinámica de Fluidos Computacional (CFD: por sus siglas en inglés Computational Fluid Dinamic); dichos métodos, se relacionan y complementan entre sí.
En lo referente a la metodología CFD, se fundamenta en lo siguiente (Carrillo, 2014)
Ecuaciones generales de conservación a una partícula de fluido, cuyos principios y formulaciones son: principio de conservación de la masa, conservación de la cantidad de movimiento o segunda ley de Newton → Ecuaciones de Navier-Stokes (si se tienen en cuenta los términos viscosos y de disipación de energía) o ecuaciones de Euler (si se desprecian) y principio de conservación de la energía.
Ecuación de estado o de comportamiento del fluido.
Ecuaciones constitutivas del medio.
En lo que respecta al procedimiento para resolver un problema de mecánica de fluidos mediante la metodología CFD, se presentan los siguientes pasos (Blanco, 2007
Especificación de la geometría del problema.
Creación del mallado o celdas en las que van a ser calculadas todas las variables.
Definición de los modelos que se van a utilizar: modelos de turbulencia.
Especificación de las propiedades del fluido: viscosidad, densidad, propiedades térmicas, etc.
Imposición de las condiciones de contorno que controlan los valores de ciertas variables en los límites del dominio.
Introducción de las condiciones iniciales.
Control de los parámetros que afectan a la resolución numérica del problema.
Proceso de cálculo.
Análisis de la solución.
Por otro lado, es importante mencionar que los modelos bidimensionales (2D) estudian la variación de la velocidad promediada en vertical, considerando un campo de velocidades horizontal, es decir, cada uno de los puntos del dominio de estudio analiza el calado y las dos componentes de la velocidad.
Por otro lado, es importante mencionar que los modelos bidimensionales (2D) estudian la variación de la velocidad promediada en vertical, considerando un campo de velocidades horizontal, es decir, cada uno de los puntos del dominio de estudio analiza el calado y las dos componentes de la velocidad.
Los modelos 2D resuelven las ecuaciones que se fundamentan en los principios básicos de la hidrodinámica, como son: la conservación de masa y la conservación de la cantidad de movimiento. Las ecuaciones implementadas en estos modelos se rigen a estos principios y toman el nombre de las ecuaciones de Saint Venant para aguas someras, las mismas que resultan de la deducción matemática de las ecuaciones tridimensionales de Reynolds; estas últimas, obtenidas a partir de las ecuaciones de Navier Stokes.
Las hipótesis consideradas son las siguientes
La escala en profundidad es mucho menor a la escala horizontal.
Distribución hidrostática de presiones.
Por lo tanto, un modelo numérico bidimensional es aplicable para el análisis de las variables hidráulicas de un cauce, cuando
El cauce presenta meandros con extensas zonas de inundación.
El cauce presenta apreciables ensanchamientos y estrechamientos de sección que provoque zonas de recirculación y consecuente reducción de la sección hidráulica.
Se evidencia un campo de velocidades heterogéneo por sección transversal.
3.2 Modelación Unidimensional
En la aplicación de la modelación unidimensional, para el cálculo de las variables hidráulicas (calado y velocidad), se utilizó el software libre HEC-RAS 5.0, desarrollado por el Hydrologic Engineering Center del US Army Corps of Engineers, el cual, nos permite simular flujos en cauces naturales o canales artificiales.
La revisión de la información cartográfica para el área de estudio, se compone de: redes hidrográficas a escala 1:50000, curvas de nivel cada 5 metros en formato .shp y .dwg, modelo digital del terreno de tamaño del pixel 3 por 3 metros en formato .tiff, curvas de nivel del cráter cada 5 y 10 metros (IGM) y orto-fotografías a escala 1:50000 (MAGAP-PRAT, SIGTIERRAS).
El drenaje norte se compone de los sistemas hidráulicos representados por los ríos Pita y El Salto, en las cercanías al volcán; y por los ríos Santa Clara y San Pedro, ubicados a media y larga distancia del volcán, respectivamente. Con este antecedente, se examina la información referente a los hidrogramas de los ríos Pita, El Salto y Santa Clara; así como también, en la zona de amenaza del lahar hacia el flanco norte, se revisa los datos de campo de los vestigios de la erupción de 1877 (EMAAP & EPN, 2004).
Escenarios eruptivos
Considerando los dinamismos eruptivos, el estado y tamaño del glaciar, la magnitud de las erupciones precedentes y el comportamiento de otros volcanes análogos al Cotopaxi; se presenta cuatro posibles escenarios en el caso de una erupción del volcán Cotopaxi, ver Tabla 1.
Escenario |
VEI |
Tipo de Actividad |
Fenómenos Eruptivos |
Espesor Promedio Afectado (metros) |
1 (pequeño) |
1-2 |
- Fuente de lava - Estromboliana |
- Caídas de ceniza - Fragmentos balísticos - Flujos piroclásticos pequeños |
0.3 |
2 (moderado) |
2-3 |
- Estromboliana - Vulcaniana |
- Caídas de ceniza - Flujos piroclásticos - Flujos de lava |
1-2 |
3 (grande) |
3-4 |
- Vulcaniana a pliniana |
- Caídas de ceniza - Flujos piroclásticos |
4 |
4 (muy grande) |
Ñ>4 |
- Pliniana |
- Caídas de ceniza - Flujos piroclásticos |
8 |
En el estudio se consideró el Escenario 3, el mismo que se encuentra asociado a un Índice de Explosividad Volcánica 4 (VEI 4: por sus siglas en inglés). Este escenario se caracteriza por presentar una alta probabilidad de ocurrencia (56%), donde el evento es explosivo con emisiones importantes de ceniza, flujos piroclásticos por el desbordamiento del magma desde el cráter y la formación de lahares con derretimiento del glaciar (EMAAP & EPN, 2004).
La cuantificación del volumen de agua generado durante una erupción, se determina en base a: la proyección del tamaño del glaciar, área igual a 3.96 km2 (Cáceres B. , 2010) y la disminución promedio del espesor del glaciar, obtenido de los estudios realizados en otros volcanes, como: el nevado del Ruiz de Colombia, con 3 metros de espesor (Thouret, 1990) y el Monte Santa Elena de Estados Unidos, con 6 metros (Brugman & Meier, 1981).
El tamaño del glaciar del Cotopaxi ha sido estudiado partiendo del análisis de las cuencas glaciares, estudiadas a partir de fotografías aéreas y controles de campo. Para lo cual, se dividió el glaciar en diecinueve secciones, que suman 21.8 km2 (año 1976). El drenaje norte se compone de las primeras seis cuencas: Sindipamba, Carnero Machay, Potrerillos, Pucahuaycu Norte, Mudadero y Cajas (Jordan, 1976)
Estudios recientes sobre el estado del casquete glaciar, presentan una disminución importante del área debido a múltiples factores como el calentamiento global (Cáceres B. , 2010), por esta razón, en un evento eruptivo se produciría una significativa reducción del caudal del flujo de escombros.
Los hidrogramas de entrada son de forma triangular con el máximo valor en el primer tercio de la duración total del evento, en el orden de los 30 a 40 minutos (EMAAP & EPN, 2004). Para el ingreso de estos datos, se consideró las secciones de entrada “Hostería Tambopaxi” para el río Salto y “Hacienda San Rafael” para el río Pita; en la siguiente figura se muestran los hidrogramas de ingreso.
Entre febrero y marzo del año 2004, el personal del Instituto Geofísico de la Escuela Politécnica Nacional realizó el levantamiento de información de los vestigios de los lahares del evento eruptivo de junio de 1877, para el Estudio de amenaza volcánica del proyecto: “EVALUACIÓN DE IMPACTOS SOBRE LA INFRAESTRUCTURA DE LA EMAAP-Q Y ESTUDIO DE FACTIBILIDAD DE LAS OBRAS DE MITIGACIÓN ANTE UNA ERUPCIÓN DEL VOLCÁN COTOPAXI”. En este trabajo de campo se levantó información de calado, perímetro mojado, velocidad media y caudal pico, en: 26 secciones del río Pita, 10 secciones del río Salto, 8 secciones del río Santa Clara y 13 secciones del río San Pedro.
En lo referente a la calibración, ésta se lleva a cabo con el modelo unidimensional en flujo permanente, en base a los datos de campo recopilados a partir de los vestigios del evento eruptivo de 1877 (Escenario 3). Los parámetros de calibración corresponden a los calados o profundidades de flujo, asociados a los caudales máximos de crecida en cada tramo de análisis.
El tramo de simulación unidimensional es de aproximadamente 50 km de longitud, desde la sección “Hostería Tambopaxi” hasta la sección “La Armenia” en el río San Pedro, este último se encuentra aguas abajo de la junta de los ríos Santa Clara y Pita en el sector de “El Triángulo” en Sangolquí, ver Figura 2.
3.3 Modelación bidimensional
La modelación numérica bidimensional se realizó en la zona densamente poblada de Sangolquí. Para esta modelación se utilizó los programas computacionales IBER 2.3 y FLOW 3D 11.1. IBER es un software libre bidimensional desarrollado por el Grupo de Ingeniería del Agua y del Medio Ambiente, GEAMA (Universidad de A Coruña, UDC) y el Instituto FLUMEN (Universitat Politècnica de Catalunya, UPC, y Centro Internacional de Métodos Numéricos en Ingeniería, CIMNE), y se aplica a modelos numéricos en hidráulica y morfología fluvial. Por otro lado, FLOW-3D es un software de simulación de fluidos computacional CFD, creado por FlowScience y aplicado principalmente al cálculo de fluidos a superficie libre, el mismo que dispone de la opción Shallow Water (aguas poco profundas), que de igual manera, realiza el cálculo como un modelo bidimensional.
El tramo de esta modelación comprende 36.36 km2, dentro de los cuales se encuentran los ríos: Pita 7.83 km, Santa Clara 5.1 km y San Pedro 3.8 km.
En las modelaciones se aplica el siguiente procedimiento
Para el ingreso de la topografía a los CFDs, se extrae el tramo a simular del modelo digital del terreno (mdt), aplicando un software de sistemas de información geográfica (SIG) y se transforma esta información a formato ASCII.
Se genera un mallado estructurado hexaédrico de 9 por 9 metros, donde cada malla presenta al final un único hidrograma de salida. Se crea tres bloques de malla con planos de ingreso de caudal en los ríos Pita y Santa Clara, y con salida de flujo en el río San Pedro de acuerdo a los drenajes presentes en el tramo.
La reología del flujo de escombros se representa a través de: el coeficiente de rugosidad de Manning (n) en el programa IBER y con la rugosidad absoluta en el programa FLOW 3D. Estos coeficientes de resistencia al movimiento además de ser un parámetro que representa la superficie del cauce, representa el comportamiento de la mezcla agua-escombros.
En el modelo IBER se utilizó un coeficiente de rugosidad de Manning (n) de 0.063, resultado de la ponderación de los coeficientes de resistencia al movimiento, según la longitud del río (EMAAP & EPN, 2004). En el programa FLOW 3D, se utiliza un valor de rugosidad absoluta de 0.216 metros, obtenido del estudio: “PREPARACIÓN DE LA BASE GEOMORFOLÓGICA PARA LA SIMULACIÓN NUMÉRICA TRIDIMENSIONAL DEL FLUJO DE LAHARES DEL VOLCÁN COTOPAXI, FLANCO NORTE", realizado en el año 2016 en la Escuela Politécnica Nacional. (BID & EPN, 2016)
El ingreso de fluido se configura como agua, con el vector de dirección en cada río. Y los hidrogramas de ingreso se obtuvieron de la modelación unidimensional.
En el análisis bidimensional se consideró cinco secciones, cuyas coordenadas se indican en la siguiente tabla
Sección |
Río |
Coordenadas UTM |
|||
Margen Izquierda |
Margen Derecha |
||||
X |
Y |
X |
Y |
||
1 |
Santa Clara |
783262.91 |
9963579.21 |
784742.28 |
9964419.89 |
2 |
Santa Clara-Pita |
782659.63 |
9965556.78 |
785281.67 |
9966198.91 |
3 |
San Pedro |
781566.51 |
9968302.79 |
784785.63 |
9967753.68 |
4 |
Pita |
786203.51 |
9963617.81 |
786479.05 |
9965408.76 |
5 |
Pita |
784765.48 |
9964692.89 |
785493.63 |
9966017.61 |
4. RESULTADOS
4.1 Modelación Unidimensional
La calibración del modelo unidimensional se alcanza representando los niveles de vestigios del evento de 1877, cuyos coeficientes de resistencia al movimiento (n) varían entre 0.01 y 0.09.
El análisis del flujo no permanente permite determinar los caudales de crecida en las secciones de entrada y salida del modelo numérico, los mismos que, debido a la laminación del flujo, van disminuyendo conforme transita el caudal a través del cauce.
De igual manera, en la siguiente tabla se presenta el tiempo en el cual se registra el caudal máximo y la duración del hidrograma. Además, se calculan los calados con el modelo unidimensional obteniéndose valores que fluctúan entre 23 y 49 metros de altura.
Tramo |
Tipo |
Caudal máximo m3/s |
Tiempo (minutos) |
Duración (minutos) |
Pita 1 |
Entrada |
23814.0 |
13 |
40 |
Pita 1 |
Salida |
22769.7 |
17 |
60 |
Salto 1 |
Entrada |
15708.5 |
10 |
30 |
Salto 1 |
Salida |
14228.0 |
20 |
60 |
Pita 2 |
Entrada |
36997.7 |
20 |
60 |
Pita 2 |
Salida |
36692.3 |
21 |
60 |
Santa Clara 1 |
Entrada |
11007.7 |
21 |
60 |
Santa Clara 1 |
Salida |
9771.1 |
40 |
80 |
Sangolquí-Ejido |
Entrada |
10306.9 |
35 |
60 |
Pita 3 |
Entrada |
25684.6 |
21 |
60 |
Pita 3 |
Salida |
22056.7 |
38 |
80 |
Colina |
Entrada |
22056.7 |
38 |
80 |
4.2 Modelación Bidimensional
Alturas Máximas
Tomando en consideración las modelaciones en IBER y FLOW 3D, en la siguiente tabla se compara para las secciones de análisis, los resultados de calados máximos; donde se obtiene un porcentaje de coincidencia de 90% en la margen izquierda y 85% en la margen derecha del río.
Sección |
Río |
Calado (IBER) |
|
Margen Izquierda |
Margen Derecha |
||
1 |
Santa Clara |
11 |
16 |
2 |
Santa Clara-Pita |
5 |
15 |
3 |
San Pedro |
19 |
17 |
4 |
Pita |
18 |
12 |
5 |
Pita |
16 |
15 |
Sección |
Río |
Calado (FLOW-3D) |
|
Margen Izquierda |
Margen Derecha |
||
1 |
Santa Clara |
11 |
15 |
2 |
Santa Clara-Pita |
4 |
12 |
3 |
San Pedro |
19 |
18 |
4 |
Pita |
15 |
7 |
5 |
Pita |
18 |
15 |
Sección |
Río |
Porcentaje de Coincidencia |
|
Margen Izquierda |
Margen Derecha |
||
1 |
Santa Clara |
100 |
94 |
2 |
Santa Clara-Pita |
80 |
80 |
3 |
San Pedro |
100 |
94 |
4 |
Pita |
83 |
58 |
5 |
Pita |
89 |
100 |
Promedio |
90 |
85 |
Velocidades Máximas
Considerando las cinco secciones de análisis, en la siguiente tabla se observa que los resultados de las velocidades máximas obtenidas en las modelaciones con los dos programas, son muy similares.
Hidrograma de Salida
En lo referente a los hidrogramas de salida, en la siguiente figura se muestran los resultados obtenidos donde se identifica el caudal pico:
Por otro lado, se verifica los volúmenes de ingreso y salida de las modelaciones. En la siguiente tabla se muestran los valores obtenidos:
4. CONCLUSIONES
El coeficiente de rugosidad “n” de Manning, que en este caso, representa las condiciones de resistencia al movimiento, tanto del lahar con el contorno, como las internas de la mezcla; es el parámetro que permitió realizar la calibración del modelo unidimensional, al comparar los calados obtenidos en la modelación, con los vestigios del evento eruptivo de 1877.
Se determinó que al utilizar un mallado estructurado hexaédrico de 9 por 9 metros, se obtienen resultados aceptables, en comparación al mallado más fino de 3 por 3 metros; lo cual permite: simplificar el modelo, reducir los recursos computacionales y los tiempos de cálculo.
Para el evento de mayor probabilidad de ocurrencia, escenario 3, al considerar el área del casquete glaciar al año 2010, se verificó que la afectación en el caso que se produzca un evento eruptivo, es menor en comparación a los resultados obtenidos en estudios anteriores.
Con base en la comparación de parámetros como: el calado, la velocidad y el área de afectación; se concluye que los resultados obtenidos al aplicar los modelos IBER y FLOW 3D, son muy similares. En general, considerando el calado de los dos modelos, se identificó una coincidencia de: 90% en la margen izquierda y 85% en la margen derecha del cauce.
Respecto a los resultados de velocidades máximas, existe un porcentaje de coincidencia de 98%, cuyo valor de velocidad máxima para la sección representativa es de 21 m/s en las dos modelaciones.
El caudal pico del hidrograma de salida es de 21.388 m/s a los 60 minutos y 23.490 a los 52 minutos, para las modelaciones con IBER Y FLOW 3D, respectivamente; alcanzando un 92% de similitud. El volumen de salida es de 43.17 millones de metros cúbicos en el modelo IBER y 44.31 millones de metros cúbicos en el modelo FLOW 3D, coincidiendo en un 97%.
En lo referente a la comparación de las áreas de afectación, estas son similares y aceptables, el área de inundación que resulta de la modelación numérica en IBER es de 8.26 km2 y en el FLOW 3D es de 7.72 km2, presentando un 93.5% de coincidencia.
Se recomienda el uso del modelo unidimensional en zonas con altas pendientes de fondo o encañonadas (encauzadas), es decir, donde predomine la dirección longitudinal del flujo del río, parar lograr la representación del modelo numérico de manera satisfactoria y optimizar recursos.
En cambio, en zonas susceptibles a inundación, como es el caso de la zona densamente poblada de Sangolquí, se recomienda el uso de un modelo bidimensional, que permita representar el flujo en dos direcciones y consiga una caracterización adecuada del fenómeno.
Con la finalidad de obtener resultados fiables en tiempos razonables, al realizar la modelación numérica bidimensional y tridimensional, se sugiere emplear computadoras de gran capacidad de procesamiento y almacenamiento. Se recomienda las siguientes características: Core I7 3.4 GHZ, 8 núcleos, 32 GB RAM, espacio disponible en disco 200 GB SSD, tarjeta gráfica de video 2 GB; o superior.