MODELO MÍNIMO UNIDIMENSIONAL DE VELOCIDADES DE LA ONDA P PARA LA CORDILLERA VOLCÁNICA DE GUANACASTE, COSTA RICA

MINIMUM 1D P WAVE VELOCITY MODEL FOR THE GUANACASTE VOLCANIC ARC, COSTA RICA

María C. Araya1, 2*, Lepolt Linkimer1, 2 & Waldo Taylor1, 3

1Red Sismológica Nacional (RSN: UCR-ICE), Apdo. 214-2060, San Pedro, Costa Rica

2Escuela Centroamericana de Geología, Universidad de Costa Rica

3Observatorio Sismológico Vulcanológico Arenal-Miravalles (OSIVAM), Instituto Costarricense de Electricidad, 0032-1000, San José, Costa Rica.

*Autora para contacto: mariacristina.araya@ucr.ac.cr

(Recibido: 25/09/2015; aceptado: 25/01/2016)

ABSTRACT: In this study a minimum 1D model for the Guanacaste Volcanic Arc (CVG) is derived from 475 local earthquakes registered between January 2006 and July 2014 by the Arenal-Miravalles Seismological and Volcanic Observatory (OSIVAM). This velocity model has six layers and covers from the surface to a depth of 80 km. The model has velocities from 3.96 to 7.79 km/s from top to bottom. The station corrections obtained range between -0.28 and 0.45. These corrections show a trend with positive values located over the volcanic arc and negative values located on the forearc, in agreement with the crustal thickness. Three main seismicity clusters were acknowledged from the earthquake relocations, which may be associated with activity on inferred faults. The obtained one dimensional velocity model displays a simplified version for the crustal structure and aim to improve the earthquake location processes performed by the OSIVAM.

Keywords: Minimum 1D velocity model, seismicity of Costa Rica, Guanacaste volcanic arc, crust structure, seismic network.

RESUMEN: En este estudio se calcula un modelo mínimo de velocidades para la Cordillera Volcánica de Guanacaste a partir de 475 sismos locales registrados por el Observatorio Vulcanológico y Sismológico Arenal Miravalles (OSIVAM) entre enero del 2006 y julio del 2014. El modelo consta de seis capas desde la superficie hasta 80 km de profundidad y posee velocidades que varían entre 3,96 y 7,79 km/s. Las correcciones para las estaciones obtenidas varían entre -0,28 a 0,45 y muestran una tendencia de valores positivos sobre el arco volcánico y negativos en el antearco, en concordancia con el espesor de la corteza. La relocalización de los sismos muestra tres grupos principales de epicentros que podrían asociarse con actividad en fallas inferidas. El MMV1D presentado en esta investigación brinda una idea simplificada de la estructura de la corteza y pretende contribuir con el mejoramiento de la localización rutinaria de sismos que realiza el OSIVAM.

Palabras clave: Modelo mínimo de velocidad, sismicidad de Costa Rica, Cordillera Volcánica de Guanacaste, Estructura de la corteza, Redes sísmicas.

INTRODUCCIÓN

Un adecuado modelo mínimo de velocidad en una dimensión (MMV1D) es un requisito para localizar sismos con buena calidad. Estos modelos son una representación simplificada de la velocidad de las ondas sísmicas al atravesar capas geológicas y son esenciales para el estudio de la estructura de la corteza y la sismicidad de una zona en particular.

Al noreste de Costa Rica se ubica la Cordillera Volcánica de Guanacaste (CVG), donde el Instituto Costarricense de Electricidad (ICE), tiene en operación proyectos de energía geotérmica, hídrica y eólica (Fig. 1). El instituto además creó el Observatorio Vulcanológico y Sismológico Arenal Miravalles (OSIVAM), para realizar estudios sismológicos y vulcanológicos en la región de influencia de estos proyectos (Barquero, 1990).

El OSIVAM inició sus funciones en 1974 con únicamente ocho estaciones sísmicas. A partir del 2004 el número de estaciones se incrementó hasta un total de 42, distribuidas en la CVG, lo cual significó el mejoramiento de la detección de sismicidad en el área. No obstante, el OSIVAM ha utilizado hasta el año 2014 un modelo de velocidades derivado a partir de datos obtenidos antes de la expansión de la red (i.e., Matumoto et al., 1977 & Taylor, 2000) y estudios internos del observatorio que no han sido publicados.

