3‐D Modeling of Vertical Gravity Gradients and the Delimitation of Tectonic Boundaries: The Caribbean Oceanic Domain as a Case Study

Geophysical data acquisition in oceanic domains is challenging, implying measurements with low and/or nonhomogeneous spatial resolution. The evolution of satellite gravimetry and altimetry techniques allows testing 3‐D density models of the lithosphere, taking advantage of the high spatial resolution and homogeneous coverage of satellites. However, it is not trivial to discretise the source of the gravity field at different depths. Here, we propose a new method for inferring tectonic boundaries at the crustal level. As a novelty, instead of modeling the gravity anomalies and assuming a flat Earth approximation, we model the vertical gravity gradients (VGG) in spherical coordinates, which are especially sensitive to density contrasts in the upper layers of the Earth. To validate the methodology, the complex oceanic domain of the Caribbean region is studied, which includes different crustal domains with a tectonic history since Late Jurassic time. After defining a lithospheric starting model constrained by up‐to‐date geophysical data sets, we tested several a‐priory density distributions and selected the model with the minimummisfits with respect to the VGG calculated from the EIGEN‐6C4 data set. Additionally, the density of the crystalline crust was inferred by inverting the VGG field. Our methodology enabled us not only to refine, confirm, and/or propose tectonic boundaries in the study area but also to identify a new anomalous buoyant body, located in the South Lesser Antilles subduction zone, and high‐density bodies along the Greater, Lesser, and Leeward Antilles forearcs. Plain Language Summary The knowledge of the density structure of the different layers that compose the solid Earth is important, for example: in the study of earthquakes, in plate tectonics reconstructions, or for the modeling of petroleum systems. These density variations affect (in small scale) the intensity of the gravity field on each point of the Earth's surface. The gravity field can be globally measured with satellites, reaching areas where the direct measurements are expensive and time consuming, such as in the ocean. In this work, we propose a new methodology with the purpose of recognizing tectonic and/or terrain limits, located in the outer most layer of the solid Earth, named crystalline crust. We calculate the gravity field of different density distributions, using four layers: seawater, sediments, crystalline crust, and mantle (a layer located below the crust), and compare the results with satellite global measurements. With our methodology it is possible to refine, confirm, and/or propose terrain limits, but additionally, we are able to estimate the average density configuration of the crystalline crust. This methodology is validated in the oceanic domain of the Caribbean, where a complex geologic history exists, due to its evolution since approximately 144 million years ago.


