Lithosphere Removal in the Sierra Nevada de Santa Marta, Colombia

The Sierra Nevada de Santa Marta (SNSM) in northwestern Colombia is one of the world's highest coastal mountains, with an elevation above 5 km. Gravity measurements show that the SNSM has a positive Bouguer anomaly (>+130 mGal), indicating that the mountain lacks a crustal root. In this work, we test the hypothesis that these observations can be explained by gravitational removal of the dense lower lithosphere. We use 2D numerical models to examine the dynamics of lithosphere removal and its effect on surface elevation and gravity anomaly. The models consist of continental lithosphere that includes a pre‐thickened crustal region, representing the SNSM. In our preferred model, the dense mantle lithosphere and crustal root are gravitationally unstable and undergo removal as local drips within ∼10 Ma from the onset of foundering. This creates an area of thinned crust (∼38 km) underlain by a buoyant sublithospheric mantle where melting and low seismic velocities are predicted. Subsequent non‐isostatic forces maintain a topography of 3.3 km with a 2D Bouguer gravity anomaly of +103 mGal. Parameter tests show that a strong lower‐crustal rheology provides greater support for the high topography and that a weak mantle lithosphere rheology produces faster removal. The models demonstrate that local lithosphere dynamics can explain the first‐order observations in the SNSM. We propose that lithosphere removal could have occurred recently (∼2 Ma), explaining the localized low seismic velocity zone below the SNSM, or at 56–40 Ma, inducing anomalous short‐lived Eocene magmatism.

gravitational foundering of the deep lithosphere as a local drip can simultaneously produce lithospheric and/ or crustal thinning, surface uplift, and symmetric high topography (e.g., Göğüş & Pysklywec, 2008;H. Wang & Currie, 2017).This mechanism has been proposed to explain localized changes in topography and trends in magmatism in various mountain ranges, such as the Sierra Nevada in California, the Puna Plateau of South America, the Wallowa Mountains in northeast Oregon and the central Anatolian plateau (e.g., DeCelles et al., 2015;Göğüş et al., 2017;Hales et al., 2005;Saleeby et al., 2012).In some cases, removal has been associated with localized Bouguer gravity anomalies of 10's of mGals (e.g., H. Wang et al., 2015;Saleeby et al., 2012, and references therein).Removal of the lower lithosphere can be produced by the growth of perturbations in the cold, dense mantle lithosphere (Houseman et al., 1981) or by gravitational instability of localized magmatic/ metamorphic eclogite (Jull & Kelemen, 2001;H. Wang et al., 2015).The SNSM region had a thick continental crust (∼64 km) in the Jurassic (Ramírez et al., 2020) and was subjected to compression during the collision of the Caribbean plateau in the Late Cretaceous-Early Paleocene (Cardona et al., 2010).Therefore, the associated crustal thickening could have fostered metamorphic eclogitization in the deep crust, while simultaneously producing a mantle lithosphere instability.Alternatively, previous stages of arc magmatism in the SNSM region (e.g., Cardona et al., 2010;Duque-Trujillo et al., 2019;Ramírez et al., 2020) could have resulted in the formation of dense residues (also called arclogites) after partial melting or fractionation in the lower crust (Ducea, Chapman, Bowman, & Triantafyllou, 2021).The subsequent removal of the dense lower crust would explain the absence of a thick crustal root in the SNSM.
In this study, we use 2D numerical models to investigate if the observations of local high surface topography and positive Bouguer gravity anomaly can be explained by drip-style lithosphere removal.We assume that the SNSM topography was produced through past tectonic processes and we investigate under which conditions this high topography can be maintained throughout a lithosphere drip.Our models explore the effect of different lithosphere rheologies on removal dynamics and the surface expressions.We also examine the effect of denser upper crustal rocks (e.g., Sanchez-Rojas & Palma, 2014) on the gravity anomaly, and the role of the adjacent vertical strike-slip faults on removal dynamics.

Tectonic/Geologic Background
The SNSM is part of a larger crustal fragment known as the Maracaibo block and is located at its northwestern corner, approximately 100 km southeast of the tectonic margin between the Caribbean and South American .BGB: Baja Guajira Basin, OF: Oca Fault, SNSM: Sierra Nevada de Santa Marta, LMB: Lower Magdalena Basin, CRB: Cesar-Rancheria Basin, PR: Perijá Range, SMBF: Santa Marta-Bucaramanga Fault.The red barbed line is the plate margin where the Caribbean plate (CAR) converges with South America.The cyan, magenta and green contours delimit respectively the three geotectonic provinces: (1) Sierra Nevada, (2) Sevilla and (3) Santa Marta.(b) Bouguer gravity anomaly map from the EIGEN-GL04C Global gravity field model (Förste et al., 2008).Profiles of topography and gravity anomaly (A-A′ and B-B′) are used to compare the models with observations.10.1029/2023JB027646 3 of 23 plates (Figure 1a).At present, the Caribbean plate is obliquely converging with South America at a rate of 10-20 mm•yr −1 toward the southeast (e.g., Freymueller et al., 1993).Regional seismological studies show that the Caribbean plate subducts (or underthrusts) the South American plate at a shallow angle (<30°) and the slab is at a depth of 80-120 km below the SNSM (e.g., Cornthwaite et al., 2021;Londoño et al., 2020;Taboada et al., 2000;Van Der Hilst & Mann, 1994;Vargas, 2020).
The triangular SNSM massif is bounded by major faults and basins on its three sides (Figure 1a).At the north, the SNSM is limited by the right lateral Oca fault and the Baja Guajira basin.At the southwest, it is bounded by the left-lateral Santa Marta-Bucaramanga fault and the Lower Magdalena basin.At the southeast, it is bounded by the Cesar-Rancheria basin and the Perijá range.The Santa Marta-Bucaramanga and Oca faults exceed 200 km in length and have experienced lateral displacements since the Paleocene of 100 and 65 km, respectively (e.g., Case & Macdonald, 1973;Tschanz et al., 1974).The Oca fault is not currently active but the Santa Marta-Bucaramanga fault system has predominantly left-lateral active kinematics (Idárraga-García & Romero, 2010).
The SNSM is commonly divided into three geo-tectonic provinces (Tschanz et al., 1974): the Sierra Nevada, Sevilla, and Santa Marta provinces (Figure 1a).According to Ramírez et al. (2020) the oldest rocks in the SNSM are in the Sevilla province (the Buriticá and Los Muchachitos gneisses) and Sierra Nevada province (Los Mangos granulite and Dibulla gneiss).These form the Precambrian metamorphic basement (e.g., Sanchez & Mann, 2015).In the Sierra Nevada province, the basement is overlain by Devonian-Carboniferous and late Triassic-Early Jurassic sedimentary rocks (e.g., Cuchilla de Carbonal and Los Indios formation, respectively).These rocks are mostly covered by outcropping Jurassic igneous bodies (Tschanz et al., 1974), which have intermediate to felsic compositions and mafic intrusions (Ramírez et al., 2020).Quandt et al. (2018) show that these are continental arc plutons associated with a Jurassic magmatic arc.In the Sevilla Province the basement is mostly overlain by Paleozoic metamorphic rocks (e.g., El Encanto orthogneiss, the Sevilla Metamorphic complex).The Santa Marta province is formed by Cretaceous metamorphic rocks (the Santa Marta metamorphic belt) which are intruded by Paleogene arc-type igneous rocks (the Santa Marta Batholith).This region was influenced by the collision of the Caribbean large igneous province starting at 75 Ma (Spikings et al., 2015), which induced compressive stresses, potentially producing crustal thickening and uplift in northwestern South America.The diachronous collision progressed toward the north during the westward displacement of the South American plate relative to the Caribbean plate (e.g., Bayona et al., 2011;Pindell et al., 2005).The collision was followed by the onset of the shallow subduction of the Caribbean plate below South America in the early Paleocene (Villagómez et al., 2011), which produced magmatism in the NW corner of the SNSM, forming the Santa Marta batholith and other arc-type plutons (Cardona et al., 2010).Therefore, the rocks of the Santa Marta province can be used to infer the tectonic evolution of northwest South America after the Late Cretaceous.According to the accepted model, during the Cenozoic the massif experienced a clockwise rotation that displaced it from the Colombian Central Cordillera to its current location.This deformation resulted in its characteristic triangular shape and its isolated location (Montes et al., 2010;Ujueta, 2003).

