Plateau reduction by drainage divide migration in the Eastern Cordillera of Colombia defined by morphometry and 10Be terrestrial cosmogenic nuclides

Catchment‐wide erosion rates were defined using 10Be terrestrial cosmogenic nuclides for the Eastern Cordillera of the Colombian Andes to help determine the nature of drainage development and landscape evolution. The Eastern Cordillera, characterized by a smooth axial plateau bordered by steep flanks, has a mean erosion rate of 11 ± 1 mm/ka across the plateau and 70 ± 10 mm/ka on its flanks, with local high rates >400 mm/ka. The erosional contrast between the plateau and its flanks was produced by the increase in the orogen regional slope, derived from the progressive shortening and thickening of the Eastern Cordillera. The erosion rates together with digital topographic analysis show that the drainage network is dynamic and confirms the view that drainage divides in the Eastern Cordillera are migrating towards the interior of the mountain belt resulting in progressive drainage reorganization from longitudinal to transverse‐dominated rivers and areal reduction of the Sabana de Bogotá plateau. Copyright © 2016 John Wiley & Sons, Ltd.


Introduction
Rift inversion and/or crustal thickening result in dramatic topography and drainage network changes in mountain belts that are mainly controlled by major tectonic structures ( Van der Beek et al., 2002). The evolution and development of drainage systems in active orogens has attracted much attention recently, particularly with regard to transverse and longitudinal drainages, and drainage divides (Willett et al., 2001Pelletier, 2004;Bonnet, 2009;Castelltort et al., 2012;Perron et al., 2012;Goren et al., 2014;Viaplana et al., 2015). Clear examples of drainage rearrangement in relation to progressive mountain building have been described in the inverted rifts of the High Atlas of Morocco (Babault et al., 2012) and Eastern Cordillera of Colombia (Babault et al., 2013;Struth et al., 2015). These studies highlight that drainage divides are dynamic features that progressively migrate and result in river capture events. Providing evidence for differential erosion rates on either side of a drainage divide adds much credence to developing and substantiating drainage evolution models for mountain belts . As such, the side of a mountain belt with the greatest erosion will progressively migrate and capture drainages from the side with the least erosion. Theoretically, a dynamic drainage network will evolve towards a steady state, maintaining a steady drainage network and stationary drainage divides (Howard, 1965).
The Eastern Cordillera of Colombia is an example of an orogen with a dynamic drainage network (Babault et al., 2013;Struth et al., 2015). This N-S oriented orogenic belt has two topographic domains: (i) an axial zone with low relief associated with longitudinal rivers that have gentle gradients (the Sabana de Bogotá); and (ii) high-relief flanks with steeper transverse rivers. These domains are separated by what we refer to as the eastern and western main drainage divides. Recent studies have suggested that the increase of regional slopes due to progressive crustal thickening results in fluvial reorganization from longitudinal to transverse dominated drainages in the Eastern Cordillera (Babault et al., 2013;Struth et al., 2015). On the basis of a morphometric analysis, field observations and a summary of paleodrainage data, these studies conclude that drainage reorganization takes place by progressive drainage divide migration toward the axial zone of the orogen by a step-by-step series of river captures.
In this study, we build on the work of Struth et al. (2015) by using 10 Be terrestrial cosmogenic nuclides (TCNs) in fluvial sands to determine catchment-wide erosion rates for the Eastern Cordillera to investigate the contrasting erosion dynamics between the Sabana de Bogotá axial plateau and the flanks of the Eastern Cordillera. TCN analysis is combined with a comparative analysis of steepness index, specific stream power and the integral of the drainage area (χ values). We also