Debido a la necesidad de contar con un MMV1D adecuado para la configuración de la red de instrumentos del OSIVAM, en este estudio se calcula un MMV1D derivado a partir de sismos locales registrados por el observatorio entre enero del 2006 y julio del 2014.

El MMV1D presentado brinda una idea simplificada de la estructura de la corteza de la CVG, pretende contribuir con el mejoramiento de la localización rutinaria de sismos de OSIVAM, permite delimitar con mejor precisión las fuentes sísmicas de la CVG y sirve de base para generar tomografías sísmicas.

MARCO TECTÓNICO

Centroamérica se encuentra situada sobre el margen occidental de la placa Caribe; esta es una placa oceánica engrosada originada durante la etapa inicial del punto caliente de Galápagos, en el Cretácico Tardío (i.e. Galli-Olivier, 1979; Hill, 1993; Von Huene et al., 2000; Barkhausen et al., 2001; Sadofsky et al., 2009). El espesor de la placa Caribe en el norte de Costa Rica es de aproximadamente 40 km (Sallarès et al., 2001). El basamento de esta placa presenta velocidades sísmicas altas y típicas de rocas máficas intrusivas (Sallarès et al., 1999) y comprende el llamado “Complejo de Nicoya” que incluye: gabros, basaltos, basaltos en almohadilla y sedimentos volcánicos y marinos (i.e. Galli-Oliver, 1979; Buffler, 1982; Flores et al., 2003., Denyer et al., 2014).

Por encima del basamento, las velocidades sísmicas y densidades de la roca son las esperadas en una corteza oceánica (Sallarès et al., 1999). En el nivel superior la placa es altamente heterogénea e incluye una secuencia de rocas sedimentarias con edades de entre el Cretácico Tardío hasta el reciente y sedimentos volcánicos porosos de vulcanismo reciente y con velocidades del orden de entre de 2 km/s y 3,3 km/s (i.e. Buffler, 1982; Sallarès et al., 1999, Sallarès et al., 2001).

A lo largo del margen oeste de la placa Caribe en la fosa Mesoamericana se subduce la placa del Coco a una velocidad entre 70 y 90 mm por año y con una edad entre 22 a 23 m.a. (Barckhausen et al., 2001, DeMets, 2001). Como consecuencia de la subducción, entre el Plioceno y Pleistoceno, se originó arco volcánico actual, paralelo a la trinchera, a una distancia media de 150 km de la fosa Mesoamericana. Este arco en la CVG consiste de noroeste a sureste de los siguientes cinco estratovolcanes activos: Orosí, Rincón de la Vieja, Miravalles, Tenorio y Arenal (Chiesa et al., 1994).

Entre la CVG y la fosa Mesoamericana existe el Bloque Antearco Centroamericano el cual se desplaza hacia el noroeste con una velocidad de entre 7 y 14 mm/año (i.e., McCaffrey , 1996; DeMets, 2001; Norabuena et al., 2004 & LaFemina et al., 2009), como resultado de la subducción oblicua de la placa del Coco y el fuerte acoplamiento entre las placas (i.e., McCaffrey, 1996; LaFemina et al., 2009).

La sismicidad en la CVG tiene dos fuentes principales: el proceso de subducción en el que se general sismos de profundidad intermedia (35 a 70 km) y las fallas de desplazamiento de rumbo con orientaciones predominantes norte-sur, noreste-suroeste y noroeste-sureste (i.e. Montero & Alvarado, 1988; Barquero, 1990) que generan sismos superficiales (0 a 15 km). Algunas de las fallas presentes en la zona de estudio son: Caño Negro, Rincón de la Vieja, Liberia, Guayabo, Bijagua, Cote, Chiripa y Cañas. Históricamente, en el área, han ocurrido siete terremotos de origen cortical en 1853 (Cañas, Mw 6,0), 1911 (Guatuso, Mw 6,5), 1935 (Bagaces, Mw 6,2), 1941 (Bagaces, Mw 6,3), 1973 (Tilarán, Mw 6,5), 2002 (Bijagua, Mw 5,4) y 2011 (Armenias, Mw 5,5).