Governing Equations and Modeling Approach
In this study we use 2D numerical models to examine whether lithosphere removal can explain the observed local positive Bouguer gravity anomaly in the elevated SNSM region.The models use the ASPECT code version 2.2.0 (Bangerth, Dannberg, Gassmoeller, Heister, et al., 2020;Heister et al., 2017).The code solves the governing equations of incompressible mass conservation, momentum conservation, and heat conservation, using the Boussinesq approximation formulation and assuming plane strain (Kronbichler et al., 2012).We model the thermal-mechanical evolution of the lithosphere and upper mantle, where the thermal and mechanical fields are coupled through temperature-dependent material properties.
We compare the model results to two key surface observables: surface topography and Bouguer gravity anomaly.The topography is determined by using a stress-free upper boundary that allows dynamic surface deformation throughout the model evolution (Rose et al., 2017).The Bouguer gravity anomaly calculation at each timestep, for each horizontal position along the profile, is based on the density differences (anomalies) with respect to the reference (non-perturbed) density structure in Figure 2b corresponding to a crustal thickness of 40 km and a total lithosphere thickness of 100 km.The contribution to the total gravity anomaly of the density difference in each element in the domain is approximated as the anomaly of a finite slab (Telford et al., 1990;details in Quiroga, 2022).

Initial Setup and Material Properties
The 2D model domain (Figure 2) has a width of 1,320 km and height of 660 km, representing the continental lithosphere and the upper mantle in the SNSM region.Our model does not include the subducting Caribbean plate because the goal is to model a lithospheric drip and its isolated contribution to surface observables.The mesh has square elements with lengths and widths of 2.58 km in the upper 200 and 10 km below.The upper region has a finer mesh to properly resolve lithospheric deformation.The side boundary conditions are free slip, the bottom boundary condition is no slip, assuming that the upper mantle is coupled to the high-viscosity lower mantle, and the top boundary has a free surface.
The continental crust is 40 km thick, with a 20 km upper crust and 20 km lower crust, underlain by a 60 km mantle lithosphere.An area of initially thick lithosphere is placed in the middle of the model (Figure 2a) to represent the SNSM.This is based on the results by Ramírez et al. (2020), who report a continental crust of about 64 km during the Jurassic arc, and also considering that the subsequent collision of the Caribbean plateau in the Late Cretaceous-Early Paleocene is inferred to have produced further thickening (Villagómez et al., 2011).In this work, we do not model shortening, but rather we look at the dynamic response after crustal thickening.In the model, the thickened region consists of an initial surface topography and crustal root, where the deflections have Gaussian shapes with half-widths of 62.5 km, resembling the width of the SNSM.At its maximum, the crust in this region is 35 km thicker than the undisturbed crust.The initial topography is set to be in isostatic balance under the assumption that equilibrium occurred right after or during shortening.The crustal thickening results in a perturbation of the lithosphere-asthenosphere boundary (LAB).The zero topography corresponds to the unperturbed surface (upper boundary) at the sides of the thickened region in the modeling setup (Figure 2a).
All materials have temperature-dependent densities and viscous-plastic rheologies using the material properties given in Table 1.The reference model (model A) is based on the setup in Figure 2a.The upper crust, lower crust, mantle lithosphere and sublithospheric mantle have viscous rheologies of wet quartzite (Gleason & Tullis, 1995), dry Maryland diabase (Mackwell et al., 1998), dry olivine and wet olivine (Karato & Wu, 1993), respectively.Frictional-plastic deformation includes linear strain weakening, where all materials have a reduction in the friction angle and cohesion of 15°-2° and 20-2 MPa, respectively, for cumulative plastic strain between 0.5 and 1.5.Other models (sets B, C, and D) examine variations in crustal density and/or lithosphere rheology, and also include the effect of lower crustal eclogitization.Table 2 contains the list of models and their respective modifications relative to model A.
The initial thermal structure is based on a conductive geotherm in the lithosphere, calculated using a surface heat flow of 54 mW•m −2 , thermal conductivity of 2.25 W•m −1 •K −1 , and radiogenic heat production in the upper and lower crust of 1 and 0.4 μW•m −3 , respectively; there is no heat production in the mantle (Table 1).This geotherm  1 for references.UC = upper crust, LC = lower crust, ML = mantle lithosphere, SLM = sublithospheric mantle, and LAB = lithosphere-asthenosphere boundary.
(Figure 2c) has an initial Moho temperature of 650°C and intersects the mantle adiabat at the LAB (100 km depth) for the unperturbed lithosphere.Below this, the mantle is adiabatic, with a potential temperature of 1292°C and adiabatic gradient of 0.4°C•km −1 .In the central region, the continental geotherm is linearly stretched following the Gaussian-shaped thickening, such that the initial Moho temperature is constant everywhere.During the model run, the model boundaries have fixed temperature conditions using the initial temperature.

Crustal Eclogitization
Eclogitization is considered to affect the lower crust in areas of thickened crust and can result in a density increase that destabilizes the crustal root, leading to its detachment (e.g., Leech, 2001).Metamorphic eclogitization can occur at temperatures and pressures greater than 600-800°C and 1.0-2.0GPa, respectively, and the phase change is promoted by the presence of water (e.g., Ahrens & Schubert, 1975;Austrheim et al., 1997;Leech, 2001).In addition, in continental arcs, eclogite facies rocks can form as dense residues after magmatic differentiation in the crustal root; these are sometimes called arclogites (e.g., Ducea, Chapman, Bowman, & Triantafyllou, 2021).It is possible that either metamorphic or magmatic eclogitization affected the SNSM, as it is an area of previous thickened crust, and because it originally formed as a continental magmatic arc in the Colombian Central Cordillera (e.g., Cardona et al., 2010;Ramírez et al., 2020).
Model A does not include eclogitization.In model sets B, C, and D, eclogitization in the lower crust occurs where its temperature is above 680°C and its pressure is above 1.2 GPa.The eclogite has a density of 3,550 kg•m −3 , representing full eclogitization of a basaltic lower crust (Austrheim et al., 1997).The rheology of eclogitized lower crust is poorly constrained.Laboratory measurements using a dry synthetic eclogite suggest relatively high strengths (e.g., Zhang & Green, 2007).On the other hand, eclogitization may be accompanied by rheological  The density (ρ) is given by: ρ where α is the thermal expansion coefficient, T is the temperature (in °C) and T 0 is the surface temperature (0°C).b The frictional-plastic yield stress (σ y ) is based on a Drucker Prager yield mechanism: σ y = P sin(ϕ) + C cos(ϕ), where ϕ is the angle of internal friction, and C is the cohesion.c Where stresses are below the plastic yield stress viscous deformation follows a dislocation creep power law, with an effective viscosity (η eff ) given by: k , where A ps is the pre-exponential factor, n is the stress exponent,   ′  is the second invariant of the deviatoric strain rate tensor, E is the activation energy, V is the activation volume, P is the pressure T k is the absolute temperature, and R is the gas constant (8.3145J mol −1 K −1 ).d The experimental uniaxial strain viscosity pre-factor A uni is scaled to plane strain A ps using a scaling factor of  3 +1 2 2 −1 .e k is the thermal conductivity, C p is the specific heat capacity, and Q p is the internal radiogenic heat production.

