Seasonal dynamics of marly badlands illustrated by field records of hillslope regolith properties, Draix – Bléone Critical Zone Observatory, South-East France

Sparsely vegetated badlands are loci of intense erosion that is sufficiently rapid to have observable effects on human timescales. Characterizing and understanding the physical weathering processes in these settings are key to predicting the temporal variability of regolith production and sediment flux, as well as their evolution under changing climate conditions. Here, we study intra-annual changes of hillslope properties and explore the relationship between sediment production and transport in steep marly badland catchments of the Draix – Bléone Critical Zone Observatory (SE France), where decades-long monitoring records show rapid morphological changes. There is evidence for seasonal dynamics of these badlands, but characterization and quantification of physical weathering processes have been lacking up to now. We explore this gap by monitoring key regolith parameters including grain-size distribution (characterized by D50), surface resistance and water content in the regolith layer (surface to (cid:1) 10 cm depth) at different locations, through repeated field surveys over a 2.5-year period. While water content appears to be directly controlled by the last previous rainfall event, the cyclic annual pattern in observed D50 suggests that loose and finely fragmented regolith is mainly produced and accumulates during the winter season, whereas sediment transport is dominant during spring – summer. This dynamic reduces regolith thickness and induces coarsening of hillslope surface material between early spring and autumn. Principal component analysis (PCA) highlights the strong correlation between resistance and D50. We therefore suggest that D50 provides the best proxy of regolith weathering in these marls. The spatial variability of the regolith was analysed through a geophysical profile, highlighting distinct behaviour depending on slope aspect. However, the distribution of slope angles is independent of aspect at the catchment scale. These results corroborate the strong annual dynamics of these catchments, where hillslopes and gullies are drained during spring and early summer high-intensity precipitation events, inducing high sediment yields.