MODELOS PREVIOS EN LA CVG

El primer modelo de velocidad de ondas sísmicas en la corteza de Costa Rica fue realizado por Matumoto et al. (1977). Este modelo se obtuvo con datos provenientes de seis estaciones cercanas al volcán Arenal y dos en la península de Nicoya. Modelos de velocidades subsecuentes para la zona norte de Costa Rica incluyen los realizados por Liaw (1981), Buffler (1982), Protti et al. (1996), Colombo et al. (1997), Sallarès et al. (1999), Sallarès et al. (2000), Taylor (2000), Sallarès et al. (2001), Quintero & Kissling (2001), Husen et al. (2003), DeShon & Shwartz (2004), DeShon et al. (2006), MacKenzie et al. (2008), Syracuse et al. (2008), Arroyo et al. (2009), Zhu et al. (2009), Harmon et al. (2013) y Arroyo et al. (2014). El detalle de algunos modelos más relevantes para este estudio se resume en el cuadro 1.

Otras investigaciones han aportado información adicional para definir la estructura de la corteza. Por ejemplo, MacKenzie et al. (2008) y Linkimer et al. (2010) por medio del análisis de funciones de receptor, determinaron la profundidad del Moho para la placa Caribe como de 37.9± 5.2 km y entre 36 y 42 km respectivamente y calcularon una razón Vp/Vs de entre 1,67 a 1,78 para Guanacaste. Además, Lücke (2012) basado en datos gravimétricos satelitales determinó que la profundidad del Moho en la CVG, varía entre 37 y 39 km bajo los volcanes Rincón de la Vieja y Arenal y es de 42 km bajo los volcanes Miravalles y Tenorio.

METODOLOGÍA

Un MMV1D es una representación simplificada de la velocidad a la que viajan las ondas símicas a través de capas geológicas inferidas como horizontales y lateralmente homogéneas bajo un punto de referencia. El MMV1D incluye además correcciones para las estaciones de forma que la heterogeneidad en la estructura es considerada. Para obtener estas velocidades es necesario registrar los sismos y para ello es esencial conocer la ubicación y características de las estaciones sísmicas y el tiempo en que los sismos son registrados por dichas estaciones. Por medio de un modelo de velocidades inicial se puede determinar el tiempo de origen de los sismos y el tiempo del viaje de las ondas sísmicas desde el hipocentro hasta la estación.

En esta investigación se recopilaron las características de las estaciones sísmicas del OSIVAM que registraron sismos durante el periodo 2006-2014, esto es latitud, longitud, elevación, tipo y sensibilidad de los sensores y frecuencia de muestreo de los digitalizadores, con el fin de obtener los datos primarios para el cálculo posterior de la magnitud momento (Mw).

La lectura de los sismos se realizó con el programa MULPLT de SeisAn (Ottëmoller et al., 2011). Un total de 950 sismos calificados como de buena calidad fueron seleccionados como la base de datos inicial. Estos sismos se caracterizaron por un mínimo de ocho arribos de la onda P, una cobertura de estaciones no menor a 200º (gap < 200º) y localizaciones iniciales ubicadas entre las latitudes 10º y 11º y las longitudes -84º y -86. Durante la identificación de los arribos, a cada lectura de la onda primaria (P), se le asignó un grado de incertidumbre (pesos) en cinco niveles de 0 a 4, siendo 0 un grado que representa un arribo muy evidente y por lo tanto fue usado en un 100% en el cálculo de la localización. El resto de pesos asignados del uno al cuatro representan un aumento en 25% de la incertidumbre del arribo, de forma tal que el peso uno es usado en un 75% en el cálculo de la localización, el peso dos en un 50%, el peso tres en un 25% y el peso cuatro no es usado. Además, se tomaron en cuenta únicamente las estaciones ubicadas a menos de 150 km del epicentro (Fig. 2).

