Erosion rate maps highlight spatio-temporal patterns of uplift and 1 quantify sediment export of the Northern Andes

Abstract


Introduction
Erosion plays a major role in controlling the large-scale topography and tectonics of mountain ranges, e.g., by redistributing mass (Beaumont et al., 1992;Wolf et al., 2022).In actively uplifting mountain ranges, erosion tends to balance rock uplift, making it possible to elucidate patterns of tectonic uplift that are otherwise hard to monitor on long (> 10 2 yr) time-scales (Kirby and Whipple, 2012).Sediments, the products of erosion, are finally transported to basins where they influence basin evolution and provide nutrients vital to ecosystem productivity (Hoorn et al., 2010).The key agents for mountain erosion and the dispersal of sediments to basins are rivers.This is a non-peer reviewed preprint submitted to EarthArxiv.
Where rivers are subjected to constant climatic and tectonic forcing, they tend to adjust their slopes so that erosion rates everywhere approximately balance rock uplift rates.In such steady-state streams, river steepness is a function of rock uplift and by extension, erosion, as well as erosional parameters such as bedrock erodibility, and climate (Howard, 1994).Spatial differences in uplift, erosional parameters, and climate lead to spatial variations in river steepness.Furthermore, temporal changes in, e.g., rock uplift lead to changes in river steepness that travel from the baselevel to the headwaters.The time-scale of landscape adjustment to new forcing conditions depends on the erosional parameters, climate, and river steepness, and typically is on the order of millions of years (Whipple et al., 2017).Therefore, temporal changes in boundary conditions can lead to disequilibrium landscapes, where different regions record different tectonoenvironmental conditions and erode at different rates.Hence, if spatial variations in river steepness, climate and erosional parameters can be constrained, one can use these to elucidate spatio-temporal patterns in erosion and thereby rock uplift.
Cosmogenic radionuclides (CRNs) measured in river sediments are commonly used to determine millennial, catchment-averaged erosion rates (Granger et al., 1996).Where rivers are well graded, CRNderived erosion rates can be used directly to infer patterns of rock uplift along mountain ranges (e.g., DiBiase et al., 2010).However, elucidating patterns of erosion over large spatial scales requires large CRN data sets, which are typically hard to generate.By analyzing only catchments with well graded river profiles, where uplift approximately equals erosion, one can separate out the effect of climate on topography (Adams et al., 2020).If the effects of climate and lithology on erosional parameters can be accounted for, the relationship between normalized river steepness and erosion rate can be used to infer erosion rates also for the disequilibrium parts of the landscape to study spatio-temporal patterns of rock uplift.
The Northern Andes of Colombia are a prime location to test the use of topography and erosion rate measurements to elucidate variations in rock uplift due to strong differences in uplift, climate, and bedrock lithology (Fig. 1).The two main mountain ranges of the Northern Andes, the Central Cordillera (CC) and the Eastern Cordillera (EC), are composed of dominantly crystalline basement rocks and clastic sedimentary rocks, respectively (Gomez and Montes, 2020) and experience significant spatial variability in precipitation (Urrea et al., 2019).Observed topographic disequilibrium has led to suggestions of pronounced spatial and temporal variations in rock uplift (Struth et al., 2017;Pérez-Consuegra et al., 2021b), but the relationship between subduction processes and inherited structures on the evolution of topography remain debated (Pérez-Consuegra et al., 2021a, 2021b).We propose that erosion rate maps could help to elucidate patterns of rock uplift through space and, if coupled with landscape evolution model predictions, through time, which can then be used to differentiate between tectonic and geodynamic models for the CC and EC.Here we present 25 new CRN-derived erosion rates from the CC, which we combine with published data from the EC to quantify erosion as a function of bedrock-related erosional parameters, climate, and topography.We developed an optimization algorithm to extrapolate our new and existing measurements and infer erosion rates from the steepness of river channels and precipitation rates at 30-m resolution for the Northern Andes.We then use these inferred erosion rate maps to interpret patterns of tectonic rock uplift in space and time and predict sediment fluxes to neighboring sedimentary basins.These analyses highlight the influence of subduction dynamics and inherited structures on the evolution of the CC and EC and quantify sediment export crucial to ecosystem productivity to the Andean foreland.

