Diagnosing climatic and tectonic controls on topography: Eastern flank of the northern Bolivian Andes

Spurred by early insights into the potential role of erosion in geodynamics (e.g., Dahlen and Suppe, 1988; Koons, 1989; Molnar and England, 1990; Beaumont et al., 1992), the tectonic geomorphology community has spent the last couple decades exploring connections between climate and tectonics through erosion. Models that assume that climate strongly affects erosional unloading show a strong link among rainfall, rock uplift rates, and resulting landscape form (e.g., Willett, 1999; Beaumont et al., 2001; Whipple and Meade, 2004; Roe et al., 2006; Whipple, 2009). Despite model predictions, firm demonstration of a tectonic response to climate has remained elusive (Whipple, 2009). Indeed, even the requisite connection between rainfall and erosional efficiency is not fully understood (e.g., Molnar et al., 2006; DiBiase and Whipple, 2011). Moreover, rainfall and relief are often correlated due to orographic enhancement of rainfall (e.g., Anders et al., 2006; Bookhagen and Burbank, 2006; Bookhagen and Strecker, 2008), and areas of high erosion rates are often in areas with both high rainfall rates and high relief. The relative contribution of relief and rainfall to localized high erosion rates remains an open question. Whereas some studies have documented a climatic influence on erosion rates (e.g., Moon et al., 2011; Bookhagen and Strecker, 2012), many have found no detectable relationship between rainfall and erosion rate (e.g., Ahnert, 1970; Riebe et al., 2001; Montgomery and Brandon, 2002; Aalto et al., 2006; Goswami et al., 2012).

Contrasts in the morphology and exhumation history between the wet Bolivian Andes north of 17°S and the dry Bolivian Andes south of 18°S have long been held as a prime example of the hypothesized climatic influence on tectonics, but no consensus has been reached. Numerous studies have argued that climate has played an important role in shaping this contrast in morphology and tectonic history (e.g., Masek et al., 1994; Horton, 1999; Montgomery et al., 2001; Strecker et al., 2007; McQuarrie et al., 2008; Norton and Schlunegger, 2011; Schlunegger et al., 2011; Barnes et al., 2012). However, many researchers remain unconvinced because differences in tectonic style and regional morphology also correlate with, and are arguably driven by, variations in the thickness of sedimentary rock packages associated with the paleogeography of Paleozoic depositional basins (e.g., Sheffels, 1995; Allmendinger and Gubbels, 1996; Baby et al., 1997; Kley, 1999). In addition, there is much debate regarding the inference of plateau surface uplift from isotopic, paleobotanical, and geomorphic records that could be biased by climate change (e.g., Barnes and Ehlers, 2009).

At a smaller scale of observation, a parallel debate has focused on both the uplift history and the relative influence of climate and tectonics on the morphology and patterns of erosion rate across strike in the northern Bolivian Andes. Whereas Barnes et al. (2012) argued that middle Miocene to recent exhumation patterns vary with rainfall in the northern Bolivian Andes, neither Safran et al. (2005) nor Insel et al. (2010) detected a climatic control on millennial-scale erosion rate patterns in this region. Safran et al. (2005) showed that erosion rates correlate best with the channel steepness index and suggested that nonuniform rock uplift drives the observed distribution of erosion rates. Similarly, Aalto et al. (2006) found that 90% of the variability in decadal-scale erosion rates could be explained by relief and lithology, and that rainfall was not significantly correlated with erosion rates. However, these authors acknowledge that the strong correlation between relief and rainfall may have inhibited their ability to isolate a climatic influence. In contrast, Schlunegger et al. (2011) have recently argued that the complex convexo-concave profiles of large rivers draining the Cordillera Real result not from the pattern of rock uplift but rather from the distribution of orographic rainfall revealed by satellite remote sensing (e.g., Bookhagen and Strecker, 2008). The unusual morphology of these rivers, however, could reflect either patterns in rock uplift or a combination of uplift and rainfall patterns, and the relative contributions of these two drivers have not been assessed.

We seek to resolve this debate for the specific example of the northern Bolivian Andes and simultaneously advance understanding of the limits of what can and cannot be deduced about the drivers of erosion from the study of landforms. Despite the divergence of prior interpretations, the details of river profiles, and their relation to the distribution of rainfall and erosion rates in the northern Bolivian Andes present an excellent opportunity to resolve the controversy, shed light on the coupling of climate and tectonics in general, and illustrate best practices for the study of tectonic geomorphology of erosional landscapes. Most critical to our approach is a systematic and theory-guided analysis of river profiles with different drainage area, landscape position relative to regional topographic and rainfall patterns, and erosion rate.

The convex to highly concave shape of the large river profiles highlighted by Schlunegger et al. (2011) is a key observation because relatively few circumstances can produce such landforms. First, the upper convexity could reflect (1) an abrupt downstream increase in rock uplift rate (e.g., Kirby and Whipple, 2001), (2) a transient response to an increase in rock uplift rate (e.g., Kirby and Whipple, 2012), and/or (3) the greater erosive power of glaciers (e.g., Brozović et al., 1997; Brocklehurst and Whipple, 2002). The first would require active deformation on either a west-vergent thrust or east-dipping normal fault at roughly the position of the 3800 m contour, but no candidate structures nor breaks in thermal histories have been recognized (e.g., McQuarrie et al., 2008; Barnes et al., 2012). The second and third alternatives speak to the question of whether there is geomorphic evidence of rapid post–10 Ma surface uplift of the northern Bolivian Andes. Second, the highly concave shape of profiles below the upper convexity could only result from: (1) a spatial gradient (decreasing from west to east) in rock uplift rate (e.g., Kirby and Whipple, 2001; Safran et al., 2005; Gasparini and Brandon, 2011), (2) a spatial gradient in erosional efficiency (increasing from west to east either due to progressively greater runoff or more erodible substrate; e.g., Roe et al., 2002; Duvall et al., 2004; Craddock et al., 2007), (3) a temporal decline or cessation of rock uplift (e.g., Whipple and Tucker, 2002; Baldwin et al., 2003), (4) oversteepening below the knickpoint associated with the upper convexity due to a shortage of abrasive tools (e.g., Whipple, 2004; Chatanantavet and Parker, 2006; Crosby et al., 2007; Gasparini et al., 2007), or (5) some combination of these. Although this list encompasses a wide range of possibility space, we will demonstrate that a combination of topographic analysis and numerical modeling in which the contrasts between adjacent rivers and their smaller tributaries are carefully examined can greatly narrow the range of possible explanations, and even reveal the dominant factor or factors controlling the landscape morphology.

After a brief description of the setting and uplift history (“Setting” section), we present detailed analyses of river profiles over a range of sizes (drainage areas of 10–1000 km²), landscape positions, and substrate lithologies in relation to modern rainfall patterns and millennial-scale erosion rates (“Topographic Analysis” section). The results of these analyses motivate numerical simulations that are used to decipher the relative contributions of climate, tectonics, and lithology (“Numerical Modeling” section). From this analysis, we demonstrate that the spatial and temporal pattern of rock uplift exerts the dominant control on landscape morphology and erosion rate patterns, and that some combination of rainfall, lithology, and transient oversteepening of transverse drainages plays a secondary role on landscape morphology (“Discussion” section). In addition, we develop a general theory on how patterns of channel steepness in rivers of a range of sizes can be compared and contrasted to decipher the dominant evolutionary driver.

We focus on a region of high relief in the northern Bolivian Andes that we call the Beni Escarpment, a topographic step of greater than 3 km that runs from north of the Rio Consata to the vicinity of the Rio La Paz drainage, where it becomes less distinct, fading away to the southeast (Fig. 1). Over a cross-strike distance of 45 km, the mean elevation drops from ∼5 to 1.5 km (Fig. 1, inset topography swath). There is a noticeable change in the elevation distribution of the landscape downslope of the foot of the escarpment, where the downslope gradient in mean, maximum, and minimum elevations abruptly decreases. These changes in topography are highlighted by changes in channel steepness (Fig. 1), local relief, and hillslope angles across the landscape (Fig. 2B).

Changes in the topography are linked to the rainfall pattern revealed in satellite data (Fig. 2A; Bookhagen and Strecker, 2008; Nesbitt and Anders, 2009). The rainfall rate begins to increase ∼20 km NE of the foot of the escarpment and rises to a peak at ∼15 km SW into the escarpment (Fig. 2A, inset), where mean elevation reaches ∼2300 m. At higher mean elevations, rainfall rates decline, but they remain above 1500 mm/yr until the mean elevation reaches ∼3300 m, some ∼25 km SW of the escarpment front (Fig. 2A). Farther SW into the high mountains, rainfall rates continue to decline, dropping below 500 mm/yr at elevations exceeding 4500 m on the crest of the Cordillera Real.

