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


Francisco Henry Cabrera-Zambrano

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).

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).
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).

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).
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).

REPRESENTACIÓN EN EL PLANO
TRIDIMENSIONAL DEL RAYO QUE UNE LA FUENTE i Y EL RECEPTOR j.
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 = UVT

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 LevenbergMarquardt 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).

MAPA
DE LAS ESTACIONES SISMOLÓGICAS CERCANAS A LA ZONA DE ESTUDIO, ESTABLECIDAS POR
LA RSNC.
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.


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).


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.


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.


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.


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.

GRÁFICA
DEL RMS OBTENIDO PARA CADA PARÁMETRO DE AMORTIGUAMIENTO Y GRUPO DE MODELOS
DURANTE LA INVERSIÓN.
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.


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.

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.
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).

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.
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.


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.


PARÁMETROS DE AMORTIGUAMIENTO TOMADOS
EN CUENTA PARA LA INVERSIÓN EN 3D.

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.
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.

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.
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.

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.
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%.

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.
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.