Holocene relative sea-level changes along the Caribbean and Pacific coasts of northwestern South America

Abstract

Predicting coastal change depends upon our knowledge of postglacial relative sea-level variability, partly controlled by glacio-isostatic responses to ice-sheet melting. Here, we reconstruct the postglacial relative sea-level changes along the Caribbean and Pacific coasts of northwestern South America by numerically solving the sea-level equation with two scenarios of mantle viscosity: global standard average and high viscosity. Our results with the standard model (applicable to the Pacific coast) agree with earlier studies by indicating a mid-Northgrippian high stand of ~2 m. The high-viscosity simulation (relevant to the Caribbean coast) shows that the transition from far- to intermediate-field influence of the Laurentide Ice Sheet occurs between Manzanillo del Mar and the Gulf of Morrosquillo. South of this location, the Colombian Caribbean coast has exhibited a still stand with a nearly constant Holocene relative sea level. By analyzing our simulations considering sea-level indicators, we argue that tectonics is more prominent than previously assumed, especially along the Caribbean coast. This influence prevents a simplified view of regional relative sea-level changes on the northwestern South American coast.

INTRODUCTION

Predicting coastal change depends upon deciphering postglacial relative sea-level (RSL) variability and its influence on modern morphology (FitzGerald et al., 2008; Passeri et al., 2015). Postglacial RSL changes control modern coastal environments, as they dictate nearshore sediment accommodation and supply (Everts, 1987; Törnqvist et al., 2020). Thus, refining our knowledge of postglacial RSL within a morphodynamic framework allows us to quantify coastal response in an increasingly human-influenced world (Nicholls and Cazenave, 2010; Kopp et al., 2015; Nerem et al., 2018; Nienhuis and van de Wal, 2021).

Here, we compare postglacial RSL change simulations with available sea-level (SL) indicators for the Caribbean and Pacific coasts of northwestern (NW) South America to decipher the spatial–temporal variability of Holocene RSL. Understanding such variability would require unraveling the glacial isostatic adjustment (GIA) response, which is critical for understanding the RSL drivers. We focus on an outstanding question: What GIA processes are responsible for RSL change—intermediate field or far field?

Considering the coastal evolution of NW South America, tectonics creates a complex deformation that prevents a simplified view of continental vertical change influencing RSL (Garrett et al., 2020). In addition, exploring postglacial RSL change informs us about modern coastal morphodynamics. Any RSL framework would imply either transgressive (intermediate-field) or regressive (far-field) conditions during the Holocene with implications for contemporary coastal change (Shadrick et al., 2022; Nienhuis et al., 2023).

Our work is organized as follows. The “Background” section includes information on postglacial relative SL changes, the geology of coastal Colombia, Holocene RSL changes along NW South America, and GIA modeling. The “Methods” section outlines details of our RSL change simulations with SELEN⁴ and the compilation of SL indicators. We then present our results in three subsections covering modeled postglacial RSL changes, modeled Holocene RSL changes, and compiled SL indicators (see the compilation in the Supplementary Material at http://dx.doi.org/10.17632/7nhpbhvfnz.2). Our “Discussion” addresses ocean versus ice loadings within GIA modeling, three-dimensional Earth structure along NW South America, the role of tectonics in RSL change, and sediment isostasy. Our conclusions are provided in the final section.

BACKGROUND

Postglacial RSL change

Coastal change in the twenty-first century will be multifaceted, depending on marine and continental processes that modify the coastal terrain (e.g., Kennedy et al., 2020; Nienhuis et al., 2020; Anfuso et al., 2021). In turn, the extent of coastal processes depends on RSL, which is determined by climatic, oceanographic, and solid Earth processes, including GIA, that control the relative position of the continental and ocean-surface levels (Pirazzoli, 1996, p. 5; Milne et al., 2009; Gregory et al., 2019).

During the Holocene (last 11.7 ka; Walker et al., 2018), GIA controlling RSL depended on the location on Earth (Lambeck et al., 2011; Kopp et al., 2015). Locally, ocean and continental surfaces varied, because ice melting redistributed ocean water; the Earth's rotation changed; and the lithosphere was deformed by water, sediment, and ice loads (Rovere et al., 2016). Therefore, ice-sheet melting produced GIA responses that depended on local conditions, that is, inner Earth structure and distance from the ice sheets (Whitehouse, 2018).

In the far field (away from the immediate influence of ice-sheet forebulge collapse), hydro-isostatic balances dominated after the ice melted and gravitational effects ceased (Clark et al., 1978). This process created Mid-Holocene high stands that implied RSL fall and marine regression during the Late Holocene. Evidence of this regression (e.g., beach ridges, strandplains, fluvial–marine terraces) is ubiquitous in the far field along tectonically inactive regions (Isla, 1989; Hein et al., 2014; Cooper et al., 2018). This regression contrasts with the transgression caused by a gradual inundation along coasts in the intermediate field, that is, controlled by the collapse of the ice-sheet forebulge that lowered the solid Earth and raised RSL. Such postglacial history implies coastline retreat with submerged evidence of Holocene transgression (Mitrovica and Peltier, 1991).

When analyzing RSL changes along the NW South American coast, as it is (geologically) divided into the Pacific and Caribbean coasts by the Panama Isthmus, we find a dichotomy in whether intermediate- or far-field GIA processes have dominated during the Holocene. Therefore, we resort to separating the Caribbean from the Pacific coasts in our analysis. Arguably, the low-latitude location makes far-field effects dominant. Also, the Caribbean coast (especially the northern region) will likely be at the limit of intermediate-field influence from the Laurentide Ice Sheet collapse (Khan et al., 2015). We, therefore, explore the dichotomy between intermediate- or far-field dominance, with the hypothesis that GIA varies latitudinally along the NW South American coast. If intermediate-field effects dominate, RSL changes will exhibit a continuous rise (e.g., Louisiana, USA; Fig. 3F in Khan et al., 2015). In contrast, if far-field effects dominate, RSL changes would exhibit a Mid-Holocene high stand, that is, an RSL higher than the present, typically at 6 ka (e.g., Rio Grande do Norte, Brazil; Khan et al., 2015; Fig. 3I). See also figure 8 in Yokoyama and Purcell (2021) for additional schematics.

Geology of NW South America and coastal Colombia

Northwestern South America exhibits a complex geological setting that frames coastal evolution and postglacial RSL changes (Fig. 1). Here, we focus on the features influencing crustal motions controlling postglacial RSL changes. For in-depth accounts of the geology and tectonics of NW South America, we suggest seminal works by Nygren (1950), Duque-Caro (1990a, 1990b), Vernette et al. (1992), Kellogg and Vega (1995), Colmenares and Zoback (2003), and IGAC and INGEOMINAS (2006), as well as more recent accounts by Cediel and Shaw (2019), Montes et al. (2019), Bustamante et al. (2016, 2024), Restrepo et al. (2021, 2023), and references therein.

Figure 1. (A) Study location geology and tectonics in South America and (B) general geology and tectonics of northern South America (black box in A), including plate boundaries, faults, fault systems (letters), and primary tectonic affinities (colors). Modified from Montes et al. (2019) and Mora-Paez et al. (2019). Plate motions relative to the north Andean block with slip rates in centimeters per year (cm/yr) after Kellogg and Dixon (1990). Block velocity relative to the stable north Andean block (1 cm/yr) from Mora-Páez et al. (2020). Coastline data (A) from the Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG) from the M_Map routine (Pawlowicz, 2020).

