A Constrained 3 DDensity Model of the Upper Crust from Gravity Data Interpretation for Central Costa Rica

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


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.

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 arcparallel gravity low with a minimum of −57 × 10 −5 m/s 2 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/s 2 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/s 2 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/s 2 on the Herradura promontory.Towards the back-arc region, the Bouguer anomaly values increase to a relative high of 15 × 10 −5 m/s 2 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/s 2 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.

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

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

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

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, lowvelocity 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 V p from Arroyo [13].As for the geometry of such heterogeneities, the integration of the data as constraints in the density model shows trenchwarddipping, 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.
International Journal of Geophysics The green box shows the extent of the modeled area.Red triangles show the locations of main Quaternary volcanic vents.Blue triangles show the locations of paleo-volcanic vents for the Neogene "Aguacate" volcanic paleo-arc, inferred from surface geology and morphology.

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.• W 8 3 • W 8 5 • W 8 4 • W 8 5 Logarithm of power spectrum (a.u.) Figure 5: Results of the power spectrum analysis of the Bouguer gravity field using 2D Fast Fourier Transform.Black lines represent the trends of the linear regression separating the spectrum in a local trend (blue squares) and a regional trend (red diamonds).
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/m 3 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/m 3 and 2.74 Mg/m 3 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/m 3 was modeled representing a 1-2 km thick layer comprised mainly of volcanoclastic Tertiary marine sediments, tephra deposits, and Quaternary alluvial sediments.

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/m 3 to 2.38 Mg/m 3 .Such bodies represent lowdensity 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 lowdensity body (2.35 Mg/m 3 ) was modeled separately for the Iraz ú-Turrialba volcanic complex; this separation concerns the basic geometry but maintains a constant density (Figures 7 and 8).

International Journal of Geophysics
A joint low-density body (2.35 Mg/m 3 ) 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/m 3 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/m 3 ), a trenchward dipping low-density (2.68 Mg/m 3 ) 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.45Mg/m 3 ) were modeled for the Neogene volcanic arc to account for the effects of paleomagmatic 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/m 3 was modeled at the location of the Cedral mountains to account for the presence of an outcropping dioritic stock.[27]).The presence of such material in the backarc 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/m 3 ).

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.38Mg/m 3 ) 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 lowdensity bodies.To account for this feature, separate lowdensity 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/m 3 density.Towards the NW a slightly higher (2.38 Mg/m 3 ) 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/m 3 for segments of the paleo-arc which have a predominantly tholeiitic composition.A density of International Journal of Geophysics   2.40 Mg/m 3 was modeled for the segments of mainly calkalkaline composition and with greater presence of pyroclastic rocks.
Along the back-arc, a 2.82 Mg/m 3 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.

InternationalFigure 1 :
Figure 1: Tectonic setting of the Costa Rican subduction zone.Major plate boundaries and tectonic features are marked by red lines.Locations of the 3D density model is indicated by a black polygon.NPDB: North Panamá Deformed Belt; PFZ: Panamá Fracture Zone; CNSC: Cocos-Nazca Spreading Center.

9 C 8 84Figure 2 :
Figure 2: Complete Bouguer anomaly map of central Costa Rica draped over the shaded SRTM [7] relief map.Red triangles show location of the main Quaternary volcanic vents.Black plus signs indicate the location of gravity stations.Red lines show locations of vertical cross sections presented on this work.

TFigure 3 :
Figure3: Simplified surface geology of the study area based on cartography by Tournon and Alvarado[9] and Denyer and Alvarado[10].The green box shows the extent of the modeled area.Red triangles show the locations of main Quaternary volcanic vents.Blue triangles show the locations of paleo-volcanic vents for the Neogene "Aguacate" volcanic paleo-arc, inferred from surface geology and morphology.

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 N ic oy a Pe ni ns ul a O sa P en .C a ri b b e a n se a

Figure 4 :
Figure 4: Semiperspective view of the modeled area.Colored dots show locations of Euler deconvolution source points relative to the modeled bodies.Main modeled bodies are shown as colored polyhedra.The yellow polyhedron represents the back arc structural high with a density of 2.82 Mg/m 3 .The red polyhedra represent the low-density zones beneath the arc with densities between 2.35 and 2.38 Mg/m 3 .The blue polyhedra represent the paleo-arc bodies with densities between 2.40 and 2.45 Mg/m 3 .Coastline is shown in black.

Figure 6 :
Figure 6: Comparison of measured (a) and calculated (b) complete Bouguer anomalies.Location and trend of model cross sections are shown as NE-SW trending black straight lines.The cross sections presented on this work are highlighted in red.Location of the main Quaternary volcanic vents are shown as red triangles.The Neogene "Aguacate" paleo-volcanic vents inferred from surface geology are shown as blue triangles.Black dots represent the location of gravity stations.Coastline and political borders are shown in bold black lines.

Figure 7 :
Figure 7: 3D density model cross section across the Poás volcano and surrounding regions.Red and black curves in the upper box show the measured and calculated gravity anomalies, respectively.Black plus-signs indicate topography derived from SRTM 3 arc sec grid which is projected on the cross section; black stars in the lower box show Euler source points.Grayscale polygons represent modeled crustal bodies with corresponding densities.Colored pixels and bars show percentage of change in V p velocities which were calculated by Arroyo [13].Red triangle indicates the location of the main Quaternary volcanic vent; blue triangle indicates the locations of inferred paleo-volcanic vent.

Figure 8 :
Figure 8: 3D density model cross section across the Iraz ú volcano and surrounding regions.Red and black curves in the upper box show the measured and calculated gravity anomalies, respectively.Black plus-signs indicate topography derived from SRTM 3 arc sec grid which is projected on the cross section; black stars in the lower box show Euler source points.Grayscale polygons represent modeled crustal bodies with corresponding densities.Colored pixels and bars show percentage of change in V p velocities which were calculated by Arroyo [13].Red triangle indicates the location of the main Quaternary volcanic vent.

Figure 9 :
Figure 9: 3D density model cross section across the Barva volcano and surrounding regions.Red and black curves in the upper box show the measured and calculated gravity anomalies, respectively.Black plus-signs indicate topography derived from SRTM 3 arc sec grid which is projected on the cross section; black stars in the lower box show Euler source points.Grayscale polygons represent modeled crustal bodies with corresponding densities.Red triangle indicates the location of the main Quaternary volcanic vent; blue triangle indicates the location of the Cedral inferred paleo-volcanic vent.