The map of complete Bouguer anomaly of Costa Rica shows an elongated NW-SE trending gravity low in the central region. This gravity low coincides with the geographical region known as the Cordillera Volcánica Central. It is built by geologic and morpho-tectonic units which consist of Quaternary volcanic edifices. For quantitative interpretation of the sources of the anomaly and the characterization of fluid pathways and reservoirs of arc magmatism, a constrained 3D density model of the upper crust was designed by means of forward modeling. The density model is constrained by simplified surface geology, previously published seismic tomography and P-wave velocity models, which stem from wide-angle refraction seismic, as well as results from methods of direct interpretation of the gravity field obtained for this work. The model takes into account the effects and influence of subduction-related Neogene through Quaternary arc magmatism on the upper crust.

1. Introduction and Tectonic Setting

A constrained 3D density model of the upper crust along the Quaternary Central American Volcanic Arc (CAVA) was carried out based on complete Bouguer anomaly data. The main focus of the study was the modeling of the fluid pathways and reservoirs in the upper crust resulting from the magmatic processes associated with the subduction of the Cocos plate beneath the Caribbean plate.

The area of interest was the portion of the CAVA known as the Cordillera Volcánica Central in Costa Rica. Throughout Central America, the volcanic arc shows a segmented disposition along the isthmus marked by gaps in Quaternary volcanism as well as sudden changes in the distance from the Middle American Trench. The geographical and morphological region known as the Cordillera Volcánica Central of Costa Rica is comprised by the Platanar, Barva, Poás, Irazú, and Turrialba volcanic edifices. It is delimited to the NW by the absence of Quaternary stratovolcanoes up to the occurrence of the Arenal-Chato volcanic complex. To the SE the arc is interrupted by a major gap in Quaternary volcanism in the Talamanca region marking the southeastern end of the portion of the CAVA related to the subduction of the Cocos plate (Figure 1).

This portion of the Costa Rican arc is characterized also by unique morphological features such as the high volume of the volcanic edifices relative to the rest of the arc (i.e., Carr et al. [1]). Special interest was put on this portion of the arc because of the presence of an elongated gravity low in the complete Bouguer anomaly map.

Until now, density modeling in the region was restricted mainly to regional 2D interpretations based on inhomogeneous gravity databases (Ponce and Case [2]), also 2D sections along seismic refraction profiles on the northwestern part of the Costa Rican arc (Sallarès et al. [3]; Gödde [4]) and across the volcanic gap in the Talamanca region (Stavenhagen et al. [5]). Density models in the Cordillera Volcánica Central were restricted to the structure of the volcanic edifices (Thorpe et al. [6]) thus accounting only for the effects of masses above the geoid. For this work, the homogenized complete Bouguer anomaly database compiled for the SFB574 was used to model in 3D the crustal structure and the effects of Neogene to Quaternary volcanism on the densities along the arc.

2. Database and Gravity Field Features

On-shore complete Bouguer anomaly maps were generated from the homogenized gravity database of the SFB574. The database includes approximately 20 000 stations and was compiled from previously existent on-shore gravity data from several government, industry, and academic institutions such as GETECH Leeds, BGI (Bureau Gravimétrique International), and ICE (Instituto Costarricense de Electricidad). The complete Bouguer anomaly map shows an arc-parallel gravity low with a minimum of −57 10-5 m/s2 along the Cordillera Volcánica Central (Figure 2).

Previous works (Ponce and Case [2], Montero et al. [8]) show a minimum of approximately −80 10-5 m/s2 for this region. However, upon further review of the gravity database, the height above sea level reported for the single station with a complete Bouguer anomaly value of −75.9 10-5 m/s2 differed in nearly 1000 m with the correspondent value of topographic height above sea level obtained from SRTM topography data. This along with other aberrant values was corrected or taken out because of lack of metadata. The corrected database was then used for the forward modeling of the density structure and mass distribution of the upper crust.

Although the main gravity low coincides in extension and trend with the main volcanic edifices of the Quaternary arc, its axis is shifted approximately 8 km towards the Middle American Trench relative to the main volcanic axis comprised by the Platanar, Poás, Barva, and Irazú volcanoes. The fore-arc region shows a relative Bouguer gravity low along the Aguacate Mountains and on-shore data along the Pacific coast show a peak of approximately 20 10-5 m/s2 on the Herradura promontory. Towards the back-arc region, the Bouguer anomaly values increase to a relative high of 15 10-5 m/s2 in the northeastern Tortuguero plains. The northwestern end of the study area shows a gradual increase in the Bouguer anomaly values along the arc with an inflection at a value of −35 10-5 m/s2 at 20 km NW of the Platanar volcano from which the values plateau until they decrease again towards the Pacific coast in SW Nicaragua. To the SE, the main gravity low shows a strong positive gradient culminating at an alignment between the port towns of Quepos and Limón.

