Hillslope denudation and morphologic response across a rock uplift gradient

Documenting the spatial variability of tectonic processes from topography is routinely undertaken through the analysis of river profiles, as 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 uplift rates. 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 5 and cosmogenic nuclides 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 Valensole Mio-Pliocene molassic basin, where a series of fold and thrust 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 hillslopes topographic properties such as hilltop curvature CHT or non-dimensional erosion rates E∗. We observed a systematic variation of these metrics coincident with the location of a major 10 underlying thrust system identified by seismic surveys. Using a simple deformation model, the inversion of the E∗ pattern allows 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 conglomerate at several hilltops locations for Be and Al concentration measurement. Calculated hilltops denudation rates range from 40 to 120 mm/ka. These denudation rates appear to be correlated with E∗ and CHT extracted from the morphological analysis, and are used to derive absolute estimates for the fault slip rate. This 15 high resolution hillslope analysis allows to resolve short wavelength variations in rock uplift, that would not be possible to unravel using commonly used channel profiles based methods. Our joint analysis of topography and geochronological data supports active thrusting at the Southwestern alpine front, and such approaches may bring crucial complementary constraints to morphotectonic analysis for slip rates on slow active faults. 1 20 https://doi.org/10.5194/esurf-2019-50 Preprint. Discussion started: 25 October 2019 c © Author(s) 2019. CC BY 4.0 License.