Una vez localizados los 950 sismos de la base de datos inicial con el uso de un modelo de velocidad inicial, se procedió siguiendo la metodología de Husen & Hardebeck (2010) para seleccionar los sismos que luego fueron usados en la inversión del MMV1D. Estos sismos seleccionados tuvieron las siguientes características: una cobertura de estaciones alrededor del evento mínima de 180° (Gap <180º), un mínimo de ocho arribos de la onda P y profundidad mínima de 0,5 km. En total 475 sismos cumplieron con estos requisitos lo que significó un total de 6491 arribos de la onda P en 79 estaciones. En promedio, estos sismos presentaron un RMS de 0,15 y un gap de 173º (Fig. 3). La estación de referencia para el cálculo del MMV1D fue CUI (Cuipilapa, Fig. 3).

Se utilizó el programa VELEST para realizar las inversiones simultáneas de los hipocentros y las velocidades con los arribos de las ondas P para cada sismo. La inversión se inició iterando posibles modelos a priori provenientes de estudios realizados en la zona (Cuadro 1). Luego de haber realizado las inversiones para cada modelo por separado, se compararon los resultados y se escogió el “mejor modelo a priori” como aquel con el menor RMS. Pequeñas modificaciones en la estructura del “mejor modelo a priori” fueron realizadas para probar la estabilidad de los resultados y evaluar el RMS final. Es importante mencionar que en estos modelos, la altura 0 es considerada como la altura del nivel del mar y las alturas negativas hacen referencia a la topografía por encima del nivel del mar.

Dos tipos de pruebas de sensibilidad fueron realizadas al modelo propuesto para evaluar su estabilidad y dependencia del modelo inicial y de los hipocentros de entrada. En el primer grupo de pruebas se realizó un cambio en la velocidad del modelo inicial, donde se aumentó y disminuyó 0,5 y 1 km/s en cada capa, con el fin de corroborar la estabilidad del MMV1D a los cambios en las velocidades iniciales. El segundo grupo de pruebas se perturbó la ubicación inicial de los hipocentros de forma tal que la ubicación inicial de los sismos fuese alterada aleatoriamente de 0 a 3 km y de 0 a 5 km en latitud, longitud y profundidad. Esta prueba tuvo el fin de comprobar la estabilidad de las localizaciones finales. Ambos grupos de pruebas fueron satisfactorios y los detalles se encuentran en Araya (2014). Por último, se relocalizaron los 950 sismos de la base de datos inicial con el MMV1D propuesto con el fin de identificar las zonas de alta sismicidad y comparar las nuevas localizaciones con respecto a las obtenidas con el modelo inicial.

MODELO MÍNIMO DE VELOCIDAD

El MMV1D fue obtenido con una densa cobertura de rayos entre hipocentros y estaciones (Fig. 4). Los sismos en profundidad se pueden separar en dos grupos: uno entre 0 y 20 km y otro con profundidades mayores a 40 km. En la base de datos usada no hay hipocentros con la calidad requerida entre 20 y 40 km de profundidad, no obstante las trayectorias de sismos con profundidad mayor a 40 km resuelven parcialmente las velocidades de las capas a estas profundidades. La región con el mejor muestreo es la más superficial, ya que aproximadamente el 80% de los hipocentros tienen una profundidad menor a 10 km.

La inversión del MMV1D para los diferentes modelos a priori estudiados (Cuadro 1) mostró un descenso importante del RMS luego de la inversión. En los casos de Liaw (1981) y Sallarès et al. (2001) el descenso del RMS con respecto al RMS inicial fue de 48% y 12%, respectivamente. El modelo con el menor error inicial y el mayor porcentaje de reducción del RMS luego de la inversión (66%) fue el modelo de seis capas de Taylor (2000), cuyo RMS inicial fue de 0,2525 y final fue de 0,1282. La profundidad de cada capa de este MMV1D es: 0, 3,2, 7,2, 16, 41 y 80 km y las velocidades obtenidas son: 3,96, 5,75, 6,22, 6,99, 7,61 y 7,79 km/s, respectivamente (Cuadro 2, Fig. 5).

Al MMV1D obtenido se le realizaron modificaciones en el espesor de las capas 1, 2 y 3 ya que estas capas presentaron un cierto grado de inestabilidad en las pruebas de sensibilidad a los cambios de la velocidad inicial en el modelo a priori. La modificaciones se realizaron a prueba y error de manera sistemática hasta encontrar el menor RMS en la inversión (Araya, 2014). De las pruebas realizadas, el MMV1D que presentó el menor RMS fue de 0,1280, el cual es muy similar al obtenido de la inversión del mejor modelo a priori. La profundidad de cada capa de este MMV1D es: -2, 3, 7,1, 16, 41 y 80 km, y las velocidades obtenidas fueron: 3,96, 5,74, 6,21, 7,02, 7,52 y 7,54 km/s, respectivamente.