3. Constraints of the 3D Density Model

The data analysis methodological approach emphasized the integration of geological and geophysical constraints into the forward modeling. For geological constraints, a simplified map (Figure 3) outlining the main superficial lithostratigraphic units was summarized from previously existing geological information (Tournon and Alvarado [9]; Denyer and Alvarado [10]) together with the integration of the inferred location of main volcanic vents related to the Neogene Aguacate arc (Alvarado [11]). The geological information was complemented by borehole stratigraphy data from Pizarro [12].

Geophysical constraining focused on direct interpretation of the gravity field (i.e., Euler deconvolution source points and power spectrum analysis) and the inclusion of local earthquake seismic tomography data generated by SFB574 collaborators (Arroyo [13]) and previously published works as well as 2D velocity models based on wide angle seismic refraction surveys (Lizarralde et al. [14]).

3.1. Euler Deconvolution

For the Euler deconvolution solutions, the software REDGER (Pašteka [15]) was used which advantages the calculus by means of regularized derivatives (Pašteka and Richter [16]). The calculation of the Euler source points is based on Euler’s homogeneity equation and results in clusters used to constrain the overall geometry of the model. In this case, a structural index of 0.01 was used to for the approximation of planar structures thereby outlining the main boundaries between bodies of contrasting density (Figure 4). The Euler source points showed clustering between geological units mainly on the southwestern and northeastern boundaries of the Quaternary volcanic arc. The clustering of the source points outlines in depth the heterogeneities in the upper crust caused by the Quaternary volcanism. The location of the source point clusters also correlates well with surface geology in the sense that they coincide with the contacts between the lithostratigraphic units that represent major events in Cenozoic volcanism as well as basement structure towards the back-arc. The distribution of the shallower Euler source points may also outline the upper boundary of heterogeneous bodies in the upper crust.

3.2. Power Spectrum Analysis

Power spectrum analysis was carried out through 2D Fast Fourier Transformation to estimate depths for the structures which cause the measured anomaly. According to this method, for example, Döring [17] and many others, the depth of a corresponding “monopole source point” is obtained from the negative slope of the linear relationship between the logarithmic power spectrum and the wavenumber of the gravity field. Results show two tendencies for the correlation between energy and wavenumber (Figure 5). The local part of the spectrum results in depths of some 11 km. The regional part of the spectrum appoints to greater depths of about 54 km. The source of this/these causing mass/masses is not yet clear. We put the main focus of the modeling on the shallower structures because they are better constrained by other geophysical data such as the seismic tomography. 3D modeling, together with the results of Euler deconvolution, constrained the local modeling in the upper 15 km of the crust. Until now the effect of deeper located structures is not investigated and may account for wider wavelengths of the regional field outside of the area of study and outside the area of interest.

3.3. Correlation with Local Earthquake Seismic Tomography

Further constraints included previously published local earthquake seismic tomography cross sections (Protti et al. [18], Colombo et al. [19], Husen et al. [20]) as well as data from SFB574 collaborators (Arroyo [13]). Based on the slab location relative to the volcanic arc as observed on the seismic tomography cross sections and the shallow results for the source of the local gravity anomaly obtained from power spectrum analysis, it is interpreted that the effects of the distribution of mass directly related to the slab do not have an effect on the gravity field at a local scale. However, low-velocity zones shown by Husen et al. [20] in the upper mantle at a depth of approximately 60 km may have a regional effect on the gravity field. Local low-velocity heterogeneities in the upper crust are present beneath the Quaternary volcanic arc in the Cordillera Volcánica Central as observed in the values of perturbation in from Arroyo [13]. As for the geometry of such heterogeneities, the integration of the data as constraints in the density model shows trenchward-dipping, low-velocity structures, originating from the slab and ascending to the upper crust beneath the Quaternary arc. Remarkable is also the presence of similar structures beneath the Neogene Aguacate arc which hints of remnant effect of paleo-volcanism on the upper crust. With regards to the issue of resolution of the seismic tomography data directly integrated as constraints (Arroyo [13]), the better resolved portions of the Quaternary arc in the upper crust domain are those beneath the Irazú-Turrialba volcanic complex as well as the Poás volcano. This is mainly due to the location of seismologic stations.

4. Results and Discussion