Introduction
The topography of the Earth evolves in response to surface processes driven by external forcings of tectonics and climatic origins (e.g.Champagnac et al., 2012;Whittaker, 2012).Tectonic uplift raises landscapes above their local base level, steepening rivers and hillslope topographic gradients.The water supply and temperature set by climatic conditions control the efficiency of weathering, erosion and transport processes, producing and moving sediment across the Earth surface.The present-day land surface morphology results from the accumulated actions of these two types of forcings through time, and a major endeavor of geomorphological research is to interpret measurable topographic properties in terms of space and time variations of either tectonic uplift or climatic conditions (e.g.Roberts et al., 2012;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 scale 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 been often motivated by the practical concern of identifying high strain zones in tectonicaly 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-profiles properties, as a direct relationship between fluvial gradient or relief and rock uplift had been identified in many incision models (Tucker and Whipple, 2002).Notably, the computation of morphometric parameters such as steepness indexes, a measure of channel gradient normalized for drainage area (Kirby and Whipple, 2012), 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 direct field measurements (e.g.Duvall, 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 (DEM : 10 to 100 m pixel size), which allow 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 or caveats when dealing with complexities of river systems such as lithological variations, transient evolution, or along-stream changes in fluvial dynamics, but are also inherently constrained by the planar structure of river network which might not always be optimal to sample the spatial 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 uplift rates.However, key elements of hillslope behavior such as the threshold stability angle for hillslope material and the non-linear relationship between sediment fluxes and topographic gradient (Roering et al., 1999(Roering et al., , 2001(Roering et al., , 2007)), implies that for fixed valley spacing, hillslope morphology can be insensitive to changes in erosion or rock uplift rates over large range 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 DEM, which are appropriate to describe river profiles developing over 1-100 km length scales, will only provide a very coarse description of hillslopes with lengths usually of the order of 100 m.
Nevertheless, this methodological limitation is currently being overcome by the growing use of LiDAR or photogrammetric techniques delivering Digital Terrain Models (DTM) with sub-metric resolutions (James and Robson, 2012;Glennie et al., 2013;Passalacqua et al., 2015).The widespread distribution of such data has spurred a great interest at new approaches to extract relevant information at the hillslope scale (e.g.DiBiase et al., 2012;Hurst et al., 2012;Grieve et al., 2016a, b;Milodowski et al., 2015).It is noteworthy that the high resolution of these DTMs allows to accurately compute the second derivative or curvature of the topographic surface, which covaries with erosion rate in linear diffusion framework, a relationship valid even for threshold hillslopes in the vicinity of the hilltop where topographic gradients are usually shallow (Hurst et al., 2012).
The possibility to access reliable proxies for erosion rate from hillslope scale metrics has lead to reconsider the potential of hillslope analysis for the assessment of denudation gradients.Notably, Roering et al. (2007) have proposed a conceptual framework, based on a formulation for non-linear 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 Radionuclides (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 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 localised change in the tectonic boundary condition is closely recorded by hillslope relief or slope angle and hilltope curvature extracted from LiDAR data, with the growth and decay phases of landscape evolution leaving a distinctive signature.
These promising results have offered critical insights into the dynamics of transient landscapes, through the point of view of hillslopes, with a spatial density of information several order of magnitude higher than what could be resolved with approaches based on fluvial network and profiles.Important methodological developments were necessary to extract the relevant information from high resolution DTMs (e.g.Grieve et al., 2016a).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., 2013aHurst et al., , 2019)).
There is therefore an urgent need to investigate further 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 confront 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 https://doi.org/10.5194/esurf-2019-50Preprint.Discussion started: 25 October 2019 c Author(s) 2019.CC BY 4.0 License.plateau in the Digne-Valensole Mio-Pliocene 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 nuclides methods we used, and we describe the corresponding datasets produced during this study.At last, we discuss the implications of these results in terms of the imprint of tectonic gradients into hillslope morphology, the constrain that can be put on tectonic structures at depth from high resolution topographic data and their relationship with denudation rates calculated from cosmogenic nuclides concentrations.

Settings
The studied area is located in the southern part of the Western European Alps, where mountain ranges result from a combination of Pyrenean (Late Cretaceous to Eocene) and Alpine (Neogene) tectonic phases (figure 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/year 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(Champagnac et al., , 2008) ) and in part to tectonic processes (Collina-Girard and Griboulard, 1990;Schwartz et al., 2017).
The flexural history of the Alps is particularly recorded in the Neogene flexural basins at the front of the Western Alps (e.g.Beck et al., 1998), such as the Digne-Valensole Basin in the southern Alps (Mercier, 1979).The Digne-Valensole Basin collected material eroded from the Alps and transported by the Durance, Bléone, Asse and Verdon rivers (figure 2).The basin is filled by marine deposits overlain by the continental Valensole conglomeratic 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 (Mercier, 1979;Clauzon et al., 1989;Dubar, 1984a;Dubar et al., 1998).The eastern edge of the basin is overthrusted by the Middle Miocene to Late Quaternary Digne Nappe (Lickorish and Ford, 1998;Roure et al., 1992;Gidon and Pairis, 1988;Hippolyte et al., 2011), and it is bordered to the West by the NNE-trending Durance seismic 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 allowed to demonstrate the Quaternary activity of several faults within the basin, and along its borders (Clauzon, 1982;Ritz, 1992;Cushing et al., 2008;Hippolyte and Dumont, 2000;Hippolyte et al., 2011).
The Digne-Valensole Basin is dissected by the Bléone and the Asse rivers which divide the area in 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 the Mediterranean sea-level dropped by about 1500 m (Ryan and Cita, 1978;Hsü et al., 1977;Clauzon, 1982;Clauzon et al., 2011).In the Valensole Basin, the Messinian paleo-canyons 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 (Dubar, 1984a;Clauzon, 1996).The Valensole-II Formation filled the Messinian canyons and covered the central and southern part of the Valensole basin (Dubar, 1984a) whose top surface presently forms the Valensole and Puimichel plateaus (figure 3).The age of the top surface of these plateaus ranges between 0.7 M.yr in the East, near the Digne Thrust (Dubar et al., 1998), and 1.7 M.yr in the West, along the Durance River (Dubar, 1984a;Dubar and Semah, 1986;Dubar, 2014).
The exceptional preservation of this recent surface allowed to demonstrate the Quaternary activity of the Lambruissier anticline (figure 1A) within the Digne-Valensole Basin (Hippolyte and Dumont, 2000).This SW-vergent fold generated a 80 m high and 5-km long morphological ridge above the Puimichel plateau surface.To the north-east of the Lambruissier anticline, an older Late Neogene fold is stratigraphically overlain by Late Pliocene horizontal deposits of the Bléone river messinian canyon (Hippolyte et al., 2011).These ages demonstrate the southward propagation of the 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 (Gabert, 1979;Dubar, 1984b), (2) the dip of the Puimichel Plateau surface which is more than three times higher (25 m/km) than the dip of the Valensole Plateau (8 m/km) (Hippolyte and Dumont, 2000) and (3) the striking dissymmetric pattern of the drainage network that may have recorded a progressive tilt of the Puimichel Plateau (Collina-Girard and Griboulard, 1990;Hippolyte and Dumont, 2000).Recent tectonic deformation of this area is also in agreement with modelling 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 (figure 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.The region is characterized by a Mediterranean climate with Mean Annual Temperature (MAT) of 13 o 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, which age is 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 globally 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 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.These characteristics leads us to consider these hillslopes as regolith-mantled and transport-limited systems, which is an important requirement of the topographic analysis described below (figure 3).

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 in the extraction of high resolution hillslope topographic metrics along a transect located at the northwestern edge of the Plateau.This analysis, and the associated parameters of interest, are based on the morphology of hillslope profiles predicted from Geomorphic Transport Laws used to describe hillslope sediment flux (Dietrich et al., 2003;Roering et al., 1999Roering et al., , 2007)).Sediment flux q s over a hillslope can be expressed as a non-linear function of hillslope gradient S, where ) and S c a critical hillslope gradient (Roering et al., 1999).One can note that for shallow slope areas, such as the vicinity of hilltops, the flux can be approximated as linear and the diffusion equation predicts a linear relationship between the erosion rate E and the second derivative of topography or hilltop curvature C HT , where β is the rock-to-regolith density ratio.Combining equation 1 with a statement for sediment conservation along the hillslope, the steady state elevation z as a function of horizontal coordinate x can be written as, With L H as the hillslope length, a reference erosion rates (Roering et al., 2007;Hurst et al., 2012;Grieve et al., 2015Grieve et al., , 2016a) ) can be defined as, and similarly a reference relief R R = S c L H which represents the maximum hillslope relief.
These two reference values allows to normalize hillslope relief R and erosion rate E into their non-dimensional equivalent as, and using equation 2, Finally, equation 3 can be expressed in non-dimensional form, 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 C HT , as well as hillslope relief R and length L H .These measurements will then allow to compute a non dimensional erosion rate E * according to equation 6.We use the approach presented by Hurst et al. ( 2012) and Grieve et al. (2016a), which we implemented into the GRASS GIS environment (Neteler et al., 2012) and R scripting language (R Core Team, 2018), as described in Godard et al. (2019) (figure 4).
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 planar curvature threshold for the definition of channel heads (e.g.Pelletier, 2013;Clubb et al., 2014).The narrow floodplains 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 range of stream orders and curvature was computed at every hilltop pixel by fitting a quadratic surface over a 30 m wide window (Hurst et al., 2012;Godard et al., 2016).Flow was routed downslope from hilltop pixels to the channels using the algorithm of Mitasova et al. (1996), and the resulting flowlines were used to compute hillsope relief and length (Hurst et al., 2012;Grieve et al., 2015Grieve et al., , 2016a)).Flowlines and associated data were grouped into patches of a least 50 contiguous hilltop pixels (Grieve et al., 2016a).
We also extracted standard metrics for the river profiles of the studied catchments (figure 5).We used χ-transformation of the river long-profiles to determine the optimal concavity for each basin (Perron and Royden, 2012).We then again χ-transformed the river profiles using a reference concavity at the scale of the study area to compute normalized steepness indexes.