Table 1 Model Parameters for the Reference Model (Model A)
10.1029/2023JB027646 6 of 23 weakening.Eclogite strength is conditioned by the extent and rate of the reaction, the amount of water available to trigger the phase change or the mechanisms by which water enters the rock (e.g., Austrheim et al., 1997;Leech, 2001).This suggests that eclogites could be weaker than their protoliths.In our models, we use the viscous rheology of eclogite of Zhang and Green (2007) with a scaling factor (f) of 0.1 which reduces its viscosity by a factor of 10 (Figure 2d), to represent a rheological change relative to the initial rheology due to the phase change (weakening relative to diabase and strengthening relative to granulite).Additional model tests show that variations in eclogite rheology produce changes in the timing of crustal root removal; weaker eclogite rheologies produce a faster root removal (and vice-versa).However, these variations do not affect the magnitudes and general behavior of surface expressions (Figures S1 and S2 in Supporting Information S1).

Reference Model (Model A)
The results of model A are shown in Figures 3 and 4 and in Movie S1.In this model, the imposed perturbation of the lithosphere-asthenosphere boundary (LAB) causes the lower part of the mantle lithosphere to undergo dripping due to its cooler temperature and resultant negative buoyancy (Figure 3).As the instability at the center founders, the mantle lithosphere is entrained from the sides and pulled toward the center.This results in thinning of the mantle lithosphere by 18 km at the sides.The drip detaches at ∼40 Ma, removing the lower 30 km of the mantle lithosphere in the central region.The shallower lithosphere is cool and therefore is too viscous to be removed (e.g., Conrad & Molnar, 1999).The crustal root is too buoyant to be removed but the drip causes a downward displacement of the Moho, thickening the crust by ∼3 km as the drip founders.
The evolution of topography above the perturbation is characterized by a total subsidence of 0.12 km produced by the pull of the dripping mantle lithosphere (Figure 4a).The drip detachment at ∼40 Ma has a minimal effect on topography, and it is followed by minor topographic changes (<0.1 km) related to continued mantle convection.Throughout the model evolution there is a gradual decrease in the Bouguer gravity anomaly of ∼20 mGal (Figure 4b).This is mostly attributed to progressive mantle lithosphere thinning at the sides of the drip as the Set C (Variations in lower-crustal rheology with a wet mantle lithosphere rheology) Set D (Variations in lower-crustal rheology with a wet mantle lithosphere rheology and a higher density upper crust in the SNSM massif) Set E g (Effect of the Santa Marta -Bucaramanga and Oca faults) Note.The mantle lithosphere rheology is varied locally (in the central region) but remains as dry olivine at the sides for all models.a   *  is the reference density in the upper crust in the perturbed region.b DMD = Dry Maryland Diabase (Mackwell et al., 1998).c DO = Dry Olivine (Karato & Wu, 1993).d LMG = Local Mafic Granulite (MG in the central region and DMD at the sides).e MG = Mafic Granulite (Y.F. Wang et al., 2012).f WO=Wet Olivine (Karato & Wu, 1993).g Model E is equivalent to model D1 but includes vertical weak zones representing the Santa Marta-Bucaramanga and Oca faults. 1 cold, dense mantle lithosphere is replaced with hot, low-density sublithospheric mantle.There may also be a minor contribution from the downward deflection of the Moho during the drip.As the drip detaches at 40 Ma, there is an abrupt gravity anomaly increase of 28 mGal produced by the rapid foundering of the cold, dense mantle lithosphere.Its effect on the surface gravity anomaly decreases rapidly as the detached mantle lithosphere sinks into the deeper mantle.

Effect of Lower Crustal Eclogitization (Model B1)
The results of model B1 are shown in Figures 4 and 5a and in Movie S2 (supplementary material).This model considers eclogitization of the deep crust; all other parameters are the same as in model A. The mantle lithosphere perturbation grows, drips, and detaches in a similar manner to model A (Figure 5a).However, the phase of instability growth prior to detachment is accompanied by progressive eclogitization of the lowermost crust.The deep crust is not initially in the eclogite stability field.The phase change conditions are induced by the pull of the foundering mantle lithosphere and the temperature increase during lithosphere removal.The negative buoyancy of both the eclogite and mantle lithosphere perturbation results in detachment of the mantle lithosphere drip at 38 Ma, which is 2 Ma earlier than in model A. Eclogitization continues after drip detachment, and the crustal root starts to founder at ∼55 Ma because of its high density.From 55 to 71 Ma, both the crustal root and underlying mantle lithosphere drip into the deeper mantle.The eclogitized crustal root detaches in two pulses (Animation 2).At 59 Ma there is detachment of the lowermost root.This leaves a ∼4 km thick eclogite layer that fully detaches from the non-eclogitized crust at 71 Ma and sinks to the deep mantle.Hereafter we refer to the time of root removal as the time when the eclogitized crustal root detaches and is no longer in contact with the non-eclogitized crust.Root removal leaves a thin crust (38 km) in the perturbed region (∼2 km thinner than the unperturbed crust), underlain by the upwelling sublithospheric mantle.
The evolution of topography in the perturbed region is influenced by the combined effects of the mantle lithosphere drip, eclogitization, and removal of the crustal root.From 0 to 71 Ma, the negative buoyancy of the growing mantle lithosphere instability, and the progressive eclogitization and foundering of the lower crust produce ∼2.4 km of surface subsidence (Figure 4a).Subsidence is not steady, but includes a period of more rapid subsidence (0.2 mm•yr −1 ) at ∼36 Ma preceding detachment of the mantle lithosphere drip.This is followed by relaxation of the downward pull that results in a subsidence rate of 0.04 mm•yr −1 by 38 Ma.The pulses of crustal root detachment at 59 and 71 Ma are also associated with temporary decreases in subsidence lasting 3 and 1 Ma, respectively.This is explained by the upwelling of hot and buoyant sublithospheric mantle as the crustal root founders, which allows for temporary support of elevated topographies of 2.5 and 2.2 km at 59 and 71 Ma, respectively.After root removal, subsidence continues as the model comes into isostatic equilibrium with the new lithosphere structure, with an elevation of ∼1 km at 100 Ma.
Initially, the central region has a Bouguer gravity anomaly of −169 mGal due to the non-eclogitized crustal root (Figure 4b).As the crustal root undergoes eclogitization, the gravity anomaly increases and has a positive value between 47 and 72 Ma.The mantle lithosphere drip produces a peak of +12 mGal at 38 Ma.At the onset of crustal root dripping, there is an increase in the Bouger gravity anomaly, reaching a maximum of +30 mGal at ∼55 Ma.During detachment of the crustal root between 59 and 71 Ma the gravity anomaly fluctuates between +4 and +32 mGal.After root detachment, the eclogitized root is replaced by buoyant sublithospheric mantle.This model shows that a positive gravity anomaly and a topographic high are simultaneously produced during crustal root removal between 46 and 72 Ma (Figure 4a).This trend is observed in the models discussed below, but the timescales and magnitudes vary with different rheologies and/or density parameters.

