Hillslope denudation and morphologic response to a rock uplift gradient

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 C_HT 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 ¹⁰Be and ²⁶Al concentration measurements. Calculated hilltop denudation rates range from 40 to 120 mm kyr⁻¹. These denudation rates appear to be correlated with E^∗ and C_HT 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; Whittaker, 2012). 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; Demoulin, 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 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 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 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 Robson, 2012; 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⁻¹ 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 Griboulard, 1990; Schwartz et al., 2017).

Figure 1 (a) Geological setting of the studied region in SE France (1∕10⁶ 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 (Mercier, 1979). 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 (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 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 (Clauzon, 1982; Ritz, 1992; Cushing et al., 2008; Hippolyte and Dumont, 2000; 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 Cita, 1978; Hsü et al., 1977; Clauzon, 1982; 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 (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), 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 (Dubar, 1984a, 2014; Dubar and Semah, 1986).

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 (Dubar, 1984b), 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 (Gabert, 1979; Dubar, 1984b), (2) the dip of the Puimichel Plateau surface which is more than 3 times higher (25 m km⁻¹) than the dip of the Valensole Plateau (8 m km⁻¹) (Hippolyte and Dumont, 2000) and (3) the striking asymmetric 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 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 Curnelle, 1978). 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

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial q_{\mathrm{s}}}{\partial x}}=\mathit{\beta }E,\end{array}$$

where x is the horizontal coordinate ([L]), q_s is the sediment flux ([L²T⁻¹]), β is the rock-to-regolith density ratio (dimensionless), and E is the erosion rate ([LT⁻¹]) 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 q_s over a hillslope as a nonlinear function of local hillslope gradient $S=\partial z/\partial x$ (Roering et al., 1999, 2007),

$$\begin{array}{}\text{(2)}& q_{\mathrm{s}}={\displaystyle \frac{DS}{\mathrm{1}-{\left(\frac{S}{S_{\mathrm{c}}}\right)}^{\mathrm{2}}}},\end{array}$$

where D is a diffusion coefficient ([L²T⁻¹]) and S_c 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 q_s=DS. Combining this expression of q_s with Eq. (1) yields a linear relationship between the erosion rate E and the second spatial derivative of topography or hilltop curvature C_HT,

$$\begin{array}{}\text{(3)}& E={\displaystyle \frac{D{C}_{\text{HT}}}{\mathit{\beta }}}.\end{array}$$

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 (C_HT) 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),

$$\begin{array}{ll}{\displaystyle }z\left(x\right)={\displaystyle \frac{D{S}_{\mathrm{c}}^{\mathrm{2}}}{\mathrm{2}\mathit{\beta }E}}& {\displaystyle }\left(\ln \left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sqrt{\mathrm{1}+{\left({\displaystyle \frac{\mathrm{2}\mathit{\beta }Ex}{D{S}_{\mathrm{c}}}}\right)}^{\mathrm{2}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\right)\right.\\ \text{(4)}& {\displaystyle }& {\displaystyle }\left.-\sqrt{\mathrm{1}+{\left({\displaystyle \frac{\mathrm{2}\mathit{\beta }Ex}{D{S}_{\mathrm{c}}}}\right)}^{\mathrm{2}}}+\mathrm{1}\right).\end{array}$$

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

$$\begin{array}{}\text{(5)}& E_R={\displaystyle \frac{D{S}_{\mathrm{c}}}{\mathrm{2}{L}_{\mathrm{H}}\mathit{\beta }}},\end{array}$$

and similarly a reference relief R_R=S_cL_H 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

$$\begin{array}{}\text{(6)}& R^{\ast }={\displaystyle \frac{R}{S_{\mathrm{c}}{L}_{\mathrm{H}}}},\end{array}$$

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

$$\begin{array}{}\text{(7)}& E^{\ast }={\displaystyle \frac{\mathrm{2}{C}_{\text{HT}}{L}_{\mathrm{H}}}{S_{\mathrm{c}}}}.\end{array}$$

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

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{R}^{\ast }& ={\displaystyle \frac{\mathrm{1}}{E^{\ast }}}(\sqrt{\mathrm{1}+(E^{\ast }{)}^{\mathrm{2}}}-\ln \left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}+\sqrt{\mathrm{1}+(E^{\ast }{)}^{\mathrm{2}}}\right)\right)\\ & -\mathrm{1}).\end{array}\end{array}$$

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 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 Team, 2018), 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⁻¹ contour curvature threshold for the definition of channel heads (e.g., Pelletier, 2013; 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⁻¹]) can be parameterized as a function of along-channel topographic gradient S and drainage area A ([L²]) as

$$\begin{array}{}\text{(9)}& I=K{A}^m{S}^n,\end{array}$$