Cosmogenic nuclides
The topographic analysis methods described in the previous section allow to identify relative spatial patterns for the intensity of surface processes such as surface denudation, and under a steady-state assumption, these pattern can be used to infer rock uplift distribution.We used in situ-produced cosmogenic nuclides concentration measurements to constrain the absolute denudation rates values.Active fluvial sediments are present along the stream network of the studied catchments (figure 2), and could be sampled to derive basin-averaged denudation rates (von Blanckenburg, 2005).However, most of the surveyed channels displayed a complex dynamics with many occurrences of aggradation, splitting of the main drain and local colluvial inputs, such that the required hypothesis of an homogeneous contribution of the whole upstream area was most likely invalid.For that reason we rather sampled material directly at the hilltops, which allows 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.leachings in concentrated HF, each removing 10% of the remaining sample mass (Brown et al., 1991).After addition of an in-house 9 Be carrier (∼100 µl at 3.025 x 10 −3 g/g (Merchel et al., 2008)), the samples were digested in concentrated HF.Al concentration in cleaned samples were in the 100-200 ppm range so no Al carrier solution was added to the samples.Be and Al were isolated for measurements using ion-exchange chromatography.
10 Be/ 9 Be and 26 Al/ 27 Al measurements were performed at the French AMS National Facility, located at CEREGE in Aixen-Provence. 10Be/ 9 Be 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).The measured 26 Al/ 27 Al 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 crosscalibrated against primary standards from a round-robin exercise (Merchel and Bremser, 2004).
Uncertainties on 10 Be and 26 Al concentrations (reported as 1σ) 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 to 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.Two process blanks were treated and measured with our samples, yielding 10 Be/ 9 Be ratios of 1.47±0.27and 0.90±0.29 ×10 −15 .It corresponds to an upper 1σ bound of 55 and 38 ×10 3 10 Be atoms for the background level in these two blanks, which is at least 20 times lower than the number of 10 Be atoms in the dissolved sample masses (30 times lower in average over our dataset).
Denudation rates were then calculated with the online calculator described in Balco et al. (2008) and the nuclide specific LSD 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).The effective attenuation length for spallation in rock is 160 g/cm 2 , and density is 2.5.No shielding correction was considered for the hilltop sites we sampled.