The constrained 3D density model was carried out by the interactive modeling with the software IGMAS (Interactive Gravity and Magnetics Application System, and Schmidt and Götze [21]). This modeling is based on the initial algorithm by Götze [22] and further enhanced by Götze [23], Götze and Lahmeyer [24], Schmidt and Götze [25]. Modeling is carried out interactively by creating cross sections from which constant density bodies are triangulated based on the location of the input vertices. Densities are directly assigned to each body or may be calculated via inversion based on the given geometry.

For this model, 22 modeling cross sections were created trending perpendicular to the volcanic arc with an NE-SW trend. Geometry is based on available constraints as well as fit between the measured and the modeled gravity field. Overall fit between measured and calculated complete Bouguer anomaly is shown on Figure 6.

4.1. Overall Context and Horizontal Discontinuities

In order to model the effects of volcanism in the upper crust, an overall background structure was modeled in accordance with available constraints. An important discontinuity in the upper crust is shown in a one-dimensional velocity model by Matumoto et al. [26] with a change on P wave velocities from 5.05 km/s to 6.2 km/s at a depth of between 8 and 10 km for northwestern Costa Rica. A discontinuity at an approximate depth of 10 km is also shown by a P wave velocity model from Gödde [4] based on wide angle refraction survey. Due to lack of constraining data for the Central region, this discontinuity was extrapolated to the area of interest and modeled as a change in density to 2.80 Mg/m3 below the 10 km depth. This density represents a mafic igneous basement assumed for the southwestern part of the Caribbean plate. A middle layer with densities between 2.70 Mg/m3 and 2.74 Mg/m3 was modeled to account for an andesitic composition of the next to uppermost crustal domain. Lateral variations were modeled to account for broad changes in the gravity field. An uppermost layer with a density of 2.60 Mg/m3 was modeled representing a 1-2 km thick layer comprised mainly of volcanoclastic Tertiary marine sediments, tephra deposits, and Quaternary alluvial sediments.

4.2. Effects of Magmatic Processes on the Upper Crust

The main features modeled beneath the Quaternary arc are elongated lower density bodies with densities ranging from 2.35 Mg/m3 to 2.38 Mg/m3. Such bodies represent low-density heterogeneities in the upper crust brought on by the effects of Quaternary volcanism. These effects may be comprised of a complex interaction between magma-derived components such as fluids and volatiles and the surrounding crust. Also, higher temperatures from higher heat flow along the active arc may play a role in lowering the overall density. In addition, the presence of melts may also be taken into account, although the relatively high volumes and broad lateral extension of the modeled bodies make it unlikely for these to be occupied entirely by melts. Thus, these are not interpreted to be magma chambers in their full extent but low-density zones in the upper crust directly related to processes such as hydrothermal alteration, higher heat flow, and the presence of melts.

As for the geometry and lateral variations of such bodies, the qualitative interpretation of the complete Bouguer anomaly map yields an indirect hint as for the disposition of these low-density zones. The main gravity low corresponds in extent and trend to the main volcanic edifices in the Cordillera Volcánica Central and appears as an elongated body constricted only in the region between the Barva volcano and the Irazú-Turrialba volcanic complex. This constriction coincides with a relatively larger gap between the main vents of the Barva and Irazú volcanoes (32 km) as compared with smaller more regular gap of about 15 km between the main vents of the Barva-Poas-Platanar volcanic edifices. This may indicate that the overlap of the effects of magmatism in the upper crust for each volcanic vent is attenuated by the greater distance represented by such gap. However, a nonconstricted gravity low may also be a product of the current sparse gravity station coverage which may merge the low wavelengths into a continuous signal. A low-density body (2.35 Mg/m3) was modeled separately for the Irazú-Turrialba volcanic complex; this separation concerns the basic geometry but maintains a constant density (Figures 7 and 8).

A joint low-density body (2.35 Mg/m3) was modeled for the upper crust beneath the Barva and Poás volcanoes (Figures 8 and 9). As for the underground structures of the Platanar-Porvenir volcanic complex, a slightly higher density of 2.38 Mg/m3 was calculated through inversion for the given geometry and assigned to the heterogeneous body resulting on a better fit between measured and calculated gravity. This is in accordance with the lower volume and lack of historical activity of the volcanic complex relative to the others, which may indicate a lower flux of magma and volatiles to the vent. Beneath the main heterogeneous low-density bodies (2.35 Mg/m3), a trenchward dipping low-density (2.68 Mg/m3) zone was modeled based on the orientation of low-velocity zones observed on the seismic tomography results. This zone is interpreted as a zone of passage of fluids and melts from the lower to the uppermost crustal domains.