The study area is underlain primarily by Ordovician, Silurian, and Devonian sedimentary rocks that have been variably metamorphosed (Martinez and Tomasi, 1978; Martinez, 1980; Guarachi et al., 2001). Little variability in rock strength at a scale that could influence landscape morphology is recognized in these Paleozoic rock packages (Safran et al., 2005; Syrek and Barnes, 2011). An important exception is the granitic and high-grade metamorphic rocks exposed in the core of the Cordillera Real (Fig. 2B), which have the highest Schmidt hammer readings (a proxy for compressive strength; Selby, 1982) in the region (Safran, 1998; Aalto et al., 2006). With the exceptions of the Rio Consata and Rio La Paz, the large channels that drain across the escarpment incise into these apparently harder rocks, a complication not considered in some previous work in the region, but that is considered here.

Overlying the Paleozoic rocks, there is a section of Neogene gravels known as the 11–7 Ma Cangalli Formation (Fig. 2B; Fornari et al., 1987; Mosolf et al., 2011). The Cangalli Formation rests unconformably over the folded and faulted Paleozoic sequence on a surface with up to 500 m relief and is associated with perched low-relief erosion surfaces that apparently grade to the top of the basin-fill level (Fornari et al., 1987). Besides constraining the timing of the exhumation of granitic rocks in the Cordillera Real and termination of shortening in this part of the orogen (Mosolf et al., 2011), the distribution of the Cangalli Formation carries important implications for the style and timing of recent uplift and deformation. Although we cannot be certain of the accuracy of the mapping of these gravels in detail, the spatial distribution of the Cangalli Formation suggests significant deformation and uplift associated with the Beni escarpment—remnants of these gravels are mapped as capping ridges up to ∼2500 m elevation some ∼20 km SW of the foot of the escarpment (Figs. 2B and 3). In addition, Fornari et al. (1987) reported that the base of the Cangalli Formation preserved along the Rio Tipuani suggests significant postdepositional tilt, as is illustrated schematically in Figure 3.

The surface uplift history of the margin of the Central Andean Plateau in Bolivia is controversial despite many independent lines of evidence that point to significant (>2 km) post–10 Ma surface uplift. Study of the Cenozoic stratigraphy of southern Bolivia indicates significant surface uplift because megafans deposited near sea level in the Eocene and Oligocene are now preserved at plateau elevations (Horton, 2005). Geomorphic evidence consisting of the elevated low-relief 12–9 Ma San Juan de Oro erosion surface, now perched above deeply incised canyons, suggests much of this uplift has been contemporaneous with deformation in the Subandean fold-and-thrust belt since ca. 10 Ma (Gubbels et al., 1993; Kennan et al., 1997; Barke and Lamb, 2006). This geomorphic and geologic evidence accords with paleobotanical evidence of 2–3 km of surface uplift since ca. 14 Ma in southern Bolivia (Gregory-Wodzicki et al., 1998; Graham et al., 2001). The paleobotanical evidence is in turn entirely compatible with two independent isotopic estimates of surface uplift of the Altiplano (Garzione et al., 2006; Ghosh et al., 2006).

Controversy arises, however, because (1) all published estimates of surface uplift carry large uncertainties, (2) the geomorphic evidence presented to date is not definitive, and (3) a combination of late Cenozoic cooling and isotopic changes in rainfall in response to gradual plateau uplift could masquerade as rapid changes in surface elevation (Barnes and Ehlers, 2009; Ehlers and Poulsen, 2009). In accordance with the observed high elevation, low-relief surfaces in the study area (Whipple and Gasparini, 2014), we will show that the details of landscape morphology in the Beni escarpment are most consistent with >3 km of rock uplift relatively recently, or over, at most, the last ∼ 15 m.y.

The goal of our topographic analysis is to illuminate trends in channel form across the Beni escarpment and to relate these trends to current rainfall, local lithology, millennial-scale erosion patterns, and a hypothesized rock uplift pattern. We use a 90 m projected Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) for the topographic analysis (USGS, 2006). We use the TRMM (Tropical Rainfall Measuring Mission) 2A25 data from years 1998–2007 to explore rainfall patterns (Nesbitt and Anders, 2009). Safran et al. (2005) used in situ cosmogenic radionuclide data to calculate watershed-averaged erosion rates in the region, and we use their measured rates. We hypothesize that there is a linear increase in rock uplift rate over a 30-km-wide zone SW of the foot of the Beni escarpment, consistent with the mapped distribution of the 11–7 Ma Cangalli Formation (Fig. 3), and that this trend, along with patterns in rainfall, can be seen in the channel morphology.

Our study focuses on patterns in both the large and small channels that flow across the Beni escarpment. The headwaters of the large channels reach back into the driest parts of the study area, and these channels flow roughly perpendicular to the escarpment front. We refer to the large channels as the transverse channels (labeled in Fig. 2). The transverse channels have an upstream length from the escarpment front of greater than 50 km (Fig. 4). The headwaters of the largest transverse channel, the Consata, have eroded farther back into the Altiplano than the other four transverse channels. We contrast the morphology of the transverse channels with that of five shorter channels that also flow roughly perpendicular to the escarpment front but that do not reach back into the driest parts of the landscape (Fig. 2). We refer to these channels as frontal channels, and they have an upstream length from the front of less than 30 km.

The transverse channels are very steep where they flow across the escarpment, and all of the channels have knickpoints or broad convexities (a downstream increase in channel slope) between elevations of ∼4 and 2.5 km (Fig. 4). The upper elevation of the convexities roughly corresponds with the lowermost snowline during the Last Glacial Maximum (Klein et al., 1999), which is illustrated by the 3800 m contour in Figures 1 and 2, and is generally above the elevation of the highest rainfall rates (blue band in Fig. 2A). The largest transverse channel, the Consata, has the smallest convexity, and below the convexity, the channel is less steep than the other four transverse channels. We explore patterns in channel form using the normalized channel steepness index.

Harder rocks and coarser bed material should result in smaller values of K, although there is not an exact mapping between rock type and K value (Stock and Montgomery, 1999), and much remains to be quantitatively understood about the role of rock strength and the properties of coarse bed material (e.g., Sklar and Dietrich, 2004). Based on observations by Safran (1998), we consider the possibility that the granitic and high-grade metamorphic rocks exposed in the core of the Cordillera Real are harder and more resistant to erosion than all of the surrounding rocks and that an influx of coarse, durable sediment from these rock units may inhibit channel incision for a considerable distance downstream.

Critical factors for discriminating among the potential causes for the highly concave regions in the channels are the relative steepness and concavity of adjacent rivers of different size, position relative to the rainfall pattern, and position relative to rock-type boundaries. Because estimates of local profile concavity are highly sensitive to the length scale over which they are measured (e.g., Wobus et al., 2006a), and because channels of different sizes crossing the same spatial gradient in rock uplift, rock strength, or rainfall will exhibit systematically different concavity indices for the same change in channel steepness across this zone (reflecting simply the different relationships between the logarithm of drainage area and distance), the spatial pattern of k_sn values is the more effective tool for quantifying variations in channel form in response to lateral gradients in boundary conditions.

We explore the k_sn patterns in five large transverse channels and five smaller frontal channels that drain across the Beni escarpment, 23 tributary or subtributary watersheds of the transverse channels in which Safran et al. (2005) measured erosion rates based on concentrations of ¹⁰Be in river sands, and an additional 35 tributary or subtributary watersheds for which we do not have erosion rate data but that are exclusively on the Beni escarpment. The tributary watersheds are smaller, and so they may be more representative of uniform uplift and are more representative of uniform rainfall conditions, in comparison with both the transverse and frontal channels, which have a wider range of rainfall rates and possibly rock uplift rates within their watersheds. We include an additional two channels south of the Rio La Paz where the Beni escarpment has largely lost definition (Fig. 2). Where the escarpment loses definition, rainfall becomes more uniform (Fig. 2A), as may the spatial pattern of rock uplift (Whipple and Gasparini, 2014). The two southern catchments are included in the study exclusively to help assess the potential role of lithologic heterogeneity.

