Articles | Volume 8, issue 2
Research article
07 Apr 2020
Research article |  | 07 Apr 2020

Hillslope denudation and morphologic response to a rock uplift gradient

Vincent Godard, Jean-Claude Hippolyte, Edward Cushing, Nicolas Espurt, Jules Fleury, Olivier Bellier, Vincent Ollivier, and the ASTER Team

Documenting the spatial variability of tectonic processes from topography is routinely undertaken through the analysis of river profiles, since a direct relationship between fluvial gradient and rock uplift has been identified by incision models. Similarly, theoretical formulations of hillslope profiles predict a strong dependence on their base-level lowering rate, which in most situations is set by channel incision. However, the reduced sensitivity of near-threshold hillslopes and the limited availability of high-resolution topographic data has often been a major limitation for their use to investigate tectonic gradients. Here we combined high-resolution analysis of hillslope morphology and cosmogenic-nuclide-derived denudation rates to unravel the distribution of rock uplift across a blind thrust system at the southwestern Alpine front in France. Our study is located in the Mio-Pliocene Valensole molassic basin, where a series of folds and thrusts has deformed a plateau surface. We focused on a series of catchments aligned perpendicular to the main structures. Using a 1 m lidar digital terrain model, we extracted hillslope topographic properties such as hilltop curvature CHT and nondimensional erosion rates E. We observed systematic variation of these metrics coincident with the location of a major underlying thrust system identified by seismic surveys. Using a simple deformation model, the inversion of the E pattern allows us to propose a location and dip for a blind thrust, which are consistent with available geological and geophysical data. We also sampled clasts from eroding conglomerates at several hilltop locations for 10Be and 26Al concentration measurements. Calculated hilltop denudation rates range from 40 to 120 mm kyr−1. These denudation rates appear to be correlated with E and CHT that were extracted from the morphological analysis, and these rates are used to derive absolute estimates for the fault slip rate. This high-resolution hillslope analysis allows us to resolve short-wavelength variations in rock uplift that would not be possible to unravel using commonly used channel-profile-based methods. Our joint analysis of topography and geochronological data supports the interpretation of active thrusting at the southwestern Alpine front, and such approaches may bring crucial complementary constraints to morphotectonic analysis for the study of slowly slipping faults.

1 Introduction

The topography of the Earth evolves in response to surface processes driven by external forcing of tectonic and climatic origins (e.g., Champagnac et al.2012; Whittaker2012). Information about tectonic uplift is sequentially transmitted through landscapes, firstly by the adjustment of river gradients, which then set the local base level of hillslopes. The water supply and temperature set by climatic conditions control the efficiency of weathering, erosion and transport processes across the Earth's surface. The present-day land surface morphology results from this accumulated actions of tectonic and climatic forcing through time, and a major endeavor of geomorphological research is to interpret measurable topographic properties in terms of space and time variations of either of these tectonic uplift or climatic conditions (e.g., Roberts et al.2012; Demoulin2012; Fox et al.2014). In particular, documenting the spatial variability of tectonic processes from topographic analysis has been a key research focus (Wobus et al.2006), as changes in topographic gradients could record variations in rock uplift rates at various scales, from regional patterns associated with crustal or lithospheric deformation (e.g., Gallen et al.2013) down to differential motion across individual faults (e.g., Boulton et al.2014). Such investigations have often been motivated by the practical concern of identifying high strain zones in tectonically active regions in order to contribute to seismic hazard assessment (e.g., Morell et al.2015).

These studies have been, by a large margin, dominated by approaches relying on the analysis of river long-profile properties, as a direct relationship between fluvial gradient and rock uplift has been identified in many incision models (Tucker and Whipple2002). Notably, the computation of morphometric parameters such as steepness indexes, a measure of channel gradient normalized for drainage area (Kirby and Whipple2012), is now a standard approach when investigating river networks, as this parameter is theoretically dependent on rock uplift and has been shown to be positively correlated with field measurements (e.g., Duvall et al.2004; Wobus et al.2006; Cyr et al.2010). The great development of these methods has been made possible by the increasing availability of medium-resolution digital elevation models (DEMs: 10 to 100 m pixel size), which allow for the reliable extraction of river profiles and accurate computation of along-stream gradients. While robust, efficient and widely used in a variety of settings, these fluvial-based approaches can encounter important limitations when dealing with complexities of river systems such as lithological variations, transient evolution or along-stream changes in fluvial dynamics, but they are also inherently constrained by the planform distribution of river networks, which might not always be optimal to sample the rock uplift patterns.

Similarly to river profiles, theoretical formulations of hillslope denudation predict a strong dependence of morphological parameters, such as slope or relief, on the rate of base level fall set by channel incision. However, key elements of hillslope behavior, such as the threshold stability angle for hillslope material and the nonlinear relationship between sediment fluxes and topographic gradient (Roering et al.1999, 2001, 2007), imply that, for fixed valley spacing, hillslope morphology can be insensitive to changes in erosion or rock uplift rates over a large range of values (Burbank et al.1996; Ouimet et al.2009). This behavior is a major limitation for the use of hillslope morphology to retrieve information about tectonic gradients. Furthermore, from a methodological point of view, the widely available intermediate-resolution DEMs, which are appropriate to describe river profiles developing over 1–100 km length scales, will only provide a very coarse description of hillslopes with typical lengths of the order of 100 m (Grieve et al.2016a).

Nevertheless, this methodological limitation is currently being overcome by the growing use of lidar or photogrammetric techniques delivering digital terrain models (DTMs) with resolutions less than 1 m (James and Robson2012; Glennie et al.2013; Passalacqua et al.2015). The widespread distribution of such data has spurred a great interest in new approaches to extract relevant information at the hillslope scale (e.g., DiBiase et al.2012; Hurst et al.2012; Grieve et al.2016b, c; Milodowski et al.2015; Clubb et al.2020). It is noteworthy that the high resolution of these DTMs allows users to compute accurate derivatives of the topographic surface. According to linear diffusion theory, the second derivative, or curvature, covaries with erosion rate. This relationship remains valid for near-threshold hillslopes in the vicinity of the hilltop where topographic gradients are usually small (Hurst et al.2012). This possibility to access reliable proxies for erosion rate from hillslope scale metrics has led researchers to reconsider the potential of hillslope analysis for the assessment of denudation gradients, which can be used to infer patterns of uplift. Notably, Roering et al. (2007) have proposed a conceptual framework, based on a formulation for nonlinear hillslope sediment flux, which highlights the links between steady-state hillslope morphology and the dynamics of erosion processes as well as the underlying tectonic or climatic forcings. Hurst et al. (2012) have built upon this formulation to construct a joint analysis of high-resolution topographic hillslope metrics and cosmogenic radionuclide (CRN) data in the Sierra Nevada (California). In particular, they used CRN-derived denudation rates to calibrate the efficiency parameter for hillslope transport processes, and they constrain the distribution of absolute erosion rates from hilltop curvature measurements. Similarly, Hurst et al. (2013a) were able to finely track transient landscape adjustment along the San Andreas Fault where long-term motion is progressively moving hillslopes in and out of a high-uplift-rate pressure ridge. This localized change in the tectonic boundary condition is closely recorded by hillslope relief or slope angle and hilltop curvature extracted from lidar data, with the growth and decay phases of landscape evolution leaving a distinctive signature.

These promising results have offered hillslope-based critical insights into the dynamics of transient landscapes, with a spatial density of information several orders of magnitude higher than what could be resolved with approaches based on the fluvial network. Important methodological developments were necessary to extract the relevant information from high-resolution DTMs (e.g., Grieve et al.2016b). However, while theoretically warranted, the applicability of these approaches to explore tectonic gradients has only been tested on a limited number of cases (e.g., Hurst et al.2013a, 2019; Clubb et al.2020). There is, therefore, an urgent need to further investigate hillslope morphological response in various types of tectonic settings to unravel the potential of such methods as an alternative or complement to the routinely applied investigation of river profiles. Additionally, even at high resolution, the analysis of hillslope or river morphological properties can only deliver estimates of the relative intensity of surface processes, the actual rates and efficiency of which can only be obtained through geochronological techniques such as cosmogenic nuclides. When facing the growing availability of such data, there is a critical need to assess whether, and under which circumstances, the high-resolution topographic properties of landscapes and their rates of evolution can both be framed into a coherent picture on the basis of available theoretical formulations for landscape evolution (Dietrich et al.2003).

The objective of our study is to investigate the changes in hillslope morphology, as observed with a lidar DTM, across a rock uplift gradient at the front of the southwestern Alps, France, and to assess what kind of information can be retrieved concerning the underlying tectonic processes. We also combine this spatial structure of the landscape with denudation rates derived from cosmogenic nuclides in order to compare the relative spatial distribution of surface processes and uplift rate inferred from the DTM analysis with absolute values. In the following, we first present the main features of our working area, the Puimichel Plateau in the Mio-Pliocene Digne-Valensole basin, with a focus on key aspects of its main structures and history that make it an interesting place to investigate the interactions between tectonic forcings and surface processes response. Then, we introduce the morphological analysis and cosmogenic nuclide methods we used, and we describe the corresponding datasets produced during this study. Finally, we discuss the implications of these results in terms of the imprint of tectonic gradients into hillslope morphology, the constraints that can be put on tectonic structures at depth from high-resolution topographic data and their relationship with denudation rates calculated from cosmogenic nuclide concentrations.

2 Setting

The studied area is located in the southern part of the European Western Alps, where mountain ranges result from a combination of Pyrenean (Late Cretaceous to Eocene) and Alpine (Neogene) tectonic phases (Fig. 1). However, the present-day tectonic activity in the area is considered to be low, and there is no significant horizontal strain rate resolved from geodetic data, despite the occurrence of earthquakes along identified tectonic structures (Jouanne et al.2001; Walpersdorf et al.2018; Nocquet et al.2016). Paradoxically, leveling measurements indicate uplift rates up to 2.5 mm yr−1 in the northwestern Alps (Nocquet et al.2016). Focal plane mechanisms show that the inner Alps are characterized by extensional stresses, whereas the external Alps, including the studied area, are still under compression (Sue et al.1999; Delacou et al.2005). Within these compressional areas, plateau surfaces at 150–400 m above the present rivers suggest active uplift of the Western Alps due in part to flexural isostatic response to Quaternary erosion (Champagnac et al.2007, 2008) and in part to tectonic processes (Collina-Girard and Griboulard1990; Schwartz et al.2017).

Figure 1(a) Geological setting of the studied region in SE France (1∕106 BRGM geological map). Black square indicates the position of inset (b). (b) Focus on the Mio-Pliocene Valensole Plateau and the main regional tectonic structures. Focal mechanisms from Nicolas et al. (1990).

The flexural history of the Alps is particularly recorded in the Neogene basins at the front of the Western Alps (e.g., Beck et al.1998), such as the Digne-Valensole basin in the southern Alps (Mercier1979). The Digne-Valensole basin collected material eroded from the Alps and transported by the Durance, Bléone, Asse and Verdon rivers (Fig. 2). The basin is filled by marine deposits overlain by conglomerates of the continental Valensole Formation, interpreted as an alluvial fan system prograding southward (e.g., Clauzon et al.1989). Deposition starts with Aquitanian marine sandstone, followed by the continental conglomerate of the Valensole Formation, which is Serravallian to Tortonian in age at its base and up to the early Pleistocene at its top (Mercier1979; Clauzon et al.1989; Dubar1984a; Dubar et al.1998). The eastern edge of the basin is overthrusted by the Middle Miocene to late Quaternary Digne Nappe (Lickorish and Ford1998; Roure et al.1992; Gidon and Pairis1988; Hippolyte et al.2011), and it is bordered to the west by the NNE-trending Durance seismically active fault, a dextral fault with a reverse west-side-up component (Roure et al.1992; Cushing et al.2008). The Digne-Valensole basin is of particular interest for the study of the interaction between tectonics and surface processes because of its structural location and its late Pliocene to Quaternary infill that allows us to demonstrate the Quaternary activity of several faults within the basin, and along its borders (Clauzon1982; Ritz1992; Cushing et al.2008; Hippolyte and Dumont2000; Hippolyte et al.2011).