De la velocidad en la capa más superficial del MMV1D se puede inferir que está compuesta por rocas poco compactadas ya que las ondas sísmicas viajan lento (menor a 5 km/s). Esto concuerda con la geología del área, la cual esta compuesta por lavas, piroclástos y rocas sedimentarias (Flores et al, 2003). La velocidad de las siguientes capas aumenta conforme incrementa la profundidad debido a que la composición de las rocas se vuelve más máfica.

Como parte del MMV1D se calcularon correcciones para las estaciones de la red del OSIVAM. Estos valores se utilizan para contrarrestar los efectos causados al asumir una geometría de capas uniformes y horizontales. Las correcciones oscilaron entre -0,28 a 0,45 (Cuadro 6-6 en Araya, 2014). Una interpolación de tipo Kriging fue realizada para los valores de las correcciones y presentado en forma de mapa (Fig. 6). Las correcciones en las estaciones muestran una tendencia con valores positivos (de 0 a 0,17) para las estaciones situadas sobre el sector sureste del arco volcánico y negativas (de -0,11 a 0) para las estaciones ubicadas principalmente en el antearco. Estos resultados podrían reflejar que las estaciones del arco volcánico están ubicadas sobre una corteza de mayor espesor que las que se encuentran en el antearco, como lo muestran otros estudios en la zona (i.e., Mackenzie 2008, Linkimer et al., 2010, Lücke, 2012). Las correcciones positivas en el arco también podrían interpretarse como el efecto del paso de las ondas sísmicas a través de cámaras magmáticas, acuíferos y zonas con anomalías geotérmicas.

La localización manual de los 950 sismos permitió distinguir que los hipocentros se concentran en los primeros 10 km del arco volcánico y que durante el periodo de tiempo analizado la sismicidad fue más alta en tres zonas ubicadas al sur volcán Rincón de la Vieja, al suroeste del volcán Miravalles y entre norte del volcán Miravalles hasta el lago Arenal (Fig. 7A).

La comparación de la ubicación epicentral entre los sismos iniciales y la relocalización de los sismos con el MMV1D obtenido muestra que los epicentros se aglomeran en las cercanías de ciertas fallas en las tres zonas (Fig. 7 A y B). La aglomeración de eventos luego de la relocalización ocurre debido a que los tiempos entre la fuente (las fallas) y la estación se ajustan mejor al usar un MMV1D que además, toma en cuenta las correcciones por estación. Esto implica que con la utilización del MMV1D los sismos se encuentran mejor localizados por lo cual aparecen aglomerados en la vecindad de las fallas que los provocó.

La ubicación de la aglomeración de eventos en algunos casos coincide con fallas conocidas o inferidas en trabajos anteriores (i.e. Denyer et al., 2003). Por ejemplo, existe un grupo de epicentros que se ajustan a una falla inferida por Denyer et al (2003) al sur del volcán Rincón de la Vieja, con un rumbo NW-SE (zona 1 en Fig. 7 A y B). Otro grupo de sismos de sismos se encuentra cercano a la falla Guayabo, con dirección NE-SW (zona 2 en Fig. 7 A y B). Otra alineación de sismos con dirección NE-SW se encuentra en la zona de transición entre el arco y el trasarco en las cercanías de la falla Bijagua (zona 3 en Fig. 7 A y B).

CONCLUSIONES

En este estudio fue derivado un MMV1D para la CVG a partir de sismos locales registrados por el OSIVAM entre enero del 2006 y julio del 2014. La base de datos inicial incluyó 950 sismos calificados como de buena calidad. Estos fueron relocalizados y seleccionados bajo criterios de cobertura y calidad para obtener un grupo de 475 sismos que fueron usados para el cálculo del MMV1D.