Introduction
The oceanic lithosphere, which includes the upper mantle, the crystalline crust and overlying sediments, is one of the least known features on outer Earth. Since the development of echosounders in the 1940s, only 10% of the sea floor has been mapped at 1 min resolution (Becker et al., 2009). Despite the current direct techniques for oceanic lithospheric mapping yield the most accurate results (e.g., refraction and reflection seismic, sidescan sonar and swath bathymetry), data acquisition is usually expensive and time consuming, due to the technical and logistical requirements involved in oceanic expeditions. Additionally, most of the tectonic structures of the oceanic crust have been mapped (inferred) from the sea-floor topography (e.g., Rosencrantz & Mann, 1991), although its spatial distribution at levels within the crystalline crust might be different. The identification of more accurate tectonic boundaries in the oceanic crust would add valuable information, for example, to test Global Positioning System (GPS) geodetic-based models, to improve 3-D configurations of gravity models, to identify boundaries between oceanic and continental crust, and as input parameters for kinematic reconstructions back in time.
Combined data sets of satellite gravity, altimetry, and terrestrial measurements make it possible to develop and test 3-D lithospheric models with homogeneous coverage and high spatial resolution (Götze & Pail, 2018), which depends on the limitations of the gravity data source: terrestrial measurements over continents and altimetry over the oceans. These lithospheric models, however, need to be tested against independent data sources to restrict the multiple options of densities, sizes, shapes and depths of bodies, which can reproduce similar gravity signals (Lowrie, 2007;Maystrenko et al., 2013).
The vertical gravity gradients (VGG or T zz -details in section 3.3) complement the interpretation of the signal given by the more commonly studied Bouguer and free-air gravity anomalies. This is due to the fact that the VGG are especially sensitive to shallow density variations, while the gravity anomalies are more useful in the study of deeper density contrasts (Álvarez et al., 2014, 2015). Recently, more researchers are giving attention to the use of satellite gravity gradients for modeling purposes (Álvarez et al., 2012Ebbing et al., 2013;Götze & Pail, 2018;Oruc, 2014;Schaller et al., 2015). However, most of the studies interpret the signal with simple (Bouguer) density corrections, without considering the heterogeneities of the different layers of the lithosphere.
Recently, Götze and Pail (2018) stated that an additional advantage of the VGG is that they can highlight with better resolution the edges of geological structures, intrusions, faults, or even the continental-oceanic transition (COT) at continental margins, where the density contrast may be higher due to the transition from continental to oceanic crust.
The definition of the COT is important for several reasons. For example, from an economic point of view, the modeling of petroleum systems can be improved with the knowledge of the crustal type, because the radiogenic heat production and thermal conductivities of the different crusts affect the hydrocarbons generation and maturation processes (Pawlowski, 2008). Furthermore, the shape of the COT can add valuable information in plate tectonics reconstructions (Götze & Pail, 2018).
In order to evaluate the use of the VGG in the delineation of tectonic boundaries using 3-D density and structural models of the lithosphere, we selected the Caribbean oceanic region (Figure 1). The crustal structure of the Caribbean is the result of a north-northeastward migration of different crustal domains (e.g., volcanic arcs, continental and oceanic realms) since Late Jurassic to Early Cretaceous times, including the igneous plateau materials that affected the oceanic crust of the former Farallon plate. Nowadays, the Caribbean plate is an active region with four subduction zones, three deformed belts, and more than 25 identified terrains (details in section 2; see Figures 1 and 2). Different regions of the Caribbean have been the target of relatively extensive seismic reflection and refraction, sonar, and drilling campaigns (e.g., Edgar et al., 1971;Kroehler et al., 2011;Mauffret & Leroy, 1997;Rosencrantz, 1990), some of them undertaken with the limitations of early seismic campaigns. Nevertheless, the coverage of these measurements is poor when compared with the complexity of most of the Caribbean morphological structures (Rosencrantz, 1990), and in other cases, the large thicknesses of sedimentary deposits have made it impossible to drill the underlying crust (e.g., Pichot et al., 2012). Therefore, its complex tectonic and terrain setting make the Caribbean one of the most interesting natural laboratories on Earth (Jiménez-Díaz et al., 2014).
Some of the topics for which no consensus has been reached in the scientific community about the tectonics of the Caribbean region include the plate boundary between the oceanic North and South American plates in the Lesser Antilles subduction (Deville et al., 2015;Van Benthem et al., 2013). Similarly, the limit between the Colombian and Venezuelan basins, which has been the topic of debate of different authors (e.g., Mattioli et al., 2014;Symithe et al., 2015) since a two-plate kinematic model was proposed for the Caribbean by Dewey and Pindell (1985). An important source of error might be the a priori assumed shape of this edge, which has been different for every kinematic model. Additionally, the crustal structure of regions such as 10.1029/2019GC008340 Geochemistry, Geophysics, Geosystems the Aves Ridge and the Lesser Antilles remains poorly constrained and data acquisition has been limited to a (relatively) small region (e.g.,: Bangs et al., 2003;Christeson et al., 2003;Evain et al., 2013;Laigle, Becel, et al., 2013, Laigle, Hirn, et al., 2013.
The main goal of this study is to explore the additional information that satellite derived VGG provide about the structure of the oceanic crust in the Caribbean domain, considering that the usefulness of these gradients is not widely recognized, and that they are especially sensitive to density contrast in the upper layers of the Earth. Specifically, we address the questions about unknown or not well-defined tectonic boundaries and the identification of crustal bodies with anomalous densities that have not been described before.
The work hypothesis is based on the analysis of the residuals resulting from the gravity gradients derived from the EIGEN-6C4 potential field model (Förste et al., 2014) downloaded from ICGEM (2018), and the gravity gradients response of 3-D density and structural models of the Caribbean lithosphere. Specifically, our hypothesis states that the residuals contain information about the structure and density heterogeneities at crustal level, as the crystalline crust is the only free parameter that we will consider in our model. Moreover, the edges of crustal fragments with a sufficiently high-density contrast should have an associated sign change on the VGG (Álvarez et al., 2012), a property that we use in the interpretation of the residuals.
Thus, a sensitivity analysis to different density solutions for the water and sedimentary layers was carried out (supporting information), with a constant density initially assumed for the crystalline crust. From all the tested density configurations, we selected the more robust 3-D model, or in other words, the model with the minimum misfits compared with EIGEN-6C4 data set. Such sensitivity analysis allowed us to minimize the error associated with density heterogeneities of the upper lithospheric layers (water and sediments).

Geochemistry, Geophysics, Geosystems
Additionally, an average crustal density field was inferred from the forward modeling of the VGG of this selected 3-D model. Based on this approach, we refine, propose, or confirm the location of plate boundaries, such as the heavily debated limit between the North and South American plates in the Lesser Antilles subduction zone, the boundary between the Colombian and Venezuelan basins toward the west flank of Beata Ridge, the limit between the Modified Saba (MS) crust and the oceanic floor of the Grenada Basin (GB), and the COT along the continental margins of the study area. We also identify a new anomalous buoyant body (volcanic arc fragment?) located in the oceanic crust of the Atlantic Ocean, close to the South Lesser Antilles subduction zone. Similarly, we recognize high-density bodies broadly distributed along the Greater, Lesser and Leeward Antilles forearcs, as well as low-density crustal fragments inside the Yucatan Basin (YB) and Cayman Trough (CT).

Subduction Zones
The oceanic floor of the Caribbean region is surrounded by four subduction zones (Figure 1). To the north, the North Lesser Antilles (NLA) subduction forms the Puerto Rico Trench (PRT - Figure 2), an area not only characterized by the south-south west dipping of the slab, but also by the transition to a transform boundary between the North American and Caribbean plates to the west of the trench (Boschman et al., 2014). As we move toward the eastern side of the Puerto Rico slab, the NLA subduction curves to the south east, dipping in a westward direction, giving as a result an amphitheatre-like geometry (Van Benthem et al., 2013).
Toward the east of the Caribbean plate, North and South American oceanic crusts subduct beneath the Caribbean at a rate of 2 cm/year (Calais et al., 2016). The tectonic edge between the North and South American plates is well known at depths of ≈400 km, where a region of low P-wave velocity is the main piece of evidence of a slab gap within 13°N and 15°N (Van Benthem et al., 2013). Despite the plate boundary at crustal level remains a topic of debate, the limit at mantle levels found by Van Benthem et al. (2013) partially overlaps with the location of the Tiburon Fault Zone (TFZ) and the Barracuda Fault Zone (BFZ; Figure 2). The deformation along the BFZ created the Barracuda Ridge, whose uplift and flexure started in the Early Pleistocene and which remains active nowadays. The deformation processes also affected the evolution of the Barracuda Trough, a deep basin north of the Barracuda Ridge that has been proposed as the  (Förste et al., 2014). Color scale saturated at −300 10 −5 m/s 2 and +300 10 −5 m/s 2 (1 mgal = 10 −5 m/s 2 ). AG = Aruba westernmost (complex and diffuse) plate boundary between the North and South American plates (Pichot et al., 2012). Based on diverse geophysical evidence, other authors have proposed this plate boundary at different locations on the crust, for example at ≈19°N (Bird, 2003 -yellow line in Figure 1) and between 14 and 17°N (Müller & Smith, 1993).
To the southern edge, the Caribbean plate has subducted beneath the continental South American plate and the Maracaibo block since the Late Cretaceous (Kroehler et al., 2011), forming the South Caribbean Deformed Belt (SCDB; Figures 1 and 2). Finally, the Middle American Trench (MAT; Figure 1) is the western most boundary of the Caribbean plate and is connected with the Panama-Chocó block, which toward the Caribbean forms the Panama Deformed Belt (PDB; Figures 1 and 2; Buchs et al., 2010;Montes et al., 2012).

Volcanic Arcs
In the Caribbean oceanic region at least four volcanic arcs, with different geological histories, have been recently described (e.g., Neill et al., 2011). These arcs include the youngest and currently active Lesser Antilles (LA; Figure 2), and three additional extinct arcs that formed the current Leeward Antilles (≈88-70 Ma), Aves Ridge (AR, ≈80-75 Ma), and the Great Arc of the Caribbean (GAC, ≈100-96 Ma).
Aves Ridge is located in a northsouth direction in the eastern Caribbean. Despite the large portion of sea floor that this ridge occupies (≈65 x 10 3 km 2 ), details about the magma sources have not been well defined yet (Neill et al., 2011). Christeson et al. (2008) described Aves as an arc of intermediate composition, although a few of the samples include granitoids and mafic rocks (Neill et al., 2011). The volcanic activity of Aves migrated to the east during the Paleocene to Early Eocene due to a rollback process of the Proto-Caribbean slab. This caused not only the opening of the Grenada Basin (GB; Figure 2) but also the formation of the currently active volcanic arc of the Lesser Antilles, and as a consequence, the Barbados Accretionary Prism (BAP; Figure 2; Aitken et al., 2011). The rollback velocity of the Proto-Caribbean slab must have been slower in the northern portion of the GB, allowing the flow of molten material and creating the Modified Saba crust (MS; Figure 2; Arnaiz-Rodríguez & Audermard, 2018). The Lesser Antilles are part of an arc of intermediate composition (Christeson et al., 2008) which has been active since the Early Eocene, with a gap in the volcanic activity of 8 to 10 Myr, due to the subduction of a buoyant body during the Late Oligocene. As a result, the axis of the volcanic arc shifted toward the west in the region north of Martinique Island (Bouysse & Westercamp, 1990). Another highly buoyant ridge has been described seaward of Guadeloupe Island by Bangs et al. (2003). This ridge accreted to the original arc crust of the Antilles in the late Miocene, displacing the backstop trenchward. Nowadays, the subduction of three nonbuoyant ridges (Santa Lucia, Tiburon and Barracuda -SLR, TFZ and BFZ; Figure 2) is broadly recognized Christeson et al., 2003;Laigle, Becel, et al., 2013).
On the other hand, Evain et al. (2013) found that the Lesser Antilles forearc is composed of two different crustal domains: an inner and an outer forearc, characterized by high and low seismic velocity gradients, respectively. Such crustal domains interact spatially with the down-going plate in different ways, for example, acting as backstops; therefore, their regional characterization will be useful in the understanding of the subduction dynamics. Nevertheless, the detailed seismic experiments carried out in the Lesser Antilles subduction have focused on a relatively narrow region between 15-18°N and 59-62°W (i.e., Bangs et al., 2003;Christeson et al., 2003;Evain et al., 2013;Laigle, Becel, et al., 2013, Laigle, Hirn, et al., 2013; thus, much of the structure of the currently active volcanic arc remains unconstrained. The origin of the Leeward Antilles involved the interaction of island arcs of depleted or mantle plume related sources with oceanic plateau processes (Neill et al., 2011). During the north-northeast migration of the Caribbean Large Igneous Plateau (CLIP), fragments of these island arcs, as well as pieces of Proto-Caribbean oceanic crust -for example picrites found in Curaçao Island (Mauffret & Leroy, 1997), were accreted along the continental margins of both the North and South American plates (Boschman et al., 2014).

Beata Ridge and the Colombian and Venezuelan Basins
The Beata Ridge (BR; Figure 2) dates from Cretaceous time and its crustal structure consists of volcanic material of plateau origin, which includes from top to bottom: basalts, original oceanic crust, a gabbroic layer and a picritic layer with mafic cumulates (see Figure 24 in Mauffret & Leroy, 1997). The ridge is an area of active Caribbean intraplate deformation since at least the Early Miocene, according to the activity of Pecos Fault Zone (Leroy & Mauffret, 1996). Beata lies between the Colombian and Venezuelan basins (CB and VB; Figure 2), and based on its deformation (which increases toward the north and west), its shape and focal mechanisms, Leroy and Mauffret (1996) suggested a relative displacement between both basins (microplates?). These authors proposed that the boundary between these (hypothetical) microplates is located at the west flank of Beata (orange line in Figure 1), dominated by left-lateral strike-slip motion with a component of compression. Similar evidence was found recently by Mattioli et al. (2014), whose best fitting inversion to GPS measurements requires a two-plate model of the Caribbean oceanic crust, which has an internal deformation of 1-3 mm/year. According to multichannel seismic profiles interpreted by Leroy and Mauffret (1996), this ridge has strong transpressive tectonic features. Therefore, the differential motion between both basins might be absorbed as a collision in southern Hispaniola Island (H; Figure 2) and transpression in the Beata Ridge area. Moreover, Symithe et al. (2015) tested this boundary locating it toward the east flank of Beata (purple line in Figure 1) on their GPS kinematic model, and found a better fitting inversion with the Caribbean acting as a single plate. Nevertheless, even though the Aruba Gap (AG; Figure 2), and especially the Pecos Fault Zone, has been suggested as the southernmost boundary between Beata and the Colombian Basin (Leroy & Mauffret, 1996), the plate boundary for the two-plate model described before remains unknown (Boschman et al., 2014).

Yucatan Basin and Cayman Trough
The YB ( Figure 2) is bounded to the north by the Yucatan borderlands and the continental platform of Cuba Island (C; Figure 2), and to the south, by the CT. Rosencrantz (1990) identified different basement domains, distinguished by three types of crust below this basin: (1) the Yucatan borderlands, an old passive margin which includes metasediments in the continental platform and slope.
(2) The Western Deep Basin, which has been interpreted as a large and inactive pull-apart system; this basin is underlain by thinner oceanic crust (less than 10 km according to Mauffret & Leroy, 1997; Figure 1) of Paleocene to Middle Eocene ages. And finally, (3) the oceanic floor of the Central Seamounts, and the Cayman Rise and Ridge (see numbers 1 to 4 in Figure 2), whose crusts are dominated by volcanics, metavolcanics, plutonics, and granodiorites dated from the Late Cretaceous (Boschman et al., 2014;Rosencrantz, 1990). These volcanic domains originated from magmatic flows during the opening of the YB, and therefore, are characterized by a heterogeneous topography.
The CT ( Figure 2) is an oceanic spreading zone dated as Early Eocene in age  which records the motion between the North American plate (Cuba segment) and the Gonave Microplate (GM; Figure 2; Boschman et al., 2014). This pull-apart basin is bounded to the south by a complex of sinistral strike-slip faults, which includes the Motagua, Swan Islands (SIFZ), Walton (WFZ) and Enriquillo-Plantain Garden faults (EFZ). To the north, the trough is bounded by the Oriente (OFZ) and the Septentrional (SFZ) fault systems (Rosencrantz & Mann, 1991). To the east of 77°W, the oceanic floor of CT is dominated by offshore rifts in a northeast-southwest orientation, dated as Paleocene -Early Eocene , and by the transition from oceanic to continental crust toward Hispaniola Island ten Brink et al., 2002).