Lower Crustal Rheology (Set B)
Lithospheric rheology strongly affects drip dynamics (e.g., Göğüş & Pysklywec, 2008;H. Wang & Currie, 2017).In model set B, we test different lower crustal rheologies (Figure 2d; Table 2).The lower crust in model B1 has the rheology of dry Maryland diabase, as in reference model A. Model B2 uses the lower crustal rheology of mafic granulite (Y.F. Wang et al., 2012) in the thickened region and has dry Maryland diabase at the sides.We refer to this as local mafic granulite, assuming it represents a locally weaker composition associated with the previous magmatic arc.In model B3, the rheology of the whole lower crust is mafic granulite, producing a weak lower crust along the full profile.
The general dynamics of the three models are similar, with an initial mantle lithosphere drip, followed by removal of the eclogitized lower crust (Figure 5).Compared to model B1, the weaker lower crust allows for a slight increase in the removal rate, where detachment occurs at 64 Ma and 60 Ma for models B2 and B3, respectively, compared to 71 Ma for model B1.In models B2 and B3, elevations are consistently lower than in model B1, suggesting less lateral support for the topography owing to the weaker crust (Figure 4a).
There are also differences between the models after crustal root detachment.In model B2 after 80 Ma, there is crustal thinning of up to 6 km induced by the upwelling of the sublithospheric mantle because of the locally weak lower crustal rheology.The crustal thinning results in a higher Bouguer gravity anomaly, resulting in a final value of −25 mGal at 100 Ma (about 25 mGal higher than in model B1; Figure 4b).
In model B3, crustal root removal and localized ascent of sublithospheric mantle in the perturbed region triggers lithosphere delamination in a similar fashion to the models of Krystopowicz and Currie (2013) (Figure S3 in Supporting Information S1).This occurs because the weak mafic granulite rheology allows for decoupling along the interface between the lower crust and mantle lithosphere.As the mantle lithosphere detaches laterally and sinks, a wide area of thinned lithosphere is produced, where the crust is thinned by ∼6 km and there is sublithospheric mantle upwelling.Consequently, Model B3 exhibits an uplift pulse of ∼0.3 km after lithosphere removal in the perturbed region and fluctuations in the gravity anomaly of up to 50 mGal.

Mantle Lithosphere Rheology (Set C)
Model set C tests the effect of mantle lithosphere rheology on drip dynamics.Model A and model set B use a dry olivine rheology for the mantle lithosphere.In set C, the mantle lithosphere in the perturbed region has the rheology of wet olivine.This is a simple approach to test a weaker mantle lithosphere (e.g., due to arc-related processes such as hydration and presence of melt).Models C1, C2, and C3 include lower crustal eclogitization and use the same modifications in the lower crust as models B1, B2, and B3, respectively (Table 2).
The results of model set C are shown in Figures 4 and 6.Model C1 is shown in Animation 3. As in model set B, lithosphere removal is induced by both the mantle lithosphere perturbation and formation of dense eclogite.In model set C, the presence of the weaker mantle lithosphere speeds up the removal process.Root detachment occurs at 25, 26, and 28 Ma for C1, C2, and C3 respectively, 35-42 Ma earlier than in the comparable set B models.In set B, the mantle lithosphere has a dry and strong rheology, which results in a slower drip, and thus the 30 km of crustal root can fully turn into eclogite before its detachment because the mantle lithosphere inhibits its descent.Conversely, in set C the mantle lithosphere drips at ∼8 Ma and eclogite starts dripping at ∼17 Ma because at that time there is no mantle lithosphere below.In set C, detachment of the eclogitized crustal root takes place in two pulses, as in set B, but the pulses occur within 5 Ma of each other, whereas they were ∼12 Ma apart in set B.
The evolution of topography and Bouguer gravity anomaly (Figure 4) in sets B and C exhibit similar magnitudes and trends, but the timing is significantly different.For set C, the subsidence associated with the removal process and the increase in the Bouguer gravity anomaly produced by eclogitization and root removal occur 30-40 Ma earlier than in set B. Because crustal root detachment occurs faster, it produces only one main peak in the gravity anomaly at ∼25 Ma, and a positive gravity anomaly occurs for only ∼5 Ma (between 23 and 28 Ma).Note that the results for set C are only shown to 50 Ma, as the crustal root is removed by 35 Ma in these models.
Models with a weak lower crust (C2 and C3) exhibit ∼6 km crustal thinning following the removal event.As in model B2, in model C2 thinning is localized and is induced by the dripping of the crustal root and the sublithospheric mantle upwelling.As in model B3, in model C3 root removal induces regional delamination of the mantle lithosphere and thinning extends over a larger area.

Comparison of Surface Observables
The effect of lithosphere rheology on surface observables for model sets B and C is summarized in Figure 7.This figure compares the surface elevation (Figure 7a) and Bouguer gravity anomaly (Figure 7b) in the central part of the model at the time of the maximum gravity anomaly (t = t(g B max )), and at 10 Ma after this (t = t(g B max ) + 10 Ma).The time of the maximum gravity anomaly corresponds to the time at which most of the crustal root has been removed through dripping; for most models this is 1-3 Ma before the time of full detachment.This allows for a comparison of the models at similar stages, and therefore an assessment of how removal of the deep lithosphere affects the surface observations.The maximum gravity anomaly for all models is +20 to +30 mGal at t(g B max ), suggesting that there is not a significant difference between the models.However, a weaker mantle lithosphere rheology results in elevations that are 1-2 km higher than those for a stronger mantle lithosphere.This is because if the removal of the crustal root is faster, there is less subsidence time prior to the gravity anomaly peak.Additionally, a weaker lower crustal rheology allows a faster subsidence rate resulting in lower elevations (blue areas in Figure 7a) compared to strong lower-crustal rheologies (red area in Figure 7a) because there is less lateral support.At t(g B max ) + 10 Ma the trend in the surface observables is the same as at t(g B max ), but for all models the magnitudes of gravity anomaly and topography are ∼10 mGal and ∼1 km lower, respectively, due to continued lower crustal heating that decreases its density, the presence of less eclogite (which has dripped), and isostatic equilibration.
Our models show that eclogitization of the crustal root can create a local area of thinned lithosphere (∼40 km).This produces an increase in the gravity anomaly while mantle upwelling and lateral strength provide support for the load of the initial topography.Models in sets B and C show that eclogitization and crustal root removal at t(g B max ) produce a positive Bouguer gravity anomaly (up to +30 mGal) together with an elevated topography (up to 4 km).In all the models the magnitude of the Bouguer gravity anomaly at t(g B max ) is significantly lower than the observed anomaly in the SNSM (>+130 mGal) suggesting that additional factors are required to explain the excess of mass.
We note that in the models presented above, the initial crustal thickness in the central mountain region prior to lithosphere removal is fairly arbitrary as there are only limited observational constraints.Additional tests using different initial crustal thicknesses in the mountain region show that greater thicknesses produce faster lithosphere removal, higher subsidence rates, and more negative final gravity anomalies (see Supporting Information S1).

Density of the SNSM Massif (Set D)
Gravity studies of northern Colombia argue that the crystalline basement rocks that underlie most of the SNSM massif have a higher density than the surrounding rocks, owing to a more mafic composition (Ceron-Abril, 2008;  Sanchez-Rojas & Palma, 2014;Sanchez & Mann, 2015).Sanchez and Mann (2015) suggest that the SNSM density is comparable to lower crust densities.Consequently, we study the effects of a locally higher density in model set D, where the upper crust in the 125 km wide perturbed region has the same reference density of the lower crust (100 kg•m −3 higher than the upper crust).Models D1, D2, and D3 use the same parameter configurations as models C1, C2, and C3, respectively, differing only in the upper crust density in the central region (Table 2).
The results of model set D are shown in Figures 8 and 9, and Model D1 is shown in Animation 4. The temporal variations in surface observables are shown in Figures 9a and 9b, together with the results of model C1 for comparison.The overall behavior of the models in set D is similar to the equivalent models in set C. The presence of the dense upper crust has little effect on the timescales of root removal.However, in set D the elevation is consistently lower (∼0.8 km lower), owing to the additional load of the dense rocks.The trend of the gravity anomaly is similar between sets C and D, but the presence of the high-density rocks results in gravity anomalies that are ∼+50 mGal higher in set D.