The tectonic framework of NW South America is a product of the junction of three tectonic plates, two oceanic (Nazca and Caribbean) and one continental (South American), which creates two subduction zones. Each subduction zone corresponds to either the Caribbean or the Pacific coast of Colombia and decisively determines the geology of each coast (Correa and Pereira, 2019). Both coasts can be considered collisional (Inman and Nordstrom, 1971). However, their coastal environments differ, ranging from spits, lagoons, mangrove swamps, and beach and dune complexes to plunging cliffs, shore platforms, and pocket beaches (Correa and Morton, 2010a, 2010b).

Geology frames six coastal fringes along the Caribbean (Correa and Morton, 2010a; Correa and Pereira, 2019; Fig. 2). The Gulf of Urabá's western margin exhibits an oceanic magmatic arc that creates plunging cliffs and pocket beaches, contrasting with intertidal Quaternary deposits (beaches–dunes, lagoons) to the southern and eastern margins (Correa et al., 2016). North of the Gulf of Urabá, the Quaternary deposits switch to fluvial–marine terraces that emerged during the last ~5 ka as a result of tectonics and (apparently) hydro-isostasy (Correa and Paniagua-Arroyave, 2016). These terraces formed out of sediments from a Cenozoic accretionary prism called the Sinú-San Jacinto Deformed Belt, spanning the coastal fringe to the Magdalena river delta, except for the Gulf of Morrosquillo, which exhibits beach–dune and lagoon deposits. North of the Magdalena River delta, there are plunging cliffs associated with the Cretaceous/Jurassic metamorphic rocks of the Santa Marta Massif. Finally, we find Quaternary deposits (beaches, dunes, and lagoons with aeolian deposits) from the Sierra Nevada of Santa Marta to the Guajira Peninsula (Correa and Pereira, 2019).

Figure 2. Coastal geomorphology, including main cities and sea-level indicators along the Caribbean coast of NW South America. Note the location of cliffs, shore platforms, and paleo-cliffs (emerged terraces) in concordance with tectonic affinities (e.g., Arboletes). Beaches, beach ridges, etc., can be found elsewhere (e.g., the Atrato River delta). Geomorphology data from figure 2 in Correa and Pereira (2019). Coastline data from the Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG) and elevation data from the ETOPO1 1 arc-minute global relief model, both from the M_Map routine (Pawlowicz, 2020).

Along the Pacific coast, the subduction of the Nazca underneath the South American plate is relatively steep, which produces a narrow continental shelf and accumulation zones (Fig. 3). The coastal zone is flanked by an oceanic plateau fragment (related to the Cali-Patía Fault system) to the east and the Garrapatas Fault to the north (González et al., 2014; Montes et al., 2019). The steep subduction has produced folding and faulting systems along regional lineaments. As a result, there are four major coastal blocks and basins: the Baudó Mountain Range, the Atrato Basin, the San Juan Basin, and the Tumaco Basin (IGAC and INGEOMINAS, 2006). The coastal morphology is characterized by Quaternary deposits related to the Patía, San Juan, Mira, and other smaller deltas that accumulate in the accommodation created by geologic structures (Correa, 1996; Correa and Morton, 2010b). The coastal environments are depositional to the south, primarily wave- and tide-dominated deltas with meso-tidal barrier islands, beaches, and mangrove swamps. To the north, an oceanic plateau fragment (the Baudó Mountain Range) is found at the ocean–continent intersection, creating plunging cliffs and related erosional environments (Correa and Pereira, 2019).