| INTRODUCTION
Soil physical properties are a critical component for most ecosystems and result from the interaction between multiple erosion and weathering processes. As the interface between bedrock, biota and air, the soil is by nature one of the main elements of the Critical Zone and records its complex interactions (e.g., Brantley et al., 2011;Milodowski et al., 2015). Here, we define soil as the weathered unconsolidated surface layer that is distinguishable from fresh rock as it is exposed to erosion and weathering processes. Since biotic and organic matter processes are not an important component of soils in denuded badlands, we consider 'soil' and 'regolith' as synonyms for our study.
This layer contains the mobile products of bedrock transformation that can be transported into the stream network. The soil evolves on short timescales and modulates the variability of landscape morphology (e.g., ). However, the diversity of soils and climates renders the relationship between soil characteristics and weathering mechanisms complex. Climate change may also affect the spatial and temporal patterns of weathering and soil production, thereby potentially perturbing sediment flux to rivers.
Determining and quantifying parameters that best describe the dynamics of the regolith in different geological contexts and under diverse climates remains a challenge. Physical soil properties can be approached by different viewpoints and need to be characterized carefully to serve the specific objectives of their study in different domains.
For agricultural purposes, the structure (packing, fragmentation, porosity, etc.) of the soil affects its fertility (Ajdadi et al., 2016;Troeh & Thompson, 2005). In this domain, surface roughness is used as an indicator of soil quality and appears to be an important property because of its control on infiltration and runoff, which affect the susceptibility of the soil to erosion (Bullard et al., 2018;Darboux et al., 2002;Vázquez-Tarrío et al., 2017). Surface roughness is also closely related to grain size, which is an essential parameter controlling sediment transport both on hillslopes and in rivers (Haddad et al., 2022;Howard, 1994;Piton et al., 2018). In hydrology, the hydraulic capacity, porosity, water content and infiltration rate are key parameters to predict water storage and overland flow. Regolith properties such as density, mechanical properties or roughness can evolve seasonally (Nadal-Romero et al., 2007;Regüés et al., 1995) and influence infiltration and detachment rates-hence the hydrosedimentary response of a catchment (Nadal-Romero & Regüés, 2010;Regüés et al., 1995).
Soil properties and transport processes on hillslopes are also important in the study of mountain hazards, as their characterization is required for understanding flash-flood and debris-flow triggering, landsliding or post-fire landscape response (Nasirzadehdizaji & Akyuz, 2022;Rengers et al., 2019). Finally, physical soil properties can affect chemical weathering, for instance by modifying porosity and thus increasing the weathering surface, or by affecting vegetation development. These interactions, initiated by the modification of soil properties, can lead to changes in the CO 2 flux between the atmosphere and soils (Soulet et al., 2021), as well as to vegetation segregation (Regüés et al., 2000). Associated with these varied fields of interest in soil properties, a wide range of techniques have been developed to characterize the regolith, ranging from laboratory analysis of field samples to highresolution imagery of the soil (Bechet et al., 2015;Loye et al., 2016;Schoeneberger et al., 2012;Tarolli, 2014) and geophysical tools to image the subsurface (Tabbagh et al., 2013;Tetegan et al., 2012).
The sparsely vegetated and strongly dissected landscape known as 'Terres Noires' in the southeastern French Alps is characterized by rapid morphological changes and interactions between hillslope and fluvial processes on short (intra-annual) timescales. Over the last 35 years, several small catchments in these marly badlands have been monitored in the framework of the Draix-Bléone Critical Zone Observatory (CZO), leading to enhanced knowledge of weathering mechanisms (Rovéra & Robert, 2005), sediment export (e.g., Antoine et al., 1995), and the relation between sediment production and export (e.g., Ariagno et al., 2022). The mountainous climate and easily erodible lithology induce a strong seasonal pattern in regolith formation and sediment export. During winter, regolith is formed mainly due to frost-cracking processes (Ariagno et al., 2022), leading to a heterogeneous surface and desiccation features, as also observed in similar Spanish badlands (Nadal-Romero & Regüés, 2010;Regüés et al., 1995). Spring flash floods induce a peak in sediment export (Ariagno et al., 2022;Mathys, 2006). Because of the ample availability of sediment and the high network connectivity (Jantzi et al., 2017), the Draix-Bléone catchments record some of the highest observed specific sediment yields worldwide, with an annual value around 14,000 t/km 2 /y for the Laval catchment (0.86 km 2 ) and sediment concentrations that can reach 800 g/L (Draix-Bléone Observatory, 2015; Le Bouteiller et al., 2021). In this context, improving the understanding of the physical processes of soil production at the surface, as well as completing the characterization of the regolith and its temporal variability, are required to quantify sediment availability on hillslopes and more accurately predict the amount of sediment that will be transmitted to downstream fluvial areas.
The aims of this study are to synthesize the spatial and temporal variability of physical soil properties, and to identify a quantitative proxy representing the regolith weathering state in the Draix-Bléone catchments. Our main hypothesis is that physical properties of the regolith vary in response to weathering and erosion processes, and can therefore be used as a proxy to track these processes. To achieve this, we combine local (plot-scale) repeat field measurements over a 2.5-year period with a hillslope-scale geophysical survey. Our approach builds on similar work in subhumid badlands in the Central Pyrenees (Nadal-Romero & Regüés, 2010;Nadal-Romero et al., 2007) that studied weathering using field measurements of bulk density, regolith moisture and surface mechanical resistance. These studies observed seasonal trends in the measured soil properties and highlighted slope aspect as a determining factor for weathering processes. Focusing on physical mechanisms, our study also addresses contrasts in slope aspect but includes grain-size measurements in the field dataset. Grain size depends primarily on the processes and degree of bedrock weathering (e.g., Neely & DiBiase, 2020). Although multiple techniques have been developed to measure grain size on hillslopes (Harvey et al., 2022), fully capturing this parameter remains complex. However, grain size is a key component in investigating sediment entrainment and motion on hillslopes (Michaelides & Martin, 2012;Sklar et al., 2017). Cohen et al. (2010) explored depth-dependent weathering functions using the grain-size variable D50 as a witness of change on hillslopes. Their modelling approach distinguished different weathering processes: 'soil-weathering' (i.e., the breakdown of large soil particles to smaller soil particles within the soil layer) and 'soil-production' (i.e., bedrock weathering into soil). Soil-production functions such as the one of Heimsath et al. (1997) are commonly used in landscape evolution models but the parameters of these laws are difficult to calibrate because of the lack of field data and supposed constant in time. Our work adds local-scale information on hillslope processes, providing qualitative and quantitative data about the temporal variability of soil production. Model calibration based on field measurements of weathering processes should improve the quality of estimates of erosional outputs and thus improve model predictions under different climate scenarios.
The specific goals of this study are: (1) to analyse the seasonal evolution of each measured physical soil property and investigate potential relationships between these variables; (2) to identify a proxy regolith property that correlates with the observed intra-annual hysteresis in sediment export (Ariagno et al., 2022) and to highlight the weathering processes associated with this parameter; (3) to assess the spatial variability in regolith thickness and structure using seismic measurements, and to discuss the relationship between aspect, slope and weathering processes on hillslopes; and (4) to calibrate models of future evolution and their implications for predictions of sediment dynamics in the marly catchments of the Draix-Bléone CZO.