Santa Marta-Bucaramanga and Oca Faults (Model E)
The SNSM is bounded by two major strike-slip faults, the Oca fault at the north and the Santa Marta-Bucaramanga fault at the southwest (Figure 1a).We include model E to test the role of these structures.Model E uses the same setup of model D1 but the adjacent faults are included using a constant internal friction angle of 2° to produce rheological weakening in a 10 km wide area that extends to the Moho at 40 km depth (Figure 10a).
The main effect of the vertical faults surrounding the thickened region is to facilitate vertical motions after the lithosphere removal episode.There are no significant differences in the dynamics or surface observables before the time of root removal (23 Ma) between models D1 and E. After removal, there is vertical motion along the faults in model E: the thickened block in the center sinks and the material at the sides rises.Figure 10b shows that the topographic profile at the time of root removal (23 Ma) is identical to model D1 but 10 Ma after removal (33 Ma) there has been 0.38 km of more subsidence in the thickened block and 0.38 km of more uplift at the sides.As the thickened block sinks along the faults, the Bouguer gravity anomaly decreases by less than 18 mGal compared to model D1, between 23 and 25 Ma (Figure 10d).Figures 10c and 10d show that the general behavior of elevation and Bouguer gravity anomaly in model E is the same as in model D1.Therefore, the effect of the faults in the dynamics is negligible and the effect in the evolution of surface expressions is minor and would only be observed ∼10 Ma after removal.

Comparison With Observations
We select model D1 as our preferred model because the magnitudes of topography and Bouguer gravity anomaly are the closest to the SNSM observations.This is a consequence of: (a) a weak local mantle lithosphere rheology (wet olivine) that produces rapid removal (within 25 Ma) and allows less time for subsidence, (b) a strong lower-crustal rheology (dry Maryland diabase) that provides more support for the load of the initial topography, and (c) the high-density upper crustal rocks in the perturbed region.In the models, a positive Bouguer gravity anomaly requires both the shallow high-density rocks of the exhumed basement and absence of a crustal root.Figure 9b shows that the initial gravity anomaly of model D1 is negative (∼−100 mGal) even with the presence of denser upper-crustal rocks, indicating that the effect of the low-density crustal root dominates.After the lithosphere removal event, the Bouguer gravity becomes positive (∼+100 mGal) and remains positive for almost 30 Ma.  between the absence of cool and dense mantle lithosphere underneath (mass deficit), the high-density shallow rocks in the SNSM (100 km•m −3 denser), and a locally thin crust (Moho depth of ∼38 km as shown in Figure 9e) with residual eclogitized material that is directly underlain by a shallow sublithospheric mantle (excess of mass near the Moho as shown in Figure 8a).Although there are uncertainties in the crustal thickness in the SNSM, the modeled value is similar to the estimates between 35 and 40 km from Sanchez-Rojas and Palma (2014), Poveda et al. (2018), andMojica-Boada et al. (2022).The observed Bouguer gravity anomaly and topography along A -A′ and B -B′ are from the EIGEN-GL04C Global gravity field model (Förste et al., 2008) and the ASTER Global Digital Elevation model V003 (NASA/ METI/AIST/Japan Spacesystems & U.S./Japan ASTER Science Team, 2019), respectively.The maximum modeled values for the Bouguer gravity anomaly and surface elevation are +103 mGal and 3.3 km, respectively, compared to the maximum observed values of +126 mGal and 5.1 km.The observed gravity shows that the positive peak is surrounded by two negative peaks that are associated with the low-density sediments of the lower Magdalena basin southwest of the SNSM and the Baja Guajira basin to the northeast.These gravity minima do not appear in the modeled gravity because the surrounding basins are not included in the numerical model.

3D Bouguer Gravity Anomaly
The 2D modeling setup and gravity anomaly calculations assume that the density structure extends infinitely along strike.However, the observed topography and gravity anomaly have a 3D nature and therefore the 2D calculations overestimate the gravity anomaly.To asses this, we use the 2D density structure of model D1 at 23 Ma to generate a 3D density structure assuming axial symmetry (Figure S6b in Supporting Information S1), considering that the SNSM is approximately axisymmetric.This is used to calculate the 3D gravity anomaly (Figure S6a in Supporting Information S1), assuming that the contribution of each element in the domain is the anomaly produced by a sphere of the same volume, following Telford et al. (1990).The resultant gravity anomaly is localized over the SNSM, with a similar width to the 2D anomaly (Figures 11b and 11d).However, the maximum value is +58 mGal.This shows that even though the gravity anomaly is still positive, the 2D calculations overestimate the magnitude of the anomaly by 45 mGal.Therefore, the 2D gravity results should be considered as first order approximations that show the behavior of the gravity anomaly during a lithosphere removal episode.

Isostatic Versus Non-Isostatic Topography
The models exhibit an overall subsidence in the central region as the result of the combined effects of the evolving density structure and mantle dynamics.To separate these effects, we calculate the isostatic topography, which is the topography that would be observed if there was isostatic equilibrium based on the density structure at each time (Figure 9c). Figure 9d shows the non-isostatic topography which is the difference between the total topography and isostatic topography.Non-isostatic effects involve the contribution to topography produced by both lateral lithosphere strength and vertical stresses from the dynamics of the deep lithosphere and underlying mantle.For models D1 and D2, the isostatic topography is significantly lower than the total topography throughout the removal episode, whereas model D3 is closer to isostatic balance throughout its evolution.This shows that the non-isostatic forces are more significant in models D1 and D2, compared to model D3.This is possibly due to the support provided by a rheologically stronger lower crust in model D1 (e.g., dry diabase), compared to the weaker crust in model D3 (e.g., mafic granulite).In model D2, the perturbed region has a mafic granulite lower crust, with dry diabase at the sides, and the role of non-isostatic forces is intermediate between models D1 and D3.
In model D1, the non-isostatic topography reaches a maximum of 3.5 km at 23 Ma (∼2 Ma before full detachment).Figures 11a and 11c shows that at 23 Ma the topography is predominantly non-isostatic.This suggests that the coexistence of a positive Bouguer gravity anomaly and high topography at this time is possible due to the delay in subsidence produced by non-isostatic effects.

Melting
Our study primarily considers the effect of lithosphere removal on surface topography and gravity.The models also make other predictions.As shown above, a lithosphere drip results in lithosphere thinning and causes upwelling of the sublithospheric mantle to fill the gap created by the drip.Decompression of upwelling mantle and simultaneous lower crustal heating may result in melting of both the mantle and crust, and therefore magmatism may accompany lithosphere removal (e.g., Kay & Mahlburg Kay, 1993;Mahlburg Kay et al., 1994).To assess this, we use the pressure-temperature conditions in model D1.Our calculations use the dry granite solidus from Elkins-Tanton (2005), and the solidus of partially hydrated peridotite with 0.15 wt% H 2 O (Katz et al., 2003), for the lower crust and mantle.The models do not include the transport of melt to the surface, nor the effects of melt on rheology and density.
Figures 12a and 12b shows model D1 at 10 Ma after root detachment (t = 33 Ma) to allow enough time for significant melt production.In the central region, melting of the sublithospheric mantle and deep crust are predicted.There is also a zone between 80 and 100 km depth outside of the main drip region where lower mantle lithosphere has been partially removed and conditions favor mantle melting.
Figure 12c shows the temporal evolution of mantle melt volume (per unit along strike) for both the full model domain and for a local 50 km wide region centered in the elevated zone.This is calculated by monitoring the regions of the model that are above the solidus, assuming a melt fraction per degree above the solidus of 0.3% for the mantle (Katz et al., 2003) and 0.75% for the crust (Annen & Sparks, 2002).The onset of melting in the central region occurs at about 23 Ma, corresponding to mantle upwelling after root detachment and crustal heating.These calculations suggest that mantle and crustal melt volumes in the central region could be significant, with approximately equal amounts of each.