Results
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.

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 (figure 6C).Calculation of average hillslope gradient shows that most hillslope have average gradients below 0.6 (figure 6B).Following (Hurst et al., 2019), we determined the value of S c for which 99% of hillslopes have a relief inferior to the maximum value (L H S c ), and obtained S c = 0.52.In the following we use S c = 0.6 for the critical gradient value.Most topographic metrics extracted at the hillslope scale display important variations along the studied profile, both in terms of basin-averaged or binned values.Non-dimensional erosion rate E * increases 2-fold from North to South (figure 7).This variation is not evenly distributed along the profile but occurs over less than 4 km of horizontal For all studied basins the river profiles display a regularly concave-up shape, and in most situation χ-transformed main trunk profiles collapses along a linear trend, as well as the main tributaries, suggesting the absence of major transient perturbation propagating through the river network (figure 5).Usual topographic indexes, such as stream concavity and normalized steepness index, were extracted from the fluvial network.Concavity ranges from 0.2 to 0.4 and a reference value of 0.3 was used in the following analysis.For both parameter no clear progressive evolution can be observed from North to South, and 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.However, it can be noted that the 4 northernmost basins display slightly higher k sn values, in average, than the rest of the section.

Cosmogenic nuclides data
Measured concentrations in our hilltop samples range from 43 to 117×10 3 at/g and 281 to 867×10 3 at/g for 10 Be and 26 Al, respectively (table 1).The corresponding 26 Al/ 10 Be ratios vary from 6.5±0.8 to 8.2±1.2.While some of these ratios are slightly higher than the theoretical production value, 1σ uncertainties ellipses are always overlapping the steady state denudation curve (figure 9).These concentrations correspond to denudation rates ranging from 42 to 102 mm/ka and from 38 to 88 mm/ka for 10 Be and 26 Al, respectively (table 2).The denudation rates calculated from both nuclides concentrations are consistent but display a small deviation from the 1:1 line, with 26 Al denudation rates slightly lower than their 10 Be equivalent.The observed 26 Al/ 10 Be 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 10 Be data for our analysis of the relationship between high resolution topography and denudation rates.

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.