Input Data
The lithospheric starting model included the thicknesses of four layers, namely from the uppermost to the lowermost: (1) seawater, derived as the difference between sea level and the General Bathymetric Chart of the Oceans (GEBCO; Weatherall et al., 2015); (2) sediments, published in the NOAA Total Sediment Thickness data set (Whittaker et al., 2013); (3) crustal thickness, computed as the difference between the top of the crystalline crust (topography in continental areas and bathymetry minus sediment thickness in oceanic domains) and Moho depth from the GEMMA model (Reguzzoni & Sampietro, 2015); and finally, (4) the mantle, which was subdivided into eight layers using the SL2013sv S-wave velocity model (Schaeffer & Lebedev, 2013), from the Moho down to 200-km depth.
An alternative model of the Moho depth is CRUST1.0 (Laske et al., 2013), but we used GEMMA instead because of four reasons: (1) GEMMA uses a uniform and homogeneous coverage of satellite measurements, while CRUST1.0 is based on seismic information, which is absent for large areas.
(2) Despite both GEMMA and EIGEN-6C4 use GOCE data, they are sufficiently independent from each other, since the latter combines measurements from diverse sources: the LAGEOS, GRACE and GOCE satellites, terrestrial information, as well as satellite altimetric data (over the oceans; Förste et al., 2014). (3) The spatial resolution of GEMMA (0.5°×0.5°) is four times smaller than CRUST1.0 (1°×1°), which provides a great improvement at the scale of our regional study. (4) Compared to CRUST1.0, GEMMA provides a more realistic Moho geometry, which allows us to minimize the uncertainty about other error sources in the forward modeling process.
The integration of the different data sets was made after a cubic interpolation to a homogeneous spatial resolution of 0.05°. Figure Figure 3c shows the Moho depths from the GEMMA model, which considers global gravity data from the GOCE satellite mission, as well as additional external information, including a model of lateral density variations in the upper mantle (Reguzzoni & Sampietro, 2015). The approach used in the computation of GEMMA allowed the calculation of Moho depths in the study area with an error standard deviation ranging between ≈1.5 km (in the deep ocean basins) and ≈7 km (in the continental platform west of YB; Figure S4). The bathymetric highs of the Nicaraguan Rise and the paleoarc of Aves Ridge are also characterized by deep Moho interfaces, between 20 and 35 km. The shallowest Moho depths are located in the YB, CT, and Colombian and Venezuelan basins, following the pattern previously described by Mauffret and Leroy (1997), where a crustal thickness of less than 10 km was reported ( Figure 1). Finally, the thicknesses of the crystalline crust are plotted in Figure 3d. In the oceanic domain of the Caribbean, the thickest crust corresponds to the paleoarc of Aves Ridge (≈35 km), followed by the North and South Nicaraguan rises where the thickness is ≈30 km.
The modeled VGG are compared with the VGG obtained from EIGEN-6C4 combined data set (Förste et al., 2014). This data has a spherical harmonic solution up to degree and order N max =2190; therefore, the smallest topographic wavelength that EIGEN-6C4 is able to resolve can be calculated as λ min = 2πR/N max ≅ 18 km. Where N max is the maximum degree and order of the spherical harmonic solution, and R is the mean Earth radius (Álvarez et al., 2012).