Seismic Velocity Structure
Our models can also be used to assess how lithosphere removal may affect the mantle seismic structure.For this, we use the temperature, pressure, and composition of model D1 to calculate the expected seismic P-and S-wave velocity structure.This uses Perple_X (Connolly, 2005) to obtain the bulk and shear moduli, assuming a pyrolite composition, with a correction for attenuation (Quiroga, 2022;van Wijk et al., 2008).The crust is excluded from this analysis, owing to uncertainties in crustal composition.
Figure 13 shows the calculated P-and S-wave velocities and their corresponding anomalies (based on the average velocity at each depth) at a time of 23 Ma.After removal of the dense crustal root, the gap created in the mantle lithosphere has a clear low velocity anomaly of up to −9.7% in P-waves and −19.1% in S-waves, with absolute values of 6.9 and 3.5 km•s −1 , respectively.Note that the calculations do not include the effects of melt. Figure 13a shows where melt is predicted and therefore the seismic wave velocities could be even lower.These calculations indicate that lithosphere removal may be associated with a significant low-velocity anomaly in the shallow mantle.This could persist for over 10 Ma following the removal event, as the lithosphere gap remains hot (Figure 13b).

Surface Heat Flow
We calculate the surface heat flow for model D1 by multiplying the vertical temperature gradient at the surface by the thermal conductivity (Figure 14).Between 0 and ∼23 Ma, the surface heat flow in the perturbed region has a decrease of <5 mW•m −2 due to the topographic changes and shallow deformation produced by crustal root removal (Figure 14b).From 23 to 50 Ma, there is a progressive increase to 81 mW•m −2 caused by upwelling of the hot sublithospheric mantle and heating of the overlying crust following root removal at 23 Ma.The largest increase in heat flow occurs after ∼30 Ma.For comparison, model A shows the surface heat flow in the absence of crustal root removal; there is an increase from 50 to 54 mW•m −2 between 0 and 50 Ma due to internal radiogenic heating and conductive heating of the crust following lithosphere thinning (Figure 14b).(Katz et al., 2003).(c) Calculations of melt volume per unit length along strike for the lower crust.Local mantle melt volume is the melt volume in a 50 km width centered in the elevated region.

Origin of High Topography in the SNSM
We have assumed that the current SNSM topography (>5 km) was produced through past tectonic processes: volcanic arc tectonics in the Jurassic (Ramírez et al., 2020) and the collision of the Caribbean plateau in the Late Cretaceous-Early Paleocene (Cardona et al., 2010).Therefore, our models start with a thick crust with initially high topography in the perturbed region.The models demonstrate that with these initial conditions a lithosphere removal episode produces subsidence.However, results show that upward forces from the upwelling sublithospheric mantle and lateral strength provide support for the pre-existing topography and delay subsidence.In the models this is only possible because the non-isostatic forces overcome the load of the initial topography (Figure 9d).

Origin of the High Gravity Anomaly in the SNSM
The models demonstrate that gravitational removal of an eclogitized crustal root can produce a local positive Bouguer gravity anomaly in an area of elevated topography.In model D1, this results in a 2D gravity anomaly of +103 mGal and topography of 3.3 km, suggesting that a lithosphere drip is a viable mechanism to explain the SNSM observations.Based on the model results and geological observations, we suggest that a drip could have occurred at two possible times: in the Eocene (56-40 Ma) or within the last 2 Ma.These possibilities are shown in Figure 15, and discussed below.However, it is also possible that there was lithosphere removal at both times.

Lithosphere Removal in the Eocene (56-40 Ma)
The first possibility is that lithosphere removal occurred in the Eocene (Figure 15a).At this time, the SNSM was located in the volcanic arc region with subduction of the Caribbean plate below South America (Cardona et al., 2010).It was at an intermediate location between the Colombian Central Cordillera and its current position, and remained there until the late Miocene (Ujueta, 2003).A removal event at 56-40 Ma would predate the slab flattening that has been proposed to explain the cessation of Paleocene-Eocene magmatic activity (e.g., Taboada et al., 2000), and therefore a low-angle Caribbean slab would not block lithospheric foundering.
We suggest that removal at this time could have been driven by metamorphic eclogitization of thickened SNSM crust.Alternatively, root densification could have been produced by magmatic processes where differentiation creates a lower layer of high-density "arclogite" rocks (Ducea, Chapman, Bowman, & Balica, 2021).This process has been shown to promote gravitational removal even in garnet-free arcs where eclogite is not produced (e.g., Jull & Kelemen, 2001), which is possibly the case in the SNSM (e.g., Duque-Trujillo et al., 2019).
There were two main Paleocene-Eocene magmatic events that were separated by a magmatic gap.The first event produced the Paleocene leucogranites (age ∼65 Ma).Duque-Trujillo et al. (2019) argue that the protolith of these rocks had a mafic composition and was within the stability field of amphibole, garnet, and rutile, which are all common minerals in arclogite assemblages (Ducea, Chapman, Bowman, & Triantafyllou, 2021).The second event created the Eocene Santa Marta Batholith (SMB) at ∼56 Ma.Duque-Trujillo et al. (2019) report that the SMB has substantial ultramafic garnet-free cumulates.Some of these have been classified as pegmatitic pyroxenites, perhaps showing evidence of a densified lower crust outside the garnet stability field that could be gravitationally unstable.Duque-Trujillo et al. (2019) suggest that the SMB magmatism does not correspond to typical Andean-type subduction magmatism because it was localized, low-volume and short term (∼10 Ma).Instead, they propose that the melting of thickened lower crust could produce the required parental magma.In our models, crustal root removal takes place on timescales of ∼10 Ma following the onset of root densification.Also, decompression and heating after removal predict lower-crustal and shallow sublithospheric mantle melting (Figure 12).Therefore, a progressive densification of the lower crust starting at ∼65 Ma (or earlier), followed by a removal event at ∼56 Ma may explain the trends in the anomalous SMB magmatism.A detailed study of the mafic intrusions is required to evaluate the possibility of a previous episode of eclogitization and root removal in the SNSM.
One problem with an Eocene removal event is that the surface expressions must persist for >40 Ma in order to explain the present-day high elevation and positive Bouguer gravity anomaly.The models in this study show only a short-lived period (∼5 Ma) of simultaneous high elevation and positive gravity anomaly, and therefore additional factors may be required, such as elastic stresses, mantle upwelling, and support from the present-day Caribbean slab.