Para la inversión del MMV1D se escogió como el mejor modelo a priori el presentado por Taylor (2000). El MMV1D obtenido tiene una geometría de 6 capas con profundidad de: -2, 3, 7,2, 16, 41 y 80 km y velocidades de: 3,96, 5,75, 6,22, 6,99, 7,61 y 7,79 km/s, respectivamente.

Las correcciones para las estaciones sísmicas obtenidas en la inversión del MMV1D muestran un patrón con valores positivos sobre parte del arco volcánico y negativos sobre el antearco. Estos resultados se pueden interpretar como el reflejo de una corteza es más espesa en la zona del arco volcánico, en concordancia con otros estudios en la zona. Las correcciones positivas en el arco volcánico podrían también manifestar los efectos del paso de las ondas sísmicas a través de cámaras magmáticas, acuíferos y zonas con anomalías geotérmicas.

La relocalización de los sismos con el MMV1D demostró que los hipocentros están concentrados principalmente en los primeros 10 km del arco volcánico. Además, hizo notar tres zonas con concentraciones de sismos que aparecen mejor definidas en las localizaciones obtenidas con el MMV1D y que podrían asociarse con actividad en una falla ubicada al sur del volcán Rincón de la Vieja y en las fallas Guayabo y Bijagua.

El MMV1D presentado en esta investigación brinda una idea simplificada de la estructura de la corteza en la zona de la CVG. El uso de este modelo puede contribuir con el mejoramiento de la localización rutinaria de sismos que realiza el OSIVAM y por consiguiente en una identificación más precisa de las fuentes sísmicas en la zona cercana a los proyectos del ICE.

AGRADECIMIENTOS

Un agradecimiento especial al personal del OSIVAM y del Área de Auscultación Sísmica y Volcánica del ICE por el aporte de los sismos utilizados para la inversión del MMV1D. Así también, al OSIVAM y la Red Sismológica Nacional (RSN) por el uso de su equipo e infraestructura. Este estudio es parte del proyecto de investigación “Propiedades físicas de la litósfera y sismogénesis en Costa Rica”, código 111-B3-263, de la Escuela Centroamericana de Geología de la Universidad de Costa Rica.

REFERENCIAS

ARAYA, M. C., 2014: Modelo mínimo de velocidades unidimensional para la Cordillera Volcánica de Guanacaste, Costa Rica.- 134 págs. Univ. de Costa Rica, San José, Costa Rica [Tesis de Lic.].

ARROYO, I. G., HUSSEN, S., FLUEH, E. R., GOSSLER, J., KISSLING, E. & ALVARADO, G. E., 2009: Three-dimensional P-wave velocity structure on the shallow part of the Central Costa Rican Pacific margin from local earthquake tomography using off- and onshore networks.- Geophys. J. Int. 179(2): 827-849, DOI: 10.1111/j.1365-246x.2009.04342.x

ARROYO, I. G., GREVENMEYER, I., RANERO, C. R. & von HUENE R., 2014: Interplate seismicity at the CRISP drilling site: The 2002 Mw 6.4 Osa Earthquake at the southeastern end of the Middle America Trench.- Geochem. Geophys. Geosyst. 15: 1-16, DOI: 10.1002/2014GC005359

BARCKHAUSEN, U., RANERO, C. R., VON HUENE, R., CANDE, S. C., & ROESER, H. A., 2001: Revised tectonic boundaries in the Cocos Plate of Costa Rica: Implications for the segmentation of the convergent margin and for Plate Tectonic Models.- J. Geophys.Res. 106(B9), 19207-19220.

BARQUERO, R., 1990: Sismicidad y tectónica de la región noreste de Costa Rica, con énfasis en la zona del Proyecto Geotérmico de Miravalles.- 117 págs. Univ de Costa Rica, San José [Tesis de Lic.].

BUFFLER, R. T., 1982: Geological structure of the forearc region off the west coast of Costa Rica in the Vicinity of the Nicoya Peninsula results of a multifold seismic reflection survey.- 55 págs. Univ. de Texas en Austin, Estados Unidos.

CHIESA, S., ALVARADO, G.E., PECCHIO, M., CORELLA, M. & ZANCHI, A., 1994: Contribution to petrological and stratigraphical understanding of the cordillera de Guanacaste lava flows, Costa Rica.- Rev. Geol. Amér. Cent. 17: 19-43.

