We studied the relationships among present-day relief, precipitation, stream power, seismic energy, seismic strain rate, and long-term exhumation rates for the Venezuelan Andes. Average long-term exhumation rates were determined for seven large catchments in the Venezuelan Andes from fission-track analysis of detrital apatite. A quantitative comparison between eight new detrital apatite fission-track (AFT) age distributions presented here and previously published bedrock AFT age patterns shows that detrital AFT ages can be used for predicting exhumation patterns across the mountain belt. Catchment-averaged exhumation rates estimated from the raw data range from 0.48 ± 0.02 km m.y.−1 to 0.80 ± 0.26 km m.y.−1 Accounting for variable sediment yield and assuming that short-term sediment production rates scale with long-term exhumation rates, these rates vary from 0.33 ± 0.07 km m.y.−1 to 0.48 ± 0.08 km m.y.−1 No variation in rates is observed between the northwestern and southeastern flanks of the mountain belt, despite a threefold increase in precipitation from the northwest to the southeast. Long-term exhumation rates are strongly correlated with relief in the different catchments, weak or negative correlations exist with precipitation data or present-day erosion indexes, while the correlation with seismic energy released by earthquakes is weak to moderate. This lack of correlation may be caused by the insufficient temporal range of the available precipitation and seismicity data, and the different time scales involved in the comparison. Long-term exhumation rates are, however, strongly correlated with seismic strain rates (which take the temporal earthquake magnitude-frequency scaling into account), suggesting that the moderate correlation with seismic energy is indeed related to the different time scales and that tectonic control on exhumation is significant. In contrast, given that precipitation patterns in the Venezuelan Andes should have been installed during Miocene times, we suggest that decoupling of relief and exhumation from present-day climate explains the lack of correlation between exhumation and precipitation.

The relative importance of tectonic and climatic control on relief development and exhumation in mountain belts remains strongly debated (e.g., Molnar and England, 1990; Burbank et al., 2003; Reiners et al., 2003; Strecker et al., 2009; Champagnac et al., 2012). Problems limiting our ability to discriminate between the driving forces for the development of topography and erosion include the limited temporal precision with which variations in exhumation rate, climate, and tectonics can be resolved and compared (e.g., Whipple, 2009), and the different temporal scales for which we have records of these processes (Burbank et al., 2003; Reiners et al., 2003; Vernon et al., 2009). While thermochronologic data record denudation on million-year time scales, seismicity (arguably a proxy for the intensity of at least brittle tectonics) and precipitation records (an admittedly reductionist climate descriptor) generally do not go back more than a few decades. Nevertheless, the comparison of spatial patterns of long-term exhumation rates and short-term seismicity and precipitation may provide some first-order insights into the processes driving exhumation, even in the absence of strong spatial correlations (Finlayson et al., 2002; Dadson et al., 2003; Koons, 2009; Vernon et al., 2009).

The Andes are important in the discussion on climate versus tectonic controls on topography and exhumation in mountain belts because of the obvious link between laterally varying topographic characteristics and climatic zonation (e.g., Montgomery et al., 2001; Lamb and Davis, 2003), as well as the strong asymmetry in orographic precipitation associated with this orogen (e.g., Bookhagen and Strecker, 2008). While many studies have concentrated on the southern and central Andes (e.g., Strecker et al., 2007; Farías et al., 2008; Thomson et al., 2010), attention has recently also focused on the Northern Andes. Mora et al. (2008), for instance, described strongly asymmetric deformation and exhumation patterns in the eastern Colombian Andes, which they linked to the strong NW-SE precipitation gradient across this part of the mountain belt. In contrast, Mora et al. (2009) and Parra et al. (2009) argued for a strong control of inherited crustal structure on spatio-temporal patterns of exhumation and relief production.

In this study, we present new detrital apatite fission-track (AFT) data from seven river catchments in the Venezuelan Andes, the northeasternmost extension of the Northern Andes. We compare the detrital AFT data to available bedrock AFT data of the Venezuelan Andes (Kohn et al., 1984; Bermúdez et al., 2010, 2011) in order to test whether they faithfully record exhumation in the catchments (e.g., Ruhl and Hodges, 2005; Brewer et al., 2006) and to spatially integrate the still relatively sparse in situ data. The objective of this study is to examine the relative controls of climate and tectonics on exhumation and relief development in the Venezuelan Andes. For this purpose, long-term exhumation rates derived from the detrital thermochronology data are compared to present-day relief, precipitation rates, seismic energy release, seismic strain rate, and short-term erosion potential predicted from the average stream power of individual drainages.