Recent Lithosphere Removal (<2 Ma)
The other possibility is that lithosphere removal occurred recently (Figure 15b).The configuration of the Caribbean slab under the SNSM is poorly constrained.A recent seismic receiver function study by Mojica-Boada et al. (2022) suggests that the Caribbean slab is 60-110 km below the Moho.This would provide enough space for lithospheric foundering between the lower crust and the Caribbean slab.
According to our models, the peak in the Bouguer gravity anomaly after lithosphere removal is short-lived (<5 Ma), and therefore the current gravity observations are compatible with recent removal.Additional evidence that favors a recent lithospheric drip comes from seismic observations.Regional seismic tomography studies show high velocity anomalies under the SNSM corresponding to the subducting Caribbean slab, but do not have suitable resolution to image the shallow (∼40 km) and local structure (e.g., Cornthwaite et al., 2021).Some local studies suggest low shear wave velocities (<3.8 km•s −1 ) in the deep crust and shallow mantle (32-35 km depth) that are partially correlated to the high topography and gravity anomaly (Poveda et al., 2018;Syracuse et al., 2016).These values are consistent with the shear wave velocities predicted for model D1 near the Moho (Figure 13a).However, other studies suggest a higher seismic velocity in the upper mantle (e.g., Lü et al., 2022).Higher resolution seismic studies are required to verify whether the current velocity structure agrees with the model predictions.
Our calculations predict that lithosphere removal would induce lower crustal and upper mantle melting, whereas the last pulse of magmatism in the SNSM was at about 54-49 Ma (Duque-Trujillo et al., 2019).Figure 12c shows that the melt volume produced by lithosphere removal in the perturbed region is only significant 10 Ma after the removal episode.This means that if there was recent lithosphere removal there is possibly melt in the deep lower crust and mantle that has not yet reached the surface.

Alternative Models
Previous studies have discussed the simultaneous occurrence of a localized high elevation and positive Bouguer gravity anomaly in the SNSM region.The challenge has been to explain how the topography is supported and the cause of the local excess of mass.
One hypothesis is that the underthrusting of the buoyant Caribbean plateau and oceanic seamounts could provide support for the SNSM topography (e.g., Case & Macdonald, 1973;Montes et al., 2005).In addition, Montes et al. (2005) propose that the Caribbean slab could have caused an upward bulge of the continental Moho, creating a thin SNSM crust and positive gravity anomaly.However, subducting/underthrusting seamounts at very shallow depths (10-30 km) generally do not generate a significant positive Bouguer gravity anomaly signature.Some examples include the subducting seamount under the Cascadia accretionary complex in North America (∼−15 mGal) (Tréhu et al., 2012) and the Papudo and San Antonio seamounts subducting under Chile, South America (<−50 mGal) (Manríquez et al., 2014).Moreover, an underthrust seamount would only explain the localized nature of the gravity anomaly and its correlation with topography if it is located right under the topographic relief and happens to have the same width.
Another possibility is that the SNSM topography is laterally supported by a rigid and strong continental lithosphere.In northwestern Colombia, the elastic thickness is 26-40 km (Arnaiz-Rodríguez & Audemard, 2014;Tassara et al., 2007).However, the results from model E show that the presence of the Oca and Santa Marta -Bucaramanga faults may weaken the lithosphere.This may limit the elastic lithosphere strength/support, especially 10 Ma after the removal episode.While the presence of a rigid lithosphere could provide mechanical support for the high topography, it does not explain the high-magnitude positive Bouguer gravity anomaly, meaning that there must be an additional factor that creates an excess of mass.

Additional Considerations
The models presented here investigate the dynamics of lithosphere removal and the associated surface expressions.This mechanism can explain the localized nature of the Bouguer gravity anomaly and its simultaneous occurrence with a high topography.However, the models do not explain their full magnitudes.There are several limitations of our models.
Our models are 2D and do not include elasticity.Previous studies show that 3D lithospheric drips occur faster, may produce greater topographic deflections (e.g., Houseman & Gemmer, 2007;Pysklywec & Cruden, 2004), and higher Bouguer gravity anomalies (e.g., Quiroga, 2022).In addition, Arnaiz-Rodríguez and Audemard (2014) suggest that elastic strength in the crust is important for supporting the SNSM massif.Kaus and Becker (2007) show that the inclusion of crustal elasticity in lithosphere drip models result in larger surface deflections and larger downwelling rates.Future 3D models including elasticity are needed to assess the implications for the general dynamics and surface expressions of lithosphere removal.Furthermore, the models address only the dynamics of the continental lithosphere and do not include the subducting Caribbean slab or regional tectonics.Seismic studies show that the SNSM region is underlain by the Caribbean slab, but there are uncertainties in the slab location and depth.It is unclear how the presence of the Caribbean plate may affect lithosphere removal, especially if the removal event occurred in the last 2 Ma.On one hand, the Caribbean plate may hinder the foundering of a crustal root; on the other hand, subduction may induce a lateral drag that promotes removal (e.g., Currie et al., 2015).Moreover, the presence of subducting seamounts could contribute to the positive Bouguer gravity anomaly if they are shallow enough to produce an excess of mass near the Moho, due to their mafic composition.Future models should include Caribbean plate subduction to assess its effect on the dynamics and structure of the overlying continent, including the SNSM.

Conclusions
This work uses numerical models of localized lithosphere removal to study the unusual surface topography and Bouguer gravity anomaly in the SNSM region of northwestern Colombia.We find that foundering of a mantle lithosphere perturbation can induce eclogitization in the crustal root, which is then gravitationally removed, leaving a locally thinner crust (∼38 km) underlain by the hot sublithospheric mantle.This results in a local excess of mass and a positive Bouguer gravity anomaly.We propose that the high gravity anomaly in the SNSM region was produced by this mechanism.
The overall topography depends on non-isostatic forces, such as lateral lithosphere strength and upwelling of the sublithospheric mantle, while the removal timescales depend on lithosphere rheology: a strong lower crustal rheology (e.g., dry diabase) provides more lateral support and sustains the initial topography for longer compared to a weak rheology (e.g., mafic granulite); a local strong mantle lithosphere rheology (e.g., dry olivine) results in removal within 70 Ma, while a weak local rheology (e.g., wet olivine) results in removal within 25 Ma allowing less subsidence time.Our models show that the SNSM observations are compatible with a strong lower crust and a locally weak mantle lithosphere (model D1), where removal produces a positive 2D Bouguer gravity anomaly of +103 mGal (∼+58 mGal in a 3D geometry), supporting a high topography (3.3 km).
The history of the SNSM is tectonically complex and previous hypotheses attribute its formation and current structure to Caribbean plate subduction and continental strength.However, these hypotheses fall short in explaining the concurrent high topography and localized positive Bouguer gravity anomaly.We propose that the observations can be explained by gravitational removal of the deep lithosphere below an area of pre-exiting high topography.The observed gravity anomaly in the SNSM (>+100 mGal) requires gravitational crustal thinning, high-density shallow crustal rocks (e.g., Sanchez & Mann, 2015), and the presence of residual eclogite in the lower crust.We suggest that lithosphere removal could have occurred either very recently (in the last 2 Ma), creating an area of low seismic velocity in the shallow mantle, or during the Eocene (56-40 Ma), also affecting regional magmatism, provided that elastic stresses and/or mantle flow provide long term topographic support.Taken together with the earlier crustal shortening and block rotation, our model provides a plausible mechanism for the enigmatic SNSM.

Figure 1 .
Figure 1.(a) Topographic map of the study region in northern Colombia (red box in inset map).BGB: Baja Guajira Basin, OF: Oca Fault, SNSM: Sierra Nevada de Santa Marta, LMB: Lower Magdalena Basin, CRB: Cesar-Rancheria Basin, PR: Perijá Range, SMBF: Santa Marta-Bucaramanga Fault.The red barbed line is the plate margin where the Caribbean plate (CAR) converges with South America.The cyan, magenta and green contours delimit respectively the three geotectonic provinces: (1) Sierra Nevada, (2) Sevilla and (3) Santa Marta.(b) Bouguer gravity anomaly map from the EIGEN-GL04CGlobal gravity field model(Förste et al., 2008).Profiles of topography and gravity anomaly (A-A′ and B-B′) are used to compare the models with observations.