Study area
The Northern Andes are formed by subduction of the Nazca and Caribbean plates beneath the South American Plate (Taboada et al., 2000).South of 5-6° N, the Nazca plate subducts steeply below the South American plate at a rate of ~ 50 mm/a (Trenkamp et al., 2002).Flat slab subduction occurs north of 5-6°N since about 6-8 Ma (Wagner et al., 2017) (Fig. 1A).There is some debate, as to whether the slab imaged north of 5-6° belongs to a previously continuous Nazca slab that tore (Chiarabba et al., 2015;Wagner et al., 2017), or whether there is overlap between the Caribbean and Nazca slabs (Taboada et al., 2000;Kellogg et al., 2019;Sun et al., 2022).
On the surface, the Northern Andes manifest as three parallel, roughly north-south striking mountain ranges, the Western, Central, and Eastern Cordillera (WC, CC, EC).The CC and EC are separated by the intermontane Magdalena Valley (Fig. 1).The CC is primarily composed of Paleozoic (meta-)granitoids, Triassic meta-sediments and meta-intrusive rocks and extensive Jurassic and Cretaceous batholiths (Villagómez et al., 2011a).South of 5°N, Plio-Quaternary volcanic rocks from current arc volcanism are found in several locations near the crest of the CC.The EC is an inverted Mesozoic rift structure that has turned into a doubly-verging fold-and-thrust belt (Cooper et al., 1995).Bedrock of the EC is mostly Mesozoic quartzose sandstones and mudstones, with some minor occurrences of Cenozoic and Paleozoic clastic sedimentary rocks and crystalline basement (Cooper et al., 1995).
The topographic evolution of the CC and EC is debated.Sedimentary and thermochronological records suggest that the CC existed throughout the Cenozoic and potentially Cretaceous (Gómez et al., 2005;Villagómez et al., 2011b).Several Cenozoic periods of rapid exhumation have been inferred, however it remains debated if the northern CC rose to high elevations (> 2000 m) by about 25 Ma (Restrepo-Moreno et al., 2019) or only within the past 5 Ma due to recent slab flattening (Pérez-Consuegra et al., 2021b).For the EC, pollen and some thermochronological studies suggest strong surface uplift and exhumation since the Late Miocene (e.g., Hooghiemstra et al., 2006;Anderson et al., 2016), while other evidence suggests mountain building commenced by or before the Oligocene (Gómez et al., 2003;Horton et al., 2010).In addition to the timing and spatial patterns of surface uplift, the drivers of mountain building in the northern Andes remain debated, especially in the EC.Some studies argue that deformation is controlled by inherited Mesozoic rift structures (Mora et al., 2006(Mora et al., , 2013)), while others suggest that slab flattening drove changes in dynamic topography and Neogene uplift (Siravo et al., 2019).

Cosmogenic 10 Be measurements
We collected 25 river sediment samples along the CC to measure in-situ, catchment-average 10 Be erosion rates.In the northern CC, most streams exhibit major knickpoints.Our aim was to focus on steady-state channel reaches; therefore, we sampled basins that were either entirely below or above major knickpoints.Sand samples were sieved to extract the 250-500 μm grain size class and quartz was isolated by magnetic separation, froth floatation, and repeated etching with hydrochloric, hydroflourosilicic, and hydrofluoric acid.The purity of the cleaned quartz was checked with an Inductively-Coupled Plasma Optical Emission Spectrometer (ICP-OES).A 9 Be carrier was added, and samples dissolved in hydrofluoric acid, before extraction of Be with ion column chromatography.The 10 Be/ 9 Be ratios were measured at the CologneAMS (Dewald et al., 2013) relative to standards KN01-6-2 and KN01-5-3.The concentrations were corrected with a 10 Be/ 9 Be blank ratio of 6.0e-15 ± 1.1e-15.
For the erosion rate calculation, pixel-based production rates were calculated from a digital elevation model of the sampled catchments, with CRONUS calculator functions published by Balco et al. (2008), Balco (2017), and the Stone (2000) scaling scheme.Average catchment production rates were used for erosion rate calculation together , assuming a bedrock density of 2.65g/cm³ and attenuation length for every basin based on the average air pressure and rigidity cutoff using the 'rawattenuationlength' function from CRONUScalc (Marrero et al., 2016).Erosion rates were determined using the CRONUS bisection method, and, following DiBiase (2018), no topographic shielding correction was applied.In the EC, Struth et al. (2017) published 23 10 Be concentrations from river sediments, for which we recalculated erosion rates in the same manner as outlined above.

Stream power incision model
To predict erosion in the Northern Andes, we build on the stream power incision law for detachment-limited rivers (Howard, 1994), which combined with mass conservation, predicts the elevation change of a river bed as (1) where channel bed elevations are raised by uplift U and lowered by river incision E, which is a function of the upstream drainage area A, channel slope S, and the dimensional erodibility coefficient K. Exponents m and n are empirical constants related to incision process, basin hydrology, and channel geometry (Whipple et al., 2000).In quasi-equilibrium conditions (   = 0), when rock uplift and incision are balanced, this equation can be rearranged to show the commonly observed power-law scaling between local channel slope and drainage area (Flint, 1974) (2)  =   *  Channel steepness   is the channel slope normalized for the downstream increase in drainage area and concavity of the channel, Ѳ.Under quasi-equilibrium conditions, incision E can be substituted for uplift U creating a direct relationship between incision rate and channel steepness.To compare the channel steepness   among multiple streams, a reference concavity   needs to be determined, which results in the normalized channel steepness   (Wobus et al., 2006).