Density Solutions for Lithospheric Layers
A sensitivity analysis to different density solutions for water and sediments was carried out (supporting information). Because the mantle contribution to the total VGG field is small ( Figure S1), we did not include it on this analysis. The density of the seawater was defined following two different approaches: (1) a 3-D distribution using the mathematical model of Gladkikh and Tenzer (2012), which predicts density based on latitude and water depth, calculated for 84 layers, 100 m thick; and (2) a constant density of 1,030 kg/m 3 , a typical value previously used in other gravity applications (e.g.: Álvarez et al., 2015;Maystrenko et al., 2013;Tenzer et al., 2010). In the 3-D mathematical model, the authors used data from the World Ocean Atlas 2009 and from the World Ocean Circulation Experiment 2004 to obtain an equation that can calculate seawater density with a relative accuracy of ≈0.25%. The density distribution ρ (kg/m 3 ) can be written in the following form (Gladkikh & Tenzer, 2012): where D is the ocean depth (m) and φ the geographical latitude (degrees). This approach considers the latitudinal contribution α(φ) in the seawater density, or in other words, the variations due to salinity and temperature as the water moves from the equator to the poles. Additionally, the change of density with depth (pressure) is included in the term β(φ)D υ(φ) . Finally, the model considers a pycnocline correction μ(φ), necessary for a better approximation of the density structure at shallow depths, where the water temperature decreases and the salinity increases with depth. The parameters described above are expressed as: The computed water densities range from 1,022.5 kg/m 3 at the ocean surface and 1,064.5 kg/m 3 at the ocean bottom ( Figure S3a).
Additionally, four different sediment-density structures were tested in the sensitivity analysis: (1) a 3-D density distribution, based on the equation of Tenzer and Gladkikh (2014), computed for 181 layers, 100 m thick each; (2) a constant value of 2,350 kg/m 3 , which is the mean value of the 3-D density model; (3) a highdensity value of 2,610 kg/m 3 , and finally, (4) the average density of marine sediments 1,700 kg/m 3 , reported by Tenzer and Gladkikh (2014). The density equation proposed by these authors (approach (1)) is the result of the analysis of 716 drill sites from the Deep Sea Drilling Project (DSDP), and can be written as where ρ s is the sediment density (kg/m 3 ), D is the ocean depth (m), d s is the sediment depth (m), and 1.66 is the nominal value of sediment density at sea level. This approximation takes into account not only the lateral density variation (second term), but also the increase of sediment density with sediment depth (compaction -third term). In general, with this approach it is possible to calculate the sediment density with an uncertainty of about 10% compared with the DSDP samples (Tenzer & Gladkikh, 2014). However, model restrictions due to limitations in the drilling depths of DSDP cores include a maximum ocean depth of 7 km and a maximum sediment thickness of 1.7 km. In the study area, the maximum water depth is approximately 8.4 km in the NLA subduction (Figure 3a), which corresponds to sedimentary deposits of less than 1 km thick ( Figure 3b). Nevertheless, the sediment thickness reaches a maximum of ≈18 km in the BAP (Figure 3b), so additional caution must be taken into consideration in order to compute the 3-D distribution of sediment densities using this approach. In our case, we applied the bedrock density contrast correction of Chen et al. (2014) to 27% of the data, which was out of the limit set by the density model. Thus, the maximum density allowed for sediments was set as 2,750 kg/m 3 , the maximum expected density for shales (Schön, 2011).
The final density distribution with depth, calculated using equation (6) and corrected due to bedrock density contrast is shown in Figure S3b. The light blue area represents the range of the computed densities at different depths. The density reaches a constant value of 2,750 kg/m 3 (bedrock contrast correction) at ≈2-km depth.
Even though the variation of the sediment density with depth strongly depends on the type of sediment (e.g., Hamilton, 1976), the sedimentary processes in the Caribbean region range from clastic-dominated (such as the Magdalena and Barbados depocenters), to biogenic-rich deposits (i.e.: in islands located in the continental platform), to dominantly pelagic (in the open ocean). This diversity does not allow an explicit modeling of the sediment composition, and we consider that a more realistic approach is given by the regression of measured core densities expressed in equation (6).
Initially, the crustal density was considered homogeneous for domains over and under the mean sea level, with a constant value of 2,810 kg/m 3 (continental) and 2,900 kg/m 3 (oceanic), respectively. These constant values do not add further disturbances to the modeled VGG at crustal level, as might be the case of a hypothetical 3-D density distribution. Therefore, the interpretation of the residuals can be undertaken under the assumption that they are a response of density heterogeneities in the crystalline crust, that were not considered in the initial set up of the starting models.
The mantle densities were computed in eight layers, from 25 to 200-km depth, using the S-wave velocity anomalies of Schaeffer and Lebedev (2013). We converted the shear wave velocities to densities following a modified approach of Goes et al. (2000), using a pressure and temperature dependent expansion coefficient (Hacker & Abers, 2004;Meeßen, 2017). The conversion requires the composition of the mantle, which was defined as the composition of oceanic upper mantle, modified from Shapiro and Ritzwoller (2004 ; Table S3). Table 1 summarizes the structural layers and the different densities tested in the starting models. Finally, Table S1 shows the configuration of the six different starting models that were considered in this study, in order to evaluate the sensitivity of the VGG to the density distribution of the water and sediments, as previously discussed.

VGG and the Identification of Tectonic Boundaries
The VGG are the second derivatives of the disturbance potential in the direction of the Earth radius, and represent the vertical component of the Marussi tensor which quantities are expressed in Eötvös (1E=10 −9 /s 2 ). In oceanic domains, the VGG signal is affected by the density contrast created by different layers at different depths, that can be written as: where T zz W accounts for the signal of water, T zz S for the sediments, T zz C for the crystalline crust, and T zz M for the mantle components. The VGG from the EIGEN-6C4 data set (T zz EIGEN − 6C4 ) represent the total signal of the VGG in the study area. On the other hand, the calculation of the total VGG based on an a priori or Typical constant density used in other studies (e.g., Álvarez et al., 2015;Maystrenko et al., 2013;Tenzer et al., 2010). c Using equation of Tenzer and Gladkikh (2014) and corrected following Chen et al. (2014). d Mean value of 3-D distribution. e A priori constant high density. f Average density of marine sediments based on DSDP cores (Tenzer & Gladkikh, 2014). g Continental crust based on Maystrenko et al. (2013). h Average density of the ocean bedrocks based on DSDP cores (Tenzer & Gladkikh, 2014).
i Velocity to density conversion using Goes et al. (2000) approximation, implemented by Meeßen (2017). j Grid resolution of 0.5°differs from model resolution, which is parameterized on a triangular grid with an ≈280-km spacing.