Figure 2 .
Figure 2. (a) Initial model setup.The inset box indicates the area shown in the following figures.(b) Unperturbed density structure.(c) Initial temperature conditions.The continental geotherm is stretched in the thickened region keeping the Moho temperature constant.(d) Reference viscosity structure at a constant strain rate of 10 −15 s −1 (black).Black dotted lines show the reference viscous rheologies above plastic yielding including: Wet Quartzite (WQ), Dry Maryland Diabase (DMD), and Dry Olivine (DO).Colored lines show other viscous rheologies tested including: Eclogite (Ec) with f = 0.1, Mafic Granulite (MG) and Wet Olivine (WO).See Table1for references.UC = upper crust, LC = lower crust, ML = mantle lithosphere, SLM = sublithospheric mantle, and LAB = lithosphere-asthenosphere boundary.

Figure 3 .
Figure 3. Evolution of the reference model (model A) at different times, showing the effect of the mantle lithosphere drip.(a) Surface topography profiles (H).(b) Bouguer gravity anomaly profiles (g B ). (c) Snapshots of the density and thermal structure.Zero topography corresponds to the surface of the modeling domain.

Figure 4 .
Figure 4. Time variation of surface observables averaged over a 50 km width centered over the initial perturbation, for the reference model (black) and model sets B (blue) and C (red).(a) Surface topography (H).(b) Bouguer gravity anomaly (g B ).Green arrows indicate the gravity response to the mantle lithosphere drip and black brackets indicate the response produced by crustal root removal.

Figure 5 .
Figure 5. Evolution of (a) model B1, (b) model B2 and (c) model B3, including profiles of surface topography (H) and Bouguer gravity anomaly (g B ), and snapshots of the density and thermal structure, at times before, during and after the main lithosphere removal episode.Model set B includes eclogitization and explores the effect of different lower crustal rheologies.

Figure 6 .
Figure 6.Evolution of (a) model C1, (b) model C2, and (c) model C3, including profiles of surface topography (H) and Bouguer gravity anomaly (g B ), and snapshots of the density and thermal structure and at times before, during and after the main lithosphere removal episode.Model set C has a weaker mantle lithosphere rheology below the perturbed region.Models include eclogitization and explore the effect of different lower crustal rheologies.

Figure 7 .
Figure 7. Magnitude of (a) surface topography (H) and (b) Bouguer gravity anomaly (g B ), averaged over a 50 km width centered on the perturbed region, at the time of the maximum gravity anomaly after crustal root removal (top) (t = t(g B max )), and 10 Ma after the peak gravity anomaly (bottom) (t = t(g B max ) + 10 Ma).The plots show the results of model sets B and C as a function of mantle lithosphere rheology (y-axis) and lower crustal rheology (x-axis): DMD = Dry Maryland diabase; LMG = Local mafic granulite; MG = Mafic granulite; WO=Wet olivine; DO = Dry olivine.

Figure 8 .
Figure 8. Evolution of (a) model D1, (b) model D2 and (c) model D3, including profiles of surface topography (H) and Bouguer gravity anomaly (g B ), and snapshots of the thermal and density structure at times before, during and after the main lithosphere removal episode.Model set D is the same as set C, but it includes high-density upper crustal rocks in the perturbed region.

Figure 11
Figure 11 compares the results of model D1 to the SNSM observations along profile lines A -A′ and B -B′ in Figure 1 at a model time of 23 Ma, corresponding to the maximum peak in the gravity anomaly (∼2 Ma before full detachment of the crustal root).At this time the local Bouguer gravity anomaly is the result of the tradeoff

Figure 9 .
Figure 9.Time variation of surface observables and Moho depth, averaged over a 50 km width centered in the perturbed region for model set D, and model C1 as a reference.(a) Surface elevation (H) (total topography).(b) Bouguer gravity anomaly (g B ). (c) Isostatic topography (green) with total topography for comparison (purple).(d) Non-isostatic topography.Blue arrow indicates time of maximum uplift rate (0.4 mm yr −1 ) in model D1.(e) Moho depth.

Figure 11 .
Figure 11.Comparison between topography and Bouguer gravity anomaly profiles of the preferred model (model D1) during root removal at 23 Ma (∼2 Ma before full detachment) with observed data along profiles A -A′ (a and b) and B -B′ (c and d); profile locations in Figure 1b.Initial topography (H) and gravity anomaly (g B ) profiles (t = 0 Ma) are shown for comparison (blue).(a) and (c) Modeled topography (red).Data from the ASTER global digital elevation model (black).The calculated non-isostatic topography (orange) shows that non-isostatic effects are localized in the elevated region.(b) and (d) Modeled Bouguer gravity anomaly (red).Data from the EIGEN-GL04C global gravity field model (black).Modeled 3D Bouguer gravity anomaly assuming a 3D axisymmetric density structure (green).SNSM: Sierra Nevada de Santa Marta, LMB: Lower Magdalena Basin, BGB: Baja Guajira Basin, CRB: Cesar-Rancheria Basin, PR: Perija Range.

Figure 10 .
Figure 10.(a) Setup of model E, showing the central region of the model (dashed area in Figure 2a).Lime dashed lines enclose areas where there is a constant friction angle (ϕ) of 2° producing rheological weakening to model the Santa Marta-Bucaramanga and Oca faults.(b) Topographic profiles of model E at the time of root removal (23 Ma), and 10 Ma after (33 Ma).Black arrows show the location of the weak zones.Time variation of (c) surface elevation (H) and (d) Bouguer gravity anomaly (g B ), averaged over a 50 km width.For comparison, the results of model D1 are shown in black.

Figure 12 .
Figure 12.(a) Density structure of the preferred model (model D1) at t = 33 Ma (8 Ma after removal).Contours enclose regions where the P-T conditions are above the wet peridotite and dry granite solidi, for the mantle (purple infill) and lower crust (white infill), respectively.(b) Temperature profile (blue line) at x = 230 km and t = 33 Ma, with the dry granite (Elkins-Tanton, 2005) and the wet peridotite solidus (Katz et al., 2003).(c) Calculations of melt volume per unit length along strike for the lower crust.Local mantle melt volume is the melt volume in a 50 km width centered in the elevated region.

Figure 13 .
Figure 13.(a) Calculated seismic structure for model D1 at t = 23 Ma, showing P-wave velocities (left), S-wave velocities (center), and vertical profiles at the side (x = 34 km) and central region (x = 230 km).Black contours enclose regions where the pressure and temperature conditions are above the peridotite solidus, which could have lower velocities.(b) P and S wave velocity anomalies at 23 Ma and profiles at the center at 23 and 33 Ma.

Figure 14 .
Figure 14.(a) Profiles of surface heatflow for model D1 at 23 and 33 Ma.(b) Time variation of surface heat flow, averaged over a 50 km width centered in the perturbed region for model D1 and model A as a reference (black).The green arrow indicates the time of root removal.

Figure 15 .
Figure 15.Proposed evolution of the lithosphere structure along profile B -B′.(a) Hypothesis of a lithosphere removal in the Eocene at 20 Ma time intervals.(b) Hypothesis of a recent removal (<2 Ma) at 15 Ma intervals.All plots show the present-day topography (H) as a reference.The present-day plots show the current Bouguer gravity anomaly (g B ), and the gray arrow indicates the location of the Caribbean-South America plate boundary.In the present-day plots, the position of the Caribbean slab and continental Moho from receiver functions (Mojica-Boada et al., 2022) are shown with orange and yellow lines, respectively; other geometry is inferred.Diagonal lines in the Caribbean slab at other times indicate that its position and geometry is not constrained.The black arrows show mantle flow induced by lithosphere removal.
Note.The reference model does not include eclogitization.The parameters for eclogite are used in all other models. a