Tomografía sísmica 3D del nido sísmico de Bucaramanga (Colombia)
3D seismic tomography in the seismic
nest of Bucaramanga (Colombia)
Francia
Jullyette Sepúlveda-Jaimes
fj.sepulveda06@gmail.com
Francisco
Henry Cabrera-Zambrano
fcabrera61@yahoo.com.br
Resumen
La tomografía
sísmica es una técnica con la cual se genera una imagen del interior de la
tierra; una de las posibilidades de datos a utilizar son los tiempos de llegada
de las ondas generadas en un sismo a las estaciones sismológicas. En el presente
trabajo se utilizó dicha técnica para generar un modelo de velocidades cuasi 3D
de la Onda P, en la zona del nido de Bucaramanga, en donde se presenta una alta
actividad sísmica, posiblemente debido a la colisión de la Placa de Nazca con
la Placa Caribe bajo la Placa Suramericana, que da lugar a la ocurrencia de
sismos a profundidades intermedias (70 – 200 km). Para el desarrollo de este
trabajo, se tuvo en cuenta un total de 1223 sismos ocurridos en el departamento
de Santander – Colombia; entre los 6° - 8° N y 72° - 74°W, en el periodo de Enero de 2012 y Junio de 2016 registrados por las estaciones
sísmicas de la Red Sismológica Nacional de Colombia (RSNC) distribuidas a lo
largo de todo el país. De estos eventos se utilizaron 1138 los cuales cumplían
con una serie de parámetros asociados al nido sísmico –tales como latitud,
longitud y profundidad – para generar un primer modelo de velocidades 1D, en el
que se asumió una homogeneidad lateral y solo se tuvo en cuenta la variación de
la velocidad a profundidad. Finalmente, se generó el modelo de velocidades
cuasi 3D, es decir, un modelo por capas horizontales a distintas profundidades,
con variación lateral en cada plano, con el cuál se determinaron dos anomalías
que se hacen más evidentes a partir de los 30 kilómetros de profundidad: la
primera es una anomalía negativa que se encuentra posiblemente asociada a la
fusión parcial del magma, y la segunda, es una anomalía positiva que se
interpreta como asociada a la cristalización del magma.
Palabras
clave: tomografía sísmica, nido sísmico de Bucaramanga,
modelo de velocidades.
Abstract
Seismic tomography is a set of
techniques used to obtain an image of velocity fields inside the Earth, using
waves measurements. In this work, this technique was used to generate a
quasi-3D velocity model of the P wave of Bucaramanga’s nest zone –a highly
seismic zone in Colombia caused by the collision of the Caribbean and Nazca
Plates under the South American plate, able to generate intermediate
earthquakes (70 – 200km). For this research, 1223 earthquakes recorded by the
seismic stations of “National Seismological Network of Colombia” were analyzed.
They all occurred in Santander, Colombia, within 6° - 8°N and 72° - 74°W
recorded from January 2012 to June 2016. Only 1138 events that met the seismic
nest parameters –like latitude, longitude and depth – were used to generate a
first 1D velocity model where only depth velocity variations were considered.
Furthermore, a quasi-3D model with both horizontal layers of different depths
and lateral variations was generated. Results show two continuous anomalies in
the study zone: the first one, a negative one interpreted as partial melt of
magma and, the second one, a positive anomaly interpreted to magma
crystallization.
Keywords: seismic tomography, Seismic Nest of
Bucaramanga, velocity model.
Forma de citar: Sepúlveda-Jaimes, F.J., y Cabrera-Zambrano, F.H. (2018). Tomografía
sísmica 3D del nido sísmico de Bucaramanga (Colombia). Boletín de Geología, 40(2),
15-33. DOI: 10.18273/revbol.v40n2-2018001.
INTRODUCCIÓN
El nido sísmico de Bucaramanga es una de las zonas de
mayor sismicidad intermedia en el mundo. Su zona focal representa una cuña de
profundidades de 120 – 220 km donde la mayor sismicidad ocurre entre los 120 –
180 km. En superficie, donde se proyecta la ubicación del nido se observa el
cambio de dirección de la cordillera oriental a través de los Sistemas de
Fallas principales de la región (Santa Marta – Bucaramanga, Suárez, Boyacá –
Belén, Gachetá – Chicamocha,
Carmen, San Vicente, Salinas) (Coral-Gómez,
1990) (FIGURA 1).
FIGURA 1
Mapa De Ubicación Del Marco
Tectónico De Colombia En Donde A. Muestra La Ubicación Del Territorio
Colombiano Y Sus Principales Características Tectónicas Debidas A La Subducción
De Las Placas De Nazca Y Caribe Bajo La Suramericana (Tomado De Zarifi Y Havskov, 2003), Las
Flechas Blancas Muestran La Dirección En La Que Las Placas De Nazca Y Caribe Se
Subducen Bajo La Placa B. Zoom Del Evidente Cambio De
Orientación Que Sufre La Cordillera Oriental Justo Dónde Comienza La Franja
Sísmica Del Departamento De Santander; Y C. Zoom De Algunos De Los Sistemas De
Falla Más Importantes De La Región Como La Falla Santa Marta – Bucaramanga Y La
Falla Suárez (Mapa Geológico De Colombia, 2015).
Este trabajo se enfoca en la técnica conocida como
tomografía sísmica, la cual tiene como objetivo principal la reconstrucción de
la distribución del campo de velocidades en el interior de la tierra a partir
de las observaciones de los tiempos de viaje de las ondas entre la fuente
sísmica (hipocentro) y los receptores localizados en diferentes posiciones
geográficas (Stein y Wysession, 2003).
NIDO SÍSMICO DE BUCARAMANGA
Un nido sísmico es una región en donde se puede
observar una concentración inusual de actividad sísmica de manera más o menos
continua: los hipocentros se encuentran concentrados a profundidades entre los
50 y 250 km, y con mayor intensidad que en las zonas adyacentes (Zarifi y Havskov, 2003).
Los nidos profundos e intermedios están relacionados con las zonas de
subducción; entre los cuáles se encuentran los nidos de Bucaramanga en
Colombia, Vrencea en Rumania, e Hindo
Kush en Afganistán (Tryggvason y Lawson, 1970).
Zarifi y Havskov (2003)
muestran una compilación de mecanismos focales en los nidos de Hindu Kush (Afganistán) y Vrancea
(Rumania) donde las soluciones son de fallas inversas, mientras que en el nido
de Bucaramanga el mecanismo focal es variable, siendo en la mayoría de casos,
mecanismos inversos. Casi todos los estudios coinciden en que la sismicidad en
Vrancea es el resultado del progreso del desprendimiento de la placa de su
parte superior (ruptura de la placa). En el nido de Bucaramanga e Hindu Kush, la fuerza motriz
responsable no se conoce exactamente. Muchos estudios sugieren que, en el nido
de Bucaramanga, el mecanismo responsable es la migración de fluidos o reacciones
de deshidratación o un complejo campo de esfuerzo concentrado, mientras que en
el nido de Hindu Kush el
mecanismo posible es la colisión de dos placas subductivas
de direcciones opuestas.
El nido sísmico de Bucaramanga está situado en la
parte nororiental del territorio colombiano, aproximadamente entre las
coordenadas 6,5° - 7° N y 73° - 73,5° W (FIGURA 1),
que corresponde a la zona con una concentración singular de actividad sísmica a
profundidades mayores a los 120 kilómetros (Perico-Martínez y Perico-Granados, 2014).
En el área se registra aproximadamente un sismo con magnitud 4,7 por mes cerca
de la ciudad de Bucaramanga (Prieto
et al.,
2012).
La actividad sísmica del nido de Bucaramanga se
concentra en un volumen mucho más pequeño (128 km3 (Schneider et al., 1987))
que los otros nidos a profundidades intermedias de la misma clase (nidos de Vrencea e Hindo Kush) (Schneider
et al.,
1987; Zarifi y Havskov, 2003) (FIGURA 2); lo cual
está asociado a la complejidad tectónica del territorio Colombiano (Zarifi y Havskov, 2003).
FIGURA 2
Modelo Esquemático De La
Interacción De Las Placas Nazca Y Caribe Que Subducen
Bajo La Placa Suramericana; A) Representa La Localización Del Nido Sísmico Sugerida
Por Van Der Hilst Y Mann (1994); B) Muestra La
Ubicación Sugerida Por Pennington (1983) Y C)
Representa El Modelo Más Actual Sugerido Por Zarifi
Et Al. (2007) (Modificado De Zarifi Et Al., 2007).
Una de las teorías que explica el origen del nido
sísmico, es la generación y migración de fluido acompañado de la fusión parcial
de la peridiotita (roca ígnea característica del
manto), causada por el escape de agua contenido en una dorsal oceánica o un
arco de isla durante la subducción (Pennington, 1981)
que a su vez es generada por el calentamiento y el cizallamiento a lo largo de
los bloques de subducción (Schneider
et al.,
1987); o también puede deberse a un complejo campo de
tensión dentro de la zona donde se ubica el nido sísmico (Zarifi y Havskov, 2003).
Debido a la complejidad del nido se han propuesto
diversos modelos tectónicos por lo que es difícil establecer a qué placa subduciente se encuentra relacionado. Pennington (1983)
sugiere que la placa del Caribe está subduciendo bajo
la placa Suramericana en dirección Sur – Este, y por ende, estaría localizado
en ésta placa. Van
der Hilst y Mann (1994)
proponen que el nido se encuentra localizado sobre la placa de Nazca, en el
segmento redefinido como el Bloque de Bucaramanga; y un tercer modelo, sugiere
que la colisión que ocurre entre ambas placas (Nazca y Caribe) es la causante
de su sismicidad (Zarifi et
al., 2007).
El registro de la sismicidad en la zona empezó a
datarse desde siglo XVI con la identificación de 111 eventos con epicentro en
el territorio santadereano hasta 1963 (Ramírez, 1969; Salcedo, 1999);
220 eventos de magnitud promedio por debajo de 4,5 entre 1963 y 1968 (Tryggvason y Lawson, 1970);
20 sismos con magnitud promedio debajo de 4,2 asociados al nido sísmico en 1976
(Pennington et
al., 1981).
Entre 1978 y 1980, con la instalación de una estación
sismológica permanente en Bucaramanga se registraron 742 sismos asociados al
nido a profundidades entre los 110 y 160 km (Gómez-Padilla, 1980).
Posteriormente en 1983, con la implementación de 12 estaciones sismológicas se
registraron 140 eventos ocurridos en un volumen de 128 km3 en la
zona (Schneider
et al.,
1987). Finalmente en 1985, las 16 estaciones sismológicas
activas registraron un total de 92 sismos ocurridos en el nido de Bucaramanga (Rivera, 1989).
En promedio la magnitud de los eventos ocurridos es
menor de 4,7; los cuales tiene origen en un volumen relativamente pequeño (128
km3 aproximadamente) y a grandes profundidades (140-180 kilómetros)
(Salcedo,
1999).
METODOLOGÍA
Uno de los problemas fundamentales en la sismología
corresponde a la localización de hipocentros, es decir, la ubicación espacio –
temporal (Xs
, Ys , Zs , ts ),
siendo Xs , Ys , Zs , las
coordenadas espaciales de la fuente (entendida esta como fuente puntual) y ts el
tiempo de ocurrencia.
Por otra parte, también está el problema de determinar
la estructura interna de la tierra, generalmente estudiada a partir de la
determinación del campo de velocidades, principalmente a partir de los primeros
arribos de la Onda P (Aki et
al., 1976).
Los dos problemas mencionados pueden ser abordados de manera conjunta (Crosson, 1976) o
de manera separada (Aki et
al., 1976).
Lo fundamental que plantean varios autores, es la posibilidad y necesidad de
desacoplar los dos problemas inversos (Kissling, 1988).
El problema simultáneo puede plantearse de la
siguiente manera: En un conjunto de estaciones en superficie E 1,
..., Ep
, (p – estaciones) se registran una serie de eventos S 1,
S
2, ..., Sq ,
(q eventos sísmicos). Para cada evento sísmico, se requiere su localización
espacio – temporal (X
1i
, X
2i
, X
3i
, X
4i
) para i=1,2…,q; siendo X → X 1,
Y
→ X
2, Z
→ X
3 y t
→ X
4; por otra parte, si la región en estudio, se divide en bloques V1 ,
V2
, ... Vl ,
donde cada bloque Vk se identifica
por el valor de velocidad de propagación de la Onda P: VK , se
requiere también el valor de la velocidad en cada uno de ellos; generados a
partir de los tiempos de tránsito desde los hipocentros i hasta las estaciones j: T
ij
, i=1, 2, …, q; j=1, 2, …, p.
En general, los tiempos de tránsito se caracterizan
por ser una función no lineal tanto de las coordenadas hipocentrales
como de los parámetros del campo de velocidad (Kissling, 1988; Crosson, 1976):
T
ij
= T ij (X 1i
, X 2i ,
..., X
1q
, X
2q
, X
3q
, X
4q
, V
1, V
2, ... Vl ,)
La expresión analítica para los tiempos de tránsito
es:
(01)
donde Hi son las coordenadas del i-ésimo hipocentro, Ej las
coordenadas de la j-ésima estaciones y V = V(X 1,
X
2, X
3) donde (X
1, X
2, X
3) es un punto típico de la curva (rayo) que une fuente – receptor (FIGURA 3).
FIGURA 3
Representación En El Plano Tridimensional
Del Rayo Que Une La Fuente I Y El Receptor J.
Si se define la “lentitud” como U(X, Y, Z) = ,
la ecuación 1 se
puede escribir como:
(02)
siendo ds el diferencial
de longitud de arco.
Los tiempos de transito dependen de forma lineal del
conjunto de parámetros U(X, Y, Z) y
de forma no lineal de los puntos de los hipocentros Hi; por lo tanto la
expresión 2 significa una dependencia no lineal del tiempo de tránsito T ij
= (X
1i
, X
2i
, ..., X
1q
, X
2q
, X
3q
, X
4q
, V
1, V
2, ... Vl ,).
Dado el carácter no lineal, el problema de hallar las coordenadas
espacio – temporal de los hipocentros (X 1i
, X 2i ,
X
3i
, X
4i
), i=1,…,
q y los parámetros V
1, V
2, ... Vl ,
tiene que ser resuelto de manera iterativa, partiendo de una estimación inicial
de las incógnitas ... Xki 0,
donde K=1,
2, 3, 4; i=1,
2,…, q y Vl 0,
donde l=1,
2, …, l, representando cada punto de la trayectoria.
Sea el tiempo estimado , el cálculo posterior de Tij podría
ser resuelto por medio de una aproximación de primer orden de la Serie de
Taylor (Aki y Richards, 1980):
(03)
o
(04)
La ecuación 4
puede escribirse de manera matricial:
(05)
con Δm = (ΔX 11,
..., ΔX 41, ΔX 12,
..., ΔX 42, ..., ΔX 4q
, ΔV 1,
..., ΔVl ).
La matriz G es la que se muestra a continuación en
forma de bloque:
(06)
donde:
(07)
k = 1, 2, …, q y:
(08)
Las submatrices S 1,
S
2, ..., Sq hacen
referencia a la parte hipocentral del problema y las submatrices M 1, M 2, ..., Mq se
refieren a la parte tomográfica del problema
simultáneo. El problema general a resolver está dado por la ecuación 5. En general, este
problema es sobredeterminado. La primera opción para
obtener una solución, es la utilización de los mínimos cuadrados, es decir,
minimizar (GΔm –
Δt)2 o lo que es
equivalente, minimizar (GΔm – Δt)
t
. (GΔm – Δt)
lo que lleva a las ecuaciones normales:
(09)
cuya solución sería:
(10)
Algunos autores, como Aster et al. (2013) y Crosson (1976)
llaman a la matriz (GTG)-1
GT
una inversa generalizada. Una característica importante desde el
punto de vista computacional, y cuando hay un gran número de eventos sísmicos, es
que la matriz GTG
es una matriz dispersa, lo que implica ventajas en cuanto al almacenamiento de
datos y también en cuanto a los cálculos a realizar. Sin embargo, para la
solución de mínimos cuadrados existe el inconveniente de que la matriz GTG
puede ser singular o también mal condicionada (Strang, 2009).
Estos inconvenientes pueden ser entendidos a partir de la llamada
descomposición en valores singulares de la matriz G.
G = U⋀VT
donde la matriz ⋀ tiene la forma:
por lo tanto:
Los valores ⋀1, ⋀2, ..., ⋀p son los valores propios de la
matriz GTG.
Cuando alguno de los valores de ⋀ i
es muy pequeño, la matriz G será mal condicionada
(Strang, 2009),
lo que significa que la solución:
(11)
tendrá problemas de estabilidad, esto significa que para pequeños errores en los Δt,
se podría obtener errores muy grandes en la soluciones Δm.
Como alternativa al problema mencionado anteriormente,
está el llamado Algoritmo de Levenberg y Marquardt, el cual modifica la solución de mínimos
cuadrados en la siguiente forma:
Minimizar {(GΔm – Δt)
T
. (GΔm – Δt) +
β 2 Δm T
. Δm}.
Minimizar ésta expresión es equivalente a resolver el
siguiente sistema de ecuaciones normales (Aki y Richards, 1980; Aster et al., 2013):
(12)
Donde β es un
coeficiente que se ajusta a los requerimientos del problema, este parámetro es
llamado por algunos autores como parámetro de regularización de Tikhonov (Aster et al., 2013). La solución de la ecuación 12 es:
(13)
Este proceso, garantiza la estabilidad del sistema de
ecuaciones a resolver; y es también llamado Método de regularización de Tikhonov.
El control de calidad de los resultados de la
inversión, puede hacerse a través de la llamada Matriz de Resolución, que para
el caso del algoritmo Levenberg – Marquardt
viene dada por la siguiente expresión (Aki y Richards, 1980; Crosson, 1976): R = HG
donde H es la matriz:
H =
(GTG
+ β 2 I)–1 GT
de manera que,
(14)
esta ecuación indica que, si la matriz R se
diferencia de una matriz identidad I, entonces habrá una discrepancia entre lo
calculado y lo real; por el contrario, si la matriz R tiende a la matriz
identidad I
(caso ideal), esto significaría que la inversión es exacta.
Para este trabajo, se utilizaron los softwares Velest
(Kissling et
al., 1994) y
Simulps14
(Thurber, 1981; Thurber, 1992),
los cuáles contienen los algoritmos que resuelven estas ecuaciones, de manera
directa y transparente para el usuario de estos softwares. El modelo de
velocidades para la Onda P en el interior del nido sísmico de Bucaramanga se
obtuvo mediante el uso de la tomografía de velocidades calculada a partir de
los tiempos del primer arribo a las estaciones (FIGURA 4, TABLA 1)
de la Red Sismológica Nacional de Colombia (RSNC) de sismos ocurridos desde
Enero de 2012 a Junio de 2016, inicialmente localizados entre 6° - 8°N y 72° -
74°W (zona nororiental del territorio colombiano) a profundidades entre los 0 y
200 kilómetros; posteriormente fueron solamente utilizados los sismos ocurridos
entre 6,5° - 7°N y 73° - 73,5°W (zona en la que se encuentra localizado el nido
de Bucaramanga).
FIGURA 4
Mapa De Las Estaciones
Sismológicas Cercanas A La Zona De Estudio, Establecidas Por La Rsnc.
TABLA 1
Características De Las Estaciones
Sismológicas Activas Cercanas A La Zona De Estudio, De Las Cuáles Se Tomaron
Los Registros De Los Primeros Tiempos De La Onda P De Los Sismos Ocurridos Y
Tomados En Cuenta En Éste Trabajo.
Para generar el modelo de velocidades 1D se utilizó
como modelo inicial el planteado por Londoño
et al.
(2010) (TABLA 2),
en un estudio realizado cerca de la zona estudiada en éste trabajo, utilizando
los tiempos de llegada de las Ondas P generadas por los sismos ocurridos, en
éste estudio, se plantea la profundidad de Moho a los 37 kilómetros
aproximadamente (Londoño
et al.,
2010).
TABLA 2
Valores Aproximados Del Modelo De
Velocidad De La Onda P, Propuesto Por Londoño Et Al. (2010).
El modelo de velocidades en 3D se generó a partir del
modelo 1D obtenido (ver TABLA 6), sin embargo, se realizó un ajuste de éste,
dado por la ecuación 15,
para obtener las velocidades entre los 50 y 200 km de profundidad, debido a que
en el modelo 1D se obtuvo un valor constante para la velocidad en el rango de
la profundidad mencionada. Este ajuste se hizo mediante la técnica de mínimos
cuadrados utilizando los datos de la TABLA 6.
(15)
RESULTADOS
Modelo de velocidades 1D
Para la inversión 1D, se utilizó el programa Velest
(Kissling et
al., 1994)
ejecutado bajo Seisan (Havskov y Ottemoller, 1999),
el cuál calcula los tiempos de llegada mediante el trazado de rayos para el
campo de velocidades propuesto; una vez conocidos estos tiempos se determina el
RMS (raíz media cuadrada de la diferencia entre los tiempos calculados y los
tiempos observados). Inicialmente, se utilizó la subrutina Select
del software Velest, para filtrar los eventos
ocurridos según una serie de parámetros establecidos previamente (TABLA 3) y se
seleccionaron 1138 sismos que cumplieron con ellos, los cuáles se utilizaron
para generar el nuevo modelo de velocidades para la Onda P en una sola
dimensión.
TABLA 3
Parámetros Con Los Cuáles Se
Escogieron Los 1138 Sismos Que Se Utilizaron Para La Inversión 1d.
*Se estableció un
máximo de 999 estaciones como referencia, para que se tomaran todas las
estaciones que registraran cada uno de los eventos.
Se tomaron en cuenta diez parámetros de
amortiguamiento (β)
(ver sección anterior) para el proceso de inversión (TABLA 4), y se
generaron sesenta modelos de velocidad divididos en seis grupos, donde el
modelo inicial a partir del segundo grupo, es el modelo de velocidad con menor
RMS obtenido en el anterior (TABLA 5).
Finalmente, se escogió el modelo con β=1,0 del
primer grupo, debido a que fue el que presentó el menor RMS general (FIGURA 5); para éste
modelo se asume una homogeneidad lateral de manera que solo se toma en cuenta
la variación de la velocidad en profundidad (Crosson, 1976),
y así el modelo obtenido consta de cuatro capas con variación de la velocidad
en cada una de ellas (TABLA 6).
TABLA 4
Parámetros De Amortiguamiento
Para Cada Una De Las Inversiones 1d Realizadas.
TABLA 5
Rms
Obtenidos Para Cada Uno De Los Modelos De Velocidad Obtenidos En Cada Grupo.
Los Números Señalados, Son Los Valores De Los Rms Más
Bajos Para Cada Grupo De Modelos Obtenidos En Cada Una De Las Inversiones.
FIGURA 5
Gráfica Del Rms
Obtenido Para Cada Parámetro De Amortiguamiento Y Grupo De Modelos Durante La
Inversión.
TABLA 6
Valores Del Modelo De Velocidades
1d Para La Onda P En La Zona Del Nido Sísmico De Bucaramanga, Generado En Éste
Trabajo.
En la FIGURA 6 se
hace una comparación entre el modelo inicial 1D utilizado (Londoño et al., 2010) y
el obtenido en este trabajo, en donde se pueden observar los cambios en las
velocidades de cada una de las capas; en la capa más superficial hay una
variación aproximada de 0,78 km/s con respecto al modelo inicial, en la capa
comprendida entre los 10 y 50 km hay una variación aproximada de 0,14 km/s,
entre los 50 y 200 km, la variación es aproximadamente de 0,29 km/s y a partir
de los 200 km de profundidad, es de 0,28 km/s.
FIGURA 6
Comparación De Los Modelos De
Velocidad; Las Líneas Punteadas Indican El Modelo Inicial Utilizado, La Línea
Continua Indica El Modelo Generado En Éste Trabajo.
El siguiente paso, es obtener el modelo tridimensional
a partir del modelo 1D obtenido en esta sección.
Modelo cuasi 3D
La inversión 3D se realizó utilizando el programa Simulps14
(Thurber, 1981; Thurber, 1992),
el cual genera un modelo de velocidades por cada uno de los sismos utilizados.
El modelo final es el promedio ponderado de todos ellos; para esto, se le da un
peso a cada uno, dependiendo de la cercanía del hipocentro con respecto a la
zona de estudio (nido sísmico de Bucaramanga, para este trabajo), es decir, a
los eventos ocurridos dentro del nido sísmico se les dio un mayor peso que a
aquellos ocurridos en los límites, o en la corteza. Inicialmente, se definió
una malla cuyo volumen es de 8694×103 km3, con divisiones
X-Y(latitud-longitud) de 23×21 con un espaciamiento de 10 y12 km
respectivamente; se establecieron, además, un total de 12 capas a diferentes
profundidades (TABLA 7),
de manera que se trabajó con un total de 6762 bloques. Durante las primeras
inversiones, se descartaron 794 sismos debido a errores en la localización o a
su cercanía a los bordes de la malla; utilizándose al final un total de 344
eventos de los 1138 iniciales, para la determinación del modelo 3D; en el
trazado de rayos de los eventos utilizados se observa una gran densidad de
datos en la zona de estudio (FIGURA 7).
FIGURA 7
Trazado De Rayos De Los 344
Eventos Utilizados Y Registrados En Cada Una De Las Estaciones Sismológicas. A.
Proyección En Planta Donde Se Aprecia La Densidad De Rayos En La Zona Del Nido
Sísmico; B. Vista Con Respecto A La Latitud, Donde La Densidad De Los Rayos Se
Acumula Entre Los 6,5° Y 7°N; Y C. Vista Con Respecto A La Longitud Donde Los
Rayos Tienen Una Mayor Densidad Entre Los 73° Y 73,5°W
.
TABLA 7
Diferentes Profundidades Y
Valores De La Velocidad, Utilizadas En El Mallado Y En La Generación Del Modelo
De Velocidades En 3d.
En la TABLA 7
también se establecen los valores de la velocidad (en km/s) para cada una de
las profundidades, obtenidas mediante el ajuste al modelo de velocidades 1D
obtenido utilizando la ecuación 15.
Al igual que con el modelo 1D, se utilizaron
diferentes parámetros de amortiguamiento (TABLA 8),
y se determinó el punto donde la varianza de los datos y la varianza del modelo
fuera compensable (Vargas
et al.,
2003); se eligió el parámetro de amortiguamiento 350 (FIGURA 8). Teniendo en
cuenta los valores del RMS en el modelo propuesto por Londoño et al. (2010) y
los obtenidos en este trabajo, se hizo un mejoramiento casi del 15,17%.
TABLA 8
Parámetros De Amortiguamiento
Tomados En Cuenta Para La Inversión En 3d.
FIGURA 8
Curva De Varianza Del Modelo Con
Respecto A La Varianza De Los Datos, Con La Que Se Determinó Que El Parámetro
De Amortiguamiento (Λ)
Más Óptimo, Fue 350.
Los resultados de la inversión 3D se muestran en la FIGURA 9 en función de
las variaciones con respecto a la velocidad media de cada una de las capas del
modelo 1D, como de los valores del campo de velocidades en cada uno de los
bloques del modelo en 3D (Aki y Lee, 1976). En este trabajo, se obtuvieron
14 secciones a diferentes profundidades definidas en la TABLA 7, en donde se
tiene una vista en planta del campo de velocidad en cada una de ellas.
FIGURA 9
Valores Del Campo De Velocidades
De Onda P (Km/S) En Cada Uno De Los Bloques Del Modelo 3d Representado En 12
Secciones A Diferente Profundidad, A. 0 Kilómetros, B. 10 Kilómetros, C. 20
Kilómetros, D. 30 Kilómetros, E. 40 Kilómetros, F. 50 Kilómetros, G. 60
Kilómetros, H. 70 Kilómetros, I. 80 Kilómetros, J. 100 Kilómetros, K. 120
Kilómetros, Y L. 150 Kilómetros; Las Variaciones Laterales De La Velocidad Con
Respecto A La Promedio En Cada Una De Las Profundidades Se Denominan “Anomalías
Positivas” Si Son Mayores Que Ésta O “Anomalías Negativas” Si Son Menores.
Dado que el nido de Bucaramanga se encuentra a
profundidad mayor de 80 km, en esta sección de resultados solo se muestra la
matriz de resolución para profundidades entre 80 y 150 km. En la FIGURA 10 se representa
espacialmente la distribución de los elementos diagonales de la matriz de
resolución (EDR) para tener una idea de cuál es la incertidumbre en los valores
obtenidos para el modelo de velocidades en la FIGURA 9. Vargas
et al.
(2003) sugieren que las zonas del modelo obtenido que tienen
una mejor resolución, son aquellas representadas en tonos claros en la
distribución espacial EDR; esto se debe a que tanto la dirección y densidad de
los datos, como el peso de las fases de Onda P ofrecen buenas perspectivas de
resolución.
FIGURA 10
Elementos Diagonales De La Matriz
De Resolución De Los Parámetros Para La Inversión Del Modelo De Velocidad De La
Onda P En 3d, Se Tienen Planos A Diferentes Profundidades, A. 80 Kilómetros, B.
100 Kilómetros, C. 120 Kilómetros Y D. 150 Kilómetros
.
Para validar el proceso de inversión y comprobar la
confiabilidad del modelo final de velocidad de este trabajo, se realizó una
prueba estadística que consistió en perturbar el modelo obtenido (dado que no
se conoce el modelo real) y se realizó la inversión tomando como referencia de
modelo inicial el modelo perturbado y se obtuvieron nuevos modelos, se observó
que en todos los casos la convergencia es hacia el modelo de referencia. Cabe
aclarar que se tuvo en cuenta la zona a partir de los 80 km de profundidad por
ser la más cercana al nido de Bucaramanga. Los resultados se muestran en la FIGURA 11. Dado que
este proceso se realizó mediante el programa Simulps14, solo fue
posible realizar una perturbación no mayor al 5%.
FIGURA 11
Se Muestra El Contraste Entre Las
Velocidades Obtenidas A Profundidades Entre Los 80 Y 150 Km Del Modelo Final De
Este Trabajo, El Modelo Perturbado Y El Generado Durante El Segundo Proceso De
Inversión. A. 80 Km De Profundidad Y La Variación Es Del 0.2%, En B. 100 Km Y
La Variación Es Del 1,04%, C. 120 Km Del 0,63% Y D. 150 Km Con Variación Del
0,24% Aproximadamente.
La distribución espacial de EDR tiene lugar en las
zonas en donde se tiene buena cobertura de rayos sísmicos (FIGURA 10), y se
observan los cambios de velocidad que ocurren con respecto a la velocidad
promedio del medio obtenida mediante la inversión 1D.
Las capas 80 y 150 kilómetros ofrecen valores de EDR
medios (0,4 – 0,65) hacia las zonas en dónde las anomalías obtenidas en el
modelo del campo de velocidad (FIGURA 10)
empiezan a juntarse hacia el centro de la región estudiada, y valores un poco
más bajos, a sus alrededores.
Los valores de EDR más altos se obtuvieron en las
capas de 100 y 120 kilómetros (0,9 – 1); esto se debe a que la mayoría de los
sismos utilizados tuvieron ocurrencias a profundidades entre los 100 y 150 kilómetros,
de manera que se tuvo una mayor densidad de datos, la cual se evidencia en la FIGURA 7 del trazado
de rayos.
Para validar el proceso de inversión y comprobar la
confiabilidad del modelo final de velocidad de este trabajo, se realizó una
prueba estadística que consistió en perturbar el modelo obtenido (dado que no
se conoce el modelo real) y se realizó la inversión tomando como referencia de
modelo inicial, el modelo perturbado, y se obtuvieron nuevos modelos; se
observó que en todos los casos la convergencia es hacia el modelo de
referencia. Cabe aclarar que se tuvo en cuenta la zona a partir de los 80 km de
profundidad por ser la más cercana al nido de Bucaramanga. Los resultados se
muestran en la FIGURA 11.
Dado que este proceso se realizó mediante el programa Simulps14, solo fue
posible realizar una perturbación no mayor al 5%.
DISCUSIÓN DE RESULTADOS
Modelo 1D
Las variaciones en los valores de la velocidad entre
el modelo obtenido en este trabajo y el propuesto por Londoño et al. (2010) se
deben en gran medida, a que los modelos no pertenecen a la misma zona, sino a
regiones cercanas.
Teniendo en cuenta que el modelo de velocidad inicial
utilizado en este trabajo presenta cuatro capas a profundidad, se obtuvieron
también cuatro capas en el modelo de velocidades 1D en el cual la velocidad va
aumentando conforme aumenta la profundidad. Las velocidades obtenidas en las
dos primeras capas del modelo (0 a 10 km, y 10 a 50 km), corresponden a
velocidades de la granodiorita, roca granítica composicional general de la
corteza y al rango de velocidades que se puede encontrar en la corteza
terrestre; las velocidades obtenidas a partir de los 50 km están relacionadas con
las de la peridiotita, que es la roca ígnea
característica del manto terrestre.
Con el modelo 1D generado en este trabajo se obtuvo un
RMS (Raíz media cuadrada) de los tiempos residuales de 0,089232 s, que representa
una reducción casi del 28% con respecto al RMS del modelo inicial obtenido por Londoño et al. (2010)
que fue de 0,32.
Modelo 3D
Teniendo en cuenta el modelo de velocidades 1D
generado y el ajuste que se le realizó para la inversión 3D, en éste trabajo se
encontró que la Onda P en la corteza posee velocidades desde los 6,52 km/s
hasta poco más de 6,85 km/s. Entre los 40 y 50 km de profundidad se presenta
una variación composicional evidenciada en el aumento de la velocidad, además
de un contraste en el tamaño de los cuerpos anómalos encontrados en el área de
trabajo; debido a esto se sugiere que la discontinuidad de Moho se encuentra en
esta zona. No es posible determinar la profundidad específica entre el rango de
40-50 km debido a que la resolución para éste trabajo es mayor de los 10 km.
De manera general, el modelo de velocidades en 3D
determinó dos cuerpos anómalos en el área de estudio: a) una anomalía positiva
(mayor velocidad que el promedio) y b) una anomalía negativa (velocidad
inferior) que se hacen más evidentes entre los 60 y 150 kilómetros de
profundidad. El primer cuerpo anómalo es interpretado como una ascensión y
cristalización magmática, mientras el segundo está relacionado con la fusión parcial
del magma por el cambio de composición. Las anomalías presentes en la zona del
nido obtenidas en este trabajo, se pueden entender como dos cuerpos en dónde su
composición comienza a variar poco antes de los 50 kilómetros de profundidad,
cambiando la velocidad de propagación de Onda P en ellos. En el manto la
velocidad promedio aumenta paulatinamente hasta alcanzar los 8,17 km/s a una
profundidad de 150 kilómetros.
El magma empieza a ascender desde profundidades
mayores a los 150 kilómetros, donde se presenta un acuñamiento
total de los cuerpos magmáticos anómalos encontrados en el área del estudio,
que se asocia a la cuña del manto siendo la zona que se localiza sobre una o
más placas que se subducen bajo la corteza terrestre
(Van
Keken, 2003), es decir
que dicho acuñamiento está posiblemente localizado
sobre el punto donde convergen las Placas de Nazca y Caribe bajo la Placa
Suramericana, siendo esta convergencia la responsable de la fusión parcial
ocurrida a esa profundidad que facilita el ascenso del magma.
El acuñamiento que se
empieza a presentar a partir de los 60 kilómetros de profundidad y que a su vez, se hace más evidente a partir de los 100
kilómetros, puede estar influenciado por la subducción de las Placas
tectónicas, delimitando los cuerpos magmáticos.
La presencia de cuerpos magmáticos de menor velocidad
al valor promedio en el manto, se deben a procesos de fusión parcial donde se
originan magmas basálticos y quedan como residuos harzburgitas
(80% olivino, 20% ortopiroxenos) y dunitas (100%
olivino); mientras que los cuerpos con mayor velocidad que el valor promedio se
deben a procesos de cristalización del magma a medida que este asciende.
El proceso de fusión parcial se hace evidente a los
100 kilómetros de profundidad en donde se presentan un cuerpo magmático en el
que la Onda P se propaga a una velocidad aproximada de 7,59 km/s (menor a la
velocidad promedio) que puede asociarse a un cuerpo de dunita (Anderson et al., 1968)
como residuo de la fusión parcial, que se encuentra rodeado de un cuerpo con
mayor velocidad asociado a la harzburgita. Por otra
parte, a la misma profundidad se encuentra un cuerpo magmático en donde la
velocidad de propagación de la Onda P es aproximadamente 7,64 km/s (mayor que
la velocidad promedio) el cual se puede asociar a la presencia de piroxeno en altas proporciones en piroxenita
y/o a la cristalización del magma y mineralización del piroxeno
(Birch, 1960).
En la corteza se tienen tres capas con diferente
velocidad promedio cada una, la primera va desde la superficie hasta los 10
kilómetros con una velocidad de 6,08 km/s; de los 10 a los 30 kilómetros se
presenta una velocidad de 6,66 km/s y finalmente de los 30 a poco antes de los
40 kilómetros con una velocidad de 6,8476 km/s. En superficie, se tiene un
medio aparentemente homogéneo dónde la velocidad promedio determinada es de
6,08 km/s; sin embargo, se encuentra una pequeña anomalía negativa (menor
velocidad) en la zona de estudio que puede ser el resultado de un cambio en el
gradiente local de temperatura, debido a que parte del fluido de la fusión
parcial ocurrida a profundidad asciende hasta llegar a superficie. En
superficie la anomalía negativa empieza a homogeneizarse con el medio general,
esto se evidencia en que la diferencia entre las velocidades es mínima.
CONCLUSIONES
·
El modelo de velocidades 1D que se generó en éste
trabajo tuvo una reducción casi del 28% en el RMS con respecto al modelo
inicial utilizado lo cual indica que se obtuvo una mayor aproximación al modelo
verdadero.
·
Se determinaron dos anomalías en el área de estudio:
una anomalía negativa que es consecuencia de la fusión parcial ocurrida en la
corteza y el manto terrestre asociada a un alto contenido de olivino; la
segunda anomalía es generada por la cristalización del magma y está asociada al
alto contenido de piroxenos.
·
La alta resolución del modelo obtenido está a partir
de los 100 y 120 km de profundidad debido a la gran densidad de rayos y datos
en ésta zona.
·
La sismicidad en el nido, también puede estar
relacionada al acuñamiento que se empieza a presentar
aproximadamente a los 120 km de profundidad, por influencia de la subducción
ocurrida, y la fusión parcial ocurrida por los fluidos acuosos provenientes de
parte de la placa subducida que recibe la cuña del
manto a los 150 km de profundidad; en dónde el adelgazamiento de los cuerpos
magmáticos anómalos se hace más evidente.
AGRADECIMIENTOS
A la Red Sismológica Nacional de Colombia (RSNC) por proporcionarnos los
datos de cada uno de los sismos utilizados en éste trabajo.
Al proyecto de investigación “Migración sísmica preapilado en profundidad por extrapolación de campos de
onda utilizando computación de alto desempeño para datos masivos en zonas
complejas” co-financiado por COLCIENCIAS-ECOPETROL,
por el apoyo para el desarrollo del trabajo.
Al Profesor Elkin de Jesús Salcedo Hurtado y todo el
equipo técnico del Observatorio Sismológico del Suroccidente Colombiano (OSSO)
por proporcionarnos los conocimientos necesarios para el manejo de los
softwares utilizados.
REFERENCIAS
Aki,
K., Christoffersson, A., and Husebye,
E. (1976). Determination of the three-dimensional seismic structure of the
lithosphere. Journal
of Geophysical Research: Solid Earth and Planets, 82(2),
277-296.
Aki,
K., and Lee, H.K. (1976). Determination of three-dimensional velocity anomalies
under a seismic array using first P arrival times from local earthquakes: 1. A
homogeneous initial model. Journal
of Geophysical Research: Solid Earth and Planets, 81(23),
4381-4399. doi: 10.1029/JB081i023p04381.
Aki,
K., and Richards, P.G. (1980). Quantitative seismology, theory and methods.
Vol. 1. New York: W.H.Freeman
& Co Ltd.
Anderson,
D.L., Schreiber, E., Liebermann, C., and Soga, N. (1968). Some elastic constant
data on minerals relevant to geophysics. Reviews of Geophysics, 6(4),
491-524. doi: 10.1029/RG006i004p00491.
Aster,
R., Borchers, B., and Thurber, C. (2013). Parameter
estimation and inverse problems. 2nd Edition. Amsterdám:
ElSevier Academic Press.
Birch,
F. (1960). The velocity of compressional waves in rocks to 10 kilobars: 1. Journal of Geophysical Research, 65(4), 1083-1102.
Coral-Gómez, C.E. (1990). La
convergencia de placas en el Noroccidente Suramericano y el origen del nido
sísmico de Bucaramanga. Revista
Académica Colombiana de Ciencias Exactas, Físicas y Naturales, 17(66),
521-529.
Crosson,
R.S. (1976). Crustal
structure modeling of earthquake data, I, Simultaneous least squares estimation
of hypocenter and velocity parameters. Journal Geophysical Research, 81(17),
3036-3046. doi: 10.1029/JB081i017p03036.
Gómez-Padilla, J. (1980).
Actividad sísmica en el departamento de Santander. Boletín de Geología,
14(28),
3-23.
Havskov, J., and Ottemoller,
L. (1999). SeisAn earthquake analysis software. Seismological
Research Letter, 70(5), 532-534. doi:
10.1785/gssrl.70.5.532.
Kissling, E. (1988). Geotomography
with local earthquake data. Review
of Geophysics, 26(4),
659-698. doi: 10.1029/RG026i004p00659.
Kissling, E., Ellsworth, W.L., Eberhart-Phillips, D., and Kradolfer,
U. (1994). Initial reference model in local earthquake tomography. Journal of
Geophysical Research: Solid Earth, 99(B10), 19635-19646. doi: 10.1029/93JB03138.
Londoño, J.M.,
Bohorquez, O.P., and Ospina, L.F. (2010). Tomografía
sísmica 3D del sector de Cúcuta, Colombia. Boletín de Geología,
32(1),
107-124.
Pennington,
W.D. (1981). Subduction of the eastern Panama basin and sismo
tectonics of northwestern South America. Journal of Geophysical Research,
86(B11),
10753-10770.
Pennington,
W.D., Mooney, W.D., Van Hissenhovenet, R., Meyer, H.,
Ramírez, J.E., and Meyer, R.P. (1981). Resultado
de un estudio de microsismos en Bucaramanga, Colombia. En: Investigaciones geofísicas sobre las
estructuras océanocontinentales del Occidente
Colombiano (pp. 42-64). Proyecto
Nariño II y III. Bogotá.
Pennington,
W.D. (1983). Role of shallow phase changes in the subduction of oceanic crust. Science, 220(4601),
1045-1047. doi: 10.1126/science.220.4601.1045.
Perico-Martínez,
N.R., y Perico-Granados, N.R. (2014). Caracterización y recurrencia sísmica del
nido de Bucaramanga. V
Congreso Internacional de Ingeniería Civil. Universidad Santo
Tomás, Tunja, Colombia.
Prieto, G.A., Beroza, G.C., Barrett, S.A,
López, G.A., and Flórez, M. (2012). Earthquake
nests as natural laboratories for the study of intermediate-depth earthquake
mechanics. Tectonophysics, 570-571,
42-56.
Ramírez, J.E.
(1969). Historia
de los terremotos en Colombia. Bogotá: Instituto Geográfico Agustin Codazzi.
Rivera, A.
(1989). Inversion du tenseur
des constraintes et des mécanismes
au foyer a partir des données
de polarités pour une population de seísmes. Application a l’Ecude du foyer de
seismicite intermeiate de
Bucaramanga – Colombia. Tésis de Doctorado. Universite Louis – Pasteur de Strasburg.
266p.
Salcedo, E.J.
(1999). Estudio de sismicidad histórica en la región de Bucaramanga (Colombia).
Revista de
la Academia Colombiana de Ciencias, Exactas, Físicas y Naturales,
23(87),
233-248.
Schneider,
J.F., Pennington, W.D., and Meyer, R.P. (1987). Microseismicity
and focal mechanisms of the intermediate-depth Bucaramanga Nest, Colombia. Journal of
Geophysical Research, 92(B13), 13913-13926.
Stein,
S., and Wysession, M. (2003). An introduction to seismology,
earthquakes, and earth structure. MA, USA: Blackwell Publishing.
Strang, G. (2009). Introduction to linear algebra.
4th Edition. MA, USA: Wellesley Cambridge Press.
Thurber,
C.H. (1981). Earth structure and earthquake locations in the Coyote Lake area,
central California. Ph.D. Thesis, Massachusetts Institute of Technology, USA.
Thurber,
C.H. (1992). Hypocenter-velocity structure coupling in local earthquake
tomography. Physics
of the Earth and Planetary Interiors, 75(1-3), 55-62. doi: 10.1016/0031-9201(92)90117-E.
Tryggvason,
E., and Lawson, J.E., Jr. (1970). The intermediate earthquake source near
Bucaramanga, Colombia. Bullentin of the
Seismological Society of America, 60(1), 269-279.
Van
der Hilst, R., and Mann, P. (1994). Tectonic
implication of tomographic images of subducted
lithosphere beneath northwestern South America. Geology, 22,
451-454.
Van
Keken, P.E. (2003). The structure and dynamics of the
mantle wedge. Earth
and Planetary Science Letters, 215(3-4), 323-338. doi: 10.1016/S0012-821X(03)00460-6.
Vargas,
C.A., Pujades, L.G., Ugalde,
A., y Canas, J.A. (2003). Tomografía
sísmica local en el territorio colombiano. Revista Internacional de Métodos
Numéricos para Cálculo y Diseño de Ingeniería, 19(3),
255-278.
Zarifi, Z., and Havskov,
J. (2003). Characteristics of dense nests of deep and intermediate depth
seismicity. Advances
in Geophysics, 46,
237-278.
Zarifi, Z., Havskov,
J., and Hanyga, A. (2007). An insight into the
Bucaramanga nest. Tectonophysics,
443(1-2),
93-105. doi: 10.1016/j.tecto.2007.06.004.