𝒌 𝒔𝒏 calculation and regression methods
We sampled channels in the CC with well graded topographic profiles, presumed to be near steady state, or within channels above knickpoint locations.To minimize the impact of varying erodibility K, we targeted catchments dominated by granites and gneisses, with a few basins containing some metamorphic mica schists.Sampling steady-state basins with homogeneous lithology enabled us to use the cosmogenicallyderived erosion rates to define the stream power parameters n and K by rearranging Eq. 4: (5)  =  *    .
We measured ksn using TopoToolbox (TT) (Schwanghart and Scherler, 2014) and the 30 m Copernicus Digital Elevation Model (DEM).Specifically, we calculated mean basin   for all sampled basins using χ -elevation regressions (TT 'chiplot' function).We employed a Bayesian optimization that linearizes stream profiles to constrain the reference concavity   (TT 'mnoptim') (Fig. 2).The concavity optimization was performed on the southern CC where stream-profiles are near equilibrium (Pérez-Consuegra et al., 2021b) to avoid a biased concavity estimate from transient profiles in the northern CC.This yielded a   = 0.5 (Fig. 2).To estimate   for the EC, the Altiplano surface was removed from the DEM, resulting in a bestfit concavity value of 0.45.
In Eq. 1, drainage area serves as a proxy for stream discharge.However, precipitation varies considerably across the study area, which could change the discharge-drainage area scaling and potentially bias   −  comparisons.To address this, we used a precipitation-corrected channel steepness  − equivalent to  − defined by Adams et al. (2020), where (6)  − = ( * )   *  with P referring to the upstream-averaged mean annual precipitation.For calculation, we used the 500 m resolution CHELSA mean annual precipitation grid (Karger et al., 2017) as weighting for the flow accumulation algorithm.To illustrate how our results would differ without this correction, we also show results based on   .