Controls on hillslopes morphology
We observe a pronounced and systematic variations of hillslope morphology along the studied transect, with hillslope curvature C HT undergoing a 2-fold increase from 7 to 10 km (figure 7B).We evaluate below the possible controls on this evolution.We first note that our study area extends over ∼10 km in length with no important changes in the range of elevation of the catchments from South to North, such that climatic conditions can be considered as constant in terms of Mean Annual Precipitation and Temperature (MAP and MAT), and can not account for the observed variations in hillslope morphology.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.
Hillslopes have been shown to record transient waves of erosion propagating through landscapes (Hurst et al., 2012(Hurst et al., , 2013a;;Mudd, 2017).The series of studied basins are directly connected to the Durance river baselevel, and the eventual propagation of incision pulses and knickpoints alongstream 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 knickpoint have been identified along the Durance river in the vicinity of our working area.Second, we note that the pattern of evolution for E * , C HT 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 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 variable rock uplift pattern associated with recent or ongoing deformation.Several studies have already pointed at long-wavelength geomorphic evidences for recent tectonic activity in the northern part of the Puimichel Plateau (Collina-Girard and Griboulard, 1990; Hippolyte and Dumont, 2000), which is confirmed by recent seismic activity (Nicolas et al., 1990), with several magnitude 4 thrust-slip events with EW strike directions (figure 1B).Notably, a major inflection of the abandonment surface occurs at 7-8 km of horizontal distance along the transect (figure 2B), suggesting post-Pliocene uplift of the northern part.The location of this surface inflection coincides with the transition area for E * , C HT and R (figure 7).Such spatial coherence between a long-term finite deformation pattern passively recorded at 1 Ma time-scale 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 located below the Northern part of the Puimichel Plateau (figure 2C and 2D).The uplift itself 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 Curnelle, 1978).A Quaternary reactivation of such faults had already been invoked to explain the surface deformation pattern of the Plateau farther North (Hippolyte and Dumont, 2000), and we propose a dominant tectonic control for the distribution of proxies for surface denudation along our transect (figure 7).
The R * vs E * relationship (with S c = 0.6) shows that the studied catchments plot slightly below the steady state line (figure 8), indicating a possible decaying dynamics of the landscape, eventually due to a change in the climatic or tectonic boundary condition Mudd (2017).All catchments appears to be similarly affected, with no obvious gradient between the northern and southern ones.For that that reason this decay is not likely to result from a decrease in the amount of differential rock uplift across the transect, but could be associated with a change in the intensity of top-down forcings, of climatic origin.However, we note that using S c = 0.5, which appears to be reasonable for the vast majority of the studied hillslopes (figure 6), reduces considerably the deviation from the steady state line.

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 Whipple, 2012), do not display as clear a pattern as hillslope metrics (figure 7D).Whereas it can be argued that normalized steepness indexes values (ks n ) are in average higher in the northern part of the transect with respect to the southern part, a lot of scatter is present in that data and, most importantly, we do not observe any clear gradual increase comparable to what is displayed by E * or C HT .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 (figure 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 Royden, 2012), and argue against the impact of local perturbation along the river profile as an origin for the observed scatter.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 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.A key difference with our study is 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 difference might be a reason for the better sampling through fluvial metrics by Hurst et al. (2019), and the clearer relationship they observe with hillslopes 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 (figure 7).Indeed, hillslope length L H remains nearly constant across the transect between 140 and 160 m, whereas for example hilltop curvature C HT , which can be considered as a proxy for erosion and rock uplift under a steady-state assumption, undergoes a nearly two-fold increase.There is no significant inverse correlation between L H and C HT as observed by Hurst et al. (2013b).The characteristic horizontal lengthscale of landscapes, which can be evaluated through different type of measurements such as drainage density, spacing of 1 st order valleys or hillslope length, has been shown to be highly sensitive to external tectonics 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 details 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/ka range corresponding to our CRN data (figure 9), and conflicts with the observation of a nearly constant L H 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. (2008bPerron et al. ( , 2009) ) as the ratio between characteristics fluvial and hillslope timescales : where K is the erodibility parameter for fluvial incision, D the hillslope diffusion coefficient, m and n the area and slope exponents of the fluvial incision law respectively, and l and ζ horizontal and vertical lengthscales 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 P e = 1 allows to retrieve l c as a characteristic lengthscale for hillslope/channel transition, which again, in the n = 1 case, does not depends on erosion or uplift rates.Therefore, the stability of L H across a two-fold erosion rate gradient, observed in our dataset, would hint toward a value of n close to unity.
While the ratio m/n 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 evidences pointing to n > 1 DiBiase and Whipple, 2011;Lague, 2014;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 Davy, 2003), 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 profiles analysis.