The Venezuelan Andes were formed by oblique convergence between the continental Maracaibo block and the South American plate, driven by eastward movement of the Caribbean plate (Fig. 1A, inset; Case et al., 1990). The orogen is characterized by high seismicity and spatially variable exhumation since Miocene times, as recorded by AFT thermochronology data (Kohn et al., 1984; Bermúdez et al., 2010, 2011). The structure of the Venezuelan Andes is controlled by reactivated faults that delineate individual tectonic blocks. These faults were inherited from early Mesozoic rifting, and possibly from earlier orogenesis (Aleman and Ramos, 2000; Pindell and Kennan, 2001). Bermúdez et al. (2010) defined at least seven tectonic blocks with contrasting exhumation and cooling histories, separated by major strike-slip and thrust faults. The orogen is bounded by two seismically active thrust belts to the northwest and southeast (Colletta et al., 1997; Fig. 1B). The most important strike-slip fault systems are the right-lateral Boconó, Central-Sur Andino, and Caparo faults and the left-lateral Icotea, Valera, and Burbusay or Carache fault systems (Fig. 1A). The Boconó fault zone extends from the border with Colombia for more than 500 km to the northeast (Fig. 1A) and divides the Venezuelan Andes symmetrically in its central part into two separate chains; the Sierra La Culata to the northwest and the Sierra Nevada to the southeast are delimited by the NW fold-and-thrust belt to the north, and the Boconó fault to the south; and the Boconó fault, other tectonic block (Caparo, inBermúdez et al., 2010) and SE fold-and-thrust belt to the south, respectively (Fig. 1B). These two blocks cooled rapidly but diachronously during the late Miocene–Pliocene (Kohn et al., 1984; Bermúdez et al., 2011). In contrast, the Caparo tectonic block is delimited by the Boconó fault system and Central-Sur Andino–Caparo faults system, and the Trujillo tectonic block is delimited by the Boconó fault system, Valera fault system, and Burbusay fault (Fig. 1A), at the southwestern and northeastern ends of the Venezuelan Andes, respectively, and both experienced slow cooling from the late Oligocene to late Miocene onward (Bermúdez et al., 2010). Both these blocks are dominated by Paleozoic and Mesozoic sedimentary rocks, whereas Proterozoic to Paleozoic gneisses and granites are exposed in the core of the mountain belt, in the Sierra Nevada and Sierra La Culata ranges. This crystalline basement is covered by Cenozoic sedimentary rocks of the Maracaibo and Barinas foreland basins to the north and south, respectively (Fig. 1).

Thermochronology of detrital minerals allows the quantification of cooling rates and exhumation processes in convergent mountain belts (e.g., Garver et al., 1999; Carter, 2007), providing a complementary record of the erosional history of a mountain belt with respect to local cooling paths deduced from in situ samples within the orogen (Brandon and Vance, 1992). Detrital thermochronology applied to modern river sediments provides a spatially integrated view of the exhumation history of an entire drainage basin, which may be the only way to characterize exhumation patterns for inaccessible areas.

Data Collection and Discrimination of Age Components

We collected samples for detrital AFT thermochronology from seven river catchments draining the Venezuelan Andes. The Coloncito, Tucaní, San Pedro, Agua Viva, and Mimbós Rivers drain the northern flank, the Chama and Chejendé Rivers drain the central range, and the Santo Domingo River drains the southern flank of this mountain belt (Fig. 2). Samples were collected along active river channel bars at all sites and prepared using standard techniques (cf. Bermúdez et al., 2010). We aimed at dating at least 100 grains per sample in order to attain statistically significant AFT age populations. This objective was attained for five out of eight samples (Table 1).

Because apatites in detrital samples are derived from different sources with variable bedrock ages within a drainage area, the grain-age distribution may contain several grain-age components (Brandon and Vance, 1992; Garver et al., 1999). Different methods have been proposed for deconstructing a fission-track grain age distribution into its age components, the most popular of which has become the binomial peak-fitting method (e.g., Galbraith and Green, 1990; Brandon 1992, 1996). Here, we use this method (as described by Stewart and Brandon, 2004) to decompose the grain-age distributions of our modern river samples. Resulting age peaks are reported in Table 1 and Figure 2. All samples except from the Tucaní catchment record a late Miocene (6–10 Ma) age peak, and all except those from the northwesternmost catchments (Agua Viva and Chejendé) also record a Miocene–Pliocene (2.5–6 Ma) age peak. Although all samples contain older (15–35 Ma) single-grain ages, only the Santo Domingo, Agua Viva, Tucaní, and Chama catchments record an older (early–middle Miocene; 14–25 Ma) age peak. The late Miocene age population is dominant (57%–100%) in most samples, except for the Santo Domingo and Coloncito catchments, where the younger (Miocene–Pliocene) age population is larger. In the Tucaní sample, 90% of the grains are grouped in a 5.9 ± 0.7 Ma age peak. The Chejendé River is a tributary of the Agua Viva and drains part of its catchment (Fig. 2); for this reason, we merged the detrital data for the corresponding samples (0507 and 0807) and applied the decomposition of the grain-age distribution as done previously. The resulting age peaks are very similar to those from the Agua Viva sample alone (Table 1).

Access to the southern flank of the Venezuelan Andes is difficult, and our data set contains only one river draining to the south. The Santo Domingo River sample shows a similar grain-age distribution and peak ages as the rivers that drain to the north (Table 1; Fig. 2), although the component ages are slightly older than corresponding age peaks in the other samples.

Although binomial peak fitting provides a convenient means of analyzing detrital thermochronology data, it may mask similarities in age structure between different samples. In order to analyze this, we use the Kolmogorov-Smirnov and Kuiper tests (Conover, 1980) to quantitatively compare the detrital age probability density functions (PDF) for the seven catchments (cf. supplementary material, Figure S1 and Table S11). These comparisons indicate that, at a 95% confidence level, the samples from four catchments located in the central part of the Venezuelan Andes (Chama, Tucaní, Santo Domingo, and San Pedro) yield a similar detrital age distribution, while the age distributions of the other samples are significantly different from each other.