Figure 3. Coastal geomorphology, including main cities, major earthquakes, and sea-level indicators along the Pacific coast of NW South America. Note the location of cliffs, shore platforms, and paleo-cliffs (emerged terraces) in concordance with tectonic affinities (e.g., Juradó). Beaches, etc., are found elsewhere (e.g., the Mira River delta). Seven earthquakes are registered with Mw > 7.0 (orange circles with black numbers; please see our Supplementary Material at http://dx.doi.org/10.17632/7nhpbhvfnz.2 for details). Earthquakes from U.S. Geological Survey's app Latest Earthquakes v. 1.3.1 (https://earthquake.usgs.gov). Geomorphology from figure 12 in Correa and Pereira (2019). Coastline data from the GSHHG Database and elevation data from the ETOPO1 1 model, both from the M_Map routine (Pawlowicz, 2020).

In NW South America, horizontal crustal deformation has been measured with GPS technology, obtaining a NE (60°) motion of the north Andean block at 8.6 mm/yr, an eastward collision of the Panama arc at 15–18 mm/yr, and a Caribbean subduction at 13 mm/yr. Due to significant uncertainties, vertical motions have not yet been published (Mora-Páez et al., 2019). In addition to faulting-related tectonics, vertical motions contributing to RSL change also depend on mud diapirism along both coasts (Martínez and López Ramos, 2010; Carvajal, 2016) and co-seismic motions linked to major earthquakes, with at least seven events of Mw > 7.0 since 1906 along the Pacific coast (González and Correa, 2001; Fig. 3).

Along both coasts, especially the Caribbean, crustal kinematics get exacerbated by mud diapirism upwarping the littoral fringe and continental shelf (Carvajal, 2016; Naranjo-Vesga et al., 2020). Mud diapirism manifests on- and offshore by dome-type features, mud, and gas emissions (Vivas-Narváez, 2019). For example, along Minuto de Dios, north of Arboletes town (Fig. 2), alongshore cliff-top elevations evidence an inclination toward the central mud volcano (Paniagua-Arroyave et al., 2018).

With this complex geologic framework in mind, we now explore the evidence of RSL changes along the NW South American coast.

Evidence of Holocene RSL changes along NW South America

Local studies have documented SL indicators in emerged terraces (Page, 1982; Martínez et al., 2010), sequences of organic deposits in coastal lagoons (González, 2017), drowned (Page and James, 1981), and raised beaches (González et al., 2014). These studies highlight tectonic motion as the most prominent process and main unknown. Despite the effort to document such SL indicators in remote locations, these records are difficult to analyze regarding postglacial RSL change, given the knowledge gap in tectonics and local GIA.

Considering the local GIA in the Caribbean, recent studies suggest that intermediate-field processes control the Holocene RSL change (González, 2017; Khan et al., 2017). These analyses used a set of available SL indicators and state-of-the-art statistical models to elucidate the RSL variability. They suggest that RSL has been rising during the past millennia due to the influence of the Laurentide forebulge collapse that induces crustal lowering with a constant ocean level. However, we hypothesize that the Caribbean coast of NW South America is within the limit of influence from such intermediate effects. We support our argument by highlighting several dubious assumptions in other authors’ work. The inconsistencies include: (1) assuming the central Caribbean coast does not exhibit tectonic deformation (González, 2017) and (2) including SL indicators from an island several hundreds of kilometers away from the Caribbean coast of NW South America (Khan et al., 2017).

We argue against the dominance of intermediate-field effects along the Caribbean, in line with recent evidence on crustal deformation. There is an ongoing discussion in favor of a “crustal block” model for the Caribbean coast of NW South America (Gómez-Álvarez, 2022), aligned with seminal contributions for the Pacific (Correa, 1996). In that vein, a recent study combined in situ instrumentation and remote sensing to propose that compressive tectonics controls RSL rise near Cartagena City by crustal subsidence (Restrepo-Ángel et al., 2021). In other words, tectonics can also lower the crust (and raise SL) along the Caribbean and Pacific coasts of NW South America. Therefore, SL indicators from coastal lagoons (e.g., Urrego et al., 2013) ought to be the submerged counterparts of SL indicators from emerged morphology, such as raised corals or marine terraces (e.g., Martínez et al., 2010).

The crustal block model can also explain RSL changes in the Pacific. Seminal works demonstrated the interdependence of hydrodynamics, sediment supply/accommodation, and local crustal structures, with coastal morphology manifesting in active and inactive (paleo) cliffs, bays, and deltas (Correa, 1996). Furthermore, co-seismic motions add to the crustal block kinematics, such that Holocene RSL indicators above the present mean SL can appear even without postglacial GIA far-field effects (Page and James, 1981; González et al., 2014). Mud diapirism along the Pacific coast should also control RSL, although the evidence is limited (Martínez and López Ramos, 2010).

We highlight two issues surrounding the current knowledge of postglacial RSL changes along the NW South American coasts. First, the results that support a dominance of intermediate-field effects depend on RSL indicators biased toward coastal lagoons irrespective of local tectonics. Second, tectonics creates a complex terrain deformation that prevents a simplified view of continental vertical change compared with local SL indicators, that is, inferring regional RSL changes from discrete locations. Therefore, we attempt to resolve these issues by simulating postglacial RSL changes with the best-suited Earth rheological models and analyze our simulations using the available SL indicators.

Modeling postglacial RSL in NW South America

Regional studies calibrated GIA models with local SL indicators for portions of the South American and Caribbean passive margins to elucidate the postglacial meltwater contribution to RSL changes (Milne et al., 2005; Milne and Peros, 2013). These studies also explained contributions from the ocean, ice, and rotational components to Holocene RSL change, with far-field effects influencing all South American coastlines and intermediate-field effects dominating the Caribbean. A critical part of these studies was finding the Earth rheology that explained the RSL change observations.

Here, we analyze the spatial variability in Holocene RSL changes along NW South America by exploring Earth rheological models in numerical simulations of postglacial RSL. We then compare the simulations with SL indicators to analyze the influence of crustal kinematics and solid Earth rheology. We emulate a recent contribution that examined Holocene RSL in the tectonically active Chilean region by reassessing SL indicators and contrasting them with GIA modeling (Garrett et al., 2020). We highlight how subduction styles and crustal deformation decisively control postglacial RSL change along the NW South American coasts.

METHODS

GIA modeling of postglacial RSL change

We modeled postglacial RSL changes along the NW South American coast (Figs. 4 and 5, Table 1) by solving the sea-level equation (SLE). The SLE, first introduced by Farrell and Clark (1976), provides a self-consistent quantitative description of the physical interactions between the cryosphere, oceans, and solid Earth in response to the melting of ice sheets. We obtained numerical solutions of the SLE with the open-source SELEN⁴ code (Spada and Melini, 2019), which is based on the pseudospectral method and considers the horizontal migration of shorelines (Peltier, 2004) and rotational feedback (Milne and Mitrovica, 1998). The numerical solution has been obtained on a global icosahedron grid (Tegmark, 1996) with resolution R = 44, corresponding to a spatial resolution of about 90 km on the Earth's surface.

Figure 4. Location of relative sea-level sites along Colombia's Caribbean Coast (red circles). Coastline data from the Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG) and elevation data from ETOPO1 1 arc-minute global relief model, both from the M_Map routine (Pawlowicz, 2020). Note that the selection of these sites is arbitrary, following the general location of sea-level indicators and landmarks.

Figure 5. Location of relative sea-level sites along Colombia's Pacific Coast (blue diamonds). Coastline data from the Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG) and elevation data from ETOPO1 1 arc-minute global relief model, both from the M_Map routine (Pawlowicz, 2020). Note that the selection of these sites is arbitrary, following the general location of sea-level indicators and landmarks.

Table 1. Relative sea-level curve locations (geographic coordinates) for SELEN⁴ simulations for the northwestern South American coast.

^a These locations are arbitrarily selected according to geographic landmarks (deltas, bays, etc.). Some sites include nearby large cities and towns (following dash and in italics) for reference (e.g., Minuto de Dios near Arboletes town).

In our solutions of the SLE, the spatial–temporal evolution of ice sheets is assumed to follow the ICE-6G_C GIA model of Peltier et al. (2015), which has been converted to a piecewise constant time history, as required by SELEN⁴, with a constant time step of 500 yr from 26 ka to 0 ka (present). The Earth's internal structure is assumed to be spherically symmetric, with incompressible, linear (Maxwell-type) mantle rheology, an elastic lithosphere, and an inviscid fluid core. Consistent with the ICE-6G_C GIA model, the radial structure is assumed to follow the VM5 rheological profile (Peltier and Drummond, 2008).

To represent the different subduction styles along the Caribbean and Pacific coasts, we explored two different mantle viscosity scenarios: a standard viscosity (that accounts for global data, VM5i from Spada and Melini [2019]; Table 2) and a high-viscosity model according to previous studies that we called “VM5h” (Milne et al., 2005; Milne and Peros, 2013; Table 3). We applied the same ice history ICE-6G_C from Peltier et al. (2015) with both models.

Table 2. Profiles of density, rigidity, and viscosity for the VM5i (standard mantle viscosity) model.^a

^a Rheological parameters after Spada and Melini (2019).

^b LT, lithosphere; UM, upper mantle; TZ, transition zone; LM, lower mantle; CO, core.

Table 3. Profiles of density, rigidity, and viscosity for the VM5h (high mantle viscosity) model.^a

^a We propose the viscosity profile after Milne and Peros (2013) and Milne et al. (2005).

^b LT, lithosphere; UM, upper mantle; TZ, transition zone; LM, lower mantle; CO, core.

Note that only our standard (VM5i) model reconciles simulations with global RSL data for a given ice thickness history (model ICE-6G_C of Peltier et al., 2015), as applied elsewhere (see Spada and Melini [2022] and others). However, changing the viscosity profile without varying the ice thickness history creates a mismatch with local RSL observations. Our high-viscosity model does not intend to reconcile simulations with global SL indicators. Instead, we want to analyze the GIA simulations with a plausible Earth rheology to compare with the available RSL indicators for the NW South American coasts.

NW South American RSL indicators

We compiled SL indicators available for the NW South American coasts, considering, whenever possible, their “fundamental characteristic attributes” as proposed in the literature: geographic location, age of formation, elevation to contemporary tidal datum, and relationship to RSL (Khan et al., 2019). SL indicators comprise coastal sediments (for correlations with marine–terrestrial transitions within the stratigraphic record) and geomorphic evidence (corals, terraces–cliffs, beach ridges), with corresponding dating and laboratory techniques to quantify ages and characterize environments. The compiled SL indicators fall into geomorphological (emerged terraces), coastal sediment, and coral reef categories (Shennan et al., 2015).

We included 62 indicators (54 for the Caribbean and 8 for the Pacific) comprising coastal lagoon sediments, beach sediments, mollusk shells, and corals, with ages given by radiocarbon and optically stimulated luminescence (OSL) dating. Given the heterogeneous methodologies used and incomplete descriptions of those methodologies, we included time and vertical location uncertainties according to the literature (Engelhart et al., 2011; Engelhart and Horton, 2012; Khan et al., 2017).

To use a standard temporal frame, we quantified the calibrated years for the 59 radiocarbon dates available using the MatCal v. 3.1 routine in MATLAB (Lougheed and Obrochta, 2016), with Marine20 and IntCal20 calibration curves, ΔR = 19 ± 23 ¹⁴C yr (Caribbean reservoir effect; Martínez et al., 2010) and R(t) = 198 ± 163 ¹⁴C yr (Pacific; Reimer and Reimer, 2001). We referenced all dates to 1950 CE as “present,” including beach sediments (three dates) at Bajo Baudó (Pacific site 5) dated by OSL from González et al. (2014).

In terms of radiocarbon ages and simulated RSL changes, we have a total of 12 dates related to Minuto de Dios (Caribbean site 2), 5 to the Gulf of Morrosquillo (Caribbean site 3), 7 to Manzanillo del Mar (Caribbean site 4), 6 to the Magdalena River delta (Caribbean site 5), 4 to the Ranchería River delta (Caribbean site 6), 2 to Utría Cove (Pacific site 6), and 5 to Solano Bay (Pacific site 7).

RESULTS

Modeled postglacial RSL changes

Figures 6 and 7 show the modeled RSL changes for NW South American (Caribbean and Pacific) coasts from the standard (VM5i) and high (VM5h) viscosity scenarios. We include postglacial (Figs. 6A and 7A) and Holocene changes (Figs. 6B and 7B).

Figure 6. Modeled relative sea-level curves from SELEN⁴ for the two mantle viscosity scenarios (VM5i, standard mantle viscosity, light red; and VM5h, high mantle viscosity, dark red) for locations along the Caribbean coast of Colombia (A) since the Last Glacial Maximum and (B) Holocene (−11.7 ka to present; Walker et al., 2018). Dashed vertical lines mark age transitions within the Holocene epoch at −8.236 ka (Greenlandian–Northgrippian) and −4.2 ka (Northgrippian–Meghalayan).

Figure 7. Modeled relative sea-level curves from SELEN⁴ for the two mantle viscosity scenarios (VM5i, standard mantle viscosity, light blue; and VM5h, high mantle viscosity, dark blue) for locations along the Pacific coast of Colombia (A) since the Last Glacial Maximum and (B) Holocene (−11.7 ka to present; Walker et al., 2018). Dashed vertical lines mark age transitions within the Holocene epoch at −8.236 ka (Greenlandian–Northgrippian) and −4.2 ka (Northgrippian–Meghalayan).

At millennial timescales, variations in time and space are insignificant on both coasts. Both viscosity models predicted RSL ~100 m below the present level at ~−26 ka. Relative SL increased from ~ −100 m to ~ −80 m between −26 and −14 ka, then increased rapidly to ~ −60 m between −14 and −13 ka, linked to the Melt Water Pulse 1A (Liu and Milliman, 2004; Liu et al., 2015). Then, both models predict an increase in RSL until −6 ka and a relatively constant RSL from −6 ka to the present.

The difference in Holocene RSL changes between models is evident from the mid-Northgrippian (6 ka) to the present. The standard model predicts a high stand and RSL fall, whereas the high-viscosity model predicts either an RSL rise or a still stand. Overall, the RSL changes for the standard model represent the response to hydro-isostatic effects along continental shelves in the far field of ice sheets (Clark et al., 1978). In contrast, the high-viscosity model results include far- and intermediate-field effects (Engelhart et al., 2009).

Modeled Holocene RSL changes

Figures 8 and 9 show the Mid- and Late Holocene RSL changes along the NW South American coast. For the Caribbean (Fig. 8A), the standard model predicts an RSL high stand of ~2 m around 6 ka at the southern locations (e.g., Minuto de Dios, site 2, and Gulf of Morrosquillo, site 3), similar to Zone VI of continental shorelines from Clark et al. (1978). This peak is on the order of centimeters at the northern locations (e.g., Caribbean site 7, Guajira Peninsula, and site 5, Magdalena River Delta). The high stand also shifts in time, from ~6 ka in the south to ~3 ka in the north. These results are consistent with extensive morphological observations that suggest coastal change depends on a regional RSL change related to postglacial hydro-isostatic effects (Correa, 1990; Correa and Vernette, 2004; Correa et al., 2005, 2007, 2016; Correa and Morton, 2010a). This scenario implies a Late Holocene RSL fall of ~2 m in 6 ka, or 0.33 mm/yr, and marine regression at the southern sites. This framework would result in Holocene SL indicators above the present mean SL, such as raised beaches and coastal paleo-cliffs (Dougherty et al., 2019).

Figure 8. Modeled relative sea-level (RSL) curves during the Late Holocene (−8.236 ka to present) for locations along the Colombian Caribbean coast: (A) standard viscosity scenario, VM5i (light red); and (B) high-viscosity scenario, VM5h (dark red). The vertical dashed line shows the Northgrippian–Meghalayan age transition at −4.2 ka.

Figure 9. Modeled relative sea-level (RSL) curves during the Late Holocene (−8.236 ka to present) for the Colombian Pacific coast: (A) standard viscosity scenario, VM5i (light blue); and (B) high-viscosity scenario, VM5h (dark blue). The vertical dashed line shows the Northgrippian–Meghalayan transition at −4.2 ka.

The high-viscosity simulations for the Caribbean (Fig. 8B) suggest that during the Meghalayan age (4.2 ka to present), RSL changes have depended on intermediate-field effects related to the Laurentide proglacial forebulge collapse, as in Zone IV of marine submergence from Clark et al. (1978). This scenario concurs with recent studies, wherein RSL indicators obtained at (presumably) tectonically stable sites suggest a gradual marine transgression during the last 4 ka (Vélez et al., 2014; González, 2017). Like the standard model, simulations predicted higher values for southern locations (e.g., Gulf of Morrosquillo, Caribbean site 3), in contrast to northern areas with RSL curves ~0.5 m below the southern zones (e.g., Guajira Peninsula, Caribbean site 7). These curves suggest a Late Holocene RSL still stand at the northern sites and a marine transgression of ~1 m in the last ~2 ka (+0.5 mm/yr) at the southern sites.

Holocene RSL curves for the Pacific coast show spatially homogeneous Mid-Holocene high stands for both models (Fig. 9). For the standard model (more representative of the Pacific subduction), we predicted an ~2 m high stand (height variability of ~0.5 m among sites), with an RSL fall of ~2 m in 6 ka (0.33 mm/yr) during the Late Holocene. On the other hand, we predicted a Mid-Holocene high stand of 0.50 m for the high-viscosity model (height variability of ~0.5 m among sites), with a variable RSL fall: −0.5 m in 6 ka (0.083 mm/yr) at the Naya River mouth (site 3) compared with −0.10 m in 2 ka (0.05 mm/yr) at the Patía River delta (site 2). We did not predict an RSL rise during the Late Holocene for the high-viscosity scenario, in contrast to the results for the Caribbean.

Compiled SL indicators

We now compare our GIA simulations with the compiled SL indicators. On the Caribbean coast, SL indicators from marine fossils are typically above modeled RSL curves for both mantle models (Fig. 10B–D). These indicators correspond to emerged marine terraces or corals from the Meghalayan age (4.2 ka to present) (Page, 1982; Correa and Paniagua-Arroyave, 2016) that result from the combination of postglacial hydro-isostasy, tectonics (Martínez et al., 2010; Restrepo-Ángel et al., 2021), and mud diapirism (Naranjo-Vesga et al., 2020).

Figure 10. Comparison between modeling results and sea-level indicators from the literature for the Caribbean coast of NW South America (Page, 1982; Martínez et al., 2010; Urrego et al., 2013; Vélez et al., 2014; González, 2017; please see our Supplementary Material at http://dx.doi.org/10.17632/7nhpbhvfnz.2). (A) Location of sea-level indicators (gray triangles) indicating the subplot in black letters; simulated relative sea-level curves for VM5i and VM5h models for (B) site 2 near Minuto de Dios; (C) site 3 near the Gulf of Morrosquillo; (D) site 4 near Manzanillo del Mar; (E) site 5 near the Magdalena River delta; and (F) site 6 near the Ranchería River delta. All dates have 1950 CE as the time 0 ka (present). Vertical dashed lines indicate the Northgrippian–Meghalayan transition at −4.2 ka. We distinguish between index points and marine-limiting date indicators following Khan et al. (2017) and references therein.

Conversely, Caribbean indicators from coastal lagoons are below the RSL predictions (Fig. 10E and F). For example, near the Gulf of Morrosquillo (Caribbean site 3 in Fig. 4A), González (2017) reported RSL values of −1.6 m at 4.3 ka and −0.4 m at 7.3 ka (Fig. 10C) that follow the marine transgression predicted by our high-viscosity model at the northern Caribbean sites (Fig. 8B).

On the Pacific coast, previous studies reported ~2 m raised beaches at Terco and Termales near the Corrientes Cape (Pacific site 5) (González et al., 2014; Fig. 11B). The active tectonics, as well as GIA, can explain the raised beaches. For the dates close to 3 ka found by colleagues, our high mantle viscosity model predicts raised beaches ~1.0 m above the current mean SL.

Figure 11. Comparison between our modeling results and sea-level indicators from the literature for the Pacific coast of NW South America (Page and James, 1981; Jaramillo and Bayona, 2000; González et al., 2014; data in our Supplementary Material at http://dx.doi.org/10.17632/7nhpbhvfnz.2). (A) Location of sea-level indicators (gray triangles) with letters indicating the subplot in black letters; simulated relative sea-level curves for VM5i and VM5h models for (B) site 5 near Bajo Baudó; (C) site 6 near Utría Cove; and (D) site 7 near Solano Bay. We distinguish between index points and marine- limiting date indicators following Khan et al. (2017) and references therein.

Finally, the SL indicators of coastal lagoons in the Pacific are below our predictions (Fig. 11C and D). Previous studies interpreted submerged SL indicators at Solano Bay as representative of co-seismic subsidence (Page and James, 1981), in contrast to the hydro-isostatic mechanisms. We argue that emergence or subsidence occurs on both coasts: the coastal fringe can either emerge or subside, with correspondent RSL fall or rise due to faulting–folding and block stacking: the “crustal block” model.

Our comparison of GIA simulations and SL indicators provides two insights. First, tectonics seems more prominent than expected, as RSL can also fall by crust subsidence. Second, GIA simulations suggest the Pacific and Caribbean differ in which GIA process dominates, with the transition from intermediate-field to far-field effects occurring on the Caribbean coast.

DISCUSSION

This work compares simulations of postglacial RSL changes with available SL indicators along the NW South American coasts to analyze the spatial variability and drivers of RSL. We applied two scenarios of mantle rheology (standard and high viscosities) and compared the results with published SL indicators. SL indicators are above or below our predictions according to compressive tectonics that produced emergence or subsidence. We now discuss our simulations regarding GIA mechanisms (intermediate- and far-field processes). Then, we discuss how tectonics influences RSL change through GIA response (solid Earth rheology and subduction styles) and solid Earth deformation (crustal block model).

Modeling postglacial RSL: ocean versus ice loadings

We consider the main uncertainties associated with GIA RSL modeling: Earth's structure and the modeling approach (Melini and Spada, 2019). From seminal works, it is well known that the postglacial ice-sheet waning influenced low-latitude coastlines by water flux from equatorial regions to zones exhibiting forebulge collapse (Mitrovica and Peltier, 1991; Mitrovica and Milne, 2002). This collapse induced a Holocene RSL fall, because the water flux from far- to intermediate-field regions drove a long-term SL fall after the instantaneous rise caused by ice-sheet melting (i.e., “ocean siphoning”). In addition, the added ocean water on continental shelves “tilted” the continents toward the sea and exposed portions of the inner shelves (i.e., “continental margin levering”) (Clark et al., 1978). Therefore, postglacial RSL models usually differentiate two GIA loadings: “ice” from the forebulge collapse and “ocean” from continental shelf hydro-isostasy.

Regional simulations with a standard solid Earth model suggest that the ocean loading happens along South America's coasts, with mid-Northgrippian (7 ka) ~2 m high stands along the southern Caribbean and northern Pacific coasts of NW South America (Milne et al., 2005). Spatially, the ocean component of these simulations includes a gradient in RSL change perpendicular to the general coastline orientation, with a zero RSL change coinciding with the coastline near the Magdalena River delta. Therefore, the model does not predict a high stand north of the Magdalena delta, because the ocean component (hydro-isostasy) is negative. Conversely, the ice-loading component is negative everywhere in northern South America (Fig. 5 in Milne et al., 2005). This contrast relates to the dichotomy we are exploring: whether the ocean or ice loadings dominated during the Meghalayan age (4.2 ka to present).

Considering the ocean–ice loading dichotomy, SL indicators suggest contrasting origins. For example, the RSL high stands related to the ocean-loading mechanism explain indicators in the southern Caribbean's marine terraces (Page, 1982; Correa et al., 2007). However, they do not concur with recent observations from coastal lagoons that suggest the dominance of the ice-loading effects (Urrego et al., 2013; González, 2017; Khan et al., 2017).

Because the SLE numerically solved by SELEN⁴ is nonlinear, we cannot separate the ocean from ice loadings. However, we can distinguish them in the resulting SL curves to assess the intermediate- to far-field influence. First, results from the standard model align with the control of the ice-loading (intermediate-field) mechanism north of our northernmost stations (e.g., Caribbean site 7, Guajira Peninsula, and site 5, Magdalena River Delta). There, the hydro-isostatic factor accounts for a mid-Northgrippian (6 ka) RSL high stand on the order of centimeters. In other words, our northernmost stations are close to the southern limit of the intermediate-field effects from the Laurentide forebulge collapse. Thus, the ocean-loading (far-field) effects are negligible in the northern Caribbean and more prominent in the south.

From our high-viscosity model, representative of the Caribbean according to Milne and Peros (2013), we propose that the transition from far- to intermediate-field effects is located between the Gulf of Morrosquillo (site 3) and Manzanillo del Mar (site 4). As this result contradicts what is currently accepted by the scientific community (Khan et al., 2017), we respectfully highlight some oversights that led to an incorrect interpretation of the Caribbean's SL indicators in previous works. Khan and colleagues used two Caribbean records to calibrate a statistical model and concluded that intermediate-field effects dominate the Caribbean (Khan et al., 2017). However, one of the records corresponds to San Andres Island, located ~700 km from continental South America (González et al., 2010). This record hardly represents the RSL variability along the Caribbean coast of NW South America. The second record corresponds to a coastal lagoon (Urrego et al., 2013), ~130 km NE from an emerged SL not considered due to tectonic “contamination” (Martínez et al., 2010). Considering the crustal block model, a question arises: Why do our colleagues consider coastal lagoons as tectonically stable sites?

Furthermore, although it was not applied in the modeling, Khan et al. (2017) discussed the Pacific RSL record by Jaramillo and Bayona (2000) as a Caribbean SL indicator. According to our estimates and recent unpublished analyses (Gómez-Álvarez, 2022), the coastal lagoon records of Jaramillo and Bayona and Urrego et al. (2013) may lie in subsiding coastal fringes like other sectors along the Caribbean coast (e.g., Cartagena Bay) (Restrepo-Ángel et al., 2021). In this case, SL indicators are below modern mean SL because of crustal block subsidence (see our discussion on tectonics: “Reconciling SL Indicators: The Role of Tectonics”).

With the standard viscosity model, we predict mid-Holocene (6 ka) high stands on the order of meters along the southern Caribbean and Pacific coasts that partially explain the emerged SL indicators. However, the model presents neither the submerged SL indicators nor their proximity to emerged indicators (e.g., Caribbean site 4 vs. 5). Exploring the response of a 3D solid Earth structure might reconcile this inconsistency, as the GIA response would be linked to subduction styles and laterally heterogeneous rheology.

3D solid Earth structure along coastal Colombia

The influence of solid Earth's rheology on postglacial GIA along the NW South America coast remains poorly understood. Optimizing solutions to the SLE for the northern Caribbean with SL indicators in Cuba resulted in relatively high mantle viscosities (Milne and Peros, 2013). Given the spherically symmetric model of SELEN⁴, we applied a similar mantle model in our VM5h simulations to represent the GIA response along the Caribbean coast (high viscosity), in contrast to the Pacific coast (standard–low viscosity) (e.g., Creveling et al., 2017). However, representing the complicated tectonic setting of NW South America may require a laterally varying, 3D solid Earth structure (Latychev et al., 2005; Hay et al., 2017; Mohammadzaheri et al., 2021). In the 3D case, the GIA response depends on lateral and vertical variations in crust thickness and mantle viscosity (Gomez et al., 2018; Thompson et al., 2023). In the far field, a 3D rheology implies variable controls of the gravitational and equatorial siphoning, with the continental levering relatively more influenced by local characteristics (Peak et al., 2022).

Considering the vertical rheology, northern South America exhibits a relatively thin crust of ~45 km (Feng et al., 2007) and three slabs with different subduction angles (Vargas and Mann, 2013; Idárraga-García et al., 2016). These features translate into low mantle viscosities along active plate boundaries with a steep subduction slab. Also, we expect high mantle viscosities along active plate boundaries with a shallower (flat) subduction slab. Global analogies include the high viscosity of the flat slab in Barbados (Austermann et al., 2013) and low viscosities in the active subduction zone of Alaska (Larsen et al., 2005; Lange et al., 2014).

A laterally variable mantle structure implies different relaxation times to surface load changes (Whitehouse, 2018). We expect such differences along the Pacific because of the latitudinally varying subduction: flat subduction south of Malpelo Island (Pacific site 3 near Naya River delta) and north of Solano Bay (Pacific site 5) and steep subduction centered at the San Juan River delta (Pacific sites 3 and 4) (Idárraga-García et al., 2016).

According to the subduction styles, the high mantle viscosity represents the Caribbean's flat slab and the Pacific subduction north of Solano Bay and south of the Naya River delta. Recent surface ice mass changes should influence the GIA response along these coasts. On the other hand, the low-viscosity model represents the steep-subduction region along the mid-Pacific coast. Emerged SL indicators because of equatorial siphoning should dominate along these coasts.

Reconciling SL indicators: the role of tectonics

Comparing RSL change simulations and SL indicators is customary, as it can shed light on the role of tectonics (Garrett et al., 2020). Seminal studies in tectonically stable sites linked emerged landforms to postglacial hydro-isostatic effects. For example, in Australia, GIA explains mid-Holocene high stands on the order of meters if observations are used to adjust a mantle viscosity model (Nakada and Lambeck, 1989). These results agree with SL indicators in South America, for example, the eastern coast of Brazil (Angulo et al., 2006) and the Rio de la Plata estuary (Prieto et al., 2017). Following the continental levering mechanism, these sites exhibit SL high stands >4 m (Mitrovica and Milne, 2002). For NW South America, previous works predicted a Holocene high stand of centimeters along the northern Caribbean coast, whereas predictions for the southern Caribbean proposed a ~3 m high stand (Clark et al., 1978; Page, 1982). Our results with the standard mantle structure concur with previous studies. On the other hand, results with the high-viscosity model predict a still stand south of the Gulf of Morrosquillo. Thus, in light of our results, tectonics are more significant than previously argued, especially on the southern Caribbean coast.

Active tectonics (including mud diapirism) shapes the NW South American Caribbean coast with the interplay of three major tectonic plates (Caribbean, Nazca, and South American) and two main crustal blocks (Panamá-Chocó and northern Andes) (Kellogg and Vega, 1995; Taboada et al., 2000; Cortés et al., 2005). The shallow subduction of the Caribbean plate beneath the South American plate at 20 mm/yr in the NE direction dominates coastal terrain deformation (Syracuse et al., 2016; Mora-Páez et al., 2019). These dynamics configure the Sinú-San Jacinto Deformed Belt in front of the subduction zone (cf. Fig. 1). Four (sites 2 to 5) of the seven analyzed sites are located within this deformational front.

The Gulf of Urabá (Caribbean site 1) is located on the eastern border of the Panamá-Chocó crustal block, at the limit of the northern Andes block. To the south, along the Atrato River valley, the Chocó block is limited by the Uramita and related faults, whereas to the north, these faults pass beneath the Gulf of Urabá and get dispersed within the Sinú Fold Belt and the North Panamá Thrust Belt. These active faults constitute thrust faults (with a minor left-lateral component) controlling the coastal zone's vertical deformation by subsidence and emergence.

Sites 2 to 5 (Minuto de Dios, Gulf of Morrosquillo, Manzanillo del Mar, and Magdalena River delta) are in the northern part of the northern Andes block, along the Sinú-San Jacinto Fold Belt. These terrains formed along folded sedimentary rocks imbricated with thrust faults and raised due to crust stacking (e.g., Vinnels et al., 2010). This deformation, which occurs along an extensive structure with many anisotropies, is poorly understood. However, the recent analyses along Cartagena Bay showed an RSL rise of ~7 mm/yr linked to coastal subsidence because of crustal block dynamics (Restrepo-Ángel et al., 2021). On the contrary, we can find the counterpart deformation that has produced coastal emergence at Manzanillo del Mar (Martínez et al., 2010; Caribbean site 4, Fig. 6D).

Considering GPS observations, preliminary analyses (Gómez-Álvarez, 2022) confirm the ~1 mm/yr of subsidence at the Gulf of Morrosquillo proposed by previous unpublished reports (Page, 1982). Assuming subsidence operated uniformly during the Meghalayan (4.2 ka to present), the SL indicators from González (2017) (Caribbean site 3) become +2.7 m at 4.3 ka and −0.4 m at 7.3 ka. These values concur with our standard model's 6 ka high stand (Fig. 8A).

More pronounced than in the Caribbean, tectonics controls the coastal morphology in the Pacific through active faulting, mud diapirism, and co-seismic subsidence. The morphology controlled by tectonics includes: (1) active cliffs and short rivers with relatively low sediment load, on which subsidence levels determine lagoon/bay morphology; (2) coastal paleo-cliffs (bluffs) and hills, flanked by faults toward the littoral fringe; (3) the San Juan River delta, deposited on an oceanward-dipping, gently sloping graben; and (4) the Patía River delta, a deltaic region of co-seismic subsidence compartments (Correa, 1996, p. 149).

Overall, the alongshore succession of thrust faults, mud diapirism (Caribbean and Pacific), and co-seismic motions (Pacific) drives RSL change as well as postglacial GIA. A variable continental level would result in a variable RSL change even for a relatively constant ocean level. We argue that these tectonic and structural responses provide the primary mechanism configuring RSL changes along the NW South American coast. In other words, contrary to what other studies suppose, none of the RSL indicators appear to result from a tectonically stable location.

What we missed: sediment isostasy

A significant unknown in RSL change along the NW South American coast relates to sediment isostasy. Sediment deposition and erosion affect RSL change as they vary the mass distribution and relative distance between the ocean and solid Earth surfaces (Dalca et al., 2013). In NW South America, this distribution depends on sediment transfer to the continental shelf in the basin–coastal zone continuum. As sediment erodes from the continent, it would reduce continental mass and imply an RSL fall (because of continental rising). In contrast, sediment accumulation on the continent leads to mass loading and RSL rise.

Despite relatively few applications of sediment isostasy, available studies highlight it as an effective mechanism of RSL change for continental shelves with appreciable sediment input. Seminal analyses at Karachi in the Arabian Sea, close to the Indus delta, proposed a postglacial RSL correction of ~7 m related to sediment isostasy. This correction implies a Late Holocene still stand instead of the high stand of ~3 m expected from ocean and ice loads (Fig. 6A in Ferrier et al., 2015).

Given the relatively large sediment load from the Caribbean and Pacific catchments, sediment isostasy promises to control RSL along NW South America. For example, fluvial sediment deposition has created arguably the most extensive delta systems along the Pacific Coast of North and South America (the Patía and San Juan river deltas) despite the narrow and high-energy shelf (Restrepo and López, 2008). Such delta progradation influences RSL through continental sediment redistribution (Dalca et al., 2013).

Another control in sediment isostasy is retention at floodplains. For example, floodplains prevent ~10% of fluvial sediment load from reaching the coastal zone in the depositional region (the Momposina Depression) within the Magdalena River catchment (Restrepo et al., 2006). Such deposition occurs on an area of ~25,000 km², translating into ~55 m of Holocene sediments (Latrubesse, 2015).

However, beyond floodplains, accumulation is complicated at NW South American deltas. For the Magdalena, jetties route fluvial sediments to the continental rise through a submarine canyon, preventing prodelta accumulation and shelf loading (Naranjo-Vesga et al., 2021). In this case, engineering structures modify delta morphodynamics, resulting in increased channel siltation and an imbalance in marine sediment fluxes (Restrepo et al., 2020; Paniagua-Arroyave and Nienhuis, 2022). This imbalance varies sediment redistribution and prodelta reworking of Late Holocene deposits, adding another source of uncertainty to RSL changes.

In postglacial timescales (tens of thousands of years), sediment deposited in NW South American floodplains would lower RSL through continental tilting. On the contrary, RSL would rise because of deposition on the continental shelf. A relatively recent publication suggests an RSL fall of ~15 m since the last interglacial (122 ka) at deltas along the Caribbean and Pacific of NW South America due to sediment isostasy (Pico, 2020). Also, crustal uplift can be associated with erosional unloading via sediment isostasy and tectonic uplift, resulting in another source of RSL fall (Ruetenik et al., 2020).

CONCLUSIONS

Based on the comparison of postglacial RSL change simulations (with a model that solves the gravitationally and topographically self-consistent SLE for a given mantle rheology) with published SL indicators for NW South America's Caribbean and Pacific coasts, we find that:

• • Far-field effects (equatorial siphoning and continental levering) have dominated along the Pacific coast during the Holocene, with high stands on the order of meters (−0.33 mm/yr of Holocene RSL fall).

• • Intermediate-field effects (related to the Laurentide forebulge collapse) were more prominent along the Caribbean coast, with a Late Holocene RSL rise in the north (+0.5 mm/yr of Holocene RSL change) and a still stand in the south (~0 mm/yr).

• • The change in influence between far- and intermediate-field effects occurs between the Gulf of Morrosquillo and Manzanillo del Mar along the Caribbean coast.

• • The lateral variability in Earth's rheology supports applying a GIA model with 3D mantle structure, including the influence of sediment isostasy.

• • According to the crustal blocks’ approximation, published SL indicators correspond to emerged or submerged sites because of faulting and folding, mud diapirism (Caribbean and Pacific), and co-seismic motions (Pacific).