Geological, morphological and climatological setting
The Eastern Cordillera of Colombia is an inverted continental rift in the northern Andes. The cordillera is composed of Precambrian-Paleozoic basement and a succession of Mesozoic and early Cenozoic sedimentary rocks within a doubly-verging thrust system (Julivert, 1970;Colletta et al., 1990;Cooper et al., 1995; Figure 1).
The Eastern Cordillera is flanked on both sides by lowlands composed of Cenozoic deposits, which infill the Magdalena Valley and the Llanos foreland basins that are at elevations of 200-300 m above mean sea level (asl; Figure 1(A)). The flanks of the Eastern Cordillera are dominated by alternating Cretaceous sandstone and shale formations. Isolated Precambrian-Paleozoic basement massifs with low-and medium-grade metamorphic rocks (mainly phyllite and schist) are present in the Quetame and Floresta massifs (Segovia, 1965;Ulloa and Rodríguez, 1979;Ulloa and Rodríguez, 1982;Parra et al., 2009a; Figure 1 Teixell et al., 2015). Total shortening along this section is 82 km, consistent with values previously reported by Tesón et al. (2013) for the segment of the Eastern Cordillera covered by the map. [Colour figure can be viewed at wileyonlinelibrary.com] large-displacement (Mora et al., 2008 and references therein; Figure 1(B)). The local relief in the Sabana de Bogotá plateau is small because up to 600 m of unconformable Pliocene-Quaternary fluviolacustrine deposits partially fill the synclinal depressions (Tilatá and Sabana formations; Julivert, 1963;Andriessen et al., 1993;Torres et al., 2005). A moderate amount of orogenic shortening, 25 to 30%, across the Eastern Cordillera has been calculated along serial transects by Tesón et al. (2013) and Teixell et al. (2015).
There is a strong precipitation gradient across the western divide, whereas no precipitation gradient exists across the eastern divide ( Figure 2). Precipitation is greatest (> 3000 mm year À1 ) near Charalá, in the eastern and western foothills and in the Quetame Massif.
Summary of the deformation and paleodrainage evolution of the Eastern Cordillera Paleogene data show a western and southwestern source for the Central Cordillera foreland basin sediments, which currently include the early successions of deposits in the Magdalena Valley and Llanos basins, and the intervening Eastern Cordillera. Sediments moved NNE in the Central Cordillera foreland basin following the regional slope and the contemporaneous structures within the basin (Cooper et al., 1995;Bayona et al., 2008;Silva et al., 2013). The drainage was longitudinal, parallel to the current Eastern Cordillera, following growing folds that thermochronological data show started to form during the Eocene (Parra et al., 2009b;Mora et al., 2010;Moreno et al., 2011;Saylor et al., 2011;Caballero et al., 2013aCaballero et al., , 2013bSilva et al., 2013).
The Magdalena and the Llanos basins became disconnected by the uplift of the Eastern Cordillera as deformation propagated eastward during the Late Oligocene to Early Miocene (Parra et al., 2009b;Mora et al., 2010). This is supported by provenance data indicating that the Eastern Cordillera developed into an effective topographic barrier that separated the Central Cordillera from the Llanos basin before the Mid to Late Miocene . The axial Eastern Cordillera was likely a closed basin south of 6°N during the late Oligocene-Miocene .
Base level started to rise in the Middle Magdalena Valley basin during the Late Oligocene to Mid Miocene, forcing its rivers to flow to the east, across the present location of the Eastern Cordillera, into the Llanos basin (Gómez et al., 2005). The drainage in the Middle Magdalena Valley basin returned to flow to the north due to the continued uplift of the cordillera during the Mid-Late and Late Miocene. Since the Late Miocene, the paleoflow pattern in the eastern foreland of the Eastern Cordillera is characterized by a transverse drainage reflecting an eastward direction  and providing sediments for the Llanos and Middle Magdalena Valley basins. Most external thrust sheets of the Cordillera are strongly deformed in recent times, and there are very young apatite fission track ages for the Quetame Massif along its eastern flank (0-3 Ma; Figure 3). This attests to strong exhumation across the Eastern Cordillera margins in the latest Neogene and Quaternary (Mora et al., 2008) . Struth et al. (2015) argue that the Eastern Cordillera is experiencing fluvial drainage network reorganization by drainage divide migration. Longitudinal drainages created in the early stages of drainage developed paralleling the structural grain of the growing tectonic orogen. The potential energy and power of transverse rivers was enhanced as the regional slope progressively increased during crustal thickening. This resulted in headward erosion causing drainage divide migration towards the Sabana plateau and capturing longitudinal streams: ultimately leading to drainage network reorganization.