Comparison of Detrital and Bedrock AFT Ages

Previous studies have shown that detrital thermochronology faithfully records source-area exhumation in a number of different settings (e.g., Bernet et al., 2004; Ruhl and Hodges, 2005; Brewer et al., 2006). Here, we compare our detrital AFT grain-age distributions of an individual drainage to previously published bedrock AFT ages. We analyze the Chama catchment, because it is the only one within the Venezuelan Andes for which sufficient in situ data (Kohn et al., 1984; Bermúdez et al., 2011) are currently available to derive a reliable bedrock age map.

Several authors (Stock et al., 2006; Vermeesch, 2007; McPhillips and Brandon, 2010) have developed quantitative methods for comparing bedrock and detrital age patterns. We adopted the approach described by Vermeesch (2007), randomly sampling a number of ages equal to the number of single-grain ages in the detrital sample from the interpolated bedrock age map. This method has the advantage that it takes into account Poisson-distributed measurement uncertainties. We sampled the bedrock age distribution and constructed the predicted age probability-density function (PDF) and cumulative-density function (CDF) 1000 times (Fig. 3A). We then used the Kolmogorov-Smirnov (K-S) (Conover, 1980) and Kuiper equality tests (Ruhl and Hodges, 2005) to compare the predicted AFT CDFs and PDFs with the detrital AFT data. Additional simulations (5000; 10,000) do not change the statistical results of our comparison. The resulting average p-values show that the inferred bedrock age distribution is significantly different from the detrital age distribution at the 95% confidence level; all simulations fail the Kuiper test, and the vast majority also fail the K-S test.

There are several possible reasons for this discrepancy. First, we do not take into account potential variations in apatite content for different lithologies in the catchment, whereas this may strongly influence the detrital age distribution (e.g., Tranel et al., 2011). However, (1) we have no information on the relative apatite abundance in different rock types of the Venezuelan Andes, and (2) the sampled Chama catchment is mainly underlain by relatively homogeneous Precambrian basement rocks (see Chama River in Fig. 1A).

Another possibility is that the catchment is not currently eroding at a uniform rate, which is an inherent assumption when randomly picking ages from the interpolated bedrock age distribution. It is relatively simple to weight the bedrock age distribution according to a model of inferred present-day erosion rates (Vermeesch, 2007; McPhillips and Brandon, 2010). Here, we assume that the short-term erosion rates in the catchment are coupled to the long-term exhumation rates, which can be estimated from the AFT ages using a simple one-dimensional (1-D) thermal model (cf. next section). We thus resampled the interpolated bedrock age distribution, letting the probability of sampling any particular age be determined by the long-term exhumation rate associated with that age, relative to the average long-term exhumation rate of the catchment. The resulting predicted CDFs/PDFs compare very favorably to the observed detrital AFT age distribution (Fig. 3B); both the K-S and Kuiper tests suggest the two distributions are similar at the 95% confidence level.

We conclude from these comparisons that the detrital AFT age distributions provide a reliable estimate of in situ bedrock ages in the sampled catchments when weighted by their associated long-term exhumation rates, implying that the apatite yield of the catchment (as recorded by our detrital samples) is controlled by the long-term exhumation rates.

Long-Term Exhumation Rates

The detrital age distributions can be used to provide predictions of average long-term exhumation rates of individual drainages, which can be compared to potential tectonic or climatic control parameters. We employed a simple 1-D steady-state thermal model developed by Brandon et al. (1998; see also Ehlers, 2005; Reiners and Brandon, 2006) to convert detrital AFT ages to exhumation rates. Because this model does not take either transient thermal effects (e.g., Rahl et al., 2007) or the three-dimensional (3-D) effects of topography (e.g., Whipp et al., 2009) into account, predicted exhumation rates should be considered as first-order estimates. Parameters used for these calculations are a surface temperature of 25 °C, layer thickness to constant temperature of 40 km, thermal diffusivity of 25 km2 m.y.−1, heat production of 10 °C m.y.−1, and a temperature at the base of the layer of 700 °C (for details, see Bermúdez et al., 2011). The last two values were constrained by inverse modeling of age-elevation relationships in the central Venezuelan Andes (Bermúdez et al., 2011) and imply a pre-exhumation surface geothermal gradient close to 25 °C km−1. The model iteratively calculates an exhumation rate (ɛ), cooling rate, closure temperature (Tc), and closure depth (zc) from the AFT ages, using the Dodson (1973) equation to solve for the closure temperature of the AFT thermochronometer and a steady-state thermal structure to relate this to the closure depth.

For each detrital sample, we used this procedure to translate single-grain ages into exhumation rates, which we then combined to predict a catchment-wide average long-term exhumation (ɛT) rate as:
where ɛj is the exhumation rate inferred for grain j, and N is the total number of grains in the detrital sample. In order to include the expected effect that rapidly exhuming areas of a catchment will contribute more sediment (and thus more datable apatite grains), as indicated by the comparison of bedrock and detrital age distributions in the previous section, we also calculated a weighted long-term catchment-averaged exhumation rate (ɛT)w:


