ESCALAMIENTO DE FUENTE SÍSMICA EN EL
NIDO DE BUCARAMANGA
GABRIEL LÓPEZ ANGARITA, PROF. GERMÁN PRIETO PHD
UNIVERSIDAD DE LOS ANDES
2011
2
Índice General
Motivación
Parte I
1. Introducción a la teoría de fuente sísmica
1.1. Ecuación de onda sísmica
1.2. Solución en términos de funciones de Green
1.3. Representación de fuente sísmica
1.4. Solución explícita de las ecuaciones elásto-dinámicas
1.5. Teoría de fuente sísmica
2. Modelos de fuente sísmica
2.1 Forma de onda y aproximación de campo lejano
2.2 Espectro de la forma de onda
2.3 Modelos específicos de fuente
2.3.1 Falla de propagación unidireccional (Haskel)
2.3.2 Falla de propagación radial (Savage)
2.3.3 Falla de propagación radial (Sato, Hirasawa)
2.3.4 Falla de propagación radial (Molnar, Tucker y Brune)
2.3.5 Falla de propagación elíptica (Dahlen)
2.3.6 Falla circular que se detiene (Madariaga)
3. Dinámica y cinemática de la fuente sísmica
3.1. Leyes de Coulomb y Amonton
3.2. Caída de estrés y velocidad de ruptura
3.3. Modelos empíricos de fuente sísmica
3.4. Caída de estrés
3.5. Escalas de magnitud y energía sísmica
3.6. Escalamiento sísmico y auto-similaridad
Parte II
1. Introducción
2. Datos
3. El nido de Bucaramanga
4. Calculo de la distribución de ruido
5. Procesamiento de datos
6. Estimación de parámetros de fuente
7. Conclusiones
8. Tabla de resultados específicos
Referencias
Agradecimientos
3
Motivación
La hipótesis de escalamiento sísmica sugiere que los terremotos pequeños son solo
“versiones miniatura” de los eventos más grandes, diferenciándose únicamente por
factor de escala. En este trabajo estudiamos la veracidad o falsedad de la hipótesis de
escalamiento de fuente sísmica en una región en el noreste de Colombia cuya
tectónica ha sido fuente de debate por parte de los investigadores por muchos años: El
nido de Bucaramanga.
Esta es una región sísmicamente muy activa que presenta terremotos a profundidades
intermedias, los cuales están concentrados en un área de menos de 13 km2. En este
estudio aplicamos herramientas de análisis espectral, estudio de ruido sísmico, y
modelamiento para determinar las propiedades de la fuente sísmica que dan origen a
las señales detectadas en la superficie.
Los parámetros de fuente se extraen directamente de los datos a partir de la
comparación directa de las formas espectrales resultantes con respecto a un modelo
teórico de fuente. Realizamos un análisis independiente tanto de las ondas P
(longitudinales) como de las ondas S (transversales) presentes en un terremoto, lo cual
nos permite estudiar parámetros de mismo, tales como la partición de la energía
sísmica, la caída de estrés aparente y la dimensión de la fuente, entre otras.
El estudio conjunto de estos parámetros de fuente nos lleva a la falsificación de la
hipótesis de escalamiento para los terremotos en el nido de Bucaramanga.
4
Parte I
En la primera parte de este documento describiremos la teoría necesaria para el
desarrollo y completa comprensión de los temas incluidos en este trabajo. No
obstante, debido a la cantidad limitada de espacio, realizaremos un paseo rápido por
las herramientas teóricas más importantes enfatizando únicamente los resultados más
relevantes para el desarrollo del documento.
En el capítulo 1 comenzaremos por estudiar las ecuaciones que describen el
movimiento de la tierra (ecuaciones elastodinámicas), derivaremos soluciones
analíticas de las mismas para el caso de un medio de propagación isotrópico y
homogéneo, y posteriormente nos centraremos en el término de fuente de las
soluciones, con el fin de describirlo y caracterizarlo.
En el capítulo 2 nos dedicaremos a repasar las características más importantes de los
modelos particulares de fuente sísmica más comunes. Los cuales se encuentran
descritos ampliamente en la literatura. Nos centraremos en las propiedades que más
nos conciernen para la caracterización de la fuente, como lo son la caída del estrés, las
frecuencias de esquina, la tasa de decaimiento a altas frecuencias, entre otros
parámetros de fuente.
Finalmente en el capítulo 3 nos concentraremos en la descripción intuitiva de la fuente
sísmica, así como de la comprensión de los procesos internos que tienen lugar en una
ruptura. Definiremos muchas de las cantidades que serán relevantes para el estudio
experimental de la fuente sísmica a través de sus parámetros propios, tales como la
caída de estrés y las frecuencias de esquina. Al mismo tiempo introduciremos muchos
de los conceptos que serán la clave para la comprensión de la segunda parte de este
trabajo.
5
Capitulo 1
Introducción a la teoría de fuente sísmica
En este capítulo estudiaremos de forma breve y compacta la mecánica del continuo
aplicada sobre un medio uniforme y homogéneo, y veremos cómo las ecuaciones de
movimiento admiten soluciones ondulatorias en el campo de desplazamiento, las
cuales se pueden expresar como productos de convolución entre la fuente generadora
de la perturbación (una dislocación en el campo de desplazamiento para el caso de un
terremoto) y el propagador asociado (función de Green). En esta descripción
utilizaremos la aproximación de respuesta de sólido elástico para el medio de
propagación, esto nos permitirá plantear relaciones que permitan establecer la
cantidad de deformación ocurrida en el medio, cuando se aplica un conjunto
determinado de esfuerzos sobre su superficie. Una vez encontrada una expresión
cuantitativa para la propagación de ondas mecánicas en el medio, estudiaremos el
campo de ondas visto por un observador lejano, lo cual nos permitirá diferenciar los
campos de onda P y S, y posteriormente aislar la contribución de la fuente. Con el
modelo de fuente planteado, basta con estudiar geometrías específicas de falla que
dan lugar a diferentes modelos de fuente sísmica, la digresión se centrará
fundamentalmente en el modelo de Madariaga (1976), el cual es uno de los modelos
más utilizados en la actualidad para describir terremotos pequeños.
Debido a la gran cantidad de tema por abordar y a la cantidad limitada de espacio
disponible para este fin, me limitaré a enfocarme únicamente en las expresiones más
importantes para nuestro propósito, dejando muchos de los detalles matemáticos y
demostraciones al lector. Para una revisión completa de estos temas invito al lector a
revisar la excelente literatura disponible para este fin1.
1.1 Ecuación de onda sísmica
Considere un elemento de volumen infinitesimal en el interior de un medio de material
uniforme y homogéneo. Para este desarrollo (y en general en toda la sismología),
utilizaremos una aproximación Lagrangiana al problema, lo cual significa que
seguiremos la evolución del centro de masa de nuestro elemento infinitesimal de
volumen a medida que este se desplaza por el medio, partiendo inicialmente de una
posición de referencia en un tiempo de referencia.
1 Véase por ejemplo: Aki & Ricards (1980), Shearer (1999), Chapman (2004), Torne
(1995), Lowrie (2007), Kostrov & Das (2005), entre otros.
6
Describiremos el desplazamiento del CM de un elemento de volumen en el medio
mediante un vector en un sistema cartesiano el cual depende de
una posición de referencia en el medio y evoluciona con un tiempo . Así el
desplazamiento en términos de la función de referencia luego de un tiempo se
escribe (t). Luego la velocidad del CM viene dada por y su
aceleración por . Desde ahora utilizaremos para describir el vector de
desplazamiento en un punto del medio.
Definidas las variables cinemáticas a utilizar, pasaremos ahora a una formulación
dinámica que nos permita establecer la ecuación de movimiento de nuestra partícula.
Las fuerzas aplicadas a nuestro cubo infinitesimal pueden ser de dos tipos: fuerzas de
cuerpo, esto es, fuerzas que actúan sobre cada partícula del volumen (por ejemplo la
gravedad) las cuales denotaremos con el símbolo , y fuerzas de superficie, que
actúan únicamente en la superficie del cuerpo estudiado (por ejemplo la presión y los
esfuerzos mecánicos) las cuales denotaremos con el símbolo , con el vector
normal a la superficie.
Puede comprobarse entonces que segunda ley de Newton requiere se cumpla
Donde es la densidad del medio, es su volumen, es su superficie. , es la
tracción (fuerza por unidad de área) sobre un elemento de área infinitesimal en la
superficie con vector normal . Esta fuerza se debe a las presiones que el medio ejerce
sobre el volumen. por su lado está relacionado con la fuerza que origina el
movimiento en el medio, y la asociaremos en general con la fuente del terremoto.
Para resolver el término de tracción en la ecuación, recurrimos al tensor de estrés el
cual se relaciona con la tracción mediante2
Recurriendo entonces al teorema de la divergencia de Gauss tenemos
2 A partir de ahora introduzco el uso de la convención de suma de Einstein para la geometría diferencial.
Dado que no nos interesa la varianza/contra-varianza de los tensores, omitiré los sistemáticamente la notación con superíndices, simplificando así la interpretación de las ecuaciones sin afectar la convención de suma.
7
Con lo cual llegamos a las ecuaciones de movimiento. Vale notar ahora la simetría en el
tensor de estrés producto de la conservación del momento angular: , lo cual
equivale al hecho de que el medio no admite rotaciones netas.
Esta descripción del sistema resulta incompleta ya que no hemos especificado las
propiedades mecánicas del material, para hacerlo, introducimos la relación lineal
estrés-deformación
Donde es el tensor de deformación, el cual cuantifica la deformación de un cuerpo
a lo largo del plano . A su vez, es el tensor elástico, el cual actúa como una
constante que determina la cantidad de la cantidad de deformación a lo largo del
plano producto de una fuerza con la dirección sobre una cara del cuerpo
orientada en la dirección .
Nótense las simetrías y , producto de y
correspondientemente. Puede mostrarse, que la forma del tensor de segundo orden
más general que admite este tipo de simetrías esta dado por (Jefreys & Jefreys, 1972)
Expresión que involucra dos parámetros independientes los cuales se conocen
como módulos de Lamé, y describen las propiedades mecánicas de un sólido
homogéneo e isotrópico. también se conoce como el módulo de corte o rigidez del
material. Utilizando esta ecuación la relación estrés-deformación se escribe
Tomando entonces estos resultados, podemos obtener finalmente las ecuaciones
diferenciales que describen el comportamiento del campo de desplazamiento.
O de forma equivalente
8
Estas ecuaciones se conocen como las ecuaciones elásto-dinámicas o ecuaciones de
onda sísmica, y juegan un rol crucial en la sismología.
1.2 Solución en términos de funciones de Green
Basta notar la complejidad de las ecuaciones diferenciales (ecuaciones elásto-
dinámicas) involucradas en el movimiento del medio, para percibir inmediatamente la
necesidad de un método poderoso de solución de ecuaciones. Recurrimos al método
de Green ya que permite obtener una expresión analítica para la solución del campo
de onda en términos de su fuente (nuestro campo de interés) y su camino de
propagación. Utilizaremos un método basado en la relación de Betti, también llamada
la relación de reciprocidad, el cual nos permitirá encontrar una solución formal a las
ecuaciones.
Para plantear la relación de Betti, supongamos un volumen limitado por una
frontera . Definamos un campo de onda producto de una fuerza , junto con
las correspondientes condiciones iníciales en y de frontera en . Tome también
un segundo campo de desplazamiento producto de una fuerza con
condiciones iníciales y de frontera en general diferentes a las de . Sea la
tracción debido al desplazamiento y la tracción devida al desplazamiento .
La relación de Betti es entonces:
Luego de un poco de manipulación algebraica, integrando sobre un rango temporal
, y omitiendo por simplicidad la dependencia en , tenemos
Por compleja que pueda parecer esta expresión, su relevancia será muy notoria, ya
que nos permitirá solucionar para el campo de onda sísmica si introducimos una
fuente de tipo impulsivo aplicada en , en la dirección de la forma
. La solución a esta fuente impulsiva, es por definición la
función de Green , la cual satisface la ecuación
Note que la función de Green depende de las posiciones espaciales y temporales de la
fuente sísmica y del punto de observación , así como de la dirección de
9
aplicación de la fuerza inpulsorial (tipo ) y la dirección de observación (recuerde que
el campo de ondas de desplazamiento es vectorial). La función de Green satisface
también relaciones de reciprocidad espacio-temporales
Con lo cual, reemplazando en la relación deducida a partir de la relación de Betti
tenemos el teorema de representación
1.3 Representación de fuente sísmica
Concluimos la sección anterior con una ecuación que mediante funciones de Green
puede darnos una expresión teórica para el campo de desplazamiento observado en
cualquier posición del medio producto de una fuente ubicada en .
Desafortunadamente esta expresión no parece contener información directa del tipo
de fuente que genera el campo de onda. Intuitivamente, esperamos que podamos
describir la falla sísmica como una distribución superficial de fuerzas a lo largo de un
plano de rozamiento.
Para una descripción adecuada de la fuente sísmica imaginemos que introducimos una
“falla”, como una superficie de lados y inmersa en el medio. Esta superficie
producirá una discontinuidad en el campo de desplazamiento
3. Note que las ecuaciones elastodinámicas no aplicarán
sobre la superficie de la falla ya que el desplazamiento es discontinuo ahí, pero sí
aplicarán sobre los laterales derechos e izquierdos de la misma. El volumen de
interés equivale entonces a las superficies .
Tome ahora como límite de frontera la condición de superficie infinita , esto
elimina la contribución de a la integración. Luego a partir de la continuidad de la
tracción a lo largo de la falla y la relación de reciprocidad tenemos
3 significa que el campo se está evaluando en .
10
lo cual se puede reescribir en términos de un producto de convolución como
En esta última ecuación la integración se realiza únicamente sobre el área de la falla.
Por lo tanto conocer el desplazamiento aquí equivale a conocer el desplazamiento en
cualquier lugar del medio. Esto era de esperarse ya que la solución debe ser
única en el interior del volumen de integración si establecemos sus condiciones
iníciales y de frontera. La función definida sobre la falla se conoce como
función de deslizamiento y se usa para describir diferentes geometrías de ruptura.
El término fuente a la derecha de la última ecuación se puede escribir efectivamente
en términos de una distribución de fuerzas sobre la falla. Se puede demostrar que el
equivalente en fuerza de cuerpo para una discontinuidad de desplazamiento sobre
es
Esto se demuestra llevando el término asociado a la discontinuidad en el teorema de
representación a una forma similar a la del primer término del mismo teorema, tal que
permita reconocer la cantidad que aparece en el integrando con una distribución de
fuerzas.
En la ecuación para vale la pena resaltar la presencia del producto de
convolución en el integrando del la forma . La aparición de
derivadas del propagador en esta ecuación es físicamente equivalente a tener un par
opuesto de fuerzas actuando en . Un producto de la forma (como el que
podría producirse si el campo de ondas se debe de forma total o parcial a una fuerza
de cuerpo ) equivaldría en este contexto a una única fuerza actuando en la posición
. En ambos casos la integración sobre , distribuye estas fuerzas a lo largo de
toda la superficie de . Estos términos corresponden entonces al -ésimo componente
del campo de desplazamiento en debido a la(s) fuerza(s) en .
El término que acompaña la función de Green en la expresión puede entenderse
entonces como la “fuerza” ejercida por determinado componente. Por ejemplo, en el
caso de la falla, el término corresponde a la intensidad de la acción de la
pareja de fuerzas. Las dimisiones de este término son de momento por unidad
de área, lo que lleva al la definición del tensor densidad de momento
Utilizando la definición de para un medio homogéneo e isotrópico, y suponiendo
un deslizamiento cortante (paralelo a ) tenemos y por lo tanto
11
Lo cual equivale a representar la falla como una configuración de dos pares de fuerzas
actuando en direcciones opuestas sobre cada punto .
El desarrollo de arriba se basa en obtener el campo de ondas como la superposición de
las contribuciones de cada uno de los puntos sobre la falla. No obstante, en la mayoría
de los casos, y si la dimensión característica de la falla es mucho más pequeña que la
distancia de observación, las ondas provenientes de las diferentes regiones de la
falla llegan aproximadamente en fase. Esto equivale a ver todo como una fuente
puntual de ondas, lo cual motiva la siguiente aproximación.
Con el tensor de momento definido según
Esta última fórmula nos permite encontrar el tensor momento asociado a una
distribución arbitraria de fuerzas sobre en como
1.4 Solución explícita de las ecuaciones elásto-dinámicas
En la sección anterior obtuvimos una expresión para el campo de desplazamiento en
todo punto del espacio dada una configuración de fuerzas de fuente, no obstante la
forma explícita de la función de Green permanece sin determinar. Para encontrar
la solución explicita a las ecuaciones nos distanciaremos ligeramente del método hasta
ahora abordado. No obstante debido a la complejidad matemática que involucra tal
solución, me limitaré únicamente a establecer los pasos a seguir para dar con el
resultado. Los detalles pueden encontrarse en Aki & Richard (1980).
La solución se consigue a partir de los siguientes pasos:
i. Introducir potenciales de Hemholtz para y :
ii. A partir del tipo de fuente presente, encontrar los potenciales de Hemholtz
asociados a la fuente . Esto se realiza solucionando la ecuación de Poisson
sobre un vector con fuente
12
Los potenciales se pueden expresar posteriormente como
iii. Tomar el rotacional y la divergencia de forma independiente en las ecuaciones
elásto-dinámicas para obtener ecuaciones de onda asociadas a los potenciales
del campo de desplazamiento y .
iv. Resolver estas ecuaciones para obtener expresiones explicitas de y . Esto
puede lograrse utilizando por ejemplo una transformada de Fourier para
convertir las ecuaciones diferenciales en ecuaciones algebraicas, despejando y
calculando la transformada inversa.
v. Componer el campo de desplazamiento a partir de y , según
Luego de realizar este procedimiento para el caso de interés de doble pareja de
fuerzas (descritas con un tensor densidad de momento), tenemos entonces la solución
explicita para el campo de desplazamiento sísmico.
Con función de Green
En esta ecuación es la distancia fuente-observador y . Los
términos hacen referencia a las ondas de campo cercano que aparecen en las
soluciones, los cuales son despreciables unas cuantas longitudes de onda fuera de la
fuente ya que sus decaimientos a una distancia de la fuente son muy rápidos
( , ).
Como vemos, obtenemos 2 campos de onda diferentes los cuales corresponden ondas
P y S. El campo de onda S tiene polarización transversal y se desplazan a una velocidad
. El campo de ondas P tiene una polarización longitudinal y se desplazan a una
velocidad . Para ambos campos, la amplitud de onda y su forma viene dada por
el cual resulta proporcional a las velocidades de las partículas promediadas
sobre el plano de falla. La amplitud viene modulada en ambos casos por patrones de
radiación que establecen planos nodales en los cuales el campo de desplazamiento es
siempre cero.
13
1.5 Teoría de fuente sísmica
Obtendremos ahora las formas de onda de desplazamiento para las ondas P y S en un
medio homogéneo e isotrópico en función de la dinámica de la fuente sísmica, cuya
información se encuentra codificada dentro de la función de deslizamiento
definida en la sección 1.3. Utilizando la expresión para junto con la expresión para la
función de Green, tenemos (usando la aproximación de campo lejano)
Realizando la diferenciación e ignorando aquellos términos que atenúan más rápido
que tenemos
Suponiendo ahora que la superficie es plana y que la discontinuidad de
desplazamiento ocurre en la misma dirección en todo , tenemos
Con una función escalar conocida como función fuente o función deslizamiento en
el caso de una falla de corte. Si la estación de observación está suficientemente lejos
de la falla, entonces y son aproximadamente constantes, por lo que podemos
sacarlos de la integral. Bajo esta línea de ideas, la ecuación de desplazamiento en el
campo lejano se reduce a
Esta es nuestra expresión final para el campo de onda sísmico en función de la
dinámica de la fuente codificada en . Note que el término dentro de la integral es el
único responsable de la forma de onda observada en el campo lejano.
El espectro del campo de onda dependerá entonces de la forma particular
que asuma la función , lo cual a su vez depende del modelo de fuente utilizado.
14
Existen varios modelos de fuente los cuales se diferencian principalmente en la
geometría de la ruptura y las condiciones de frontera que se imponen a ella. Los
modelos particulares de fuente se describirán en el siguiente capítulo.
15
Capitulo 2
Modelos de fuente sísmica
En este capítulo estudiaremos los diferentes modelos de fuente sísmica que pretenden
dar una descripción cuantitativa de la dinámica presente en la fuente de ondas
sísmicas. De todas las posibles fuentes del campo de onda sísmico (explosiones,
transformaciones de fase rápidas, etc.), nos concentraremos en aquellas de origen
tectónico. Esto implica un plano de falla sometido a estrés tectónico sobre el cual se
desencadena eventualmente un desplazamiento súbito paralelo al plano de falla. Este
proceso comienza inicialmente en un punto del plano y se denomina nucleación. Dada
una falla, el punto particular en el cual se desencadenará el movimiento es bastante
difícil de predecir, y depende generalmente de la distribución de estrés sobre el plano,
así como de su relativa debilidad material. Una vez el movimiento se genera en un
punto de la falla, se propaga sobre esta hasta que finalmente su energía se disipa y el
movimiento se detiene. Durante todo este proceso, el plano de ruptura genera ondas
sísmicas que luego pueden detectarse a grandes distancias de la fuente.
En el capitulo anterior obtuvimos una expresión que determina el desplazamiento
observado en el campo lejano debido a una ruptura de área y de función escalar de
fuente , esta ultima determina la forma y evolución temporal de la discontinuidad
de desplazamiento sobre la ruptura.
2.1 Forma de onda y aproximación de campo lejano
A partir de la ecuación de desplazamiento a campo lejano, podemos separar el
término de fuente de los términos de patrón de radiación y atenuación geométrica,
con este fin definimos
Donde es la velocidad de propagación de onda ( para P y para S) y es la
forma de onda de desplazamiento en el campo lejano. Nótese que es la
distancia desde la fuente infinitesimal , centrada en ; el punto de observación
es . Para describir de forma más simple el término , nótese que
16
Donde es la distancia del punto de observación al origen y es el vector unitario apuntando al punto de observación. Para valores grandes de con respecto a las dimensiones de la falla podemos aproximar según
Esta aproximación resulta válida siempre y cuando el error cometido (del orden del término más grande omitido en la expansión) sea mucho menor a un cuarto de la longitud de onda dominante emitida. En términos de la longitud de la falla esto equivale a
Bajo esta condición la forma de onda de desplazamiento queda entonces
Note que la forma de onda depende más fuertemente de la dirección de observación (posición en la esfera focal) que de la distancia .
2.2 Espectro de la forma de onda
Tomando la transformada de Fourier de tenemos
Esto implica que podemos ver el espectro de la forma de onda (con un desplazamiento
de fase de ) como una superposición de ondas planas. Note además que bajo la
identificación , la ecuación toma la forma de una transformada doble de
Fourier
Desafortunadamente a partir de la observaciones a campo lejano, podemos obtener
únicamente la proyección de al plano de falla , y dado que es un vector
unitario tenemos . Lo cual implica que no podemos ver detalles de la
17
fuente sísmica con escalas de longitudes menores a la menor longitud de onda
observada.
Con estas herramientas teóricas, podemos pasar a describir modelos particulares de
fuente, comenzando con el modelo de propagación unidireccional.
2.3 Modelos específicos de fuente
2.3.1 Falla de propagación unidireccional (Haskell)
Suponga un plano de falla con longitud y ancho . La ruptura comenzará en uno de
los extremos de la falla y se propagará con un frente plano a través de ella a velocidad
, la dirección de propagación será . Estableciendo el sistema de coordenadas
sobre un plano de falla rectangular y tomando la función de desplazamiento
como
dentro de la falla y cero fuera de ella. Definiendo como el ángulo entre el eje de
propagación y el vector hacia el punto de observación, se obtiene el siguiente
resultado para el espectro de la forma de onda (note una transformada de Fourier en
)
con
donde la envolvente representa el efecto del tamaño finito de la falla y
genera ceros en el espectro sobre cada uno de sus nodos ubicados en
Dado que representa la evolución temporal del desplazamiento sobre la falla,
deberíamos esperar que fuese inicialmente cero (antes de la ruptura), creciese
linealmente durante un tiempo finito T (a medida que se fractura el material), y
finalmente se detuviera en un desplazamiento total D (correspondiente a un
desplazamiento neto entre las caras opuestas de la falla). Escribimos entonces:
Con lo cual la magnitud del espectro de la forma de onda queda
18
Existen algunas cosas importantes que vale la pena resaltar para este resultado:
1. La envolvente produce ceros en el espectro y actúa de forma más
fuerte cuando se observa la dirección opuesta a su propagación ( ). Esto
introduce directividad en el mismo, ya que implica un menor contenido en
altas frecuencias para . El espectro de la forma de onda decae como
a altas frecuencias4.
2. La duración T del desplazamiento neto durante el sismo (o tiempo de
ascensión) se atenúa como a altas frecuencias.
3. Los efectos combinados de y T finito hacen que el espectro decaiga
como a altas frecuencias.
4. Los parámetros de fuente: describen completamente la forma de
onda en el modelo de Haskell.
5. El espectro es plano a bajas frecuencias y proporcional al momento sísmico
6. El espectro puede pensarse como el producto de un factor de finitud
(ya que contiene la dimensión de la ruptura) y un espectro de fuente puntual.
Ambos decaen según a altas frecuencias.
2.3.2 Falla de propagación radial (Savage)
A pesar de que el modelo de fuente de Haskell presenta muchas de las
características que observamos en los estudios de fuente sísmica, resulta poco
realista ya que el frente de ruptura es plano y la falla es perfectamente rectangular.
Resulta más convincente en cambio permitir que la ruptura inicie en un punto y
luego se propague radialmente, bajo un frente de onda circular , hasta
concluir el movimiento en una superficie arbitraria , donde son
coordenadas cilíndricas en el plano de falla y es la velocidad de ruptura.
Pondremos la falla centrada en el eje y propagándose sobre los ejes
Supondremos que la ruptura se propaga desde el origen en todas las direcciones
del plano de falla .
Savage (1966) modeló la discontinuidad de desplazamiento con una distribución de
Heaviside con un valor final de . Con esto la función de desplazamiento es
4 El decaimiento a altas frecuencias es importante, ya que como veremos a continuación resulta una propiedad medible en datos espectros reales. Además proporciona una forma de caracterizar los diferentes modelos de fuente.
19
Lo cual implica que solo ocurre deslizamiento si y , en este modelo
el deslizamiento ocurre instantáneamente en el frente de onda, esto es, salta de un
valor cero antes de la llegada del frente de onda, a un valor despues.
Para la ecuación de la forma de onda tenemos
Donde hemos usado coordenadas cilíndricas para expresar la posición del punto de
observación y cilíndricas para expresar la posición del
elemento infinitesimal de área en la fuente . Continuando
con la integración se encuentra:
donde , lo cual implica que la ruptura es
subsónica . Existen algunas cosas importantes que resaltar de este
modelo:
1. La forma de onda (en el caso de aproximadamente constante) presenta
una evolución como función de rampa hasta que la
señal alcanza el perímetro de la falla
2. La discontinuidad de desplazamiento en equivale a un
desplazamiento instantáneo de las partículas, lo cual introduce velocidades
y aceleraciones de tipo Delta de Dirac en estos puntos.
3. Se puede observar en el espectro una fase de nucleación con decaimientos
a altas frecuencias proporcionales a y una fase de detención con
decaimientos a altas frecuencias proporcionales a . A altas
frecuencias, la fase de detención domina la fase de nucleación, por lo que el
espectro como un todo decae según .5 La fase de nucleación propaga
con el frente de onda la información sobre el inicio de la ruptura. La fase de
detención propaga con un frente de onda reflejado la información sobre la
detención de la ruptura en los bordes de la falla; esta información se lleva al
interior de la falla luego de la reflexión del frente de onda en los bordes de
la misma.
5 De hecho la variación de la fase de detención es propia de cualquier frente de onda de superficie suave.
20
2.3.3 Falla de propagación radial (Sato, Hirasawa)
El modelo de Savage (1966) discutido en la sección anterior resulta irrealista ya que
el deslizamiento no resulta consistente con la solución estática y ocurre de forma
instantánea en todo punto. Para sortear la primera de estas dificultades Sato y
Hirasawa (1973) propusieron la siguiente función de desplazamiento
Con . Modelo que describe una falla circular bajo estrés
cortante uniforme . Inyectando este modelo en la ecuación para la forma de
onda tenemos
Donde y . Algunas cosas para notar de este
modelo son:
1. La forma de onda inicial presenta una variación en
contraposición al modelo inicial de Savage.
2. Espectralmente, la fase de nucleación presenta una variación a altas
frecuencias, y la de detención una variación , correspondiente al
decaimiento de la forma de onda completa; lo cual demuestra la
dominancia de la fase de detención sobre la fase de nucleación.
A pesar de ser una mejora al modelo de Savage, el modelo de Sato-Hirasawa todavía
presenta dificultades, ya que el movimiento se detiene al mismo tiempo sobre todo el
plano de falla. Esto es opuesto a la realidad, en la cual, el deslizamiento se detiene
sobre la ruptura únicamente luego de que el frente de onda se refleja sobre los bordes
de la falla dando lugar al frente de detención el cual detiene el deslizamiento de la roca
a medida que se propaga hacia adentro de la ruptura. Una detención inmediata sobre
todos los puntos de la falla implicaría que estos no están conectados casualmente.
21
2.3.4 Falla de propagación radial (Molnar, Tucker y Brune)
Para superar las dificultades en los modelos de Savage y Sato-Hirasawa, Molnar,
Tucker y Brune (1973) propusieron el siguiente modelo cinemático de falla circular
donde es el radio del área circular rota y es la velocidad relativa de las partículas,
la cual es constante sobre la falla. En este modelo la ruptura sufre nucleación en el
centro, crece radialmente en todas las direcciones a velocidad hasta el radio , y
luego se contrae de nuevo al centro a la misma velocidad. Aquí el deslizamiento ocurre
sobre el frente de onda, se propaga hasta la frontera y luego se refleja hacia el
interior.
Existen algunos detalles a resaltar en este modelo
1. La forma de onda inicial presenta una variación en
contraposición al modelo inicial de Savage.
2. Espectralmente la fase de detención presenta un decaimiento a altas
frecuencias. En la forma de onda completa el decaimiento a altas frecuencias
va a depender del ángulo de observación barriendo desde hasta
(Aki & Richards, 1980).
2.3.5 Falla de propagación elíptica (Dahlen)
Para generalizar los modelos de ruptura circular existentes, Dahlen (1974) propuso un
modelo de fuente cuyo frente de onda elíptico se mantiene creciendo sin perder la
forma. El utilizó la siguiente función de deslizamiento
donde es la velocidad relativa alrededor del centro de la ruptura y , son las
velocidades de propagación de ruptura en los ejes y respectivamente. Este
modelo da lugar a una caída de estrés uniforme sobre la falla. Nótese que si , la
función de deslizamiento se reduce a la del modelo de Sato-Hirasawa para .
La amplitud espectral de este modelo en el campo lejano es
22
donde los ángulos , corresponden a la posición angular del punto de observación
en coordenadas esféricas. Nótese que a diferencia de los modelos hasta ahora
estudiados, el modelo de ruptura elíptica no presenta en su espectro ningún
parámetro que haga referencia al tamaño final de la ruptura, por lo cual el modelo de
fuente resulta invariante bajo escalamiento, es decir, auto-similar. La auto-similaridad
del modelo de Dahlen está fuertemente relacionada con la dominancia de la fase de
nucleación sobre la fase de detención, y se destruiría si este orden se invirtiese.
Algunos detalles a notar de este modelo son:
1. De nuevo la forma de onda inicial presenta una forma parabólica en función del
tiempo.
2. La fase de nucleación presenta un decaimiento espectral de tipo . La fase
de nucleación resulta dominante sobre la fase de detención, la cual a diferencia
de los modelos anteriores, no ocurre de forma súbita, sino de forma pausada, a
medida que el frente de onda se adentra en regiones de mayor fricción o
menor estrés tectónico. Esto se puede derivar directamente de la función de
deslizamiento del modelo .
2.3.6 Falla circular que se detiene (Madariaga)
Hasta ahora nos hemos limitado a describir la fuente mediante modelos cinemáticos
de la misma, esto implica que nos limitamos únicamente a determinar una forma
apropiada para la función de deslizamiento tal que describa la cinemática que
queremos reproducir, y luego sencillamente la inyectamos dentro de la ecuación para
la forma de onda, con lo cual obtenemos una solución para el desplazamiento (y su
espectro) en el campo lejano. Desafortunadamente, los modelos desarrollados hasta
ahora nos dan una descripción de cómo se propaga la ruptura en la fuente, pero no
nos describen de forma satisfactoria el porqué se propaga de esta forma. Además, por
buenos que parezcan, los modelos cinemáticos generalmente resultan ser físicamente
insatisfactorios, y conducen a consecuencias que contradicen la intuición física.
Para una descripción más apropiada de la fuente sísmica es necesario adoptar un
modelo dinámico de propagación. En un modelo dinámico, la propagación misma se
describe en términos de la condiciones iníciales de estrés sobre la falla.
Consideraremos ahora el modelo dinámico propuesto por Madariaga (1976) ya que es
uno de los más satisfactorios desde el punto de vista físico para la descripción de la
fuente, además de ser el modelo más ampliamente utilizado en la actualidad para este
tipo de estudios. El modelo describe una falla circular que luego de propagarse a
velocidad constante, comienza a detenerse hasta detenerse al entrar en una región de
mayor fuerza cohesiva. Debido a que este modelo tiene en cuenta las complejas
reflexiones que ocurren una vez alcanzado el borde de la región de velocidad
23
constante, no existe solución analítica y debemos adoptar una aproximación numérica.
Madariaga solucionó este problema mediante una aproximación de diferencias finitas
con red escalonada.
Asumiremos un modelo de fuente circular de radio . Para esto posicionaremos la
ruptura a lo largo del plano con el origen en el centro de la misma. Tomaremos
perpendicular a la falla, será el campo de desplazamiento, y el tensor de
estrés. Supondremos que la caída de estrés en la ruptura ocurre únicamente en el
componente , tomaremos en el plano , y afuera del
área fracturada sobre el plano . Las condiciones de frontera de este problema
son entonces:
En coordenadas cilíndricas, las ecuaciones a solucionar son 3 ecuaciones de
movimiento producto de la relación y 6 ecuaciones que relacionan el
estrés con la deformación en un medio homogéneo e isotrópico según
, lo cual nos da lugar a ecuaciones independientes. Como incógnitas tenemos las
3 componentes de velocidad y 6 componentes del estrés, para un total de 9
ecuaciones y 9 incógnitas, con lo cual el problema tiene solución. Definimos
Con lo cual el sistema de ecuaciones a solucionar es
24
En el modelo de Madariaga (1976) la fase de nucleación decae a altas frecuencias
como y la fase de detención decae como . La fase de detención domina el
decaimiento asintótico, por lo cual es de esperar un decaimiento a altas
frecuencias.
La importancia del modelo de Madariaga radica en que es uno de los modelos más
acertados cuando se pretende describir terremotos pequeños. Como veremos en
capítulos posteriores, el nido de Bucaramanga presenta terremotos en su mayoría de
pequeña magnitud, con una tasa de terremotos grandes muy por debajo de lo
esperado. Este será nuestro modelo preferido al momento de describir la fuente
sísmica. De hecho, como se verá en la parte II, la tasa de decaimiento es la que
mejor se adapta a nuestros datos experimentales, mostrando concordancia con el
modelo de Madariaga.
La figura 3.1 muestra algunos espectros de campo lejano para las fases P y S con
velocidad de ruptura (con la velocidad de las ondas S). Puede observarse el
efecto de la directividad de la fuente sísmica si se observa el corrimiento del borde de
la figura (o frecuencia de esquina como se definirá en los siguiente capítulos) con
respecto al ángulo entre el observador y la línea de falla.
Figura 3.1. Espectro de campo lejano para las ondas P y S en el modelo de Madariaga. Tomado
de Madariaga (1976)
25
Capitulo 3
Dinámica y cinemática de la fuente
sísmica
Dedicaremos ahora nuestros esfuerzos a una comprensión un poco más intuitiva del
formalismo desarrollado en los capítulos 2 y 3. El objetivo de este capítulo consiste en
brindar al lector las herramientas necesarias para abordar la segunda parte de esta
obra. Para entender el proceso de deslizamiento en la falla conviene dividirlo en 3
etapas: (1) formación de la grieta inicial producto de la acumulación progresiva de
estrés sobre la roca. La grieta se forma cuando el estrés alcanza el punto máximo de
fractura que puede soportar el material. (2) Crecimiento de la zona de deslizamiento.
Esto equivale a la propagación del frente de ruptura hacia el exterior. (3) Terminación
del proceso de ruptura. La cual ocurre debido a que el frente de ruptura alcanza una
zona de mayor fricción, la cual ya no puede romper dado que su energía ha sido
disipada en forma de ondas sísmicas y atenuación interna (calor), además de ser
disipada por la atenuación geométrica propia de un frente en expansión.
Comenzaremos estudiando el primer proceso mediante la teoría de la mecánica de
falla, cuyas bases recaen en las leyes Coulomb y Amonton.
3.1 Leyes de Coulomb y Amonton
Las leyes de Coulomb y Amonton se encuentran en las raíces de la mecánica de fallas y
describen de forma simple, y precisa la mecánica de fuerzas internas al momento de la
fractura inicial en un terremoto.
La teoría de Coulomb establece que el estrés cortante máximo sobre una roca
justo antes de la fractura equivale a una constante denominada cohesión de la roca
más el producto del estrés normal sobre la roca y una constante denominada
coeficiente de fricción interna, más la presión de poros (en el caso de tener una roca
porosa con un líquido en su interior).
La constante está relacionada con la fuerza de ligadura de las fuerzas interatómicas
las cuales mantienen la forma y estructura del material. La constante tiene una
interpretación similar a , esta vez cuantificando la respuesta de un estrés externo
sobre el material. Esta ecuación se denomina el criterio de fractura de Coulomb.
26
Una ecuación similar cuantifica el estrés friccional sobre una falla existente una vez
que se ha excedido el estrés máximo soportado por el material y se ha
establecido un deslizamiento entre las caras internas de la falla
Aquí se conoce como el coeficiente de fricción cinética, y está relacionado con las
protrusiones de las caras internas de la falla conocidas como asperezas. Esta última
ecuación se conoce como la segunda ley de Amonton. (La primera ley de Amonton
establece que las fuerzas fricciónales son independientes del tamaño de la superficie
de la falla).
Se ha realizado trabajo experimental extensivo en este tópico, el cual ha podido
determinar que las leyes de Amonton son aplicables a en un rango amplio de valores
de estrés. En efecto el trabajo de Byerlee (1978) demostró que
Lo cual se conoce como la ley de Byerlee. Sorprendentemente, esta ley resulta casi
completamente independiente del tipo de roca en consideración.
En un terremoto, la grieta inicial se forma al comienzo debido a la acumulación de
estrés sobre el material hasta alcanzar un estrés de falla , posteriormente el
material comienza a deslizar extendiendo la grieta para formar el plano de falla.
Durante su deslizamiento siente un estrés friccional debido a la rugosidad de las
superficies en contacto. En general, entre más grande sea la fuerza de fricción, menor
será la velocidad a la cual ocurrirá el desplazamiento; esto se conoce debilitamiento
de la velocidad.
3.2 Caída de Estrés y velocidad de ruptura
Una vez el material ha generado una grieta que se ha propagado para dar lugar a una
falla, el plano de falla establece una región de debilidad en el material, por lo que
continuará deslizándose bajo la aplicación de estrés en el futuro. Con el tiempo, la falla
acumula un historial de deslizamientos, y la cantidad de estrés acumulada necesaria
para producir un deslizamiento se hará cada vez menor, por lo que se deslizará con
menor dificultad. Este proceso se conoce como debilitamiento de la falla.
La teoría descrita hasta ahora resulta útil para comprender la estática de la fuente
sísmica, pero dice muy poco acerca de la dinámica de la fuente. Dinámicamente, la
ruptura se genera en el centro de la falla y se propaga a lo largo de esta mediante un
frente de ruptura, que no es otra cosa que un frente de onda de estrés sobre la falla.
La propagación del frente de la ruptura sobre el plano de falla puede así entenderse
27
como un proceso de acumulación/liberación de estrés, tal y como lo ilustra la figura
3.1. Inicialmente un punto en el material se encuentra sujeto a un estrés que se ha
venido acumulando con el tiempo. A medida que el frente de ruptura se aproxima el
estrés sobre el punto comienza a aumentar hasta alcanzar un punto de falla , en este
punto comienza el desplazamiento. A medida que el frente de ruptura se aleja, el
estrés sobre el punto cae a un valor , donde el desplazamiento finalmente se
detiene. Luego de esto el nivel de estrés se ajusta a un nivel un poco más arriba o
debajo de , dependiendo de si ha ocurrido debilitamiento o endurecimiento6 de la
velocidad. La caída de estrés del proceso equivale entonces a
Figura 3.1. Estrés en un punto sobre la superficie de fala a medida que llega y se va el frente de ruptura.
La caída de estrés sobre el punto es a la diferencia entre los estados inicial y final de estrés en todo el
proceso . Adaptado de Lay & Wallace (1995).
Uno de los temas todavía por resolver en la dinámica de la ruptura es la duración del
deslizamiento en cada punto sobre la falla. En la mayoría de modelos de fuente, el
deslizamiento continúa incluso luego de que pase el frente de ruptura. En estos
modelos, el deslizamiento se mantiene hasta la llegada de la información sobre la
ruptura del frente de onda. Esta información proviene de un frente de sanación
producto de la reflexión de la energía sísmica al interior del área fracturada, una vez el
frente de onda comienza a desvanecerse. El deslizamiento en cada punto se detiene a
medida que el frente de sanación se propaga hacia adentro. Por ejemplo, en el caso de
una falla circular, la energía sísmica se refleja en los extremos de la falla propagándose
de nuevo al centro. De esta forma el punto central es el primero en deslizarse y el
último en detener su deslizamiento.
Sabemos que para cada punto sobre la falla, el deslizamiento comienza luego de que
pasa la punta del frente de onda, y no se detiene hasta que llega el frente de sanación.
Por lo tanto, podemos modelar el deslizamiento de cada punto sobre la falla en
función del tiempo en términos de una función de rampa, tal y como se explicó en el
modelo de Haskell. De esta forma la derivada temporal del deslizamiento, o la tasa de
deslizamiento tendrá la forma de una función boxcar (figura 3.2).
6 El endurecimiento de la velocidad es el efecto opuesto al debilitamiento de la velocidad.
28
Figura 3.2. Deslizamiento (der) y tasa de deslizamiento (izq) en un punto sobre la falla.
Sabemos que la densidad de momento sísmico es proporcional al deslizamiento
según
donde es la rigidez, es el área de la falla y es el deslizamiento. Luego el
deslizamiento total promedio es directamente proporcional al momento sísmico ,
según
Luego podemos concluir que el momento sísmico sobre toda la falla es la suma del
deslizamiento de todas las partículas sobre la falla; así, la tasa de momento es la
convolución de la historia de la tasa de deslizamiento en un punto y la historia de la
propagación del frente de ruptura a lo largo de la falla.
3.3 Modelos empíricos de fuente sísmica
A pesar de la complejidad intrínseca de los modelos teóricos, desde el punto de vista
práctico el espectro a campo lejano de una señal sísmica se puede describir
básicamente con 3 parámetros: (1) El nivel a baja frecuencia, debe ser proporcional al
momento sísmico y permite la determinación de la magnitud de un evento sísmico. (2)
La frecuencia de esquina, la cual se puede definir como la intersección de las asíntotas
de bajas y altas frecuencias en una grafica log-log de densidad espectral contra
frecuencia. La presencia de la frecuencia de esquina es una consecuencia de las
dimensiones finitas de la fuente. (3) La potencia de la asíntota a bajas frecuencias,
también llamada tasa de decaimiento. Esta potencia nos da información sobre el tipo
de fuente sísmica que tenemos, en particular, determina la fase dominante durante el
proceso de ruptura (ya sea la fase de detención o la fase de nucleación).
En el modelo de ruptura bidireccional de Savage (muy parecido al modelo presentado
en la sección 2.4), la nucleación ocurre en el centro de una falla rectangular, y se
propaga mediante dos frentes de onda laterales planos. En este caso no obstante, la
asíntota a altas frecuencias decae como . Asumiendo (con la velocidad
29
de onda S) y , los promedios angulares de las frecuencias de esquina de las
ondas P y S respectivamente, Savage (1972) encontró:
Con y la longitud y ancho de la falla rectangular. Su razón es entonces
asumiendo un sólido de Poisson (véase Shearer, 1999).
En el modelo de Sato-Hirosawa tenemos en cambio
Donde es el radio de la fuente y . Su razón es
En el modelo dinámico de Madariaga, con tenemos
Como vemos la predicción del valor de varia de modelo a modelo
dependiendo de la geometría de la falla, por lo cual puede utilizarse para determinar
aproximadamente la forma de la misma. En la mayoría de modelos, al igual que en la
gran mayoría de observaciones tenemos . Es importante notar también que
existen muchos modelos de fuente que predicen dos frecuencias de esquina en el
espectro, una a altas y otra a medianas frecuencias. No obstante, debido a las
dificultades en la observación de tales frecuencias, comúnmente se supone una única
frecuencia de esquina, ya que el ruido o la banda limitada de frecuencia entre otros
efectos, comúnmente impiden la observación de mayores detalles en el espectro
Para describir empíricamente la forma del espectro en un plano log-log, recurrimos a
modelos empíricos de fuente, el modelo general propuesto por Abercrombie (1995)
es:
30
Donde es el espectro de la fuente, es la amplitud a bajas frecuencias
(proporcional al momento sísmico), es la frecuencia de esquina, es la tasa de
decaimiento a altas frecuencias y es una constante. La variación en generalmente
se limita a los valores y , donde la principal diferencia entre estas radica en
que los espectros con valores mayores de presentan bordes más puntiagudos.
En las figuras 3.3 y 3.4 se pueden observar algunos ejemplos de espectros con
diferentes valores de parámetros , y . Nótese como la intersección aproximada de
las asíntotas de bajas y altas frecuencias coincide aproximadamente con la frecuencia
de esquina. Observe también como la variación del parámetro hace más puntiagudo
el espectro.
Figura 3.3. (Izquierda) Espectro con , , . (Derecha) Espectro con , , .
Figura 3.4. (Izquierda) Espectro con , , . (Derecha) Espectro con , , .
Dado que los espectros observados para las ondas P y S se propagan desde la fuente
por la tierra antes de alcanzar las estaciones, el espectro observado no corresponda
exactamente con el espectro de la fuente. La razón de esta discordancia radica en las
complejidades del medio de propagación y efectos tales como:
focalización/defocalización de las ondas, directividad, atenuación cercana a la
superficie, atenuación intrínseca, respuesta del instrumento, efectos de camino y
1.00.5 5.0 10.0 50.0Frequency Hz
1.000
0.500
0.100
0.050
0.010
0.005
0.001
u f
0
1.00.5 5.0 10.0 50.0Frequency Hz
0.10
1.00
0.50
0.20
0.30
0.15
0.70
u f
0
1.00.5 5.0 10.0 50.0Frequency Hz
1.000
0.500
0.100
0.050
0.010
0.005
0.001
u f
0
1.00.5 5.0 10.0 50.0Frequency Hz
0.10
1.00
0.50
0.20
0.30
0.15
0.70
u f
0
1.00.5 5.0 10.0 50.0Frequency Hz
0.10
1.00
0.50
0.20
0.30
0.15
0.70
u f
0
1.00.5 5.0 10.0 50.0Frequency Hz
0.10
1.00
0.50
0.20
0.30
0.15
0.70
u f
0
1.00.5 5.0 10.0 50.0Frequency Hz
0.10
1.00
0.50
0.20
0.30
0.15
0.70
u f
0
1.00.5 5.0 10.0 50.0Frequency Hz
0.10
1.00
0.50
0.20
0.30
0.15
0.70
u f
0
1.00.5 5.0 10.0 50.0Frequency Hz
1.000
0.500
0.100
0.050
0.010
0.005
0.001
u f
0
1.00.5 5.0 10.0 50.0Frequency Hz
1.000
0.500
0.100
0.050
0.010
0.005
0.001
u f
0
1.00.5 5.0 10.0 50.0Frequency Hz
1.000
0.500
0.100
0.050
0.010
0.005
0.001
u f
0
1.00.5 5.0 10.0 50.0Frequency Hz
1.000
0.500
0.100
0.050
0.010
0.005
0.001
u f
0
31
efectos de sitio. En general el espectro de amplitud observado se puede describir
como
donde es el espectro de la fuente, representa la influencia de la tierra a lo
largo del camino de propagación, representa la respuesta de la corteza,
incluyendo posibles reverberaciones y representa la respuesta del instrumento.
De todas las correcciones, la atenuación intrínseca es tal vez el efecto más importante
a considerar a la hora de estudiar el espectro. Esta atenuación ocurre debido a la
estructura anaelástica de la tierra o de forma equivalente a procesos de fricción
interna durante la propagación de las ondas. Este efecto se cuantifica empíricamente
mediante el parámetro de calidad . Este parámetro establece una medida de que tan
eficiente es la tierra como medio de propagación. Lo interesante de este tipo de
atenuación es que su efecto se hace principalmente fuerte en las altas frecuencias, las
cuales son fuertemente atenuadas, mientras que las bajas frecuencias pasan casi
inadvertidas por medios de alta atenuación intrínseca. Para tener en cuenta este
efecto en el espectro proponemos el espectro observado como
Donde es el espectro fuente y es el tiempo de vuelo (el tiempo que le toma al
frente de onda de cierta fase (P o S) viajar de la fuente al recibidor).
3.3 Caída de Estrés
Uno de los parámetros empíricos más importantes que se puede determinar con datos
reales es la caída del estrés . En la sección anterior definimos la caída del estrés
sobre un punto como la diferencia entre el estrés inicial y el estrés final
. Hay que tener en cuenta no obstante que tal medida puede variar
de punto a punto sobre la falla. Definimos entonces la caída de estrés estática como
Con el área de la falla. Podemos relacionar ahora la deformación
con el estrés aplicado según la ley de Hooke
32
Donde es la dimensión característica de la falla (el radio en el caso de una falla
circular, o la longuitud en el caso rectangular) y es una constante adimensional que
depende de la geometría de la falla. En la tabla 3.1 se pueden encontrar algunos
valores característicos de para algunas geometrías particulares, también se
encuentran las relaciones complementarias en términos del momento sísmico. Note
que .
Tabla 3.1 Caída de estrés y momento sísmico para algunas geometrías de falla
Falla circular
Falla de transformación
Falla de subducción
Donde , son las dimensiones de la falla rectangular, es el radio de la falla circular
y son los parámetros de Lamé. Como vemos en la tabla 3.1 la caída del estrés varía
según , lo cual introduce gran dispersión en los datos debido al gran error en
la estimación de .
Es posible estudiar el escalamiento de la caída de estrés como función del momento si
calculamos el área de la falla, la cual se puede estimar a partir de los primeros pulsos
de llegada en P y S (función fuente-tiempo) en el sismograma.7 Muchos de los estudios
muestran que es esencialmente independiente de (Prieto, 2004, Yamada,
2007, Ide y Beronza, 2001); esta idea establece la base de la hipótesis de escalamiento,
a cual exploraremos posteriormente.
3.4 Escalas de magnitud y energía Sísmica
En general resulta importante tener una forma de cuantificar el “tamaño” de un
terremoto. Historicamente se ha utilizado la “magnitud” de un sismo como una forma
de cuantificar su tamaño, existen varias escalas de magnitud en la literatura. La
magnitud local introducida en 1930 por Charles Ritcher (también conocida como
escala Ritcher) cuantifica el tamaño de los eventos en función de la amplitud más
grande de le señal detectada según
7 La función de fuente temporal debe tener una forma trapezoidal (para una fuente de Haskell) a partir
de la cual se puede medir el tiempo de ruptura , donde es el ángulo entre la falla y la estación. Si suponemos , podemos encontrar el área de ruptura en función de . El area bajo la función de fuente temporal debe ser proporcional a .
33
Donde es la amplitud máxima medida en el sismograma del evento y es la
amplitud máxima medida en un evento de referencia a la misma distancia.
De forma similar se puede definir la magnitud de onda de cuerpo , proporcional a
la amplitud máxima detectada en las señales de las ondas de cuerpo (P y S); y la
magnitud de onda superficial , proporcional a la amplitud máxima detectada en las
señales de las ondas de ondas superficiales (Rayleigh y Love). Véase (Shearer, 1999).
Una medida más exacta del tamaño de un terremoto es por supuesto el momento
sísmico el cual ya he introducido, y su escala de magnitud asociada, la magnitud de
momento , que se define según la relación de Kanamori (1977)
con dado en Nm. En una primera aproximación tenemos .
Una forma alternativa de medir el tamaño de un terremoto, es mediante la energía
sísmica liberada. Existen varias formas para realizar esto, algunas más precisas que
otras. Esta energía puede calcularse a partir de la magnitud de ondas superficiales
o la magnitud de ondas de cuerpo según lo encontraron Gutenberg y Richter
También es posible relacionar la energía con el momento sísmico o la caída del estrés
según
La forma más precisa de hacerlo es calcular la energía directamente del espectro de
velocidad observado según la relación obtenida por Boatwright y Fletcher (1984)
donde es la densidad en la fuente, es la distancia estación-fuente, y es el
factor de corrección por atenuación intrínseca del cual hablaremos posteriormente.
34
3.5 Escalamiento sísmico y auto-similaridad
En general los parámetros de fuente en un terremoto parecen estar ligados. Al
aumentar el momento sísmico (o el deslizamiento total) aumenta la duración de la
ruptura ; de igual forma el tiempo de ruptura y deslizamiento total son
proporcionales a la longitud de la ruptura . Estas relaciones entre
parámetros de fuente se denominan relaciones de escalamiento, y su validez descansa
en la suposición de que es constate sobre . Las relaciones de escalamiento son:
Donde es el tiempo de ruptura (tiempo que demora en fracturarse la falla), es el
tiempo de ascensión (tiempo que le toma al deslizamiento ir de cero a su valor final),
son las dimensiones de la falla, es la velocidad de ruptura (típicamente
) y son constantes de proporcionalidad.
Note que con estas relaciones el momento se escala según la tercera potencia de la
longitud de la falla , lo cual indica que si doblamos la longitud de la falla,
el momento aumenta en un factor de 8.
La condición de razón de aspecto constante no necesita demostración, ya que se
satisface siempre para cualquier falla rectangular. La condición de deformación
constante es equivalente a la condición de estrés constante, ya que
. La condición de similaridad dinámica establece , pero dado que
, llegamos a constante lo cual equivale a la condición de deformación
constante.
Físicamente la auto-similaridad implica que no existe una escala de longitud o tiempo
intrínseca a la falla, por lo que la dinámica y evolución de de la fuente no serán
funciones de la ninguna dimensión intrínseca a la misma. En los modelos estudiados en
el capitulo anterior, solamente el modelo de de falla elíptica de Dahlen y el de modelo
de falla circular de Madariaga son auto-similares. La auto-similaridad en un modelo
teórico implica que todas sus variables espacio-temporales pueden expresarse de
forma adimensional, mediante el cambio de coordenadas apropiado.
Una vez establecidas las relaciones de escalamiento, los espectros correspondientes a
eventos de diferentes magnitudes forman una única familia, cuyo único parámetro
libre es el momento sísmico . Una vez establecido la frecuencia de esquina
35
queda fija, por lo que la forma del espectro queda completamente determinada
(suponiendo que conocemos la tasa de decaimiento a altas frecuencias).
Figura 3.5. Formas espectrales de eventos de diferentes magnitudes asumiendo escalamiento sísmico en la fuente,
las frecuencias de esquina siguen la líneas (Negro). (Rojo) , (Verde) , (Azul) ,
(Naranja) , (Cyan) .
La figura 3.5 muestra los espectros de terremotos de diferentes magnitudes. Note
como se desplaza la frecuencia de esquina a las bajas frecuencias a medida que la
magnitud del terremoto se hace mayor. Aquí para terremotos con ,
para terremotos con , para terremotos con ,
para terremotos con .
En general las frecuencias de esquina siguen la línea en el diagrama log-log. Esto
se puede hacer intuitivo si consideramos un pequeño ejemplo: Imaginémonos dos
terremotos de fuente idéntica y de longitudes , y funciones fuente en
tiempo , . El área de ruptura será , , el deslizamiento
será , . Por lo que el momento será , . Lo
cual implica a que la funciones fuente-tiempo se escalarán según
Tomando la transformada de Fourier tenemos entonces
Con lo cual
De lo cual se deduce que las frecuencias de esquina siguen la línea en el diagrama
log-log (Prieto, 2004). Véase Figura 3.5.
1.00.5 5.0 10.0 50.0Frequency Hz
108
106
104
uf
1.00.5 5.0 10.0 50.0Frequency Hz
0.10
1.00
0.50
0.20
0.30
0.15
0.70
u f
0
1.00.5 5.0 10.0 50.0Frequency Hz
1.000
0.500
0.100
0.050
0.010
0.005
0.001
u f
0
36
37
Parte II
La segunda parte de este trabajo está dedicada al estudio particular que ocupa el
cuerpo principal de esta tesis. En este estudio evaluaremos las propiedades de
escalamiento de los parámetros de fuente de un área sismológicamente muy activa y
compacta conocida como el nido de Bucaramanga (“BCN” por sus siglas en inglés)
ubicado en el noreste de Colombia. Los datos provienen de la Red Sismológica
Nacional de Colombia (RSNC), de los cuales escogemos únicamente las estaciones
banda-ancha que se encuentren lo más próximo posible al nido.
Utilizamos el algoritmo multi-ventanas para la estimación espectral de las señales de
tiempo en las fases P y S. Con el fin de corregir en los datos las estructuras asociadas a
la presencia de ruido, calculamos el espectro de la ventana pre-P (de decir la ventana
de tiempo justo antes de la llegada de la onda P) para todos los eventos en el conjunto
de datos de la estación. Luego calculamos la distribución de probabilidad de ruido en
cada punto del plano amplitud espectral contra frecuencia. Mediante el ajuste de esta
distribución a una distribución Gaussiana, obtenemos los intervalos de confianza para
una estimación de amplitud 90% libre de ruido en cada una de las frecuencias
presentes en el espectro. Esto establece un límite inferior de corte para los espectros P
y S, tal que estos se encuentren bien arriba de las amplitudes presentes en el ruido.
Luego, después de remover las estructuras espectrales relacionadas con el ruido,
podemos ajustar los espectros resultantes para las fases P y S a un modelo de fuente
de Brune, con un decaimiento asintótico a altas frecuencias de tipo , y extraer de
este modelo las frecuencias de esquina, las cuales serán la clave para el cálculo
posterior de los parámetros de fuente. La energía sísmica se obtiene por integración
directa del modelo de fuente espectral.
Los resultados arrojan una relación lineal entre la caída de estrés y el momento
sísmico, indicando un comportamiento no auto-similar en los terremotos provenientes
del nido de Bucaramanga.
38
1. Introducción
Este estudio evalúa los parámetros de fuente y la hipótesis de escalamiento de un
grupo de terremotos de profundidades intermedias altamente concentrado conocido
como el Nido de Bucaramanga (Colombia). Actualmente existe una falta de consenso
en la comunidad científica con respecto a la hipótesis de auto-similaridad. Algunos
investigadores han encontrado que el estrés aparente se escala con el momento
sísmico , lo cual implicaría que los eventos más grandes son radiadores de
energía sísmica más eficientes que los eventos pequeños (Kanamori, 1993,
Abercrombie, 1995, Mayeda & Walter, 1996, Walter 2006). Esta conclusión rompe la
idea intrínseca de simetría presente en la auto-similaridad, la cual implica imaginar a
los eventos pequeños sencillamente como versiones escaladas de los eventos grandes;
esta conclusión lleva también a la búsqueda del mecanismo físico responsable por este
comportamiento. Otros estudios que apoyan la idea de auto-similaridad, en la cual el
estrés aparente aparece constante sobre el rango de magnitud analizado (típicamente
), son por ejemplo: Ide & Beronza (2001), Prieto (2004), Yamada (2007).
Es importante notar, no obstante, que el rompimiento de la auto-similaridad puede ser
atribuido a las diferencias en la mecánica de falla de los terremotos, tal como lo
sugiere Malagnini (2008). La ruptura de la hipótesis de escalamiento también puede
atribuirse a fallas en los modelos, tales como la subestimación de la energía (Ide y
Beronza, 2003). Somos cuidadosos de evitar este tipo de causas de error en este
estudio.
Los métodos más comunes para la estimación de parámetros de fuente basados en las
llegadas de onda directas (P y S) son el método de constante y el método de razón
espectral.
El método de constante, descansa en la suposición de que el factor de calidad es,
a primera aproximación, independiente de la frecuencia; o por lo menos resulta
independiente para la banda de frecuencia de interés en la sismología (Abercrombie,
1995). La ventaja de este método radica en que puede aplicarse para cualquier a
cualquier terremoto independientemente de la banda de momento presente en el
conjunto, además proporciona una medida absoluta de las amplitudes espectrales a
largo periodo.
Los métodos de razón espectral no hacen ningún tipo de suposición con respecto al
comportamiento del factor de calidad sobre la frecuencia. Su principal ventaja radica
en que pueden tener en cuenta las desviaciones en los sismogramas producto de los
términos cercanos a la estación y los términos de camino; para esto, el método de
razón espectral se basa en la escogencia de pares de terremotos con hipocentros
cercanos dentro del mismo conjunto de eventos. Una vez escogidos los pares
apropiados se realiza una razón espectral entre el espectro del terremoto más grande
39
y espectro del terremoto más pequeño, esperando que muchas de las complejidades
propias del camino de propagación (las cuales típicamente se describen con factores
multiplicativos en los espectros) se eliminen. Dado que se quiere poder resolver las
frecuencias de esquina de los eventos escogidos, la diferencia entre las frecuencias de
esquina de los terremotos debe ser grande, de modo que los eventos deben estar
separados en magnitud. Típicamente debe dejarse un orden de magnitud entre ellos
con el fin de poder realizar una estimación correcta de los parámetros de fuente.
Desafortunadamente, la aplicación de este tipo de análisis a nuestros datos resulta
imposible debido a la cantidad comprometedora de ruido observada en los eventos de
magnitudes pequeñas. Escoger eventos de magnitudes para eliminar las
complejidades de camino en los eventos más grandes eliminaría la mayor parte de los
eventos claros disponibles en nuestras bases de datos ( [3.0, 5.0]), dejándonos
únicamente alrededor de 80 eventos para el análisis. Además, el intervalo de
frecuencia disponible para el estudio de estos eventos sería muy pequeño, debido a la
reducida banda de frecuencia libre de ruido en los terremotos de menor magnitud.
Otros métodos para estudios de fuente han sido también desarrollados a partir del
análisis de las ondas coda (Baltay, 2009), las cuales son ondas secundarias producto de
múltiples reflexiones en el medio de propagación. Sus llegadas son posteriores a las
llegadas de las ondas directas y se presentan con amplitudes mucho menores, por lo
que se requiere una buena resolución instrumental para estudiarlas.
Las ondas coda proveen un método robusto para la estimación de parámetros
espectrales, permitiendo la corrección empírica de los efectos de sitio y de camino, con
efectos reducidos de directividad y patrón de radiación (Mayeda & Walter, 1996,
Baltay et al., 2009,2011). Estas propiedades, que son siempre un problema en los
métodos convencionales, las hacen idóneas para el estudio de los parámetros de
fuente sísmica.
Desafortunadamente, este método no puede ser aplicado a nuestro conjunto de datos,
dado que las señales coda no tienen la suficiente claridad para identificarlas de forma
fidedigna en nuestras señales de tiempo.
En este estudio adoptamos el método de factor de calidad constante (Abercrombie,
1995) para la estimación de parámetros de fuente sísmica.
40
2. Datos
Utilizamos las señales provenientes de la red sismológica nacional, las cuales se
obtienen de estaciones banda-ancha asentadas sobre lechos rocosos, lo más cerca
posible del nido de Bucaramanga. Las características de las estaciones utilizadas son8:
Tabla 2.1. Características de las estaciones utilizadas
Nombre Estación
Ubicación Depto Longitud (Wº)
Latitud (Nº)
Altitud (msnm)
Tipo Digitalizador Sensor
BRR Barrancabermeja Santander -73,71 7,11 137 Banda Ancha
QUANTERRA Q330
STRECKEISEN STS-2
HEL Santa Helena Antioquia -75,53 6,19 2815 Banda Ancha
QUANTERRA Q330
GURALP CMG-3T
RUS Cerro La Rusia Boyacá -73,08 5,89 3697 Banda Ancha
GURALP DM 24
GURALP 3T
CHI Páramo de Chingaza
Cundinamarca -73,73 4,63 3140 Banda Ancha
QUANTERRA Q330
STRECKEISEN STS-2
PRA Represa de Prado
Tolima -74,89 3,71 468 Banda Ancha
QUANTERRA Q330
STRECKEISEN STS-2
En general estas estaciones demostraron la obtención de buenas razones señal/ruido,
con una frecuencia de muestreo de 100 Hz ( ) en medidas de velocidad.
Además de una cantidad de eventos detectados en BCN suficiente para el
procesamiento.
A pesar de que la banda de frecuencia de los sismómetros va de a , en
la práctica solamente una fracción de esta banda de frecuencia es útil para el estudio.
La razón de este inconveniente: posibles deficiencias en la instalación y asentamiento
de las estaciones. Desafortunadamente para casi todos los eventos provenientes de la
RSNC, existen estructuras espectrales a bajas frecuencias las cuales terminan por
apoderarse de la señal nublándola completamente bajo una estructura de escalón. En
general el ruido afecta con mayor fuerza a los eventos pequeños en la base
de datos. La figura 2.1 muestra claramente la presencia de estas estructuras alrededor
de los .
Adicionalmente, la caída del espectro a altas frecuencias lo obliga a atravesar el nivel
esperado del ruido, haciendo la señal inservible en el rango .
El efecto conjunto de estas dos limitantes deja un rango limitado de frecuencia de
aproximadamente. En la práctica la banda de frecuencia útil varía de
evento en evento en proporción con su magnitud.
La figura 2.1 muestra los espectros de velocidad calculados para las ondas P y S de 7
eventos en la estación BRR durante el 2008. Las líneas rojas muestran el espectro del
ruido presente en la estación justo antes de la llegada de las primeras ondas de
8 Datos obtenidos a partir de los registros de la Red Sismológica Nacional de Colombia. http://seisan.ingeominas.gov.co/RSNC/
41
cuerpo, y las líneas negras representan los espectros de las ondas P y S de los
diferentes eventos. Obsérvese como los eventos más pequeños están muy próximos al
nivel de ruido, o incluso por debajo del ruido en buena parte de la banda de
frecuencia. Nótese la estructura escalonada en el espectro a bajas frecuencias tanto en
el espectro de la señal, como en el espectro del ruido. Esto es evidencia directa de que
tales estructuras no son de origen sísmico; sino más bien, producto de efectos propios
a la estación, y por lo tanto deben ser apropiadamente removidos con el fin de evitar
sesgos en el análisis.
Figura 2.1. Espectros de velocidad para 7 eventos en el nido de Bucaramanga, ocurridos en el 2008 y
detectados en la estación BRR. La señal se muestra en negro mientras que el ruido aparece en rojo. Aquí
se muestran tanto la fase P (izquierda) como la fase S (derecha).
Únicamente los eventos en el interior del nido de Bucaramanga son considerados para
el procesamiento. Definimos la ubicación del nido de Bucaramanga entre
de latitud, los de longitud y los 120-180 km de profundidad.
3. El nido de Bucaramanga y la sismicidad profundidad
intermedia
El nido de Bucaramanga (BCN) es una región ubicada en el departamento de Santander
sísmicamente muy activa, y altamente localizada en la cual ocurren continuamente
terremotos cuyos hipocentros se encuentran localizados a distancias intermedias
(>60km) por debajo de la superficie. Con alrededor de 15 terremotos por día
concentrados a profundidades de 140 km en una región de alrededor de 15 km2, el
nido de Bucaramanga es único en el mundo; existiendo pocas zonas en el mundo que
se puedan comparar a su alta tasa de actividad y gran compacidad. Este nido de
terremotos se encuentra ubicado alrededor de los , 160 km de
profundidad (ver figura 3.1).
Estudios tomográficos en esta región (Pennington, 1981) han demostrado que el nido
de Bucaramanga es el producto de la colisión de la placa de Maracaibo y la placa de
42
Bucaramanga, donde esta última subduce la primera (Taboada et al, 1998, Zarifi,
2007).
Los mecanismos focales en el nido siguen un claro patrón acorde parcialmente con su
geología, mostrando subducción al oeste en el eje principal de alineación de la presión
(eje P) y subducción al este en el eje principal de alineación de la tensión (eje T).
Aunque el eje T está sobre el plano de falla, el eje P no lo está. Además mediante el
estudio de sus tensores de momento, se ha demostrado que corresponden a fuentes
de doble pareja casi perfectas (Frohlich, 1995, Cortés, 2005).
Se han encontrado tasas de deformación muy altas en esta región, de alrededor de
a , o entre y en todo el nido. Dado que
las tasas de deformación en el nido de Bucaramanga son de uno a tres órdenes de
magnitud mayores que las presentes en los terremotos superficiales, no es
sorprendente que se desarrolle una tasa alta de actividad en esta región. La dimensión
característica de la fuente ha demostrado ser pequeña en comparación con la de los
terremotos superficiales, alrededor de los 4 km (Frohlich, 1995).
Con respecto a la distribución estadística de sus terremotos, el nido de Bucaramanga
muestra una clara escasez de terremotos con , lo cual se traduce en una
estadística de Gutenberg-Richter anómala (Eysteinn & Edward, 1970). Estudios
posteriores (Frohlich, 1995) arrojan un valor de para la pendiente de
la distribución momento-frecuencia, lo cual sin duda equivale a un valor alto en
comparación con estadísticas provenientes de otros sitios. Este valor de representa
una dimensión fractal de 2.5 a 2.8 correspondiente a una geometría volumétrica9.
La figura 3.1 (b) muestra la posición relativa de 2000 eventos sísmicos detectados en el
nido de Bucaramanga durante el 2008. Puede verse claramente lo pequeña y
compacta que resulta la geometría del nido, así como su clara concentración alrededor
de los 140-160 km de profundidad. En la figura 3.2 se puede apreciar también un corte
transversal del nido.
9 equivale a una falla lineal, equivale a una falla planar, y
corresponde a una geometría volumétrica de la fuente (Frohlich, 1995).
43
Figura 3.1. (a) Localización del nido de Bucaramanga (círculos negros) junto con las estaciones de la Red
Sismológica Nacional de Colombia. (b) Localizaciones de los hipocentros dentro del nido de
Bucaramanga. Aquí mostramos alrededor de 2000 eventos sísmicos.
Figura 3.2. Corte transversal del nido de Bucaramanga sobre los 7º de latitud. Aquí mostramos alrededor
de 2000 eventos sísmicos.
La sismicidad profunda presente en el nido de Bucaramanga y en otras regiones a lo
largo del mundo resulta todavía intrigante para la comunidad científica. Sabemos que
los terremotos profundos se caracterizan por ocurrir a profundidades mayores a 60
km. En contraste con los terremotos superficiales, estos suelen tener magnitudes
menores, caídas de estrés más grandes y funciones fuente-tiempo más impulsivas, así
como muy pocas réplicas. Aún así, existen ejemplos de terremotos profundos que
contradicen estas ideas.
Consideremos momentáneamente una falla a 140 km de profundidad. Para producir
deslizamiento en esta falla es necesario que el estrés cortante en la falla supere las
fuerzas de fricción sobre la misma, la cuales se oponen al deslizamiento. No obstante si
hacemos lo cálculos, y multiplicamos la presión a tales profundidades por el
coeficiente de fricción esperado para las rocas, obtendremos una medida de estrés
gigantesca, que supera con creces el estrés que soportan las rocas. Todo indica que
tales eventos sencillamente no deberían ocurrir. Sin embargo suceden.
(a) (b)
44
Es posible que el coeficiente de fricción no sea el mismo a tales presiones y
temperaturas. Otra alternativa que abriría la posibilidad de la ocurrencia de
terremotos profundos en un marco conceptual es la presencia de rocas porosas en el
plano de falla, tal que sus poros están llenos de algún tipo de líquido. Este líquido
puede ser simplemente agua o algún fluido de naturaleza magmática. En este caso la
fase líquida establece una presión negativa que lleva a la posibilidad de deslizamientos
debido a una fricción reducida.
Sea cual sea el caso, no hay duda de que los terremotos profundos continúan
representando uno de los campos activos de investigación en la sismología. La
comprensión de la mecánica de su fuente será de gran utilidad a la hora de contrastar
los modelos teóricos que se pueda proponer con datos experimentales precisos.
4. Cálculo de la distribución de ruido
Una estimación precisa de parámetros de fuente recae en una estimación robusta de
las frecuencias de esquina por medio de ajustes del espectro a modelos teóricos. Los
métodos comunes de ajuste involucran algoritmos de mínimos cuadrados o mínimos
cuadrados robusto, los cuales frecuentemente introducen la suposición de que la
desviación de los datos con respecto al modelo está causada por ruido Gaussiano o
casi Gaussiano. Este no es precisamente el caso para las señales de tiempo sísmicas, en
las cuales el espectro del ruido tiene estructuras grandes alrededor de los 8s producto
de patrones de interferencia entre las ondas oceánicas, alrededor de los 15s producto
de las olas rompiendo en las costas, y alrededor de las 12h/24h producto de las
mareas.
Para una estimación más confiable de la distribución del ruido en nuestro conjunto de
datos, adoptamos un esquema similar al utilizado por McNamara & Bunland (2004).
Bajo este esquema, tomamos para cada evento una ventana de tiempo de 16s de
longitud antes de la llegada de la onda P y calculamos su espectro utilizando el método
multitaper (Park, 1987).
La figura 4.1 muestra una señal de tiempo típica proveniente del nido de
Bucaramanga. La longitud y posición de las ventanas pre-P, P y S en tiempo se
muestran en las líneas de colores rojo, verde y azul, respectivamente. Las ventanas de
tiempo P y S serán de utilidad posteriormente en nuestro análisis
45
Figura 4.1. Señal de tiempo típica detectada en la estación BRR en el 2008 para un evento de magnitud
4.7. Se muestran las ventanas de tiempo pre-P (rojo), P (verde), y S (azul). La amplitud se mide en
conteos.
A partir de estos registros calculamos la frecuencia de ocurrencia de ruido en cada
punto del plano Log(Frecuencia) - Log(Amplitud Espectral). Para realizar este cálculo el
plano es dividido en celdas de área tal que una distribución de
frecuencias suave se obtiene para ambas variables. Dado que nuestro interés se basa
en la distribución de ocurrencia de ruido en función de la frecuencia, calculamos la
frecuencia de ocurrencia de ruido por cada línea de frecuencia.
La figura 4.2 muestra el plano discretizado mediante una malla que lo
recubre. Cada punto del espectro se aproxima a alguno de los nodos de la malla,
superponemos todos espectros a la red, asegurándonos de contar cuantas veces se
aproxima un punto del espectro a un determinado nodo de la red, cada punto del
espectro equivale a un conteo. Posteriormente se cuentan la cantidad total de conteos
sobre cada punto de la malla. Fijándonos en una línea de frecuencia constante, es
posible graficar la cantidad de conteos obtenidos sobre todo el rango de amplitud.
Podemos convertir estos conteos a probabilidad de ocurrencia si normalizamos los
conteos obtenidos en cada intervalo de amplitud por la cantidad total de conteos
sobre toda la banda de amplitud.
tiempo (s)
Amplitud
46
Figura 4.2. Diagrama esquemático de la discretización del plano . El ancho de las celdas
se ha exagerado para facilitar su visualización. Adicionalmente se muestran algunos 42 espectros de
ruido pre-P en la estación RUS, los cuales son la base para la estimación de la densidad de amplitud de
ruido en esta estación.
Figura 4.3. Diagramas de frecuencia de ocurrencia de amplitud de ruido para bandas de frecuencia
centradas en (a) , (b) y (c) y (d) . Note el cambio en la
distribución del ruido a medida que aumenta la frecuencia. (azul) datos, (rojo) ajuste gaussiano.
Una vez calculada la distribución de probabilidad de ocurrencia, hacemos una
aproximación al continuo para derivar de esta la densidad de probabilidad de
ocurrencia. Para esto es necesario asegurar que la integral de la densidad sobre todo el
rango de amplitudes sea siempre 1. Es necesario realizar este procedimiento para
todas las frecuencias disponibles, generando así una curva densidad de probabilidad
por frecuencia disponible.
En general, la probabilidad se muestra más dispersa a bajas frecuencias y tiene a
concentrarse mucho más a medida que nos movemos a altas frecuencias. Esto se
muestra en la figura 4.3, en la cual se muestran varios cortes de frecuencia contante.
La figura 4.4 muestra la densidad de probabilidad para el espectro de ruido la estación
BRR. Esta figura se obtiene mediante la agrupación de todas las curvas densidad de
probabilidad.
(a) (b)
(c) (d)
47
Figura 4.4. Densidad de probabilidad de ruido en la estación RUS. Se muestra una sección del plano
con el fin de observar el comportamiento del ruido en función de la frecuencia. La escala
de amplitud de conteos por segundo y la de frecuencia es Hertz.
Realizamos un ajuste a una densidad de probabilidad Gaussiana para cada línea de
frecuencia constante. Para esto utilizamos un algoritmo de mínimos cuadrados que
nos permita obtener la media y la desviación estándar de cada una de las
densidades de probabilidad. La línea de confianza definida por el conjunto de todos los
puntos en el plano corresponde a la línea de amplitudes libres de ruido con
un 90% de probabilidad.
La figura 4.5 muestra las distribuciones de probabilidad estimadas, junto con las líneas
correspondientes a la media (rojo) y los intervalos de confianza (verde)
correspondientes a las líneas . A modo de comparación, se grafican también
los espectros de las ondas directas. Note que en algunos puntos los espectros se
encuentran por debajo de la región “libre de ruido” .
Los intervalos de confianza nos permiten recortar los puntos espectrales para los
cuales la condición de estar libre de ruido a un 90% de probabilidad no está asegurada.
Estos puntos están ubicados por debajo de la línea de confianza , y son
removidos de los datos. Este método permite también la remoción de estructuras
espectrales asociadas al instrumento y posibles resonancias de sitio, dejando
únicamente señales de datos reales con poca probabilidad de incidencia de ruido.
48
(a) (b)
(c) (d)
(e) (f)
49
(g) (h)
Figura 4.5. Densidad de probabilidad de ruido para las estaciones (a,b) BRR (c,d) HEL (e,f) RUS (g,h) CHI
en las ondas P y S. Se pueden ver claramente las líneas de confianza (verde) y la línea de la
media (rojo) de la distribución. Se grafican también los espectros de las ondas de cuerpo. Toda región
del espectro por debajo de la línea de confianza será recortada.
5. Procesamiento de datos
5.1 Eliminación de estructuras indeseadas
Ventanas de tiempo disjuntas son utilizadas para extraer de los registros las ondas P y
S de cada evento. Utilizamos una longitud de ventana constante de 16s de largo con el
fin acceder a las frecuencias más bajas posibles. Ventanas más grandes o más
pequeñas (entre 7s y 16s) no parecen afectar significativamente la razón señal/ruido
de nuestros datos. Tomamos ventanas de tiempo que comienzan 0.5 s antes de las
llegadas de las ondas directas para evitar efectos de borde en la estimación espectral.
Figura 5.1. Espectros de onda P (izquierda) y onda S (derecha) para un evento de magnitud 3 detectado
en la estación HEL. Nótese las estructuras en forma de escalón a bajas frecuencias incompatible con la
estructura de cuña esperada para un espectro de velocidad.
50
La densidad espectral de potencia para los registros de onda P, y S (en sus 3
componentes SE, SN y SZ) son calculados a partir del algoritmo multitaper (Park &
Vernon, 1987), el cual es conocido por presentar la mayor concentración de energía
sobre cualquier banda de frecuencia, sin efectos considerables de filtración espectral o
sesgo. Las amplitudes espectrales son calculadas luego de la raíz cuadrada del espectro
de potencia. La amplitud espectral total presente en la onda S se calcula a partir de la
suma vectorial de las amplitudes en sus componentes , y mediante la suma
vectorial . En contraste, la amplitud espectral de la onda P se
calcula únicamente a partir del componente ya que en nuestros datos más del 90%
de la energía se concentra aquí.
Antes de comenzar con el procesamiento, necesitamos remover las estructuras
indeseables en nuestros datos. La primera de estas correcciones corresponde a la
remoción de las estructuras no reales en el espectro de bajas frecuencias, las cuales
provienen de posibles defectos en la instalación de las estaciones, niveles altos de
ruido sísmico en Colombia, o ruido proveniente de la cordillera oriental. Estas
estructuras presentan típicamente unas desviaciones en forma de escalón en los
espectros, los cuales siguen una línea recta en la grafica
, luego de agrupar varios de estos espectros en la misma grafica. Solo
los puntos espectrales por debajo de esta línea son eliminados de los datos (véase fig.
5.2).
Figura 5.2. Línea de corte a bajas frecuencias (rojo) para los registros de densidad espectral (negro) en la
estación BRR fase P. Solamente se dejan los datos por encima de la línea de corte, los cuales
corresponden a la forma espectral esperada para medidas de velocidad. Los valores de corte y de
pendiente que definen estas líneas en los ejes log-log se encuentran en la Tabla 5.1.
51
Debido a que la línea de recorte puede variar de estación a estación, es necesario
calcularla para cada caso particular. La tabla 5.1 muestra los valores de la pendiente
y el corte de cada una de las estaciones y cada una de las fases P y S .
Tabla 5.1. Valores de las pendientes y los interceptos verticales de las líneas de corte en cada estación.
Las líneas se calculan en los ejes log-log de densidad espectral contra frecuencia.
Estación
HEL -4.77 0.66 -8.59 -0.56
BRR -4.26 1.00 -4.32 1.43
CHI -4.19 0.78 -6.11 0.26
RUS -3.30 0.85 -7.52 0.58
Luego de eliminar las estructuras no reales a bajas frecuencias, se hace también
necesario limitar la señal únicamente hasta los 20Hz de frecuencia. Esto debe hacerse
así con el fin de evitar muchos de los efectos de sitio presentes en el espectro de altas
y medianas frecuencias. Estas contribuciones provienen del lugar en particular en el
cual está ubicado el sismógrafo, y se deben generalmente a resonancias y/o bandas de
absorción debido a la geología del sitio en particular donde se encuentra la estación.
Utilizamos también la línea de confianza al 90% de probabilidad calculada
anteriormente para recortar del espectro los puntos ruidosos. Este procedimiento nos
asegura que cada punto que permanezca en el espectro tiene un 90% de probabilidad
de representar la señal y no el ruido.
Es necesario tener cuidado aquí ya el recorte sobre la línea de ruido puede producir en
algunos casos, puntos aislados los cuales no siguen la tasa de caída promedio del
espectro, especialmente en las frecuencias más altas. Dejar estos puntos en el
espectro puede sesgar el ajuste espectral posterior. Para evitar este problema
aplicamos a nuestros datos un algoritmo de eliminación de puntos de baja densidad,
con el objetivo de eliminar aquellos puntos aislados del espectro. Estos puntos se
caracterizan por su baja densidad de número en un intervalo pequeño de frecuencia
alrededor de ellos. Luego de la eliminación de estos, quedan en nuestros espectros
solamente los puntos altamente precisos y prácticamente libres de ruido. El algoritmo
elimina los datos en los cuales no existe una suficiente densidad de puntos alrededor
para que estos sean considerados como conexos con el resto de la secuencia.
La figura 5.3 muestra estas líneas de recorte para los eventos detectados en la estación
CHI.
52
Figura 5.3. Líneas de recorte para los eventos detectados en la estación CHI en los años 2007-2010. La
grafica muestra las regiones eliminadas en colores, las regiones de color negro corresponden a aquellas
por debajo del 90% de confianza, las de color verde corresponden a aquellas que fueron consideradas
“error instrumental”, las regiones de color azul se eliminaron porque contenían posiblemente efectos
de sitio (nótese el fuerte efecto de sitio en esta estación). Dejando entonces las regiones de color
morado para el procesamiento.
Finalmente debemos corregir en el espectro el sesgo a altas frecuencias producto de la
atenuación intrínseca. Como vimos en la sección 3.3 de la Parte I, esta atenuación se
puede modelar de la forma
Donde es el espectro observado de velocidad y es el espectro fuente de
velocidad. La corrección del espectro, para obtener el espectro fuente, consiste
únicamente en multiplicar el espectro observado por . Con el tiempo de
vuelo y el factor de calidad. Explicaremos la escogencia de este valor de
en la próxima sección.
La figura 5.4 muestra el resultado de la corrección de atenuación para un evento de
magnitud detectado en la estación RUS. Nótese el cambio significativo en la
región de altas frecuencias.
53
Figura 5.4. Comparación entre los espectros sin corrección por atenuación (negro) y con corrección por
atenuación (azul); tanto para ondas P (a) como para ondas S (b). Ambos espectros provienen de un
único evento de magnitud detectado en RUS. Si no se corrigen los efectos de la atenuación,
podrían obtenerse espectros con frecuencias de esquina menores a las reales. Aquí =1317.
5.2 Escogencia del factor de calidad
La escogencia del factor de calidad resulta de fundamental importancia para la
determinación de los parámetros de fuente. Una escogencia errónea de este
parámetro desviará considerablemente el espectro en la región de altas frecuencias, lo
cual llevaría a una estimación incorrecta de las frecuencias de esquina.
Existen varios estudios que determinan el valor de justo al momento de realizar el
ajuste al modelo (Gök & Hutchings, 2009). En este estudio evitamos esta aproximación,
ya que existe una relación de intercambio entre y . Es posible reproducir de forma
aproximada un mismo espectro ya sea aumentando su factor de calidad y
disminuyendo su frecuencia de esquina o aumentando su frecuencia de esquina y
disminuyendo su factor de calidad. La implementación de esta relación de intercambio
lleva a valores variables, aún para espectros estudiados en una única estación. Esto
es contradictorio, ya que se espera que para caminos de propagación similares el
factor de calidad sea aproximadamente igual. Además de esto se podrían encontrar
frecuencias de esquina inusualmente más grandes que las esperadas intuitivamente.
Para evitar este inconveniente decidimos fijar el valor de desde el inicio, y utilizar el
mismo valor para corregir todos los eventos. Para obtener el valor especifico que se
debe usar, imaginémonos un único terremoto observado en dos estaciones diferentes.
Los espectros de velocidad , recibidos para una determinada fase (P
o S) en ambas estaciones serán entonces (Boatwright, 1980):
(a) (b)
54
Donde las fórmulas para y son las amplitudes a bajas frecuencias; , son los
tiempos de vuelo hasta cada estación, y es la frecuencia de esquina del evento. La
razón espectral será
Despejando tenemos
Para resolver , utilizamos las fórmulas de la sección 6.2 y suponemos una distancia
estación-evento dada por y , con la velocidad de onda.
A partir de este tipo de procedimiento realizado sobre varias estaciones y sobre todos
los eventos disponibles obtenemos una familia de curvas . Estas curvas suelen
presentar un comportamiento constante sobre , con algunos picos muy pronunciados
debido a las imprecisiones numéricas producto de la división . Para
eliminar la contribución fuertemente localizada de estos picos, tomamos la mediana
de cada una de las curvas . Luego calculamos la media de todas las medianas para
obtener un valor de . Este valor representa el factor de calidad promedio
para la región de interés, y lo utilizamos para corregir todos los espectros. Resulta
coherente desde el punto de vista físico, ya que se espera una valor alto de este
parámetro en el manto superior , el cual es débilmente atenuante. En
contraste con los resultados obtenidos para terremotos superficiales, en los cuales en
promedio (Ide, 2003, Tomic, 2009).
55
6. Estimación de parámetros de fuente
6.1 Ajuste al modelo
La fuente sísmica puede ser descrita en términos de parámetros de fuente tales como
las frecuencias de esquina , , la dimensión de la fuente , el momento sísmico
, la energía sísmica liberada , y la caída de estrés . Con el modelo apropiado de
fuente, todas estas cantidades nos brindan una descripción física de la mecánica de la
fuente. Esto puede brindar información importante sobre el comportamiento de los
terremotos, en particular en el nido de Bucaramanga.
Utilizamos aquí el análisis de constante para el modelamiento de los parámetros de
fuente (Boatwright, 1980). El espectro de velocidad con correcciones previas de
atenuación intrínseca se modela entonces como
en donde es la frecuencia de esquina, es el factor de calidad, es el tiempo de
vuelo de la onda directa (ya sea la onda P o la onda S), y es la amplitud a bajas
frecuencias.
Utilizamos aquí el modelo de decaimiento asintótico propuesto por Brune (1970),
ya que es el que mejor se ajusta a los datos, los cuales muestran una caída ligera a
altas frecuencias y una frecuencia de esquina suave. Este tipo de modelo está también
de acuerdo con la caída a altas frecuencias esperada según el modelo de fuente
circular de Madariaga (1976), el cual utilizamos para describir nuestros datos.
Para la estimación de y , los datos espectrales de las ondas P y S de los eventos
son ajustados al modelo teórico por medio de un algoritmo de mínimos cuadrados. El
tiempo de vuelo de cada una de las fases se obtiene directamente del catálogo.
Una buena estimación de parámetros de fuente se basa en una estimación precisa de
las frecuencias de esquina. Para evitar desviaciones en el cálculo de las frecuencias de
esquina, es necesario asegurarse de realizar el ajuste al modelo teórico en el eje lineal
de frecuencia; en contraste con los ejes log-log estándar que se usan al mostrar los
espectros. Forzar el algoritmo de regresión a trabajar con el logaritmo de la frecuencia,
le daría más resolución a las bajas frecuencias y mucha menos resolución a las altas. Es
importante tener en cuenta este detalle, ya que una regresión sobre un eje de
frecuencia logarítmico resultaría en frecuencias de esquina gigantes para algunos de
los eventos, producto de un error numérico mucho mayor. (ver figura 6.1)
56
Figura 6.1. Distribución simétrica de puntos sobre un eje lineal (arriba) y su correspondiente
representación en un eje logarítmico (abajo). Como se puede ver, el eje logarítmico le da especial
relevancia a las bajas frecuencias (ej: la región [0.1, 1]) mientras que comprime las frecuencias más
altas (ej: la región [10, 25]). Hacer una estimación sobre un eje logarítmico generará imprecisión en el
algoritmo de interpolación para ciertos intervalos de frecuencia. Esto implica, por ejemplo, que en el eje
logarítmico una frecuencia de esquina de se resolverá con mayor resolución que una
frecuencia de esquina de
Para evitar imprecisiones que pueden resultar luego en grandes errores en los
parámetros de fuente, optamos por realizar la interpolación en los ejes log-lin.
La figuras 6.2 y 6.3 muestran algunos espectros fuente obtenidos para las ondas P y S
junto con sus respectivos modelos teóricos. Para una mejor comparación visual se
muestran los espectros de un mismo evento detectado en las 5 estaciones. Nótese la
presencia de un efecto de sitio notorio en la estación BRR (figura 6.2 e y j, figura 6.3 e y
j ). A pesar de este efecto, las frecuencias de esquina obtenidas en esta estación no
resultan demasiado diferentes a las frecuencias de esquina obtenidas en otras
estaciones, como se puede comprobar visualmente en las figuras a continuación.
57
Figura 6.2. Se muestran varios espectros de onda (a) P y (b) S correspondientes a un único evento junto
con sus modelos teóricos más afines en las 5 estaciones. La frecuencia de esquina estimada del espectro
se muestra en cada grafica. La magnitud local del evento es de 3.7
(a) (b) (c)
(d) (e)
(f) (g) (h)
(i) (j)
58
Figura 6.3. Se muestran varios espectros de onda (a) P y (b) S correspondientes a un único evento junto
con sus modelos teóricos más afines en las 5 estaciones. La frecuencia de esquina estimada del espectro
se muestra en cada grafica. La magnitud local del evento es de 2.9.
Hasta ahora hemos realizado el análisis de datos de forma independiente para todas
las estaciones. No obstante no podemos basar nuestras conclusiones en resultados
individuales a cada estación, ya que existen muchos efectos que hacen que las señales
sísmicas detectadas en una estación puedan resultar diferentes a las señales
provenientes de otra estación. Entre estos se encuentran efectos como la directividad,
el patrón de radiación, los efectos de sitio, o los efectos de camino de propagación.
Para una mejor estimación de las frecuencias de esquina, buscando al mismo tiempo
minimizar la influencia de los efectos ajenos a la fuente, optamos por realizar un
estudio conjunto a todas las estaciones (5 estaciones). Para esto, escogemos
(a) (b) (c)
(d) (e)
(f) (g) (h)
(i) (j)
59
únicamente los eventos que sean comunes a todas las estaciones y luego de calcular su
frecuencia de esquina en cada estación, tomamos la mediana de las frecuencias de
esquina disponibles para cada evento. Esto nos produce una única frecuencia de
esquina por evento, la cual será utilizada para los cálculos posteriores. Realizamos el
mismo procedimiento para las fases P y S independientemente.
La figura 6.4 muestra las frecuencias de esquina obtenidas para cada evento en
función de su momento sísmico. Se muestran también las líneas de caída de estrés
constante correspondientes a , y . Como
puede verse, las frecuencias de esquina no parecen seguir las líneas de caída de estrés
constante, como lo supondría la hipótesis de auto-similaridad. Un ajuste lineal a las
curvas muestra una dependencia de la forma y . Lo
cual equivale a una desviación de la auto-similaridad de
para las ondas P y para las ondas S.
Figura 6.4. Frecuencias de esquina promedio para las ondas P (a) y S (b) de los eventos comunes a todas
las estaciones. Las barras de error se obtienen a partir del estimado de la varianza en cada punto. La
línea morada muestra la regresión lineal de las frecuencias de esquina. Las líneas rojas muestran la
tendencia esperada si las frecuencias de esquina siguiesen las líneas de caída de estrés constante.
La razón de las frecuencias de esquina promedio obtenidas para las ondas P y S
se aproxima al valor esperado según el modelo de Madariaga, el cual
predice (Madariaga, 1976). Nuestros resultados muestran
mostrando un porcentaje de diferencia10 del 11% con respecto al
modelo teórico. Es importante resaltar que esta diferencia no necesariamente
corresponde a una falla del modelo o del análisis realizado. Corresponde a un
resultado real, que puede distanciarse de los modelos a menudo idealizados utilizados
para describirlo. Esta diferencia podría corresponder por ejemplo a una fuente cuya
geometría no es perfectamente circular, o cuya velocidad de ruptura no es
necesariamente .
10
(a) (b)
60
La figura 6.5 muestra la relación entre las frecuencias de esquina para los eventos
analizados
Fig 6.5. Relación entre y para los 80 eventos analizados. La línea roja muestra el valor
estimado por Madariaga para la razón de las frecuencias de esquina correspondientes a una falla
circular, de velocidad de ruptura , y cuyas frecuencias de esquina se encuentran promediadas
sobre toda la esfera focal. La línea azul corresponde a la tendencia observada en los datos.
Por otro lado, la dimensión de la fuente puede ser determinada inmediatamente
luego de la estimación de la frecuencia de esquina, para lo cual Madariaga (1976)
predice
donde es una constante ( , ), es la velocidad de
las ondas S, y es la frecuencia de esquina. Para la estimación de la dimensión de la
fuente usamos ambos estimadores del tamaño de la fuente y promediamos sus valores
La figura 6.6 muestra los resultados para en como función del momento sísmico. Los
resultados demuestran concordancia con los valores esperados para este tipo de
terremotos. Teniendo los eventos dimensiones de 100m a 300m dependiendo de su
magnitud. Las rupturas de longitudes más grandes están asociadas a mayor energía
sísmica liberada y mayor deslizamiento total, luego es de esperar que exista una
proporcionalidad entre el tamaño del terremoto y su dimensión característica. Como
vemos, esta relación entre y es confirmada por los datos.
61
Figura 6.6. Dimensión de la fuente en función del momento sísmico . Se puede apreciar una
tendencia lineal la cual implica que los eventos más grandes presentan mayores momentos sísmicos.
La incertidumbre en está directamente asociada a la estimación en y a la
estimación en la velocidad de onda cortante . Por su parte, se calcula suponiendo
una velocidad de onda P en la fuente de y luego (Ojeda &
Havskov, 2001).
El momento sísmico se calcula aproximando la magnitud local de catálogo a la
magnitud de momento según y utilizando la relación de Kanamori (1977)
Una vez estimadas las frecuencias de esquina, la caída de estrés puede también ser
calculada directamente como lo demostró Elsheby (1957) para una falla circular
La figura 6.7 muestra la caída de estrés calculada. De la figura vemos que la caída de
estrés parece escalarse con el aumento del momento sísmico, en concordancia con la
hipótesis de no auto-similaridad. Las barras de error corresponden a la propagación de
la incertidumbre presente en y . Calculamos mediante una interpolación lineal
(línea roja) la pendiente de escalamiento, la cual nos arroja una relación .
Lo cual muestra una clara desviación del comportamiento auto-similar ( constante).
Los intervalos de confianza al 95% en la interpolación se muestran con líneas moradas
punteadas.
62
Figura 6.7. Caída del estrés en función del momento sísmico para 80 eventos comunes a todas las
estaciones. La línea roja corresponde a una interpolación lineal de los datos, a partir de la cual se puede
obtener la pendiente de escalamiento.
Las incertidumbres presentes en la estimación de resultan mayores que las
presentes en la estimación de , y son producto de los errores en cuya precisión
está limitada a la precisión del catálogo, y a la propagación de la incertidumbre11 en
(Prieto et al, 2007).
6.2 Calculo de la energía sísmica
Resulta también una buena práctica contrastar la claridad de los resultados en con
el cálculo del estrés aparente , el cual debe reflejar una estructura de
puntos similar a la de , con una constante de escalamiento vertical relacionada con
la eficiencia sísmica del proceso de ruptura. Esto nos debe dar una medida de que
tanta energía presente al momento de la ruptura es convertida a ondas sísmicas.
Para realizar un test de la auto-similaridad en la fuente sísmica es común analizar la
variación de la caída de estrés aparente en función del momento .
Donde por supuesto, la energía también se escala con el momento sísmico. Para que la
hipótesis de auto-similaridad se verifique es necesario que sea constante sobre el
momento sísmico tal que todos los terremotos radien con la misma eficiencia sísmica,
independientemente de su tamaño.
11
La incertidumbre se calcula mediante propagación de errores, según la fórmula general:
con las incertidumbres en y , y la
incertidumbre en .
63
Podemos calcular la amplitud a bajas frecuencias en el espectro mediante la fórmula
estándar (Aki & Richards, 1980)
Donde es la densidad, es la velocidad sísmica (ya sea para ondas
P o para ondas S), es el patrón de radiación promedio (0.52 para ondas P y
0.63 para ondas S) y es la distancia de la fuente al receptor (calculada a partir del
tiempo de vuelo y la velocidad de onda sísmica). Estimamos a partir de y no del
modelo más aproximado a los datos, ya que hemos observado variaciones sistemáticas
en los registros de las amplitudes a bajas frecuencias entre años consecutivos. Estas
variaciones nos llevan a imaginar posibles cambios de instrumentación y/o re-
calibración de los instrumentos existentes en algunas de las estaciones. Haciendo
dudosas las dimensiones de amplitud en los registros de la RSNC.
A partir del modelo de pareja doble, la componente (P o a S) de la energía sísmica se
encuentra dado por (Boatwright & Fletcher, 1984)
Para este cálculo integramos sobre el modelo que mejor se ajusta a cada uno de los
espectros, asegurándonos de extrapolar el modelo a altas frecuencias para asegurar
una cobertura completa sobre una banda ancha de frecuencia. De esta forma evitamos
el sesgo en la estimación de la energía, resaltado por Ide & Beronza (2001). Este tipo
de integración asegura que se tengan en cuenta las regiones del espectro por encima
de la frecuencia de esquina, en las cuales se concentra la mayor parte de la energía.
Con los estimados de la energía a nuestra disposición, calculamos ahora la caída de
estrés aparente en cada uno de los terremotos detectados. Para cada
evento común a todas las estaciones, calculamos luego la mediana de estos estimados,
a partir de la cual deducimos el comportamiento global del estrés aparente.
La figura 6.8 muestra vs . Una relación lineal se observa directamente en los
datos (puntos negros) demostrando la no auto-similaridad en los terremotos presentes
en el nido de Bucaramanga. Realizamos una regresión lineal de mínimos cuadrados
(rojo) para determinar la magnitud del escalamiento. Los intervalos de confianza al
95% de la regresión se muestran a los lados de esta línea (morado).
A partir de la pendiente de la regresión se puede obtener una relación de escalamiento
, lo cual contradice de nuevo la hipótesis de auto-similaridad que implica
constante. Al mismo tiempo, este resultado verifica la magnitud del escalamiento
64
obtenido en , siendo el porcentaje de diferencia de estos del 1.75%. Nótese la clara
semejanza de estos resultados, aún cuando las medidas de y se realizaron
mediante métodos completamente independientes12.
Figura 6.8. Relación entre el estrés aparente y el momento sísmico para 51 eventos comunes a
todas las estaciones. La línea roja corresponde a una interpolación lineal de los datos, a partir de la cual
se puede obtener la pendiente de escalamiento. La línea morada establece los intervalos de confianza al
95% de la interpolación. La estimación de las barras de error para resulta muy difícil aquí.
Existen varias fuentes de error en este cálculo, la más importante de ellas es el cálculo
de a partir de los estimados de en el catálogo. Otras fuentes de error pueden
provenir de la incertidumbre en la estimación de las frecuencias de esquina y las
velocidades de onda y supuestas para la fuente. Aún cuando las velocidades son
tomadas de estudios tomográficos en la región (Ojeda & Havskov, 2001), su
dependencia con la quinta potencia en la energía conduce a errores notorios,
aún para incertidumbres pequeñas en .
Para hacer más completo nuestro estudio sobre la energía de las ondas sísmicas, vale
la pena también calcular la partición de energía sísmica
La cual determina cuanta energía adicional tiene la onda S por encima de la onda P, y
sirve para estudiar qué fracción de la energía total se le asigna a cada fase. Boatwright
& Fletcher (1984) realizaron el estudio teórico de esta expresión y llegaron a
12 De hecho las figuras de estrés aparente y de caída de estrés en función de momento sísmico son tan parecidas, que con el escalamiento vertical adecuado, podrían fácilmente superponerse. Lo cual indica una relación lineal entre estas variables, la cual hemos calculado como .
65
Este valor nos sirve a modo de parámetro de calibración para probar la auto-
consistencia de nuestro método, ya que si alimentamos la fórmula anterior con el valor
promedio obtenido para deberíamos ser capaces de reproducir ,
donde y se calculan por integración directa del modelo.
La figura 6.9 muestra los valores de obtenidos para cada evento a partir del cálculo
de las energías y . La línea roja corresponde al valor teórico, suponiendo
, tal y como se obtuvo anteriormente.
Fig. 6.9. Razón entre las energías S y P obtenidas. Los resultados demuestran la auto-congruencia del
método utilizado, validado los cálculos de y como equivalentes.
A pesar de la dispersión de los datos en la figura 6.9, la mediana de estos puntos
concuerda con el valor esperado según las relaciones teóricas. Obtuvimos
, lo cual concuerda muy bien con el valor de calculado a partir de
las frecuencias de esquina . Estos resultados corresponden a un error
relativo13 alrededor del 1% y a una desviación estándar de 2.87. Comprobando así la
auto-consistencia del método.
7. Conclusiones
Con el fin de determinar las propiedades de escalamiento sísmico en el nido de
Bucaramanga, estudiamos 80 eventos provenientes del nido, los cuales fueron
detectados durante los años 2007- 2010 por la Red Sismológica Nacional de Colombia
(RSNC) en las estaciones banda ancha BRR, CHI, PRA, RUS y HEL.
13
.
66
Calculamos el espectro de las ondas P y S presentes en cada uno de los terremotos,
corregimos el espectro y limitamos la banda de frecuencia de los datos con el fin de
obtener el espectro de la fuente sísmica.
Luego, mediante la comparación con un modelo de fuente teórico calculamos la
frecuencia de esquina asociada a cada evento. Paralelamente, calculamos también la
energía sísmica liberada en cada terremoto mediante la integración del espectro de
velocidad sobre toda la banda de frecuencia.
Con los resultados de la frecuencia de esquina en cada una de las estaciones,
limitamos nuestra atención únicamente a los eventos comunes a las estaciones
estudiadas (lo cual elimina los eventos más pequeños). Calculamos la mediana de las
frecuencias de esquina de cada evento, y a partir de las frecuencias de esquina
resultante, calculamos la caída del estrés y la dimensión de la fuente . Por otra
parte, el cálculo del estrés aparente asociado a cada evento se realizó tomando la
mediana de los valores calculados para cada una de las estaciones. De esta forma nos
aseguramos de obtener un estimado confiable para .
Observamos un escalamiento en las frecuencias de esquina producto de una relación
no auto-similar con respecto al momento sísmico, la desviación del comportamiento
auto-similaridad fue de para las ondas P y para
las ondas S. Los intervalos de confianza al 95% para los estimados son
y dependiendo del componente escogido (P
o S).
Los resultados obtenidos para parecen distanciarse un poco del modelo
teórico utilizado para su modelación, lo cual puede indicar diferencias en la geometría
de la fuente sísmica, variaciones en su velocidad de ruptura, o sencillamente un
cubrimiento acimutal insuficiente de estaciones. La razón de las frecuencias de
esquina P y S es según Madariaga (1976) , el valor obtenido en este
estudio es de , incurriendo en un porcentaje de diferencia del 11%
y una desviación estándar de 0.46.
La partición de energía entre las fases P y S ( ) resulta en valores particulares
de mucho más dispersos, pero con una tendencia precisa a replicar el resultado
esperado según las frecuencias de esquina calculadas. Calculando el valor de
mediante energías y mediante frecuencias de esquina obtenemos valores de medias
muy aproximadas (1% de diferencia entre ellos), con una desviación estándar de 2.87.
Lo cual comprueba la validez del método utilizado y demuestra la equivalencia entre
los indicadores de escalamiento y .
Los resultados de la variación de y en función del momento sísmico
muestran dependencias de la forma y , lo cual se distancia de
la dependencia funcional esperada en el caso de cumplirse la hipótesis
67
de auto-similaridad. Estos resultados nos llevan a refutar el comportamiento auto-
similar para los terremotos en el nido de Bucaramanga.
Una de las limitantes importantes de este estudio se relaciona con las bases de datos
utilizadas. Desafortunadamente son pocos los eventos que se pueden extraer de las
bases de datos actuales de la RSNC tal que permitan un estudio detallado y confiable
de la fuente. Algunos de los obstáculos, junto con su respectiva solución implementada
son:
Tabla 7.1 Obstáculos e Implementaciones en los datos de la RSNC
Obstáculo Implementación Cantidad limitada de estaciones banda ancha lo suficientemente cerca para detectar terremotos pequeños y profundos (algunas estaciones con pocos eventos son: POP2, FLO2, TUM, MON)
Ampliar la cantidad de datos al máximo incluyendo todos los eventos posibles
Fuertes efectos de sitio encontrados en algunas estaciones. (ej: BRR y CHI)
Limitar la banda de frecuencia dependiendo de los efectos observados
Presencia predominante del ruido en muchos de los registros. Esto es, baja razón señal/ruido para muchos de los eventos (en especial aquellos con magnitudes
).
Calculo del error y recorte de las regiones por debajo de la línea de confianza al 90%. Limitarse a eventos con
Una tasa de muestreo insuficiente para muchas de las estaciones. Lo cual se traduce en una banda de frecuencia limitada. (ej: ROSC)
Evitar estaciones con tasas de muestreo menores a 100 Hz
Historia limitada de muchas de estas estaciones, pues presentan eventos viables solamente después del 2008/2009. Ya sea porque la estación es nueva o porque sus eventos no son viables.
Buscar incluir la mayor cantidad de estaciones posible
Eventos sin llegadas P y S demarcadas Omitir tales eventos
Espectros observados en una misma estación (ej: CHI, RUS) con densidades espectrales de amplitud mucho menores que eventos de igual magnitud. Producto de re-calibraciones o cambios de instrumentación en la estación.
Omitir aquellos eventos problemáticos y dejar solo los eventos que concuerden con los valores esperados según las tablas de instrumentación actuales en la RSNC.
Fuimos especialmente cuidadosos de tomar todas las medidas necesarias para evitar o
combatir cada uno de estos limitantes, con el fin de basarnos únicamente en los
eventos más confiables para la estimación de los parámetros de fuente.
En particular, tuvimos que limitar la banda de frecuencia de los terremotos hasta un
máximo de 20 Hz, cortar las regiones a bajas frecuencias con efectos notorios de error
y limitar el análisis a los terremotos con magnitudes locales mayores o iguales a 3
(debido a una baja SNR).
68
8. Tabla de resultados específicos
Por completitud, incluyo una tabla con los parámetros de fuente obtenidos para cada
uno de los eventos estudiados. Las unidades utilizadas son: Hz para las frecuencias, Pa
para la caída de estrés y el estrés aparente y m para el radio.
frec. P frec. S Estrés Ap_Estres radio Es/Ep Magnitud
7.21 4.63 4.01E+06 4.69E+04 201.75 7.00 3.11
8.42 2.89 1.23E+08 2.44E+06 248.63 1.08 4.49
12.71 8.21 1.12E+07 1.30E+05 114.05 7.16 2.91
10.12 5.14 4.37E+07 5.75E+05 162.78 3.49 3.70
5.94 5.31 3.03E+06 3.27E+04 209.86 18.74 2.91
16.03 9.23 1.59E+07 1.93E+05 95.99 5.10 2.91
10.69 5.22 2.57E+08 3.48E+06 157.60 3.10 4.20
8.77 5.47 5.24E+07 6.20E+05 168.33 6.42 3.70
13.77 7.93 1.42E+07 1.73E+05 111.75 5.08 3.01
3.44 2.16 2.90E+08 3.45E+06 427.00 6.55 4.99
10.25 6.69 6.06E+06 7.01E+04 140.71 7.36 2.91
12.71 5.91 1.67E+07 2.33E+05 136.37 2.69 3.31
8.68 5.95 1.20E+07 1.37E+05 162.16 8.49 3.21
9.81 5.47 6.60E+06 8.21E+04 159.67 4.59 3.11
7.22 4.00 7.29E+06 9.14E+04 217.60 4.51 3.40
4.83 2.63 2.06E+06 2.62E+04 328.53 4.24 3.40
7.38 5.05 4.12E+07 4.72E+05 190.91 8.42 3.70
6.47 2.46 3.78E+07 6.58E+05 303.22 1.45 4.30
16.17 6.16 3.75E+07 6.40E+05 121.13 1.49 3.50
8.15 5.15 4.38E+07 5.15E+05 179.89 6.67 3.70
11.50 6.40 1.06E+07 1.31E+05 136.34 4.58 3.11
7.64 4.88 6.64E+06 7.78E+04 190.80 6.90 3.21
7.41 4.65 1.15E+07 1.35E+05 198.48 6.54 3.40
7.84 4.91 4.79E+06 5.65E+04 187.85 6.51 3.11
9.82 6.66 5.98E+06 6.83E+04 144.10 8.24 2.91
12.76 6.86 1.30E+07 1.65E+05 125.20 4.13 3.11
13.50 7.57 8.77E+06 1.08E+05 115.64 4.69 2.91
18.45 8.57 1.80E+07 2.50E+05 94.05 2.70 3.01
14.17 7.82 9.69E+06 1.20E+05 111.09 4.49 2.91
8.72 4.68 2.07E+06 2.64E+04 183.39 4.08 2.91
9.69 6.51 7.88E+06 9.04E+04 146.71 8.01 3.01
8.80 6.51 5.58E+06 6.22E+04 154.25 10.68 2.91
8.12 4.38 9.59E+06 1.22E+05 196.20 4.16 3.40
9.40 4.19 2.11E+06 3.08E+04 189.05 2.36 3.01
13.31 7.58 1.25E+07 1.52E+05 116.29 4.92 3.01
10.06 5.25 1.17E+07 1.51E+05 161.40 3.78 3.31
7.16 4.62 1.12E+07 1.31E+05 202.65 7.10 3.40
10.35 4.70 4.21E+06 6.03E+04 169.83 2.50 3.11
8.73 4.30 1.28E+07 1.72E+05 192.23 3.17 3.50
5.00 3.77 1.53E+06 1.71E+04 269.03 11.23 3.01
69
12.83 6.50 7.87E+06 1.03E+05 128.61 3.47 3.01
6.19 3.08 1.87E+07 2.52E+05 269.36 3.27 3.90
11.75 5.37 2.49E+07 3.55E+05 149.05 2.55 3.50
9.31 5.09 2.68E+06 3.37E+04 169.86 4.34 2.91
8.28 5.18 2.81E+06 3.31E+04 178.03 6.46 2.91
12.04 4.41 3.45E+06 6.24E+04 166.86 1.31 3.11
5.95 3.98 2.86E+07 3.31E+05 239.24 7.89 3.80
11.36 6.50 7.87E+06 9.63E+04 135.85 4.98 3.01
9.81 5.41 2.54E+07 3.18E+05 160.61 4.44 3.50
11.86 6.91 6.67E+06 8.08E+04 128.95 5.24 2.91
10.64 4.91 4.79E+06 6.78E+04 163.62 2.62 3.11
5.25 2.89 5.51E+06 6.95E+04 299.88 4.41 3.60
13.80 7.86 1.39E+07 1.70E+05 112.19 4.92 3.01
12.14 6.17 5.32E+07 6.99E+05 135.79 3.49 3.60
6.77 2.92 3.17E+07 4.80E+05 268.00 2.12 4.10
10.04 5.62 2.86E+07 3.55E+05 155.57 4.66 3.50
9.07 4.48 2.88E+07 3.87E+05 184.76 3.19 3.70
9.71 4.06 5.40E+06 8.38E+04 190.24 1.95 3.31
13.60 7.58 1.76E+07 2.17E+05 115.15 4.61 3.11
9.31 2.69 4.45E+06 1.17E+05 252.32 0.65 3.60
14.25 4.67 5.81E+06 1.22E+05 151.51 0.95 3.21
13.43 4.70 5.92E+06 1.13E+05 153.94 1.15 3.21
11.47 4.09 1.38E+06 2.59E+04 178.12 1.21 2.91
8.12 3.36 8.59E+06 1.35E+05 229.17 1.87 3.60
9.25 3.70 2.05E+06 3.34E+04 205.22 1.71 3.11
11.24 4.67 2.92E+06 4.55E+04 164.94 1.92 3.01
10.24 5.07 2.64E+06 3.54E+04 163.22 3.23 2.91
7.84 3.52 8.82E+05 1.28E+04 225.85 2.40 2.91
6.12 2.44 1.87E+07 3.06E+05 310.56 1.69 4.10
9.33 3.45 1.32E+07 2.36E+05 213.89 1.35 3.70
6.73 3.45 5.24E+07 6.89E+05 243.72 3.56 4.10
11.85 4.69 4.67E+07 7.68E+05 161.33 1.65 3.80
6.50 3.36 8.63E+06 1.13E+05 251.05 3.66 3.60
4.86 2.44 9.27E+06 1.24E+05 341.63 3.33 3.90
11.82 5.57 6.97E+06 9.68E+04 145.62 2.79 3.11
6.89 3.10 8.53E+05 1.24E+04 256.46 2.41 3.01
14.95 4.94 2.44E+06 5.07E+04 143.57 0.97 2.91
11.96 4.94 3.44E+06 5.40E+04 155.69 1.88 3.01
4.94 1.81 5.33E+06 9.74E+04 406.94 1.29 4.00
11.19 4.04 5.30E+06 9.78E+04 181.15 1.25 3.31
70
Agradecimientos
Este trabajo es el resultado de alrededor de dos años de constante aprendizaje e
investigación. Por supuesto, nada de esto hubiese sido posible sin el constante apoyo
de Germán Prieto, quien ha sido mi guía e inspiración en los últimos años. Le
agradezco a él por confiar en mí y por darme las herramientas que necesito para seguir
adelante.
Mi familia ha sido también un apoyo moral enorme, ya que son quienes siempre me
han motivado para continuar adelante en mis proyectos y en mi vida; ellos
indirectamente, fueron de vital importancia para este trabajo.
Ningún análisis habría podido realizarse sin los datos de la Red Sismológica Nacional de
Colombia, los cuales nos cedieron amablemente para el desarrollo de esta y futuras
investigaciones sismológicas en Colombia.
Quiero agradecer finalmente a la universidad de los Andes por financiar durante un
año completo este proyecto. Esta ayuda fue muy importante, ya que me permitió
concentrar mis esfuerzos en el estudio completo de la sismología teórica y
experimental, dándome al mismo tiempo la capacidad de tener un buen rendimiento
en mis obligaciones académicas.
Siempre estaré orgulloso de mi universidad y feliz de verla progresar. Veo con mucha
expectativa el desarrollo del departamento de Geociencias en los Andes, el cual
permitirá la presencia de grandes científicos al servicio del desarrollo del país, así como
una nueva generación de profesionales ampliamente capacitados en las ciencias de la
tierra.
71
Referencias
Abercrombie, R. E. (1995), Earthquake source scaling relationships from -1 to 5 Ml
using seismogram recorded at 2.5-km depth, J. Geophys. Res., 100, 24, 036.
Aki, K., and P.Richards (1980), Quantitative Seismology, 932pp., W. H. Freeman, New
York.
Baltay, A., G. Prieto, y G. C. Beroza (2010), Radiated seismic energy from coda
measurements and no scaling in apparent stress with seismic moment, J. Geophys.
Res., 115, B08314.
Boatwright J. A. (1980), A spectral theory for circular seismic sources: Simple estimates
of source dimesion, dynamic stress drop and radiated seismic energy, Bull. Seis. Soc.
Am., 70, 1-27.
Byerlee J. D. (1978), Friction of rocks, Pure Appl. Geophys., 116, 615-626.
Boatwright, J. A., y J. B. Fletcher (1984), The partition of radiated energy between P and
S waves, Bull. Seismol. Soc. Am., 74, 361-373.
Brune J. N. (1970), Tectonic stress and seismic shear stress from earthquakes, J.
Geophys. Res., 75, 4997-5009.
Chapman C. (2004), Fundamentals of seismic wave propagation, Cambridge Uni. Press,
Cambridge UK.
Cortés M. y Angelier C. (2005), Current states of stress in the northern Andes as
indicated by the focal mechanisms of earthquakes, Tectonophysics, 29-58.
Dahlen F. A. (1974), On the ratio of P and S-waves corner frequencies for shallow
earthquake sources, Bull. Seis. Soc. Am., 64, 1159-1180.
Elsheby, J. D. (1957), The determination of the elastic field of an ellipsoidal inclusion
and related problems, Proc. R. Soc. London A, 241, 376-396.
Eysteinn T. y Lawson J. E. (1970), The intermediate earthquake source near
Bucaramanga, Colombia, Bull. Seis. Soc. Am., 60, 269-273.
Frohlich et. al (1975), A reexamination of the Bucaramanga, Colombia, earthquake
nest, Bull. Seis. Soc. Am., 85, 1622-1634.
Gök R. y Hutchings L. (2009), Source Parameters for 1999 North Anatolian Fault Zone
Aftershocks, Pure and Applied Geophysics.
Hirasawa T. y Stauder W. (1973), On the seismic waves from a finite moving source,
Bull. Seis. Soc. Am., 55, 1811-1842.
72
Hough, S. E. (1997), Empirical Green's function anlaisis: Taking the next step, J. of
Geophys. Res., 102, 5369-5384.
Ide S., y G. C. Beronza (2001), Does apparent stress vary with earthquake size?,
Geophis. Res. Lett., 28, 3349-3352.
Ide S., y G. C. Beronza, S. G. Prejean, and W. L. Ellworth (2003), Apparent break in
eathquake-scaling due to path and site effects on deep bore hole recordings, J.
Geophys. Res., 108(B5), 2271.
Jeffreys H., Jeffreys B. S., Methods of Mathematical Physics, 3 ed, Cambridge Univ.
Press, Cambridge UK.
Kanamori, H. (1977), The energy realese in great earthquakes, J. Geophys. Res., 82,
2981-2987.
Kanamori, H., J. Mori, E. Hauksson, T. H. Heaton, L. K. Hutton, and L. M. Jones (1993),
Detemination of earthquake energy realese and ML uning TERRAscope, Bull. Sismmol.
Soc. Am., 83, 330-346.
Kostrov B.V., Das S. (1998), Principles of earthquake source mechanics,Cambridge Univ.
Press, Cambridge UK.
Madariaga, R. (1976), Dynamics of an expaning circular crack, Bull. Seimolog. Soc. Am.,
66, 639-666.
Malagnini (2008), Strong evidence of non-similar earthquake source scaling in central
Italy, Geophysical Research Letters, Vol 35.
Malagnini, L., S. Nielsen, K. Mayeda, y E. Boschi (2010), Energy radiation from
intermediate to large magnitude earthquakes: Implications for dynamic fault
weakening, J. Geophys. Res., 115, B06319, doi:10.1029/2009JB006786.
Mayeda, K., y W. R. Walter (1996), Moment, Energy, Stress drop and source spectra of
Western United States earthquakes from regional coda envelopes, J. Geophys. Res.,
101, 11, 195-11, 208.
McNamara D. E, y Buland R. P., Ambient Noise Levels in the Continental United States,
Bull. Seis. Soc. of Am., 94, 4, 1517-1527.
Molnar P., Tucker B. E., Brune J. N. (1973), Corner frequencies of P and S waves and
models of earthquake source, Bull. Seis. Soc. of Am., 63, 2091-2104.
Ojeda A., Havskov (2001), Crustal Structure and Local Seismicity in Colombia, Journal of
Seismology, Kluwer Academic Publishers.
73
Park, J., C. R., y F. L. Vernon (1987), Multitaper spectral analysis of high frequency
seismograms, J. Geoph. Res., 92, 12,675-12,648.
Pennington W. (1981), Subduction of the eastern panama basin and seismotectonics of
northwestern South America, Jour. Gephys. Res., 86, 10753-10770.
Prieto G. A., P. M. Shearer, F. L. Vernon, y D. B. Kilb (2004), Earthquake source scaling
and self-similarity estimation from stacking P and S spectra, J. Geophys. Res., 109,
B08310.
Prieto, G. A., D. J. Thomson, F. L. Vernon, P. M. Shearer y R. L. Parker (2007),
Confidence intervals of earthquake source parameters, Geophys. J. Int., 168, 1227-
1234.
Savage J. C. (1966), Radiation from a realistic model of faulting, Bull. Seis. Soc. Am., 56,
577-592.
Shearer P. (1999), Introduction to seismology, Cambridge Univ. Press, Cambridge UK.
Taboada A., Dimate C., Fuenza A. (1998), Sismotectónica de Colombia: deformación
continental activa y subducción, Física de la Tierra, 10, 111-147.
Thorne L. (1995), Modern Global Seismology, Academic Press, San Diego CA.
Tomic J. et. al (2009), Source parameters and rupture velocity of small m<=2.1 resevior
induced earthquakes, Geophy. J. Int., 179, 1013-1023.
Van der Hilst R., y Mann, P. (1994), Tectonic implications of tomographic images of
subducted lithosphere beneath northwestern South America, Geo. Soc. Am., 22, 451-
454.
Yamada T., J. J. Mori, S. Ide, R. E. Abercrombie, H. Kawata, M. Nakatami, Y. Iio, and H.
Ogasawara (2007), Stress drops and radiated seimic energies of microearthquakes in a
South Sfrican gold mine, J. Geophys. Res., 112, B03305.
Walter, W. R., K. M. Mayeda, R. Gok, y A. Hofstetter (2006), The scaling of seismic
energy with moment: simple models compared with observations, in Earthquakes:
Radiated energy and the physics of fauling, Geophys. Mon. Series., 170, 22-41.
Zarifi Z., Havskov J., Hanyga A. (2006), An insight into the Bucaramanga nest,
Tectonophysics, 443, 93-105.
Top Related