Constraints on tectonic processes from hillslope morphology
Observations of absolute or relative spatial 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).As discussed in the previous section, the hypothesis that the evolution of hillslope morphology along the studied transect (figure 7) results from a spatial variations 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 a simple elastic dislocation model (Okada, 1985) to predict surface deformation distributions and compare them with our observed E * evolution along the profile presented on figure 7.This type of modelling approach is usually associated with the study of upper crustal deformation at the scale of the seismic cycle, but has also been successfully applied to interpret surface deformation patterns associated with blind thrusts over longer timescales (Benedetti et al., 2000;Myers et al., 2003).We consider a single planar dislocation embedded in an elastic medium and vary the horizontal and vertical https://doi.org/10.5194/esurf-2019-50Preprint.Discussion started: 25 October 2019 c Author(s) 2019.CC BY 4.0 License.
positions of its upper limit, its dip angle and length.We compare the predicted surface deformation profile to the observed E * pattern.We do not fit the absolute value of E * which has no signification with respect to the elastic dislocation model, but rather the wavelength of the predicted deformation and simply scale its amplitude to that of the observed E * profile.
We perform a MCMC driven sampling of the 4-parameters space and the corresponding marginal distributions are presented on figure 10.We observe that the horizontal position of the top of the dislocation is the best constrained parameter with a most likely value around 8 km (same horizontal reference frame as the profile on figure 7).The most likely depth for the upper limit of the dislocation is between 2 and 4 km and the suggested dip of the structure is >50 o .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 (figure 10B).The high E * northern part of the studied transect roughly corresponds to the Mées structure identified from seismic surveys and drilling (Dubois and Curnelle, 1978), which is a large anticline inherited from the late-Cretaceous-Eocene compression of the Pyrenean phase (figure 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 Curnelle, 1978), and that the Cenozoic formations are submitted to 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 whose Quaternary reactivation 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 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 Dumont, 2000).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(Hurst et al., , 2019)), our results are the first to illustrate the use of such data to infer the geometry of tectonic structures at depth.