COLOMBO, D., CIMINI, G. B. & de FRANCO, R., 1997: Three-dimensional velocity structure of the upper mantle beneath Costa Rica from teleseismic tomography study.- Geophys. J. Int. 131: 189-208.

DEMETS, C., 2001: A new estimate for present-day Cocos-Caribbean plate motion: Implications for slip along the Central American volcanic arc.- Geophys. Res. Lett. 28: 4043-4046.

DENYER, P., AGUILAR, T. & MONTERO, W., 2014: Cartografía geológica de la Península de Nicoya, Costa Rica: Estratigrafía y tectónica.- 207 págs. Editorial UCR, San José.

DENYER, P., MONTERO, W. & ALVARADO, G. E., 2003: Atlas tectónico de Costa Rica.- 79 págs. Editorial de la Universidad de Costa Rica.

DESHON, H. R. & SCHWARTZ, S. Y., 2004: Evidence for serpentinization of the forearc mantel wedge along the Nicoya Peninsula, Costa Rica.- Geophys. Res. Lett. 31: 1-4.

DESHON, H. R., SCHWARTZ, S. Y., NEWMAN, A. V., GONZÁLES, V., PROTTI, M., DORMAN, L. M., DIXON, T. H., SAMPSON, D. E. & FLUEH, E. R., 2006: Seismogenic zone structure beneath the Nicoya Peninsula, Costa Rica, from three-dimensional local earthquake P- and S-wave tomography.- Geophys. J. Int. 164: 109-124, DOI: 10.1111/j.1365-246X.2005.02809.x

FLORES, K., DENYER, P. & AGUILAR, T., 2003: Nueva propuesta estratigráfica: geología de las hojas Matambú y Talolinga, Guanacaste, Costa Rica.- Rev. Geol. Amér. Cent. 28: 131-138.

GALLI-OLIVER, C., 1979: Ophiolite and island-arc volcanism in Costa Rica.- Bull. Geol. Soc. Am. 90(5), 444-452, doi 10.1130/0016-7606

HARMON, N., SALAS, M., RYCHERT, C. A., ABERS, G. & FISHER, K., 2013: Crustal and mantle shear velocity structure of Costa Rica and Nicaragua from ambient noise and teleseismic Rayleigh wave tomography.- Geophys. J. Int. 195: 1300-1313, DOI: 10.1093/gii/ggt309

HILL, R. I., 1993: Mantle plumes and continental tectonics.- Lithos, 30: 193-206.

HUSEN, S., QUINTERO, R., KISSLING, E., & HACKER, B., 2003: Subduction-zone structure and magnetic processes beneath Costa Rica constrained by local earthquake tomography and petrological modelling.- Geophys. J. Int. 155: 11-32.

HUSEN, S. & HARDEBECK, J. L., 2010: Earthquake Location Accuracy.- Community Online Resource for Statistical Seismicity Analysis, DOI: 10.5078/cosrssa-55815573

LAFEMINA, P., DIXON T. H., GOVERS, R., NORABUENA, E., TURNER, H., SABALLOS A., MATTIOLI G., PROTTI, M. & STRAUCH, W., 2009: Fore-arc motion Cocos Ridge collision in Central America.- Geochem. Geophys. Geosyst. 10(5): 1-21, DOI: 10.1029/2008GC002181

LIAW, H. B., 1981: Modeling velocity from an ensemble of earthquakes. - 135 págs. Univ. de Texas en Dallas, Estados Unidos [Tesis de Ph.D.].

LINKIMER, L., BECK, S. L., SCHWARTZ, S. Y., ZANDT, G. & LEVIN, V., 2010: Nature of crustal terrenes and the Moho in northern Costa Rica from receiver function analysis.- Geochem. Geophys. Geosyst. 11(1): 1-24, DOI: 10.1029/2009GC002795

LÜCKE, O. H., 2012: Moho structure of Central America based on three-dimensional lithosphere density modeling of satellite-derived gravity data.- Int. J. Earth Sci. DOI: 20.1007/s00531-012-0787-y

McCAFFREY, R., 1996: Estimates of modern arc-parallel strain rates in fore arcs.- Geology. 24(1): 27-30, DOI: 10.1130/0091-7613

