INTRODUCTION
Predicting coastal change depends upon deciphering postglacial relative sea-level (RSL) variability and its influence on modern morphology (FitzGerald et al., Reference FitzGerald, Fenster, Argow and Buynevich2008; Passeri et al., Reference Passeri, Hagen, Medeiros, Bilskie, Alizad and Wang2015). Postglacial RSL changes control modern coastal environments, as they dictate nearshore sediment accommodation and supply (Everts, Reference Everts, Nummedal, Pilkey and Howard1987; Törnqvist et al., Reference Törnqvist, Jankowski, Li and González2020). 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, Reference Nicholls and Cazenave2010; Kopp et al., Reference Kopp, Hay, Little and Mitrovica2015; Nerem et al., Reference Nerem, Beckley, Fasullo, Hamlington, Masters and Mitchum2018; Nienhuis and van de Wal, Reference Nienhuis and van de Wal2021).
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., Reference Garrett, Melnick, Dura, Cisternas, Ely, Wesson, Jara-Muñoz and Whitehouse2020). 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., Reference Shadrick, Rood, Hurst, Piggott, Hebditch, Seal and Wilcken2022; Nienhuis et al., Reference Nienhuis, Kim, Milne, Quock, Slangen and Törnqvist2023).
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 SELEN4 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., Reference Kennedy, Oliver, Tamura, Murray-Wallace, Thom, Rosengren, Ierodiaconou and Augustinus2020; Nienhuis et al., Reference Nienhuis, Ashton, Edmonds, Hoitink, Kettner, Rowland and Törnqvist2020; Anfuso et al., Reference Anfuso, Postacchini, di Luccio and Benassai2021). 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, Reference Pirazzoli1996, p. 5; Milne et al., Reference Milne, Gehrels, Hughes and Tamisiea2009; Gregory et al., Reference Gregory, Griffies, Hughes, Lowe, Church, Fukimori and Gomez2019).
During the Holocene (last 11.7 ka; Walker et al., Reference Walker, Head, Berkelhammer, Björck, Cheng, Cwynar and Fisher2018), GIA controlling RSL depended on the location on Earth (Lambeck et al., Reference Lambeck, Antonioli, Anzidei, Ferranti, Leoni, Scicchitano and Silenzi2011; Kopp et al., Reference Kopp, Hay, Little and Mitrovica2015). 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., Reference Rovere, Stocchi and Vacchi2016). Therefore, ice-sheet melting produced GIA responses that depended on local conditions, that is, inner Earth structure and distance from the ice sheets (Whitehouse, Reference Whitehouse2018).
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., Reference Clark, Farrell and Peltier1978). 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, Reference Isla1989; Hein et al., Reference Hein, FitzGerald, de Menezes, Cleary, Klein and Albernaz2014; Cooper et al., Reference Cooper, Meireles, Green, Klein and Toldo2018). 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, Reference Mitrovica and Peltier1991).
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., Reference Khan, Ashe, Shaw, Vacchi, Walker, Peltier, Kopp and Horton2015). 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., Reference Khan, Ashe, Shaw, Vacchi, Walker, Peltier, Kopp and Horton2015). 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., Reference Khan, Ashe, Shaw, Vacchi, Walker, Peltier, Kopp and Horton2015; Fig. 3I). See also figure 8 in Yokoyama and Purcell (Reference Yokoyama and Purcell2021) 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 (Reference Nygren1950), Duque-Caro (Reference Duque-Caro1990a, Reference Duque-Caro1990b), Vernette et al. (Reference Vernette, Mauffret, Bobier, Briceno and Gayet1992), Kellogg and Vega (Reference Kellogg and Vega1995), Colmenares and Zoback (Reference Colmenares and Zoback2003), and IGAC and INGEOMINAS (2006), as well as more recent accounts by Cediel and Shaw (Reference Cediel and Shaw2019), Montes et al. (Reference Montes, Rodriguez-Corcho, Bayona, Hoyos, Zapata and Cardona2019), Bustamante et al. (Reference Bustamante, Archanjo, Cardona and Vervoort2016, Reference Bustamante, Cardona, Restrepo, Zapata, Beltrán-Triviño, Bustamante and Valencia2024), Restrepo et al. (Reference Restrepo, Bustamante, Cardona, Beltrán-Triviño, Bustamante, Chavarría and Valencia2021, Reference Restrepo, Bustamante, Cardona, Beltrán-Triviño and Valencia2023), and references therein.
(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. (Reference Montes, Rodriguez-Corcho, Bayona, Hoyos, Zapata and Cardona2019) and Mora-Paez et al. (Reference Mora-Páez, Kellogg, Freymueller, Mencin, Fernandes, Diederix and LaFemina2019). Plate motions relative to the north Andean block with slip rates in centimeters per year (cm/yr) after Kellogg and Dixon (Reference Kellogg and Dixon1990). Block velocity relative to the stable north Andean block (1 cm/yr) from Mora-Páez et al. (Reference Mora-Páez, Kellogg, Freymueller, Gómez and Pinilla–Pachon2020). Coastline data (A) from the Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG) from the M_Map routine (Pawlowicz, Reference Pawlowicz2020).
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, Reference Correa, Pereira, Cediel and Shaw2019). Both coasts can be considered collisional (Inman and Nordstrom, Reference Inman and Nordstrom1971). 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, Reference Correa, Morton and Bird2010a, Reference Correa, Morton and Bird2010b).
Geology frames six coastal fringes along the Caribbean (Correa and Morton, Reference Correa, Morton and Bird2010a; Correa and Pereira, Reference Correa, Pereira, Cediel and Shaw2019; 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., Reference Correa, Prüssmann-Uribe, Garrido-Escobar, Blanco-Libreros and Londoño-Mesa2016). 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, Reference Correa, Paniagua-Arroyave and Hermelin2016). 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, Reference Correa, Pereira, Cediel and Shaw2019).
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 (Reference Correa, Pereira, Cediel and Shaw2019). 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, Reference Pawlowicz2020).
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., Reference González, Shen and Mauz2014; Montes et al., Reference Montes, Rodriguez-Corcho, Bayona, Hoyos, Zapata and Cardona2019). 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, Reference Correa1996; Correa and Morton, Reference Correa, Morton and Bird2010b). 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, Reference Correa, Pereira, Cediel and Shaw2019).
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 (Reference Correa, Pereira, Cediel and Shaw2019). Coastline data from the GSHHG Database and elevation data from the ETOPO1 1 model, both from the M_Map routine (Pawlowicz, Reference Pawlowicz2020).
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., Reference Mora-Páez, Kellogg, Freymueller, Mencin, Fernandes, Diederix and LaFemina2019). 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, Reference Martínez and López Ramos2010; Carvajal, Reference Carvajal and Hermelin2016) 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, Reference González and Correa2001; Fig. 3).
Along both coasts, especially the Caribbean, crustal kinematics get exacerbated by mud diapirism upwarping the littoral fringe and continental shelf (Carvajal, Reference Carvajal and Hermelin2016; Naranjo-Vesga et al., Reference Naranjo-Vesga, Ortiz-Karpf, Wood, Jobe, Paniagua-Arroyave, Shumaker, Mateus-Tarazona and Galindo2020). Mud diapirism manifests on- and offshore by dome-type features, mud, and gas emissions (Vivas-Narváez, Reference Vivas-Narváez2019). 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., Reference Paniagua-Arroyave, Correa, Anfuso and Adams2018).
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, Reference Page1982; Martínez et al., Reference Martínez, Yokoyama, Gomez, Delgado, Matsuzaki and Rendon2010), sequences of organic deposits in coastal lagoons (González, Reference González2017), drowned (Page and James, Reference Page and James1981), and raised beaches (González et al., Reference González, Shen and Mauz2014). 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, Reference González2017; Khan et al., Reference Khan, Ashe, Horton, Dutton, Kopp, Brocard and Engelhart2017). 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, Reference González2017) and (2) including SL indicators from an island several hundreds of kilometers away from the Caribbean coast of NW South America (Khan et al., Reference Khan, Ashe, Horton, Dutton, Kopp, Brocard and Engelhart2017).
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, Reference Gómez-Álvarez2022), aligned with seminal contributions for the Pacific (Correa, Reference Correa1996). 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., Reference Restrepo-Ángel, Mora-Páez, Díaz, Govorcin, Wdowinski, Giraldo-Londoño, Tosic, Fernández, Paniagua-Arroyave and Duque-Trujillo2021). 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., Reference Urrego, Correa-Metrio, González, Castaño and Yokoyama2013) 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., Reference Martínez, Yokoyama, Gomez, Delgado, Matsuzaki and Rendon2010).
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, Reference Correa1996). 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, Reference Page and James1981; González et al., Reference González, Shen and Mauz2014). Mud diapirism along the Pacific coast should also control RSL, although the evidence is limited (Martínez and López Ramos, Reference Martínez and López Ramos2010).
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., Reference Milne, Long and Bassett2005; Milne and Peros, Reference Milne and Peros2013). 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., Reference Garrett, Melnick, Dura, Cisternas, Ely, Wesson, Jara-Muñoz and Whitehouse2020). 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 (Reference Farrell and Clark1976), 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 SELEN4 code (Spada and Melini, Reference Spada and Melini2019), which is based on the pseudospectral method and considers the horizontal migration of shorelines (Peltier, Reference Peltier2004) and rotational feedback (Milne and Mitrovica, Reference Milne and Mitrovica1998). The numerical solution has been obtained on a global icosahedron grid (Tegmark, Reference Tegmark1996) with resolution R = 44, corresponding to a spatial resolution of about 90 km on the Earth's surface.
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, Reference Pawlowicz2020). Note that the selection of these sites is arbitrary, following the general location of sea-level indicators and landmarks.
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, Reference Pawlowicz2020). Note that the selection of these sites is arbitrary, following the general location of sea-level indicators and landmarks.
Relative sea-level curve locations (geographic coordinates) for SELEN4 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. (Reference Peltier, Argus and Drummond2015), which has been converted to a piecewise constant time history, as required by SELEN4, 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, Reference Peltier and Drummond2008).
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., Reference Milne, Long and Bassett2005; Milne and Peros, Reference Milne and Peros2013; Table 3). We applied the same ice history ICE-6G_C from Peltier et al. (Reference Peltier, Argus and Drummond2015) with both models.
Profiles of density, rigidity, and viscosity for the VM5i (standard mantle viscosity) model.a
a Rheological parameters after Spada and Melini (Reference Spada and Melini2019).
b LT, lithosphere; UM, upper mantle; TZ, transition zone; LM, lower mantle; CO, core.
Profiles of density, rigidity, and viscosity for the VM5h (high mantle viscosity) model.a
a We propose the viscosity profile after Milne and Peros (Reference Milne and Peros2013) and Milne et al. (Reference Milne, Long and Bassett2005).
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., Reference Peltier, Argus and Drummond2015), as applied elsewhere (see Spada and Melini [Reference Spada and Melini2022] 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., Reference Khan, Horton, Engelhart, Rovere, Vacchi, Ashe, Törnqvist, Dutton, Hijma and Shennan2019). 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., Reference Shennan, Long and Horton2015).
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., Reference Engelhart, Horton and Kemp2011; Engelhart and Horton, Reference Engelhart and Horton2012; Khan et al., Reference Khan, Ashe, Horton, Dutton, Kopp, Brocard and Engelhart2017).
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, Reference Lougheed and Obrochta2016), with Marine20 and IntCal20 calibration curves, ΔR = 19 ± 23 14C yr (Caribbean reservoir effect; Martínez et al., Reference Martínez, Yokoyama, Gomez, Delgado, Matsuzaki and Rendon2010) and R(t) = 198 ± 163 14C yr (Pacific; Reimer and Reimer, Reference Reimer and Reimer2001). 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. (Reference González, Shen and Mauz2014).
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).
Modeled relative sea-level curves from SELEN4 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., Reference Walker, Head, Berkelhammer, Björck, Cheng, Cwynar and Fisher2018). Dashed vertical lines mark age transitions within the Holocene epoch at −8.236 ka (Greenlandian–Northgrippian) and −4.2 ka (Northgrippian–Meghalayan).
Modeled relative sea-level curves from SELEN4 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., Reference Walker, Head, Berkelhammer, Björck, Cheng, Cwynar and Fisher2018). 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, Reference Liu and Milliman2004; Liu et al., Reference Liu, Milne, Kopp, Clark and Shennan2015). 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., Reference Clark, Farrell and Peltier1978). In contrast, the high-viscosity model results include far- and intermediate-field effects (Engelhart et al., Reference Engelhart, Horton, Douglas, Peltier and Törnqvist2009).
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. (Reference Clark, Farrell and Peltier1978). 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, Reference Correa1990; Correa and Vernette, Reference Correa and Vernette2004; Correa et al., Reference Correa, Alcántara-Carrió and González R2005, Reference Correa, Acosta and Bedoya2007, Reference Correa, Prüssmann-Uribe, Garrido-Escobar, Blanco-Libreros and Londoño-Mesa2016; Correa and Morton, Reference Correa, Morton and Bird2010a). 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., Reference Dougherty, Thomas, Fogwill, Hogg, Palmer, Rainsley and Williams2019).
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.
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. (Reference Clark, Farrell and Peltier1978). 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., Reference Vélez, Escobar, Brenner, Rangel, Betancourt, Jaramillo, Curtis and Moreno2014; González, Reference González2017). 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, Reference Page1982; Correa and Paniagua-Arroyave, Reference Correa, Paniagua-Arroyave and Hermelin2016) that result from the combination of postglacial hydro-isostasy, tectonics (Martínez et al., Reference Martínez, Yokoyama, Gomez, Delgado, Matsuzaki and Rendon2010; Restrepo-Ángel et al., Reference Restrepo-Ángel, Mora-Páez, Díaz, Govorcin, Wdowinski, Giraldo-Londoño, Tosic, Fernández, Paniagua-Arroyave and Duque-Trujillo2021), and mud diapirism (Naranjo-Vesga et al., Reference Naranjo-Vesga, Ortiz-Karpf, Wood, Jobe, Paniagua-Arroyave, Shumaker, Mateus-Tarazona and Galindo2020).
Comparison between modeling results and sea-level indicators from the literature for the Caribbean coast of NW South America (Page, Reference Page1982; Martínez et al., Reference Martínez, Yokoyama, Gomez, Delgado, Matsuzaki and Rendon2010; Urrego et al., Reference Urrego, Correa-Metrio, González, Castaño and Yokoyama2013; Vélez et al., Reference Vélez, Escobar, Brenner, Rangel, Betancourt, Jaramillo, Curtis and Moreno2014; González, Reference González2017; 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. (Reference Khan, Ashe, Horton, Dutton, Kopp, Brocard and Engelhart2017) 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 (Reference González2017) 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., Reference González, Shen and Mauz2014; 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.
Comparison between our modeling results and sea-level indicators from the literature for the Pacific coast of NW South America (Page and James, Reference Page and James1981; Jaramillo and Bayona, Reference Jaramillo and Bayona2000; González et al., Reference González, Shen and Mauz2014; 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. (Reference Khan, Ashe, Horton, Dutton, Kopp, Brocard and Engelhart2017) 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, Reference Page and James1981), 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, Reference Melini and Spada2019). 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, Reference Mitrovica and Peltier1991; Mitrovica and Milne, Reference Mitrovica and Milne2002). 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., Reference Clark, Farrell and Peltier1978). 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., Reference Milne, Long and Bassett2005). 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., Reference Milne, Long and Bassett2005). 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, Reference Page1982; Correa et al., Reference Correa, Acosta and Bedoya2007). However, they do not concur with recent observations from coastal lagoons that suggest the dominance of the ice-loading effects (Urrego et al., Reference Urrego, Correa-Metrio, González, Castaño and Yokoyama2013; González, Reference González2017; Khan et al., Reference Khan, Ashe, Horton, Dutton, Kopp, Brocard and Engelhart2017).
Because the SLE numerically solved by SELEN4 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 (Reference Milne and Peros2013), 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., Reference Khan, Ashe, Horton, Dutton, Kopp, Brocard and Engelhart2017), 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., Reference Khan, Ashe, Horton, Dutton, Kopp, Brocard and Engelhart2017). However, one of the records corresponds to San Andres Island, located ~700 km from continental South America (González et al., Reference González, Urrego, Martínez, Polanía and Yokoyama2010). 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., Reference Urrego, Correa-Metrio, González, Castaño and Yokoyama2013), ~130 km NE from an emerged SL not considered due to tectonic “contamination” (Martínez et al., Reference Martínez, Yokoyama, Gomez, Delgado, Matsuzaki and Rendon2010). 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. (Reference Khan, Ashe, Horton, Dutton, Kopp, Brocard and Engelhart2017) discussed the Pacific RSL record by Jaramillo and Bayona (Reference Jaramillo and Bayona2000) as a Caribbean SL indicator. According to our estimates and recent unpublished analyses (Gómez-Álvarez, Reference Gómez-Álvarez2022), the coastal lagoon records of Jaramillo and Bayona and Urrego et al. (Reference Urrego, Correa-Metrio, González, Castaño and Yokoyama2013) may lie in subsiding coastal fringes like other sectors along the Caribbean coast (e.g., Cartagena Bay) (Restrepo-Ángel et al., Reference Restrepo-Ángel, Mora-Páez, Díaz, Govorcin, Wdowinski, Giraldo-Londoño, Tosic, Fernández, Paniagua-Arroyave and Duque-Trujillo2021). 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, Reference Milne and Peros2013). Given the spherically symmetric model of SELEN4, 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., Reference Creveling, Mitrovica, Clark, Waelbroeck and Pico2017). However, representing the complicated tectonic setting of NW South America may require a laterally varying, 3D solid Earth structure (Latychev et al., Reference Latychev, Mitrovica, Tromp, Tamisiea, Komatitsch and Christara2005; Hay et al., Reference Hay, Lau, Gomez, Austermann, Powell, Mitrovica, Latychev and Wiens2017; Mohammadzaheri et al., Reference Mohammadzaheri, Sigloch, Hosseini and Mihalynuk2021). In the 3D case, the GIA response depends on lateral and vertical variations in crust thickness and mantle viscosity (Gomez et al., Reference Gomez, Latychev and Pollard2018; Thompson et al., Reference Thompson, Creveling, Latychev and Mitrovica2023). 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., Reference Peak, Latychev, Hoggard and Mitrovica2022).
Considering the vertical rheology, northern South America exhibits a relatively thin crust of ~45 km (Feng et al., Reference Feng, van der Lee and Assumpção2007) and three slabs with different subduction angles (Vargas and Mann, Reference Vargas and Mann2013; Idárraga-García et al., Reference Idárraga-García, Kendall and Vargas2016). 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., Reference Austermann, Mitrovica, Latychev and Milne2013) and low viscosities in the active subduction zone of Alaska (Larsen et al., Reference Larsen, Motyka, Freymueller, Echelmeyer and Ivins2005; Lange et al., Reference Lange, Casassa, Ivins, Schröder, Fritsche, Richter, Groh and Dietrich2014).
A laterally variable mantle structure implies different relaxation times to surface load changes (Whitehouse, Reference Whitehouse2018). 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., Reference Idárraga-García, Kendall and Vargas2016).
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., Reference Garrett, Melnick, Dura, Cisternas, Ely, Wesson, Jara-Muñoz and Whitehouse2020). 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, Reference Nakada and Lambeck1989). These results agree with SL indicators in South America, for example, the eastern coast of Brazil (Angulo et al., Reference Angulo, Lessa and Souza2006) and the Rio de la Plata estuary (Prieto et al., Reference Prieto, Mourelle, Peltier, Drummond, Vilanova and Ricci2017). Following the continental levering mechanism, these sites exhibit SL high stands >4 m (Mitrovica and Milne, Reference Mitrovica and Milne2002). 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., Reference Clark, Farrell and Peltier1978; Page, Reference Page1982). 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, Reference Kellogg and Vega1995; Taboada et al., Reference Taboada, Rivera, Fuenzalida, Cisternas, Philip, Bijwaard, Olaya and Rivera2000; Cortés et al., Reference Cortés, Angelier and Colletta2005). 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., Reference Syracuse, Maceira, Prieto, Zhang and Ammon2016; Mora-Páez et al., Reference Mora-Páez, Kellogg, Freymueller, Mencin, Fernandes, Diederix and LaFemina2019). 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., Reference Vinnels, Butler, McCaffrey and Paton2010). 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., Reference Restrepo-Ángel, Mora-Páez, Díaz, Govorcin, Wdowinski, Giraldo-Londoño, Tosic, Fernández, Paniagua-Arroyave and Duque-Trujillo2021). On the contrary, we can find the counterpart deformation that has produced coastal emergence at Manzanillo del Mar (Martínez et al., Reference Martínez, Yokoyama, Gomez, Delgado, Matsuzaki and Rendon2010; Caribbean site 4, Fig. 6D).
Considering GPS observations, preliminary analyses (Gómez-Álvarez, Reference Gómez-Álvarez2022) confirm the ~1 mm/yr of subsidence at the Gulf of Morrosquillo proposed by previous unpublished reports (Page, Reference Page1982). Assuming subsidence operated uniformly during the Meghalayan (4.2 ka to present), the SL indicators from González (Reference González2017) (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, Reference Correa1996, 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., Reference Dalca, Ferrier, Mitrovica, Perron, Milne and Creveling2013). 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., Reference Ferrier, Mitrovica, Giosan and Clift2015).
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, Reference Restrepo and López2008). Such delta progradation influences RSL through continental sediment redistribution (Dalca et al., Reference Dalca, Ferrier, Mitrovica, Perron, Milne and Creveling2013).
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., Reference Restrepo, Kjerfve, Hermelin and Restrepo2006). Such deposition occurs on an area of ~25,000 km2, translating into ~55 m of Holocene sediments (Latrubesse, Reference Latrubesse2015).
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., Reference Naranjo-Vesga, Paniagua-Arroyave, Ortiz-Karpf, Wood, Jobe, Galindo, Shumaker and Mateus-Tarazona2021). In this case, engineering structures modify delta morphodynamics, resulting in increased channel siltation and an imbalance in marine sediment fluxes (Restrepo et al., Reference Restrepo, Orejarena-Rondón, Consuegra, Pérez, Llinas, Otero and Álvarez2020; Paniagua-Arroyave and Nienhuis, Reference Paniagua-Arroyave and Nienhuis2022). 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, Reference Pico2020). 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., Reference Ruetenik, Ferrier, Creveling and Fox2020).
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).
Acknowledgments
JFP-A acknowledges funding from the Vice-Presidency of Science, Technology, and Innovation at EAFIT University (award 952-000015), Florida State University Department of Earth, Ocean and Atmospheric Science, Utrecht University Department of Physical Geography, the Polar Earth Observing Network (POLENET), and the International Association of Cryospheric Sciences (IACS, to attend the 2019 Glacial Isostatic Adjustment Training School in Gavle, Sweden). The authors recognize supercomputing resources made available by Centro de Computación Científica Apolo at EAFIT University (http://www.eafit.edu.co/apolo), with the gracious support of Juan G. Lalinde-Pulido, Laura Sánchez-Córdoba, and Jacobo Monsalve-Guzmán. The authors sincerely thank Glenn A. Milne for sharing his modeling results and guiding the comparison with SL indicators. JFP-A acknowledges Iván D. Correa for inspiration and guidance regarding discussions of RSL changes along coastal Colombia. The authors acknowledge thoughtful reviews by Barbara Mauz, Colin Murray-Wallace, and Juan L. González on an earlier version of this article and editorial work by Lewis Owen, Karin E. Perring, Lucie Taylor, and Mary Safford Curioli. The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. This work contributes to IGCP Project 725 “Forecasting Coastal Change.”
Data Availability Statement
Supplementary material and data can be accessed at http://dx.doi.org/10.17632/7nhpbhvfnz.2.
Open access