Results are reported in Table 2. Catchment-averaged exhumation rates vary between 0.48 ± 0.02 km m.y.−1 for the Santo Domingo catchment and 0.80 ± 0.26 km m.y.−1 for the Mimbós catchment. As expected, weighted-average rates are lower, and all weighted long-term exhumation rates in the central Venezuelan Andes overlap within error, varying between 0.33 ± 0.07 km m.y.−1 (Agua Viva–Chejendé) and 0.48 ± 0.08 (Coloncito) km m.y.−1 The northeastern Agua Viva–Chejendé catchment shows the lowest long-term exhumation rates (<0.35 km m.y.−1). The central Chama and Santo Domingo catchments show intermediate rates (0.35–0.40 km m.y.−1), and the catchments located on the northwestern flank of the range show the highest exhumation rates (>0.40 km m.y.−1). In the following discussion, we will focus on the weighted-average catchment exhumation rates, which we consider to be the most reliable estimates.

In order to explore the potential topographic, tectonic, or climatic controls on the exhumation rates determined herein, we derive quantitative measures for these controls (relief, precipitation, seismicity, and stream power) in the following sections.

Elevation and Relief

The Venezuelan Andes rise from sea level at Lake Maracaibo, to the north of the mountain belt, to close to 5000 m elevation at Pico Bolívar in the Sierra Nevada block (Figs. 1 and 4A). This part of the Andes most probably emerged as an orographic barrier at ca. 8 Ma (Hoorn et al., 1995; Bermúdez et al., 2011). The central Sierra Nevada and Sierra la Culata ranges (Fig. 1B) include more than 60 peaks with elevations of 4300 m and above (Fig. 4A). These highest parts of the Venezuelan Andes were glaciated during the Quaternary. Moraines of the Last Glacial Maximum are encountered at 3400–3600 m elevation (Schubert, 1984), with older moraines occurring at elevations several hundred meters lower. The mean elevation of the range is significantly lower to the northeast and southwest, varying between 1160 and 1695 m (Bermúdez et al., 2010).

Relief was calculated as the maximum elevation difference within a variable radius, between 1 and 15 km, for every pixel in the digital elevation model (e.g., Montgomery and Brandon, 2002). There is a power-law relationship between average relief of the different catchments and the radius at which it was calculated, with scaling exponents around 0.5–0.6, characteristic for low-latitude mountain belts (Champagnac et al., 2012). The present-day relief of the Venezuelan Andes for a 5 km radius is shown in Figure 4A. Highest relief values of close to 3500 m are found along the central Chama River valley, which follows the Boconó fault. Other areas of high relief are located on the northwestern and southeastern flanks of the mountain belt. Relief is much lower within the central Sierra Nevada and Sierra la Culata ranges and in the Caparo and Trujillo blocks in the southwest and northeast of the Venezuelan Andes (Figs. 1 and 4A). The former are characterized by extensive low-relief summit surfaces that show widespread evidence for glacial beveling through cirque retreat (e.g., Mitchell and Montgomery, 2006; Foster et al., 2008), whereas the latter show subdued relief at lower mean elevations with widespread soil mantling of slopes.