10.1029/2019GC008340
Geochemistry, Geophysics, Geosystems starting 3-D lithospheric model gives the T zz Modeled , and includes all the components described before, whose structural and density information was constrained using up-to-date high spatial resolution geophysical databases (see sections 3.1 and 3.2).
Taking into account that our goal is to analyze crustal heterogeneities, the density distribution of the crystalline crust was assumed constant in the starting models (section 3.2 and the supporting information). This approach allows us to hypothesize that the residuals of the VGG (T zz Residuals , equation (8)) can be interpreted as evidence of density heterogeneities (density contrasts) in the crystalline crust, which were not initially considered in the starting models.

Modeling Approach
The VGG were computed in spherical coordinates, at a calculation height of 1 meter above the topography or sea level, on a regular grid of 0.05°, using the software Tesseroids (Uieda et al., 2016). This spherical approximation is valid in regions close to the equator, where the radius of the sphere (6,378,137 m) approximates the semimajor axis of the WGS84 reference ellipsoid (Li & Götze, 2001). Therefore, the comparison between EIGEN-6C4 (referenced to the ellipsoid) and the modeled results (referenced to the sphere) is valid in the Caribbean domain. The modeled area is located between 5-25°N and 55-95°W; however, in order to avoid boundary effects, an internal buffer of 2°was established for the study area ( Figure 1). The 3-D density and structural models were converted to tesseroids: elementary bodies delimited by two meridians, two parallels, and two concentric spheres (Uieda, 2015). The advantage of using Tesseroids lies in the fact that this software allows the computation of the gravity derivatives in spherical coordinates. This avoids the use of the flat Earth approximation, which in a large area like our study domain (≈2,000 km wide by ≈4,200 km long) can be a significant source of error, as demonstrated for example, by Álvarez et al. (2012). Additionally, a high-pass Gaussian filter was applied to the modeled results, using GMT software (Wessel et al., 2013), with the purpose of homogenizing the minimum wavelength of the model with the minimum wavelength that can be resolved by EIGEN-6C4 data set (≈18 km). Figure S5 shows the comparison between the non-filtered residuals and the results after applying the Gaussian filter.

VGG Residuals and the Identification of Tectonic Boundaries
After the sensitivity analysis was carried out, we selected the density solution of the lithosphere that showed the best fit between the modeled VGG and the VGG from EIGEN-6C4 data set, considering a constant density for the crystalline crust. In this case, it corresponds to the density distribution of the starting model SM6, which includes 3-D density distributions for both the water and the sediments (supporting information). Thus, based on the VGG residuals obtained from this starting model, we attempted to identify heterogeneities in the crystalline crust, following our work hypothesis. We analyzed these residuals considering the previously described geologic setting of the Caribbean (section 2), and the provinces defined by Case and Holcombe (1980).
The edges of crustal fragments with a sufficiently high-density contrast should have an associated sign change of the VGG residuals, from positive to negative (Álvarez et al., 2012). We use this property for refining, proposing or confirming major tectonic boundaries in the Caribbean oceanic domain.

Forward Modeling of the VGG and the Inferred Average Crustal Density Field
To test our work hypothesis, we forward modeled the VGG of the best fit starting model (SM6, supporting information). In this model, the density assigned to the oceanic crystalline crust was 2,900 kg/m 3 (Table  S1). Thus, the sign of the residuals is directly related with a deficit or excess with respect to this reference value. Specifically, positive residuals indicate that the real crustal density is higher compared with the reference (and vice versa with the negative ones). Therefore, an iterative modeling process of the VGG was carried out, in which we modified the crustal density ± 10 kg/m 3 on each step (with respect to the reference density) considering the sign of the residuals on each cell of the grid (grid searching). Additionally, we also considered a representative density for the different provinces previously identified in the Caribbean region. Such as: 2,800 kg/m 3 for island arcs of intermediate composition (e.g., Aves Ridge, Leeward and Lesser Antilles), 3,000 kg/m 3 for picrites (present in Curaçao and the west flank of Beata Ridge), 2,810 kg/m 3 for continental crust (such as the North Nicaraguan Rise), and 2,800 kg/m 3 for basaltic flows (CLIP). In each iteration, the residuals of the VGG were computed, and again the crustal density at each cell was modified aiming to reduce the residuals. Eventually, the iterations converged to a density model which includes a heterogeneous density distribution for the crystalline crust, and whose residuals had the minimum standard deviation and root-mean-squared error. We called this final model SM-F (Table S1). Figure 4 shows the general methodological workflow used in this study for the analysis of tectonic boundaries, and for the inference of the average distribution of crystalline crust densities, using the VGG.

VGG in the Caribbean Region
From the EIGEN-6C4 field, some structural features of the different Caribbean oceanic terrains can already be identified (Figure 5a). The minimum values (negative gradients) are generally distributed along the subduction zones (e.g., MAT and NLA), but also along the Oriente and Septentrional faults, in the northern edge of GM. The maximum values (positive gradients), on the other hand, appear mainly over bathymetric highs such as the Leeward and Lesser Antilles, Beata and Aves ridges, and some broadly distributed highs on Cayman, North and South Nicaraguan rises and along the continental platforms of Bahamas and Central America. It is worth mentioning that the Barracuda Ridge has a high, positive signal, whereas its associated deep basin, the Barracuda Trough, is characterized by strong, negative VGG.
The modeled and filtered results (see section 3.3.1) obtained with the density structure of the starting model SM6 (Figure 5b and Table S1) share most of the characteristics previously described for the EIGEN-6C4 data set. The negative gradients are associated with bathymetric lows such as trenches, troughs and small subbasins, but in this case, the extreme changes in bathymetric gradients, as for example those that occur across the continental slopes (e.g., Yucatan and Venezuelan basins), are especially highlighted by negative values. In general, the continental platforms are characterized by positive modeled gravity gradients, in contrast with EIGEN-6C4 data, where short wavelengths of positive and negative signals are present. Even though both EIGEN-6C4 and the modeled and filtered VGG include spherical harmonic solutions up to degree and order 2190 (λ≈18 km), the calculated gradients do not show small density contrasts that can be broadly evidenced throughout the study area in the EIGEN-6C4 data set. The smallest and well resolved features on the model correspond with seamounts concentrated on the Pacific Ocean, especially over the Cocos plate.

