Introducción
El vulcanismo es uno de los fenómenos geológicos más devastadores a los que se ha enfrentado el hombre. En particular, por la naturaleza y duración de las erupciones volcánicas, las cuales pueden llevar a veces a evacuar y relocalizar poblaciones enteras, con todas las implicaciones sociales, políticas y ambientales que esto conlleva. En Colombia, se conocen, hasta la fecha, 23 volcanes activos, localizados principalmente en la cordillera Central y cerca a varias de las regiones más pobladas y estratégicas del país, como el Eje Cafetero, la zona del Cauca y Nariño, entre otras.
Entre los volcanes más explosivos y peligrosos de Colombia, se encuentra el Volcán Cerro Machín (VCM), un volcán con composición geoquímica félsica, el cual ha presentado actividad explosiva en los últimos 5 ka. Su última erupción importante se registró hace 900 años (Méndez etal., 2002). Se estima que su ciclo de actividad eruptiva violenta se produce cada 900 años, razón por la cual es de vital importancia estudiar este volcán en detalle, ya que una erupción en la actualidad del VCM podría causar una gran tragedia ambiental, social y económica para Colombia. También se estima que el occidente de Colombia quedaría aislado de la parte central y oriental del país, debido a que la vía terrestre conocida como La Línea sufriría enormes afectaciones. Por otra parte, el municipio de Cajamarca, con cerca de 15.000 personas, podría desaparecer en minutos si una gran erupción del VCM ocurriera. Además, gran cantidad de municipios del Tolima y algunos de Cundinamarca se verían afectados por lahares que se generarían durante y después de la erupción. Este panorama de riesgo volcánico indica que es de gran importancia conocer en el mayor grado posible el comportamiento del VCM desde diferentes puntos de vista y disciplinas.
Desde que Gutenberg y Richter (1954) establecieron la relación frecuencia-magnitud (FM) de los sismos (Log10N=a-bM; donde N= número de sismos acumulados con magnitud mayor o igual a M, a= productividad sísmica, b= pendiente o valor b), esta ha sido usada ininterrumpidamente para obtener información sobre la actividad sísmica de una región. Con el pasar del tiempo, se pudo determinar que esta relación brinda más información que la simple actividad sísmica de una región. Varios estudios han demostrado que cambios en las condiciones del medio como la presión de poros, variaciones en la temperatura, acumulación de esfuerzos, entre otros, hacen que el parámetro conocido como valor b, o la pendiente de la relación FM, varíe (Scholz, 1968; Warren y Latham, 1970; Wyss et al., 1997; Wiemer y Wyss, 2000, entre otros). El valor de 1 para b ha sido considerado el valor normal de la relación entre sismos grandes y pequeños, mientras que valores por encima o debajo de 1 se consideran anómalos. La distribución espacial del valor b ha sido usada como ayuda para determinar zonas donde este parámetro presenta valores anómalos. Estas zonas permiten identificar parte de la estructura interna de una región. En particular, en zonas volcánicas, el mapeo o distribución espacial del valor b ha permitido definir zonas con valores anómalos mayores a 1, que posiblemente están asociadas con acumulación de magma o fluidos (Wyss et al., 1997; Wiemer y Wyss, 1997; Wyss et al., 2001; McNutt, 2005; Sánchez et al., 2005; Londoño y Rodríguez, 2013; Chiba y Shimizu, 2018). Por su parte, se han reportado cambios temporales en el valor b antes de la ocurrencia de sismos (Smith, 1986; Enescu e Ito, 2001; Tsukakoshi y Shimazaki, 2008; Mallika et al., 2013; Zheng y Zhou, 2014; Chen y Zhu, 2020) y erupciones volcánicas (Patané et al., 1992; Nishimura et al., 2016). Esto demuestra la utilidad del valor b como herramienta de monitoreo sísmico y volcánico, aunque hasta la fecha no se ha definido oficialmente por la comunidad científica como un premonitorio de sismos o erupciones volcánicas.
Algunos estudios sísmicos se han realizado para el VCM (Figura 1), el cual ha registrado sismicidad desde que se inició su monitoreo sísmico continuo en 1989 (Londoño, 2004, 2012). Sin embargo, a la fecha no se ha realizado un estudio de la variación espacio-temporal del valor b en el VCM. Este trabajo busca realizar un análisis de la variación espacio-temporal del valor b en el VCM y establecer relaciones de esas variaciones con su actividad para el período 2007-2020.
Actividad reciente del VCM
El VCM localizado en la cordillera Central de Colombia (Figura 1) es catalogado como uno de los volcanes más explosivos del país (Cepeda y Murcia, 2000). Se han podido identificar al menos seis erupciones importantes en los últimos 5000 años, la última hace 900 años (Rueda et al., 2005). El monitoreo continuo del VCM se inició en 1989, y posteriormente en 2004 se diversificó la red de monitoreo, con diferentes estaciones instaladas en campo, que cubren disciplinas como sismología, geodesia, geoquímica y geofísica de campos potenciales.
En los últimos años, el VCM ha presentado incrementos en su actividad, particularmente en la actividad sísmica, deformación y geoquímica, lo cual evidencia cambios en los sistemas hidrotermal y magmático (Inguaggiato et al., 2017). Se han observado cambios en la deformación en cercanías al domo principal, con una subsidencia o deflación hacia el SSW, mientras que el resto de zonas han permanecido estables sin signos de deformación.
Entre 2008 y principios de 2013, el volcán mostró la mayor actividad sísmica registrada hasta la fecha en términos de energía y número de sismos; se destaca el enjambre sísmico ocurrido en noviembre de 2008 hacia el NE y E del domo principal (INGEOMINAS, 2008), en el cual se registró el sismo de mayor magnitud observado hasta la fecha, ocurrido el 9 de noviembre (ML 4,9, prof=3,0 km, a 1 km al E del domo principal). Así mismo, durante ese período, se registró sismicidad distal (4-10 km) hacia el SE del domo principal, particularmente entre 2010 y principios de 2013 (julio 2010 a abril de 2011 y octubre 2012 a febrero de 2013), con profundidades entre 9 y 18 km, la cual no se había registrado anteriormente.
Así mismo, en enero de 2013, se pudo detectar geoquímicamente la entrada de un posible nuevo pulso de magma en la parte profunda del volcán, relacionada con un cambio en la firma isotópica en la relación R/Ra vs log(C/3He). Esta entrada de material magmático se evidenció al cambiar la relación R/Ra de 3,6 en noviembre de 2011 a 5,87 en enero de 2013. Simultáneamente, la relación log(C/3He) pasó de 57% a 80% en esas mismas fechas, lo que indica un origen mantélico de dicho incremento (Inguaggiato et al., 2017).
En febrero de 2018, ocurrió un enjambre sísmico hacia el SW del domo principal, entre 3 y 5 km de profundidad, con el mayor número de sismos diarios registrado hasta la fecha, aunque no tan energéticos como en el período 2008-2012. La sismicidad de tipo volcano-tectónica, asociada con fracturamiento y agrietamiento de las rocas al interior del volcán, es la que predomina en el VCM (Figura 2), aunque en el período 2018-2019 se han registrado unas pocas señales sísmicas que podrían estar asociadas a actividad de fluidos dentro del volcán.
Metodología
Para calcular el valor b (b), se utilizó el método de máxima verosimilitud (Shi y Bolt, 1982; Aki, 1965; Utsu, 1965), mediante la expresión:
Donde es la magnitud promedio y M m¡n es la magnitud mínima de una muestra de datos, que puede expresarse en función de la magnitud de completitud, Mc. Así mismo, Bender (1983) y Utsu (1999), entre otros, establecieron un AMcorrespondiente a un "bin" o rango de magnitudes agrupadas o acumuladas, que por lo general tiene un valor de 0,1, modificando la fórmula de Aki (1965) y Utsu (1965) así:
De acuerdo con Aki (1965), el valor b obtenido por la fórmula de Utsu (1965) (ecuación 1) es la estima con máxima verosimilitud; este autor también propuso los límites de confianza de dicha estimación para diferentes muestras, estableciendo que:
Donde b'=b/log10e, de=1,96 para una probabilidad (e) del 95% o un nivel de confianza del 95%. Aki (1965) define además los valores de de/√n para diferentes niveles de confianza y muestras.
Para determinar la porción de la curva FM, donde esta toma una forma lineal, y sobre la cual en la práctica se calcula el valor de la magnitud de completitud, Mc, de la que depende en gran medida el valor b, existen diferentes métodos (Woessner y Wiemer, 2005; Mignan y Woessner, 2012). En este trabajo se utilizó el método de máxima curvatura, MaxC (Wiemer y Wyss, 2000), que se basa en cambios en la distribución de FM o cambios en la pendiente de dicha curva. Este método proporciona un valor de Mc confiable de manera rápida, y busca el máximo punto de curvatura, el cual se define como el valor máximo de la primera derivada de la curva FM, que coincide con el bin de magnitud que contiene la más alta frecuencia de eventos en la curva FM no acumulada. Este método, por lo general, produce un valor subestimado de Mc, ya que puede omitir curvaturas graduales en la curva FM debido a las posibles heterogeneidades espaciales o temporales. No obstante, según Mignan y Woessner (2012), si se eliminan dichas heterogeneidades, la técnica de MaxC no subestima los valores de Mc.
Para mapear la distribución espacial del valor b, se utilizó el método de Wiemer y Wyss (1997), el cual consiste en discretizar la zona de estudio mediante una retícula, donde en cada nodo se define un cilindro vertical y se calcula el valor b usando un radio variable hasta alcanzar un número determinado de sismos (muestras) o mediante un radio y un número fijo de sismos. En este caso, se utilizó una retícula de 0,5 x 0,5 x 0,5 km, con un radio variable hasta encontrar 300 sismos para la distribución del valor b en perfil o profundidad, y se usaron estos mismos criterios, pero sin importar la profundidad del sismo, para el mapa del valor b en planta. Debido a la diferencia de calcular el valor b en planta y perfil, es posible observar discrepancias en cercanías a la superficie entre los valores de b en planta y perfil.
Para la variación temporal del valor b, se tomaron ventanas de tiempo de 90 días, con traslape de 30 días, y en cada una de ellas se calculó el valor b y su error estándar. Se elaboraron gráficas y mapas de distribución temporal y espacial del valor b para el VCM.
El catálogo sísmico fue facilitado por el Observatorio Vulcanológico y Sismológico de Manizales (OVSM) del Servicio Geológico Colombiano (SGC). Dicho catálogo de sismos abarca el período enero de 2007 a diciembre de 2019 y contiene un total de 12.155 sismos. Los sismos fueron localizados con el algoritmo HYPO71PC (Lee y Valdez, 1985) y un modelo de corteza propuesto por Londoño (2004). La Figura 3 muestra los hipocentros usados para el estudio. La magnitud local se calculó a partir de una fórmula derivada para la zona de estudio por Londoño (2016).
Resultados
Antes de analizar los resultados obtenidos del cálculo espacial del valor b es necesario validarlos. Para ello, se elaboraron mapas del valor b superponiendo contornos de errores, número de sismos por nodo y de radio necesario (en km) para incluir 300 sismos en un nodo. Las Figuras 4, 5 y 6 muestran los resultados de la distribución espacial del valor b, su error, número de sismos y radio para cada nodo. En estas figuras se puede observar que los errores de los valores de b, en general, son bajos, en su mayoría menores a 0,07, si se toma como referencia un valor de 1,0 (Wyss et al., 2001), aunque se pueden observar algunos valores más altos hasta máximo 0,23 hacia la parte SW del domo principal (Figura 4).
Por su parte, la distribución del número de sismos por nodo muestra que todos los nodos que cubren la zona del volcán presentan 300 sismos o más (Figura 5) en un radio menor a 1,5 km, tomado desde las coordenadas del nodo. Para los extremos NE y SW del volcán, no se obtuvieron buenos resultados del valor b, por lo que se consideran zonas con mala resolución, con menos de 200 sismos por nodo y radios mayores a 3 km para obtener un mínimo de 300 sismos ( Figura 6). Estos resultados demuestran que los valores obtenidos del valor b para la zona que cubre el edificio del VCM son confiables y realistas y no se deben a artefactos matemáticos.
Con miras a validar la confiabilidad e independencia de los valores b calculados espacialmente, se seleccionaron aleatoriamente tres nodos a los cuales se les calculó el valor b y se compararon las muestras tomadas para dicho cálculo. La Figura 7 muestra los resultados. De esta figura es posible concluir que los tres valores son independientes y que las muestras usadas para el cálculo también lo son.
En la Figura 4 se puede observar que, en general, los valores de b para el VCM son bajos (<0,7) para las zonas del domo principal y hacia el SE del mismo, mientras que se observan valores relativamente altos para el SW del domo principal y para la zona SE más distal del edificio, cerca al domo de Tapias. La Figura 8 muestra un perfil NW-SE cruzando el domo principal y el domo de Tapias, donde se observan los valores de b en profundidad. Para el cálculo de este perfil se usó la misma retícula de 0,5 x 0,5 x 0,5 km. De esta figura es posible observar que el VCM presenta tres zonas anómalas de valores b en profundidad; una localizada entre -1 y 1 km de profundidad al SW del domo (b>1,3±0,13); otra debajo del domo central entre 0 y 1 km (b>1,3±0,13) y otra entre 8 y 13 km de profundidad debajo del domo Tapias (1,2±0,1).
La Figura 9 muestra la variación temporal de la magnitud de completitud (Mc) y del valor en el VCM entre 2007 y principios de 2020. Mc, en general, varió entre -1 y 0,7, con promedio cercano a cero. Se puede observar también que el valor b presentó cambios temporales; en particular, tres cambios importantes se pudieron detectar. El primero (1 en Figura 9), en el período diciembre 2007 a febrero de 2008, que varió de 0,8 antes de ese período hasta 1,4 durante ese período, luego descendió hasta 0,55 antes del sismo de ML 4,9 el 9 de noviembre de 2008. El segundo (2 en Figura 9), ocurrido en el período octubre a la primera mitad de noviembre de 2010, varió de 0,8 antes de ese período a 1,5 durante ese período, justo después de la aparición de sismicidad profunda (9-18 km) localizada entre 4 y 10 km hacia el SE del domo principal. El tercero (3 en Figura 9), ocurrido en el mes de febrero de 2018, se asocia con un aumento importante de la sismicidad del VCM (Figura 2), variando los valores de b de 1,0 a finales de enero de 2017 a valores de 2,5 a mediados de febrero de 2018, luego de lo cual retornó a valores de 1,0 a finales del mismo mes. Entre noviembre de 2013 y noviembre de 2017, el valor b mostró inestabilidad, con una gran variación oscilatoria (0,6 a 1,5), sin mostrar una tendencia clara. Así mismo, también es posible observar en la Figura 9 que, a partir de octubre de 2013, el valor b aumentó su promedio, mientras que la Mc permaneció similar, por lo que es posible separar en dos grandes períodos los valores de b en el VCM para la ventana de tiempo analizada: entre enero de 2007 y septiembre de 2013, y entre octubre de 2013 y enero de 2020.
Análisis e interpretación
Para este estudio, no se realizó desagrupamiento (declustering) del catálogo sísmico, dado que la sismicidad del VCM básicamente se produce a manera de enjambres sísmicos, fenómeno que es intrínseco a un posible proceso de intrusión que puede estar presentando el VCM, y que no está relacionado con la historia sísmica de la zona (Van Stiphout et al., 2012).
Por otra parte, de acuerdo con los resultados de las pruebas realizadas sobre la independencia de los valores de b obtenidos en cada nodo, se asume que los valores de b obtenidos en este trabajo son independientes en cada nodo, lo que permite observar variaciones espaciales que se pueden originar en los cambios reales en el valor b, más que en efectos estadísticos o de repetición de datos.
Los resultados obtenidos en este trabajo permiten concluir que el valor b en la zona del VCM ha experimentado cambios tanto espaciales como temporales en el lapso estudiado (2007-enero 2020). Estos cambios podrían obedecer a variaciones en la actividad del VCM. Las tres zonas con valores anómalos altos de b (>1,2) en profundidad, localizadas SW (-1, 1 km), debajo del domo principal (0, -1 km) y al SE del domo central (8-13 km), permiten suponer que es posible que estén asociadas a parte de la estructura interna del volcán y sus alrededores (Figura 8). La zona al SW más superficial podría estar asociada con material fluido y caliente, posiblemente perteneciente al sistema hidrotermal del VCM, y que es en parte responsable de la actividad fumarólica y de fuentes termales con relativa alta temperatura que se encuentran en cercanías del domo central (Inguaggiato et al., 2017). Por su parte, la zona anómala al SE, más profunda, podría estar relacionada con parte de uno de los reservorios magmáticos del VCM. La zona anómala debajo del domo central se podría relacionar con parte del sistema magmàtico superficial del VCM, posiblemente con parte de la estructura interna del domo central. Londoño (2012) realizó una tomografia sísmica 3D de velocidad para el VCM y definió una zona de baja velocidad entre 10 y 12 km de profundidad, coincidente con la zona anómala de valores de b hacia el SE del VCM, en cercanías al domo Tapias.
Por otra parte, las variaciones temporales del valor b estarían indicando cambios en la actividad del volcán a través del tiempo. Los tres rangos de fechas en los que el aumento del valor b es más claro, descritos anteriormente y observados en la Figura 9, corresponden a períodos de cambios en la actividad del volcán, lo que evidencia la utilidad del valor b como herramienta rutinaria de monitoreo volcánico.
En el primer período (rotulado 1 en la Figura 9), en noviembre de 2008, como se mencionó anteriormente, ocurrió el sismo de mayor magnitud registrado hasta la fecha en el VCM. Este sismo ocurrió dentro de un enjambre sísmico importante localizado en el domo principal y sus alrededores (Figura 10), que registró el segundo pico mayor en número de sismos hasta la fecha en el VCM. El sismo se localizó a 1,5 km al NE del domo principal y a 3 km de profundidad.
Es interesante observar, inicialmente, el aumento importante del valor b (>1,4), en diciembre de 2007 y enero de 2008, es decir, entre 10 y 11 meses antes del enjambre, y posteriormente la disminución importante (<0,6) justo un mes antes del enjambre (Figura 9), lo que podría estar indicando una posible pequeña intrusión de magma entre abril y mayo de 2008, que posiblemente pudo iniciar su ascenso sin que se registrara sismicidad en grandes proporciones ni deformación (Figuras 2 y 10); esta posteriormente hizo que se generara un aumento en la presión de poros dentro del edificio volcánico en la parte más superficial, lo cual disminuyó drásticamente el valor b (<0,6) y generó luego el fuerte enjambre sísmico que incluyó el sismo de mayor magnitud registrado hasta la fecha, el cual tuvo manifestaciones superficiales importantes como la salida de fuentes de agua en lugares donde normalmente no se presentaba, agrietamiento del terreno y varios deslizamientos de tierra.
Así mismo, se detectaron cambios geoquímicos antes del enjambre (INGEOMINAS, 2008). Por lo tanto, es posible suponer que este primer cambio importante del valor b estuvo asociado a una intrusión magmática, posiblemente en forma de dique o silo, sugerida además por la distribución de los hipocentros (ver más adelante).
El segundo aumento importante del valor b sucedió aproximadamente dos meses después de la ocurrencia de una sismicidad profunda (12-18 km), que inició a mediados del mes de agosto de 2010 y finalizó en diciembre de 2010, localizada entre 4 y 10 km al SE del domo principal, cubría parte del domo de Tapias y la cual no se había presentado antes en el VCM desde que se tiene registro instrumental, el cual empezó en 1988 (Figura 10). Si bien es cierto que los valores de b eran menores a 1,0 antes del incremento, no se observó una tendencia a disminuir de manera importante, como sí sucedió en el primer incremento discutido anteriormente, razón por la cual probablemente no se generó un sismo de tamaño similar al del primer incremento.
Es posible que este segundo incremento temporal del valor b pueda estar asociado a un nuevo ascenso de un pulso de magma, tal vez más profundo que el primero, pero ciertamente localizado hacia el SE del VCM a diferencia del primero que estaría más cerca del domo principal, o a un movimiento de una falla enterrada que no se observa en superficie. La primera opción es más factible que la segunda, ya que la sismicidad inició de manera distal y con el tiempo fue acercándose al domo principal. Incluso de forma simultánea ocurría sismicidad tanto en el domo principal como en la parte distal (Figura 10). Este comportamiento sísmico y distribución de la sismicidad, la cual cambia o rebota de una zona a otra, es muy común en procesos de intrusión de diques en zonas volcánicas (Savage y Cockerham, 1984; Morishita et al., 2016; Woods et al., 2019; Bonaccorso y Giampiccolo, 2020; Koike y Nakamichi, 2021).
Adicionalmente, después de la ocurrencia de esta sismicidad distal de 2010, se observó un cambio geoquímico importante con un aumento en isótopos en la relación R/Ra vs log(C/3He), asociado con una inyección de magma profunda, entre 2011 y 2013, como se mencionó anteriormente (Figura 9) (Inguaggiato et al., 2017); esto soporta la idea de una intrusión magmática como la responsable de este cambio en el valor b. Es notable que a esa profundidad se observa una zona de valores altos de b (1,1-1,3; Figura 4).
Por último, es interesante observar que justo después del cambio en la firma isotópica geoquímica mencionada anteriormente, el promedio del valor b cambia, aumentando su valor. Al comparar los valores antes y después de esa fecha (2013), se observa una variación del promedio general (Figura 9). Este comportamiento del valor b podría indicar el efecto de la inyección magmática profunda tanto en la sismicidad como en el medio del VCM, la cual posiblemente produjo un aumento en la temperatura al interior del volcán. Otra posibilidad es un aumento en la presión de poros debido a dicha intrusión magmática, lo que genera un aumento en el valor b (Wiemer et al., 1998). Así mismo, se observó un cambio en la pendiente de la curva de deformación de uno de los inclinómetros electrónicos instalados en el domo principal (estación CIMA en Figura 2C). Esta estación empezó a mostrar cambios importantes a partir de junio de 2012 en la pendiente de la curva, asociados con una deflación hacia el SSW (Figura 2C).
El tercer incremento del valor b se dio en medio de un enjambre con el mayor número de sismos que se haya registrado en un día en el VCM hasta la fecha (enero 2022), con más de 1000 sismos en un día, ocurrido entre el 17 y 27 de febrero de 2018 (SGC, 2018). Este evento fue similar al ocurrido en noviembre de 2008, aunque un poco menos energético (Figura 3 ). Esto lleva a suponer que muy posiblemente también se haya tratado de una intrusión de un pulso de magma ubicado en el mismo sector que en 2008, pero localizado un poco más hacia el SW del domo principal, definiendo una posible zona de intrusión superficial para el VCM con dirección SW-NE y que cruza el domo principal (Figura 10). Este aumento del valor b se dio en medio de un proceso de deformación en la zona del domo principal, el cual se observa en la Figura 2C, luego del cual la curva del vector resultante del inclinómetro electrónico CIMA muestra un cambio en la pendiente. Esta observación puede tener implicaciones en la evaluación de la amenaza volcánica del VCM, ya que la posible presencia de una intrusión de un dique o un silo relativamente superficial puede aumentar la probabilidad de una erupción en el VCM.
Si bien los valores de b presentan cierta incertidumbre intrínseca, es posible usarlo como una herramienta rutinaria de monitoreo volcánico, a la vez que su distribución espacial ayuda a delimitar posibles zonas relacionadas con reservorios magmáticos o fluidos hidrotermales. En este trabajo se muestra la utilidad del cálculo del valor b como apoyo a la evaluación de la evolución de la actividad del VCM.
Conclusiones
Se llevó a cabo un estudio de la variación espacio-temporal del valor b en el VCM para el período 20072020. Se pudieron detectar cambios espaciales y temporales de dicho valor, los cuales se han asociado con cambios en la actividad del VCM. En particular se detectaron tres episodios temporales de aumento del valor b que pueden estar asociados a tres posibles intrusiones de pulsos de magma en 2008, 2010 y 2018.
Adicionalmente, se pudieron identificar dos zonas con valores anómalos altos del valor b (>1,2) hacia el SW del domo principal, entre -1 y 1 km de profundidad, y en la parte SE distal del domo principal, en cercanías al domo de Tapias, entre 8 y 13 km de profundidad, las cuales pueden estar asociadas, la primera, con la presencia de fluidos hidrotermales y, la segunda, con un reservorio de magma.
El cálculo rutinario del valor b sirve como herramienta de monitoreo para el VCM, que ayuda a entender mejor su comportamiento y a realizar un mejor diagnóstico de su estado de actividad.