MACKENZIE, L. G., ABERS, G. A., FISHER, K. M., SYRACUSE, E. M., PROTTI, J. M., GONZALES, V. & STRAUCH, W., 2008: Crustal Structure along the southern Central American Volcanic front.- Geochem. Geophys. Geosyst. 9(8): 1-19.

MATUMOTO, T., OHTAKE, M., LATHAM, G. & UMAÑA, J., 1977: Crustal Structure in Southern Central America.- Bull. Seism. Soc. Am. 67(1): 121-134.

NORABUENA, E., DIXON, T. H., SCHWARTZ, S., DESHON, H., NEWMAN, A., PROTTI, M., GONZALEZ, V., DORMAN, L., FLUEH, E.R., LUNDGREN, P., POLLITZ, F. & SAMPSON, D., 2004: Geodetic and seismic constraints on some seismogenic zone processes in Costa Rica.- J. Geophys. Res. 109: 1-25, DOI: 10.1029/2003JB002921

OTTERMÖLER, L., VOSS, P.,& HAVSKOV, J., 2011: SEISAN: the Earthquake Analysis Software for Windows, Solaris, LINUX, and MACOSX, version 9.0.1.- 361 págs. Univ. de Bergen, Noruega.

PROTTI, M., SCHWARTZ, S. Y. & ZANDT, G., 1996: Simultaneous Inversion for Earthquake Location and Velocity Structure Beneath Central Costa Rica.- Bull. Seism. Soc. Am., 86(1): 19-31.

QUINTERO, R. & KISSLING, E., 2001: An improved P-wave velocity reference model for Costa Rica.- Geophys. Int. 40(1): 3-19.

SADOFSKY, S., HOERNLE, K. & DUGGENG, S., 2009: Geochemical variations in the Cocos Plate subducting beneath Central America: implications for the composition of the arc volcanism and the extent of the Galápagos Hotspot influence on the Cocos oceanic crust.- Int. J. Earth Sci. 90: 901-913.

SALLARÈS, V., DAÑOBEITIA, J., FLUEH, E.R. & LEANDRO, G., 1999: Seismic velocity structure across the Middle American land bridge in northern Costa Rica.- J. Geodynamics, 27: 327-344.

SALLARÈS, V., DAÑOBEITIA, J. J. & FLUEH, E. R., 2000: Seismic tomgraphy with local earthquakes in Costa Rica.- Tectonophysics, 329: 61-78.

SALLARÈS, V., DAÑOBEITIA, J. J. & FLUEH, E. R., 2001: Lithospheric structure of the Costa Rican Isthmus: Effects of subduction zone magmatism on an oceanic plateau.- J. Geophys. Res. 106(B1): 621-643.

SYRACUSE, E. M., ABERS, G. A., FISHER, K., MACKENZIE, L., RYCHERT, C., PROTTI, M., GONZÁLES, V. & STRAUCH, W., 2008: Seismic tomography and earthquake locations in the Nicaraguan and Costa Rican upper mantle.- Geochem. Geophys. Geosyst. 9(7):1-22, DOI: 10.1029/2008GC001963

TAYLOR, W., 2000: Shear-wave splitting and improved upper crustal velocity model at the Miravalles Geothermal field, Guanacaste-Costa Rica.- 74 págs. Univ de Bergen, Noruega [Tesis de M.Sc.].

VON HUENE, R., RANERO, C. R., WEINREBE, W. & HINZ, K., 2000: Quaternary Convergent Margin Tectonics of Costa Rica, Segmentation of the Cocos Plate, and Central American Volcanism.- Tectonics, 19(2): 314-334.

ZHU, J., KOPP, H., FLUEH, E. R., KLAESCHEN, D., PAPENBERG, C. & PLANERT, L., 2009: Crustal structure of the central Costa Rica subduction zone: implications for basal erosion from seismic wide-angle data.- Geophys. J. Int. 178: 1112-1131, DOI: 10.1111/j.1365-246X.2009.04208.x

ARAYA, C, LINKIMER, L. & TAYLOR, W., 2016: Modelo mínimo unidimensional de velocidades de la onda P para la cordillera volcánica de Guanacaste.- Rev. Geol. Amér. Central, 54: 179-191, doi: 10.15517/rgac.v54i0.23283