| Study site
The Draix-Bléone CZO is part of the French network for the study of the critical zone OZCAR (Gaillardet et al., 2018) and is specifically dedicated to the study of erosion and sediment transport in a mountainous region. Hydrological and sedimentary fluxes have been monitored in several catchments of this CZO since 1983(Draix-Bléone Observatory, 2015. The Draix-Bléone CZO is situated in the Alpine foothills, at elevations between 850 and 925 m above sea level, 12 km northeast of the town of Digne-les-bains in southeast France (Figure 1). The Draix catchments are drained by the Bouinenc stream, a tributary of the Bléone, which is itself one of the main tributaries of the Durance River. The bedrock geology is characterized by Mesozoic marls and limestones that were folded and faulted during the Alpine orogeny (Lemoine et al., 1986;Lickorish & Ford, 1998). The Bouinenc catchment has a mountainous and Mediterranean climate due to its relatively high elevation (>800 m above sea level). The mountainous influence is responsible for cold winters with frequent frost, whereas the Mediterranean influence leads to dry summers interspersed with intense thunderstorms. Annual rainfall is about 900 mm/yr (varying between  This study focuses on the Moulin catchment of the Draix-Bléone CZO (Figure 1), which has a drainage area of 0.10 km 2 , a vegetation cover of around 46%, and is entirely underlain by thick Middle Jurassic black marl series. The dominant bedding dip is about 35 to the east in the Draix area. The catchment morphology is characterized by high drainage density and deeply incised gullies; the catchment is wide in its upper part and becomes narrower and concentrated around the main meandering channel in its lower part. Sediment transport occurs through gravitational processes on hillslopes, minor landslides (<1 m 3 ) and debris flows in the upper drainage network, and fluvial transport as bedload and suspended load in the main drainage network. Considering only the unvegetated parts of the catchment as contributing to the sediment yield and the measured sediment density of 1700 kg/ m 3 , the average erosion rate in the Moulin catchment is around 8 mm/yr (Mathys, 2006).
Based on previous studies of regolith characteristics Mallet et al., 2018;Mathys, 2006), we chose specific locations for recording regolith properties over several years. Slope aspect and bedding orientation have been identified as important factors controlling hydrological processes (e.g., infiltration and runoff; Esteves et al., 2005;Mathys et al., 2005). We test the influence of these factors on regolith production and evolution by choosing four measurement plots of 1 m 2 each, all situated in the lower part of the catchment, encompassing a specific combination of these two factors: 1. P1: east-facing slope (36 ) and surface parallel to the bedding; 2. P2: west-facing slope (32 ) and surface perpendicular to the bedding; 3. P3: north-facing slope (54 ) perpendicular to the bedding; 4. P4: south-facing slope (49 ), perpendicular to the bedding. and therefore added another plot with a lower slope angle (P3bis; 51 ), which is placed in the meander upstream of the one containing P3 and P4.
The temporal analysis can be combined with quantitative measurements of the spatial variability in regolith thickness along a seismic profile. The geophysical transect is located between P3 and P4 through the meander ridge ( Figure 1).

| Plot measurements
Data were collected four times per year to record the seasonal evolution and the impact of the temperature and precipitation regimes on regolith characteristics. We started the field measurements in November 2019 and repeated the campaigns in February/March, April and August for 2.5 years until March 2022 (a total of 10 measurement campaigns, with one additional campaign for grain-size measurements; Supporting Information Table S1).
For each of the four plots and each campaign, we measured the soil resistivity in the field and collected samples for water content, density and grain-size measurements in the lab. The small-scale variability in water content and density was characterized by taking three different samples on each plot during each campaign, using metallic cylinders of 2.5 cm depth and 5 cm diameter. Because of the dryness of the marl surface, it was frequently difficult to catch the samples in the cylinder. Setting the cylinder into the soil was also challenging for several campaigns (especially in spring and summer). Cylinder samples were stored in a watertight plastic bag for transport and weighed before being placed in an oven at 70 C.
They were weighed again 3 and 6 days later to ensure obtaining a steady dry weight.
Samples for grain-size measurements were collected from the upper 2.5 cm of soil (i.e., a similar height to the cylinder) at undisturbed locations adjacent to the plots. We used a trowel to avoid breaking particles and collected samples in a single move. Samples were stored in aluminium baskets for transport, dried for 2 days in the oven and weighed. Samples were subsequently sieved into eight size fractions: 8, 4, 2, 1, 0.4, 0.2 and 0.08 mm. Sieving was performed under water by hand, without shaking the sieves to avoid fragmentation of the grains. All sieves, each containing a part of the sample, were dried in the oven for 2 days and weighed; the grain-size distribution was obtained by subtracting the weight of each sieve from the measurements. Because it was difficult to capture the fine fraction of the samples (<0.08 mm), we subtracted the sum of each dry fraction from the initial dry mass to obtain the weight of the fine fraction. All samples were sieved by the same operator to avoid uncertainty from the sieving process as much as possible. In order to assess the reproducibility of our grain-size measurements we sampled plot P1 five times in May 2020. We also collected a vertical profile for grain-size analysis at each plot during the May 2020 campaign. To this end, we collected two or three cylinder samples at each depth. After each sample collection at a particular depth, we homogenized and cleaned the new surface (previous surface À thickness of the cylinder) before collecting the following sample.
A mini-penetrometer was used to obtain resistance measurements from each plot. We randomly collected 20 measurements in each plot for each campaign to characterize the local-scale spatial variability in resistance.
The cumulative rainfall in the last 10 days before each campaign was obtained by summing the detailed event records, measured with a tipping-bucket rain gauge located at the outlet of the neighbouring Laval catchment, <500 m from the Moulin outlet ( Figure 1).

| Geophysical survey
We collected geophysical data in December 2018 along a 24 m long profile crossing the ridge separating plots P3 and P4 ( Figure 2). We generated seismic waves in the subsurface with the impact of a hammer on a high-density polyethylene plate placed on the ground. We repeated these shots at 25 locations spaced 1 m apart along the transect. For each of these shots, we recorded the seismic-wave propagation using 96 4.5 Hz geophones spaced 0.25 m apart and connected to four 24-channel Geometrics Geode seismic recorders. Total station surveying was carried out to accurately retrieve the locations and elevations of shots and geophones along the profile.

| Data processing: Water content, density, grain size and resistance
We computed the volumetric water content (ϑ) using the measured cylinder weights: with m h and m d the mass of the humid and dry samples respectively, ρ s the dry sample density and ρ w water density (1 g cm À3 ).
The reported water content is the mean of the three samples collected for each plot/campaign. The uncertainty is the highest value between either the mean of the uncertainty calculated for each sample, or the standard deviation of the water content of the three samples.
From the same samples, we computed the dry sample density ρ s (g/cm 3 ): The value of D50 was computed by linear interpolation between the two cumulative masses (associated with sieve size) that bracket the D50 mass. The uncertainty associated with the D50 value is set by taking Δm d = 0.1 g for each sieve (implying Δm d_total = 0.8 g). From the five samples collected at the same plot and time we infer that the D50 value is reproducible to within 0.3 mm (Supporting Information Figure S1).
To obtain the resistance value for each plot/campaign we computed the mean of the 20 measurements. We assigned a value of 5 to measurements that were out of the mini-penetrometer bounds (maximum = 4.5 units). In winter and early spring, the regolith was sometimes crusty at the surface, with loose and fragmented material underneath. With the pressure on the crust, the mini-penetrometer pierced the crust suddenly; we therefore recorded the value required to pierce the crust. The uncertainty associated with this measurement is the standard deviation of the 20 measured values.
The cumulative rainfall (D-10) was calculated because we observed that volumetric water content appeared to be significantly influenced by short-term weather conditions preceding the campaigns. Several observation periods preceding the campaigns, varying between 1 and 15 days, were tested, and the 10-day period was chosen based on the best correlation with volumetric water content.
We performed statistical analyses on this dataset, including a multivariate linear regression, using the sklearn-PCA Python package.

| Geophysical and morphological data processing
Collected seismic data were processed using seismic-refraction tomography to estimate the P-wave velocity (Vp) in the subsurface.
We initially picked first arrival times on each of the 96 traces of the 25 recorded shot gatherers with a sufficiently high signal-to-noise ratio. These first arrival times were then inverted using the seismicrefraction tomography module of the pyGIMLi Python library (Rücker et al., 2017), where the inversion domain corresponds to a triangular mesh with constant velocity cells through which rays are traced using a shortest-path algorithm. The velocity in each cell is solved with a generalized Gauss-Newton inversion, which starts with an initial model consisting of a velocity field that increases linearly with depth.
The subsurface velocity distribution is updated at each iteration in order to decrease the discrepancy between estimated and observed travel times. We followed the inversion strategy proposed by Pasquet  (Vallauri, 1997). We used a vegetation cover shapefile to obtain a clipped DEM of bare parts of the catchment, which we used for the morphological analyses.

| RESULTS
3.1 | Relation between field measurements and state of the regolith Supporting Information Table S2  Multivariate linear regression analysis shows that around 55% of the variability in resistance is explained by D50 and water content (Supporting Information Table S3). Figure 4A shows the positive linear trend between resistance and D50. This figure also shows some variability between the different plots: P3 and P4 are situated in the extreme parts of the regression, whereas measurements from P1 and P2 are spread out along the trend.
A two-factor analysis of variance (ANOVA) test was run for each regolith property (density, water content, resistance and D50) to check for significant differences between seasons and between plots.
A significant effect of the season factor (p-value < 0.01) was detected for all properties. A significant effect of the plot factor was detected only for the D50 and resistance.
The temporal evolutions of the mean D50 and the mean resistance values present a similar sinusoidal trend ( Figure

| Geophysical results
The near-surface geophysical survey provides a complementary view of the regolith to the field measurements, as it has lower resolution but images a deeper layer over a larger area. Following Flinchum et al.

| Morphological analysis of the Moulin catchment
Using only the unvegetated area of the Moulin catchment (0.038 km 2 ), we compared the slope frequency for four slope-aspect classes (north-, east-, south-and west-facing; Figure 7).

| Uncertainties related to field measurements
The characteristics of the marly lithologies, including their softness and flaky grain shape, lead to uncertainties during both sampling and laboratory analyses. We aimed to estimate these uncertainties by taking multiple samples (cf. Sections 2 and 3). Controlling the cylindersample volume was particularly difficult in the non-cohesive slope material, as the strong preferred orientation of the marl grains induced an irregular surface at the base of the cylinder.
Grain-size samples incorporate aggregate grain sizes with up to centimetre-sized particles that disaggregate upon contact with water.
Our sample processing used a single operator and a similar procedure in order to reduce uncertainties associated with grain-size determinations. Sampling necessarily induced some fragmentation, the effect of which on the measured grain-size distribution is impossible to quantify. We therefore took care to sample the same quantities of regolith in the same manner during each campaign, both at the surface and at depth (cf. Section 2). Samples were not subsampled to avoid segregation bias. Results from the reproducibility assessment validate the protocol and quantify the uncertainties (Supporting Information Figure S1).
The resolution of geophysical measurements is closely linked to the spacing between geophones and the seismic wavelength (e.g., Scott et al., 2017). For seismic refraction analysis, resolution is highest at the surface and decreases with depth. Lateral resolution is about half the geophone spacing close to the surface (i.e., 0.125 m) and vertical resolution is about a tenth of the wavelength (i.e., less than 0.4 m in the shallow, low-velocity materials). Although mean saprolite depths on both south-and north-facing slopes are only slightly higher than the vertical resolution limit, they are sufficiently different to interpret them in terms of spatially variable weathering intensity.

| A proxy of regolith weathering
Looking at the dependency between measured field variables is a common analysis strategy, and several studies have noticed a positive correlation between water content and resistance or between density and D50 (e.g., R. S. Anderson et al., 2013;Nadal-Romero et al., 2007;Regüés et al., 1995). The PCA shows that $50% of the variability in the dataset is related to the resistance. Although this variable is associated with significant uncertainties, due to the sensitivity of the micro-penetrometer to small spatial variations in the regolith surface, it is well correlated with both D50 and water content, which are independent of each other (Figure 3). Water content itself is correlated with D-10 cumulative rainfall. As shown by the depth profiles, D50 can be related to the state of fragmentation of the regolith. Therefore, the resistance of the regolith is explained by its state of weathering (D50) and humidity (water content and cumulative rainfall). The multivariate regression confirms that D50 and water content have a significant impact on the resistance and explain more than half of its variability. As D50 shows a clear seasonal variability ( Figure 4C), the temporal variability of the dataset can be explained by seasonal weathering dynamics combined with high-frequency variability of the humidity, linked to precipitation. variation in resistance and D50 (Figure 4) may also shed light on the relationship between weathering and sediment export. Ariagno et al. (2022) have shown recently that frost-cracking during winter is the main control on sediment production in the Draix catchments. This study also highlighted an annual hysteresis cycle in sediment export from these catchments, with a peak in sediment export for a moderate amount of rainfall in June, contrasting with a peak in cumulative rainfall linked to moderate sediment export during autumn. The sinusoidal trends observed in D50 and resistance provide physical support for these inferences: the regolith is more fragmented, finer and less resistant during the winter; spring rainfall events will thus easily mobilize this regolith and export it as sediment in a transport-limited regime.
As the easily mobilizable top layer of regolith is removed, the surface becomes coarser and more resistant, as indicated by our measurements at the end of the summer. This situation corresponds to a supply-limited regime, in which sediments are more difficult to detach and transport. However, the transition between both regimes is not sharp, as illustrated by comparing our D50 measurements over time with the depth profiles (Figures 4 and 5): the surficial regolith at the end of summer is still significantly finer grained than the material at depth. It is therefore weathered material that corresponds to the removal of <2.5 cm from the surface (depth of the sample cylinder), consistent with estimates of sediment export from the catchment (Ariagno et al., 2022). The timing of the transition between these two regimes is also controlled by the distribution of rainfall events over spring and summer.
Moreover, two anomalies need to be highlighted and discussed: (1) the particularly high D50 and resistivity values of the P3 samples; and (2) the anomalous measurements during winter 2020 (Feb20) compared to the following winters (Feb21 and, for D50, Feb22).
In comparison to plots P1, P2 and P4, which show a clear seasonality in D50 and resistance (Figure 4), the P3 plot shows less variability, a less predictive pattern in changes and overall higher values of resistance and D50. In contrast to the other plots, this slope also shows significant sediment deposition at its base. A possible explanation for these anomalies could be the slope angle of $54 , which is steeper than that of the three other plots and may prevent sediment accumulation on this hillslope. Thus, the lack of regolith could be due to dry-ravel transport operating on this relatively steep slope. However, regolith formation has been observed elsewhere on hillslopes with similar aspect and slope angles. An additional explanation could be that, in addition to the steep slope on the P3 plot, steeper slopes above it might contribute to enhancing dry ravel at this location and continuously remove any regolith produced. A reduced intensity of weathering processes on this north-facing slope could be an alternative explanation. However, we have no field arguments to support this hypothesis and it goes against previous studies showing that weathering processes are enhanced on north-facing slopes (S. P. Anderson et al., 2014;Nadal-Romero et al., 2007).
To avoid the potential specific characteristics of this site, we added a second sample on a north-facing slope (P3bis) to the dataset from November 2020. Data from this plot are much more coherent with those from plots 1, 2, and 4, confirming that weathering processes are not significantly different on north-facing slopes. In light of these observations, we choose not to take into account the records of P3 in the calculated mean values but to include data from P3bis when available. values of D50 and resistance at the end (Feb20) compared to the beginning (Nov19) of the winter season is interpreted as a lack of weathering of the regolith. As frost-cracking has been inferred to be the main physical weathering process in this catchment (Ariagno et al., 2022), we computed the frost-cracking intensity indicator (FCI) for both winters, following the approach of Ariagno et al. (2022).
Between the campaigns of November 2019 and February 2020, we found an FCI of about À11 C min/cm on south-facing slopes, whereas for the same period during winter 2020-2021 and 2021-2022 the FCI is about À411 C min/cm and À866 C min/cm, respectively. These large differences imply significantly less intense frost-weathering processes in the first winter of our campaign (Nov19-Feb20) compared to the other two.
Moreover, because precipitation above a threshold is the main driver of erosion in these badlands (Ariagno et al., 2022;Mathys, 2006), we also looked at the distribution of rainfall over these two winters. In the 2019-2020 winter season, nine rainfall events (as defined by Mathys, 2006) were recorded, with a cumulative rainfall of 278 mm. In the following winter (2020-2021), seven events were recorded, with a cumulative rainfall of 113 mm. These records suggest that, in addition to reduced sediment production, the 2019-2020 winter was also characterized by significant rainfall, which would have driven slope erosion during this winter season. Together, these observations explain the anomalous value of D50 found between November 2019 and February 2020: limited frost-weathering combined with enhanced entrainment of regolith by rainfall induced an absence of regolith fining during this season.

| Can we identify a weathering front in black marls?
In the denuded marly badlands of Draix-Bléone, where 8 mm/yr of erosion is recorded and a supply-limited regime occurs during part of the year (Ariagno et al., 2022;Bechet et al., 2016), estimating the regolith thickness is determinant to predicting sediment export.
We observed a rapid increase of D50 in the first centimetres below the surface, indicating rapidly diminishing fragmentation of the regolith with depth. This trend can be interpreted in two different ways.
The trend of D50 values is consistent with an exponential increase in weathering intensity towards the surface. Taking D50 as a proxy, such a trend can be expressed as where D inf is the D50 value of unweathered rock [mm], h is depth in the soil [mm] and h* is the characteristic soil decay depth [mm].
The exponential decline of weathering intensity with depth ('soil production function') has been well documented (e.g., Gilbert, 1877;Heimsath et al., 1997) but there are fewer studies quantifying soil-weathering rates (Wells et al., 2006;Yoo & Mudd, 2008). Cohen et al. (2010) explored the relation between the bedrock-weathering (soilproduction) function and the soil-weathering function. They showed that the soil-weathering function controls the distribution of the soil grading down the profile, but that soil-surface properties were controlled by the integrated degree of weathering and thus by the soilproduction function.
However, our results can also be interpreted as a linear increase in D50 with depth in the top 7 cm of the soil. This depth would indicate a boundary below which the regolith is characterized by another weathering state, illustrated by a constant value of D50. From a penetrography survey, Maquaire et al. (2002) also found a first boundary between 6 and 12 cm depth, depending on the location of the resistance profile. The impossibility to collect grain-size samples below 7 cm depth in the P2 and P4 plots can also be interpreted as indicating a resistance boundary at that depth. In both cases, these results confirm the presence of a fragmented regolith at the surface at the end of winter/early spring (i.e., in May 2020) with a first weathering front around 7 cm depth.
Additionally, the geophysical survey provides continuous data and images a deeper weathering front, whereas the field measure- In a more temperate climate, the Gordon Gulch watershed (Colorado) has been well studied, and contrasts in vegetation, rock weathering and moisture content have been found between north-and south-facing slopes, with enhanced weathering on north-facing slopes subject to long periods of snow cover (S. P. Anderson et al., 2014).
Such an opposite trend could be related to differences in climate and insolation conditions (e.g., Pelletier et al., 2018;Smith & Bookhagen, 2021) that would deserve further investigations.
Previously, we developed arguments that could explain the particularly thin regolith layer on the P3 plot (Section 4.2). However, our results between P4 and P3bis are still opposite to previous studies, although the distinct behaviour between aspects is reduced. Further field measurements on other north-facing slopes would be needed to compare our study site with other locations. These aspect distinctions are also found at a larger scale, as demonstrated by the geophysical survey: the weathering front of the saprolite is deeper on the southfacing (MST = 1.2 m) than on the north-facing (MST = 0.6 m) slope at this site. Although we know that the north-facing slope site is probably influenced by its high slope angle, this complementary continuous method confirms the difference in regolith and saprolite thickness between these two slopes of opposite aspect.
To extend the spatial analysis from the hillslope to the catchment scale and put the observed contrast between opposite-facing slopes into context, we analysed the morphology of the Moulin catchment.
Within the Moulin catchment ( Figure 1B), the main drainage is not located symmetrically within the catchment; a wider vegetated area is present on the north side of the channel (i.e., dominated by southfacing slopes). This distribution of vegetation is consistent with the observations of Regüés et al. (2000) and Nadal-Romero and Regüés (2010) that vegetation was more stable and persistent on south-facing slopes. However, the distribution of slope angles throughout the catchment does not highlight a significant difference between slope aspects ( Figure 6). Thus, the morphology of the ridge that we sampled in a meander of the Moulin catchment, which was clearly steeper on its north-facing than on its south-facing slope, cannot be related to a more general morphological trend over the catchment (Supporting Information Figure S3). 4.5 | Perspectives: Integrating regolith characteristics in landscape morphology analysis Although our field measurements show consistent trends, match with the observed sediment-export cycle ( Figure 7A) and confirm the results of previous studies on weathering processes in these badlands, field observations need to be continued and extended to better understand which processes drive the spatial heterogeneity of the regolith and what are the consequences of this heterogeneity.
Increasing the measurement frequency, although time consuming, could also help tracking how regolith properties react to event-scale precipitation and erosion.
First, to more fully assess the potential contrast between weathering on north-and south-facing slopes, a larger dataset of D50 on north-facing slopes would be required. Although slope angle is undoubtedly one of these key parameters, slope aspect may also play a role by influencing soil temperature, humidity or exposure to solar energy.
Tracking the location of sediment deposits could also help to explain the observed asymmetry of the position of the main channel in the Moulin catchment ( Figure 1B) as well as the adjacent Laval catchment, with the main channel of both catchments clearly shifted toward the south. Johnstone et al. (2017) showed that 'the asymmetry in the position of the thalweg within mainstream valleys is positively correlated with asymmetry in the length of fans'. Further study of these catchments is required to establish a potential relationship with sediment production and transport, or vegetation development, as a function of slope aspect.
At smaller scales, Bechet et al. (2015Bechet et al. ( , 2016 investigated the change in hillslope and regolith characteristic using terrestrial laser scanning. Similar studies could improve the quantification of regolith surface roughness over time, which could provide another proxy of the weathering state.
The results of the geophysical survey complement our field measurements and provide support for active erosional removal of regolith on the hillslope where the P3 plot is located, leading to a thinner saprolite layer. Performing additional surveys on longer and steeper slopes could further clarify the impact of aspect on regolith thickness.
In particular, it would be interesting to confirm the particular behaviour of the P3 site by looking at other north-facing slopes, in order to test whether slope angle or aspect is the prime parameter affecting regolith thickness.

| Implications for modelling landscape evolution in badlands
With the rapid growth of digital elevation models and numerical analysis tools, landscape-evolution models have become key tools to investigate the drivers of spatial and temporal variability in erosion rates, combine and test hypotheses and explore quantitative predictions (e.g., Tucker & Hancock, 2010). The development of numerical models of badland evolution provides significant promise for a better understanding of the rapid erosion processes operating in them, the interaction with vegetation and the long-term evolution of erosion rates (e.g., Gallart et al., 2013

| CONCLUSIONS
Based on our field measurements and statistical analyses, and accounting for the inevitable uncertainties in our dataset, we show that D50 is a good proxy to characterize the state of regolith weathering in the Draix badlands.
Both soil resistance and D50 follow seasonal weathering dynamics: small values of these parameters are found during the winter, when the regolith is particularly fragmented and crumbly, while high values are recorded at the end of the summer, when mobile regolith has been washed from the hillslopes by the intense spring/early summer rainfall events. This cyclicity of D50 is consistent with the intra-annual hysteresis in sediment export observed by Ariagno et al. (2022). A small D50 value implies availability of easily detached sediment on hillslopes for transport in the drainage network. In contrast, the highest values imply a less-weathered regolith that was exposed to the surface due to the removal of the more mobile sediment above it. Moreover, the frost-weathering intensity indicator, shown to predict sediment production in the Draix catchments (Ariagno et al., 2022), can explain field measurements of D50 that do not conform to the seasonal trend. In contrast to the seasonal cyclicity, the dynamics of the water content show a high-frequency (daily) variability. Independent of D50, this variable is highly correlated to the cumulative rainfall in the previous 10 days.
Although further studies on aspect influence are needed to confirm the results obtained here, both discrete (field samples) and continuous (geophysical survey) measurements highlight the spatial variability of the weathered layer thickness and show a distinct behaviour of the regolith between the south-facing and north-facing slopes studied. However, this trend is not visible at the catchment scale, as the morphological analysis of the denuded area of the Moulin catchment shows a similar slope distribution independent of aspect.
Our field measurements over a 2

DATA AVAILABILITY STATEMENT
The cumulative rainfall data used in this paper are available on the BDOH data portal at https://bdoh.irstea.fr/DRAIX/ and referenced under https://doi.org/10.17180/OBS.DRAIX (Draix-Bléone Observatory, 2015). All the plot data collected for this study is available in the Supplementary Material and will be uploaded to an online database upon acceptance of the manuscript.