Mauricio A Bermúdez
email: mauricio.bermudez@ing.ucv.ve

Peter A Van Der Beek

Matthias Bernet

Strong tectonic and weak climatic control on exhumation rates in the Venezuelan Andes

Keywords: Venezuelan Andes, detrital apatite fission-track, exhumation, erosion rates, relief, climate

We have studied the relationships between 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 to 0.80±0.26 km Myr -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 to 0.48±0.08 km Myr -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, but no 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 timescales 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.

INTRODUCTION

The relative importance of tectonic and climatic control on relief development and exhumation in mountain belts remains strongly debated (e.g., [START_REF] Molnar | Late Cenozoic uplift of mountain ranges and global climate change: chicken or egg?[END_REF][START_REF] Burbank | Decoupling of erosion and precipitation in the Himalayas[END_REF][START_REF] Reiners | Coupled spatial variations in precipitation and long-term erosion rates across the Washington Cascades[END_REF][START_REF] Strecker | Does the topographic distribution of the central Andean Puna Plateau result from climatic or geodynamic processes?[END_REF][START_REF] Champagnac | Tectonics, climate, and mountain topography[END_REF].

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., [START_REF] Whipple | The influence of climate on the tectonic evolution of mountain belts[END_REF], and the different temporal scales for which we have records of these processes [START_REF] Burbank | Decoupling of erosion and precipitation in the Himalayas[END_REF][START_REF] Reiners | Coupled spatial variations in precipitation and long-term erosion rates across the Washington Cascades[END_REF][START_REF] Vernon | Spatial correlation between longterm exhumation rates and present-day forcing parameters in the western European Alps[END_REF]. While thermochronologic data record denudation on million-year timescales, 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 [START_REF] Finlayson | Spatial coincidence of rapid inferred erosion with toung metamorphic massifs in the Himalayas[END_REF][START_REF] Dadson | Links between erosion, runoff variability and seismicity in the Taiwan orogen[END_REF][START_REF] Koons | On the implications of low spatial correlation of tectonic and climate variables in the western European Alps[END_REF][START_REF] Vernon | Spatial correlation between longterm exhumation rates and present-day forcing parameters in the western European Alps[END_REF].

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., [START_REF] Montgomery | Climate, tectonics, and the morphology of the Andes[END_REF]Lamb and Davies, 2003), as well as the strong asymmetry in orographic precipitation associated with this orogen (e.g., [START_REF] Bookhagen | Orographic barriers, high-resolution TRMM precipitation, and relief variations along the eastern Andes[END_REF]. While many studies have concentrated on the southern and central Andes (e.g., [START_REF] Strecker | Tectonics and climate of the southern central Andes[END_REF][START_REF] Farías | Late Miocene high and rapid surface uplift and its erosional response in the Andes of central Chile[END_REF][START_REF] Thomson | Glaciation as a destructive and constructive control on mountain building[END_REF], attention has recently also focused on the northern Andes. [START_REF] Mora | Climatic forcing of asymmetric orogenic evolution in the Eastern Cordillera of Colombia[END_REF], for instance, describe strongly asymmetric deformation and exhumation patterns in the eastern Colombian Andes, which they link to the strong NW-SE precipitation gradient across this part of the mountain belt. In contrast, [START_REF] Mora | The role of inherited extensional fault segmentation and linkage in contractional orogenesis: a reconstruction of Lower Cretaceous inverted rift basins in the Eastern Cordillera of Colombia[END_REF] and [START_REF] Parra | Episodic orogenic front migration in the northern Andes: Constraints from low-temperature thermochronology in the Eastern Cordillera[END_REF] argue 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 eight river catchments in the Venezuelan Andes, the north-easternmost extension of the Northern Andes.

We compare the detrital AFT data to available bedrock AFT data of the Venezuelan Andes [START_REF] Kohn | Mesozoic-Pleistocene fission track ages on rocks of the Venezuelan Andes and their tectonic implications[END_REF][START_REF] Bermúdez | Spatial and temporal patterns of exhumation across the Venezuelan Andes: Implications for Cenozoic Caribbean geodynamics: Tectonics[END_REF][START_REF] Bermúdez | Asynchronous Miocene-Pliocene exhumation of the central Venezuelan Andes[END_REF] in order to test whether they faithfully record exhumation in the catchments (e.g., [START_REF] Ruhl | The use of detrital mineral cooling ages to evaluate steady state assumptions in active orogens: An example from the central Nepalese Himalaya: Tectonics[END_REF][START_REF] Brewer | Downstream development of a detrital cooling-age signal[END_REF] and to spatially integrate the still relatively sparse in-situ data. The objective of this study is to examine the relative control 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.

GEOLOGIC SETTING

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; [START_REF] Case | Geology of the Northern Andes; an overview: The Geology of North America[END_REF]. The orogen is characterized by high seismicity and spatially variable exhumation since Miocene times, as recorded by AFT thermochronology data [START_REF] Kohn | Mesozoic-Pleistocene fission track ages on rocks of the Venezuelan Andes and their tectonic implications[END_REF][START_REF] Bermúdez | Spatial and temporal patterns of exhumation across the Venezuelan Andes: Implications for Cenozoic Caribbean geodynamics: Tectonics[END_REF][START_REF] Bermúdez | Asynchronous Miocene-Pliocene exhumation of the central Venezuelan Andes[END_REF]. The structure of the Venezuelan Andes is controlled by reactivated faults that delineate individual tectonic blocks (Fig. 1B). These faults were inherited from Early Mesozoic rifting, and possibly from earlier orogenesis [START_REF] Aleman | Northern Andes[END_REF][START_REF] Pindell | Kinematic evolution of the Gulf of Mexico and Caribbean[END_REF]. [START_REF] Bermúdez | Spatial and temporal patterns of exhumation across the Venezuelan Andes: Implications for Cenozoic Caribbean geodynamics: Tectonics[END_REF] 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 (Coletta et al., 1997). 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. 1B). The Boconó fault zone extends from the border with Colombia for more than 500 kilometers to the northeast (Fig. 1) 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. These two blocks cooled rapidly but diachronously during the late Miocene-Pliocene [START_REF] Kohn | Mesozoic-Pleistocene fission track ages on rocks of the Venezuelan Andes and their tectonic implications[END_REF][START_REF] Bermúdez | Asynchronous Miocene-Pliocene exhumation of the central Venezuelan Andes[END_REF]. In contrast, the Caparo and Trujillo blocks (Fig. 1B), at the southwestern and northeastern ends of the Venezuelan Andes respectively, experienced slow cooling from the late Oligocene -late Miocene onward [START_REF] Bermúdez | Spatial and temporal patterns of exhumation across the Venezuelan Andes: Implications for Cenozoic Caribbean geodynamics: Tectonics[END_REF]. 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).

DETRITAL APATITE FISSION-TRACK THERMOCHRONOLOGY

Thermochronology of detrital minerals allows quantifying cooling rates and exhumation processes in convergent mountain belts (e.g., [START_REF] Garver | Exhumation history of orogenic highlands determined by detrital fission track thermochronology[END_REF][START_REF] Carter | Heavy minerals and detrital fission-track thermochronology[END_REF], providing a complimentary record of the erosional history of a mountain belt with respect to local cooling paths deduced from in-situ samples within the orogen [START_REF] Brandon | Tectonic evolution of the Cenozoic Olympic subduction complex, Washington State, as deduced from fission track ages for detrital zircons[END_REF]. 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 characterise exhumation patterns for inaccessible areas.

Data collection and discrimination of age components

We collected samples for detrital AFT thermochronology from eight 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. [START_REF] Bermúdez | Spatial and temporal patterns of exhumation across the Venezuelan Andes: Implications for Cenozoic Caribbean geodynamics: Tectonics[END_REF]. We aimed at dating at least 100 grains per sample in order to attain statistically significant AFT age populations. This objective was attained for 5 out of 8 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 [START_REF] Brandon | Tectonic evolution of the Cenozoic Olympic subduction complex, Washington State, as deduced from fission track ages for detrital zircons[END_REF][START_REF] Garver | Exhumation history of orogenic highlands determined by detrital fission track thermochronology[END_REF]. Different methods have been proposed for decomposing 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;[START_REF] Brandon | Decomposition of fission-track grain-age distributions[END_REF][START_REF] Brandon | Probability density plots for fission-track grain age distributions[END_REF]. Here, we use this method (as described by [START_REF] Stewart | Implications for late Cenozoic evolution of the Cascadia subduction wedge[END_REF] to decompose the grain-age distributions of our modern river samples.

Resulting age peaks are reported in Table 1 andFigure 2. All samples except 1707 from the Tucaní catchment record a Late Miocene (6-10 Ma) age peak, and all except those from the northwestern-most catchments (Agua Viva and Chejendé) also record a Mio-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 to 100%) in most samples, except for the Santo Domingo and Coloncito catchments where the younger (Mio-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. 1); for this reason we merged the detrital data for the corresponding samples (0507 and 0807) and proceed to apply the decomposition of the grain-age distribution as above. 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 dataset contains only one river draining to the south. The Santo Domingo River sample shows a similar grainage distribution and peak ages as the rivers that drain to the north (Table 1, Figure 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 these, we use the Kolmogorov-Smirnov and Kuiper tests [START_REF] Conover | Practical Nonparametric Statistics: 2 nd ed[END_REF] to quantitatively compare the detrital age probability density functions (PDF) for the seven catchments (cf. Supplementary Material; Figure S1 and Table S1). 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. [START_REF] Bernet | Fundamentals of detrital zircon fission-track analysis for provenance and exhumation studies with examples from the European Alps[END_REF][START_REF] Ruhl | The use of detrital mineral cooling ages to evaluate steady state assumptions in active orogens: An example from the central Nepalese Himalaya: Tectonics[END_REF][START_REF] Brewer | Downstream development of a detrital cooling-age signal[END_REF]. Here, we compare our detrital AFT grain-age distributions of an individual drainage to previously published bedrock AFT ages. We analyse the Chama catchment, as it is the only one within the Venezuelan Andes for which sufficient in-situ data [START_REF] Kohn | Mesozoic-Pleistocene fission track ages on rocks of the Venezuelan Andes and their tectonic implications[END_REF][START_REF] Bermúdez | Asynchronous Miocene-Pliocene exhumation of the central Venezuelan Andes[END_REF] are currently available to derive a reliable bedrock age map.

Several authors [START_REF] Stock | Where does sediment come from? Quantifying catchment erosion with detrital apatite (U-Th)/He thermochronometry[END_REF][START_REF] Vermeesch | Quantitative geomorphology of the White Mountains (California) using detrital apatite fission track thermochronology[END_REF][START_REF] Mcphillips | Using tracer thermochronology to measure modern relief change in the Sierra Nevada, California[END_REF] have developed quantitative methods for comparing bedrock and detrital age patterns. We adopt the approach described by [START_REF] Vermeesch | Quantitative geomorphology of the White Mountains (California) using detrital apatite fission track thermochronology[END_REF], randomly sampling a number of ages equal to the number of single-grain ages in the detrital sample from the interpolated bedrock age map (see Glotzbach et al., in review for details). This method has the advantage that it takes into account Poisson-distributed measurement uncertainties. We sample the bedrock age distribution and construct the predicted age probability-density function (PDF) and cumulative-density function (CDF) 1000 times (Fig. 3a). We then use the Kolmogorov-Smirnov (KS) [START_REF] Conover | Practical Nonparametric Statistics: 2 nd ed[END_REF] and Kuiper equality tests [START_REF] Ruhl | The use of detrital mineral cooling ages to evaluate steady state assumptions in active orogens: An example from the central Nepalese Himalaya: Tectonics[END_REF] to compare the predicted AFT CDFs and PDFs with the detrital AFT data. Additional simulations (5000; 10000) 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 KS 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., [START_REF] Tranel | Spatial variation of erosion in a small, glaciated basin in the Teton Range, Wyoming, based on detrital apatite (U-Th)/He thermochronology[END_REF]. However, [START_REF] Bookhagen | Orographic barriers, high-resolution TRMM precipitation, and relief variations along the eastern Andes[END_REF] we have no information on the relative apatite abundance in different rock types of the Venezuelan Andes, (2) the sampled Chama catchment is mainly underlain by relatively homogeneous Precambrian basement rocks (see Chama River in Fig. 1B).

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 [START_REF] Vermeesch | Quantitative geomorphology of the White Mountains (California) using detrital apatite fission track thermochronology[END_REF][START_REF] Mcphillips | Using tracer thermochronology to measure modern relief change in the Sierra Nevada, California[END_REF]Glotzbach et al., in review). 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 1D thermal model (cf. next section). We thus resample 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 CDF's/PDF's compare very favourably to the observed detrital AFT age distribution (Fig. 3B); both the KS and Kuiper tests suggest the two distributions are similar at the 95% confidence level.

We conclude from the above 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 employ a simple 1-D steady-state thermal model developed by [START_REF] Brandon | Late Cenozoic exhumation of the Cascadia accretionary wedge in the Olympic Mountains, northwest Washington State[END_REF] see also [START_REF] Ehlers | Crustal thermal processes and the interpretation of thermochronometer data[END_REF][START_REF] Reiners | Using thermochronology to understand orogenic erosion[END_REF] to convert detrital AFT ages to exhumation rates. As this model does not take either transient thermal effects (e.g. [START_REF] Rahl | Quantifying transient erosion of orogens with detrital thermochronology from syntectonic basin deposits[END_REF] or the 3-D effects of topography (e.g., [START_REF] Whipp | Effect of exhumation kinematics and topographic evolution on detrital thermochronometer data[END_REF] 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 km 2 My -1 , heat production of 10 °C My -1 and a temperature at the base of the layer of 700 °C (see Bermúdez et al., 2011 for details). The last two values were constrained by inverse modeling of age-elevation relationships in the central Venezuelan Andes [START_REF] Bermúdez | Asynchronous Miocene-Pliocene exhumation of the central Venezuelan Andes[END_REF] 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 (T c ) and -depth (z c ) from the AFT ages, using the [START_REF] Dodson | Closure temperature in cooling geochronological and petrological systems[END_REF] equation to solve for the closure temperature of the AFT thermochronometer and a steadystate thermal structure to relate this to the closure depth.

For each detrital sample, we use the above procedure to translate single-grain ages into exhumation rates, which we then combine to predict a catchment-wide average long-term exhumation (ε T ) rate as: [START_REF] Bookhagen | Orographic barriers, high-resolution TRMM precipitation, and relief variations along the eastern Andes[END_REF] 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 calculate a weighted long-term catchment-averaged exhumation rate (ε T ) w :

(2)

Results are reported in Table 2. Catchment-averaged exhumation rates vary between 0.48±0.02 km Myr -1 for the Santo Domingo and 0.80±0.26 km Myr -1 for the Mimbós catchments. 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 (Agua Viva-Chejendé) and 0.48±0.08 (Coloncito) km Myr -1 . The northeastern Agua Viva/Chejendé catchment shows the lowest long-term exhumation rates (< 0.35 km Myr -1 ). The central Chama and Santo Domingo catchments show intermediate rates (0.35-0.40 km Myr -1 ) and the catchments located on the northwestern flank of the range show the highest exhumation rates (> 0.40 km Myr -1 ). In the following discussion, we will focus on the weighted-average catchment exhumation rates, which we consider to be the most reliable estimates.

POTENTIAL CONTROLLING PARAMETERS

In order to explore the potential topographic, tectonic or climatic controls on the exhumation rates determined above, we derive quantitative measures for these controls (relief, precipitation, seismicity, 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 (Fig. 1B). This part of the Andes most probably emerged as an orographic barrier at about 8 Ma [START_REF] Hoorn | Andean tectonics as a cause for changing drainage patterns in Miocene northern South America[END_REF][START_REF] Bermúdez | Asynchronous Miocene-Pliocene exhumation of the central Venezuelan Andes[END_REF]. The central Sierra Nevada and Sierra la Culata ranges (Fig. 1C) 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 [START_REF] Schubert | The Pleistocene and recent extent of the glaciers of the Sierra Nevada de Merida, Venezuela[END_REF], 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 [START_REF] Bermúdez | Spatial and temporal patterns of exhumation across the Venezuelan Andes: Implications for Cenozoic Caribbean geodynamics: Tectonics[END_REF].

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., [START_REF] Montgomery | Topographic controls on erosion rates in tectonically active mountain ranges[END_REF]. 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 [START_REF] Champagnac | Tectonics, climate, and mountain topography[END_REF]. 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 (Fig. 1B,4A). The former are characterized by extensive low-relief summit surfaces that show widespread evidence for glacial bevelling through cirque retreat (e.g., [START_REF] Mitchell | Influence of a glacial buzzsaw on the height and morphology of the Cascade Range in central Washington State[END_REF][START_REF] Foster | Small valley glaciers and the effectiveness of the glacial buzzsaw in the northern Basin and Range[END_REF], whereas the latter show subdued relief at lower mean elevations with widespread soil mantling of slopes.

Precipitation

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 centre 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 about 8 Ma, in response to tectonically driven surface uplift [START_REF] Hoorn | Andean tectonics as a cause for changing drainage patterns in Miocene northern South America[END_REF][START_REF] Bermúdez | Asynchronous Miocene-Pliocene exhumation of the central Venezuelan Andes[END_REF]. We have used precipitation data from 30 meteorological stations in and around the Venezuelan Andes to elaborate a mean annual precipitation map for the past twenty years.

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) NOAA World Temperature -Precipitation dataset (http://bonnet19.cs.qc.edu:7778/pls/rschdata/rd_start.main); (2) published data [START_REF] Stansell | Last glacial maximum equilibrium-line altitude and paleo-temperature reconstructions for the Cordillera de Mérida[END_REF][START_REF] Naranjo | Estimación de la oferta de agua superficial y conflictos de uso en la cuenca alta del Río Chama[END_REF] Espejo (MARS). The present-day precipitation map shown in Figure 4B was compiled from these data using nearest-neighbor interpolation [START_REF] Arya | An optimal algorithm for approximate nearest neighbor searching in fixed dimensions[END_REF] and 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 (TRMM) precipitation dataset synthesized by [START_REF] Bookhagen | Orographic barriers, high-resolution TRMM precipitation, and relief variations along the eastern Andes[END_REF] shows overall similar patterns, both qualitatively and quantitatively.

Seismicity

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 km 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 about 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 [START_REF] Dewey | Seismicity and tectonics of western Venezuela[END_REF][START_REF] Niu | Receiver function study of the crustal structure of the southeastern Caribbean plate boundary and Venezuela[END_REF]. A seismic zone of intermediate depth (~160 km) occurs towards 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 [START_REF] Schneider | Microseismicity and focal mechanisms of the intermediate depth Bucaramanga Nest, Colombia[END_REF] 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, leftlateral 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 [START_REF] Colmenares | Stress field and seismotectonics of northern South America[END_REF][START_REF] Corredor | Seismic strain rates and distributed continental deformation in the northern Andes and three-dimensional seismotectonics of northwestern South America[END_REF][START_REF] Cortés | Current states of stress in the northern Andes as indicated by focal mechanisms of earthquakes[END_REF] Andes, together with regional tectonic reconstructions, suggest that the current tectonic regime was installed during Pliocene-Quaternary times [START_REF] Backé | Spatial and temporal relationships between compression, strike-slip and extension in the Central Venezuelan Andes: Clues for Plio-Quaternary tectonic escape[END_REF][START_REF] Egbue | Pleistocene to Present North Andean "escape[END_REF].

In order to quantify the effects of seismicity across the Venezuelan Andes, we compiled a seismicity record over the last 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 Venezuelan Foundation for Seismological Investigation (FUNVISIS; http://www.funvisis.gob.ve). The data were filtered using two criteria: (1) epicenters of earthquakes located between latitudes 72.25°W -70.00°W and longitudes 7.75°N -10.00°N; (2) only data with reported local magnitudes M l were considered. Figure 5A shows a summary of earthquakes that occurred in the Venezuelan Andes during the last century. We calculated the released seismic energy (Se) from the local magnitudes using the classical expression of [START_REF] Gutenberg | Seismicity of the Earth and Associated Phenomena[END_REF]:

(3)

The values of the parameters a and b are 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 [START_REF] Gutenberg | Seismicity of the Earth and Associated Phenomena[END_REF]. Released seismic energy values were cumulated 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 M l > 5 earthquakes that occurred on the Burbusay, Valera and Icotea fault systems.

The record of seismic energy release spans only 100 years, significantly shorter than the return period of major earthquakes in the Venezuelan Andes (300 years for M ≥ 7 earthquakes; [START_REF] Audemard | Holocene and historical earthquakes on the Boconó fault system, southern Venezuelan Andes: Trench confirmation[END_REF]. 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. [START_REF] Holt | On the determination of selfconsistent strain rate fields within zones of distributed continental deformation[END_REF]. 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 timescales longer than the observation interval, using the observed earthquake magnitude-frequency (Gutenberg-Richter) relationship. To achieve this, we use the method described by [START_REF] Braun | Constraints on the current rate of deformation and surface uplift of the Australian Continent from a new seismic database and low-T thermochronological data[END_REF] and calculate seismic/brittle strain rate as:

(4) in which the parameters a and b are defined and calculated as explained above; M max is the maximum observed magnitude; µ is elastic shear modulus; Δ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 years).

The depth of the maximum magnitude earthquakes (33-163 km) generally exceeds the depth of brittle-ductile transition (~15-20 km). Because we are using 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, r 2 < 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 2001/12/21, M w = 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 centre 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 [START_REF] Finlayson | Spatial coincidence of rapid inferred erosion with toung metamorphic massifs in the Himalayas[END_REF]. 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 [START_REF] Wilson | Terrain Analysis: Principles and Applications[END_REF][START_REF] Wobus | Tectonics from topography: procedures, promise and pitfalls[END_REF]. The analysis is based on a prediction of bedrock incision rate as a function of stream power [START_REF] Finlayson | Spatial coincidence of rapid inferred erosion with toung metamorphic massifs in the Himalayas[END_REF][START_REF] Tucker | Topographic outcomes predicted by stream erosion models: Sensitivity analysis and intermodel comparison[END_REF]:

(5

)
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:

EI = A m S n (6)
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:

€ EI p = A p ( ) P [ ] m S n ∑ (7)
where (A p ) is the pixel area, P is the local precipitation and the summation sign implies summing the along the flow-lines within the catchment in order to calculate the flow accumulation. We term EI p 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 [START_REF] Finlayson | Spatial coincidence of rapid inferred erosion with toung metamorphic massifs in the Himalayas[END_REF][START_REF] Tucker | Topographic outcomes predicted by stream erosion models: Sensitivity analysis and intermodel comparison[END_REF]. In the case where incision is controlled by total stream power (TSP), m = n = 1. For incision controlled by stream power 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 (Figure 4B), with maximum erosion potential shifting to the southeastern flank of the mountain belt (Figure 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).

DISCUSSION: RELATIONS BETWEEN TECTONICS, CLIMATE AND EROSION IN THE VENEZUELAN ANDES

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 [START_REF] Montgomery | Climate, tectonics, and the morphology of the Andes[END_REF], Washington Cascades [START_REF] Reiners | Coupled spatial variations in precipitation and long-term erosion rates across the Washington Cascades[END_REF], while others argued for decoupling between long-term exhumation and precipitation, such as in the central Himalaya [START_REF] Burbank | Decoupling of erosion and precipitation in the Himalayas[END_REF] or in the European Alps [START_REF] Vernon | Spatial correlation between longterm exhumation rates and present-day forcing parameters in the western European Alps[END_REF]. Recently, [START_REF] Mora | Climatic forcing of asymmetric orogenic evolution in the Eastern Cordillera of Colombia[END_REF] have 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 which timescales? 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 presentday precipitation or seismicity data that was collected over the past 100 years or less.

In order to gain insight into the potential controls on erosion and exhumation rates in the Venezuelan Andes, we have 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.88) between long-term exhumation rates and 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 [START_REF] Wittmann | The relation between rock uplift and denudation from cosmogenic nuclides in river sediment in the Central Alps of Switzerland[END_REF][START_REF] Vernon | Spatial correlation between longterm exhumation rates and present-day forcing parameters in the western European Alps[END_REF] or in the San Bernardino Mountains of California [START_REF] Binnie | Tectonic uplift, threshold hillslopes, and denudation rates in a developing mountain range[END_REF]. 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 [START_REF] Montgomery | Topographic controls on erosion rates in tectonically active mountain ranges[END_REF] is less important with respect to these other mountain belts. In effect, a large proportion of slopes close to the threshold for landsliding is only encountered in the central Sierra Nevada and Sierra la Culata blocks of the Venezuelan Andes [START_REF] Bermúdez | Spatial and temporal patterns of exhumation across the Venezuelan Andes: Implications for Cenozoic Caribbean geodynamics: Tectonics[END_REF].

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 (EI p ). 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 [START_REF] Mora | Climatic forcing of asymmetric orogenic evolution in the Eastern Cordillera of Colombia[END_REF].

However, two caveats need to be recalled before concluding from this data that climate and exhumation are decoupled in the Venezuelan Andes. First, the relatively limited spatial extent of our sampled catchments, for which the wet southern flank is clearly under-represented, and second the timescale problem mentioned previously. Considering the first, it is clear that our dataset would benefit from inclusion of additional remote south-flank catchments. As for the problem of comparing measurements on strongly varying timescales, 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 [START_REF] Hoorn | Andean tectonics as a cause for changing drainage patterns in Miocene northern South America[END_REF][START_REF] Díaz De Gamero | The changing course of the Orinoco River during the Neogene: a review: Palaeogeography[END_REF]; and (2) Quaternary glaciations, which may strongly modify erosion patterns in some mountain belts (e.g., [START_REF] Burbank | Decoupling of erosion and precipitation in the Himalayas[END_REF][START_REF] Gabet | Modern erosion rates in the High Himalayas of Nepal[END_REF][START_REF] Vernon | Spatial correlation between longterm exhumation rates and present-day forcing parameters in the western European Alps[END_REF] have been relatively limited in the Venezuelan Andes [START_REF] Schubert | The Pleistocene and recent extent of the glaciers of the Sierra Nevada de Merida, Venezuela[END_REF][START_REF] Stansell | Last glacial maximum equilibrium-line altitude and paleo-temperature reconstructions for the Cordillera de Mérida[END_REF]. It has been suggested, both in the Andes and elsewhere (e.g., [START_REF] Bookhagen | Late Quaternary intensified monsoon phases control landscape evolution in the northwest Himalaya[END_REF][START_REF] Abbühl | El Niño forcing on 10 Be-based surface denudation rates in the northwestern Peruvian Andes?[END_REF], that spatial patterns of erosion are strongly 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 use cumulative seismic energy release and seismic strain rate as proxies for tectonic forcing (Tables 2 and3). We estimate the cumulative seismic energy released within the different catchment 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 normalising 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 catchmentaveraged exhumation rate or 10 km-radius relief, with r-values varying between 0.40 and 0.64 (Table 3). As noted previously, the pattern of seismic energy release is strongly controlled by large (M l > 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 10 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 Myr -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 of precipitation. It also suggests 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., [START_REF] Whipple | Bedrock rivers and the geomorphology of active orogens[END_REF][START_REF] Lague | Reduction of long-term bedrock incision efficiency by short-term alluvial cover intermittency[END_REF]. We have previously suggested an overall tectonic control on relief and exhumation rates in the Venezuelan Andes as the orogen can be divided in fault-bounded blocks that each have their characteristic relief and exhumation history [START_REF] Bermúdez | Spatial and temporal patterns of exhumation across the Venezuelan Andes: Implications for Cenozoic Caribbean geodynamics: Tectonics[END_REF], as also observed in other obliquely convergent mountain belts (e.g., [START_REF] Spotila | Patterns of bedrock uplift along the San Andreas fault and implications for mechanisms of transpresion[END_REF]. It would be interesting to map out erosion rates throughout the Venezuelan Andes using cosmogenic nuclides in stream sediments in order to analyse whether the same patterns emerge on shorter time scales.

CONCLUSIONS

Detrital AFT analysis of modern river sediment provides an efficient tool for studying longterm 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 [START_REF] Bermúdez | Spatial and temporal patterns of exhumation across the Venezuelan Andes: Implications for Cenozoic Caribbean geodynamics: Tectonics[END_REF]2011) and can thus be used to infer exhumation rates in areas where in-situ data is absent. Our results suggest that using weighted average exhumation rates, in order to include the effect that rapidly exhuming areas of a catchment have higher sediment yields 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 Myr -1 . Rates are higher (>0.4 km Myr -1 ) on the northwestern flank of the mountain belt (Mimbós, San Pedro, Tucani and Coloncito rivers), whereas they are lowest (0.31 km Myr -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 Myr -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 for driving exhumation and relief development in the Venezuelan Andes, or that the shortterm 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-year 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 studying relationships between long-term exhumation rate, relief and tectonics. Our study suggests that the climatic Sierra Nevada de Mérida (Modified from [START_REF] Colletta | Tectonic inheritance, crustal architecture, and contrasting structural styles in the Venezuela Andes[END_REF].

Figure captions

A                                                                      !  !  !  !  !  !  !  !  Lake Maracaibo
!  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  !  < 1.

  Figure 6. These maps emphasize distinct zones of high erosion potential associated with the

Fig. 1 .

 1 Fig.1. A) Shaded digital elevation model(1-km resolution;[START_REF] Jarvis | Hole-filled SRTM for the globe Version 4[END_REF] showing

Fig

  Fig.2. Shaded 90-m resolution digital elevation model of the Venezuelan Andes(Jarvis et al.,

Fig. 4 .

 4 Fig. 4. (A) 5-km radius relief across the Venezuelan Andes, draped over a shaded relief

Fig. 5 .

 5 Fig. 5. (A) Major active faults systems (modified from[START_REF] Audemard | Map and Database of Quaternary Faults in Venezuela and its Offshore Regions[END_REF] and seismicity

Fig. 6 .

 6 Fig. 6. Average logarithmic stream power index maps across the Venezuelan Andes. Left

Fig. 7 .

 7 Fig. 7. Correlation plots for: A) 5-km relief, B) present-day precipitation, C) seismic energy

  Figure 3. Bermúdez et al.

Figure 5 .

 5 Figure 5. Bermúdez et al.

  Figure 6. Bermúdez et al.

  www.cecalc.ula.ve/redbc/colecciones/colecciones_datos.html); and (4) data provided by Gerard Kopp (pers. comm.) of the Institute for Meteorology and Climate Research, University of Karlsruhe (Germany) for the Mérida Atmospheric Research Station at Pico

; (3) Bio-climatic net stations of Mérida (http://

  . 2. Shaded 90-m resolution digital elevation model of the Venezuelan Andes[START_REF] Jarvis | Hole-filled SRTM for the globe Version 4[END_REF] with location of bedrock apatite fission-track ages (white circles) and detrital samples

	(star symbols). Thick dark gray polygons indicate catchment boundaries. Histograms show
	distribution of detrital AFT ages for the different samples (note logarithmic age scale); the
	grain-age distributions were decomposed into peak age components using binomial peak
	fitting (Stewart and Brandon, 2004). Peak ages (P1, P2, P3) are indicated, with the number of
	grains (percentage in parentheses) that composes each population.
	Fig. 3. Comparison of measured (continuous black line; dashed lines indicate 95% confidence
	limits) and predicted bedrock (grey lines indicate 1000 random samples) cumulative density
	age functions (CDF) for the Chama catchment; inset shows probability-density functions
	(PDF) for comparison. (A) shows a comparison for random (un-weighted) sampling of
	bedrock age distribution; (B) shows results for a weighted age distribution where the
	probability of sampling a particular bedrock age is a function of the associated long-term
	exhumation rate relative to the average exhumation rate (see text for discussion). Numbers
	show average Kolmogorov Smirnov (KS) and Kuiper (KT) statistics, with associated p-
	values, which test the null hypothesis (H 0 : the distributions are indistinguishable). The number
	of simulations that fail the KS and KT test at a significance level p = 0.05 is also indicated.

ACKNOWLEDGEMENTS

This study was supported by the CDCH de la Universidad Central de Venezuela (UCV), Project Number PI 08-00-6219-2006 and ECOS-Nord project V08U01. We thank Juan Flores and Antonio León for support during fieldwork. We acknowledge Dr. Gerhard Kopp of the Institute for Meteorology and Climate Research, University of Karlsruhe, Germany, for providing the meteorological data of 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 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.

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,b;[START_REF] Kirchner | yr, 10 k.y., and 10 m.y. time scales[END_REF]. without considering errors bars, and r w is the correlation coefficient calculated for the [START_REF] Williamson | Least-squares fitting of a straight line[END_REF] regression line.

Supplementary Material


Florent Lafarge

Xavier Descombes

Josiane Zerubia

Marc Pierrot-Deseilligny

Automatic Building Extraction from DEMs using an Object Approach and Application to the 3D-city Modeling

Keywords: building extraction, high resolution satellite data, Digital Elevation

In this paper, we present an automatic building extraction method based on an object approach from Digital Elevation Models. First, a rough approximation of the building footprints is realized by a method based on marked point processes: the building footprints are modeled by rectangle layouts. Then, these rectangular footprints are regularized by improving the connection between the neighboring rectangles and detecting the roof height discontinuities. The obtained building footprints are structured footprints: each element represents a specific part of an urban structure. Results are finally applied to a 3D-city modeling process.

Introduction

Three dimensional models of urban areas are very useful for many kinds of applications. Urban planning, radiowave reachability tests for wireless communications, disaster recovery or creation of virtual realities require to be able to automatically build 3D-models with connected planar facets of towns.

However the 3D building reconstruction is a difficult problem, mainly due to the complexity of the scenes. Urban environments are very dense and composed of many types of buildings. Therefore, it makes their analysis difficult.

Many methods have been developed in various contexts which depend on the automaticity of the methods and the kind of data. [START_REF] Flamanc | Evaluation of 3D city model production from PLEIADES HR satellite images and 2D ground maps[END_REF] and [START_REF] Gruen | CC-Modeler : a topology generator for 3D city models[END_REF] propose semi-automatic processes which require human operator interventions. However, the 3D building reconstruction demands automaticity which is the case of the methods mentioned below. In [START_REF] Lin | Building detection and description from a single intensity image[END_REF], the input is a single view image -building shadows are used to extract 3D information. Multiple view aerial images are the most common input. Such data are rich in terms of 3D information and allow to extract 3D primitives. For example, [START_REF] Scholze | A probabilistic approach to building roof reconstruction using semantic labelling[END_REF] and [START_REF] Kim | Automatic description of complex buildings from multiple images[END_REF] use 3D lines in their 3D reconstruction methods. In [START_REF] Jibrini | Automatic building reconstruction from very high resolution aerial stereopairs using cadastral ground plans[END_REF], planes are extracted and then arranged in a polyhedral approach. [START_REF] Baillard | A plane-sweep strategy for the 3D reconstruction of buildings from multiple images[END_REF] and [START_REF] Taillandier | Automatic Builings Reconstruction from Aerial Images : a Generic Bayesian Framework[END_REF] propose generic models by combining several kinds of primitives. Laser scanning is also very popular since the decrease of the acquisition cost and the 2 accuracy of the measures. Some interesting models have been proposed using laser data such as [START_REF] Maas | Two algorithms for extracting building models from raw laser altimetry data[END_REF][START_REF] Haala | Extraction of buildings and trees in urban environments[END_REF][START_REF] Früh | Constructing 3D city models by merging aerial and ground views[END_REF]. Several methods, such as [START_REF] Vinson | Multiple rectangle model for buildings segmentation and 3D scene reconstruction[END_REF][START_REF] Weidner | Towards Automatic Building Reconstruction from High Resolution Digital Elevation Models[END_REF], are based on Digital Elevation Models (DEMs).

The use of DEMs is motivated by the fact that they already provide a geometric description of a scene from aerial imagery or airborne scanner data.

Nowadays, this problem can be tackled by another type of data: the submetric satellite images. The future PLEIADES satellites are especially well adapted to deal with 3D building reconstruction through the high resolution of images (0.7 meter) and the stereoscopic characteristics (three views -B/H=0.2). The main advantages of satellite data compared to aerial images are a high swath width and ground coverage. However, such data have a "relatively low" resolution and a "low" Signal to Noise Ratio (SNR) to deal with 3D reconstruction problems compared to aerial data. For example, PLEIADES simulations have 2 pixels per square meter density contrary to aerial images used in [START_REF] Taillandier | Automatic Builings Reconstruction from Aerial Images : a Generic Bayesian Framework[END_REF] or [START_REF] Fuchs | Efficient subgraph isomorphism with a-priori knowledge. Application to building reconstruction for cartography[END_REF] which have about 25 pixels per square meter. In this paper, we aim at automatically extracting buildings by exploiting satellite images, and more precisely PLEIADES simulations.

A robust process is necessary to face up the "low" quality of such data. Basic 3D-primitive extraction is not especially adapted to our data quality. For instance, the 3D-segment extraction by using [START_REF] Taillandier | Reconstruction of 3d linear primitives from multiple views for urban areas modelisation[END_REF] does not provide convincing results. This context leads us to propose a method based on important prior knowledge concerning urban structures. An object approach is well adapted since it allows the introduction of prior knowledge on the object layout. Most of the methods using such an approach need building localization masks in order to initialize the footprint extraction algorithm, as for example active contours based methods in [START_REF] Oriot | Statistical snakes for building extraction from stereoscopic aerial images[END_REF][START_REF] Bailloeul | Spatio-temporal prior shape constraint for level set segmentation[END_REF].

We propose an object approach which does not need localization information and only use a single input well adapted to geometric description: a DEM. DEMS are generated from 3-view satellite images (more precisely from PLEIADES simulations) by an algorithm detailed in [START_REF] Pierrot-Deseilligny | A multiresolution and optimization-based image matching approach : an application to surface reconstruction from SPOT5-HRS stereo imagery[END_REF] and based on a multi-resolution implementation of a Cox and Roy optimal flow matching image algorithm [START_REF] Roy | A maximum-flow formulation of the n-camera stereo correspondence problem[END_REF]. PLEIADES simulations were derived from aerial images which have been degraded by the French Space Agency, taking into account various PLEIADES satellite characteristics such as the image resolution and the SNR. In this paper, an approximation of the global shape of buildings is first extracted by a method based on marked point processes: the building footprints are modeled by rectangle layouts. Then, these rectangular building footprints initiate a process leading to a polygonal description of the footprints while detecting the roof height discontinuities. The obtained building footprints are structured footprints: each element represents a specific part of an urban structure. Then, these results drive a 3D-city modeling process.

2 Rectangular approximation of the building footprints by marked point processes

The first step consists in extracting a rough description of the global shape of buildings. High density of urban areas and complexity of human made objects make the building extraction process difficult to achieve. We aim at developing an automatic method which is robust with respect to the DEM noise and does not need building localization maps or initialization states. To do so, a rectangular footprint extraction algorithm based on marked point processes is used [START_REF] Ortner | Processus ponctuels marqués pour l'extraction automatique de caricatures de bâtiments à partir de modèles numériques d'élévation[END_REF][START_REF] Ortner | Building outline extraction from Digital Elevation Models using marked point processes[END_REF].

Principle

This method is based on an object approach: it provides rectangle layouts which represent the building footprints. Faced with the complexity and the diversity of the building forms, such an approach is well adapted since it provides a modeling by simple geometric objects (i.e rectangles) and allows the introduction of a prior knowledge concerning the object layout and doesn't need any initialization. An energy is associated to each configuration. The global minimum of this energy is then found by applying a Reversible Jump Markov Chain Monte Carlo (RJMCMC) sampler embedded into a simulated annealing scheme [START_REF] Green | Reversible Jump Markov Chains Monte Carlo computation and Bayesian model determination[END_REF].

Marked point process

Let us consider a point process X living in K = [0, X max ] × [0, Y max ] which represents an image. X is a measurable mapping from a probability space (Ω, A, P) to configurations of points of K:

∀ω ∈ Ω, x i ∈ K, X(ω) = {x 1 , ..., x n } (1)
A marked point process in S = K × M is a point process where each point is associated with a mark which allows to define an object. M is the mark space which allows, in our case, to describe random configurations of rectangles:

M = [- π 2 , π 2 ] × [L min , L max ] × [l min , l max ] (2)
We denote C the set of all finite configurations of objects of S. 6

Density formulation

We consider a marked point process X defined by an unnormalized density h(.) from C to R and a reference intensity measure ν (defining a reference

Poisson point process distribution).

The unnormalized density h(.) derives from three energies presented in the following. More details concerning these energy terms can be found in [START_REF] Ortner | Building outline extraction from Digital Elevation Models using marked point processes[END_REF][START_REF] Ortner | Processus ponctuels marqués pour l'extraction automatique de caricatures de bâtiments à partir de modèles numériques d'élévation[END_REF].

For a configuration x ∈ C, the density is given by:

h(x) = exp -[U ext (x) + U int (x) + U excl (x)]
(3)

External energy U ext

This term measures the quality of a configuration given the data, i.e. the DEM. It can be decomposed as a sum of energies U d (.) per object such as

U ext (x) = u∈x U d (u).
For each object u ∈ x, a local data energy U d (u) is computed by using a method which compares the DEM discontinuities to the borders of the object.

Internal energy U int

The internal energy U int (x) allows to give a particular spatial structure to the configuration x. It can be seen as a regularizing term composed of three kinds of interactions between objects which favor alignment between rectangles, completion and some paving behavior. Fig. 3 represents these three types of interactions.

Exclusion energy U excl

This term avoids redundant objects by penalizing the intersection of parallel rectangles. Such an exclusion energy allows to ensure that the attractive interactions do not provoke object accumulations.

Optimization

Once the model is defined, we want to find the configuration maximizing the density h(.). We search for the Maximum A Posteriori (MAP) estimator x M AP :

x M AP = arg max x h(x) (4)
This is a non convex optimization problem in a high and variable dimension space C. A Reversible Jump Monte Carlo Markov Chain sampler [START_REF] Geyer | Simulation and likelihood inference for spatial point processes[END_REF][START_REF] Green | Reversible Jump Markov Chains Monte Carlo computation and Bayesian model determination[END_REF][START_REF] Van Lieshout | Markov point processes and their applications[END_REF] embedded in a simulated annealing scheme is particularly well adapted to this problem.

RJMCMC

It consists in simulating a discrete Markov Chain (X t ) t∈N on C, having π as an invariant measure (specified by the density h(.)), which performs small jumps between spaces of variable dimensions. One of the main advantages of such a sampler is that the chain asymptotically converges towards π for any initial configuration X 0 . It means that the initialization steps or building localization masks are not needed. The transitions of this chain are managed by some proposition kernels (denoted Q i ), which are moves (or jumps) proposed from the current configuration.

The RJMCMC algorithm is described below:

Algorithm 1 . At iteration k, with X k = x k :
• Choose uniformly one of the proposition kernels Q i (x k , .).

• According to Q i , propose a new state y.

• Compute the acceptance ratio of this move

α i (x k , y) = min π(y)Q i (y, x k ) π(x k )Q i (x k , y) , 1 (5)
• take X k+1 = y with probability α i (x k , y), and X k+1 = x k otherwise.

Proposition kernels

The proposition kernels allow to define the transformations applied to objects in the current configuration. The kernels we used are birth, death, rotation, dilatation and translation (see Fig. 4). Birth and death kernels allow to respectively add / remove an object from / in the current configuration. These two transformations, which correspond to jumps in spaces of higher (birth) / lower (death) dimension, guaranty that the Markov chain visit the whole configuration space. However, it is important to define relevant moves in order to speed up the convergence of the Markov chain. That is why we introduce other moves (rotation, dilatation and translation) which correspond to perturbations of the current configuration (i.e. without jumping in a space of different dimension).

The computation of these kernels are specified in [START_REF] Ortner | Building outline extraction from Digital Elevation Models using marked point processes[END_REF].

Simulated annealing

In practice, the RJMCMC sampler is embedded into a simulated annealing scheme: the density h(.)

1 D t is substituted to h(.)
where D t is a sequence of temperatures which tends to zero as t tends to infinity. At the beginning of the algorithm (i.e. when the temperature is high), the process is not really selective: it allows to explore the density modes. When the temperature decreases, configurations which have a high density are favored. Fig. 5 shows an example of this optimization process.

The simulated annealing theoretically ensures the converge to the global optimum for any initial configuration x 0 by using a logarithmic temperature decrease (see [START_REF] Salamon | Facts, conjectures, and improvements for simulated annealing[END_REF]). In practice, we prefer to use faster decreases which provide an approximate solution close to the optimal one. More precisely, we use a geometrical decrease, detailed in [START_REF] Van Laarhoven | Simulated Annealing : Theory and Applications[END_REF], and defined as follows:

D t = D 0 .α t (6)
where α < 1 and close to 1 and D 0 is the initial temperature. In practice,

we take α = 0.99999. The initial and ending temperatures D 0 and D f inal are estimated through the variations of the energy by [START_REF] White | Concepts of scale in simulated annealing[END_REF].

Comments

Fig. 6 presents two results from PLEIADES simulations on buildings of the city of Amiens (France). The results are satisfactory considering a satellite and automatic context. However, applying directly a 3D reconstruction process from these rectangular building footprints engenders artifacts in the 3Dresults as it is shown in (Lafarge et al., 2006a). These artifacts are mainly due to non optimal connections between neighboring rectangles (see Fig. 7).

The proposed approach is fully automatic but requires heavy computation. To reduce the computation time, we had proposed to only estimate a caricature of the building footprints consisting of a collection of rectangles. To reach a 3D reconstruction, these rectangle layouts need to be regularized, especially the linking up between neighboring rectangles. We address this problem in the following section.

3 Regularization: towards structured footprints

The rectangle layouts are regularized by improving firstly the connections between neighboring rectangles (through a local fusion of rectangles) and, secondly, the position of facades (through a roof height discontinuity detection). At the end, we obtain structured building footprints whose elements are correctly connected: each one represents a specific part of an urban structure.

Connection between neighboring rectangles

We aim at improving the connection between neighboring rectangles. Fig. 8-(a) represents a typical configuration of two neighboring rectangles that we want to connect.

Generalizing 3D volumes, such as in [START_REF] Kada | 3D building generalization based on half-space modeling[END_REF], could allow to improve the linking up between neighboring rectangles, nevertheless some building details would be lost. Applying small perturbations on the rectangle parameters could improve the connection. However, such a method cannot be used since the information related to the general layout of rectangles will be lost. The chosen solution consists in merging the neighboring rectangles, i.e. morphing two rectangles into polygons which are perfectly connected as we can see on Fig.

8-(b).

(a) (b)

Fusion

We aim at merging two neighboring rectangles (R 1 , R 2 ), which means that we It exists many relative positions of two neighboring rectangles which depend on the width of the rectangles, their orientation and their distance. Faced with the high number of possibilities, an exhaustive description of all the fusion configurations will be very heavy. We have chosen to restrict the number of possibilities to 16 fusion configurations. 16

Cost function

In order to select the most adapted fusion configuration, we define a cost function which measures the coherence between the proposed configuration and the data represented by both the DEM and the two initial rectangles.

This cost is composed of three terms:

C = C DEM + ω 1 C recovering + ω 2 C contour (7)
where ω 1 and ω 2 represent weights which are tuned by trial and error. In practice, we choose ω 1 = 0.5 and ω 2 = 0.8. These values are stable for all scenes and various ranges of DEMs since they weight geometric criterion related to rectangle layouts provided by the same algorithm.

DEM cost

The DEM term allows to define the quality of the proposed configuration, given the DEM. To do so, we estimate the rate of high Z values which are inside the proposed footprint (i.e. pixels which correspond to pieces of building). The higher the rate, the lower the cost. The high intensity pixels are detected by thresholding: the threshold value corresponds to half a floor height with respect to the ground altitude. Let us consider N, the number of these high intensity pixels inside the proposed footprint, and N t , the area of the footprint (i.e. the total number of pixels). Then, the DEM cost can be written as follows:

C DEM = 1 - N N t (8)

Recovering cost

This term provides a area criterion: it corresponds to the recovering rate between the proposed configuration and the two rectangles. Let us consider R 1 and R 2 , the two initial rectangles, and P , the union of polygons defined by the proposed fusion. Then, the recovering term is given by:

C recovering = | Area(P ) -Area(R1 ∪ R2) | Area(R1 ∪ R2) (9)
For example, the configuration of Fig. 10-(b) is favored (in term of recovering cost) contrary to the one of Fig. 10-(a).

Contour cost

This term allows to estimate the contour coherence between the proposed configuration and the two initial rectangles. This cost penalizes complex fusions as we can see on Fig. 11. Let us consider L in (respectively L out ), the total length of the inside sector polygons P (respectively outside sector) and L 1 and L 2 , the length of the two rectangles. The contour cost is defined as follows:

C contour = L out + L in 2 × (L 1 + L 2 ) (10)

Algorithm

Rectangular footprints are morphed into polygonal footprints whose elements are perfectly connected. To do so, the fusion process between two rectangles is applied for each couple of neighboring rectangles which are positioned in series (a rectangle is positioned in series if it has at most two neighbors). These fusions "in series" are well adapted to these rectangle layouts since a large majority of rectangles is positioned in series. It is due to the internal field of the rectangular footprint extraction method (see section 2.3.2) which allows to regularize the rectangles, and especially favors their alignments.

The algorithm consists in merging all the neighboring rectangles which are positioned "in series". However, some buildings owning complex junctions cannot be considered as a single sequence of quadrilaterals since some neighboring rectangles are positioned "in parallel" as we can see in Fig. 12. In this case, the result corresponds to several sequences of quadrilaterals where the main sequence is realized by considering the longest train of fusions (i.e. the longest sequence of rectangles). Consequently, instead of considering parts of buildings represented by rectangles, we work on complete buildings. Each building b is associated with a set of connected polygons P b . Fig. 13 presents the result from the rectangular footprints of Fig. 6 (left). This result is satisfactory: the artifacts of rectangle connections have disappeared.

Roof height discontinuity detection

The second improvement consists in detecting the roof height discontinuities in the building footprints. In fact, several roof heights can exist inside a footprint as we can see in Fig. 14-(a).

To do so, a rooftop profile is computed for each footprint (see Fig. 14-(b)).

This profile is obtained by estimating the average roof height (i.e. Z) in a sliding window. Then, the high gradients are detected by applying a thresholding on the profile derivative (see Fig. 14-(c)). The threshold is a parameter which has been set up to 3 (this value corresponds to the Z-variation / X-variation ratio) w.r.t. a reference building containing various roof height discontinuities of one floor height. Finally, the selected gradients are accumulated (Fig. 14-(d)) in order to obtain the roof height discontinuities presented on Fig. 14-(e).

Fig. 15 shows the result from the footprints of Fig. 6. This result represents "structured building footprints" since each element represents a specific part of an urban structure.

Results

Two sets of parameters exist. They correspond to weights which respectively tune energy terms of the global shape extraction step (equation 3) and cost function terms of the regularizing step (equation 7). They are chosen by trial and error. The calibration of these parameters is made for a range of DEMs.

In the future, it could be interesting to work on the estimation of these parameters in order to have a fully-automatic process. rate is satisfactory with respect to the satellite context (9.7%). This rate can be improved by introducing a vegetation mask in the process in order to prevent the detection of trees (see Fig. 16). The error rates are detailed in Tab.

1. The computing time is quite high. For example, 3 hours are necessary to obtain the result of Fig. 17 (Fig. 19-c) This section presents an application of the previous results: the 3D-city modeling. The structured building footprints provide useful knowledge concerning the localization of buildings, the change of rooftop direction and the roof height discontinuities. It means that these footprints help to do the 3D interpretation of the scenes. We propose a simple 3D-reconstruction process using these structured footprints. It is based on a parametric approach which has been chosen since it is particularly robust. However, this approach deals with restricted 3D-models.

Skeleton process of the footprints

The proposed parametric method is based on a skeleton process of the structured footprints which allows to model the rooftops of buildings. It is a very simple process: roofs are constrained to be symmetric two-planes. Such a process has been used in [START_REF] Haala | Interpretation of urban surface models using 2d building information[END_REF][START_REF] Flamanc | Evaluation of 3D city model production from PLEIADES HR satellite images and 2D ground maps[END_REF] for instance. In other words, for each building b, we consider the median axis of the connected polygon set P b as the rooftops. The roof forms are fully determined by the structured building footprints: it allows to reduce the computing time compared to other approaches such as (Lafarge et al., 2006a;[START_REF] Jibrini | Automatic building reconstruction from very high resolution aerial stereopairs using cadastral ground plans[END_REF][START_REF] Taillandier | Automatic Builings Reconstruction from Aerial Images : a Generic Bayesian Framework[END_REF].

Setting of the building heights

The different heights of buildings constitute the parameters of the model. For each polygon of P b , a roof gutter height and a rooftop height are estimated by a method developed in a previous work [START_REF] Lafarge | A parametric model for automatic 3D building reconstruction from high resolution satellite images[END_REF]. These heights are then regularized on all the building b by using the "K-means" algorithm including a regularizing term. A similar algorithm using an entropic term has been proposed in [START_REF] Palubinskas | An unsupervised clustering method using the entropy minimization[END_REF]. The advantage of such an algorithm is that the number of classes does not need to be known a priori. It allows to favor the rooftop continuity and reduces the presence of artifacts.

Let us consider a set of N connected polygons P b representing the footprint of a building b (see, for example, Fig. 13) and Hg i , the estimated height of the roof gutter of the i th polygons (i ∈ [1, N]). The proposed "K-means" algorithm including a regularizing term consists in searching for c j , the centroid of the j th class, which minimizes the criterion J defined as follows:

J = C j=1 N i=1 d 2 (Hg i , c j ) -αC (11)
where C represents the number of classes (which is unknown).

The criterion J is composed of two terms. The first one corresponds to the classical criterion of the "K-means" algorithm. The second one, which is weighted by the parameter α ∈ R +⋆ , represents a priori knowledge on the classes. This term linearly favors a low number of classes. Each height Hg i is associated with its centroid c j which represents the final regularized height of the roof gutter. The parameter α is tuned by trial and error for a range of DEMs (in our case, we choose α = 100). If the value of α decreases, models will have a bad regularization. On the contrary, if α is too high, the roof height discontinuities will disappear. The same process is applied for the rooftop heights. 29

Results

Fig. 21 presents an example of building 3D-reconstruction (a) associated with the structured building footprints (b) and the 3D-ground truth (c). This result is satisfactory compared to the 3D-ground truth. The proposed modeling does not allow to represent some details. However, the global shape is respected and the generalization level is acceptable for satellite data in an automatic context.

Fig. 22 shows the 3D-results on three Amiens areas used in section 4. These results are satisfactory compared to the associated 3D-ground truths which can be found in (Lafarge et al., 2006b) and have been semi-automatically established by the French Mapping Agency (IGN) from aerial images (0.2 meter resolution) and ground maps. Even if the 3D-modeling is simple (the roof forms are constrained to be symmetric two-planes), a majority of urban structures is close to reality. Moreover, these results are aesthetically satisfactory compared, for example, to results presented in [START_REF] Lafarge | A parametric model for automatic 3D building reconstruction from high resolution satellite images[END_REF]. The computing time is negligible with respect to the footprint extraction process (11 seconds for the result presented in Fig. 22-b).

Conclusion

A method for the building footprint extraction has been proposed in this paper. This method is automatic: building localization maps or initialization states are not needed. It provides structured footprints where each element represents a specific part of a building and is perfectly connected to the neighboring ones. The obtained results are satisfactory and the method seems to be efficient on dense urban areas using high resolution satellite data.

An extension to the 3D building reconstruction has been proposed: it gives satisfactory 3D city models. Even if this method is very simple (the roof forms are constrained to be symmetric two-planes), a majority of urban structures is close to reality and results are aesthetically correct thank to the lack of artifacts in the structured building footprints.

In the future, we should extend the roof form possibilities in order to have more realistic models. To do so, we could use a grammar of 3D-models which include more complex roof forms such as dissymmetric two-planes or gambrel roofs. Moreover, it could be interesting to work on parameter estimation in order to set a fully-automatic process.

Fig. 1 -Fig. 1 .

 11 Fig. 1. A realization of a point process in K (a), a realization of a marked point process in S (b), an element of S -(X c , Y c ) ∈ K, the center of the rectangle -(θ, L, l) ∈ M , the orientation of the rectangle, its length and its width (c).

Fig. 2 -

 2  shows a part of a DEM , a rectangle (in white) and slices evenly disposed. We aim at extracting some points of interest from each slice and checking the coherence between these points and the rectangular shape of the object. Points of interest are detected using a profile simplification algorithm detailed in[START_REF] Ortner | Automatic 3D land register extraction from altimetric data in dense urban areas[END_REF] -they represent DEM discontinuities on the slices (see Fig.2-(b)). Then, three criteria of coherence are used to compute U d (u). The first criterion measures a volume rate (see Fig.2-(c) -grey lines represent segments used to compute the volume rate). The second one corre-

Fig. 2 .

 2 Fig. 2. External energy: a part of a DEM with a proposed rectangle (in white) and its slices (a), detected points of interest (b), volume rate (c), moment rate (d), localization rate (e).

Fig. 3 .

 3 Fig. 3. Internal energy: alignment interaction (a), paving interaction (b), completion interaction (c).

Fig. 4 .

 4 Fig. 4. Mixture of the kernels.

Fig. 5 .

 5 Fig. 5. Example of the optimization process -Evolution of the configuration (top) associated with the energy decrease and temperature decrease schemes (bottom).

Fig. 6 .

 6 Fig. 6. Examples of building extraction results.

Fig. 8 .

 8 Fig. 8. Two neighboring rectangles (a), example of a fusion (b).

  search for connecting both points A and B (outside sector), and points C and D (inside sector) of Fig. 8-(a) in order to create polygons (P i ) i∈N ⋆ which are perfectly connected (see Fig. 8-(b)).

Fig. 9 -Fig. 9 .

 99 Fig. 9-(a) presents 4 connection configurations of the outside sector (it means

Fig. 10 .

 10 Fig. 10. A fusion configuration (o-3 "outside sector -configuration 3"/ i-3 "inside sector -configuration 3") in three different rectangle relative positions.

Fig. 11 .

 11 Fig. 11. Contour cost: configuration penalized by the contour cost (left), favored configuration (right).

Fig. 12 .

 12 Fig. 12. Building footprints represented by a single sequence of rectangles (left) -Building footprints represented by two sequences of quadrilaterals (right).

Fig. 13 .

 13 Fig. 13. Rectangular footprints (left), footprints with fusion of rectangles (right).

Fig. 14 .

 14 Fig. 14. Roof height discontinuity detection : initial footprint (a), rooftop profile (b), gradient selection (c), gradient accumulation (d), resulted footprint with roof height discontinuities (e).

Fig. 15 .

 15 Fig. 15. Building footprints with roof height discontinuities.

Fig. 16 ,

 16 Fig. 16, 17 and 18 show results of structured building footprints (c) associated with PLEIADES simulations (a) and the rectangular footprints (b). These results are satisfactory since the footprints are now regularized: rectangle connections have been improved and roof height discontinuities have successfully been detected.Fig.19represents the error maps. The True Negative (TN) rate is quite high (15.3%). However, TN mainly correspond to low flat buildings of inner courtyards (no more than two floor height) that the proposed method cannot detect since these buildings have low DEM discontinuities. Without taking into account these low flat buildings, the TN rate is 4.5%. The False Positive (FP)

Fig. 16 .

 16 Fig. 16. PLEIADES simulation on a piece of Amiens (a), rectangular building footprints (b), structured building footprints (c).

Fig. 17 .

 17 Fig. 17. PLEIADES simulation on Amiens downtown (a), rectangular building footprints (b), structured building footprints (c).

Fig. 18 .

 18 Fig. 18. PLEIADES simulation on a piece of Amiens (a), rectangular building footprints (b), structured building footprints (c).

Fig. 19 .

 19 Fig. 19. Evaluation maps of Fig. 16, 17 and 18 (respectively (a)&(b)&(c)).

Fig. 20 Fig. 20 .

 2020 Fig. 20 presents examples of a reconstructed building. The first one (a) is obtained without regularizing building heights: we can see the presence of artifacts through the bad adjustments of the rooftop heights and roof gutter heights. The second building (b) is regularized.

  Fig. 21. 3D-result of Amiens city hall (a), structured building footprints (b), 3D-ground truth (c).

Fig

  Fig. 22. 3D-results of Fig. 16, 17 and 18 respectively (respectively (a)&(b)&(c)).

Table 1

 1

Acknowledgments

The first author would like to thank the French Mapping Agency (IGN) and the French Space Agency (CNES) for partial financial support during his PhD.

The authors thank the French Space Agency (CNES) for providing PLEIADES simulations.