Figure 2(a) Location of the studied basins and sampling sites at the western edge of the Puimichel Plateau (Fig. 1). White line is the location of the profile used in Fig. 7. The numbers denote distance along the profile in kilometers and match the horizontal distance coordinate of Figs. 7 and 11. (b) Projection of the summit surface along the profile. Horizontal coordinates are consistent with labels of the cross section in panel (a). (c) Red contour lines of the two-way travel times (ms) to the top of Oxfordian black marls from a seismic survey synthesis included in the scientific report of drilling BSS002DWDJ available in the subsurface BSS (Banque du Sous-Sol) database of BRGM. Brown circle indicates the location of the drill site. Dark red dashed lines delineate possible geological structures. Dashed white line is parallel to the main cross section and used as an additional section line in Fig. 11. (d) Simplified geological section across the Puimichel Plateau adapted from Dubois and Curnelle (1978). Note that the orientation of this cross section is slightly different from the geomorphological transect studied here.

Figure 3(a) Panorama of the western edge of Puimichel Plateau and the studied basins, viewed from Ganagobie Abbey (see Fig. 1b for location). (b) Hillslope flank covered with regolith clasts. (c) Hilltop sampling site K (see Fig. 2a for location). (d) Hilltop sampling site M.


The Digne-Valensole basin is dissected by the Bléone and the Asse rivers, which divide the area into three parts: the Valensole Plateau to the south; the Puimichel Plateau in the middle, which is partly eroded and incised by the Rancure river; and the mountains of the Duye valley to the north. These rivers also dissected the Valensole basin during the Messinian crisis, when the Mediterranean sea level dropped by about 1500 m (Ryan and Cita1978; Hsü et al.1977; Clauzon1982; Clauzon et al.2011). In the Valensole basin, the Messinian paleocanyons are mostly buried under Pliocene and early Quaternary sediments. A few sections of these paleocanyons could be mapped near Oraison (Dubar et al.1998) and near Digne (Hippolyte et al.2011), where the Messinian erosional surface separates the Valensole I and II formations (Dubar1984a; Clauzon1996). The Valensole II formation filled the Messinian canyons and covered the central and southern part of the Valensole basin (Dubar1984a), with its top surface presently forming the Valensole and Puimichel plateaus (Fig. 3). The age of the top surface of these plateaus ranges between 0.7 Ma in the east, near the Digne Thrust (Dubar et al.1998), and 1.7 Ma in the west, along the Durance river (Dubar1984a, 2014; Dubar and Semah1986).

This surface was used by Champagnac et al. (2008) as a passively deformed marker to identify long-wavelength tilting of the Alpine foreland, in part as a response to erosional unloading. At shorter wavelengths, the exceptional preservation of this surface allowed Hippolyte and Dumont (2000) to demonstrate the Quaternary activity of the Lambruissier anticline (Fig. 1a) within the Digne-Valensole basin. This SW-verging fold generated an 80 m high and 5 km long morphological ridge above the Puimichel Plateau surface. Younger terraces have been mapped in the area (Dubar1984b), but due to active erosion and poor preservation of their surfaces it is not possible to use them as reliable benchmarks to measure finite deformation. To the northeast of the Lambruissier anticline, an older, late Neogene fold is stratigraphically overlain by horizontal late Pliocene deposits of the Bléone river Messinian canyon (Hippolyte et al.2011). These ages demonstrate the southward propagation of deformation within the Digne-Valensole basin, which would be further confirmed if the Mées ramp anticline (Dubar et al.1978), located south of the Lambruissier anticline, was active. This activity is suggested by (1) elevation anomalies in the Mindel terrace of the Durance river (Gabert1979; Dubar1984b), (2) the dip of the Puimichel Plateau surface which is more than 3 times higher (25 m km−1) than the dip of the Valensole Plateau (8 m km−1) (Hippolyte and Dumont2000) and (3) the striking asymmetric pattern of the drainage network that may have recorded a progressive tilt of the Puimichel Plateau (Collina-Girard and Griboulard1990; Hippolyte and Dumont2000). Recent tectonic deformation of this area is also in agreement with modeling of the thermal history of the northern tip of the Digne-Valensole basin (Schwartz et al.2017).

Our study is specifically focused on the western edge of the Puimichel Plateau, which is dissected by a series of small basins (Fig. 2) draining directly into the Durance river. The eastern limit of these basins corresponds to the post-Pliocene abandonment surface of the summit of the plateau. Seismic surveys and drilling have identified an important uplifted basement structure below this plateau (Mées Structure, Fig. 2d) of Upper Cretaceous to Eocene age, with possible Alpine Miocene reactivation (Dubois and Curnelle1978). The region is characterized by a Mediterranean climate with mean annual temperature (MAT) of 13 C and mean annual precipitation of 700 mm (data at Saint-Auban Météo-France weather station over the 1981–2010 period). The dominant lithology is the Mio-Pliocene Valensole conglomerate, with an age either pre- or post-Messinian depending on the geometry of erosional surface, which has not been continuously mapped in this area. However, bedrock geology is uniform, with mostly clast-supported conglomerates, which weather primarily by destruction of the sandy matrix over a few tens of centimeters below the surface. The clasts (5 to 10 cm in size) are set loose but remain interlocked with little vertical movement and mixing inside the regolith profile. Once the clasts reach the surface, they are free to move downslope.

For the hillslope domain, the existence of this extensive mobile regolith cover implies that hillslopes are mostly under transport-limited regime, which is an important requirement of the topographic analysis described below (Fig. 3). Concerning the fluvial network, we did not observe any large-scale alluviation but only a thin cover of sediment with many occurrences of outcropping unweathered conglomerate bedrock, which leads us to consider that these streams are mostly under detachment-limited dynamics.

3 Methods

3.1 Topographic analysis