Geochemistry, Geophysics, Geosystems
Finally, the modeled signal over Barracuda and Tiburon ridges and Barracuda Trough is not as intense as in the original data.
The residuals of the VGG calculated based on the starting model SM6 are shown in Figure 5c. Here lowdensity crustal domains such as the continental platforms and terrains with continental affinity composition (e.g., NNR) are characterized by negative residuals; however, these domains are interbedded with some medium to small scale bodies with positive signal, which become especially important along the Leeward and Lesser Antilles forearcs. In Figure 5c it is also possible to identify a wide body with negative residuals, which appears on the oceanic crust of the South American plate, close to the Lesser Antilles subduction. This body is located in the northern portion of the polygon which denotes the location of the thick BAP. On the other hand, in the location of the tectonic boundary between North and South American plates in the Lesser Antilles subduction, published by Bird (2003) and Case and Holcombe (1980), it is not possible to observe any density contrast, or at least not any one that can be traceable by the VGG. Nevertheless, it is possible to recognize an eastwest oriented strip with strong positive and negative anomalies, in the region of the BFZ, which can be traced until it merges with the Lesser Antilles subduction zone.
The Beata Ridge has associated strong negative residuals, especially concentrated toward its western flank. They define an elongated and irregular body in a northeastsouthwest direction on its northern portion, but that slowly curves toward the southeast, south of ≈15°W. The strongest and widest negative anomalies are located north of 14°W, but it is possible to follow weaker and narrower northwest-southeast residuals, which connect Beata with the SCDB in the Aruba Gap region.
Furthermore, the signal of the residuals over the GB suggests two different crustal domains. The northern portion, which corresponds to the MS crust, is dominated by negative residuals, while the oceanic floor of the GB is controlled by positive ones. Similarly, negative residuals dominate the volcanics of the YB (Central Seamounts province and Cayman Rise and Ridge), while positive values are present on the subbasins. Finally, the spatial behavior of the VGG residuals is not homogeneous along the CT, where negative values are present at ≈2°around the axis of the spreading centre. Toward the eastern portion of GM, close to Hispaniola Island, the residuals alternate positive and negative values. As can be observed in Figure 5c, the residuals change their sign from the continental platforms (negative) to the oceanic basins (positive), suggesting the location of the COT.

Inferred Average Density Field of the Crystalline Crust
The inferred density distribution for the crystalline crust is shown in Figure 6, and is the result of the iterative process described in section 3.4. This average crustal density has been tested in the 3-D lithospheric starting model named SM-F (Table S1). The modeled VGG based on this starting model has a better performance compared with EIGEN-6C4 data set (Table S2), with the minimum values of standard deviation and rootmean-squared error of the associated residuals (9.79 × 10 −9 /s 2 and 23.94 × 10 −9 /s 2 , respectively). Figure  S2 shows the spatial distribution of the residuals from the different starting models. In general, the residuals

Contributions to the Identification of Terrain Boundaries
Considering the potential of the VGG residuals for sensing density contrasts in the crystalline crust, we evaluate their signal over regions where poorly understood tectonic boundaries exist, as well as their performance in the definition of the COT, and in the identification of crustal bodies with anomalous densities.

Boundary of the North and South American Plates in the Lesser Antilles
The tectonic boundary at crustal level between North and South American plates, in front of the Lesser Antilles arc, is still a topic of debate (Pichot et al., 2012;Van Benthem et al., 2013). The Fifteen Twenty Fault Zone (Figure 1) is the well-known boundary, which can be traced from the Middle Atlantic Ridge (MAR), but becomes diffuse near the Antilles subduction zone. Patriat et al. (2011) and Pichot et al. (2012) suggested that the relative motion between these two plates must be accommodated, toward the west, by a heterogeneous and broad fracture zone, which consists of north-south compressional structures, including Tiburon and Barracuda ridges. Figure 7 shows a zoom in the residuals obtained in the NLA region. Figure 7b includes the main fracture zones, as well as the ages of the crust at the subduction zone (Cogné & Humler, 2004) The crustal ages range from 83 Ma in the north of the Lesser Antilles subduction, to 105 Ma in the south. The northernmost fault associated with the magnetic anomalies draws the attention, because of its almost perfect correspondence with the negative residuals at the Barracuda Through, north of BFZ. The crust is 83 Ma old north of this fault and 91 Ma old toward the south. This "jump" in crustal ages is the most pronounced (8 Ma) along the Lesser Antilles, suggesting not only a differential subduction rate, but also a possible decoupling zone between both plates.

Boundary Between the Colombian and Venezuelan Basins
The possible tectonic boundary between the Colombian and Venezuelan basins around the Beata Ridge area has been heavily debated (Dewey & Pindell, 1985;Leroy & Mauffret, 1996;Mattioli et al., 2014;Symithe et al., 2015). The west flank of Beata is a highly deformed region, which is bounded by a large westward-dipping normal fault, where the basement offset reaches up to ≈3,750 m (i.e., see Figure 10 in . According to Driscoll and Diebold (1999), the compressional and extensional-unloading processes that Beata has been subjected to, are recorded in the uplifted and tilted crustal fragments present on its western flank, where the flexural rebound of the lithosphere must have been a relevant factor on its formation. In Figure 5c, negative residuals are concentrated along the western flank of Beata. The inferred average density structure on this part of the ridge ( Figure 6) includes low-density crustal fragments (<2,850 km/m 3 ), which in this specific case, can be related with a prevalence of basaltic and/or gabbroic material in the crust, considering the composition previously reported by Mauffret and Leroy (1997). This dominance of lighter components in the crustal structure may have played an important role in the tectonic uplifting that this ridge has experienced Mauffret & Leroy, 1997). Considering that the negative residuals are located in the most uplifted, deformed, and faulted zone of the ridge, we suggest that they are delineating a region of transition between different crustal domains. In this case, they might be indicating a possible decoupling zone between the Colombian and Venezuelan basins, supporting the tectonic model of Leroy and Mauffret (1996; see orange line in Figure 1). In this area, the negative residuals can be followed from southern Hispaniola, until they merge with the continental margin in the SCDB, along the Aruba Gap.
It is important to mention that the observed signal of the residuals over this gap is not present either in the EIGEN-6C4 data set or in the modeled gradients; therefore, only the analysis of the VGG residuals makes this interpretation possible. Nevertheless, the role of this boundary as an actual kinematic feature needs to be proved considering additional geophysical evidence. The integration and acquisition of new GPS kinematic measurements and seismic data are important pieces in this tectonic "puzzle." For example, contrary to what Leroy and Mauffret (1996) suggested, the seismic results found by Driscoll and Diebold (1999) support the one-plate model, arguing that the majority of the deformation occurred during the early history of the Caribbean plate prior to the sedimentation in this region. Figure 7 provides a closer look at the VGG residuals in the GB, where two different domains can be identified. The northern portion of the basin, which corresponds to the MS crust, is controlled by negative residuals. In this region, the inferred crustal densities (Figure 6) range between 2,700 and 2,850 km/m 3 , suggesting a light crust probably affected by the Aves Ridge and/or the Lesser Antilles magmatism, which might have continued in this region during the slab roll-back (Arnaiz-Rodríguez & Audermard, 2018). To the south of the basin, the residuals are positive and the inferred densities (2,900 to 3,050 km/m 3 ) are in agreement with oceanic crust formed during the extensional process that opened the basin (Aitken et al., 2011); at this location, the highest density areas might suggest the presence of localized underplated material. The differences in the extensional processes that affected the evolution of the currently named GB let us consider that the characterization of its crustal structure may be improved, and that the definition of a boundary between these two sub-basin domains should be considered.

Continental-Oceanic Transition (COT)
Considering that one of the strongest density contrasts might occur along the continental margins, the sign change of the residuals of the VGG along these regions (Figure 5c) is interpreted as the COT. However, it is possible to recognize that in most of the cases (e.g., southern margin of the Caribbean plate and Muertos Thrust Belt) such change does not coincide with the COT boundaries defined by Bird (2003) and Case and Holcombe (1980). Furthermore, the inferred crustal densities of continental-like material along the continental platforms are in general lower than 2,850 km/m 3 ; however, some high-density intrusions (>2,900 km/m 3 ) can be identified, especially along the Leeward and Lesser Antilles forearcs. These intrusions are discussed in section 5.2.1. Lucia, positive residuals dominate this axis and can be followed until they merge with the continental platform of Venezuela. Additionally, a generalized trenchward pattern that alternates two along-arc stripes of positive and negative residuals is present in the forearc, from eastern Hispaniola to offshore Grenada Island.