where K is an erodibility coefficient ([L^1−2mT⁻¹]), and m and n are empirical exponents (Howard, 1994; Whipple and Tucker, 1999). Under steady-state conditions, river incision I equals rock uplift U and Eq. (9) can be reorganized as

$$\begin{array}{}\text{(10)}& S=\left|{\displaystyle \frac{\mathrm{d}z}{\mathrm{d}x}}\right|={\left({\displaystyle \frac{U}{K}}\right)}^{\frac{\mathrm{1}}{n}}{A}^{-\frac{m}{n}}.\end{array}$$

Parameters $k_{\mathrm{s}}={\left(\frac{U}{K}\right)}^{\frac{\mathrm{1}}{n}}$ and $\mathit{\theta }=\frac{m}{n}$ are referred to as the steepness index and channel concavity (Wobus et al., 2006; Kirby and Whipple, 2012), 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 m∕n value is chosen and used in the regression in order to obtain a normalized steepness index k_sn.

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 Royden, 2012), using $m/n=\mathrm{0.23}$. Inset shows the evolution of the R² of the linear regression of elevation against χ for a range of m∕n 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 x_b, to x as

$$\begin{array}{}\text{(11)}& z\left(x\right)-z\left(x_{\mathrm{b}}\right)=\underset{x_{\mathrm{b}}}{\overset{x}{\int }}{\left({\displaystyle \frac{U}{K{A}^m}}\right)}^{\frac{\mathrm{1}}{n}}\mathrm{d}x.\end{array}$$

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

$$\begin{array}{}\text{(12)}& z\left(x\right)=z\left(x_{\mathrm{b}}\right)+{\left({\displaystyle \frac{U}{K{A}_{\mathrm{0}}^m}}\right)}^{\frac{\mathrm{1}}{n}}\mathit{\chi },\end{array}$$

where A₀ is a reference drainage area, with

$$\begin{array}{}\text{(13)}& \mathit{\chi }=\underset{x_{\mathrm{b}}}{\overset{x}{\int }}{\left({\displaystyle \frac{A_{\mathrm{0}}}{A}}\right)}^{\frac{m}{n}}\mathrm{d}x.\end{array}$$

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 m∕n value (concavity) leading to the best linearization of the χ-transformed river profile (Perron and Royden, 2012). We fixed the reference value according to the mean of observed m∕n value across all basins. We then again χ-transformed the river profiles using this reference m∕n 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 Blanckenburg, 2005). 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 ¹⁰Be and ²⁶Al 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 H₂SiF₆ and submitted to vigorous mechanical shaking until pure quartz was obtained. Decontamination from atmospheric ¹⁰Be 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 ⁹Be carrier, the samples were digested in concentrated HF. Be and Al were isolated for measurements using ion-exchange chromatography. ¹⁰Be∕⁹Be and ²⁶Al∕²⁷Al 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 1 Cosmogenic nuclide ¹⁰Be results.