The topography of the Venezuelan Andes generates a well-developed orographic precipitation pattern. Moist air masses from equatorial South America in the south need to rise in order to cross the Venezuelan Andes to the north, causing heavy precipitation on its southern flank, whereas the center and northern flank of the belt are more arid. It has been proposed that the present-day precipitation and drainage system of the Venezuelan Andes started developing at ca. 8 Ma, in response to tectonically driven surface uplift (Hoorn et al., 1995; Bermúdez et al., 2011). We used precipitation data from 30 meteorological stations in and around the Venezuelan Andes to elaborate a mean annual precipitation map for the past 20 yr. Although these measurements do not record snowfall as accurately as rainfall, the vast majority of precipitation is rainfall in the low-latitude Venezuelan Andes. The database was compiled from different sources: (1) National Oceanic and Atmospheric Administration (NOAA) World Temperature- Precipitation data set (2011); (2) published data (Stansell et al., 2006; Naranjo and Duque, 2004); (3) bioclimatic net stations of Mérida (http://www.cecalc.ula.ve/redbc/colecciones/colecciones_datos.html); and (4) data provided by Gerard Kopp (2007, personal commun.) of the Institute for Meteorology and Climate Research, University of Karlsruhe (Germany) for the Mérida Atmospheric Research Station at Pico Espejo (MARS). The present-day precipitation map shown in Figure 4B was compiled from these data using nearest-neighbor interpolation (Arya et al., 1998), and it highlights the strong precipitation gradient between the northern and southern flanks of the Venezuelan Andes, with precipitation ranging from 1.6 to 3.4 m yr−1 on the southern flank as opposed to 0.09–1.1 m yr−1 on the northern flank. A particular feature is the arid patch that expands into the orogen between El Vigía and Mérida: Hot dry winds from the north exploit the Chama River valley to penetrate deep into the orogen, generating a unique arid climate in this area (Fig. 4B). A remotely sensed (Tropical Rainfall Measuring Mission [TRMM], NASA) precipitation data set synthesized by Bookhagen and Strecker (2008) shows overall similar patterns, both qualitatively and quantitatively.


Historical and instrumental seismicity in the Venezuelan Andes is strongly concentrated along the Boconó fault system. Smaller seismic events are scattered within a band of several tens of kilometers width adjacent to this fault zone, indicating that many of the fault branches are also active (Fig. 5A). However, most of the seismicity occurs along the main trace of the fault at an average depth of ∼15 km. Larger seismic events tend to occur at greater depth to the northwest (Lake Maracaibo basin) and southeast (Barinas basin) of the surface trace of the Boconó fault, reaching more than 40 km depth (Dewey, 1972; Niu et al., 2007). A seismic zone of intermediate depth (∼160 km) occurs toward the southwestern boundary of the Venezuelan Andes and continues below the Eastern Cordillera of Colombia. This significant concentration of events is known as the Bucaramanga seismic nest (Schneider et al., 1987) and occurs within the subducted Caribbean slab below northern South America (van der Hilst and Mann, 1994). More scattered, but relatively high-magnitude (M > 5) earthquakes are associated with the north-south–oriented Icotea, Valera, and Burbusay fault systems to the north of the Venezuelan Andes (Fig. 5A).

Focal mechanisms suggest predominantly right-lateral faulting along the Boconó fault, left-lateral faulting along the north-south–trending faults to the north of the Venezuelan Andes, and orthogonal thrusting along the northwest and southeast foreland fold-and-thrust belts (Colmenares and Zoback, 2003; Corredor, 2003; Cortés and Angelier, 2005). These focal mechanisms indicate a compressional stress regime with σ1 oriented approximately WNW-ESE across the central Venezuelan Andes, evolving toward a strike-slip regime with a NW-SE–directed σ1 axis to the northeast. Studies of active tectonic landforms in the Venezuelan Andes, together with regional tectonic reconstructions, suggest that the current tectonic regime was installed during Pliocene–Quaternary times (Backé et al., 2006; Egbue and Kellogg, 2010).

In order to quantify the effects of seismicity across the Venezuelan Andes, we compiled a seismicity record over the past century (from January 1911 to January 2011), using data available from the digital library of the Geophysics Laboratory of Universidad de Los Andes, Mérida (http://lgula.ciens.ula.ve/), and data provided by the Fundación Venezolana de Insvestigaciones Sismológicas (FUNVISIS; http://www.funvisis.gob.ve). The data were filtered using two criteria: (1) epicenters of earthquakes located between latitudes 72.25°W and 70.00°W and longitudes 7.75°N and 10.00°N; (2) only data with reported local magnitudes Ml were considered. Figure 5A shows a summary of earthquakes that occurred in the Venezuelan Andes during the past century. We calculated the released seismic energy (Se) from the local magnitudes using the classical expression of Gutenberg and Richter (1954):


The values of the parameters a and b were estimated by a least-squares fit of cumulative magnitude-frequency relationships constructed from subsamples of the seismic database for each catchment in 0.5° × 0.5° cells and are equivalent to the intersect and slope, respectively, of the Gutenberg-Richter relationship (Gutenberg and Richter, 1954). Released seismic energy values were summed within circles with radius of 25 km around the epicenter of each earthquake. Figure 5B shows the resulting map of seismic energy obtained from this procedure. The rate of seismic energy release is highest in the central Venezuelan Andes (Chama catchment) and decreases toward the northeast and southwest. However, the calculated pattern is strongly affected by the Ml > 5 earthquakes that occurred on the Burbusay, Valera, and Icotea fault systems.

The record of seismic energy release spans only 100 yr, which is significantly shorter than the return period of major earthquakes in the Venezuelan Andes (300 yr for M > 7 earthquakes; Audemard, 1997). Consequently, the pattern of released seismic energy may be locally underestimated. However, the seismicity records brittle deformation of the upper crust; the spatial distribution and frequency of earthquakes are intuitively related to the rate of brittle deformation (e.g., Holt et al., 2000). Thus, we can use the compiled seismic database to estimate the present-day distribution of brittle strain rate and extrapolate the total amount of seismic strain over time scales longer than the observation interval, using the observed earthquake magnitude-frequency (Gutenberg-Richter) relationship. To achieve this, we used the method described by Braun et al. (2009) and calculated seismic/brittle strain rate as:
in which the parameters a and b are defined and calculated as explained already; Mmax is the maximum observed magnitude; μ is elastic shear modulus; and ΔV is the volume of the crust (that is, the moving 0.5° × 0.5° cell area multiplied by depth of the maximum-magnitude earthquake) in which the earthquakes were observed over a period of time Δt (in this case, Δt = 100 yr). The depth of the maximum-magnitude earthquakes (33–163 km) generally exceeds the depth of the brittle-ductile transition (∼15–20 km). Because we used the seismic strain rate as a proxy for brittle deformation, we recalculated seismic strain rates using only earthquakes with hypocentral depths less than 20 km. However, we did not find any significant differences between the two approaches. The maximum-magnitude bin used in the least-squares fit always had at least one earthquake in it. Cells with a correlation coefficient r2 < 0.95 (mostly due to an incomplete catalogue) were not considered. In order to check the robustness of the results, we tested the effect of removing the minimum- and maximum-magnitude bins on the a and b values and the effect of forcing b to be exactly 1, but neither caused the seismic strain rate to change significantly.

Figure 5C shows the resulting map of seismic/brittle strain rate from subsamples (sliding 0.5° × 0.5° cells) of the seismic database for the Venezuelan Andes. The resulting seismic strain-rate map (Fig. 5C) shows maximum deformation occurring in the central and northwest Venezuelan Andes, with a second maximum occurring to the south of the mountain belt in the area of the 21 December 2001, Mw = 5.6 earthquake. Seismic strain rates vary between ∼10−17 s−1 at the northeastern (Agua Viva–Chejendé catchment) and southwestern extremities of the belt, the latter influencing average seismic strain rates in the Chama catchment, and ∼10−15 s−1 in the center of the belt (Mimbós, San Pedro, and Coloncito catchments).

Short-Term Erosion Potential

The long-term exhumation rates estimated from detrital AFT ages can be compared to model predictions of the intensity of short-term erosion using the erosion-index approach of Finlayson et al. (2002). The erosion index (EI) can be calculated in different ways as a function of stream power, which is the rate of potential energy expenditure by flowing water and has been used extensively in studies of erosion, sediment transport, and geomorphology as a measure of the erosive power of rivers and streams (Wilson and Gallant, 2000; Wobus et al., 2006). The analysis is based on a prediction of bedrock incision rate as a function of stream power (Finlayson et al., 2002; Tucker and Whipple, 2002):
where e is the local incision rate, A is upstream drainage area (used as a proxy of discharge), S is local slope, and m, n, and k are constants. The parameter k is mainly related to bedrock erodibility. Given the similar geologic and geographic conditions of the studied catchments and the lack of information on relative erodibility of the different lithologies outcropping in the Venezuelan Andes, we take k to be uniform throughout the study area. Setting k to unity, the predicted erosion index becomes:


Spatial variations in precipitation P can be incorporated into the prediction of EI in order to study their influence on spatial variability of the erosion potential:
where (Ap) is the pixel area, P is the local precipitation, and the summation sign implies summing along the flow-lines within the catchment in order to calculate the flow accumulation. We term EIp the precipitation-modulated erosion index.

Different m and n values can be employed in Equations 6 and 7, depending on the control of river incision rates by total stream power, stream power per unit channel width, or shear stress (Finlayson et al., 2002; Tucker and Whipple, 2002). In the case where incision is controlled by total stream power (TSP), m = n = 1. For incision controlled by stream power per unit channel width (USP), m = 1/2 and n = 1. If incision is controlled by fluvial shear stress power (SSP), m = 1/3 and n = 2/3.

Logarithmic erosion index maps for the Venezuelan Andes predicted by the six possible models (TSP, USP, and SSP for uniform or spatially variable precipitation) are shown in Figure 6. These maps emphasize distinct zones of high erosion potential associated with the steepest terrain and largest discharge of the main rivers. If uniform precipitation is considered, the most important zones of high erosion potential of the Venezuelan Andes are located in the central Chama and Santo Domingo valleys, with secondary regions of high erosion potential on both the northern and southern flanks. Consequently, under this assumption, the Chama and Santo Domingo Rivers, together with the small catchments (Tucaní, San Pedro, and Mimbós) on the northwest flank of the orogen, have high relative erosion indexes (whether considering TSP, USP, or SSP), with the peripheral Coloncito and Agua Viva–Chejendé River catchments showing lower erosion indexes (Table 2).

When the precipitation data are incorporated in the erosion index predictions, the results differ significantly because of the strong orographic effect (Fig. 4B), with maximum erosion potential shifting to the southeastern flank of the mountain belt (Fig. 6). In particular, the predicted high relative erosion index in the central Chama River valley disappears when precipitation is taken into account, due to the pronounced aridity of this valley. Consequently, the southeast-draining Santo Domingo catchment is predicted to have the highest relative erosion index (whether considering TSP, USP, or SSP), followed by the Chama River catchment. Incorporating spatially variable precipitation lowers the relative erosion index values for the small northwestern Tucaní, San Pedro, and Mimbós catchments, which become indistinguishable from those for the peripheral Agua Viva–Chejendé and Coloncito catchments (Table 2).

Whether long-term exhumation rates in mountain belts are controlled to first order by tectonic or climatic factors remains a matter of debate. Some authors have argued for strong coupling between long-term exhumation and precipitation, for instance, for the Central Andes (Montgomery et al., 2001) and Washington Cascades (Reiners et al., 2003), while others have argued for decoupling between long-term exhumation and precipitation, such as in the central Himalaya (Burbank et al., 2003) or in the European Alps (Vernon et al., 2009). Recently, Mora et al. (2008) argued for coupling between precipitation and exhumation rates in the Eastern Cordillera of Colombia, which is contiguous to the Venezuelan Andes. Can long-term exhumation rates in the Venezuelan Andes be correlated to present-day precipitation rates? If they are uncorrelated with precipitation, what could be the controlling factor driving exhumation? The underlying question is, how strongly are tectonic and climatic processes coupled and on what time scales? Table 2 summarizes our observations on long-term exhumation, predicted short-term erosion potential, and potential controlling factors for the seven studied catchments in the Venezuelan Andes. The challenge lies in comparing exhumation rates over several millions of years, derived from detrital AFT ages, with present-day precipitation or seismicity data that were collected over the past 100 yr or less.

In order to gain insight into the potential controls on erosion and exhumation rates in the Venezuelan Andes, we calculated Pearsonian correlation coefficients between each of the measures reported in Table 2. Table 3 summarizes the results of this analysis. Given that we compare all variables for seven catchments, correlations are statistically significant (at a 95% confidence level) for a Pearson correlation coefficient r ≥ 0.7.

Figure 7 shows how long-term exhumation rates in the Venezuelan Andes correlate with different potential control parameters (present-day relief, precipitation, seismic energy release, seismic strain rate, and stream power). A first observation is the very strong positive correlation (r = 0.92) between long-term exhumation rates and 5-km radius relief, implying that the present-day relief of the Venezuelan Andes is adapted to the long-term exhumation rates. The correlation we find here is much stronger than that reported, for instance, in the European Alps (Wittmann et al., 2007; Vernon et al., 2009) or in the San Bernardino Mountains of California (Binnie et al., 2007). The reason for this is probably the comparably lower exhumation rates and relief in the Venezuelan Andes, so that bedrock landsliding on threshold slopes, which tends to decouple relief and exhumation rates (Montgomery and Brandon, 2002), is less important with respect to these other mountain belts. In effect, large proportions of slopes close to the threshold for landsliding are only encountered in the central Sierra Nevada and Sierra la Culata blocks of the Venezuelan Andes (Bermúdez et al., 2010).

In contrast, long-term exhumation rates appear to be decoupled from present-day precipitation rates; they appear even (weakly) negatively correlated with present-day precipitation (i.e., highest exhumation rates occur on the dry northwest flank of the orogen). Also, whereas short-term erosion indexes (EI) are not correlated to long-term exhumation rates for any of the three models (i.e., TSP, SSP, or USP), the correlations become even weaker or negative when spatially variable precipitation rates are taken into account in the calculation of these measures (EIp). The lack of correlations suggests that the present-day precipitation pattern does not strongly control long-term exhumation rates in the Venezuelan Andes, in contrast to what has been inferred for the contiguous Eastern Cordillera of Colombia (Mora et al., 2008). However, two caveats need to be recalled before concluding from these data that climate and exhumation are decoupled in the Venezuelan Andes: (1) the relatively limited spatial extent of our sampled catchments, for which the wet southern flank is clearly under-represented, and (2) the time scale problem mentioned previously. Considering the first, it is clear that our data set would benefit from inclusion of additional remote south-flank catchments. As for the problem of comparing measurements on strongly varying time scales, whereas we obviously acknowledge that we cannot readily extrapolate the present-day precipitation measurements several million years into the past, we do note that (1) the Neogene sedimentary record of northern South America suggests that the topographic relief of the Venezuelan Andes, and therefore its associated orographic imprint, has been in place since late Miocene times (Hoorn et al., 1995; Díaz de Gamero, 1996); and (2) Quaternary glaciations, which may strongly modify erosion patterns in some mountain belts (e.g., Burbank et al., 2003; Gabet et al., 2008; Vernon et al., 2009), have been relatively limited in the Venezuelan Andes (Schubert, 1984; Stansell et al., 2006). It has been suggested, both in the Andes and elsewhere (e.g., Bookhagen et al., 2005; Abbühl et al., 2010), that spatial patterns of erosion are strongly influenced by the occurrence of extreme events, leading to decoupling between patterns of erosion and present-day average climate. This could be the case in the Venezuelan Andes, but sufficient resolution in the regional climate record to address this question is currently lacking.

In order to assess potential tectonic controls on long-term exhumation rates, we used cumulative seismic energy release and seismic strain rate as proxies for tectonic forcing (Tables 2 and 3). We estimated the cumulative seismic energy released within the different catchments in two different ways: by averaging the values of Figure 5B within each catchment, or by simply summing the seismic energy released by each earthquake within the catchment boundaries and normalizing by catchment area. The difference between the two methods is the extent to which earthquakes occurring outside the catchment boundaries are taken into account.

The inferred seismic energy release correlates weakly to moderately with either catchment-averaged exhumation rate or 5-km-radius relief, with r values varying between 0.41 and 0.64 (Table 3). As noted previously, the pattern of seismic energy release is strongly controlled by large (Ml > 5) strike-slip earthquakes located outside the Venezuelan Andes, which contribute significantly to accumulated seismic energy in the low-exhumation Agua Viva–Chejendé catchment, for instance. The seismic energy release only integrates a century of data and, similar to the precipitation records, may therefore not be directly comparable to long-term exhumation rates and relief. However, there is very little correlation either between measures of seismic energy release and stream-power–based erosion indexes, suggesting that short-term erosion rates (if these are correctly predicted by the stream-power models) are decoupled from seismic energy release.

In contrast to the pattern of released seismic energy, the inferred seismic strain rate correlates very strongly with both catchment-averaged exhumation rates (r = 0.81) and 15-km-radius relief (r = 0.85). Catchments on the northwest flank of the orogen that show the largest seismic strain rates (Mimbós, San Pedro, Coloncito; ∼10−15 s−1) are also characterized by the highest relief (>1650 m) and exhumation rates (>0.4 km m.y.−1). In contrast, the northwestern Agua Viva–Chejendé catchment, which has the lowest relief and exhumation rates, also shows very low seismic strain rates (2.5 × 10−17 s−1).

Our data thus suggest that active tectonic deformation of the Venezuelan Andes provides a much stronger control on long-term exhumation rates throughout the orogen than the pattern weighted by precipitations. They also suggest that the widely used stream-power–based erosion potential does not adequately model erosion rates throughout the mountain belt. The reason for this may be that the simple stream-power expression does not take into account the effects of sediment flux, stochasticity of discharge, and incision thresholds (e.g., Whipple, 2004; Lague, 2010). We have previously suggested an overall tectonic control on relief and exhumation rates in the Venezuelan Andes because the orogen can be divided into fault-bounded blocks that each have their characteristic relief and exhumation history (Bermúdez et al., 2010), as also observed in other obliquely convergent mountain belts (e.g., Spotila et al., 2007). It would be interesting to map out erosion rates throughout the Venezuelan Andes using cosmogenic nuclides in stream sediments in order to analyze whether the same patterns emerge on shorter time scales.

Detrital AFT analysis of modern river sediment provides an efficient tool for studying long-term exhumation rates in complex areas with difficult access, such as the Venezuelan Andes. We have shown that the technique efficiently records average catchment-wide long-term exhumation rates when compared to in situ samples. The detrital AFT ages of seven major catchments of the Venezuelan Andes reproduce and extrapolate the exhumation patterns reported by Bermúdez et al. (2010, 2011) and can thus be used to infer exhumation rates in areas where in situ data are absent. Our results suggest that using weighted average exhumation rates, in order to include the effect of higher sediment yields in rapidly exhuming areas of a catchment, improves both the fit to observed detrital age distributions and the correlations with potential control parameters.

Catchment-averaged long-term exhumation rates in the Venezuelan Andes are fairly constant at 0.31–0.48 km m.y.−1 Rates are higher (>0.4 km m.y.−1) on the northwestern flank of the mountain belt (Mimbós, San Pedro, Tucani, and Coloncito Rivers), whereas they are lowest (0.31 km m.y.−1) in the northeastern Agua Viva–Chejendé catchment. The large Chama and Santo Domingo catchments draining the central Venezuelan Andes show intermediate rates of 0.37–0.38 km m.y.−1 Average exhumation rates correlate very well with relief of the catchments as well as with seismic strain rates, but correlations with both stream-power–based predictions of short-term erosion potential and precipitation are weak or negative.

The lack of correlation between long-term exhumation rates and relief, on the one hand, and precipitation, on the other, suggests that either precipitation is not the main controlling factor driving exhumation and relief development in the Venezuelan Andes, or that the short-term record of present-day precipitation is not representative of long-term average values. Alternatively, erosion may be effectively controlled by extreme climate events that are not captured by the record. A moderate correlation with seismic energy release indicates that the 100 yr seismic record does not fully capture the current deformation, because of the long recurrence time of large seismic events. However, we show that calculating the seismic strain rate and its extrapolation over geologic time is feasible and allows us to study relationships among long-term exhumation rate, relief, and tectonics. Our study suggests that the climatic control on exhumation rates in the active orogens of northern South America may have been overstated. Climatic variations appear to only have a second-order control on relief, erosion, and long-term exhumation rates, as has been suggested elsewhere (e.g., Riebe et al., 2001a, 2001b; Kirchner et al., 2001).

This study was supported by the Consejo de Desarrollo Científico y Humanístico (CDCH) de la Universidad Central de Venezuela (UCV), project number PI 08-00-6219-2006 and ECOS-Nord project V08U01. We thank engineers Juan Flores and Antonio León for support during field work. We acknowledge Gerhard Kopp of the Institute for Meteorology and Climate Research, University of Karlsruhe, Germany, for providing the meteorological data from the Mérida Atmospheric Research Station for Pico Espejo (MARS). María Elena Naranjo of the Universidad de Los Andes (ULA) provided precipitation data of the following stations: Páramo de Mucuchies, Páramo Pico El Aguila, Los Plantíos, Mucubají, Valle Grande, and Tabay. FUNVISIS (Fundación Venezolana de Investigaciones Sismológicas) provided access to the seismicity data. Finally, we thank Jean Braun and Christoph Glotzbach for discussion and assistance during the writing of this manuscript, and Bodo Bookhagen and an anonymous reviewer for constructive comments.

1GSA Data Repository Item 2012271, Fig S1: comparison among cumulative density functions (CDF) of single-grain ages for all the catchments; Table S1: Results of Kolmogorov-Smirnov equality test for the CDF’s, is available at www.geosociety.org/pubs/ft2012.htm, or on request from editing@geosociety.org, Documents Secretary, GSA, P.O. Box 9140, Boulder, CO 80301-9140, USA.