Comparison of cosmogenic nuclides data with high resolution topography metrics
Our understanding of landscape dynamics relies on the formulation of Geomorphic Transport Laws (GTL), such as equation 1, 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 model requires to obtain actual measurements for these spatial and temporal descriptions of the landscape.Spatial properties of landscapes are usually derived from the analysis of Digital Elevation Models at various resolutions, upon which topographic gradient For example, the confrontation of Catchment-Wide Denudation Rates (CWDR) calculated from measured 10 Be concentrations in river sediments with topographic gradient extracted from Digital Elevation Models has provided critical tests of GTL 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 properties 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 suffer 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., 2012Hurst et al., , 2013b;;Godard et al., 2016Godard et al., , 2019;;Neely et al., 2019).
Our dataset allows to carry out such confrontation 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 formulation such as equation 2 relating hilltop curvature C HT 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 C HT and erosion rates are compatible with D in the 0.003 to 0.006 m 2 /a range (figure 11).The very high denudation rate observed at site P, and resultant high diffusion coefficient, can be a consequence of a recent anthropogenic disturbance observed at the hilltop and 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 this 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 can not be invoked as possible control for the ∼2 fold variations for the estimates of D in our dataset, which in any case does not present a clear spatial pattern or clustering.Similarly, bedrock geology is homogenenous between the different sites, with conglomerates of Miocene and Pliocene ages releasing clasts of an homogeneous 5-10 cm size group.We note that, over our dataset, the amplitude of variability for our estimates of D at individual hilltops is similar to the uncertainty range for D from studies based on CWDR (e.g.Hurst et al., 2013b), rather than local estimation of denudation rates as in our case.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 basin-wide denudation rates rather than local estimations based on bedrock CRN inventories.
We use the estimated distribution for D to convert the C HT pattern along the transect into denudation rates using equation 2.
We repeatedly sample the cumulative distribution for D and use C HT aggregated for hillslope patches.The resulting values for denudation rates are binned at 1 km interval along the transect (figure 11).The interquartile ranges of the bins are largely overlapping, but a systematic increase can still be observed, consistently with the underlying variation in C HT .The estimated median denudation values are ∼30 and ∼50 mm/ka in the southern and northern part, respectively, which is comparable to denudation rates reported on surrounding landscapes (Siame et al., 2004;Thomas et al., 2017Thomas et al., , 2018)).Considering a dip angle of 60 o (figure 10E) for the hypothetical fault responsible for this differential uplift rate, and assuming topographic steady state, this would imply a ∼20 mm/ka slip rate on this blind thrust.This slip rate is similar to observations on slow motion faults of Western Provence, where deformation usually can not be resolved from geodetic data, and proposed long-term slip rates are usually < 100 mm/ka.For example, on the Trevaresse ridge fault, which produced the last major earthquake in metropolitan France (1909) with an estimated M w = 6 (Baroux et al., 2003), trenches have yielded a Late Quaternary slip rate in the 50 to 300 mm/ka range (Chardon et al., 2005), whereas accumulated deformation since the late Miocene indicates a slip rate of 30±20 mm/ka (Chardon and Bellier, 2003).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/ka has been proposed over the last Ma (Siame et al., 2004).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 environments.

Conclusions
Our analysis of hillslope properties along a transect into the Digne-Valensole basin, at the Southern Alpine front, allows 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 globally 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 concavity or steepness index, do not show similar variations along the studied profile.This absence of response suggests that such short-wavelength variations in rock uplift in slowly deforming settings are poorly resolved with these approaches.
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 hillslopes morphology, with estimation of denudation rates from cosmogenic nuclides concentration measurements, to evaluate the diffusion coefficient for hillslope sediment transport.In addition to the investigation of relative rock uplift patterns allowed by the approaches described above, the combined use of cosmognenic nuclides derived denudation rates at selected hilltop sites can be used to evaluate the amount of differential rock uplift across the transect.Such direct confrontation of the spatial (high resolution morphological analysis) and temporal (cosmogenic nuclides) Table 2. Calculated denudation rates using the online calculator from (Balco et al., 2008) and the LSD scaling scheme (Lifton et al., 2014).
Reported uncertainties are ±1σ and combine the internal and external uncertainties.τ10 is 10 Be denudation rates integration timescale.

Samples for 10
Be and 26 Al concentrations measurements were collected at 10 hilltop sites, by amalgamating 30 to 40 individual sandstone clasts derived from the bedrock conglomerate (figure2A).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 H 2 SiF 6 and submitted to vigorous mechanical shaking until pure quartz was obtained.Decontamination from atmospheric 10 Be was achieved by a series of three successive https://doi.org/10.5194/esurf-2019-50Preprint.Discussion started: 25 October 2019 c Author(s) 2019.CC BY 4.0 License.
https://doi.org/10.5194/esurf-2019-50Preprint.Discussion started: 25 October 2019 c Author(s) 2019.CC BY 4.0 License.distance.The hilltop curvature C HT pattern closely mimics that of E * , increasing from 0.01 to 0.02 m −1 .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 show an increase from ∼40 to ∼60 m.On the contrary no clear systematic changes in hillslope length can be observed along the profile, with values ranging from 140 to 160 m.
https://doi.org/10.5194/esurf-2019-50Preprint.Discussion started: 25 October 2019 c Author(s) 2019.CC BY 4.0 License.and curvature can be computed.On the other hand, geochronological techniques, such as cosmogenic nuclides concentration measurements in bedrock or sediment samples, can allow to constrain the rate of lowering of the topographic surface though time, provide the framework to evaluate the temporal part of the landscape evolution problem.