Analysis of the digital drainage network
We used digital elevation model (DEM) SRTM90v4 that has a horizontal resolution of 90 m (Jarvis et al., 2008) to analyze the drainage network for the Eastern Cordillera. The DEM was corrected in narrow areas that had low resolution using elevations from Instituto Geográfico Agustín Codazzi (IGAC) 1:100 000 topographic maps. Various geomorphic parameters provide information about erosion in channel networks and therefore information about the dynamism of the drainage. The most important parameters that we examine include stream power, steepness index and χ. We calculated these parameters to characterize the contrast between the topographic domains of the Sabana de Bogotá plateau and its flanks, and to identify which one best correlates with the TCN-derived erosion rates. We propose that a correlation between geomorphic parameters and TCN-derived erosion rates may allow a first estimation of erosion rates in areas devoid of TCN data.Stream power relates to the energy rate per unit distance along a river (Bagnold, 1966) and reflects the incision power of a river into bedrock under detachment-limited conditions (Howard and Kerby, 1983;Howard et al., 1994;Whipple and Tucker, 1999;Kirby and Whipple, 2001). River incision is calculated according to the stream power model expressed as: where U represents the bedrock uplift at x, K is the erosional efficiency, A is local drainage area and S is the channel slope. The positive exponents m and n describe the relative dependency of stream erosion rates on A and S.
Analysis of slope-area has frequently been used to reveal erosion trends in channel networks defined as steepness index or Ksn (Kirby and Whipple, 2001;Kirby, 2003;Snyder et al., 2003;Wobus et al., 2006;Dibiase et al., 2010). However, local slope calculations undertaken for regions that have low DEM resolution such as the Eastern Cordillera of Colombia provide scattered and noisy results. For this reason, we calculate the channel slopes using the χ gradient referred to as Mx Mudd et al., 2014). To calculate Mx, we use the river profile elevation instead of the slope as the dependent variable, against χ as the independent variable. This approach produces more reliable results. We plot values of Mx, i.e. the slope in χ-elevation space (Mudd et al., 2014), a parameter related to the ratio between erosion rate and erodibility , as a proxy to identify the distribution and magnitude of erosion. We defined the best concavity of the river profiles based on AICc-collinearity tests and χ-plots following the method of Mudd et al. (2014) to calculate Mx values. The AICc is a statistical method that selects a model that balances goodness of the fit against model complexity (Akaike Information Criterion, AICc;Akaike, 1974; Hurvich and Tsai, 1989; Burnham and Anderson, 2002). We extract the AICccollinearity test and the χ plots for each basin based on iteration through a range of concavity values for the main channel and tributaries to define the best-fit concavity. The concavity with the minimum AICc value will be the best fitting one using this method. The mean calculated concavity for all the basins is 0.45 (see SD1). We compared the Mx values between different catchments by fixing the obtained concavity. Stream power theory predicts that river profiles will have a linear χ-plot and Mx proportional to the erosion rates, assuming, first, a steadystate condition, where rock uplift is balanced by erosion, and second, that erosion, erodibility and uplift are constant through time and space (Royden and Taylor Perron, 2013;Mudd et al., 2014). In reality, U and K can be variable in space and time, which are dependent upon the tectonic and climatic history, and rock type. We calculated the concavity for all the basins following the method of Mudd et al. (2014) by extracting the steepness index, knickpoints and localize all the known active structures as well as highly erodible lithologies to solve this problem. In the case of stepped channel profiles associated with spatial or temporal changes in uplift rates, we used the colinearity test to identify the best concavity for each catchment (Mudd et al., 2014). This technique allows the magnitude and distribution of the erosion rates to be identified. We extract the χ-map following the method of Willett et al. (2014), and with a same base level determined by the elevation of the most external tectonic structures (300 m asl) and with a critical drainage area of 1 km 2 to extract information about the dynamics of the catchments. The main focus of this analysis rests on mapping differences in χ-coordinate values across drainage divides. Similar χ values on both sides of a drainage divide would suggest that the region is in equilibrium, while large differences in χ values across the drainage divide imply that river networks are in disequilibrium where divide migration or river capture is likely to occur . Drainage divides generally migrate towards the higher χ values to achieve equilibrium, and hence, catchments with high values are prone to capture and may eventually disappear . We extract the χ-map for the entire central segment of the Eastern Cordillera as a proxy for the dynamics of the drainage divides using a modified version of χ that makes a correction factor for the precipitation (Yang et al., 2015) such that the χ value is defined by the following equation: where P 0 and A 0 are arbitrary scaling factors for the precipitation rate and drainage area, respectively, P is precipitation rate, A is the upstream drainage area and m and n are empirical and non-integer constants.
We also calculated an averaged specific stream power (SSP, Equation (3), Knighton, 1999) for each catchment following the method described in Godard et al. (2012) such that: where d is density, g is acceleration due to gravity, Ksn is the steepness index slope, Q is discharge and W is the channel width.
Terrestrial cosmogenic nuclides analysis Sediment samples were collected for 10 Be TCN analysis to investigate erosional contrasts between the Sabana plateau and its flanks, including the main streams and major tributaries in the Guayuriba and Turmequé (representing the eastern flank of the Cordillera), Bogotá (Sabana de Bogotá plateau) and upper Suárez catchments. We followed the sampling methods of Carretier et al. (2015) who collected fluvial sand from the subcatchments to minimize issues associated with local lithologic variations. Rivers in the eastern flank and plateau have their headwaters on the same drainage divide, which allows us to make direct comparison between them. We also analyzed four 10 Be TCN samples in the lower Suárez and Chicamocha catchments ( Figure 3, Table I) to examine the effect of a different vegetation cover, bedrock erodibility and climate. These four samples are not included in the comparative analysis between the plateau domain and the eastern flank.
Larger sampled catchment areas have higher likelihoods of including mixed fluvial signatures. In such cases, the erosion rate obtained may not be representative of all the reaches of that river, which may have experienced a diverse reorganization history.
Catchments including a relict flat area in their upper parts (Chicamocha, Suárez, Bogotá and Turmequé) reflect two different erosion regimes separated by a knickpoint (white stars in Figure 6). As such the upper part the catchment, above the mean knickpoints, should have lower erosion rates (related to lower relief, gentler slopes and lower Mx values) compared with those downstream. The Turmequé catchment illustrates this well, with an ancient plateau domain above the main knickpoint where the erosion rates are low (LUCN-45, 26) and higher erosion rates in its lower reaches 27,28). We only sampled the upper part of the Bogotá catchment 19,30,31,33) above the Tequendama Falls (see Figure 4 for location) to define the amount of erosion in it lowest erosion rate domain. We sampled its gently sloping upper reaches 50), and the steeper lower reaches below the knickpoint (LUCN-53, 55) in the Suárez catchment. Only one sample (LUCN-53) was collected in the Chicamocha valley catchment, which was located downstream of the main knickpoint.
We collected~1 kg of sand for each sample. Catchment areas for each sample were large (>100 km 2 ) as recommended by Niemi et al. (2005) and Yanites et al. (2009) to provide a good representation of the 10 Be TCN inventory and the erosion rates. We use the erosion values obtained in the upstream reaches of the rivers that flank the Eastern Cordillera to compare the erosion rates between the rivers of the Sabana and the flanks of the cordillera. These flank rivers and those located in the Sabana, flow through alternating sandstone and shale formations in proportions that are similar in each different catchment. This similarity reduces the possibility of lithological biasing in sampling between our chosen catchments. The erosion rates determined by 10 Be reflect erosion on timescales of tens of thousands of years (Von Blanckenburg, 2005). In addition we sampled rivers across the western flank of the Eastern Cordillera, but we could not determine 10 Be concentrations in this region because we were unable to extract quartz from samples which were dominated by argillaceous sediment.
We also collected five quartzite samples for 10 Be TCN surface exposure dating from a well-preserved river terrace in the Guayuriba Basin (Figure 3(C), Table II). Such landforms are rare in our study area. Our sampled river terrace, which lies at~406 m above the current stream, is not deformed, and it shows no evidence of erosion or landsliding. We collected 500 g rock from the upper surfaces of fresh quartzite boulders inset into the river terrace. Topographic shielding was calculated by measuring the inclination to the skyline for every cardinal and intercardinal direction following the approach by Balco et al. (2008).
Quartz isolation, purification, dissolution and preparation of BeO were undertaken in the geochronology laboratories at the University of Cincinnati following the methods of Kohl and Nishiizumi (1992) and described in detail in Dortch et al. (2009). All river sediment and river terrace boulder samples were sieved and crushed to 250-500 μm for the analysis. Samples were cleaned using HNO 3 , HCl and HF and passed through a Frantz magnetic separator. Density separation was undertaken using lithium heteropolytungstate followed by an additional HF leach (Brown et al., 1991;Kohl and Nishiizumi, 1992;Cerling and Craig, 1994). For each sample, 15-30 g of clean quartz grains was dissolved in HF and HNO 3 with 350 mg of 9 Be carrier. Beryllium was separated using anion and cation exchange columns. Beryllium hydroxide was obtained after fuming with HClO 4 acid and passing through anion and cation exchange columns (Bourlès, 1988;Brown et al., 1992). The Be(OH) 2 was heated in an oven at 750°C to form BeO and then loaded in steel target mixed with of Niobium (Nb) powder. 10 Be/ 9 Be ratios were measured by Accelerator Mass Spectrometry in the Purdue Rare Isotope Measurement (PRIME) Laboratory at Purdue University, Indiana, USA.
We took the standard approach and assumed that the amount of 9 Be in the prepared quartz was negligible. Sometimes this may not be the case, and quartz may contain some 9 Be. Such cases are rare, and almost invariably occur when quartz is derived from beryl-bearing granites or pegmatites, so that traces of Be and/or small amounts of beryl exist in the supposedly pure quartz sample, as documented in some areas of the Himalaya (Portenga et al., 2015). However, pegmatites and beryl-bearing granites are not present in our study area. Any significant native 9 Be is therefore highly unlikely in our sampled quartz. If native 9 Be is present then denudation rates calculated from 10 Be concentrations will be overestimates, and the calculated erosion rates should be considered as apparent; however these rates can still be used for comparisons of variation across a region of similar lithologies (Corbett et al., 2013;Portenga et al., 2015).
To model catchment-wide erosion rates from 10 Be results we assumed that: (i) the sediment volume is proportional to the erosion rate, i.e. the catchment is close to steady-state; (ii) all the sediment collected at the sample locality was well mixed; (iii) the contribution of quartz was homogeneous in the catchment; (iv) that there is an isotopic equilibrium within the catchment, i.e. TCNs production in the catchment equals the transport of TCNs out of the catchment; and (v) the erosional timescale was significantly larger than the sediment transfer through the catchment (Lal and Arnold, 1985;Granger et al., 1996;von Blanckenburg, 2005).
The corrected catchment-averaged production rates for 10 Be were calculated for each catchment using the SRTM90 digital elevation model following the method of Dortch et al. (2011) using MATLAB v. 2008 and the scaling factors provided in Lal (1991) and Stone (2000). With this method, the production rate for each pixel is calculated accounting for shielding to make an average for all the pixels in the catchment to obtain a spatially average production rate for the entire catchment (Lal, 1991). We use a 10 Be half-life of 1.36 AE 0.07 Ma (Nishiizumi et al., 2007), the scaling model of Lal (1991) and Stone (2000) and a sea-level high-latitude production rate of 4.49 AE 0.39 10 Be atoms/g SiO 2 /a, for catchment erosion rates and river terrace dating. 10 Be ages for the river terrace samples were calculated using the CRONUS calculator (http://hess.ess.washington.edu/math/; Balco et al., 2008). These ages will be minimum values and hence the resulting incision rates are maximum values. We are aware that the production rates and scaling models for 10 Be are being refined Marrero et al., 2016). Borchers et al. (2016) recently revised the production rate for sea level at high latitudes to 4.01 10 Be atom/g SiO 2 /a for the Lal (1991) / Stone (2000) scaling model. However, we prefer to use the established scaling model and production rates in Balco et al. (2008) until there is community-wide agreement on the appropriate production and scaling. The new production rate results in an~10% difference in the erosion rate from our preferred values.

Digital drainage analysis
The axial zone, which includes the Sabana de Bogotá and the southern part of the Sogamoso basin, is characterized by the lowest SSP values (Figure 4). In contrast, higher SSP values occur on the Eastern Cordillera flanks and in the northern areas of the Sogamoso basin, which argue for higher erosion rates than in the axial zone of the cordillera. We extend the Mx map for the Eastern Cordillera of Colombia of Struth et al. (2015) that was limited to the area between 4°N and 5°30′N to between 4°N and 7°N. The χ values in the axial region of Sabana de Bogotá of the Eastern Cordillera ( Figure 5) are higher than those for its flanks in the area of the drainage divide. A strong differentiation in χ values exists within the same domain: while the Sabana de Bogotá has low values near the central divide, the Sogamoso basin shows high values near the divide (Figure 5(B1)). In addition, χ values indicate disequilibrium between the Turmequé (with low values) and the Chicamocha (with higher values near the divide) basins (Figure 5(B2)).
Low SSP, Mx and precipitation, and high χ values in the headwaters of the rivers characterize the plateau, while high SSP, Mx and precipitation, and low χ values in the headwaters of the rivers characterize the flanks. In summary, the geomorphic indices clearly differentiate the plateau and flank domains within the Eastern Cordillera.

Landscape evolution rates
The large erosional contrast between the plateau and its flanks is evident from the variation in erosion rates derived from 10 Be analysis ( Figure 6). The 10 Be analysis shows that for similar drainage areas, the erosion rates are lower in the Sabana de Bogotá plateau region (< 20 mm/ka) than on its flanks (>400 mm/ka (Table III). In more detail, the upper half of the Guayuriba catchment is eroding at a rate of~75 mm/ka (samples LUCN-01, 03, 04, 08; Figure 6), whereas downstream the erosion rates are~100, 479 and 670 mm/ka (LUCN-07, LUCN-05 and LUCN-06). Two samples (LUCN-05 and LUCN-06) show that erosion rates are an order of magnitude higher than the rest of the flank samples and provide the highest erosion rates (see discussion below). In the Turmequé catchment, the erosion rates increase progressively downstream from 5.5 mm/ka in the upper part of the catchment in the longitudinal tract (LUCN-45) to 15.7, 48, 53 and 59 mm/ka (LUCN-26, 28, 27 and 25) in the transverse tract ( Figure 6).
The Sogamoso catchment drains the axial plateau of the Eastern Cordillera to the north, and is eroding at a rate of 167 mm/ka (LUCN-54; Figure 6). This catchment is divided into the Suárez catchment in the west and the Chicamocha catchment in the east, which are eroding at a rate of 64 and 228 mm/ka (LUCN-55 and 56), respectively. The upper part of the Suárez catchment is eroding at a rate of 8 and 10 mm/ka (LUCN-46 and 50). The erosion rate for the middle part of the Sogamoso catchment is 46 mm/ka (LUCN-53), and is similar to the erosion rate for the Suárez catchment.
TCN exposure ages for terrace boulders in the Guayuriba basin, assuming zero erosion and considered as minimum ages, range from 257 to 697 ka (Figure 7). However, the ages will be significantly older if we estimate an erosion rate using the oldest boulder age (LUCN-12) by applying the methods of Lal (1991). The erosion rate obtained from the oldest boulder is 1 mm/ka and if we assume that all the samples erode at this rate, the initial ages with zero erosion will increase by up 50% for ages up to~700 ka. However, since we sampled large boulders with no apparent weathering, we can assume that a correction factor for erosion is not really necessary for our ages.
The ranges of ages show significant scatter, but the three oldest samples (LUCN-11 to LUCN-13) clustered particular well given their antiquity, with a mean age of 656 AE 69 ka. We argue that the two youngest ages probably reflect exhumation of the boulders from the terrace and so we do not use those ages in our incision rate calculation. Given the mean elevation of the three oldest samples is 2129 AE 3 m asl and the elevation in the adjacent valley is~1723 m asl, we use a height difference of 406 AE 5 m to calculate an incision rate of~62 mm/ka.   This is comparable with the catchment-averaged erosion rates calculated from river sediment from the flank streams.

Discussion
Calculated erosion rates for the Eastern Cordillera using 10 Be and values of SSP, Mx and χ show clear erosional contrasts between the axial plateau and its flanks. Intra-domain and intra-basin differences are also apparent between regions.

Drainage network dynamics
Rivers across the flanks of the Eastern Cordillera need more erosion capacity than the longitudinal rivers in the Sabana plateau for river capture to occur (Struth et al., 2015). Higher erosion rates along the flanks of the Eastern Cordillera (>50 mm/ka) compared with the Sabana plateau (<20 mm/ka) suggests that retreat of the main divides and capture of the longitudinal plateau rivers by the transverse flank rivers is occurring. The Guayuriba and the Bogotá catchments provide examples of flank catchments with a drainage area of 200 km 2 yielding an erosion rate of 79 mm/ka (LUCN-01), whereas the plateau catchments with comparable drainage areas of 163 and 245 km 2 yield erosion rates of 18 (LUCN-30) and 4 mm/ka (LUCN-31), respectively (Figure 6). For a plateau catchment with an area of 240 km 2 , an erosion rate of 15 mm/ka (LUCN-15) was determined, which is lower than a similar size catchment (293 km 2 ) on the eastern flank that has an erosion rate of 77 mm/ka (LUCN-08). Comparisons between the Turmequé and Sogamoso catchments also illustrate this contrast; sampled catchments for Suárez, with areas of 313 and 1186 km 2 , provided rates of 7 and 10 mm/ka (LUCN-46 and 50). This is markedly different from the Turmequé flank samples from drainage areas of 460 and 1291 km 2 that have higher erosion rates of 48 and 59 mm/ka (LUCN-28 and 25), respectively. These results confirm that the main drainage divides in the Eastern Cordillera will migrate towards the axial plateau because rivers in the flanks have more energy (indicted by SSP values) and erosion capacity (reflected in the Mx values).
A special case is illustrated by comparing the Sabana catchment where sample LUCN-19 was collected (with a drainage area of 265 km 2 and an erosion rate of 11 mm/ka) with the Turmequé catchment flank where samples LUCN-45 (65 km 2 , 6 mm/ka) and LUCN-26 (346 km 2 , 15 mm/ka) were collected. This sampled area round these samples was interpreted by Struth et al. (2015) as a drainage capture zone on the basis of topography, featuring a depressed topographic profile and a mapped reentrant toward the west of the main divide, and of the occurrence of knickpoints upstream of a fluvial elbow (sharp change in the river channel direction) in map view. Similar erosion rates on both sides of the drainage divide in the Turmequé area argue in favor of drainage capture. In addition, the χ-values map shows that the upper part of the Turmequé catchment was part of the plateau area in the past and is now incorporated into the flank domain ( Figure 5(B1)).
As documented by Struth et al. (2015), the Eastern Cordillera of Colombia has clear geomorphic differences between its flanks and axial zone, which includes the Sabana de Bogotá and the southern low-relief part of the upper Sogamoso catchment. Comparison of the geomorphic indices and the erosion rate results with the drainage divide features argue for capture processes. This confirms the view of Struth et al. (2015) and Babault et al. (2013) who proposed drainage divide migration and a longitudinal to transverse drainage rearrangement in the Eastern Cordillera.
Correlation of erosion rates with Mx values (Figure 9) is stronger than with the SSP (Figure 8). The SSP is, at a first order, a function of precipitation, slope and drainage area. Calculation of the SSP from a series of raster measurements and the low-resolution data of the DEM or the precipitation provide equivocal results. Mx is solely a function of the elevation and distance of the river profile, providing a more realistic value to describe the steepness of the river reaches.

Calibration of geomorphic parameters
The links between erosion rates and SSP or Mx have been examined in previous studies (Kirby and Whipple, 2012;Perron and Royden, 2013;Safran et al., 2005;Bookhagen and Strecker, 2012). Our data reveal high erosion rates derived from samples located on the flanks of the Eastern Cordillera (which show high SSP) and in the northern parts of the Sogamoso basin, and low rates for the axial plateau (Sabana de Bogotá, upper part of Turmequé basin and southern part of the Sogamoso basin). Geomorphic indices cluster into four groups (A, B, and C in Figure 8(B)). Group A comprises Sogamoso catchment samples from the most incised northern part of that catchment 54,55,56). Group B is associated with flanks that are characterized by high relief, and group C includes the lowest relief area samples. Group C is subdivided into C′, including samples with flank characteristics but with lower relief 27,28 and 01), and C″, with plateau characteristics 26,31,33,45,46,50).
Erosion rates and the catchment-wide average Mx values correlate positively (Figure 9). Despite the four samples having different lithologies, climate and uplift conditions (dashed box in Figure 9), a polynomic equation is defined for the rest of samples (Trend A = red line in Figure 9). Trend A reflects samples with similar lithologies and in areas of similar precipitation and uplift. This trend is linear and can be expressed in the form: where a 1 is 12.28. Equation (4) might be used as a proxy for estimating erosion rates when TCN data are not available or cannot be obtained. However, the applicability of Equation (4) for the shale-rich western flank of the Eastern Cordillera of Colombia might be problematic. This is because the western flank is a shaledominated area with more erodible rocks while the eastern flank is composed of alternating sandstone and shale formations that are more resistant; hence application of Equation (4) for the western side will provide underestimates of the erosion rates.
Local climatic, bedrock and tectonic effects Classifying our samples into four groups, according to their erosion rates, provides an additional means of differentiating between the axial plateau and the cordillera flanks ( Figure 10): Group I represents the plateau domain, with the lower erosion rate values (<20 mm/ka). This area is a relatively little tectonically deformed belt of the cordillera (Mora et al., 2008, and is associated with longitudinal, structure-controlled fluvial drainage. Samples LUCN-26 and 45 are not located on the plateau, but the low relief of these catchments and the longitudinal trend of the rivers suggest they were part of the plateau in the past. Group II samples represent an incised domain, including some of the samples located in the flanks and in the Sogamoso catchment (Suárez River). The samples located in this group are from high-energy rivers that are transverse to the chain, and cut across the structural grain.
Group III includes only two samples, LUCN-54 (167 mm/ka, Sogamoso river) and LUCN-56 (228 mm/ka, Chicamocha river). The difference between erosion rates for these samples and the adjacent catchments, 46 mm/ka (LUCN-53) compared with 64 mm/ka (LUCN-55), correlates with a difference in climate, vegetation cover and rock strength. Basement composed of low-and medium-grade metamorphic rocks is exposed in the upper and medium part of the Chicamocha basin (LUCN-56), and the region has little vegetation. In contrast, the Suárez river basin (LUCN-55) mainly traverses alternating Cretaceous sandstone and shale formations. In terms of erodibility, and according to Castro (1992) and González and Jiménez (2015), basement rocks are more erodible than Cretaceous sedimentary rocks. Moreover, erosion is enhanced because of the lack of vegetation cover in the upper and middle part of the Chicamocha basin. The Chicamocha valley has an arid microclimate unlike any other part of the Eastern Cordillera. The valley is very narrow and the obtained data ( Figure 2) from the few pluviometer stations do not reflect the real conditions inside it, resulting in overestimates of precipitation. The Cocuy-Santander massif to the east and the Floresta massif to the west bound the Chicamoca valley (Figure 1). These massifs form an orographic barrier and result in a rain shadow zone within the Chicamocha valley. The Chicamocha river is characterized by a very high mean Mx (~8.13) with steep hillslopes and an erosion rate of 228 AE 32 mm/ka. The Suárez river (sample LUCN-54) mainly flows across Cretaceous formations and Mx values are low (~5.92), with an erosion rate of 63.8 AE 8.6 mm/ka. The Chicamocha catchment, with higher erosion rate, has lower precipitation than the Suárez catchment, which is similar to the Sabana de Bogotá. If high precipitation results in large erosion rates, then higher erosion rates would be expected in the wetter Suárez catchment than in the Chicamocha catchment; but this is not the case. We propose that the high erosion rate in the Chicamocha valley is related to the local 'semiarid' conditions, low vegetation cover and more erodible bedrock, all of which will enhance the erosion capacity. We underline the importance of rock type in determining the magnitude of the erosion within similar topographic domains, as shown in other studies in dynamic settings where higher denudation rates that have been related to high erodibility of the rocks (Salgado et al., 2008;Chadwick et al., 2013;Bierman et al., 2014;Pupim et al., 2015). Safran et al. (2005) showed that the rock type might play a secondary non-dominating role in catchment-averaged erosion rates in the Bolivian Andes, where high erosion rates were found in catchments that have weak or resistant bedrock. That work proposed that uplift had a first-order effect on erosion rates, with rock type playing a second-order effect, and climate had little if any influence on erosion. However, rock type, vegetative cover and climate have a first-order effect on erosion in the Eastern Cordillera of Colombia.
Group IV is composed only by two samples, LUCN-05 (479 mm/ka) and LUCN-06 (670 mm/ka), which yielded strongly different erosion rates compared with the rest of the Guayuriba catchment samples, as for example the sample LUCN-03 with an erosion of 73 mm/ka. Group IV samples are located in the Quetame basement massif (that is easily eroded), characterized by higher precipitation values and younger exhumation ages (Mora et al., 2008Parra et al., 2009a) than in the upper part of the catchment, as indicated by greater erosion and incision (see Figure 2). The combination of greater precipitation and exhumation defines a closed feedback erosive cycle, in accordance with the high erosion rates determined for this area.
Up to 150 mm/ka differences in erosion rates are observed in catchments within the same topographic domain (Figure 6). The difference between the Suárez (64 AE 9 mm/ka) and Chicamocha (228 AE 32 mm/ka) catchments is probably due to a combination of precipitation and lithologic effects. Erosion rates in the eastern flank catchments are not homogeneous, as is, reflected by the Guayuriba basin showing higher erosion rates, relief and incision values than the Turmequé catchment (Figures 2 and 11).
The contrast in erosion rates described above argues for different fluvial dynamics between the Guayuriba and Turmequé catchments. The upper part of the Turmequé catchment has clearly been captured, as suggested by geomorphic analysis in Struth et al. (2015). Our current work adds new analysis using the χ and SSP values suggesting a rearrangement of the axial longitudinal drainage by transverse rivers that traverse the flanks of the cordillera. The disequilibrium in the χ values for both sides of the eastern and central divides suggests drainage divide migration toward the inner part of the plateau's interior and to the north, respectively. This drainage divide migration is more evident for the Turmequé catchment than for the Guayuriba catchment. This supports the view of Struth et al. (2015) who suggested that the Turmequé area was a captured reentrant.
The different geomorphic characteristics of the eastern drainage divide and between the two catchments are interpreted as the result of the competition between uplift and erosion. In the Guayuriba catchment, the erosive potential results from the high contrast between slopes during mountain building and the active tectonics in the basin that has young thermochronologic ages (Mora et al., 2008Parra et al., 2009a; Figure 2), arguing for high erosive potential to compensate the high uplift in the foothills of the flanks. In this way, the upstream propagation of erosion and then, the capture and divide migration is less probable in a catchment with low exhumation such as the Turmequé catchment. Struth et al. (2015) interpret the different river dynamics between the flanks and the axial plateau of the Eastern Cordillera as a product of mountain building and drainage development. As such there was progressive increase in the regional slope caused by the accumulation of crustal shortening and thickening, as documented for the Moroccan High Atlas (Babault et al., 2012). Essentially, the contrast in regional slopes results in different erosion rates and topographic dynamics across the Eastern Cordillera. Variations in local precipitation, tectonics and bedrock within a basin may also influence its erosion dynamics, and may have played a secondary role in the dynamism of divide migration and landscape evolution.

Conclusions
A smooth axial plateau flanked by steep topographic belts characterizes the Eastern Cordillera of Colombia. New 10 Be TCN data reveal high erosion rates for catchments along the highrelief flanks of the Eastern Cordillera, with a mean value of 70 AE 10 mm/ka (exceeding 400 mm/ka in some catchments). In contrast, the mean erosion rate is 11 AE 1 mm/ka for the low-relief axial plateau. This argues for erosional contrasts between the two domains and a migration towards the plateau of the N-S oriented plateau-flanks drainage divides. Results of digital morphometric analysis, including specific stream power, steepness index and χ values confirm the view that the drainage divide in the Eastern Cordillera is asymmetric and is moving by processes of river capture. Drainage reorganization from longitudinal to transverse, dominated by means of a series of river capture events, will lead to a progressive reduction of the extension of the axial plateau of the Eastern Cordillera. The erosional contrast between the two morphologic domains of the Eastern Cordillera was primarily driven by the increase in the orogen regional slope by progressive accumulation of crustal shortening and thickening. Local climate, tectonics and rock type play a secondary role in controlling the erosion rates and the basin dynamics at a local scale.
Comparison of the TCN-derived erosion rates with the geomorphic digital parameters show positive correlations with SSP and/or Mx. Using our derived Equation (4) the Mx-erosion rate plot for the Eastern Cordillera can be applied to help acquire first-order estimates of erosion rates in areas where there are similar lithologic (layered sedimentary rocks in alternating competent and incompetent formations) and pluviometric characteristics (but where no TCN data are available