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