As for the near fore-arc, local earthquake seismic tomography data (Arroyo [13]) show low-velocity zones beneath the paleo-arc. The location of these zones also shows a relative gravity low trending NW-SE which itself coincides with the surface geology units corresponding to the Neogene Aguacate volcanic arc. Constrained also by Euler deconvolution source points obtained for this work as well as the inferred location of paleo-volcanic vents (Alvarado [11]), low-density bodies (2.40–2.45 Mg/m3) were modeled for the Neogene volcanic arc to account for the effects of paleo-magmatic processes on the upper crust. The lower densities may be the result of pervasive hydrothermal alteration which can be recognized from surface geology. Also, the presence of granitic intrusions (Denyer and Alvarado [10]) which crop out in limited extent along the paleo-arc may contribute to this signal. A limited extent crustal body with a density of 2.80 Mg/m3 was modeled at the location of the Cedral mountains to account for the presence of an outcropping dioritic stock.

4.3. Back-Arc Structure

Towards the back-arc, a less constrained body with a density of 2.82 Mg/m3 (Figures 7, 8 and 9) was modeled as a basement high located beneath the Caribbean plains. Its structure follows the trend of a Bouguer gravity high in the back-arc and represents a shallow mafic basement in accordance with outcroppings of serpentinized peridotite near the Costa Rica-Nicaragua border (Tournon et al. [27]). The presence of such material in the back-arc at depths of less than 2 km corresponds with borehole stratigraphy logs (Pizarro [12]). The basement high marks the boundary between the San Carlos sedimentary basin related to the Nicaragua graben and the Northern Limón sedimentary basin towards the Caribbean passive margin. The sedimentary fill of the San Carlos basin was modeled as thickening of the uppermost layer (2.6 Mg/m3).

5. Conclusions

A total of 22 cross sections were designed to construct a constrained 3D density model for the upper crust in the central region of Costa Rica. Power spectrum analysis of the gravity field results in a maximum depth of approximately 11 km for sources of the local along-arc gravity anomaly modeled thus suggesting a shallow depth for the main heterogeneities beneath the Quaternary volcanic arc. Such heterogeneities were modeled by low-density (2.35–2.38 Mg/m3) zones which are interpreted as the effects of arc magmatism derived from the subduction of the Cocos plate beneath the Caribbean plate. The geometry of these bodies was constrained by seismic tomography data and Euler deconvolution source points. Along the Cordillera Volcánica Central, the low of the Bouguer anomaly is segmented towards the SE where shorter wavelengths are observed suggesting a change in the geometry of the low-density bodies. To account for this feature, separate low-density bodies were modeled for the Irazú-Turrilaba volcanic complex and the Barva and Poás volcanoes. However, these features consist of a unique 2.35 Mg/m3 density. Towards the NW a slightly higher (2.38 Mg/m3) density was modeled beneath the Platanar-Porvenir volcanic complex.

Low-velocity zones from local earthquake seismic tomography results from the work of Arroyo [13], and quantitative interpretation of the gravity field with Euler deconvolution source points suggests the presence of shallow heterogeneous low-density bodies beneath the Neogene “Aguacate” volcanic arc. These bodies were modeled with densities of 2.45 Mg/m3 for segments of the paleo-arc which have a predominantly tholeiitic composition. A density of 2.40 Mg/m3 was modeled for the segments of mainly calk-alkaline composition and with greater presence of pyroclastic rocks.

Along the back-arc, a 2.82 Mg/m3 body was modeled and interpreted as a basement high. This density accounts for the presence of a shallow mafic basement which consists of serpentinized peridotite which crops out in northeastern Costa Rica and is present in borehole logs at a shallow depth of approximately 2 km towards the back-arc.


The 3D density modeling presented here was supported by the Collaborative Research Centre “SFB574”: Volatiles and Fluids in Subduction Zones: Climate Feedback and Trigger Mechanisms for Natural Disasters which is financed by the Deutsche Forschungsgemeinschaft, Bonn, Germany, and the Christian-Albrechts-Universität Kiel, Germany. The authors sincerely acknowledge the support of colleagues in Kiel and San José: I. Arroyo for providing important seismological data for use as constraints within the model, S. Schmidt for her collaboration with the modeling software IGMAS, and T. Damm for his help with the 3D visualization and other SFB574 collaborators in Kiel and at the University of Costa Rica for providing valuable geological and geophysical insights into the investigated area. This is publication no. 188 of SFB574 and a contribution to the SPP 1257 Mass transport and mass distribution.