Inferred erosion rate maps
Provided that n and K are known and can be regarded as constant within each mountain range, Eq. 5 makes it possible to convert a  − map to an inferred erosion rate map.To do this, we smoothed river bed elevations using a constrained regularized smoothing (TT 'crs') with a tau value (elevation quantile) of 0.25 (to account for positive DEM errors in valley bottoms; (Schwanghart and Scherler, 2017) and smoothing value of 10.To calculate  − , we used a critical drainage area for stream initiation of 5 km², based on our observations of a systematic drop in  − at lower drainage areas with the smoothing parameters applied (Fig. S1). − was calculated using the best fit concavity values of 0.5 and 0.45 for the CC and EC, respectively.Subsequently,  − values were projected from streams onto the hillslopes by reversing the flow routing, ensuring that no smoothing occurs across drainage divides.Eq. 5 was used together with n and K to convert the  − to an erosion inferred erosion rate map.Projecting river incision values onto the hillslopes assumes that incision and hillslope erosion are coupled, which is an implicit assumption of all studies comparing   to cosmogenic nuclide erosion rates (e.g., DiBiase et al., 2010;Adams et al., 2020).To limit the variability in K, we estimate erosion rates in the CC only within areas underlain by crystalline basement rocks, including minor Quaternary volcanic exposures.Similarly, in the EC, we only estimate erosion rates within areas predominantly covered by siliciclastic sediments.Based on the range of mean  − of CRN-sampled catchments in the CC and EC, in both cases about 75% of the erosion rate map is interpolated between measured erosion rate data and 25% are extrapolated.
To derive the individual n and K parameters for the inferred erosion map, we develop an optimization algorithm that minimizes the misfit between measured and predicted erosion rates in the two cordilleras separately.Traditionally, power law regressions of   or  − and erosion rate data are used to determine the values of K and n.However, our approach has several advantages over using standard power law regressions: (1) Xiao et al. (2011) showed that selection of a proper regression method (non-linear model in linear space vs. linear regression in log-space) depends on the data distribution; (2) no mean basin  calculation is involved, making it viable to include non-equilibrium catchments, thereby increasing the number of available data points, (3) the consistency of erosion rate map and measured data is directly tested, and ( 4) predicting an erosion rate by projecting   onto the hillslopes and converting every pixel to an erosion rate ensures that ksn is weighted according to hillslope area.In other words, if the common assumption of equilibrium between hillslope erosion and channel incision rates is justified, then the erosion rate prediction has to consider hillslope area and their relation to local channel values, instead of using only channel-based   -regressions to define n and K.We tested a range of reasonable n (0.5 to 4.5) and K (1e-14 to 1e-6) values.For every parameter combination, the algorithm calculates an inferred erosion rate map and uses it to predict erosion rates for each sampling location, by taking the average of all upstream pixel values.Subsequently, a weighted and non-weighted misfit function are applied to define the best fit models, where   refers to the modelled erosion rate, dE to the reported erosion rate uncertainty, and n to the number of samples i.We apply two different misfit functions to investigate the effect of uneven erosion rate distribution on the optimization parameters.a normalized to the standards KN01-6-2 and KN01-5-3 with a nominal 10 Be/ 9 Be value of 5.35e-12 and 251 6.320e-12.Subtracted average blank ratio for corrections is 6. 0e-15 ± 1.1e-15. 252 Our 25 CRN samples were taken from different geomorphic subregions of the CC (Fig. 3).Perez-Consuegra 253 et al. (2021b) showed that the southern CC comprises high relief topography with steep river profiles (Fig.

Erosion rates
in the west and in the east gently slope down to towards the Magdalena River (Fig. 3B).Therefore, we separate the erosion rates from the CC into four geomorphic regions: southern CC, northern CC low relief, northern CC high relief, and the east sloping surface of the northern CC.
CRN-derived erosion rates in the CC are between 0.007 and 0.394 mm/a (Tab.1).Above the flab slab subduction, samples on the low-relief-high-elevation surfaces have low average erosion rates of 0.017 ± 0.004 mm/a (Fig. 2).Samples taken of the east-sloping part of the low-relief surface have higher average erosion rates of 0.045 ± 0.007 mm/a.Erosion rates in the high relief canyons of the northern CC are significantly higher and average 0.138 ± 0.022 mm/a including one sample from the WC.The highest measured erosion rates are found south of the slab tear in the deeply incised valleys of the southern CC with an average of 0.316 ± 0.039 mm/a (Fig. 3).A Mann-Whitney test shows that the differences in erosion rates from geomorphic domains are statistically significant (Tab.S1).

Erosion rate -𝒌 𝒔𝒏 relationship
Both our newly determined erosion rates from the CC and the ones published by Struth et al. (2017) from the EC, follow power law relationships between the measured erosion rate and   / − (Fig. 4).In the CC, the regression parameters from non-linear models yield n-values of 1.22-1.32and K of 7e-8 -1.1e-7 (yr -1 ).Values for n in the EC are higher (3.45-4.41)with a K of 2e-14 -1.8e-12 (m 0.1 yr -1 ).In contrast, the linear regression in log-space returns lower n and higher K values.Values of n and K depend on the regression method (non-linear model versus linear model in log-space), and the use of   versus  − .
We performed the same analysis using only quasi-equilibrium catchments (see supplement) and found similar results (Fig. S2), verifying that we sampled equilibrium parts of the landscape as intended.In the CC,   and  − regression parameters are similar, whereas in the EC,  − regressions yield lower n and higher K-values.This suggests that precipitation gradients can explain some non-linearity between channel steepness and erosion in the EC.

Inferred erosion rate maps
We used our new erosion rates from the CC and the ones published by Struth et al. (2017) for the EC to calculate inferred erosion rate maps for both cordilleras with our optimization approach.The distribution of measured erosion rates is skewed, with many low rates and fewer high rates.Therefore, the weighted misfit function is biased towards fitting the lower end of erosion rates where most data points exist.For the CC and EC data, this leads to an underestimation of high erosion rates (Fig. 5 A&D).In contrast, the nonweighted misfit function provides a better fit to the higher erosion rates, but underestimates lower erosion rates in the CC, and overestimates both intermediate and low erosion rates in the EC (Fig. 5 B & E).To minimize the influence of uneven data distribution on the regression, we attempted to bin the data into even spaced erosion rate bins before performing our optimization.However, erosion rates in the EC were too unevenly distributed to cluster into erosion rate bins with more than three samples per bin.Therefore, we develop an alternative misfit function.For the tested parameter combinations of n and K, we converted the weighted and non-weighted misfit matrices to percentiles.We added the two percentile matrices to determine the global minimum from both simulations.We refer to this method as 'semi-weighted'.Note, that despite some underprediction for high erosion rates in the EC case, the semi-weighted fit performs best in the sense of fitting low and high erosion rates reasonably well in both Cordilleras (Fig. 5 C & F).In the supplement we show the n-K-parameter space misfits for all fitting methods, which highlight the trade-offs between n and K (Fig. S3-S5).We take the parameters from our semi-weighted optimization as the best-fit model and use them to convert a  − map into an inferred erosion rate map (Fig. 6).To estimate uncertainties, we calculate maps of minimum and maximum erosion rate end members by overlaying inferred erosion rate maps from the three different misfit functions and taking the minimum and maximum pixel values across all methods (Fig. 6 B & C).
We used inferred erosion rate maps to estimate the total volume of eroded material for both Cordilleras (Tab 2).Our preferred model (semi-weighted misfit) indicates 5.4 km³/ka of erosion for the CC, with other models indicating a range from 4.1 to 8.7 km³/ka.In the EC, a larger spread can be observed with precipitation-corrected values ranging from 10.3 to 35.9 and a best-fit estimate of 20.5 km³/ka.Eroded volumes predicted by   optimizations are similar (Tab.S2).Based on our best-fit models the total sediment export from the EC is nearly four times larger than from the CC.

Limitations of erosion rate maps
Inferred erosion rate maps present a novel and versatile tool to study landscape evolution.Previous attempts used a coarse (5 km radius) moving window to map erosion values from the stream pixels to the hillslopes (Adams et al., 2020;Clementucci et al., 2022) or lacked the data to infer both n and K values (Clementucci et al., 2022).In contrast, our approach accounts not only for climatic-gradients (Adams et al., 2020) but also for variations in major rock type between cordilleras and maintains the original flow routing to not smooth values across drainage divides.This is especially important when investigating landscapes with drainage reorganization and allowed us to test the inferred erosion rate map for internal consistency by forward modelling catchment-average erosion rates.
A limitation of erosion rate maps is that the fitting parameters depend on the regression or optimization method used to derive them.For our data, linear regression models applied in log-space consistently returned lower power-law exponents compared to non-linear fits in linear space.Therefore, if n and K are derived from bivariate regressions, a careful error analysis should be conducted to choose the appropriate fitting method (Xiao et al., 2011).Typically, only quasi-equilibrium catchments are used to derive parameter predictions from regressions (Adams et al., 2020).Our optimization method allows us to include more data points, yielding more robust results.Moreover, the optimization method ensures that  − values are weighted proportionally to the size of the neighboring hillslopes.Due to an uneven distribution of erosion rates, we found different best-fit parameters using weighted versus non-weighted optimization misfit functions.However, it may be possible to bin more uniformly distributed erosion rate data sets in other locations prior to optimization to generate erosion rate maps following our procedure.In cases where the data distribution does not allow sufficient binning, we show that our semi-weighted approach yields a satisfactory regression between high and low erosion rate and normalized channel steepness measurements.
In any case, observed erosion rates should be compared to modelled ones to elucidate potential biases.
Another assumption of our analysis is that the pattern of rainfall has been stable during the integration time of the catchment-averaged erosion rates.Most of our erosion rates integrate until the mid-Holocene, with the slowest rates integrating into the Pleistocene.Paleo-precipitation models for the mid-Holocene and Last Glacial Maximum (Fick and Hijmans, 2017) suggest no major shifts in the patterns of precipitation across the study area, supporting this assumption.
The inferred erosion rate map assumes that channels and hillslopes are coupled, and local channel incision is in balance with hillslope erosion.This is a common assumption (e.g., Adams et al., 2020) and studies with high-resolution topographic data have shown a strong coupling between channels and hillslopes across large gradients of uplift rate (Hurst et al., 2019), despite potential time lags in areas of recent uplift rate change (Clubb et al., 2020).Hence, we assume that most hillslopes in the CC and EC are tightly coupled to the local channel gradient, allowing us to extrapolate incision rates onto the hillslopes.

Tectonic implications of the erosion rate map
In active mountain ranges, erosion tends to balance rock uplift (e.g., Brandon et al., 1998), and therefore our inferred erosion rate map can be used to elucidate tectonic signals.In equilibrium parts of the landscape (U = E), uplift rates should directly equal our inferred erosion rates.In non-equilibrium regions, incision dynamics make it possible to infer spatio-temporal patterns in uplift from our erosion rate maps.In particular, landscape evolution models indicate that areas close to base-level tend to be equilibrated with recent tectonic conditions, whereas upstream areas typically record past tectono-environmental conditions.
This pattern is consistent with our results in the northern CC and geochronologic evidence of a recent acceleration in uplift rates; our inferred erosion rate map indicates faster erosion of the downstream flanks of the Antioqueno Plateau in the northern CC compared to slow erosion of the relic upland low-relief plateau surface (Fig. 6).Rapid erosion of the plateau flanks has been previously inferred from 6-7 Ma Apatite (U-Th-Sm)/He ages (AHe) (Pérez-Consuegra et al., 2022) in the Cauca Canyon (Fig. 6).Our analyses support the interpretation that upstream plateau areas in the northern CC formed during a past period of slower tectonic uplift and that a recent acceleration in uplift led to the incision of deep canyons.
Perez An area with lower inferred erosion rates near the crest of the CC has recently been linked to glacial planation based on the close correlation between the extent of this area and the altitude of moraines (Fig. 6) (Pérez-Consuegra et al., 2021b).This highlights that the interpretation of inferred erosion rate maps should be limited to regions where fluvial erosion is the dominant erosion process on the time-scale of landscape adjustment.
In the EC, erosion rate variability is high and three important patterns can be observed: (1) the flanks of the EC erode rapidly and rim a slowly eroding central high plateau, (2) the eastern flank erodes more rapidly than the western flank, and (3) significant along-flank variation exists along the eastern EC, where erosion rates are high in the south and decrease towards the central part of the flank before increasing north towards the Cocuy range.A space-for-time-substitution suggests that uplift rates in the EC were slow in the past and accelerated more recently, leading to transient landscape adjustment similar to that observed in the northern CC (Struth et al., 2015(Struth et al., , 2017)).Several thermochronologic and geologic studies also suggest an increase in exhumation rates and mountain building since the Late Miocene (Mora et al., 2008(Mora et al., , 2013;;Siravo et al., 2019), which agrees with findings from pollen studies near Bogota indicating a rise from lowland elevations to > 2 km within the same time period (Hooghiemstra et al., 2006).We note though that this pollen-based paleoaltimetry is debated (Molnar and Perez-Angel, 2021).Furthermore, a faster eroding and therefore faster uplifting eastern flank is consistent with the earthquake distribution in figure 1A that shows a higher density of shallow earthquakes recorded along the eastern flank.Moreover, a balanced cross section from the southern EC (Mora et al., 2013) shows substantially deeper exhumation of rocks in on the eastern flank, where inferred erosion rates are highest (Fig. 7F).This suggests elevated erosion rates on the southeastern EC flank have persisted long enough to create asymmetrical unroofing of ~4 kilometers compared to the western flank (Fig. 7F).
In the EC, two regions stand out as erosional hotspots, the Cocuy Range in the northern EC and the southeastern flank of the EC.Both loci of erosion coincide with locations of rapid exhumation with AFT ages < 10 Ma and AHe ages < 5 Ma (Mora et al., 2008(Mora et al., , 2015;;Siravo et al., 2018;Pérez-Consuegra et al., 2021).The central portion of the eastern EC flank has lower inferred erosion rates and coincides with an area where the orogen widens, and several parallel thrust systems are active simultaneously (Jimenez et al., 2013) (Fig. 7C).This suggests that slower erosion rates in this region may be related to more distributed exhumation, occurring primarily along new thrust faults in the former foreland basin.In this region, balanced cross-sections also indicate more distributed unroofing (Fig. 7C).The activity of several thrust sheets leads to a transient slope reduction of the flank until the main deformation front has migrated basinward.The consistent spatial patterns between our erosion rate maps and tectonic data suggests that erosion rate maps can be used to identify transient adjustments of thrust tectonics during orogen growth.
In the southern part of the Middle Magdalena Valley, inferred erosion rates in the CC are higher than in the corresponding western flank of the EC (Fig. 7).Independent support for this observation can be found in the current geographic distribution of Plio-Quaternary alluvial fans.The Magdalena River is currently flowing on the eastern side of its intermontane valley.This is likely a consequence of the higher sediment flux from the CC compared to the EC, which produces larger alluvial fans that deflect the river eastwards (Fig. 7D).showing high inferred erosion rates on the eastern EC flank, independently supported by the balanced cross-section in (F).

Variations in climate and erosion parameters
Our optimization results show a large difference in exponent n, with n > 3in the EC, and < 2 in the CC (Tab.2).This means that in the CC, subtle differences in  − would suggest only small differences in erosion, whereas subtle differences in  − in the EC correspond to more substantial differences in erosion rate due to the larger exponent.Power law fits with  − also return lower n-values compared to models with   (Fig. 4).
A comparison of erosion rate maps based on  − and   suggests that only maps based on  − predict erosional hotspots in the southeastern EC and the Cocuy (Fig. 8A).Independent geochronologic data also support the existence of these erosional hotspots.This highlights the importance of including rainfall into   calculations (Adams et al. 2020) and underscores the significant influence of climate on relationships between erosion and topography in the Northern Andes.
By accounting for spatial variations in erosion parameters and climate gradients simultaneously, our analyses reveal important local variations in tectonic forcing that would otherwise be obscured in traditional river profile analyses.For instance, variations in   along the eastern flank of the EC are minor (Fig. S5A), however our erosion rate maps and independent thermochronology and neotectonic data suggest substantial variations in uplift rate along strike (Mora et al., 2010).This implies that the combined effects of large exponent n and precipitation patterns in the EC lead to the situation where subtle differences in   along the eastern EC translate into substantial gradients in inferred erosion rates and rock uplift, highlighting the importance of moving beyond   or  − analyses.We used our inferred erosion rate maps to predict sediment fluxes exported from the Northern Andes to neighboring sedimentary basins (Fig. 8B).The respective volumes per unit time were calculated by summing the pixel values of erosion rates on either side of the main drainage divides.This approach assumes that long-term sedimentary sinks in the CC and EC are negligible.Even though we have shown that sediment export from the southern CC is greater than from the western EC flank in the southern Middle Magdalena Valley (Fig. 7D), the total export from the EC into this basin is nearly three times higher.The largest volume of sediments is exported to the Llanos foreland basin with 9.4 km³/ka, whereas only about 0.6 km³/ka are being supplied to the intermontane Cauca Valley.High sediment export to the Llanos basin may have contributed to the exceptionally high biodiversity in the Northern Andean foreland by providing a large amount of nutrients (Hoorn et al., 2010).

Predicting sediment flux to foreland basins and coeval tectonic processes
This is a non-peer reviewed preprint submitted to EarthArxiv.
These predicted sedimentary fluxes agree well with the relative stratigraphic thicknesses deposited since the Neogene in the adjacent foreland basins.The highest thicknesses and volumes of preserved Neogene to recent sedimentary units (ca. 4 km, see Reyes-Harker et al., 2015) have been documented in the Llanos foreland basin adjacent to the eastern flank of the EC (Hermeston and Nemčok, 2013;Silva et al., 2013;Reyes-Harker et al., 2015).The second highest volumes of preserved Neogene to recent sedimentary units have been documented in the Middle Magdalena intramountain Basin with a thickness of up to 2 km (Moreno et al., 2013;Tesón et al., 2013;Reyes-Harker et al., 2015).In contrast, the preservation of Neogene to recent sedimentary units is minor in the Cauca Valley west of the Central Cordillera (Suter et al., 2008) --consistent with our predictions.
In the CC, sediment fluxes are almost identical between the eastern and western flank, suggesting symmetrical uplift of the orogen.In contrast, in the EC, asymmetrical erosion and sediment export suggest higher rates of tectonic uplift along the eastern flank of the EC.This interpretation is consistent with the off-center location of the main drainage divide separating the Upper Magdalena Valley and Llanos Basin (Fig. 8).In mountain ranges experiencing asymmetric uplift, the main drainage divide should move towards the mountain side with higher uplift rates (He et al., 2021).Higher uplift rates along the eastern flank of the EC may have congruently shifted the main drainage divide to the east.
Several structural studies have concluded that deformation along the western flank of the EC commenced in the Eocene/Oligocene (Gómez et al., 2003;Horton et al., 2010), which suggests that the locus of deformation and uplift has moved eastward since the Late Miocene (Mora et al., 2013;Siravo et al., 2018).
It remains unclear, however, what drove eastward migration of deformation and Late Miocene to recent surface uplift.It has been hypothesized that the location of widening north of 4.5°N may be related to inherited structures.It has also been postulated that faster exhumation rates in the eastern EC flank are linked to a deeper rifting in this region (Mora et al., 2015;Pérez-Consuegra et al., 2021a).However, this does not explain the temporal shift in uplift from west to east revealed by our erosion maps and other geochronologic data.Based on the low crustal thickness of the northern EC and the symmetric widening of the EC north of ~4.5°N together with the Late Miocene timing of increased mountain building, we speculate that Late Miocene slab flattening not only caused dynamic uplift (Siravo et al., 2019) but also shifted deformation eastward to the eastern flank of the EC, consistent with the predictions of geodynamic models (Martinod et al., 2020).This would suggest that changes in subduction geometry over the Cenozoic were the main driver of topographic evolution of the Northern Andes with a potential overprint by inherited tectonic structures.

Summary and Conclusions
We determined CRN erosion rates and used them together with published data and interpolation methods to generate an inferred erosion rate map of the Northern Andes.Our main findings are: (1) Subduction geometry exerts first-order control on spatial and temporal patterns of erosion in the northern Andes.CRN-derived erosion rates are highest in the southern CC above the normal slab subduction; erosion rates on the plateau surfaces in the northern CC were slower, with faster rates below major knickpoints.The pattern most likely reflects an acceleration of uplift rates in the northern CC in response to Late Miocene slab flattening.
(2) Topographic signatures of lanscape equilibirum and transience provide surface evidence of changes in subduction geometry along the northern Andean margin.The southern CC above the normal slab segment exhibits a steady-state topography which suggests that there were no major changes in subduction geometry during the Neogene, whereas the EC shows a pronounced topographic disequilibrium.
(3) Fast erosion of the EC flanks and low erosion of the interior Altiplano suggest strong topographic growth on the time-scale of landscape adjustment.
(4) Faster erosion rates in the eastern EC compared to the western flank show asymmetric tectonic uplift, in agreement with thermochronometric and geological data.
(5) Along strike variations in erosion rate on the eastern flank of the EC are most likely linked to different stages of wedge growth and fault stepping.
(6) Spatial differences in climate and erosional parameters highlight the importance of moving from   analysis to inferred erosion rate maps, when trying to use topography to study erosion and tectonics.
(7) Sediment flux from the EC is nearly four times higher than from the CC.The Llanos foreland basins receives the highest sediment flux of all foreland basins in Colombia, providing a large amount of nutrients; potentially related to the northern Andean foreland biodiversity hotspot.

Figure1:
Figure1: Overview of study area and CRN sampling locations.(A) Regional subduction zones and earthquake hypocenter depth (> Mw 4) (U.S. Geological Survey, 2020).Inset overview map of study area with slab depth contour lines from Wagner et al. (2017).Blue shading indicates region of flat slab subduction.(B) Catchment-average erosion rates in the CC and WC from this study (black outlines) and

Figure 2 :
Figure 2: Concavity optimization for the CC (A) and EC (B) showing the objective function from the Bayesian optimization with error bars.

Figure 3 :
Figure 3: (A) Catchment average erosion rates of the CC with approximate location of the slab tear at depth (for location see Fig. 1D).The catchment outlines are colored by geomorphic region with the color code shown in (E).(B) and (C) swath profiles across the northern and southern CC, respectively.Swath locations are shown in (A).Colors indicate different geomorphic domains.Note the low standard deviation of elevation in the northern CC highlighting the low-relief surfaces.(D) Erosion rate histogram grouped by geomorphic region.(E) Mean erosion rates for the geomorphic regions, computed by drawing normally

Figure 4 :
Figure 4: Erosion rate versus   and  − for the CC (left) and EC (right).Left panels show non-linear fit in linear space, right panels linear regression in log space.First row   , second row  − .Note that due to different concavities the units for   and K vary between the CC and EC.

Figure 5 :
Figure 5: Predicted versus measured erosion rates for the CC and EC for  − and different misfit functions.(A & D) Weighted misfit, (B & E) non-weighted misfit function, (C & F) semi-weighted misfit function.Data in the second row are plotted in log space to better visualize the fit of erosion rate data across scales.Colored arrows highlight biases in the optimization.Note that the semi-weighted approach partially mitigates the biases across the entire erosion rate range.

Fig. 6 :
Fig. 6: (A) Inferred erosion rate map for CC and EC with from the semi-weighted optimization.Minimum (B) and maximum (C) erosion rates computed from the range in values of the three estimates.Color range for (B) and (C) is the same as in (A).

Fig. 7 :
Fig. 7: Comparison between inferred erosion rate maps and balanced cross-sections.(A) Inferred erosion rate map highlighting the locations of swath profiles (B) and (E).(B) Swath profile across the northern CC and EC with elevation (black) and erosion rate (blue).Due to the high variability a spline was fit throughthe erosion rate data.(C) Balanced x-section(Mora et al., 2013) showing several thrust systems causing exhumation on the eastern flank of the EC.(D) Higher erosion rates in southern CC compared to western EC, create alluvial fans (white outlines) that shift the course of the Magdalena to the eastern side of theMagdalena Valley.Map location highlighted in (A).(E) Swath profile across the southern CC and EC

Fig. 8 :
Fig. 8: (A) Difference in erosion rate between best-fit semi-weigthed inferred erosion rate map based on  − and   .Positive values indicate higher erosion rates predicted by  − , and vice versa.Note positive values along the southeastern and northeastern EC flank highlight that only  − is predicts these independently documented erosional hotspots.(B) Volumes of eroded material supplied from the CC and EC to neighboring sedimentary basins from different sides of the main drainage divides.Sedimentary basins are shaded gray.Different colors highlight the areas that drain to different sedimentary basins.
Fig. S1: Median   and  − versus drainage area for the CC (left) and EC (right).Note the steep drop at drainage areas smaller ~5 km².

Figure S2 :
Figure S2: Erosion rates from quasi-equilibrium basins versus   and  − in the CC (left) and EC (right).Left panels show non-linear fit in linear space, right panels linear fit in log space.First row   , second row  − .A few of the low gradient plateau streams in the CC exhibit minor knickpoints due to small steps in the landscape.Therefore, we employed a r² threshold between the χ-elevation data and the linear χ-elevation fit to objectively define quasi equilibrium basins.We applied a r² criterion of 0.75, which is similar to previous studies (Hilley et al., 2019) and visually matches χ-elevation plot expectations.The r² criterion removes 6 out of 25 basins in CC and 7 out of 23 in the EC.Another 3 samples in the EC are nested catchment samples within close proximity of each other and were therefore removed, leaving 13 samples for the EC.

Figure S3 :
Figure S3: Misfit (unitless) for weighted misfit function and the tested n and K parameters combinations for   and  − .Best fit model indicated by the star.The misfit increases rapidly away from the best-fit solution and therefore we set the maximum of the color map to the misfit of the top 2% solutions.

Figure S4 :
Figure S4: Misfit (m²/a²) for non-weighted misfit function and the tested n and K parameters combinations for   and  − .Best fit model indicated by the star.The misfit increases rapidly away from the best-fit solution and therefore we set the maximum of the color map to the misfit of the top 2% solutions.

Figure S5 :
Figure S5: Misfit for semi-weighting misfit function.Best fit model indicated by the star.The color indicates the percentile of the misfit from 0 to 100%, with 0% being the lowest misfit of all parameter combinations.

Fig. S5 :
Fig. S5:   (A) and  − (B) map of the Northern Andes.The color bars for the CC and EC are adjusted differently with the upper limit at the 98% percentile, because the units differ between cordilleras, due to the different concavities.

Table 1 :
AMS data and erosion rates from the CC.250

Table 2 :
Best fit optimization parameters and eroded volumes based on  − .
-Consuegra et al., (2021b)ascribe this acceleration in surface uplift to dynamic uplift due to Late Miocene/Pliocene slab flattening.South of the slab tear, erosion rates in the CC are substantially higher and less variable compared to the northern CC (Figs.6, 7B,E).The lower erosion rate variability suggests that this part of the mountain range is close to a steady-state topography.These observations corroborate the hypothesis byPerez-Consuegra et al. (2021b)that slab flattening since ~ 6 Ma led to increased uplift in the northern CC, whereas the steady-state topography of the southern CC suggests that no major changes in subduction geometry occurred over the time-scale of landscape adjustment.

Table S1 :
p-value matrix for Mann-Whitney test.All geomorphic regions show statistically significant differences in erosion rate (2-sigma).Only samples from the eastern slope of the northern CC are not statistically different from other populations due to the small sample size (n=3).

Table 2 :
Best fit optimization parameters and eroded volumes based on   .