^a Dissolved pure quartz mass. ^b In-house carrier mass, ∼ 150 µL at 3.025 × 10⁻³  g g⁻¹ (Merchel et al., 2008). ^{c 10}Be∕⁹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⁻¹¹ (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 ¹⁰Be∕⁹Be ratios of 1.47 ± 0.27 and 0.90 ± 0.29 × 10⁻¹⁵. 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 ¹⁰Be 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⁻² for the effective attenuation length for spallation in rock and a density of 2.5 g cm⁻³. No shielding correction was considered for the hilltop sites we sampled, which were selected on nearly horizontal ridgelines.

Table 2 Cosmogenic nuclide ²⁶Al results and comparison with ¹⁰Be 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 ²⁶Al∕²⁷Al ratios were normalized to the in-house standard SM-Al-11 whose ratio of 7.401 ± 0.064 × 10⁻¹² (Arnold et al., 2010) has been cross-calibrated against primary standards from a round-robin exercise (Merchel and Bremser, 2004). ^c Uncertainties are reported at the 1σ level. ^d An analytical blank yielded a ratio of 7.27 ± 2.17 × 10⁻¹⁵. ^e Procedures for that calculation of denudation rates are identical to the ones described in Table 1. ^f Integration timescales for ¹⁰Be and ²⁶Al denudation rates (von Blanckenburg, 2005).

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 (Okada, 1985) 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 Valensise, 1994; 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 10⁶ in length with a 10⁴ length burn-in phase. The multivariate potential scale reduction factor is 1.004, suggesting that convergence was achieved (Gelman and Rubin, 1992).

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⁻¹ (Fig. 6a). Calculation of average hillslope gradient shows that most hillslopes have average gradients below 0.6 m m⁻¹ (Fig. 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_HS_c) and obtained S_c=0.52 m m⁻¹. In the following we use S_c=0.6 m m⁻¹ for the critical gradient value (Fig. 6c). This S_c 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 C_HT pattern closely mimics that of E^∗, increasing from 0.01 to 0.02 m⁻¹ (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 C_HT 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 7 Projection 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^{\ast }=\mathrm{2}{L}_{\mathrm{H}}{C}_{\text{HT}}/S_{\mathrm{c}}$ (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 m∕n 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 m∕n ratio and normalized steepness index (k_sn), were extracted from the fluvial network. The m∕n 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 m∕n 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 k_sn values than the rest of the section, and that basin averaged E^∗ and k_sn are significantly positively correlated (Fig. 8).

Figure 8 Evolution of basin averaged E^∗ as a function of normalized steepness index (k_sn). 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 km²). Solid and dashed lines correspond to a linear fit and its 95 % confidence interval (R²=0.62 and $p<{\mathrm{10}}^{-\mathrm{3}}$).

4.2 Cosmogenic nuclide data

Measured concentrations in our hilltop samples range from 43 to 117×10³ atoms g⁻¹ and 281 to 867×10³ atoms g⁻¹ for ¹⁰Be and ²⁶Al, respectively (Tables 1 and 2). The corresponding ²⁶Al∕¹⁰Be 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⁻¹ and from 38 to 114 mm kyr⁻¹ for ¹⁰Be and ²⁶Al, 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 ²⁶Al denudation rates slightly lower than their ¹⁰Be equivalents. The observed ²⁶Al∕¹⁰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 ¹⁰Be 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 ²⁶Al∕¹⁰Be ratio for steady-state denudation and constant exposure histories, respectively. (b) Comparison of denudation rates calculated from measured ¹⁰Be and ²⁶Al 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 C_HT 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; Mudd, 2017). 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^∗, 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 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 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 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^∗, C_HT 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 Curnelle, 1978). Quaternary reactivation of such faults has been invoked to explain the surface deformation pattern of the plateau farther to the north (Hippolyte and Dumont, 2000), 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 S_c=0.6 m m⁻¹) 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 (Mudd, 2017). 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 S_c=0.5 m m⁻¹, 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 10 Basin averaged R^∗ vs. E^∗ plot (Roering et al., 2007; Grieve et al., 2016b). Two values for the critical hillslope gradient S_c are tested. Open circles correspond to S_c=0.6 m m⁻¹, colored according to the position of the basins along the transect (Figs. 2 and 7). Small dark filled circles correspond to S_c=0.5 m m⁻¹. Pale yellow symbols are averages over hilltop patches, for S_c=0.6 m m⁻¹ (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 (Okada, 1985), 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 Curnelle, 1978). (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 Curnelle, 1978), 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 Curnelle, 1978), 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 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, 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 Whipple, 2012), do not display as clear a pattern as hillslope metrics (Fig. 7d). Normalized steepness index values (k_sn) 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 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 (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 Royden, 2012), 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 k_sn across the inferred rock uplift gradient. We note, however, that in such a situation, stream concavity (m∕n) should display some sensitivity to changes in rock uplift (Wickert and Schildgen, 2019), 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, L_H, remains nearly constant across the transect between 140 and 160 m, whereas hilltop curvature, C_HT, 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 L_H and C_HT 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⁻¹ range corresponding to our CRN data (Fig. 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. (2008b, 2009) as the ratio between characteristic fluvial and hillslope timescales:

$$\begin{array}{}\text{(14)}& P_e={\displaystyle \frac{K{l}^{\mathrm{2}(m+\mathrm{1})-n}}{D{\mathit{\zeta }}^{\mathrm{1}-n}}},\end{array}$$

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 P_e=1 allows us to retrieve l_c 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 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 pieces of evidence 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 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 ¹⁰Be 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) ¹⁰Be 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⁻¹ 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 C_HT in Fig. 7b, because sampling sites where selected on relatively high-curvature ridges, higher than 0.02 m⁻¹ in most cases, whereas spatially averaged values along the transect are on average lower than 0.02 m⁻¹ (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 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² a⁻¹ 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 C_HT pattern along the transect into denudation rates using Eq. (3). 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 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 C_HT. The estimated median denudation values are ∼ 30 and ∼ 50 mm kyr⁻¹ 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 C_HT at these sites are ≥0.02 m⁻¹ (Fig. 12a), whereas continuously averaged C_HT along the transect are on average below 0.02 m⁻¹ (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⁻¹ 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⁻¹, which implies that slip rate is most likely < 50 mm kyr⁻¹ 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⁻¹. For example, on the Trévaresse 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 kyr⁻¹ range (Chardon et al., 2005), whereas accumulated deformation since the late Miocene indicates a slip rate of 30 ± 20 mm kyr⁻¹ (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 kyr⁻¹ 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.