The topographic analysis carried out in this study relies on a 1 m resolution airborne lidar digital terrain model (DTM) acquired in 2014 as part of the RGE ALTI® database from IGN (Institut National de l'information Géographique et Forestière) and covering our study area, as well as most of the Valensole and Puimichel plateaus. The core of our analysis consists of the extraction of high-resolution hillslope and channels topographic metrics along a transect located at the northwestern edge of the plateau in order to identify short-wavelength variations in the distribution of surface processes (Fig. 2), as opposed to the long-wavelength deformation investigated by previous studies (Champagnac et al.2008).

3.1.1 Hillslope morphology

We present here the theoretical background supporting the interpretation of hillslope-scale morphological parameters. Mass conservation across a steady-state 1-D hillslope profile can be expressed as

(1) q s x = β E ,

where x is the horizontal coordinate ([L]), qs is the sediment flux ([L2T−1]), β is the rock-to-regolith density ratio (dimensionless), and E is the erosion rate ([LT−1]) which is equal to the rock uplift rate under steady-state conditions. Equation (1) can be combined with a geomorphic transport law (GTL; Dietrich et al.2003), describing sediment flux qs over a hillslope as a nonlinear function of local hillslope gradient S=z/x (Roering et al.1999, 2007),

(2) q s = D S 1 - S S c 2 ,

where D is a diffusion coefficient ([L2T−1]) and Sc is a critical hillslope gradient (Roering et al.1999). For gentle slope areas, such as in the vicinity of hilltops, Eq. (2) can be linearly approximated as qs=DS. Combining this expression of qs with Eq. (1) yields a linear relationship between the erosion rate E and the second spatial derivative of topography or hilltop curvature CHT,

(3) E = D C HT β .

Equation (3) will be central in the interpretation of our results, as it allows us to combine the two types of data acquired at hilltop sites during this study: high-resolution morphometric measurements (CHT) and denudation rates (E) derived from cosmogenic nuclide concentrations.

Combining Eqs. (1) and (2) and integrating yields the steady-state elevation profile z associated with a spatially uniform erosion rate E (Roering et al.2001),


With LH as the hillslope length (horizontal distance from hilltop to channel), a reference erosion rate (Roering et al.2007) can be defined as

(5) E R = D S c 2 L H β ,

and similarly a reference relief RR=ScLH which represents the maximum hillslope relief can also be defined.

These two reference values allow us to normalize hillslope relief R and erosion rate E into their nondimensional equivalents as

(6) R = R S c L H ,

and, dividing Eq. (3) by Eq. (5),

(7) E = 2 C HT L H S c .

Finally, Eq. (4) can be expressed in nondimensional form:

(8) R = 1 E ( 1 + ( E ) 2 - ln 1 2 1 + 1 + ( E ) 2 - 1 ) .

The purpose of our analysis of the DTM is to extract these various metrics characterizing the relief structure and erosion of hillslopes, and in particular the hilltop curvature CHT, as well as hillslope relief R and length LH. These measurements will then allow us to compute a nondimensional erosion rate E according to Eq. (7). We use the approach presented by Hurst et al. (2012) and Grieve et al. (2016b), which we implemented into the GRASS GIS environment (Neteler et al.2012) and R scripting language (R Core Team2018), as described in Godard et al. (2019) (Fig. 4).

Figure 4(a) Hillshade image from a 1 m IGN RGE ALTI® lidar digital terrain model (DTM). Thick brown lines indicate the hilltops extracted from the DTM. Light blue polygons are the floodplains extracted using the approach by Clubb et al. (2017). Thin purple lines are flow lines routed from the hilltop toward the floodplain. (b) Corresponding orthophotography (IGN BD ORTHO®).

First, the DTM was filtered to remove short-wavelength noise using the despeckling algorithm of Sun et al. (2007) and its application to DTM filtering as described in Stevenson et al. (2010). The river network was extracted using a geometric approach by defining a 0.2 m−1 contour curvature threshold for the definition of channel heads (e.g., Pelletier2013; Clubb et al.2014). The narrow floodplains that are present in our working area were delineated following the approach proposed by Clubb et al. (2017). Hilltops were then identified as the intersecting margins of basins over all ranges of stream orders, and curvature was computed at every hilltop pixel by fitting a quadratic surface over a 30 m wide window, which is large enough to filter out short-wavelength surface roughness and small enough to avoid perturbation from the ridge-and-valley topographic signal (Hurst et al.2012; Godard et al.2016; Grieve et al.2016c; Lashermes et al.2007). Flow was routed downslope from hilltop pixels to the edge of floodplains using the algorithm of Mitasova et al. (1996), and the resulting flow lines were used to compute hillslope relief and length (Hurst et al.2012; Grieve et al.2016a, b). Flow lines and associated data were grouped into patches of a least 50 contiguous hilltop pixels (Grieve et al.2016b).

3.1.2 Channel morphology

We also extracted standard metrics from the river profiles of the studied catchments (Fig. 5). River incision I ([LT−1]) can be parameterized as a function of along-channel topographic gradient S and drainage area A ([L2]) as

(9) I = K A m S n ,

where K is an erodibility coefficient ([L1−2mT−1]), and m and n are empirical exponents (Howard1994; Whipple and Tucker1999). Under steady-state conditions, river incision I equals rock uplift U and Eq. (9) can be reorganized as

(10) S = d z d x = U K 1 n A - m n .

Parameters ks=UK1n and θ=mn are referred to as the steepness index and channel concavity (Wobus et al.2006; Kirby and Whipple2012), and they are often determined by regression in a slope–area diagram. Steepness indexes are of particular interest in tectonic studies, due to their direct dependence upon the rock uplift U, and, under the assumption of constant erodibility K, they can be used to decipher relative spatial variation in U (e.g., Kirby et al.2003). In order to allow for a meaningful comparison between channels of different concavities, a reference mn value is chosen and used in the regression in order to obtain a normalized steepness index ksn.

Figure 5(a) River profile for the main trunk and tributaries of the Moureisse catchment (see Fig. 2 for location). (b) χ transform of the longitudinal profile (Perron and Royden2012), using m/n=0.23. Inset shows the evolution of the R2 of the linear regression of elevation against χ for a range of mn values.


Calculating S by differentiating the river long-profile often yields noisy data. Following Perron and Royden (2012), Eq. (10) can be integrated from base level, at position xb, to x as

(11) z ( x ) - z ( x b ) = x b x U K A m 1 n d x .

Under the assumption that U and K are spatially constant, Eq. (11) becomes

(12) z ( x ) = z ( x b ) + U K A 0 m 1 n χ ,

where A0 is a reference drainage area, with

(13) χ = x b x A 0 A m n d x .

Equation (12) implies that in a z vs. χ diagram a steady-state channel profile should plot a straight line, with a slope proportional to steepness index.

For every investigated basin, we searched the mn value (concavity) leading to the best linearization of the χ-transformed river profile (Perron and Royden2012). We fixed the reference value according to the mean of observed mn value across all basins. We then again χ-transformed the river profiles using this reference mn to compute normalized steepness indexes.

3.2 Cosmogenic nuclides

The topographic analysis methods described in the previous section allow us to identify relative spatial patterns for the intensity of surface processes such as surface denudation, which can be used to infer rock uplift distribution. We used in situ-produced cosmogenic nuclide concentration measurements to constrain the absolute denudation rate values. Active fluvial sediments are present along the stream network of the studied catchments (Fig. 2a), and they could be sampled to derive basin averaged denudation rates (von Blanckenburg2005). However, most of the surveyed channels displayed complex dynamics with localized occurrences of aggradation, splitting of the main channel and local colluvial inputs, such that the required hypothesis of a homogeneous contribution of the whole upstream area was most likely invalid. For that reason, we rather sampled material directly at the hilltops, which allows for a clear identification of the origin of the sampled material and a straightforward comparison with the topographic metrics extracted from the high-resolution DTM. As noted earlier, weathering principally affects the sandy matrix of the conglomerate, liberating the clasts which remain interlocked until they reach the surface, such that vertical mixing within the regolith is minimal at the hilltop.

Samples for 10Be and 26Al concentration measurements were collected at 10 hilltop sites, by amalgamating 30 to 40 individual sandstone clasts derived from the bedrock conglomerate (Fig. 2a). Samples were crushed and sieved to extract the 250–1000 µm fraction, which was submitted to sequential magnetic separation. The remaining fraction was leached with 37 % HCl to remove carbonate fragments. The samples were then repetitively leached with H2SiF6 and submitted to vigorous mechanical shaking until pure quartz was obtained. Decontamination from atmospheric 10Be was achieved by a series of three successive leaching processes in concentrated HF, each removing 10 % of the remaining sample mass (Brown et al.1991). After addition of an in-house 9Be carrier, the samples were digested in concentrated HF. Be and Al were isolated for measurements using ion-exchange chromatography. 10Be∕9Be and 26Al∕27Al measurements were performed at the French AMS National Facility, located at CEREGE in Aix-en-Provence. All results and technical characteristics of the measurements and calculations are provided in Tables 1 and 2.

Table 1Cosmogenic nuclide 10Be results.

a Dissolved pure quartz mass. b In-house carrier mass,  150 µL at 3.025 × 10−3  g g−1 (Merchel et al.2008). c 10Be∕9Be ratios were calibrated against the National Institute of Standards and Technology standard reference material 4325 by using an assigned value of 2.79 ± 0.03×10−11 (Nishiizumi et al.2007). d Uncertainties are reported at the 1σ level. e Uncertainties on isotopic ratios are calculated according to the standard error propagation method using the quadratic sum of the relative errors and include a conservative 0.5 % external machine uncertainty (Arnold et al.2010), the uncertainty on the certified standard ratio, a 1σ uncertainty associated with the mean of the standard ratio measurements during the measurement cycles, a 1σ statistical error on counted events and the uncertainty associated with the chemical and analytical blank corrections. f Two process blanks were treated and measured with our samples, yielding 10Be∕9Be ratios of 1.47 ± 0.27 and 0.90 ± 0.29 × 10−15. It corresponds to an upper 1σ bound of 55 and 38 × 103 10Be atoms for the background level in these two blanks, which is at least 20 times lower than the number of 10Be atoms in the dissolved sample masses (30 times lower on average over our dataset). g Denudation rates were then calculated with the online calculator described in Balco et al. (2008) and the nuclide-specific LSD (Lifton–Sato–Dunai) scaling scheme (Lifton et al.2014), using the CRONUS-Earth calibration dataset (Borchers et al.2016) for the calculation of spallation production rates and muon production rates according to Balco (2017). We use 160 g cm−2 for the effective attenuation length for spallation in rock and a density of 2.5 g cm−3. No shielding correction was considered for the hilltop sites we sampled, which were selected on nearly horizontal ridgelines.

Download Print Version | Download XLSX

Table 2Cosmogenic nuclide 26Al results and comparison with 10Be results from Table 1.

a Naturally occurring Al measured by ICP-OES (inductively coupled plasma optical emission spectrometry). No Al carrier solution was added to the samples. b The measured 26Al∕27Al ratios were normalized to the in-house standard SM-Al-11 whose ratio of 7.401 ± 0.064 × 10−12 (Arnold et al.2010) has been cross-calibrated against primary standards from a round-robin exercise (Merchel and Bremser2004). c Uncertainties are reported at the 1σ level. d An analytical blank yielded a ratio of 7.27 ± 2.17 ×10−15. e Procedures for that calculation of denudation rates are identical to the ones described in Table 1. f Integration timescales for 10Be and 26Al denudation rates (von Blanckenburg2005).

Download Print Version | Download XLSX

3.3 Surface deformation modeling

Under the assumption that the spatial distribution in erosion rates results from tectonic forcing, observations of absolute or relative variations in the intensity of surface processes are commonly used to derive information on rock uplift patterns and the geometry of associated structures (Wobus et al.2006; Scherler et al.2014; Le Roux-Mallouf et al.2015). While the key hypothesis of a tectonic origin for the variability in erosion proxies will need to be discussed in detail, we present here a modeling approach which can be used to interpret surface erosion patterns, which are considered proxies for rock uplift, in terms of the associated structures at depth.

We use a simple elastic dislocation model (Okada1985) to predict surface deformation distributions and compare them with our observed E evolution along the investigated profile (Fig. 2). This type of modeling approach is usually associated with the study of upper crustal deformation during the seismic cycle, but it has also been successfully applied to interpret surface deformation patterns associated with blind thrusts over longer timescales (Ward and Valensise1994; Benedetti et al.2000; Myers et al.2003). As noted by Myers et al. (2003), who used elastic dislocation modeling to study folding in the Los Angeles Basin, it allows us to keep track, at first order, of the displacement of material associated with the activation of the fault, independently of the mechanical parameters. We consider a single planar dislocation embedded in an elastic medium. We define its geometry with four parameters: the horizontal position and depth of its upper limit (varying from 0 to 20 km along the profile and from 0 to 10 below the surface, respectively), dip angle (varying from 0 to 90) and length (varying from 0 to 4 km). This set of four parameters allows us to predict a surface deformation profile which is compared to the observed E pattern. We do not fit the absolute value of E, which has no significance with respect to the elastic dislocation model, but rather we fit the wavelength of the predicted deformation and simply scale its amplitude to that of the observed E profile.

We use a standard Monte Carlo Markov Chain (MCMC) approach to move through this parameter space and estimate the posterior distributions (Metropolis et al.1953). The consecutive displacements during the sampling procedure were driven by the Metropolis–Hastings algorithm, with an acceptance rate of 20 %. We ran 16 independent chains, each 106 in length with a 104 length burn-in phase. The multivariate potential scale reduction factor is 1.004, suggesting that convergence was achieved (Gelman and Rubin1992).

4 Data

In this section we present the primary data acquired during this study and the associated direct observations. The interpretations built on these datasets are developed and discussed in the next section.

4.1 Topographic analysis

Most hillslope profiles display a clear convexity and progressive downward increase in slope, with almost linear portions at the bottom of the hillslope and gradients close to 0.6 m m−1 (Fig. 6a). Calculation of average hillslope gradient shows that most hillslopes have average gradients below 0.6 m m−1 (Fig. 6b). Following Hurst et al. (2019), we determined the value of Sc for which 99 % of hillslopes have a relief inferior to the maximum value (LHSc) and obtained Sc=0.52m m−1. In the following we use Sc=0.6m m−1 for the critical gradient value (Fig. 6c). This Sc value is lower than what was found in other settings (Grieve et al.2016b), but it is close to the natural angle of repose in many granular materials ( 30). We note that two particularities of our study area are the nature of the regolith, which is mostly constituted of highly mobile conglomerate-derived clasts moving downslope by both creep and dry ravel, and the isolated vegetation providing little cohesion (Fig. 3b), as opposed to the finer-grained soils supporting denser root networks observed in many other similar studies. Most topographic metrics extracted at the hillslope scale display important variations along the studied profile, both in terms of basin averaged or binned values. Nondimensional erosion rate E increases 2-fold from south to north (Fig. 7a). This variation is not evenly distributed along the profile but occurs over less than 4 km of horizontal distance. The hilltop curvature CHT pattern closely mimics that of E, increasing from 0.01 to 0.02 m−1 (Fig. 7b). For both metrics, basin averaged values are highly consistent from one catchment to the other and delineate a clear trend along the section, with the exception of one outlying small catchment at  4 km distance. The evolution of hillslope relief is less pronounced but also shows an increase from  40 to  60 m (Fig. 7c). On the contrary, no clear systematic changes in hillslope length can be observed along the profile, with values ranging from 140 to 160 m (Fig. 7c).

Figure 6(a) Selected hillslope profiles from the studied basins (Fig. 2). Grey lines indicate topographic gradient values. (b) Distributions of hillslope average topographic gradient for all the studied basins (solid line) and five northernmost basins (dashed line) where relief R, hilltop curvature CHT and nondimensional erosion rate E are the highest (Fig. 7). (c) Joint distribution of hillslope relief and length for the flow lines extracted from the studied basins (Fig. 4). Grey lines indicate topographic gradient values.


Figure 7Projection of hillslope and fluvial parameters along the profile of Fig. 2. The horizontal distance values used here match the distance measured along the profile in Fig. 2. For all panels, open symbols refer to basin averaged values (location of basins in Fig. 2) and closed symbols to 1 km length bin averages along the profile (error bars are ±1σ). (a) Evolution of nondimensional erosion rate calculated as E=2LHCHT/Sc (Roering et al.2007). (b) Hilltop curvature computed over 15 m radius window. (c) Hillslope length and relief from the flow-line patches (Fig. 4). (d) Normalized steepness index and mn ratio (concavity, reference value of 0.25) extracted from channel profiles of the studied basins.


For all studied basins, the river profiles display a regular concave-up shape, and in most situations χ-transformed main trunk profiles, as well as the main tributaries, collapse along a linear trend, suggesting the absence of a major transient perturbation propagating through the river network (Fig. 5). Small tributaries usually show higher dispersion due to changes in processes in small colluvial valleys. Usual topographic indexes, such as mn ratio and normalized steepness index (ksn), were extracted from the fluvial network. The mn ratio ranges from 0.2 to 0.4, with an average of 0.24 ± 0.06, and a reference value of 0.25 was used in the following analysis. While this mn value is lower than the often reported 0.4 to 0.5 ratio, it is within the range of observations from Harel et al. (2016) for high erodibility lithologies, such as the setting we consider here. These fluvial metrics display a larger amount of scatter from one basin to another when compared with the patterns extracted from the hillslope morphology analysis (Fig. 7d). However, it can be noted that the four northernmost basins display higher ksn values than the rest of the section, and that basin averaged E and ksn are significantly positively correlated (Fig. 8).

Figure 8Evolution of basin averaged E as a function of normalized steepness index (ksn). Circles are colored according to the position of the corresponding basins along the transect (Figs. 2 and 7). Circle radius is a function of catchment size (ranging from 0.5 to 4 km2). Solid and dashed lines correspond to a linear fit and its 95 % confidence interval (R2=0.62 and p<10-3).


4.2 Cosmogenic nuclide data

Measured concentrations in our hilltop samples range from 43 to 117×103atoms g−1 and 281 to 867×103atoms g−1 for 10Be and 26Al, respectively (Tables 1 and 2). The corresponding 26Al∕10Be ratios vary from 6.5 ± 0.8 to 8.2 ± 1.2. While some of these ratios are slightly higher than the theoretical value, confidence ellipses are always overlapping the steady-state denudation curve (Fig. 9). These concentrations correspond to denudation rates ranging from 42 to 115 mm kyr−1 and from 38 to 114 mm kyr−1 for 10Be and 26Al, respectively (Tables 1 and 2). The denudation rates calculated from both nuclide concentrations are consistent but display a small deviation from the 1:1 line, with 26Al denudation rates slightly lower than their 10Be equivalents. The observed 26Al∕10Be values argue against a significant contribution of an inherited CRN inventory from the history of the clasts prior to their deposition inside the Valensole conglomerate, in particular if they derived from the Valensole II formation. In the following we only consider 10Be data for our analysis of the relationship between high-resolution topography and denudation rates, due to their lower uncertainty.

Figure 9(a) Two nuclides plot for the sampled sites, with 68 % confidence ellipses (see Fig. 2 for location and Tables 1 and 2 for data). Solid and dashed brown lines indicate the predicted 26Al∕10Be ratio for steady-state denudation and constant exposure histories, respectively. (b) Comparison of denudation rates calculated from measured 10Be and 26Al concentrations using the LSD scaling scheme (Lifton et al.2014) and calculation procedure from Balco et al. (2008). Blue line and envelope are linear fit and its 95 % confidence interval, respectively.


5 Discussion

We now discuss the morphological and geochronological data presented in the previous section in terms of the main controls on hillslope morphology, the relationship between hillslope and channel properties, the geometry of the underlying tectonic structures, and the comparison between morphological observations and denudation rates at the same sites.

5.1 Interpretation of the spatial variability in hillslope morphology

5.1.1 Possible controls on hillslope morphology

We observe a pronounced and systematic variation of hillslope morphology along the studied transect, with hillslope curvature CHT undergoing a 2-fold increase from S to N (between 7 and 10 km in Fig. 7b). We evaluate below the possible controls on this evolution. We first note that our study area extends over  10 km in length, with catchment average elevation and relief ranging from 460 to 620 m and from 100 to 300 m, respectively. These limited changes in elevation imply that climatic conditions can be considered constant in terms of mean annual precipitation and temperature (MAP and MAT), and they cannot account for the observed variations in hillslope morphology. Vegetation cover is also homogeneous over the western flank of the Puimichel Plateau, with a forest dominated by Quercus pubescens and occurrences of Quercus ilex and Pinus sylvestris. Similarly, the investigated basins are eroding into Mio-Pliocene conglomerates with no major changes in the nature and properties of the bedrock or regolith material. We note that this homogeneity of geological, climatic and biological properties over the transect is a specificity of our studied area, and this might not be warranted in other settings, where eventual disparities in these properties might complicate the interpretation of fluvial and hillslopes morphologies.

Hillslopes have been shown to record transient waves of erosion propagating through landscapes (Hurst et al.2012, 2013a; Mudd2017). The series of studied basins are directly connected to the Durance river base level, and the eventual propagation of incision pulses and along-stream knickpoints might impact the network of tributaries and adjacent hillslopes, inducing a differential hillslope response (Hurst et al.2012). However, several lines of evidence argue against such control on the observed distribution of hillslope properties. First, no major knickpoints have been identified along the Durance river in the vicinity of our working area. Second, we note that the pattern of evolution for E, CHT and R along the transect is characterized by almost constant values for 6 km, followed by a rapid increase over less than 4 km, rather than the progressive variations that would be expected to result from the propagation of a knickpoint in front of the western edge of the plateau. At last, we note that the propagation of an incision wave would result in a pattern with higher erosion areas in the southern part connected to the adjusted or adjusting landscape downstream of the propagating incision pulse, and slower erosion in the northern yet unaffected areas, which is exactly contrary to what we observe here. Therefore, we propose that the observed pattern is unlikely to result from the transient adjustment of the landscape to the propagation of a wave of incision along the Durance river.

Another possibility to generate the observed distribution of hillslope parameters would be a sustained differential rock-uplift pattern associated with recent or ongoing deformation. Several studies have already pointed at geomorphic evidence for recent tectonic activity in the northern part of the Puimichel Plateau (Collina-Girard and Griboulard1990; Hippolyte and Dumont2000), which is confirmed by recent seismic activity (Nicolas et al.1990), with several magnitude-4 thrust-slip events with E–W strike directions (Fig. 1b). Notably, a major inflection of the abandonment surface occurs at 7–8 km of horizontal distance along the transect (Fig. 2b), suggesting post-Pliocene uplift of the northern part. The location of this surface inflection coincides with the transition area for E, CHT and R (Fig. 7). Such spatial coherence between a long-term finite deformation pattern passively recorded at 1 Myr timescale by the summit surface, and the shorter-term active erosional response of the landscape observed through hillslope morphometric indices, argue for a control by differential uplift rates across the transect.

Additionally, we note that this transition zone is also coincident with the southern flank of a major basement uplift, identified by seismic surveys, located below the northern part of the Puimichel Plateau (Fig. 2c and d). This uplifted basement is a long-lasting structure formed during the Pyrenean orogeny and is associated with basement thrusts on its southern edge, at the location of the observed geomorphic transition (Dubois and Curnelle1978). Quaternary reactivation of such faults has been invoked to explain the surface deformation pattern of the plateau farther to the north (Hippolyte and Dumont2000), and we propose a dominant tectonic control for the distribution of proxies for surface denudation along our transect (Fig. 7).

The R vs. E relationship (with Sc=0.6m m−1) shows that the studied catchments plot slightly below the line denoting steady state predicted by Eq. (8) (Fig. 10), indicating possible decaying dynamics of the landscape toward lower relief, likely due to a change in the climatic or tectonic boundary condition (Mudd2017). All catchments appear to be similarly affected, with no obvious gradient between the northern and southern ones. For that reason, this decay is not likely to result from a decrease in the amount of differential rock uplift across the transect, but it could be associated with a regional change in the intensity of the top-down forcing of climatic origin, which would modify the value of the diffusion coefficient D. However, we note that using Sc=0.5m m−1, which appears to be reasonable for the vast majority of the studied hillslopes (Fig. 6), reduces considerably the deviation from the line denoting steady state.

Figure 10Basin averaged R vs. E plot (Roering et al.2007; Grieve et al.2016b). Two values for the critical hillslope gradient Sc are tested. Open circles correspond to Sc=0.6m m−1, colored according to the position of the basins along the transect (Figs. 2 and 7). Small dark filled circles correspond to Sc=0.5m m−1. Pale yellow symbols are averages over hilltop patches, for Sc=0.6m m−1 (see text for details). Thick grey line corresponds to the steady-state relationship between R and E predicted by Eq. (8).


5.1.2 Constraints on tectonic structures from hillslope morphology

As discussed in the previous section, the hypothesis that the evolution of hillslope morphology along the studied transect (Fig. 7) results from a spatial variation of rock uplift is supported by several lines of evidence and notably the presence of a coincident major warping of the Pliocene abandonment surface of the plateau. Here we develop further the interpretation of this surface deformation pattern in order to put constraints on the geometry of the associated structures. We use the dislocation modeling approach presented above to constrain the geometry of a fault whose displacement could explain the observed change in E along the transect (Fig. 7) through an inversion procedure. The parameters characterizing this fault geometry are the horizontal position and depth of the upper tip of the fault, its dip angle, and its length (Fig. 11).

Figure 11(a) Nondimensional erosion rate evolution from hilltop patches (mean and standard deviation binned every 1 km). Red curve is the result of the optimization of a simple dislocation model (Okada1985), with amplitude adjusted to the range of E values. See text for details. (b) Projection of the deformed summit surface of the Puimichel Plateau (Fig. 2b). (c) Two-way travel times to the top of Oxfordian black marls, interpolated from seismic surveys across the studied area, as indicated in the report of drilling BSS002DWDJ available in the subsurface BSS database of BRGM. Data are projected onto the section indicated as a dashed white line in Fig. 2c. (d) Geometry of the dislocation (black line) used to compute the surface deflection on panel (a) (red curve). Dashed lines (position and depth of the fault upper end) and light red surface (fault dip angle) correspond to 68 % density intervals from the marginal distributions of panel (e). Thin dark lines indicate the limits of the geological units from the cross section presented in Fig. 2d (Dubois and Curnelle1978). (e) Marginal distributions from a MCMC exploration for the parameters of the elastic dislocation model.


We observe that the horizontal position of the top of the fault is the best constrained parameter with a most likely value around 8 km (same horizontal reference frame as the profile in Fig. 7). The most likely depth for the upper limit of the dislocation is around 2 km, and the suggested dip of the structure is > 50. The length of the dislocation is not well constrained by our inversion.

Interestingly, the suggested horizontal position for the upper tip of the dislocation is close to the major deflection of the Pliocene abandonment surface, imaged by the lidar data (Fig. 11b). The high-E northern part of the studied transect roughly corresponds to the Mées Structure identified from seismic surveys and drilling (Dubois and Curnelle1978), which is a large anticline inherited from the Late Cretaceous–Eocene compression of the Pyrenean phase (Fig. 2b). We note that the suggested horizontal position of the dislocation is also coincident with the southern flank of the structure and the complex of thrusts responsible for the folding. Some of these thrusts correspond to reactivated high-angle basement structures, which is in agreement with the inferred dip angle of the dislocation. The depth of at least 2 km for the top of the dislocation is also consistent with the observation that the Mio-Pliocene reflectors do not display major offsets in the available seismic data (Dubois and Curnelle1978), and that the Cenozoic formations underwent long-wavelength folding rather than localized faulting. Overall, the geometry of the structure constrained by our simple model is compatible with the activity of steep south-verging inherited structures, affecting the basement and Mesozoic series, and the Quaternary reactivation which induced a long-wavelength warping of the Mio-Pliocene cover and differential uplift along our transect. We note that several other recently active structures have been documented farther to the north, corresponding to similarly oriented south-verging thrust-and-fold systems, such as the Lambrussier anticline which affects the northern edge of the Puimichel Plateau (Hippolyte and Dumont2000). The amount of finite deformation accommodated by these folds progressively decreases southward, such that the structure we identified could correspond to the most recently activated as an in-sequence system. Finally, while high-resolution hillslope morphology analysis has already been used to constrain rock uplift patterns in a limited number of studies (Hurst et al.2013a, 2019; Clubb et al.2020), our results are the first to illustrate the use of such data to infer the geometry of tectonic structures at depth.

5.2 Hillslopes and channel dynamics

Interestingly, parameters extracted from the analysis of river long-profiles such as steepness indexes, which are commonly used to decipher tectonic patterns in erosional landscapes (Kirby and Whipple2012), do not display as clear a pattern as hillslope metrics (Fig. 7d). Normalized steepness index values (ksn) are on average higher in the northern part of the transect with respect to the southern part and positively correlated with E (Fig. 8), but the data are scattered and we do not observe a clear progressive increase comparable to what is displayed by E or CHT. In each of the studied basins, the main trunk and its tributaries display regular concave-up profiles, which collapse along a single trend in χ plots (Fig. 5). This observation suggests that, at the scale of each catchment, the river network is globally equilibrated with respect to a common rock uplift rate (Perron and Royden2012), and argues against the impact of local perturbation along the river profile as an origin for the observed scatter. An underlying assumption of our river profile analysis is that the streams behave as purely detachment-limited systems. While our field survey did not allow us to identify thick alluvial cover along the stream network, local and intermittent shifts toward transport-limited behavior could explain the apparent subdued response of ksn across the inferred rock uplift gradient. We note, however, that in such a situation, stream concavity (mn) should display some sensitivity to changes in rock uplift (Wickert and Schildgen2019), which is not observed here (Fig. 7d).

A key difference between the two approaches is the very high spatial density of the metrics extracted for the hillslope dataset, which is several orders of magnitude denser than the evaluation of steepness index and concavity at 18 basins. This comparison illustrates the resolving power of high-resolution hillslope morphology analysis, which allows us to document short-wavelength patterns of erosion and uplift that are undersampled by scarcely distributed fluvial metrics. We note that Hurst et al. (2019) observed a clear response of both hillslope and channel metrics across the tectonic gradient they studied in the vicinity of the San Andreas Fault. Two important differences compared to our study are the existence of transient channel adjustment along the Bolinas Ridge and the dimension of the section, which is  30 km long in their case compared to the  10 km of our profile. Combined with a longer wavelength of the underlying tectonic signal, this latter difference might be a reason for the better sampling through fluvial metrics by Hurst et al. (2019), and the clearer relationship they observe with hillslope properties.

One prominent feature of the hillslope evolution across the rock-uplift gradient is the lack of significant changes in hillslope length, contrasting with the other extracted metrics (Fig. 7). Hillslope length, LH, remains nearly constant across the transect between 140 and 160 m, whereas hilltop curvature, CHT, which can be considered a proxy for erosion and rock uplift under a steady-state assumption, undergoes a nearly 2-fold increase. There is no significant inverse correlation between LH and CHT as observed by Hurst et al. (2013b). The characteristic horizontal length scale of landscapes, which can be evaluated through different types of measurements such as drainage density, spacing of first-order valleys or hillslope length, has been shown to be highly sensitive to external tectonic and climatic forcings, as well as internal parameters controlling erosion processes (e.g., Perron et al.2008a; Pelletier et al.2016; Clubb et al.2016; Hurst et al.2019). For example, Clubb et al. (2016) studied in detail the relationship between denudation rate and drainage density with analytical and numerical models as well as high-resolution topographic and cosmogenic nuclide data. They show a sensitivity of drainage density to erosion rates, which is very pronounced in the 50–100 mm kyr−1 range corresponding to our CRN data (Fig. 9) and conflicts with the observation of a nearly constant LH along our transect. However, this relationship is highly dependent on the parameters used for the fluvial and hillslope erosion laws, as shown by the formulation of the landscape Péclet number proposed by Perron et al. (2008b, 2009) as the ratio between characteristic fluvial and hillslope timescales:

(14) P e = K l 2 ( m + 1 ) - n D ζ 1 - n ,

where K is the erodibility parameter for fluvial incision; D the hillslope diffusion coefficient; m and n are the area and slope exponents of the fluvial incision law, respectively; and l and ζ are horizontal and vertical length scales for the landscape. In the special case where the slope exponent n of the stream power formulation for river incision is equal to 1, this Péclet number becomes independent from relief ζ, and hence uplift rate. Considering Pe=1 allows us to retrieve lc as a characteristic length scale for hillslope or channel transition, which again, in the n=1 case, does not depend on erosion or uplift rates. Therefore, the stability of LH across a two-fold erosion rate gradient, observed in our dataset, would hint toward a value of n close to unity. While the ratio mn can usually be constrained from the measured concavity of river profiles, the absolute values of the slope and area exponents of the stream power description for fluvial incision are debated, with many pieces of evidence pointing to n>1 (DiBiase and Whipple2011; Lague2014; Harel et al.2016). However, it is noteworthy that values closer to unity have often been reported for high-erodibility sediments (Harel et al.2016) or small catchments affected by colluvial processes (Lague and Davy2003), which are both notable characteristics of the area we investigate. In any case, our study clearly illustrates the potential of high-resolution hillslope morphological properties to resolve short-wavelength variations in rock uplift that are difficult to capture from the conventional methods based on fluvial profile analysis.

5.3 Comparison of cosmogenic nuclide data with high-resolution topography metrics

Our understanding of landscape dynamics relies on the formulation of geomorphic transport laws (GTLs), such as Eq. (2), as the foundation of landscape evolution models (Dietrich et al.2003). In most situations such models can be expressed as a differential equation involving spatial and temporal derivatives of the topographic surface elevation. More precisely, these equations will often relate the spatial structure of the landscape involving topographic slope or curvature with its rate of evolution as, for example, hillslope denudation rate or fluvial incision. Evaluating the relevance of such models requires obtaining actual measurements for these spatial and temporal descriptions of the landscape. Spatial properties of landscapes are usually derived from the analysis of DEMs at various resolutions, from which topographic gradient and curvature can be computed. On the other hand, geochronology techniques, such as cosmogenic nuclide concentration measurements in bedrock or sediment samples, allow for constraining the rate of lowering of the topographic surface though time and provide the framework to evaluate the temporal component of the landscape evolution problem.

For example, the comparison of catchment-wide denudation rates (CWDRs) calculated from measured 10Be concentrations in river sediments with topographic gradient extracted from digital elevation models has provided critical tests of GTLs for hillslope sediment flux (Ouimet et al.2009). While not often explicitly formulated in terms of the evaluation of a GTL, this kind of connection between some spatial property of landscapes (slope, relief, steepness index, etc.) and their rates of evolution is now standard in CWDR studies. However, it is to be noted that the interpretation of CRN concentrations in river sediments in terms of spatially averaged denudation rates suffers from important limitations resulting from the heterogeneity and stochasticity of erosion processes in space and time (e.g., Yanites et al.2009). Furthermore, the overwhelming majority of these studies rely on the evaluation of topographic metrics derived from medium-resolution DEMs (pixel size > 10 m) for which the computation of the spatial derivatives controlling erosion rates and sediment fluxes is inaccurate at the scale of hillslopes. Only a handful of studies have actually attempted to reconcile CRN-based denudation data with topographic metrics extracted from high-resolution DTMs into a GTL-based physical framework (e.g., DiBiase et al.2012; Hurst et al.2012, 2013b; Godard et al.2016, 2019; Neely et al.2019).

Figure 12(a) 10Be denudation rate against hilltop curvature. Dashed lines correspond to Eq. (3) for different values of the diffusion coefficient D. (b) Computed diffusion coefficient according to Eq. (3) at the various sites and corresponding probability density function. (c) Orange bars indicate the evolution of denudation rates along the profile from Fig. 7, calculated by sampling the distribution from the previous inset and applying Eq. (3) to individual hillslope patches. White squares are median values and orange bars indicate the interquartile range. Dashed horizontal lines delineate the  20 mm kyr−1 differential denudation rate across the transect. For comparison, light blue circles indicate the denudation rates actually measured at hilltop sites (excluding site P). These rates are systematically higher than those computed from Eq. (3) and the evolution of CHT in Fig. 7b, because sampling sites where selected on relatively high-curvature ridges, higher than 0.02 m−1 in most cases, whereas spatially averaged values along the transect are on average lower than 0.02 m−1 (Fig. 7b).


Our dataset allows for carrying out such a comparison of CRN denudation rates determined at hillslope sites with high-resolution hillslope topographic properties. In particular, we test the consistency of our dataset with prediction of simple hillslope diffusion formulations such as Eq. (3) relating hilltop curvature CHT with denudation rate. We observe that no single diffusion coefficient can explain the distribution of our data, but that, at most sites, the values of CHT and erosion rates are compatible with D in the 0.003 to 0.006 m2 a−1 range (Fig. 12). The very high denudation rate observed for sample P, and resultant high diffusion coefficient, can be a consequence of a recent anthropogenic disturbance, with the presence of small walls made from collected cobbles and possible shallow excavation in the vicinity of the sampling site. This data point is not further considered in the following analysis. The range of observed diffusion coefficients is consistent with values reported by available compilations for similar lithologies and climate (mean annual precipitation of 700 mm) (Hurst et al.2013b; Richardson et al.2019). There is no climatic gradient over the limited extent of our study area, so this cannot be invoked as a possible control for the  2-fold variations for the estimates of D in our dataset, which in any case do not present a clear spatial pattern or clustering. Similarly, bedrock geology is homogeneous between the different sites, with conglomerates of Miocene and Pliocene ages releasing clasts of an homogeneous 5–10 cm size group. We amalgamated up to 40 clasts at each sampling site and carefully selected well-defined near-horizontal ridges with negligible topographic shielding, such that we do not consider that this spread in D can arise from a systematic sampling bias. We note that the amplitude of variability for our estimates of D at individual hilltops is similar to the uncertainty range for D from studies based on CWDRs (e.g., Hurst et al.2013b). We consider the distribution of the observed D values to be mostly controlled by the natural variability of the hillslope processes and persistent internal transience at the  100 m wavelength. Such variability, observed here at individual hilltop locations, is usually averaged out in studies relying on CWDR rather than local estimations based on bedrock CRN inventories.

We use the estimated distribution for D to convert the CHT pattern along the transect into denudation rates using Eq. (3). We repeatedly sample the cumulative distribution for D and use CHT aggregated for hillslope patches. The resulting values for denudation rates are binned at 1 km intervals along the transect (Fig. 12). The interquartile ranges of the bins are largely overlapping, but a systematic increase can still be observed, consistent with the underlying variation in CHT. The estimated median denudation values are  30 and  50 mm kyr−1 in the southern and northern parts, respectively, which is comparable to denudation rates reported in surrounding landscapes (Siame et al.2004; Thomas et al.2017, 2018). We note that the denudation rate values measured at the various sites (blue circles in Fig. 12c) display an overall similar increasing trend but with systematically higher values. This deviation results from a sampling bias toward high-curvature ridges due to better field conditions at such sites. Indeed, measured CHT at these sites are ≥0.02m−1 (Fig. 12a), whereas continuously averaged CHT along the transect are on average below 0.02m−1 (Fig. 7b). Transport-limited conditions were prevailing at all surveyed sites such that we have no reason to question the validity of Eq. (3). For that reason the predominance of high-curvature/high-denudation sites in the dataset should not introduce a systematic deviation in the calibration of the diffusion coefficient D.

Under an assumption of steady-state conditions, we can propose that the observed 21 ± 17 mm kyr−1 differential denudation rate between the southern (distance < 5 km in Fig. 12c) and northern (distance >10 km) parts of the transect reflects a similar differential rock uplift rate across the transect. Considering a dip angle of 64 ± 11 (Fig. 11e) for the hypothetical fault responsible for the uplift pattern, it converts into a slip rate of 23 ± 20 mm kyr−1, which implies that slip rate is most likely < 50 mm kyr−1 on this blind thrust. This slow slip rate estimation is consistent with observations on slowly slipping faults of western Provence, where deformation usually cannot be resolved from geodetic data, and proposed long-term slip rates are < 100 mm kyr−1. For example, on the Trévaresse ridge fault, which produced the last major earthquake in metropolitan France (1909) with an estimated Mw=6 (Baroux et al.2003), trenches have yielded a late Quaternary slip rate in the 50 to 300 mm kyr−1 range (Chardon et al.2005), whereas accumulated deformation since the late Miocene indicates a slip rate of 30 ± 20 mm kyr−1 (Chardon and Bellier2003). For the Middle Durance Fault, which is located directly west of our study area, across the Durance river, a slip rate ranging from 10 to 70 mm kyr−1 has been proposed over the last Ma (Siame et al.2004). At last, the uplift rate we calculated is of the same order of magnitude as the one reported for the Lambrussier anticline by Hippolyte and Dumont (2000), directly north of our study area. In these slow-tectonics landscapes, erosion processes often outpace uplift rates, eroding passive deformation markers. While some aspects of the investigated area, such as the lithology, size and orientation of our basins, present favorable characteristics, our study illustrates the potential of tracking active erosion processes through hillslope morphology analysis to retrieve tectonic information in this type of environment.

6 Conclusions

Our analysis of hillslope properties along a transect into the Digne-Valensole basin, at the southern Alpine front, allows us to identify an important systematic variation across a short horizontal distance (< 10 km). In the absence of any climatic, lithologic or vegetation gradient, the observed increase in hilltop curvature, hillslope relief and normalized erosion rate points to a coincident increase in rock uplift. Hillslope lengths appear to be constant along the studied transect and thus unaffected by this major change in the tectonic boundary condition. Usual metrics associated with channel profile geometry, such as normalized steepness index, only capture the first-order change along the profile. Our results illustrate the utility of high-resolution analysis of hillslope morphology in low-uplift areas. Such techniques have the advantage to provide a dense spatial coverage, whereas approaches based on the fluvial metrics are inherently limited by the geometry and distribution of the river network.

We also demonstrate that, using simple deformation models, the relative uplift pattern resolved from high-resolution hillslope morphology analysis can be used to constrain the geometry and activity of underlying tectonic structures. We merged the morphological information obtained for hillslope morphology with an estimation of denudation rates from cosmogenic nuclide concentration measurements to evaluate the diffusion coefficient for hillslope sediment transport. In addition to the investigation of relative rock uplift patterns that are allowed by the approaches described above, the combined use of cosmogenic-nuclide-derived denudation rates at selected hilltop sites can be used to evaluate the amount of differential rock uplift across the transect. Such direct comparison of the spatial (high-resolution morphological analysis) and temporal (cosmogenic nuclides) aspects of landscape evolution has only been explicitly addressed by a limited number of studies and holds tremendous potential to decipher the response of various landscape elements to tectonic and climatic forcings.

Data availability

The DEM used in this study is part of the RGE ALTI® database from IGN (Institut National de l’information Géographique et Forestière;; IGN2017). This dataset is available under license from IGN.

Team list

Georges Aumaître (Aix-Marseille Univ., CNRS, IRD, INRAE, Coll France, CEREGE, Aix-en-Provence, France), Didier L. Bourlès (Aix-Marseille Univ., CNRS, IRD, INRAE, Coll France, CEREGE, Aix-en-Provence, France), and Karim Keddadouche (Aix-Marseille Univ., CNRS, IRD, INRAE, Coll France, CEREGE, Aix-en-Provence, France).

Author contributions

VG carried out the topographic analysis, processed the CRN samples and developed codes used in the interpretation. VG and JCH conducted field work. ASTER Team performed the AMS measurements. JCH, EC, NE, JF, OB and VO contributed to the understanding of the tectonics and geomorphology of the study area. VG prepared the article with contributions from all co-authors.

Competing interests

The authors declare that they have no conflict of interest.


The 10Be and 26Al measurements were performed at the ASTER AMS national facility (CEREGE, Aix en Provence), which is supported by the INSU/CNRS, the ANR through the “Projets thématiques d'excellence” program for the “Equipements d'excellence” ASTER-CEREGE action and IRD. We thank Franck Thomas for assistance in the field. Some of the computing for this project was performed on the OSU Pythéas HPC cluster. Insightful reviews by Martin D. Hurst, Marta Della Seta, Peter van der Beek and associate editor Veerle Vanacker allowed us to greatly improve our article.

Financial support

This research has been supported by the ECCOREV Federation (Transval). This work is also a contribution to the Labex OT-Med (ANR-11-LABX-0061) funded by the French Government “Investissements d'Avenir” program of the French National Research Agency (ANR) through the AMIDEX project (ANR-11-IDEX-0001-02).

Review statement

This paper was edited by Veerle Vanacker and reviewed by Martin D. Hurst, Peter van der Beek, and Marta Della Seta.


Arnold, M., Merchel, S., Bourlès, D. L., Braucher, R., Benedetti, L., Finkel, R. C., Aumaître, G., Gottdang, A., and Klein, M.: The French accelerator mass spectrometry facility ASTER: Improved performance and developments, Nucl. Instrum. Meth. B, 268, 1954–1959,, 2010. a, b

Balco, G.: Production rate calculations for cosmic-ray-muon-produced 10Be and 26Al benchmarked against geological calibration data, Quat. Geochronol., 39, 150–173,, 2017. a

Balco, G., Stone, J. O., Lifton, N. A., and Dunai, T. J.: A complete and easily accessible means of calculating surface exposure ages or erosion rates from 10Be and 26Al measurements, Quat. Geochronol., 3, 174–195,, 2008. a, b

Baroux, E., Pino, N. A., Valensise, G., Scotti, O., and Cushing, M. E.: Source parameters of the 11 June 1909, Lambesc (Provence, southeastern France) earthquake: A reappraisal based on macroseismic, seismological, and geodetic observations, J. Geophys. Res.-Sol. Ea., 108, 2454,, 2003. a

Beck, C., Deville, E., Blanc, E., Philippe, Y., and Tardy, M.: Horizontal shortening control of Middle Miocene marine siliciclastic accumulation (Upper Marine Molasse) in the southern termination of the Savoy Molasse Basin (northwestern Alps/southern Jura), Geol. Soc. Spec. Publ., 134, 263–278,, 1998. a

Benedetti, L., Tapponnier, P., King, G. C. P., Meyer, B., and Manighetti, I.: Growth folding and active thrusting in the Montello region, Veneto, northern Italy, J. Geophys. Res.-Sol. Ea., 105, 739–766,, 2000. a

Borchers, B., Marrero, S., Balco, G., Caffee, M., Goehring, B., Lifton, N., Nishiizumi, K., Phillips, F., Schaefer, J., and Stone, J.: Geological calibration of spallation production rates in the CRONUS-Earth project, Quat. Geochronol., 31, 188–198,, 2016. a

Boulton, S., Stokes, M., and Mather, A.: Transient fluvial incision as an indicator of active faulting and Plio-Quaternary uplift of the Moroccan High Atlas, Tectonophysics, 633, 16–33,, 2014. a

Brown, E. T., Edmond, J. M., Raisbeck, G. M., Yiou, F., Kurz, M. D., and Brook, E. J.: Examination of surface exposure ages of Antarctic moraines using in situ produced 10Be and 26Al, Geochim. Cosmochim. Ac., 55, 2269–2283,, 1991. a

Burbank, D. W., Leland, J., Fielding, E., Anderson, R. S., Brozovic, N., Reid, M. R., and Duncan, C.: Bedrock incision, rock uplift and threshold hillslopes in the northwestern Himalayas, Nature, 379, 505–510,, 1996. a

Champagnac, J., Molnar, P., Anderson, R., Sue, C., and Delacou, B.: Quaternary erosion-induced isostatic rebound in the western Alps, Geology, 35, 195,, 2007. a

Champagnac, J.-D., van der Beek, P., Diraison, G., and Dauphin, S.: Flexural isostatic response of the Alps to increased Quaternary erosion recorded by foreland basin remnants, SE France, Terra Nova, 20, 213–220,, 2008. a, b, c

Champagnac, J.-D., Molnar, P., Sue, C., and Herman, F.: Tectonics, climate, and mountain topography, J. Geophys. Res., 117, B02403,, 2012. a

Chardon, D. and Bellier, O.: Geological boundary conditions of the 1909 Lambesc (Provence, France) earthquake: structure and evolution of the Trévaresse ridge anticline, B. Soc. Geol. Fr., 174, 497–510,, 2003. a

Chardon, D., Hermitte, D., Nguyen, F., and Bellier, O.: First paleoseismological constraints on the strongest earthquake in France (Provence) in the twentieth century, Geology, 33, 901–904,, 2005. a

Clauzon, G.: The Messinian Rhone canyon as a definite proof of the “desiccated deep-basin model”, B. Soc. Geol. Fr., S7-XXIV, 597–610,, 1982. a, b

Clauzon, G.: Sequence boundaries and geodynamic evolution, Geomorphologie, 2, 3–21, 1996. a

Clauzon, G., Aguilar, J. P., and Michaux, J.: Investigation of the time-sedimentation relation thanks to examples taken from the French Mediterranean Neogene, B. Soc. Geol. Fr., V, 361–372,, 1989. a, b

Clauzon, G., Fleury, T. J., Bellier, O., Molliex, S., Mocochain, L., and Aguilar, J. P.: Morphostructural evolution of the Luberon since the Miocene (SE France), B. Soc. Geol. Fr., 182, 95–110,, 2011. a

Clubb, F. J., Mudd, S. M., Milodowski, D. T., Hurst, M. D., and Slater, L. J.: Objective extraction of channel heads from high-resolution topographic data, Water Resour. Res., 50, 4283–4304,, 2014. a

Clubb, F. J., Mudd, S. M., Attal, M., Milodowski, D. T., and Grieve, S. W.: The relationship between drainage density, erosion rate, and hilltop curvature: Implications for sediment transport processes, J. Geophys. Res.-Earth, 121, 1724–1745,, 2016. a, b

Clubb, F. J., Mudd, S. M., Milodowski, D. T., Valters, D. A., Slater, L. J., Hurst, M. D., and Limaye, A. B.: Geomorphometric delineation of floodplains and terraces from objectively defined topographic thresholds, Earth Surf. Dynam., 5, 369–385,, 2017. a, b

Clubb, F. J., Mudd, S. M., Hurst, M. D., and Grieve, S. W.: Differences in channel and hillslope geometry record a migrating uplift wave at the Mendocino triple junction, California, USA, Geology, 48, 184–188,, 2020. a, b, c

Collina-Girard, J. and Griboulard, R.: The deep-structure of the Valensole Plateau (French Alps). Contribution of network valley and topographic analysis, Géologie méditerranéenne, 17, 153–171, (last access: 4 January 2017), 1990. a, b, c

Cushing, E. M., Bellier, O., Nechtschein, S., Sébrier, M., Lomax, A., Volant, P., Dervin, P., Guignard, P., and Bove, L.: A multidisciplinary study of a slow-slipping fault for seismic hazard assessment: the example of the Middle Durance Fault (SE France), Geophys. J. Int., 172, 1163–1178,, 2008. a, b

Cyr, A. J., Granger, D. E., Olivetti, V., and Molin, P.: Quantifying rock uplift rates using channel steepness and cosmogenic nuclide-determined erosion rates: Examples from northern and southern Italy, Lithosphere, 2, 188–198,, 2010. a

Delacou, B., Sue, C., Champagnac, J.-D., and Burkhard, M.: Origin of the current stress field in the western/central Alps: role of gravitational re-equilibration constrained by numerical modelling, Geol. Soc. Spec. Publ., 243, 295–310,, 2005. a

Demoulin, A.: Morphometric dating of the fluvial landscape response to a tectonic perturbation, Geophys. Res. Lett., 39, 1–5,, 2012. a

DiBiase, R. A. and Whipple, K. X.: The influence of erosion thresholds and runoff variability on the relationships among topography, climate, and erosion rate, J. Geophys. Res., 116, F04036,, 2011. a

DiBiase, R. A., Heimsath, A. M., and Whipple, K. X.: Hillslope response to tectonic forcing in threshold landscapes, Earth Surf. Proc. Land., 37, 855–865,, 2012. a, b

Dietrich, W. E., Bellugi, D. G., Sklar, L. S., Stock, J. D., Heimsath, A. M., and Roering, J. J.: Geomorphic transport laws for predicting landscape form and dynamics, in: Prediction in Geomorphology, edited by: Wilcock, P. R. and Iverson, R. M., vol. 135 of Geophysical Monograph Series, pp. 103–132, American Geophysical Union, Washington, D.C., USA,, 2003. a, b, c

Dubar, M.: Chronology and significance of the late Neogene deposits of Riez-Valensole basin (Alpes-de-Haute-Provence, France), B. Soc. Geol. Fr., S7-XXVI, 971–978,, 1984a. a, b, c, d

Dubar, M.: Fluviatile terraces in South of the Alps, Bulletin de l'Association française pour l'étude du quaternaire, 21, 134–138,, 1984b. a, b

Dubar, M.: Correlation of the villafranchian mammalian sites of Puimoisson (valensole basin, Alpes de Haute-Provence, France) and others dated sites (palaeomagnetism, radionuclides) of Southern Europe (France, Italy and Spain), Quaternaire, 25, 85–90, 2014. a

Dubar, M. and Semah, F.: Paleomagnetic Data Bearing on the Age of High Terrace Deposits (Durance Sequence) in Alpine Valleys of Southeastern France, Quaternary Res., 25, 387–391,, 1986. a

Dubar, M., Guerin, C., and Heintz, E.: Les nouveaux gisements villafranchiens du ravin de Cornillet (Moustiers Sainte-Marie, Alpes de Haute Provence, France) et leur contexte géologique, Geobios, 11, 367–381,, 1978. a

Dubar, M., Aguilar, J. P., Chaline, J., Michaux, J., and Semah, F.: Chronological data (mammalian faunas, paleomagnetism) on the Plio-Pleistocene deposits from the top of the Valensole Basin: morphodynamic implications, Géologie de la France, 1, 57–68, 1998. a, b, c

Dubois, P. and Curnelle, R.: Résultats apportés par le forage les Mées no. 1 sur le plateau de Valensole (Alpes de Haute-Provence), Comptes-Rendus Sommaires des séance de la Société Géologique de France, 4, 181–184, 1978. a, b, c, d, e, f

Duvall, A., Kirby, E., and Burbank, D. W.: Tectonic and lithologic controls on bedrock channel profiles and processes in coastal California, J. Geophys. Res., 109, F03002,, 2004. a

Fox, M., Goren, L., May, D. A., and Willett, S. D.: Inversion of fluvial channels for paleorock uplift rates in Taiwan, J. Geophys. Res.-Earth, 119, 1853–1875,, 2014. a

Gabert, J.: Les terrasses de la Moyenne Durance, Bulletin de l'Association française pour l'étude du quaternaire, 16, 101–108,, 1979. a

Gallen, S. F., Wegmann, K. W., and Bohnenstiehl, D. W. R.: Miocene rejuvenation of topographic relief in the southern Appalachians, GSA Today, 23, 4–10,, 2013. a

Gelman, A. and Rubin, D. B.: Inference from Iterative Simulation Using Multiple Sequences, Stat. Sci., 7, 457–472,, 1992. a

Gidon, M. and Pairis, J.-L.: An example of the influence of coevel sedimentation on structures developed at the front of a thrust sheet moving in free-air: the structure of the region around Digne (SE France), Comptes rendus de l'Académie des Sciences, 307, 1283–1288, 1988. a

Glennie, C. L., Carter, W. E., Shrestha, R. L., and Dietrich, W. E.: Geodetic imaging with airborne LiDAR: the Earth's surface revealed, Rep. Prog. Phys., 76, 086801,, 2013. a

Godard, V., Ollivier, V., Bellier, O., Miramont, C., Shabanian, E., Fleury, J., Benedetti, L., and Guillou, V.: Weathering-limited hillslope evolution in carbonate landscapes, Earth Planet. Sc. Lett., 446, 10–20,, 2016. a, b

Godard, V., Dosseto, A., Fleury, J., Bellier, O., Siame, L., and ASTER Team: Transient landscape dynamics across the Southeastern Australian Escarpment, Earth Planet. Sc. Lett., 506, 397–406,, 2019. a, b

Grieve, S. W., Mudd, S. M., and Hurst, M. D.: How long is a hillslope?, Earth Surf. Proc. Land., 41, 1039–1054,, 2016a. a, b

Grieve, S. W. D., Mudd, S. M., Hurst, M. D., and Milodowski, D. T.: A nondimensional framework for exploring the relief structure of landscapes, Earth Surf. Dynam., 4, 309–325,, 2016b. a, b, c, d, e, f, g

Grieve, S. W. D., Mudd, S. M., Milodowski, D. T., Clubb, F. J., and Furbish, D. J.: How does grid-resolution modulate the topographic expression of geomorphic processes?, Earth Surf. Dynam., 4, 627–653,, 2016c. a, b

Harel, M.-a., Mudd, S., and Attal, M.: Global analysis of the stream power law parameters based on worldwide 10Be denudation rates, Geomorphology, 268, 184–196,, 2016. a, b, c

Hippolyte, J.-C. and Dumont, T.: Identification of Quaternary thrusts, folds and faults in a low seismicity area: examples in the Southern Alps (France), Terra Nova, 12, 156–162, 2000. a, b, c, d, e, f, g, h

Hippolyte, J.-C., Clauzon, G., and Suc, J.-P.: Messinian-Zanclean canyons in the Digne nappe (southwestern Alps): tectonic implications, B. Soc. Geol. Fr., 182, 111–132,, 2011. a, b, c, d

Howard, A. D.: A detachment-limited model of drainage basin evolution, Water Resour. Res., 30, 2261–2285,, 1994. a

Hsü, K. J., Montadert, L., Bernoulli, D., Cita, M. B., Erickson, A., Garrison, R. E., Kidd, R. B., Mèlierés, F., Müller, C., and Wright, R.: History of the Mediterranean salinity crisis, Nature, 267, 399–403,, 1977. a

Hurst, M. D., Mudd, S. M., Walcott, R., Attal, M., and Yoo, K.: Using hilltop curvature to derive the spatial distribution of erosion rates, J. Geophys. Res., 117, F02017,, 2012. a, b, c, d, e, f, g, h, i

Hurst, M. D., Mudd, S. M., Attal, M., and Hilley, G.: Hillslopes record the growth and decay of landscapes, Science, 341, 868–871,, 2013a. a, b, c, d

Hurst, M. D., Mudd, S. M., Yoo, K., Attal, M., and Walcott, R.: Influence of lithology on hillslope morphology and response to tectonic forcing in the northern Sierra Nevada of California, J. Geophys. Res.-Earth, 118, 832–851,, 2013b. a, b, c, d

Hurst, M. D., Grieve, S. W., Clubb, F. J., and Mudd, S. M.: Detection of channel-hillslope coupling along a tectonic gradient, Earth Planet. Sc. Lett., 522, 30–39,, 2019. a, b, c, d, e, f

IGN: Composante altimetrique du RGE®, available at:, last access: 1 September 2017. a

James, M. R. and Robson, S.: Straightforward reconstruction of 3D surfaces and topography with a camera: Accuracy and geoscience application, J. Geophys. Res., 117, F03017,, 2012. a

Jouanne, F., Hippolyte, J. C., Gamond, J. F., and Martinod, J.: Current deformation of the Digne Nappe (southwestern Alps) from a comparison between triangulation and GPS data, Geophys. J. Int., 144, 432–440,, 2001. a

Kirby, E. and Whipple, K. X.: Expression of active tectonics in erosional landscapes, J. Struct. Geol., 44, 54–75,, 2012. a, b, c

Kirby, E., Whipple, K. X., Tang, W., and Chen, Z.: Distribution of active rock uplift along the eastern margin of the Tibetan Plateau: Inferences from bedrock channel longitudinal profiles, J. Geophys. Res., 108, 2217,, 2003. a

Lague, D.: The stream power river incision model: Evidence, theory and beyond, Earth Surf. Proc. Land., 39, 38–61,, 2014. a

Lague, D. and Davy, P.: Constraints on the long-term colluvial erosion law by analyzing slope-area relationships at various tectonic uplift rates in the Siwaliks Hills (Nepal), J. Geophys. Res., 108, 2129,, 2003. a

Lashermes, B., Foufoula-Georgiou, E., and Dietrich, W. E.: Channel network extraction from high resolution topography using wavelets, Geophys. Res. Lett., 34, L23S04,, 2007. a

Le Roux-Mallouf, R., Godard, V., Cattin, R., Ferry, M., Gyeltshen, J., Ritz, J.-F., Drupka, D., Guillou, V., Arnold, M., Aumaître, G., Bourlès, D. L., and Keddadouche, K.: Evidence for a wide and gently dipping Main Himalayan Thrust in western Bhutan, Geophys. Res. Lett., 42, 3257–3265,, 2015. a

Lickorish, W. H. and Ford, M.: Sequential restoration of the external Alpine Digne thrust system, SE France, constrained by kinematic data and synorogenic sediments, Geol. Soc. Spec. Publ., 134, 189–211,, 1998. a

Lifton, N., Sato, T., and Dunai, T. J.: Scaling in situ cosmogenic nuclide production rates using analytical approximations to atmospheric cosmic-ray fluxes, Earth Planet. Sc. Lett., 386, 149–160,, 2014. a, b

Merchel, S. and Bremser, W.: First international 26Al interlaboratory comparison – Part I, Nucl. Instrum. Meth. B, 223-224, 393–400,, 2004. a

Merchel, S., Arnold, M., Aumaître, G., Benedetti, L., Bourlès, D., Braucher, R., Alfimov, V., Freeman, S., Steier, P., and Wallner, A.: Towards more precise 10Be and 36Cl data from measurements at the 10–14 level: Influence of sample preparation, Nucl. Instrum. Meth. B, 266, 4921–4926,, 2008. a

Mercier, H.: Le Néogène et le Pléistocène inférieur duranciens, Géologie Alpine, 55, 111–132, 1979. a, b

Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., and Teller, E.: Equation of State Calculations by Fast Computing Machines, J. Chem. Phys., 21, 1087–1092,, 1953. a

Milodowski, D. T., Mudd, S. M., and Mitchard, E. T. A.: Topographic roughness as a signature of the emergence of bedrock in eroding landscapes, Earth Surf. Dynam., 3, 483–499,, 2015. a

Mitasova, H., Hofierka, J., Zlocha, M., and Iverson, L. R.: Modelling topographic potential for erosion and deposition using GIS, Int. J. Geogr. Inf. Syst., 10, 629–641,, 1996. a

Morell, K. D., Sandiford, M., Rajendran, C. P., Rajendran, K., and Sanwal, J.: Geomorphology reveals seismogenic structures in the central Himalayan seismic gap, Lithosphere, 7, 247–256,, 2015. a

Mudd, S. M.: Detection of transience in eroding landscapes, Earth Surf. Proc. Land., 42, 24–41,, 2017. a, b

Myers, D. J., Nabelek, J. L., and Yeats, R. S.: Dislocation modeling of blind thrusts in the eastern Los Angeles basin, California, J. Geophys. Res.-Sol. Ea., 108, 2443,, 2003. a, b

Neely, A. B., DiBiase, R. A., Corbett, L. B., Bierman, P. R., and Caffee, M. W.: Bedrock fracture density controls on hillslope erodibility in steep, rocky landscapes with patchy soil cover, southern California, USA, Earth Planet. Sc. Lett., 522, 186–197,, 2019. a

Neteler, M., Bowman, M. H., Landa, M., and Metz, M.: GRASS GIS: A multi-purpose open source GIS, Environ. Modell. Softw., 31, 124–130,, 2012. a

Nicolas, M., Santoire, J., and Delpech, P.: Intraplate seismicity: new seismotectonic data in Western Europe, Tectonophysics, 179, 27–53,, 1990. a, b

Nishiizumi, K., Imamura, M., Caffee, M. W., Southon, J. R., Finkel, R. C., and McAninch, J.: Absolute calibration of 10Be AMS standards, Nucl. Instrum. Meth. B, 258, 403–413,, 2007. a

Nocquet, J.-M., Sue, C., Walpersdorf, A., Tran, T., Lenôtre, N., Vernant, P., Cushing, M., Jouanne, F., Masson, F., Baize, S., Chéry, J., and van der Beek, P. A.: Present-day uplift of the western Alps, Sci. Rep.-UK, 6, 28 404,, 2016. a, b

Okada, Y.: Surface deformation due to shear and tensile faults in a half-space, B. Seismol. Soc. Am., 75, 1135–1154, 1985. a, b

Ouimet, W. B., Whipple, K. X., and Granger, D. E.: Beyond threshold hillslopes: Channel adjustment to base-level fall in tectonically active mountain ranges, Geology, 37, 579–582,, 2009. a, b

Passalacqua, P., Belmont, P., Staley, D. M., Simley, J. D., Arrowsmith, J. R., Bode, C. a., Crosby, C., DeLong, S. B., Glenn, N. F., Kelly, S. a., Lague, D., Sangireddy, H., Schaffrath, K., Tarboton, D. G., Wasklewicz, T., and Wheaton, J. M.: Analyzing high resolution topography for advancing the understanding of mass and energy transfer through landscapes: A review, Earth-Sci. Rev., 148, 174–193,, 2015. a

Pelletier, J. D.: A robust, two-parameter method for the extraction of drainage networks from high-resolution digital elevation models (DEMs): Evaluation using synthetic and real-world DEMs, Water Resour. Res., 49, 75–89,, 2013. a

Pelletier, J. D., Nichols, M. H., and Nearing, M. A.: The influence of Holocene vegetation changes on topography and erosion rates: a case study at Walnut Gulch Experimental Watershed, Arizona, Earth Surf. Dynam., 4, 471–488,, 2016. a

Perron, J. T. and Royden, L.: An integral approach to bedrock river profile analysis, Earth Surf. Proc. Land., 38, 570–576,, 2012. a, b, c, d

Perron, J. T., Dietrich, W. E., and Kirchner, J. W.: Controls on the spacing of first-order valleys, J. Geophys. Res., 113, F04016,, 2008a. a

Perron, J. T., Kirchner, J. W., and Dietrich, W. E.: Spectral signatures of characteristic spatial scales and nonfractal structure in landscapes, J. Geophys. Res., 113, F04003,, 2008b. a

Perron, J. T., Kirchner, J. W., and Dietrich, W. E.: Formation of evenly spaced ridges and valleys, Nature, 460, 502–505,, 2009. a

R Core Team: R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria,, 2018. a

Richardson, P. W., Perron, J. T., and Schurr, N. D.: Influences of climate and life on hillslope sediment transport, Geology, 47, 1–4,, 2019. a

Ritz, J. F.: Recent tectonics and seismotectonics in the southern Alps: stress analysis, Quaternaire, 3, 111–124,, 1992. a

Roberts, G. G., Paul, J. D., White, N., and Winterbourne, J.: Temporal and spatial evolution of dynamic support from river profiles: A framework for Madagascar, Geochem. Geophy. Geosy., 13, 1–23,, 2012. a

Roering, J. J., Kirchner, J. W., and Dietrich, W. E.: Evidence for nonlinear, diffusive sediment transport on hillslopes and implications for landscape morphology, Water Resour. Res., 35, 853–870,, 1999. a, b, c

Roering, J. J., Kirchner, J. W., and Dietrich, W. E.: Hillslope evolution by nonlinear, slope-dependent transport: Steady state morphology and equilibrium adjustment timescales, J. Geophys. Res., 106, 16499–16513,, 2001. a, b

Roering, J. J., Perron, J. T., and Kirchner, J. W.: Functional relationships between denudation and hillslope form and relief, Earth Planet. Sc. Lett., 264, 245–258,, 2007. a, b, c, d, e, f

Roure, F., Brun, J.-P., Colletta, B., and Van Den Driessche, J.: Geometry and kinematics of extensional structures in the alpine foreland basin of southeastern France, J. Struct. Geol., 14, 503–519,, 1992. a, b

Ryan, W. B. and Cita, M. B.: The nature and distribution of Messinian erosional surfaces – Indicators of a several-kilometer-deep Mediterranean in the Miocene, Mar. Geol., 27, 193–230,, 1978. a

Scherler, D., Bookhagen, B., and Strecker, M. R.: Tectonic control on 10Be derived erosion rates in the Garhwal Himalaya, India, J. Geophys. Res.-Earth, 119, 1–23,, 2014. a

Schwartz, S., Gautheron, C., Audin, L., Dumont, T., Nomade, J., Barbarand, J., Pinna-Jamme, R., and van der Beek, P.: Foreland exhumation controlled by crustal thickening in the Western Alps, Geology, 45, 139–142,, 2017. a, b

Siame, L., Bellier, O., Braucher, R., Sébrier, M., Cushing, M., Bourlès, D., Hamelin, B., Baroux, E., de Voogd, B., Raisbeck, G., and Yiou, F.: Local erosion rates versus active tectonics: cosmic ray exposure modelling in Provence (south-east France), Earth Planet. Sc. Lett., 220, 345–364,, 2004. a, b

Stevenson, J. A., Sun, X., and Mitchell, N. C.: Despeckling SRTM and other topographic data with a denoising algorithm, Geomorphology, 114, 238–252,, 2010. a

Sue, C., Thouvenot, F., Fréchet, J., and Tricart, P.: Widespread extension in the core of the western Alps revealed by earthquake analysis, J. Geophys. Res.-Sol. Ea., 104, 25611–25622,, 1999. a

Sun, X., Rosin, P., Martin, R., and Langbein, F.: Fast and Effective Feature-Preserving Mesh Denoising, IEEE T. Vis. Comput. Gr., 13, 925–938,, 2007. a

Thomas, F., Godard, V., Bellier, O., Shabanian, E., Ollivier, V., Benedetti, L., Rizza, M., Espurt, N., Guillou, V., Hollender, F., and Molliex, S.: Morphological controls on the dynamics of carbonate landscapes under a mediterranean climate, Terra Nova, 29, 173–182,, 2017. a

Thomas, F., Godard, V., Bellier, O., Benedetti, L., Ollivier, V., Rizza, M., Guillou, V., Hollender, F., Aumaître, G., Bourlès, D. L., and Keddadouche, K.: Limited influence of climatic gradients on the denudation of a Mediterranean carbonate landscape, Geomorphology, 316, 44–58,, 2018. a

Tucker, G. E. and Whipple, K. X.: Topographic outcomes predicted by stream erosion models: Sensitivity analysis and intermodel comparison, J. Geophys. Res., 107, 2179,, 2002. a

von Blanckenburg, F.: The control mechanisms of erosion and weathering at basin scale from cosmogenic nuclides in river sediment, Earth Planet. Sc. Lett., 237, 462–479,, 2005. a, b

Walpersdorf, A., Pinget, L., Vernant, P., Sue, C., and Deprez, A.: Does Long-Term GPS in the Western Alps Finally Confirm Earthquake Mechanisms?, Tectonics, 37, 3721–3737,, 2018. a

Ward, S. N. and Valensise, G.: The Palos Verdes terraces, California: Bathtub rings from a buried reverse fault, J. Geophys. Res.-Sol. Ea., 99, 4485–4494,, 1994. a

Whipple, K. X. and Tucker, G. E.: Dynamics of the stream-power river incision model: Implications for height limits of mountain ranges, landscape response timescales, and research needs, J. Geophys. Res., 104, 17661–17674,, 1999. a

Whittaker, A. C.: How do landscapes record tectonics and climate?, Lithosphere, 4, 160–164,, 2012.  a

Wickert, A. D. and Schildgen, T. F.: Long-profile evolution of transport-limited gravel-bed rivers, Earth Surf. Dynam., 7, 17–3,, 2019. a

Wobus, C., Whipple, K. X., Kirby, E., Snyder, N. P., Johnson, J., Spyropolou, K., Crosby, B., and Sheehan, D.: Tectonics from topography: Procedures, promise, and pitfalls, Geol. S. Am. S., 398, 55–74,, 2006. a, b, c, d

Yanites, B. J., Tucker, G. E., and Anderson, R. S.: Numerical and analytical models of cosmogenic radionuclide dynamics in landslide-dominated drainage basins, J. Geophys. Res., 114, F01007,, 2009. a

Short summary
Slow-slipping faults are often difficult to identify in landscapes. Here we analyzed high-resolution topographic data from the Valensole area at the front of the southwestern French Alps. We measured various properties of hillslopes such as their relief and the shape of hilltops. We observed systematic spatial variations of hillslope morphology indicative of relative changes in erosion rates. These variations are potentially related to slow tectonic deformation across the studied area.