We calculated k_sn across all channels in our study area that have a drainage area greater than 10 km², and we display the results from all the channels, but we focus primarily on the patterns from our 70 study channels and tributary watersheds. We used a freely available tool that calculates k_sn over a moving window in all channels greater than the prescribed minimum size (http://www.geomorphtools.org). Channel slopes are smoothed over a 2 km moving window, and k_sn is calculated over a 2 km distance at every point in the channel. When illustrating only the data from the Bolivian channels, θ_ref is set to 0.45, a value commonly used in other studies and therefore convenient for comparisons with other locations (Kirby and Whipple, 2012). When comparing the Bolivian channels with landscape evolution model data, θ_ref is set to 0.5 because the intrinsic concavity (m/n) of the model channels is 0.5. Because there is a linear relationship between k_sn values calculated with θ_ref = 0.45 and θ_ref = 0.5, patterns in k_sn with distance are the same regardless of the reference concavity value (Fig. 5).

In the 23 watersheds that Safran et al. (2005) sampled, we compare the relationships among k_sn, distance from the escarpment front, watershed-averaged rainfall rate, and erosion rate. Safran et al. (2005) investigated trends in erosion rate with normalized channel steepness index calculated using a slightly different method than we use here, and we repeat their analysis using our measurements of k_sn. We use a subset of the data collected by Safran et al. (2005). We only consider watersheds that are tributaries, or subtributaries, to a transverse channel because they integrate erosion rates over an area in which rainfall and erosion rates are likely less variable. We also do not include data from watersheds with a drainage area less than 10 km².

For the analysis of the Safran et al. (2005) watersheds and the additional 35 watersheds, we follow Ouimet et al. (2009) and DiBiase et al. (2010) by using the mean k_sn value across all reaches in the watershed with drainage area greater than 1 km² to represent the normalized channel steepness of the watershed. In the 35 additional watersheds, we only include those stretches of the tributary channel that are below 3800 m in order to ensure that we sample the adjusted or adjusting portions of the landscape that we can confidently interpret to be fluvially dominated. The mean annual rainfall across all of the tributary watersheds is also calculated.

All of our data analysis and modeling are based on the modern rainfall data. What matters most for our study is the pattern of rainfall, which is dictated by orographic precipitation mechanics, and this pattern is unlikely to change over the time scale of our erosion rate data. As we will show, our analysis indicates a secondary dependence of topography on rainfall, which minimizes the sensitivity of our analysis to possible paleoclimate change. Further, we fully acknowledge that climate and orographic precipitation enhancement will change at the million-year time scale of surface uplift. By using the modern rainfall pattern in the numerical modeling scenarios, we maximize the predicted influence of rainfall.

There is a zone of high channel steepness values (the highest values, k_sn > 500 m^0.9) roughly parallel to the escarpment front (Fig. 1). In four of the five transverse channels, the highest k_sn values are just upstream of (or reach slightly into) the highest rainfall zone (Fig. 2A) and just downstream of the 3800 m contour (Figs. 1 and 2). However, the Consata, which also flows through the zone of high rainfall, does not have a similar reach of high k_sn values. Many of the tributaries to the transverse channels (including tributaries to the Consata; Fig. 1) also have k_sn values above 500 m^0.9 in this zone. Channel steepness values decline toward the escarpment front in all of the transverse channels and in their tributaries. Most of the frontal channels have k_sn values below 350 m^0.9; there are only a few short reaches with k_sn values above 350 m^0.9, and none with k_sn values above 500 m^0.9 (Fig. 2). Downslope of the escarpment front, the channel steepness index is generally less than 150 m^0.9 (Fig. 1).

Channel slope versus drainage area data from the transverse channels have some similar trends (Fig. 6). All of the transverse channels have a notable convexity, identified as a region in which slope increases with drainage area. In the Consata, k_sn values increase from ∼100 to ∼350 m^0.9 across the convexity (Fig. 6A), whereas in the Tipuani, Challana tributary, Challana, and Zongo Rivers (referred to herein as the high k_sn transverse channels), k_sn values increase from ∼100 to ≥500 m^0.9 (Figs. 6B–6E). The high k_sn transverse channels incise into regions with granite and high-grade metamorphic rocks, whereas the Consata does not (Fig. 2B). In all of the transverse channels, the convexity occurs in a region where rainfall rates are less than 1000 mm/yr. In the Zongo and Consata Rivers, this convexity occurs well upstream of the point where the channel enters the region where rainfall rates are above 1000 mm/yr, but in the other three transverse channels, the downstream end of the convexity appears to roughly coincide with the location where the channel enters the region of rainfall rates above 1000 mm/yr. In all of the transverse channels, the convexity occurs between drainage area values of 10 and 100 km² (see also Fig. 4). All of the channels are primarily concave downstream of this convexity. They have a region in which the concavity remains close to the reference value of 0.45 and k_sn is relatively uniform, and then downstream of this region, the concavity increases, and k_sn values decrease back to ∼100 m^0.9 at the foot of the Beni escarpment.

The five frontal channels do not have the convexity that is present in the transverse channels (Fig. 6F). They have a relatively constant slope (zero concavity) until ∼5 km², a feature common to steep mountainous drainages that is interpreted as a dominance of debris-flow processes at small drainage area (e.g., Montgomery and Foufoula-Georgiou, 1993; Lague and Davy, 2003; Stock and Dietrich, 2003). The short interval over which k_sn increases downstream on these channels is artifactual, simply reflecting the low concavity in this reach (less than the reference value of 0.45). At drainage areas exceeding 5 km², the frontal channels have a relatively uniform, high concavity. Because the channel concavity is higher than the reference value, k_sn decreases downstream.

The patterns in the slope-area data are illustrated in the k_sn trends with distance from the escarpment front (Fig. 7). The four high k_sn transverse channels have a region in which k_sn increases toward the escarpment front (roughly between 40 and 30 km), a region in which k_sn remains high or decreases (roughly between 30 and 20 km), and a region in which k_sn decreases (roughly between 20 and 0 km) (Fig. 7A). These three regions correspond to the convex, concave, and high-concavity regions, respectively, illustrated in the slope-area data (Figs. 6B–6E and 6G). Even though k_sn increases and then decreases from ∼40 to 15 km, the rainfall rate increases throughout this zone toward the escarpment front. From ∼15 km to the escarpment front, k_sn and rainfall rates both decrease. Downslope of the escarpment, k_sn values remain low, and rainfall rates continue to decrease with distance from the escarpment. The dark-gray band in Figures 7A and 7B highlights the pattern in the four high k_sn transverse channels.

Similar to the high k_sn transverse channels, k_sn decreases toward the escarpment front in the five frontal channels on the Beni escarpment (gray symbols in Fig. 7A). However, k_sn values in the frontal channels are lower than those in the high k_sn transverse channels. The five frontal channels on the Beni escarpment do not incise into any areas of mapped granite, and only their headwaters reach into the highly metamorphosed rocks. The light-gray band in Figures 7A and 7B highlights the k_sn pattern across the escarpment in the frontal channels. Downslope of the escarpment, k_sn values in the frontal channels are similar to the nearly uniform values in the transverse channels.

The k_sn pattern in the Consata, a transverse channel that does not incise into any of the mapped areas of granite or high-grade metamorphic rocks (Fig. 2B) and has a larger drainage area than the high k_sn channels, differs from the other transverse channels (Fig. 7B). There is an increase, leveling off, and decrease in the k_sn values in the Consata, but the locations of these changes are farther away from the escarpment front than in the high k_sn channels. Further, k_sn never goes above ∼400 m^0.9 in the Consata, in comparison with maximum values above 600 m^0.9 in the high k_sn transverse channels.

The two channels south of the Beni escarpment drain an area with nearly uniform rainfall, are immediately adjacent to each other, are of similar size, but one drains an area of mapped granite, and the other does not (Figs. 2 and 7C). Although there is a great deal of scatter, channel steepness values do not appear to vary systematically with distance from the reference line in these two channels. However, at a given distance from the reference line, the channel that flows across an area of mapped granite has a higher k_sn value than the channel that does not drain any of the mapped granite. On average, the granite channel has a k_sn value that is 1.6–1.9 times higher than the nongranite channel (depending on the area over which k_sn is averaged in the granite channel), which can be taken as an indication of the magnitude of the potential lithologic effect.

Mean channel steepness indices, mean annual rainfall, and erosion rates in the sampled tributary watersheds from Safran et al. (2005) all have noticeable trends with distance from the reference line along the foot of the escarpment (see Fig. 8 for watershed locations and Fig. 9 for the data). Mean k_sn and erosion rate in the tributary watersheds generally increase with distance upslope (SW) of the reference line (Figs. 9A and 9C). The value of k_sn remains low downslope of the reference line, where erosion rates are also low but rainfall rates are high (Fig. 9B). Rainfall rates are lowest at a distance of 30 km (and greater) upslope of the reference line, where the highest erosion rates are measured.

Several interesting trends are observed among tributary watershed k_sn values, mean annual rainfall, and erosion rate. Converse to expectation and common assumptions, erosion rates generally decrease with increasing watershed-averaged annual rainfall (Fig. 9D). Erosion rates do, however, generally increase with mean k_sn (Fig. 9E), as predicted by Equation 3. Although there are not many points to generate trends, it is worth highlighting that, for a given mean k_sn, erosion rates are higher in the driest watersheds (illustrated by the fit lines in Fig. 9E). This is counter to what is predicted by the theoretical relationship (Eq. 3). We note that the rainfall pattern varies southeast of the Beni escarpment, and the rock uplift pattern in this area likely also varies (Whipple and Gasparini, 2014), and this could contribute to some of the scatter in the trends.

Mean channel steepness values from the 35 tributary watersheds on the Beni escarpment could be indicative of a spatial rock uplift gradient, or a change in k_sn with lithology, but not both (Fig. 10). Although channel steepness is generally higher in watersheds underlain by granitic or high-grade metamorphic rocks (290 and 420 m^0.9) than in those underlain by weaker sedimentary and low-grade metasediments, the weaker rocks experience greater rainfall rates (Fig. 10C) and exhibit a well-defined increase in mean k_sn from ∼150 to ∼300 m^0.9 over a 20 km distance upslope of the escarpment front (Fig. 10B). The higher channel steepness values in the stronger rocks fall on the same trend of increasing k_sn with distance from the range front, reaching ∼400 m^0.9 at 30 km upslope (Fig. 10B), suggesting a common response to an uplift gradient as suggested by the pattern in erosion rates and deformation of the Cangalli Formation (Figs. 9C and 3, respectively). However, the difference in mean channel steepness between stronger and weaker rocks is commensurate (approximately a factor of 2) with the contrast seen in the two channels south of the Beni escarpment, which appears to be attributable to a lithologic effect (Fig. 7C), and the difference in means in the tributary watersheds is supported by a one-way analysis of variance (ANOVA) test (p < 0.0001). Interestingly, the tributary channels entirely within granitic and high-grade metamorphic rocks have significantly lower k_sn values than the main-stem transverse drainages in the same landscape position, even where the mean rainfall rate is less than 1000 mm/yr (Fig. 10). This suggests that the very high k_sn on the transverse channels is unlikely to be purely a lithological effect.

The relationship between rainfall rates and k_sn in the tributary watersheds is similarly complex. As seen in the subset of watersheds studied by Safran et al. (2005), k_sn values decline with average annual rainfall (Fig. 10C), although another possible interpretation is that there is a threshold at ∼2000 mm annual rainfall, resulting in channels with generally lower channel steepness values above this threshold. However, the wettest watershed has a similar k_sn value (308 m^0.9) as the driest watershed (316 m^0.9), and the apparent threshold at ∼2000 mm/yr coincides with the lithologic boundary discussed already. Erosion rate patterns and the distribution of the Cangalli Formation (Figs. 9 and 3, respectively) suggest an underlying tectonic driver for the spatial distribution of tributary channel steepnesses. The tributary morphometric data are arguably consistent with this interpretation, but not uniquely so. The strong spatial correlations among rainfall rate, lithology, channel steepness, and erosion rate decidedly complicate interpretation. Fortunately, careful, systematic investigation with a landscape evolution model can help us to determine the scenario that best matches all available data.

Based on our topographic analysis, we make the following initial interpretations:

• (1) The foot of the Beni escarpment marks an uplift boundary. The channels and tributary watersheds downslope of the escarpment front have uniformly lower k_sn values and erosion rates in comparison with the channels on the escarpment.

• (2) The change in rock uplift rate at the base of the escarpment is likely an uplift gradient (increasing to the west) and not an abrupt step change in rock uplift rate as would be expected from an active thrust fault (no fault is mapped). The erosion rate and k_sn patterns, along with the slope on the surface of the Cangalli Formation, support an uplift gradient.

• (3) The upper convexity in the transverse channels is unlikely to record a downstream (eastward) increase in rock uplift rate: (a) Knickpoint positions are defined by drainage area and elevation (∼3800 m) and do not trace a likely structural boundary (Fig. 1); (b) deformation within the Eastern Cordillera is thought to have ceased by late Miocene time (e.g., Barnes et al., 2012; Lease and Ehlers, 2013); and (c) neither cosmogenic nor thermochronometric data indicate a sudden change in erosion rate at this position (Fig. 9; Safran et al., 2005; Insel et al., 2010). Thus, these upper convexities are either transient features, or the low k_sn regions at high elevations have reduced slopes as a result of glacial erosion, or both. These upper convexities do not appear to be controlled by the rainfall pattern as suggested by Schlunegger et al. (2011) because there is no consistent relationship between the position of the convexities and the modern rainfall pattern (Fig. 6); this hypothesis will be more rigorously tested in our numerical simulations.

(4) The difference in k_sn values between transverse and frontal channels that incise into the granite/high-grade metamorphic rocks and those that do not suggests that lithology may play a significant role in landscape morphology. However the relationships among channel steepness, lithology, rainfall rate, and erosion rate in tributary channels are complex, leaving uncertain the relative roles of lithology, tectonics, and rainfall rate (Figs. 9 and 10). Additional analyses are required to make further progress.

We use a numerical model to more quantitatively assess the dominant drivers of landscape evolution across the Beni escarpment and to extract a general understanding of how, and to what limits, the dominant drivers can be inferred from study of topography and erosion rate patterns. We use the models to test the following questions: (1) Has the landscape reached steady state, or is the morphology indicative of a transient response (or in other words, has there been significant uplift since the late Miocene?)? (2) Could rainfall patterns alone cause a large convexity in these channels, as suggested by Schlunegger et al. (2011)? (3) Can an uplift gradient alone cause the observed channel morphology? (4) How much of the variation in channel steepness values between the frontal and transverse channels could be explained by lithology, or variations in rock strength?

We use the CHILD landscape evolution model (Tucker et al., 2001a, 2001b) to explore the topographic signatures of different uplift and rainfall scenarios and to explore the potential role of lithologic variations. We perform a number of numerical experiments in which we systematically vary the spatial rock uplift pattern, spatial rainfall pattern, or both the rainfall and uplift patterns, in order to illuminate the signature that each of these forces could have on landscape evolution. The two different uplift patterns that we use are designed to test whether there is a gradual or step change in uplift rate at the foot of the escarpment. Our experiments use two different rock strength values to explore the possible influence of variable lithology. Only fluvial processes are modeled, or in other words, we do not resolve hillslopes in our model, and glacial processes are not simulated. Because we do not simulate glacial erosion, we are only modeling scenarios in which the upper convexity reflects the rainfall gradient, a transient adjustment to a change in rock uplift rate, or both. The potential role of glacial processes is discussed in section on “Role of Tectonics, Rainfall, and Lithology on the Morphology of the Beni Escarpment.” We include both detachment-limited and transport-limited fluvial processes in all of the model simulations, which allow us to capture sediment deposition in channels with large drainage areas that may be inundated with large volumes of sediment during transient conditions (Whipple and Tucker, 2002; Baldwin et al., 2003).

The geometry of the modeled landscape is designed to explore only the 45 km upslope of the escarpment front, which would encompass approximately the drainage divide at the edge of the plateau to the foot of the escarpment for most of the transverse channels (upslope/southwest portion of the swath topography in Fig. 1). All of the modeled landscapes have the same dimensions, resolution, and boundary conditions (Fig. 11A). The spacing between points on the landscape is ∼200 m.

The experiment details are given in Table 1. Two of the experiments (X1 and X2) are designed to contrast the effects of rainfall patterns on steady-state morphology, but our topographic analysis suggests that the landscape is transient. These two experiments start from the same initial nearly flat topography with 0.5 m of white noise. We present six transient experiments (X3–X8), all of which start from the topography resulting from X1, a low-uplift, higher-erodibility, low-relief, steady-state topography. Because all of the simulated landscapes evolve from the same initial white-noise surface, they all have the same channel network. Four of the six transient runs are focused on exploring the four high k_sn transverse channels between the Consata and La Paz Rivers (Fig. 1) that largely define the overall morphology of the landscape. These channels incise into the granitic and highly metamorphosed units, and the simulations use lower erodibility (K and K_t) values (X3–X6). An additional two experiments with higher K and K_t values explore whether a lithologic difference may affect the difference in steepness values between the frontal and transverse channels (X7 and X8) by comparing the morphology of the frontal channels from these runs with that of the transverse channels from comparable runs with lower erodibility (X5–X6).

The transient experiments (X3–X8) are designed to explore three of the four plausible conditions discussed earlier that are necessary to create the broad convexity and general k_sn patterns that are observed in the channels on the Beni escarpment. We discuss, but do not attempt to model, the potential role of a shortage of abrasive tools in generating the very high steepness index values along the transverse drainages given the limitations of current tools-and-cover models to capture the observed transient behavior of many large rivers (e.g., Gasparini et al., 2007; Kirby and Whipple, 2012). In our transient experiments, we systematically vary the spatial uplift pattern, temporal uplift pattern, rainfall pattern, and rock erodibility in order to quantify the effects that each of these variables have on landscape morphology.

In the model, rainfall gradients drive discharge gradients, which may alter both the transient erosion pattern (Han et al., 2014; Colberg and Anders, 2014) and the steady-state channel slope (e.g., Craddock et al., 2007). We explore two different rainfall patterns (Figs. 11C and 11D) in the model experiments, uniform and nonuniform, based on the observed pattern across the Beni escarpment (Fig. 2A). In the transient experiments that include nonuniform rainfall, the rainfall pattern is implemented at the start of the transient simulation. Ehlers and Poulsen (2009) found that the current orographic rainfall pattern was not established until the topography reached 75% of its current elevation; therefore, if anything, we likely overestimate the effects of the rainfall pattern on the channel morphology. The equations for the nonuniform pattern are given in the supplemental material (see ^{footnote 1}).

We explore two different rock uplift patterns in the model experiments: uniform and a ramp-uplift pattern. The ramp-uplift pattern is a linear increase in rock uplift rates for the first 30 km into the escarpment, at which point the rock uplift rate becomes uniform. The equations describing the uplift pattern are given in the supplemental material (see ^{footnote 1}). All of the transient simulations explore the consequences of a total of 3.3 km of rock uplift; this amount is dictated by the ∼3.8 km elevation of the upper convexity in transverse channels, ∼3.3 km above the ∼0.5 km elevation of rivers exiting the foot of the Beni escarpment. There is a direct trade-off between rock uplift rate and erosional efficiency (K, K_t), and we calibrate K and K_t, based on the uplift values, to best match observed channel steepness and relief. In other words, the outcomes of our numerical experiments do not and cannot constrain either the rate of rock uplift or the duration of the simulated young phase of uplift. For simplicity, we set the maximum rock uplift rate to 1 mm/yr, and consequently we simulate a 3.3 m.y. duration of rock uplift. Note that the simulation that explores the effect of a decline in uplift rates through time (X3) is run for slightly longer after the initial 3.3 km of uplift, when we impose a decrease in rock uplift rate, and the topography begins to decay.

When comparing between the model and Bolivian data, we use a reference concavity value of 0.5. The model channels have a known reference concavity value (m/n), and the results are easier to interpret when this value is used. Because the channel steepness index values are sensitive to the reference concavity, for comparison with model output, we recalculate the Bolivian data with the same reference concavity used in the model (θ_ref = 0.5). As noted earlier, the trends in channel steepness in the Bolivian data are the same regardless of the reference concavity value (Fig. 5).

The topographies of the two steady-state, uniform-uplift landscapes, one with uniform rainfall (X1) and the other with nonuniform rainfall (X2), are very different, as evidenced simply by looking at the topographic maps (Figs. 11A and 11B). In the landscape with nonuniform rainfall, less of the topography is at the highest elevations (less red and yellow areas in the topographic map). The differences in elevation distribution across the landscapes are clearly illustrated in the topographic swaths (Figs. 11E and 11F). In the landscape with uniform rainfall, the maximum elevation gradually decreases from the back of the mountain toward the front (Fig. 11E). In contrast, in the landscape with nonuniform rainfall, there is an abrupt decline in the maximum elevation of the landscape at the point where the rainfall rate begins to increase toward the mountain front (at 40 km) (Fig. 11F). In both cases, the modeled elevation distributions are noticeably different from that observed in the Beni escarpment.

The rainfall pattern drives a region of increased concavity in the modeled steady-state transverse channel (Fig. 12, black solid line and squares) that bears some similarity to the observed region of increased concavity in the transverse channels on the Beni escarpment (Craddock et al. [2007] observed similar patterns in the Marsyandi River in Nepal). In the area where rainfall rates are nonuniform in the modeled landscape (between 0 and 40 km from the mountain front), the upstream spatially averaged rainfall rate increases toward the mountain front. As a result, k_sn declines downstream (Eq. 3), and channel concavity is greater than the reference value, even in the region where rainfall rates decline toward the mountain front in the transverse channel. This pattern matches the observed pattern, but it is more subdued. Models predict only a factor of 2 decline in k_sn along the transverse drainages, i.e., far less than observed. Similarly, the nonuniform rainfall model reproduces some key attributes of the frontal channels but not others. The higher mean rainfall rate in the modeled frontal channels yields frontal channels that are less steep than the transverse channels, as observed. However, the upstream spatially averaged rainfall rate in the frontal channels decreases toward the mountain front, and this results in a slight downstream increase in k_sn (Fig. 12, black dashed line and x symbols), i.e., the opposite of what is observed on the Beni escarpment (Fig. 7). The clear implication is that a rainfall gradient alone cannot explain the morphology of the Beni escarpment, but that the rainfall gradient may contribute to differences in the k_sn values (at a given distance from the escarpment front) between the frontal and transverse channels.

Neither of the steady-state, uniform-uplift-rate scenarios results in landscapes with much resemblance to the upper reaches of the Beni escarpment. Both scenarios result in transverse channels with no knickpoints or convexities (Fig. 12). Contrary to the argument put forward by Schlunegger et al. (2011), the rainfall pattern alone cannot create a convexity at steady state, nor does one emerge during the evolution toward this steady state. Moreover, any enhancement of rock uplift in the high rainfall zone, such as that suggested by Schlunegger et al. (2011), would act to damp the high concavity of the modeled transverse drainages, which is one of the ways in which the nonuniform rainfall simulation actually bears a resemblance to the topography of the Beni escarpment. Because there are no convexities in the transverse channels in the modeled steady-state landscapes, no evidence for an active structure that coincides with observed convexities, and a general correspondence of the convexities with the 3800 m contour, we conclude that the convexity in the transverse channels on the Beni escarpment is either caused by a transient response to an increase in rock uplift rate or glacial erosion. Because we do not simulate glacial erosion, the rest of the numerical experiments explore different transient scenarios to determine possible drivers responsible for the channel morphologies observed on the Beni escarpment.

Experiment X3 isolates the effects of a temporal reduction in a uniform rock uplift rate and has some similarities to the Beni escarpment topography but also some notable differences. The trend in the minimum elevation of the modeled topography is similar to the Bolivian topography, but the mean and maximum modeled elevations are too high across the entire landscape (Fig. 13A). The decline in uplift rate over time does result in a reduction in k_sn values toward the escarpment front (green data in Fig. 14A). However, this decline is very localized in the modeled channels, and as a result, the channels are more concave in this region than the Bolivian channels (supplemental Fig. 1 [see ^{footnote 1}]). In addition, the larger, transverse channels respond more quickly than the frontal channels, and the slope in the transverse channels near the mountain front quickly declines toward the transport slope necessary to move the sediment delivered from upstream (Whipple and Tucker, 2002). This is evident in the nearly constant k_sn values near the mountain front in the transverse channel (green squares in Fig. 14A). The frontal channels remain steeper because they do not respond as quickly as the transverse channels (green x symbols in Fig. 14A). As a result, near the mountain front, there is a region in which the frontal channels have higher k_sn values than the transverse channels, which is the opposite of the trend observed in the channels on the Beni escarpment. These results suggest that a temporal decline in rock uplift is not responsible for the morphology of the Beni escarpment.

Experiment X4 isolates the effects of the rainfall pattern on a transient landscape, and it also has noticeable contrasts with the real landscape. The mean, peak, and minimum elevations are higher than the normalized Beni escarpment topography near the mountain front (Fig. 13B). The rainfall pattern does result in a decrease in k_sn values near the mountain front in the transverse channels (blue squares in Fig. 14A). However, the decline in k_sn is not as steep in the modeled transverse channel as it is in the Beni escarpment transverse channels, and as a result, the modeled channels are not as concave as the Beni escarpment channels (supplemental Fig. 2 [see ^{footnote 1}]). The rainfall pattern causes the frontal channels to have lower k_sn values than the transverse channels at the same perpendicular distance from the mountain front, but k_sn increases toward the front in the modeled frontal channels (same patterns as observed in steady-state experiment X2). Both the transverse and frontal modeled channels have an upper convexity. The upper convexity in the modeled transverse channel is much more discrete (occurs over a shorter channel length) than the convexities in the Beni escarpment transverse channels. Unlike the modeled frontal channel, the Beni escarpment frontal channels do not have an upper convexity. Overall, a transient landscape resulting from a uniform increase in rock uplift and the onset of enhanced orographic rainfall does not appear to match the morphology of the Beni escarpment.

The ramp-uplift scenario with uniform rainfall (X5) has a number of similarities with the Beni escarpment. Peak elevations increase linearly upslope of the escarpment front and then almost level off (Fig. 13C). A similar pattern is observed across the Beni escarpment, although the location at which peak elevations level off differs between the two landscapes. Mean elevations rise more linearly and are slightly higher in the modeled landscape in comparison with the Beni escarpment. The simulated and observed transverse channels are generally similar in form, and the location of the channel convexity in the modeled transverse channel is nearly the same as in the Challana channel (Supplemental Fig. 3 [see ^{footnote 1}]). The general trend in k_sn is similar between the modeled and Beni escarpment channels; that is, both modeled and real channels have a linear rise in k_sn into the escarpment (red data in Fig. 14A). However, we highlight some notable differences between the modeled and real data. The k_sn values of the frontal and transverse channels in the model are the same near the downstream edge where they have adjusted to the uplift pattern, whereas the Beni escarpment frontal channels have smaller k_sn values than the Beni escarpment transverse channels at a given distance from the escarpment front. Further, the upstream portion of the modeled frontal channel has a notable convexity (increase in k_sn) that records the same transient stage of channel profile development as the convexities in the transverse channels, but this is not observed in the frontal channels of the Beni escarpment. This scenario appears to capture many, but not all, of the quantified morphological details of the Beni escarpment.

The simulation with the ramp-uplift scenario and nonuniform rainfall pattern (X6) has even greater similarity with the morphology of the Beni escarpment than the previous experiment (X5). Mean, maximum, and minimum elevations of the model topography rise at a slightly lower gradient into the escarpment than they do in the ramp uplift–only simulation (cf. Figs. 13D and 13C). The pattern in k_sn in the modeled transverse channel nearly matches the pattern in the Challana River (compare blue squares and gray squares in Fig. 14B), and this is also evident in the slope-area and channel profiles (supplemental Fig. 4 [see ^{footnote 1}]). In contrast with the previous experiment, the decline in k_sn toward the mountain front is not linear, and this is because the spatially averaged rainfall rate does not increase linearly toward the mountain front. The modeled frontal channels also have slightly lower k_sn values than the transverse channels, although the difference between k_sn values in the transverse and frontal channels is much less in the modeled topography than in the Beni escarpment. The upper convexity in the frontal channels is much smaller in this experiment (in contrast with X5). The addition of the rainfall pattern with the ramp uplift improves the match with the real data; however, the relative difference between the k_sn values in the transverse and frontal channels is not as large as it is in the real channels.

Runs X7 and X8 were designed to explore the possibility that lithology is responsible for, or contributes to, the distinctly different k_sn values between the transverse and frontal Beni escarpment channels. Although where side-by-side the transverse channels (higher k_sn) are cutting through the same rocks as the frontal channels, the transverse channels are transporting coarse, durable cobbles and boulders derived from the granitic and high-grade metamorphic rocks upstream and may be steeper as a consequence. Rather than attempting to recreate this complex scenario in a model, we simply compare the morphology of the frontal channels from a set of runs with enhanced erodibility (X7–X8) with that of the transverse channels from earlier, comparable runs with lower erodibility (X5–X6). Naturally, the larger erodibility values in experiments X7 and X8 result in lower k_sn values in both the frontal and transverse channels (red and green data in Fig. 14B) in comparison with the earlier ramp-uplift experiments with either uniform rainfall (X5; red data in Fig. 14A) or with nonuniform rainfall (X6; blue data in Fig. 14B). Both X7 and X8 have a ramp-uplift pattern, and X8 also includes the rainfall pattern. The rainfall pattern results in lower k_sn values near the mountain front (green data in Fig. 14B) in comparison with the uniform rainfall case (red data in Fig. 14B). In both X7 and X8, k_sn values in the frontal channels are within the range observed in the Beni escarpment frontal channels, but in the lowest part of the observed range. The transverse channels are less steep than the Beni escarpment transverse channels (see also supplemental Figs. 5 and 6 [see ^{footnote 1}]) but interestingly match the steepness of the Rio Consata—the one transverse drainage that does not tap into the granitic or high-grade metamorphic rocks, breaching the Cordillera Real where there is a gap in the outcrop of these stronger rocks.

It is tempting to assume that the concentration of high rainfall rates on the Beni escarpment must play a major role in setting the morphology of the landscape (e.g., Schlunegger et al., 2011; Barnes et al., 2012), but we have shown that the rainfall pattern alone cannot drive the patterns we observe in the topography of the escarpment and the morphology of the channels flowing across it. Here, we discuss all the lines of evidence that suggest that the morphology of the Beni escarpment is primarily driven by an uplift gradient, and that rainfall patterns play a secondary role in controlling landscape morphology.

Because rainfall gradients affect fluvial incision patterns only indirectly by changing the downstream distribution of fluvial discharge, more extreme rainfall gradients would be required for rainfall to be the sole driver of the observed morphology of the transverse channels. The local rainfall rate does not directly impact the channel slope (Eq. 3). Rather, incision rates are controlled by the fluvial discharge, which is a function of the upstream spatially averaged rainfall rate within the catchment, and not the local rainfall rate. At steady state with uniform rock uplift rates, modeled transverse channels from the nonuniform rainfall simulation (X2) do not have any areas in which slope increases downstream (no upper convexities), and modeled channel concavities are weaker than those in the Beni escarpment (Fig. 12). Similarly, during a transient adjustment to uniform faster rock uplift, modeled transverse channels subjected to nonuniform rainfall (X4) do exhibit the upper convexity as expected; they do have an ∼25-km-wide zone of enhanced concavity, but again the simulated concavities are far weaker than those observed in the Beni escarpment (with k_sn declining by only a factor of ∼2 compared to a factor of ∼5). Further, the morphology of the frontal channels on the Beni escarpment does not match the predicted and modeled morphology of frontal channels without an uplift gradient; i.e., the Beni escarpment frontal channels also express enhanced concavities despite the downstream decrease in rainfall rate in these smaller catchments (Fig. 2).

Schlunegger et al. (2011) suggested that the convexity in the transverse channels between 10 and 100 km² results from high uplift rates driven by high rainfall rates, but this is not consistent with either observations or models. Although the downstream end of the convexity does roughly correspond to the zone of high rainfall rates (>1000 mm/yr) in three of the transverse channels, the convexity is well upstream of the high rainfall zone in the Zongo and Consata channels (Fig. 6). Moreover, we can find no mechanism by which this can occur and conclude that this convexity rather reflects either a transient channel adjustment to an increase in rock uplift rate as in our simulations, or possibly glacial scour at and above the Last Glacial Maximum equilibrium line altitude (e.g., Brozović et al., 1997; Brocklehurst and Whipple, 2002), or most likely both (glacial scour enhancing a transient convexity). This follows because the locations of the upper convexity in all five transverse channels (including the Consata) are most consistent with the transient adjustment to uplift scenario, as summarized in the following text. The fact that the convexity is in the same approximate drainage area range in all of the transverse channels supports our interpretation that this convexity was driven by a wave of increased incision that was spurred by increased uplift rates. Such a wave of incision will migrate upstream at a rate proportional to the local drainage area and would be expected to be at roughly the same drainage area (Whipple and Tucker, 1999) and elevation (Niemann et al., 2001) in all of the transverse channels of the same size. The upstream convexities in modeled transverse channels with low K values, ramp uplift, and either with or without the rainfall gradient (X5 and X6; supplemental Figs. 3 and 4 [see ^{footnote 1}]) have similar forms to those observed in the high k_sn transverse channels on the Beni escarpment, with the upper convexities occurring at consistent elevations and drainage areas. Although we cannot rule out an important role of glacial scour above ∼3800 m elevation, the fact that the upper convexity in the Consata River occurs at the same elevation and approximate drainage area suggests that the primary control is the amount of transient surface uplift. This follows because the headwaters of the Consata rise onto the low-relief and weakly eroded Central Andean Plateau that lacked the powerful alpine glaciers of the Cordillera Real.

Interpretation of the frontal channels, erosion rate data, tributary morphometric data, and the mapped distribution of the 7–11 Ma Cangalli Formation all further support the hypothesis that an uplift gradient is primarily responsible for the morphology of the Beni escarpment. Channel steepness values in the frontal channels decrease toward the escarpment front (enhanced concavity; Fig. 7A), suggesting that the uplift gradient has a greater impact on channel form than the rainfall gradient, which would be associated with a downstream increase in channel steepness in this case. Millennial-scale erosion rate data collected by Safran et al. (2005) vary systematically with distance from the escarpment front (Fig. 9C), whereas the trend with rainfall rate is counter to what would be expected if higher rainfall rates were driving higher rates of erosion (Fig. 9D). Finally, channel steepness values from 35 additional tributary watersheds also increase with distance from the escarpment front (Fig. 10B). Channel steepness values in the wettest tributaries are generally smaller than those in the driest tributaries, which supports the idea that rainfall rates influence channel morphology (Fig. 10C) but does not support the idea that erosion rates increase with rainfall rate. However, it is difficult to discern from the tributary watershed data alone whether the rainfall pattern or hypothesized uplift pattern is the primary driver of channel morphology, because these two variables likely covary.

Although we argue that the uplift pattern dominates the channel morphology, rainfall does appear to affect both channel steepness values in the larger channels and the overall morphology of the landscape. The transverse channels drain the driest parts of the landscapes, so for a given distance from the escarpment front, the upstream averaged rainfall rate in the transverse channels is lower than in the frontal channels. In agreement with theory (Eq. 3) and model simulations (X2, X4, X6, and X8; Figs. 12 and 14), the Beni escarpment data show that the transverse channels have higher channel steepness values than the frontal channels at a given distance from the escarpment (Fig. 7A). The model results also show that the rainfall gradient reduces mean elevation values near the escarpment front, resulting in a modeled gradient in mean elevation near the escarpment front that more closely resembles that of the real landscape (Fig. 13D) in comparison to the model simulation that does not include the rainfall gradient (Fig. 13C).

Adding the rainfall gradient alone to the transient ramp-uplift scenario (X6), however, does not fully explain important details of the observed landscape morphology: (1) The modeled frontal streams are comparatively slightly too steep, (2) modeled frontal streams exhibit transient upper convexities (blue x symbols in Fig. 14B) not present in the Beni escarpment frontal channels, and (3) there is no mechanism in this model to explain the observation that the steepness of the Consata is more similar to the frontal channels than to the other transverse channels (Fig. 7B). To best match the contrast in steepness between the frontal and high k_sn transverse channels, either (1) the rainfall influence needs to be stronger than incorporated into our model (Eq. 3), or (2) the lithologic difference between these drainages is important (X8, red and green symbols in Fig. 14B), or (3) the contrast in channel steepness reflects a transient oversteepening at and below the upper convexity such as is predicted under some circumstances with tools-and-cover models (e.g., Chatanantavet and Parker, 2006; Gasparini et al., 2007; Crosby et al., 2007).

As noted earlier, available erosion-rate data suggest a weaker relationship with runoff than assumed in our models, so appealing to a stronger rainfall influence does not appear to be a likely explanation for the slight misfit between model predictions and data (Fig. 14B, blue symbols). Moreover, a stronger rainfall influence would be expected to reduce the concavity of the frontal channels, inconsistent with our observations. Similarly, a lithologic impact on channel morphology could explain the difference between the high k_sn transverse channels and the frontal channels (Fig. 14B, red and green symbols); however, a lithologic signal is not consistently observed in all of the data. We do see an approximate twofold difference in k_sn similar to that between the frontal and high k_sn transverse channels on the Beni escarpment in two otherwise similar frontal channels southeast of the Beni escarpment with and without significant granitic rocks (Fig. 7C). However, the k_sn data from the 35 tributary watersheds (Fig. 10B) either reveal the uplift gradient or they are consistent with a lithologic explanation for the disparate k_sn values in the transverse and frontal channels, but not both, and many lines of evidence support the presence of an important uplift rate gradient. Finally, we do not have the requisite data to quantify differences in rock strength, nor field observations to confirm that transverse channels are armored with granitic boulders downstream of granite outcrops, as would be required for the presence of harder rocks in the upper reaches to affect the steepness of the transverse channels all the way down to the foot of the Beni escarpment as observed. Consequently, although we cannot rule it out, we cannot conclude that there is a clear lithologic influence at work.

A transient oversteepening is predicted to occur when the upstream sediment load does not supply enough tools for local incision rates to keep pace with downstream incision rates that have increased due to an increase in rock uplift rate. The oversteepening occurs downstream of a critical area that is a function of the original uplift rate, the factor of increase in uplift rate, and the K and K_t values (see Eq. 32 in Gasparini et al., 2007). The analytical relation in Gasparini et al. (2007) suggests that the parameter values that we use in this study are compatible with a transient oversteepening in the transverse channels—a possibility substantiated by the observation of numerous fluvial hanging valleys (an extreme case of transient oversteepening) southeast of the Beni escarpment (Whipple and Gasparini, 2014). If transient oversteepening downstream of the convexity is contributing to the high k_sn values in four of the five transverse channels, this would be consistent with the lack of lithologic signal in the tributary watershed data (assuming an uplift gradient) and could explain the slight misfit between models and data (Fig. 14B). However, this mechanism is inconsistent with the Rio Consata (Fig. 7B), which lacks granitic outcrops and has much lower steepness values than the other four transverse rivers (Fig. 7B) despite a pronounced upper convexity (Figs. 4 and 6)—it ought to express the same transient oversteepening phenomenon, but it does not.

Thus, we can find no single model or mechanism that can explain all the data in full detail. Consequently, we cannot precisely quantify the relative roles of rainfall and lithology on landscape evolution in the Beni escarpment. Nonetheless, we can conclude with certainty that the dominant influence on channel profiles and the overall morphology of the Beni escarpment is the pattern and history of rock uplift. Moreover, the model with the ramp-uplift pattern plus the rainfall gradient comes fairly close to explaining the data (Fig. 14B, blue symbols), and transient oversteepening is a plausible explanation of the misfit in four of the five transverse channels, without invoking any influence of lithology.

Several studies have suggested that uplift of the Bolivian Altiplano has slowed or stopped (e.g., Garzione et al., 2006; Ghosh et al., 2006; Hoke and Garzione, 2008), but it is difficult to reconcile the landscape morphology of the Beni escarpment with the results from the model simulation (X3) of uplift slowing and cessation for a significant period of time (2 m.y.). The minimum elevations across the modeled landscape have a similar distribution as those across the Beni escarpment; however, the mean and maximum elevations remain too high (Fig. 13A). With time, the mean and maximum elevations would decline, but by that time, the minimum elevations in the center of the modeled escarpment would be lower than those in the Beni escarpment. Further, a cessation in uplift leads not only to channels with nearly uniform slope in the downstream reaches (supplemental Fig. 1 [see ^{footnote 1}]; Fig. 14A, green data), rather than the highly concave channels as are observed in the Beni escarpment, it also leads rapidly to a configuration where transverse rivers have lower k_sn values than the adjacent frontal streams—the opposite of what is observed (Fig. 14).

In contrast, our results suggest that ∼3 km of surface uplift has occurred in at most the last ∼15 m.y., and that if uplift has ceased, the landscape has not yet had time to respond to the cessation in uplift. The transient simulations had 3.3 km of maximum rock uplift in 3.3 m.y., but the simulations constrain only the total amount of uplift, not the uplift rate and duration. Cosmogenic isotope data suggest that the maximum rock uplift rate (>30 km SW from the foot of the Beni escarpment) is ∼0.7 mm/yr (Safran et al., 2005), and the longer-term mean is 0.4 ± 0.29 mm/yr, determined from apatite fission-track thermochronometry (Insel et al., 2010). This lower maximum rock uplift rate suggests a duration of 5–16.5 m.y. of rock uplift, consistent with previous interpretations that the recent phase of rapid rock uplift began ca. 10–12 Ma (e.g., Gillis et al., 2006; Barnes et al., 2012). The age of the Cangalli Formation (11–7 Ma) supports the observation that the uplift is relatively young; otherwise, these deposits would not be deformed and would likely no longer be preserved on ridgelines within the Beni escarpment. The preservation of low-relief surfaces at ∼4000 m (Whipple and Gasparini, 2014) and the age of the Cangalli Formation (11–7 Ma) support the observation that the uplift is relatively young. If the uplift were older, the Cangalli Formation deposits would not be deformed and would likely no longer be preserved on ridgelines within the Beni escarpment. Although we cannot constrain the exact timing of the uplift, we can say that the uplift began after 15 Ma, and if it has stopped, the channels have not yet had time to respond.

The numerical simulations and the topography of the Beni escarpment suggest that this landscape has not fully adjusted to the pattern of rock uplift. The upstream convexity in the transverse channels links a channel in a relatively lower local relief landscape with the downstream portion of the channel that is actively incising into a high local relief landscape. The low-relief parts of the landscape above the upstream convexity are currently above ∼3800 m (Figs. 1, 2, and 4), and this part of the landscape may represent preserved patches of the former low-relief landscape that once stretched across the Beni escarpment before it was an escarpment, as evidenced by the Cangalli Formation (Mosolf et al., 2011). This would also be consistent with the elevated low-relief surfaces that have been observed south of the study area (e.g., Barke and Lamb, 2006) and NE of the Beni escarpment (Whipple and Gasparini, 2014). These low-relief patches are also roughly coplanar with the low-relief erosional margin of the Central Andean Plateau, which is well-developed at ∼4000 m elevation north of the Cordillera Real in southern Peru (Lease and Ehlers, 2013). However, without additional data, we cannot rule out the possibility that these low-relief parts of the landscape were shaped by glaciers, especially in light of the fact that the equilibrium line altitude during the Last Glacial Maximum was at ∼3800 m in this region (Klein et al., 1999).

Our theoretical and numerical results suggest a number of methods that when used in concert can best determine the relative roles of climate, tectonics, and lithology in landscape evolution. We advocate for a combined approach that uses DEM analysis, numerical modeling, and erosion rate data along with geologic maps and hydrometeorologic data as available. When erosion or incision rate data are available, they are particularly helpful for deciphering whether high rainfall rates drive high rock uplift and erosion rates. Patterns in erosion rates, at steady state, should always follow the uplift pattern, and previous studies support that the channel steepness index can be used as a proxy for erosion rate (Safran et al., 2005; Harkins et al., 2007; Ouimet et al., 2009; Cyr et al., 2010; DiBiase et al., 2010). If rainfall and erosion rates also follow the same patterns, this could mean that rainfall rates are driving higher erosion rates, or it may simply mean that rainfall rates covary with uplift rates. However, because rainfall and uplift rates affect topography in different ways, deciphering whether there is a rainfall signal atop of the uplift signal requires topographic metrics along with erosion rate data. In such cases, numerical models are especially useful because it is not easy to quantitatively predict patterns in landscape morphology in response to diverse drivers from theory alone, especially during transient adjustment.

We find that all of the different climatic and tectonic scenarios explored in this study have a unique signature on channel morphology, but only when both transverse and frontal channels are considered together. This is predicted by theory (Eq. 3), confirmed by our numerical experiments, and summarized in Figure 15. When only rock uplift rate varies spatially, k_sn varies directly with the uplift rate and should be the same in both transverse and frontal channels (Fig. 15A). In contrast, when only rainfall rate varies spatially, k_sn decreases downstream where rainfall rates both increase and decrease downstream in the channels that drain the driest parts of the landscape (transverse channels), and k_sn increases downstream in the channels that drain only the wettest parts of the landscape (frontal channels) (Fig. 15B). Further, at a given distance from the mountain front, transverse channels have higher k_sn values than frontal channels in this scenario. Although it may be unlikely to find either of these idealized scenarios in nature, they are useful for interpreting the more complex case of spatially varying rock uplift and rainfall rates (Fig. 15C). Similar to the rock uplift gradient–only case, in scenarios with both uplift and rainfall gradients, k_sn in both transverse and frontal channels declines toward the mountain front where the gradients in forcing influence channel morphology, indicating that the uplift pattern has a stronger control on channel morphology than the rainfall gradient. However, the rainfall gradient results in frontal channels that have lower k_sn values than transverse channels at a given distance from the mountain front. Further, the modeled decline in k_sn in the steady-state channel reaches is not linear due to the nonlinear variation in spatially averaged rainfall rate. A temporal decline in uplift rate also results in a nonlinear decline in k_sn downstream near the mountain front in both transverse and frontal channels (Fig. 15D). However, in this case, at a given distance from the mountain front, the frontal channels have a higher k_sn value than the transverse channels in the model. Theory also predicts that channel steepness downslope of the escarpment can be used to decipher the relative roles of climate and tectonics. Comparisons of channels that drain the driest parts of the escarpment with those that originate downslope of the escarpment (assuming this area is relatively wetter) can also be helpful for quantifying the local influence of climate (Fig. 15B).

The patterns in the Beni escarpment channels, theory, and our numerical results lead to the following observations that can be used to quantitatively compare the influence of tectonics, climate, and lithology in other settings:

• (1) In order to explore the influence of rainfall rates on fluvial discharge and incision rates, it is critical to compare k_sn values in channel reaches that have different upstream averaged rainfall rates. In a landscape with a dominant rainfall gradient in one direction, small watersheds oriented perpendicular to this gradient will each experience relatively uniform rainfall rates. Different rainfall rates among such small watersheds arranged along the rainfall gradient will be particularly useful for deciphering the influence of rainfall on channel morphology. Equally diagnostic will be similarities or differences in profile concavity (best quantified and illustrated on plots of k_sn vs. distance, not on classic slope-area plots) among channels that drain areas with rainfall that is variously uniform, increasing downstream, and decreasing downstream. The power of such comparisons is exemplified on the Beni escarpment by the similar rate of decrease in k_sn with distance in transverse and frontal channels despite their contrasting rainfall patterns.

• (2) Once channels have reached steady state, and assuming all other variables are uniform, k_sn values should always increase with rock uplift rate (Kirby and Whipple, 2012). However, the relationship between k_sn and rock uplift rate is affected by other variables. Ideally, to decipher an uplift pattern, k_sn data would be measured in channels from different locations in the landscape that have the same upstream averaged rainfall rate and lithology. However, if there are other confounding variables, or cross-correlations among variables, a comparison between k_sn values in reaches of large and small channels that have the same hypothesized uplift rates (as diagramed in Fig. 15) should aid in identifying uplift patterns. In all four of the proposed scenarios presented in Figure 15, the combination of the k_sn data from transverse and frontal channels is diagnostic of both the uplift and rainfall pattern, whereas the data from just a transverse or frontal channel are not always diagnostic, even in these four relatively simple cases. In many cases, deciphering the relative influences of different variables on k_sn patterns in large and small catchments will require numerical modeling in which variables are systematically varied, especially when patterns are more complex than those presented in Figure 15.

• (3) Lithology may impact channel steepness indices both locally and downstream. This follows because exposure of harder rocks that produce an abundance of coarse, durable gravels and boulders may increase channel slopes, not just locally, but for some distance downstream. Therefore, when interpreting climate and tectonics from topography, it is important to consider rock strength distributions both locally and upstream.

Based on the combined analysis of DEM data, erosion rate data, and numerical model output, we find that tectonic processes exert the dominant control on the morphology of the Beni escarpment. The morphology of the Beni escarpment suggests that uplift rates increase linearly into the escarpment, but whether the base of the Beni escarpment is a fault, shear zone, or fold and the way in which it relates to deeper structure (décollement ramp, duplex, or basement thrust sheet) are not known. Orographic rainfall in response to surface uplift does appear to play a secondary role in shaping fluvial channel profiles and trends in topography, but we find no measurable evidence for a feedback between rainfall rates and rock uplift rates. We find that the four larger, transverse channels draining harder lithologies (granites and high-grade metamorphic rocks) have steeper profiles in comparison with the other observed transverse or frontal channels, all else being equal. However, assuming a gradient in rock uplift rates, lithology does not appear to impact the morphology of smaller tributary channels. An alternative explanation for the higher channel steepness in four of the transverse channels in comparison with the frontal channels is a transient oversteepening due to a lack of abrasive sediment tools downstream of the convexity (Gasparini et al., 2007). However, this explanation is not consistent with the morphology of the Consata River.

We also demonstrate that diagnosing the relative roles of climate and tectonics in landscape evolution requires a comparison of river profile morphology across drainages with a range in sizes, locations, and orientations relative to structural and climatic gradients. Moreover, a combined analysis of landscape morphology (DEMs), geologic maps, and hydrometeorologic data, along with erosion rate patterns, and output of landscape evolution models, proves to be a powerful approach, especially if models are used to test alternative hypotheses suggested by analysis of field data.

We thank R. Aalto, J. Barnes, S. Nesbitt, and E. Safran for generously sharing data. We acknowledge H. Hoey and J. Adams, who helped on some of the data analyses. We are grateful for thoughtful and thorough critiques from Editor E. Kirby and reviewers B. Bookhagen and S. Miller.