Identification of Crustal
In Figure 7, a zoom in the Lesser Antilles shows more details about the structure of these residuals. Figure 7b includes a compilation of the backstop locations (limit of the overriding crystalline crust) in the NLA, inferred by Laigle, Becel, et al. (2013) and Evain et al. (2013), plotted on top of the VGG residuals. The importance of characterizing the backstop position in subduction systems lies on the fact that they are a widely used proxy for the up-dip limit of the seismogenic zone (i.e., Laigle, Becel, et al., 2013), and because they could play an important role in the stress field between the overriding and down-going plates, affecting the structure of the respective forearc basins (Noda, 2016). The backstop locations identified in the Lesser Antilles partially overlap with a region of transition from negative to positive residuals (transition from overriding to down-going plates, respectively), especially in front of Guadeloupe and Dominica Islands, where the outer forearc acts as backstop .
It is important to consider that north of Guadeloupe Island, the location of the backstop identified by Laigle, Becel, et al. (2013) is not accurate, because the sediments of the forearc basin (normally used as a proxy for the location of the backstop) are highly deformed, making their differentiation from the accretionary prism difficult. Exactly in this region, the negative residuals do not fit the spatial location of the backstop defined by these authors. Toward the south of Dominica, the outer forearc does not act as a backstop ; however, the negative residuals continue seaward, defining an extensive (≈31 x 10 3 km 2 ) anomalous low-density body (<2,850 km/m 3 ; Figure 6) east of 60°W. This body is discussed in section 5.2.2. Evain et al. (2013) proposed two hypotheses in order to explain the origin of the forearc structure in front of Guadeloupe, Dominica and Martinique Islands: (1) the outer forearc is an eastward extension of the inner forearc, which has been subjected to processes of alteration and erosion due to the subduction of the Atlantic crust, and (2) the outer forearc represents an accreted allochthonous body, especifically, a fragment of a volcanic ridge. Based on the spatial distribution of the negative gravity residuals in the North and South Lesser Antilles (Figure 5c), we cannot rule out a composite origin and evolution for this forearc, which includes both hypotheses. Nevertheless, considering that a structure similar to that of the outer forearc (outermost negative residuals) is present and continuous from the northeast of Hispaniola to Martinique Islands, it is unlikely that allochthonous terrains alone could have created this low-density forearc of the Antilles.
Similarly, we interpret the anomalous high-density bodies located in the forearc of the Leeward Antilles, in the SCDB, as fragments of either underplated and/or mafic intrusions, or fragments of picrites accreted to the continental margin, such as those reported by Mauffret and Leroy (1997) in Curaçao Island. It is important to note that this interpretation of the crustal structure along the Antilles forearcs only can be done based on the residuals, because they highlight the signal of crustal features that have not been included on the starting model. Nevertheless, the results over the Lesser Antilles might also be affected by the Moho error ( Figure S4), which ranges from 1 to 7 km in this specific region.

Anomalous Low-Density Body in the Atlantic Ocean
Based on the low density of the anomalous body present in the Atlantic Ocean, in front of Dominique and Martinique Islands, we suggest that its origin might be associated with a volcanic arc of intermediate composition. Similar findings and hypotheses have been previously described by other authors in the Lesser Antilles (i.e.: Bangs et al., 2003;Bouysse & Westercamp, 1990;Christeson et al., 2003). This body is present in all the residuals of the different starting models, as can be seen in Figure S2. However, because it is partially buried under the thick layer of sediments of the BAP, the location of its edges depends on the density distribution used for the modeling of the sedimentary layer.
Moreover, Figure 7c shows the seismicity from 1980 to 2016 (International Seismological Centre, 2019) in the Lesser Antilles, plotted on top of the VGG residuals. In a regional scope, this body seems to be the boundary between two different clusters of seismicity in the Lesser Antilles: a clear Wadati-Benioff zone toward the north, and a more diffuse structure toward the south. The oceanic crust is ≈92 Ma old over the anomalous body, and it is partially delimited toward the north and south by two fracture zones (Figure 7b). The northern fracture zone can be followed from TFZ until it merges with the Lesser Antilles forearc. Southward of the anomalous body, the  Rosencrantz (1994;in Mann, 2007)): A = 0-10 Ma. B1 and B2 = 10-20 Ma. C1 and C2 = >20 Ma. Terrain acronyms as depicted in Figure 2.
crust becomes 98 Ma old; thus, the second biggest "jump" (6 Ma) in the crustal ages along the trench occurs. Furthermore, this body partially overlaps with the slab gap found at ≈400-km depth, reported by Van Benthem et al. (2013) between 13 and 15°N, and that has been interpreted as the North American -South American plate boundary at these mantle depths. Therefore, considering the evidence that supports a northward migration of this plate boundary (probably from TFZ to BFZ -Pichot et al., 2012), a possible interpretation is that the slab gap at ≈400-km depth represents the paleo-boundary between both plates, which nowadays has moved toward the north. In this hypothetical scenario, the anomalous buoyant body could have played an important role, locking the subduction and favoring the northward migration of transpression.

Yucatan Basin and Cayman Trough
The pattern of the residuals over the YB indicates at least two crustal domains below the basin (Figures 8a  and 8b). In the north and northwest, positive residuals dominate the oceanic sub-basins. On the other hand, the southern portion is characterized by negative residuals, connecting the volcanics of the Central Seamounts province, the Cayman Rise and Cayman Ridge. Rosencrantz (1990) proposed a common origin for these volcanic features; therefore, a similar crustal density distribution could be expected for these terrains in YB. In our case, the inferred crustal densities for these domains range from 2,750 to 2,950 kg/m 3 , with the lightest crust concentrated over the western half portion of Cayman Ridge, where Perfit and Heezen (1978) described a crust formed by Cretaceous amphibolites and Late Cretaceous to Paleocene plutonics, whose densities are in the range of our estimates (Schön, 2011).
On the other hand, the crust of the CT shows a transition from oceanic to highly extended continental crust (from the axis of the spreading centre to the edges of the rift), with serpentinized sections in the middle of the trough (ten Brink et al. (2002)), as was also reported by Perfit and Heezen (1978). A serpentinized mantle can change the behavior of the lower crust decreasing its density. In fact, the density of 40-50% serpentinized rocks ranges between 2,800 and 2,900 kg/m 3 (ten Brink et al., 2002). The inferred crustal densities in CT ( Figure 8c) include a region with relative low-density crust within approximately 2°from each side of the ridge, which partially overlaps the serpentinized crust identified by ten Brink et al. (2002). This low-density crust connects the Oriente Fault Zone with the central part of the ridge, close to the limit of the 10-to 20-Maold crust defined by Rosencrantz (1994) (see regions B 1 and B 2 in Figure 8c). The densities there range from

10.1029/2019GC008340
Geochemistry, Geophysics, Geosystems 2,750 to 2,900 kg/m 3 , the latter being the dominant value. Thus, our results might be indicating the extension of the serpentinized material (where it dominates the crustal structure), as a result of the hydrothermal circulation along the different faults that bound this trough (ten Brink et al., 2002). However, we could not assure that other portions of crust (or upper mantle) in CT have not been affected by serpentinization.
Moreover, the outer edges of CT do not have a common crustal density structure, as can be seen from the gradient residuals and the corresponding inferred densities. The eastern edge, which belongs to the GM, shows a clear transition from dense (≈3,050 kg/m 3 ) and more homogeneous crust -with oceanic affinity-(see region C 2 in Figure 8c), to a more heterogeneous, extended-continental crust (2,750 to 3,050 kg/m 3 ), which extends toward Hispaniola Island. The boundary between both domains approximately corresponds with the >20-Ma-old crust limit reported by Rosencrantz (1994;in Mann, 2007). Considering the tectonic evolution proposed for this area, the CT migrated from a location near the continental platform of Belize toward the northeast, until its current location (Leroy et al., 2000). This tectonic model implies that the east side of the trough should have moved in a northeast direction, experiencing more extension compared with the west side. Our results indicate that the western edge of CT is dominated by a more homogeneous density distribution ( Figure 8c). Thus, our inferred crustal densities are in agreement with the tectonic model discussed by Leroy et al. (2000).
Taking into account the previous discussions and the compiled terrain boundaries of the Caribbean domain published by Case and Holcombe (1980), we propose or confirm major tectonic or terrain boundaries in the study area, in the regions where the density contrast is high enough to be highlighted by VGG residuals. These boundaries include the following: 1. the limit between the North and South American plates, in the Lesser Antilles subduction; 2. the possible boundary between the Colombian and Venezuelan basins; 3. the location of the edge of the MS crust in the GB; and finally, 4. the COT.
Additionally, we identify widespread high-density (inner) and low-density (outer) forearc domains in the Lesser and Leeward Antilles; an anomalous low-density body in the Atlantic Ocean; the edges of the volcanic provinces inside the YB; and regions of low-density material, possibly related with serpentinized crustal fragments in the CT. These boundaries are plotted on top of the VGG residuals in Figure 9. It is important to mention that the limits of the CT and North and South Nicaraguan rises remain the same as in Case and Holcombe (1980), because they have been well recognized along main fault zones (see Figure 2).
The location of the COT along the north-western side of the Panama Deformation Belt and along the north of the South Lesser Antilles remains unclear based on the VGG residuals (question marks and magenta lines in Figure 9). However, the negative residuals in the north SLA indicate the presence of low-density crust which origin might be related with volcanic-arc processes. Figure S5 includes the boundaries based on this study plotted on top of the free-air anomalies.

Summary and Conclusions
The VGG are gravity derivatives especially sensitive to high-density contrasts in the upper layers of the lithosphere. Using the VGG from EIGEN-6C4 data set in the Caribbean region, we evaluated the results of six different lithospheric density models and their performance in comparison with observations. The model which includes 3-D density solutions for water and sediments (starting model SM6; Table S1) showed the best performance (Table S2), even though an a priori constant density was used for the crystalline crust. Based on the VGG residuals (equation (8)) of this lithospheric model, we refined, proposed and/or confirmed tectonic (or terrain) boundaries in the Caribbean oceanic region. Additionally, the forward modeling of the VGG allowed us to infer an average density field for the crystalline crust, which provides information about the dominant density out of the different layers that compose the crystalline crust. Supplementary data sets with the results and the scripts to reproduce the workflow are available at Gómez-García et al. (2019a) and Gómez-García et al. (2019b).
Based on our approach, the main conclusions that we would like to highlight are: 1. The analysis of the VGG residuals allowed us to recognize the COT in the study area. The identification of this boundary is important, for example, in the modeling of petroleum systems, but can also add valuable information in plate tectonic reconstructions.

10.1029/2019GC008340
Geochemistry, Geophysics, Geosystems 2. The residuals suggest the existence of an anomalous buoyant body (AB; Figure 9) of ≈31 × 10 3 km 2 in the oceanic crust of the South American plate, in front of Dominique and Martinique islands. Its inferred density is <2,850 kg/m 3 , and considering the tectonic history of the Lesser Antilles subduction, this body is most likely related with a fragment of a volcanic arc of intermediate composition. Even though this body seems to be the boundary between two different clusters of seismicity in the Lesser Antilles (Figure 7c), its existence needs to be proved using more direct geophysical techniques (i.e.: drill cores, seismic reflection), able to overcome the thick sedimentary layer (up to 18 km) of the BAP, which overlies the crystalline crust in this area. 3. Based on the diverse geophysical information collected in this study, we propose that the current main boundary between North and South American plates, in the Lesser Antilles, is located at the Barracuda Trough. At that location, strong negative residuals can be followed until the subduction, whose spatial distribution coincides with a transform fault. Additionally, a "jump" of 8 Ma in the crustal ages suggests a decoupling zone (Figure 7b). 4. The boundary of the inner and outer forearc domains, previously identified in the north central portion of the Lesser Antilles by Evain et al. (2013) and Laigle, Becel, et al. (2013), coincides with regions where there is a change in the residuals sign ( Figure 7b). The inner forearc is characterized by positive residuals and high crustal densities (Figure 6), as was expected from the high-seismic velocities reported by Evain et al. (2013). Similarly, the outer forearc is associated with negative residuals and lower crustal densities, supporting the low-seismic velocities reported by these authors. Such interpretation was extrapolated to the entire forearc domain of the North and South Lesser Antilles, as well as to the Leeward Antilles, where clear high and low-density bodies (positive and negative residuals, respectively) are present ( Figure 9). The location of these forearc domains could be considered as up-dip limits in the modeling of the seismogenic zones in these subduction systems. 5. The boundary between the Colombian and Venezuelan basins corresponds with negative VGG residuals, that can be followed from the Aruba Gap (Pecos Fault) to the south of Hispaniola Island, toward the west flank of Beata Ridge (Figure 9). Spatially, these residuals coincide with a westward dipping fault . We suggest that this boundary should be considered in GPS based models in order to evaluate the two-plate hypothesis proposed for the Caribbean plate. 6. The residuals of the VGG highlight the boundary between two different crustal domains in GB. The low-density volcanic material of the MS crust dominates the northernmost third of the basin. On the other hand, the positive residuals dominate the rest of the basin, which according to the inferred densities is made of oceanic crust with few fragments of volcanic material, probably associated with the extension process of Aves Ridge. Our findings support the hypothesis of Arnaiz-Rodríguez and Audermard (2018), whose interpretations include a differential roll-back speed in the NLA subduction during the opening of the GB, allowing the flow of molten material in the northern portion of the basin. 7. The inferred average density structure of the crystalline crust ( Figure 6) is in agreement with the previously known geology of the different Caribbean provinces. For example, the low-density material of the continental platforms has an inferred density ranging from 2,700 to 2,800 kg/m 3 . Similarly, the density structure that characterizes the NNR is related with its continental and island arc fragments; as well as the structure of Aves Ridge, which can be associated with an island arc of intermediate composition. Finally, the residuals of the VGG and the low-density fragments that can be inferred inside the Cayman Through let us to propose the possible extension of serpentinized crustal fragments, besides the spatial distribution of the volcanic provinces in YB (Figure 8).