# Morphotectonic Evolution of an Alluvial Fan: Results of a Joint Analog and Numerical Modeling Approach

Clément Garcia-Estève, Yannick Caniven, Rodolphe Cattin, Stéphane

Dominguez, Romain Sylvain

► To cite this version:

Clément Garcia-Estève, Yannick Caniven, Rodolphe Cattin, Stéphane Dominguez, Romain Sylvain. Morphotectonic Evolution of an Alluvial Fan: Results of a Joint Analog and Numerical Modeling Approach. Geosciences, 2021, 11 (10), pp.412. 10.3390/geosciences11100412. hal-03761787

HAL Id: hal-03761787

https://hal.science/hal-03761787v1

Submitted on 26 Aug 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Article

# Morphotectonic Evolution of an Alluvial Fan: Results of a Joint Analog and Numerical Modeling Approach

Clément Garcia-Estève *, Yannick Caniven, Rodolphe Cattin, Stéphane Dominguez and Romain Sylvain

<figure><img src="image_2.png" /></figure>

check for
updates

**Citation:** Garcia-Estève, C.; Caniven, Y.; Cattin, R.; Dominguez, S.; Sylvain, R. Morphotectonic Evolution of an Alluvial Fan: Results of a Joint Analog and Numerical Modeling Approach. *Geosciences* **2021**, *11*, 412. https://doi.org/10.3390/geosciences11100412

Academic Editors: Nicolas Loget, Julien Babault and Jesus Martinez-Frias

Received: 30 June 2021
Accepted: 21 September 2021
Published: 1 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Géosciences Montpellier, Université de Montpellier, CNRS, 34090 Montpellier, France; yannick.caniven@umontpellier.fr (Y.C.); rodolphe.cattin@umontpellier.fr (R.C.); stephane.dominguez@umontpellier.fr (S.D.); romain.sylvain@etu.umontpellier.fr (R.S.)

* Correspondence: clement.garcia-esteve@umontpellier.fr

**Abstract:** Surface topography results from complex couplings and feedbacks between tectonics and surface processes. We combine analog and numerical modeling, sharing similar geometry and boundary conditions, to assess the topographic evolution of an alluvial fan crossed by an active thrust fault. This joint approach allows the calibration of critical parameters constraining the river deposition-incision laws, such as the settling velocity of suspended sediments, the bed-rock erodibility, or the slope exponent. Comparing analog and numerical models reveals a slope-dependent threshold process, where a critical slope of ca. 0.081 controls the temporal evolution of the drainage network. We only evidence minor topographic differences between stable and stick-slip fault behavior localized along the fault scarp. Although this topographic signature may increase with the slip rate and the return period of slip events, it remains slight compared to the cumulated displacement along the fault scarp. Our results demonstrate that the study of morphology cannot be used alone to study the slip mode of active faults but can be a valuable tool complementing stratigraphic and geodetic observations. In contrast, we underline the significant signature of the distance between the fault and the sediment source, which controls the degree of channels incision and the density of the drainage network.

**Keywords:** alluvial fan; analog and numerical modeling; inversion; fault slip modes

## 1. Introduction

Geomorphological approaches consider the well-established principle that landscape morphology is controlled by the interactions between tectonics, sedimentation, and erosion [1–3]. During the last decades, many field methods have been developed to characterize these interactions by documenting surface processes such as soil transport, river incision, or basin-wide denudation [4]. Local landforms, such as fault scarps, offset alluvial fans, or deformed river terraces are associated with well-known geomorphic processes. Thanks to cosmogenic isotopes, processes such as incision or denudation are quantified [1]. They are commonly used to estimate a mean fault slip rate or an average incision rate (e.g., [5–8]). However, most erosion parameters remain poorly constrained by field observations due to their multivariate and interrelated dependencies on climate, lithology, and vegetation cover (e.g., [9,10]). At the same time, analog and numerical models have been built to explore erosion laws and the associated parameters (e.g., [11–17]). They underline the complexity of geomorphic processes, including non-linearity, thresholds, feedbacks, and equifinality.

Comparing modeling results with natural landscapes is challenging. A few previous studies combined analog and numerical approaches for studying the surface and tectonic processes at the crustal scale (e.g., [18]). This method is currently used to estimate erosion and sedimentation processes according to the granular composition [19] and the erosion influence on tectonics in a context of crustal shortening [20]. Following this approach,
we investigated the evolution of a geomorphologic marker by comparing analog and numerical models.

In this paper, we study the morphotectonic evolution of an alluvial fan crossed by an active reverse fault. We focus on a short time scale of 100 ka to explore poorly studied mechanisms, such as the topographic signature of the fault slip mode. We choose an alluvial fan because this geomorphological object has a high deposition rate and forms rapidly, which corresponds to the study's time scale. Moreover, it is simple to characterize in terms of its dynamics and geometry. The context of thrust fault allows characterizing several markers, such as the fault scarp and syntectonic deposits. Together, these markers in the morphology are used to study the morphological signature of the fault slip behavior. Last, this setting associated with an alluvial fan crossed by a thrust fault is frequently observed in the field. Here, we assume active tectonic with rapid surface processes, which can be compared to those observed along the Tian Shan's or the Himalayas' front.

After defining the experimental setup and the numerical methods, we describe the general assumptions of the approach. Next, we present our results, devoted to focusing on (1) the repeatability of the laboratory experiments; (2) the calibration of the numerical simulation setup from the analog results; and (3) the analysis of the general features of the modeled morphotectonic evolutions.

We paid specific attention to the fault slip mode associated with stick-slip or stable behavior. This parameter is critical for improving seismic hazard assessment. The most commonly used approaches focus on seismic coupling calculation derived from geodetic measurements (e.g., [21]). Weak coupled areas are associated with stable sliding areas of low-stress accumulation and, therefore, potentially less prone to trigger major earthquakes. However, the main limitation of this approach is the short time coverage of observations, which only document a few decades. Another method is the paleoseismology approach conducted by performing trenches and outcrops analyses (e.g., [8,22,23]). The observations then come from stratigraphy, lithology, grain size, and dating. This method provides critical information on the seismic potential by completing the regional seismicity catalogs with the inferred paleo-earthquakes. However, although these paleoseismology studies are critical to collecting missing data necessary to assess the seismic hazard, they remain often limited by the local character of measurements and their time-consuming implementation. Here, we propose a complementary approach to assess the fault slip mode from the analysis of topography.

Finally, based on numerical simulations, we discuss the influence of other parameters, such as the fault slip rate, the return period of slip events, and the fault location relative to the upstream outlet, on the morphological signature of the fault slip behavior.

## 2. Methods

We developed a joint approach coupling analog and numerical models to address the impact of surface and tectonic processes on topography. Although these models have their specificities associated with their implementation, particular attention was paid to sharing an identical geometry, initial and boundary conditions, and similar external and internal forcings.

### 2.1. Analog Modeling

#### 2.1.1. Experimental Setup

Graveleau et al. (2008, 2011) [12,13] initially developed the experimental device to perform morphodynamic experiments at the scale of an active foreland. Here, we adapted this device to investigate the morphologic evolution of a simpler model with a single alluvial fan fed by a single upstream outlet (Figure 1). The model comprised three main parts from upstream to downstream:

1. A reservoir, equipped with rotating blades, containing a mix of granular material and water.

2. A transfer zone used to reduce and control the flow.
3. A sediment deposit area where the fan grows.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Experimental device used to perform the analog models. (a) General picture of the setup. On the left, the sediments go through a regulator flux then are weighed on the scale. On the right, the alluvial fan is deposed on an area of 60 × 60 cm. The alluvial fan boundary is outlined by an orange line, and the fault by a red line. (b) Schematic diagram of the experimental setup showing the circulation of sedimentary flux. The depositional area has a small slope of 0.8° and a reverse fault with a dip of 43°.</figcaption></figure>

A hand control sluice gate allowed the mix of sediments and water to get out of the reservoir into a 1 × 1 × 10 cm aluminum channel. Before the depositional area, the flow went through a mass control device used to monitor the flow transfer rate during the experiments (Figure S1). Next, the transfer zone enabled a slow water flow of 4 drops/s to avoid deposits that could restrain the sediment flow. Finally, the flow reached a 60 × 60 cm deposit area, tilted 0.8° downstream (Figure 2).

A reverse fault crossed the deposit area perpendicular to the stream direction with a dip of 43° upstream (Figure 2b). In our experiments, the fault trace was located 10 cm from the upstream outlet. The hanging wall compartment could be lifted by incremental displacements using a step-by-step computerized motor. This system could reproduce a stable or stick-slip fault behavior depending on the set size and duration of the rigid increments.

We performed five sets of experiments, each simulating a stable and a stick-slip fault behavior and sharing an average fault slip rate of 1 cm/h and a sedimentation rate of ca. 1 cm/h. This sedimentation rate was measured at the top of the alluvial fan formed during the phase with no tectonics. Stick-slip models were characterized by fault slip increments of 500 microns every 180 s (Table 1), whereas stable fault slip behavior was simulated with fault increments of 6.2 microns every 2.23 s.

### 2.1.2. Material Properties

The analog material composition was derived from the water-saturated mixture developed initially by Graveleau and Dominguez (2008) [12], refined by Graveleau et al. (2011) [13] and Strak et al. (2011) [17]. Its solid part comprised several granular materials with micrometric grain sizes such as plastic powder (PVC), silica powder, and micro-beads. The main interests of such a mixture are to reproduce natural deposit sorting processes thanks to contrasted densities and wide grain-size distribution, inducing differential transport distances for each material component. Indeed, it was shown that this type of material
is suitable to reproduce realistic geomorphological objects such as river channels, drainage networks, and alluvial fans, sharing strong similarities with their natural equivalents in size, shape, and evolution [10,14,24].

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Main geometric and tectonic model set-up. (a) General scheme of the device. The coordinate system is centered in the middle of the fault, given by the red line. The closed box boundaries are represented with hatching lines. Input flux is represented with two arrows: a blue arrow for the total flux (water and sediments) and an orange arrow for sedimentary flux only. (b) Close-up on the fault zone. The topography is represented with an orange line. $\alpha$ is the fault angle and the red arrow the fault slip direction $V_{tecto}$. (c) Temporal evolution of the fault displacement, including a fan building stage with no tectonics and a tectonic forcing stage. Two slip modes are modeled: stable sliding and stick-slip behavior. $i$ represents the increment of displacement per earthquake and $T$ the interseismic duration.</figcaption></figure>

<table><caption>**Table 1.** Scaling parameters between experimental and numerical reference models.</caption><thead><tr><th>Parameter</th><th>Experimental Model</th><th>Numerical Model</th></tr></thead><tbody><tr><td>Length</td><td>1 mm</td><td>10 m</td></tr><tr><td>Time</td><td>1 s</td><td>10 yr</td></tr><tr><td>L</td><td>60 cm</td><td>6 km</td></tr><tr><td>T<sub>fan</sub></td><td>1.5 h</td><td>54 ka</td></tr><tr><td>T<sub>tecto</sub></td><td>1 h</td><td>36 ka</td></tr><tr><td>v<sub>ref.</sub></td><td>1 cm.h<sup>-1</sup></td><td>2.8 mm.yr<sup>-1</sup></td></tr><tr><td>T<sub>ref.</sub></td><td>180 s</td><td>1.8 ka</td></tr><tr><td>q</td><td>1.4 × 10<sup>-4</sup> L.s<sup>-1</sup></td><td>140 × 10<sup>3</sup> m<sup>3</sup>.yr<sup>-1</sup></td></tr><tr><td>q<sub>s</sub></td><td>1.4 × 10<sup>-5</sup> L.s<sup>-1</sup></td><td>14 × 10<sup>3</sup> m<sup>3</sup>.yr<sup>-1</sup></td></tr><tr><td>S<sub>c</sub></td><td>0.3</td><td>0.3</td></tr></tbody></table>

The mixture used in our experiments was composed of 90 % water and 10% granular material, including 30% silica, 51% PVC, 6% anthracite, 10% glass micro-beads, and 3% pumice stone (see Supplementary Materials S2 for physical properties). Compared to the previous material, the main differences were an increase in the percentage of PVC, a decrease in micro-beads, and the addition of pumice powder. Such modifications were made to decrease the density of the material favoring a sediment transport consistent with a single kilometric alluvial fan.
### 2.1.3. Scaling

Since the pioneer paper of Hubbert (1937) [25], many studies have focused on analog models scaling (e.g., [19,20,26–28]). Here, the scaling of our experiments was done by taking into account dynamic, geometric, and kinematic similarity criteria between model parameters and field observations.

#### Dynamic Scaling

The scaling of morphotectonic models is complex to achieve. Different processes act at different time scales, leading to a wide variance from few mm.yr⁻¹ for tectonic slip rate and incision–deposition rate to few m/s for flow circulation. Moreover, water has a dual role in transport and erosion processes.

In our experiments, the tectonic deformation was concentrated along the fault scarp. Assuming that this deformation can be neglected, the dynamic similarity criterion is constrained by the hydraulic forces balance only [29,30]. This hypothesis can be tested by calculating the dimensionless Reynolds (Re) and Froude (Fr) numbers [30,31]:

$$Re = \frac{u \times l}{v} \quad \text{and} \quad Fr = \frac{u}{\sqrt{g \times l}} \qquad (1)$$

where *u* is the flow velocity, *l* is a characteristic length as the flow thickness, *v* is the fluid's kinematic viscosity, and *g* is the gravity acceleration.

While some portions of natural rivers are characterized by turbulent flows with a high Reynolds number [29], in analog experiments, the water flow is laminar with a low Reynolds number [13,16]. This discrepancy is mainly related to the use of water in analog experiments. Therefore, the analog models and nature share the same fluid's kinematic viscosity, *v*, and a similar flow velocity. The Reynolds numbers are then essentially controlled by *l* (see Equation (1)), which significantly differs from models to nature.

In our experiments, the flow velocity was estimated to be a few mm.s⁻¹–cm.s⁻¹, and the flow thickness to 0.1–1 mm. These values lead to a Froude number of $Fr = 0.02–0.6$, consistent with $Fr = 0.01–1$ obtained from natural rivers [32]. This good agreement indicates that the flow-inertia-to-gravity ratio is well respected in our experiment. Thus, although rigorous dynamic scaling is difficult and elusive to achieve, the conservation of the Froude number implies that the same force equilibrium controls the dynamic of both the analog models and nature.

#### Geometric Scaling

The spatial scaling factor $L^*$, defining the ratio between a characteristic length in the analog model and nature, is imposed by the size of the apparatus ($L = 60$ cm, Table 1) and the studied natural object dimensions (5–10 km). We obtained $L^*$ in the range of 6–12 × 10⁻⁵, i.e., 1 cm in the analog model upscaled to between 80 m and 170 m in nature. In this paper, to facilitate extrapolating model results to nature, we considered that 1 cm upscales to 100 m (Table 1).

#### Kinematic Scaling

To evaluate the time scaling, we assumed that accelerations of the simulated geomorphic processes are small enough to be neglected. The time ratio $T^*$ is then obtained by comparing average long-term velocities between the model and nature using the following relationship:

$$T^* = \frac{L^*}{V^*} \qquad (2)$$

where $V^*$ is the model-to-nature velocity ratio. Here, $V^*$ is defined from the average sedimentation rate. In our experiments, this rate is ca. 1–10 mm/h, whereas field observations range between 1 and 10 mm.yr⁻¹ [33], leading to a $V^*$ of ca. 10³–10⁵. These values and the previously calculated $L^*$ give a time ratio $T^*$ of 6 × 10⁻¹⁰–12 × 10⁻⁸. Thus, 1 s in the
model upscales to 0.3–50 years in nature. In this paper, we considered that 1 s upscales to 10 years (Table 1).

The time and velocity scaling factors were then used to constrain the fault slip rate in our models. As we intend to simulate a natural fault slip rate of 1–5 mm.yr⁻¹, we set the device for an analog fault slip rate of 1 cm.h⁻¹, equivalent to 2.8 mm.yr⁻¹, in our reference model (Table 1).

### 2.1.4. Model Monitoring and Data Processing

Two camera acquisition systems monitored the model:

I. a CCD Panasonic Gx80 camera located above the deposit area took one frame every 30 s. This acquisition gave a collection of 270 high-resolution pictures (4600 × 3500 px) for each experiment. These pictures were used to make a time-lapse video and analyze the fan building phase qualitatively (see movies in Supplementary Materials S3).

II. The final morphology was quantified by photogrammetry from 9 high-resolution pictures taken with a Sony Alpha7R2 camera (42 Mpx). After processing the data with Micmac [34] and CloudCompare software [35], we obtained a georeferenced digital elevation model (DEM) composed of about 10 million points with an average spatial resolution of 0.16 mm and an accuracy of 0.6 mm. The DEMs were analyzed using the library Matplotlib in Python.

### 2.2. Numerical Modeling

We also investigated evolution models numerically using the Landlab2.0 modeling toolkit [11,15]. To account for both surface processes and tectonics, we considered the following governing equation of mass conservation for a land surface elevation $z(x, y, t)$

$$ \frac{\partial z}{\partial t} = \left(\frac{\partial z}{\partial t}\right)_{\text{hill}} + \left(\frac{\partial z}{\partial t}\right)_{\text{fluv}} + v_z \quad (3) $$

where $\partial z/\partial t$ is the variation in elevation with time and $v_z$ is the tectonic uplift rate.

Hillslope processes modeling uses the TransportLengthHillslopeDiffuser component of Landlab, derived from the formulation proposed by Carretier et al. (2016) [36]:

$$ \left(\frac{\partial z}{\partial t}\right)_{\text{hill}} = q_s \times \frac{1 - \left(\frac{s}{S_c}\right)^2}{dx} - \kappa_{\text{hill}} S \quad (4) $$

where $\kappa_{hill}$ is an erodibility coefficient, $S = \partial z/\partial x$ is the local slope, $q_s$ is the sediment flux, $S_c$ is a critical slope, and $dx$ is the model node spacing. This law requires to know a priori two parameters: $S_c$ and $\kappa_{hill}$. A critical slope $S_c$ of 0.3 is set for all the numerical models to be consistent with the internal friction coefficient $\varphi = 20^{\circ} \pm 10^{\circ}$ of analog material (see Supplementary Materials S2). Estimating $\kappa_{hill}$ is less evident because this parameter exhibits a wide variance due to its dependency on climate, rock type, and vegetation. Soil transport efficiency is often described by a soil transport coefficient $D = \kappa_{hill} \times dx$, which is a diffusivity-like coefficient. In their recent compilation, Richardson et al. (2019) [10] propose a parameter $D$ ranging between $10^{-5}$ and $10^{-1}$ m²/yr. Here, we assumed a model node spacing $dx$ of 50 m. This value leads to a range for $\kappa_{hill}$ between $2 \times 10^{-7}$ and $2 \times 10^{-3}$ m.yr⁻¹.

Fluvial erosion and deposition were calculated from the Davy and Lague (2009) [37] formulation. A mass balance approach, taking into account transport of both topographic material and river sediment content, was coupled with an explicit representation of the sediment transport distance. We used the Erosion.Deposition component of Landlab [38]

$$ \left(\frac{\partial z}{\partial t}\right)_{\text{fluv}} = v_s \times \frac{q_s}{q} - \omega \quad (5) $$
where $v_s$ is the net settling velocity, and $q$ is the water surface flux. In their approach, Davy and Lague (2009) [37] defined a parameter $\xi = \frac{q}{d^* \times v_s}$, which represents the average travel distance of sediment grains from their erosion to their deposition. The parameter $d^*$ is the ratio between the mean sediment concentration in the stream and at the bed interface. They demonstrated that detachment-limited and transport-limited models can be encompassed using $\xi \to +\infty$ (i.e., small $v_s$) and $\xi \ll 1$ (i.e., large $v_s$), respectively. This study considered $v_s$ between $10^{-2}$ and $10^2$ m.yr$^1$, a range in agreement with the $\xi - q$ relationships tested by Davy and Lague (2009) [37].

The erosion $\omega$ was calculated from the classical power-law type equation

$$ \omega = \kappa_{river} q^m S^n \qquad (6) $$

$m$ and $n$ are two dimensionless exponents, which depend on catchment hydrology and the nature of the dominant erosional process [39,40]. The values of these two parameters are still debated. Based on the compilation presented by Harel et al. (2016) [9], we assumed a concavity index $m/n = 0.5$ and a slope exponent $n$ ranging between 0 and 4. The erodibility of the river bedrock $\kappa_{river}$ is far less constrained. This coefficient mostly depends on climate, lithology, vegetation cover, and the values of the exponents $m$ and $n$. In the following, we considered the range $10^{-5}-10^0$ according to Stock and Montgomery (1999) [41] and Harel et al. (2016) [9].

An uplift rate $v_z$ was imposed on the hanging wall's grid nodes to simulate surface tectonic displacements due to a reverse fault slip (Figure 2). In contrast to the experimental device, the elevation change due to horizontal tectonic advection was not considered. This assumption is equivalent to considering in the numerical models; a vertical fault with a slip equals the vertical component of the analog model fault slip.

## 2.3. Common Modeling Assumptions

Our approach relies on the comparison between the analog and numerical modeling results. The parameter values used in the numerical models come from studies of natural morphological landscapes (see Section 2.1.4). In this paper, the spatial and temporal scales used in the numerical models are therefore identical to those associated with natural processes. Hence, the scaling coefficients presented above define the geometry, the initial and boundary conditions, and the time scale of the numerical models.

### 2.3.1. Geometry

For the experimental approach, we considered a square model with a side length $L$ of 60 cm. Its initial topography consisted of an inclined plane of $0.8^\circ$ leading to a change in elevation of 8.5 mm when $y$ varies from 0 to $L$ (Figures 1 and 2). A reverse fault crossed the model at $y = 0$ (Figure 2).

The numerical model experiments were performed on a regular square lattice with a similar geometry. Using the scaling factor given in Table 1, we considered a $6000 \times 6000$ m model with a node spacing of 50 m. This discretization led to a mesh size of $120 \times 120$, which was a good compromise between computation time and model resolution of the calculated surface.

### 2.3.2. Initial and Boundary Conditions

Initially, the models did not undergo any external or internal forcing. Top ($y = L/6$), left ($x = -L/2$), and right ($x = L/2$) boundaries were closed, i.e., no water or sediment output flux was permitted. An outflow of water and sediment was only allowed through the bottom face ($y = -5L/6$). We simulated the formation and evolution of an alluvial fan by adding a flow of water and sediment coming from the top-center of the model (Figure 2).

In the analog approach, a constant water flux of $1.4 \times 10^{-4}$ L.s$^{-1}$ with a sediment concentration of 10% was brought from the reservoir described in Section 2.1.1. The water flux was upscaled to $140 \times 10^3$ m$^3$.yr$^{-1}$ in the numerical approach and imposed
from the Landlab variable runoff\_rate, assigned to the grid field water\_unit\_flux\_in. The functionality to add an external input sediment flux is not implemented in the current Landlab version. The ErosionDeposition component of Landlab was then modified to provide an external sediment flux of 14 × 10³ m³.yr⁻¹ as a source term to the node with coordinates (x = 0, y = L/6).

Concerning tectonic forcing, the uplift pattern was controlled by a slip rate *v* along a *α* = 43° dipping thrust fault (Figure 2b) with no flexural response associated with surface loading. In the reference analog model, *v*<sub>ref</sub> equaled 1 cm.h<sup>-1</sup>, leading to *v*<sub>ref</sub> = 2.8 mm.yr<sup>-1</sup> in numerical models (Table 1). Two fault slip mode end-members were tested (Figure 2c). First, we assumed a stable fault behavior with a constant uplift rate *v*<sub>z</sub> = *v* × sin(*α*) imposed at all time steps. Next, we assumed a fault stick-slip behavior for which the cumulated slip over a return period *T* was released during a single time step. The associated coseismic vertical displacement *u*<sub>z</sub> was then defined as the product *v*<sub>z</sub> × *T*. In the reference analog model, we assumed a return period *T*<sub>ref</sub> = 180 s, leading to *T*<sub>ref</sub> = 1800 yr for the numerical reference model.

2.3.3. Time Scenario

Laboratory experiments evolved through two main stages. The first stage lasted 1.5 h and was dedicated to fan building without tectonics (Figure 2c). Sediments then accumulated on the table so that a sufficient thickness (nearly 10 mm) avoided a too early complete incision of the fan when, in a second phase, the fault was set active during the so-called “tectonic forcing” stage, which lasted 1 h (Figure 2c).

To be consistent with the analog modeling scaling procedure, the total duration of numerical simulations was set to 90 ka with a time step of 10 yr. These simulations were also divided into two stages: (1) a sedimentation stage of 54 ka without tectonics, where the transport and deposition of sediments contributed to forming an alluvial fan, (2) a stage of 36 ka where all surface and tectonic processes acted together.

**3. Results**

In this section, we first tested the analog experiments’ repeatability. Next, we calibrated
the numerical erosion law parameters. Finally, we explored the effect of tectonics and
surface processes on alluvial fan morphological evolution.

3.1. *Repeatability of Analog Experiments*

Using the same operators, equipment, and boundary conditions, eight experiments
were performed to test the robustness of the analog results. This repeatability test was
performed using the conditions presented in Section 2.3. for the first modeling stage, i.e.,
without tectonic forcing.

Regardless of the experiment, an alluvial fan formed within the first ten minutes and
then grew due to the avulsion phenomenon (e.g., [42]). After one hour, the resulting fans
had an identical form with a quasi-radial pattern elongated in the y-direction (Figure 3).

Taking advantage of the specific color of deposits, we detected the fan boundaries
from the top-view pictures taken after one hour (Figure 3). We then assessed the fan length
in the x-direction, Lx, and in the y-direction, Ly. We also used the fan DEM obtained from
Micmac processing to calculate the elevation of the fan apex ha.

Based on the experiment dataset, we estimated the mean value and the standard
deviation of these three geometric parameters. Our results demonstrate a good repeatability
with L<sub>x</sub> = 30.7 ± 3.7 cm, L<sub>y</sub> = 36.7 ± 2.2 cm, h<sub>a</sub> = 18.5 ± 2 mm. These values were used
in the following section to calibrate the erosion law parameters of the numerical models.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Experimental alluvial fan geometry obtained after the 1-hour first modeling phase without tectonics. (a–h) Top-view pictures of alluvial fans for each experiment. The alluvial fan boundary is orange. $L_x$ is for width, $L_y$ is for length, and $h_a$ for apex height. (i) Alluvial fans geometry of all analog experiments. The coordinate system is centered in the middle of the fault. Orange lines give the fan boundaries' locations. The black box represents the deposition area.</figcaption></figure>

## 3.2. Calibration of Numerical Erosion Law Parameters from the Alluvial Fan Building Phase

The erosion law parameters depend on climate, lithology, or vegetation cover, and some of them are interdependent, such as $\kappa_{hill}$, $m$, and $n$. Therefore, assessing these parameters from natural landscapes results in a wide variance, limiting the use of erosion laws and thus the relevance of numerical models.

We calibrated these parameters by performing a numerical inversion, which compared the geometry of alluvial fans obtained from analog and numerical modeling. We imposed consistent geometry, initial, and boundary conditions for both modeling approaches, verifying the estimated scaling factors. This comparison was performed based on a systematic exploration of the erosion law parameters $\kappa_{hill}$, $v_s$, $n$, and $\kappa_{river}$ in the ranges of values
presented in Section 2.2, resulting in a collection of nearly 10,000 numerical models. We compared the calculated normalized $L_x/L$, $L_y/L$, and $h_a/L$ with those obtained from analog models for each of them. $L_x$ and $L_y$ were obtained (1) by calculating the difference between the final topography and the initial elevation associated with a 0.8° inclined plane and (2) by assuming that the fan boundary is defined by a threshold elevation of 2 mm. One can note that $L = 60$ cm for analog models and $L = 6$ km for numerical models.

We defined the cost function as the weighted root mean squared error,

$$cost = \sqrt{\frac{1}{3} \sum_{i=1}^{3} \left[ \left( \frac{obs_i - calc_i}{sigma_{obs_i}} \right)^2 \right]} \quad (7)$$

The two vectors *obs* and *sigmaobs* included $L_x/L$, $L_y/L$, $h_a/L$, and their standard deviation from analog modeling. The vector *calc* is associated with the same parameters calculated from numerical modeling.

The best model inferred from our inversion was associated with the minimum cost function. It was obtained for $\kappa_{river} = 10^{-5}$, $v_s = 1$ m.yr$^{-1}$, $\kappa_{hill} = 8.10^{-4}$ m.yr$^{-1}$, and $n = 2.5$ (Figure 4). The quasi-radial pattern elongated in the y-direction observed in the analog model can thus be explained by a low river bedrock erodibility and an equal contribution of incision and sediment transport. The horizontal geometry of the best-fitting simulated fan was characterized by normalized lengths $L_x/L$ of 0.57 and $L_y/L$ of 0.57, consistent with $L_x/L = 0.51 \pm 0.06$, $L_y/L = 0.61 \pm 0.04$ obtained from analog models. Finally, the numerical and analog approaches gave a very similar topography apex with $h_a/L = 0.32$ and $h_a/L = 0.31 \pm 0.03$, respectively. Next, we defined the likelihood function as

$$likelihood = exp(-cost) \quad (8)$$

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Alluvial fan geometry obtained after one hour of fan building without tectonics. The color scale gives the calculated elevation of the best-fitting numerical model with $\kappa_{river} = 10^{-5}$, $v_s = 1$ m.yr$^{-1}$, $\kappa_{hill} = 8.10^{-4}$ m.yr$^{-1}$, and $n = 2.5$. Orange lines give the location of the alluvial fans' boundaries obtained in the analog experiments (see Figure 3). The thick black line is the fan boundary obtained from the numerical approach. Distance and elevation are normalized by the model size $L$, which is equal to 60 cm and 6 km in analog and numerical models, respectively. In the analog experiments, the fan apex height $h_a$ is 1.85 cm, which corresponds to an elevation/L ratio of 0.031.</figcaption></figure>
The likelihood distribution was normalized with the likelihood of the best-fitting model presented above. Hence, the calculated normalized likelihood ranged between 0 (low agreement between analog and numerical results) and 1 (best agreement).

According to Figure 5, we showed that the erosion parameters did not control the likelihood distribution in the same way. Not surprisingly, the likelihood distribution did not depend on the assumed $\kappa_{hill}$ values. Indeed, alluvial fan building was associated with slopes <5°, which prevent hillslope processes. Our calibration approach was thus not relevant to constrain $\kappa_{hill}$.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Normalized likelihood distribution of numerical modeling parameters calculated from alluvial fan geometries obtained in eight experimental models (see Figure 3). $\kappa_{river}$ is the river bedrock erodibility. Its unit depends on the water discharge exponent $m$, which is defined as $m = 0.5 n$ in this study. $n$ is the slope exponent. $v_s$ is the net settling velocity. $\kappa_{hill}$ is the erodibility coefficient associated with the hillslope process. A red cross is associated with the best-fitting model with $\kappa_{river} = 10^{-5}$, $v_s = 1 \text{ m.yr}^{-1}$, $\kappa_{hill} = 8.10^{-4} \text{ m.yr}^{-1}$ and $n = 2.5$.</figcaption></figure>

In contrast, fluvial erosion and deposition significantly influenced fan formation by controlling channel avulsion and migration to a new path radiating from the apex fan. Our result suggests a net settling velocity $v_s$ of ca. 1 m.yr⁻¹ (Figure 5). This value was associated with an intermediate case between detachment-limited and transport-limited models [13]. The more likely models also suggested a slope exponent $n$ between 0.5 and 3 and a river bedrock erodibility between $10^{-5}$ and $10^{-2}$. In more detail, the obtained likelihood distribution showed a semi-logarithmic relationship between $n$ and $\kappa_{river}$: a low $\kappa_{river}$ is associated with a high $n$, and vice versa (Figure 5). This relationship can
be explained by the power incision law itself (Equation (6)), for which a decrease in the exponent *n* can compensate for a logarithmic increase in $κ_{hill}$.

Together, these first results demonstrate the robustness of both approaches and allow to put forward their complementarity to study the physical processes controlling landscape evolution. In the following section, we will use this joint approach to study the morphotectonic evolution of an alluvial fan.

### 3.3. Morpho-Tectonic Evolution

In the previous section, we modeled an alluvial fan formation only due to water and sediment flux external input. Here, we studied the morphological evolution of a fan affected by additional tectonic forcing, all other conditions remaining unchanged. We assumed a stable fault with a constant uplift rate of 0.7 cm.h⁻¹ and 1.9 mm.yr⁻¹ for analog and numerical models, respectively. This new boundary condition was applied for 1 h for experimental models, leading to an equivalent duration of 36 ka for the numerical models.

The morphology obtained with the analog models exhibited a more complex pattern compared to the surface without tectonic forcing (Figure 6). Two major landforms can be observed from the hanging wall topography: (1) an inherited alluvial fan associated with the fan building stage and (2) a central valley. This valley was narrow at the point source of sediments and widened until it crossed the fault (Figure 6b). On both sides of this central channel, one can mention very narrow abandoned channels suggesting a temporal evolution of the channel network. These abandoned channels represent the direction variations of the sedimentary flux at the outlet. The footwall was marked by a neo-formed fan, highlighted by its semicircular shape and visible drains (Figure 6a). The location of this new fan at the outlet of the central valley underlines the complex tectonic-erosion interactions between the two fault blocks. While the sediments deposited on the footwall were carried in the hanging wall central valley, the footwall fan elevation controlled the incision rate in this valley and its sediment transport capacity.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Final stage of an analog experiment. (a) General view. The reverse fault is outlined with a red line and triangles, and the contour of the alluvial fan appears in orange. The central valley is limited by red dashed lines and the flux entry with an orange arrow. (b) Close-up picture of the central valley. The minimum width outlined by a white arrow is the area where the flux is active. The maximum width includes old, abandoned river channels.</figcaption></figure>
The numerical simulations were performed using the 200 best-fitting models obtained in the previous calibration section, i.e., those having a normalized likelihood exceeding ca. 0.5. The numerical models most in line with the results of the analog approach suggested a river bedrock erodibility of 10⁻⁴, a settling velocity of 1 m.yr⁻¹, and a slope exponent of 2. Hillslope erodibility remained the least constrained parameter. A range between 5 × 10⁻⁷ m.yr⁻¹ and 2 × 10⁻³ m.yr¹ gave almost identical results. These best-fitting models exhibited a normalized elevation whose deviation from the surface of analog models did not exceed 0.002 (Figure 7). For both approaches, the normalized topographies were thus quite similar with (1) an elevation of up to 0.04, (2) a lateral extension of the inherited fan on the hanging wall of 0.8, and (3) a neo-formed fan on the footwall. The main difference was the absence of a wide central valley. Indeed, in numerical models, the width of this valley was limited to the mesh size used in the numerical models. However, we can see other narrow abandoned valleys on the hanging wall, suggesting a non-static river network with lateral movements of central channels during the simulation (Figure 7b).

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Final topography obtained from both experimental and numerical approaches. (a) Experimental modeling result. Color scale gives the normalized elevation associated with a fan building stage of 1.5 h followed by a 1h fault activity at a slip rate of 1 cm/h. Dashed lines show the location of profiles plotted in Figure 6c,d. (b) Same as (a) for numerical modeling with $\kappa_{river} = 10^{-4}$, $v_s = 1 \text{ m.yr}^{-1}$, $\kappa_{hill} \in [5 \times 10^{-7}; 2 \times 10^{-3}] \text{ m.yr}^{-1}$, and $n = 2$, assuming a fan building stage of 54 ka and a fault slip rate of 2.8 mm.yr⁻¹t during 36 ka. (c,d) show swath profiles obtained from experimental (color circles) and numerical (black lines) approaches.</figcaption></figure>
To test this hypothesis, we analyzed the temporal evolution of the modeled landscapes in more detail. We identified four distinct stages from the numerical modeling results (Figure 8). During the fan building stage, one can observe a radial network of channels diverging from the apex fan and spreading laterally over the entire model's width (Figure 8a). At the onset of tectonic forcing, the surface deformation associated with cumulated fault slip is small compared to the geometry inherited from the fan formation. The geometry of the drainage network, therefore, remains unchanged (Figure 8b). Then, during a transition period, the water and sediment flow becomes channeled. A wide central valley forms, where the drains can move laterally (Figure 8c). Finally, when the cumulative displacement on the fault becomes large enough, the central valley deepens, and the channel remains localized along a single drain in the hanging wall (Figure 8d). We observed a similar scenario for analog modeling. A central channel forms by incising the hanging wall alluvial fan after nearly 30 min of tectonic forcing. The drainage network becomes more channeled, forming a valley of ca. 5 cm width at the end of the experiment.

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Landscape evolution of numerical modeling over a time duration of *t*<sub>end</sub>. The color scale shows the calculated elevation. Blue lines are associated with the channel drainage over a normalized time window of 1/27. Red circles give the location of the river outlet along the fault. (a) *Time* / *t*<sub>end</sub> = 0.17—Fan building period. (b) *Time* / *t*<sub>end</sub> = 0.42—Beginning of the tectonic forcing. (c) *Time* / *t*<sub>end</sub> = 0.63—The tectonic forcing causes the channel drainage localization on the hanging wall within a central valley. A new alluvial fan forms on the foot wall. (d) *Time* / *t*<sub>end</sub> = 0.88—Final phase of tectonic loading. The central valley deepens, and the channel remains localized along a single drain in the hanging wall.</figcaption></figure>

To further compare the two modeling approaches, we estimated the temporal evo-
lution of the zone where the channels and the fault intersect (Figures 6b and 8 and
Supplementary Materials S3). Thanks to the drainage network, the width of this zone
was defined as the largest distance between the channels. While during the fan building
period, a width *L* was obtained, at the end of the experiment, the width corresponded to
the lateral extension of the central valley. Our results demonstrate that the mean channel
slope was the key parameter in this temporal evolution. During the fan building period, the channel slope increased gently in response to sediment deposition (Figure 9). This low slope allowed the drains to move laterally along the entire model. The application of tectonic forcing generated a rapid slope increase and a continuous uplift rate. As long as the critical slope of ca. 0.081 was not reached, there was no channelization. Beyond this threshold, one can observe the channelization of the central valley with a pseudo-stable state of the mean channel slope. This slope reached a stable value of 0.85 before the end of the simulations, highlighting an equilibrium state.

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Time evolution of morphological metrics over a time period of $t_{end}$, which includes the last stage of the fan building period over $t_{end}/3$ and tectonic forcing over the last $2t_{end}/3$. Round wheat boxes give the number of figures associated with the four stages described Figure 8. (a) Time variation of the channel slope along the main central valley. (b) Mean accumulated vertical displacements along the fault scarp, which includes both tectonic uplift and fan deposition. (c) Red circles give the location of channel outlets along the fault scarp obtained from numerical modeling. Dashed blue lines are the envelope defined by the local maxima and minima observed in the experimental approach (see Figure 6).</figcaption></figure>

Our results confirm the relevance of our joint approach to explaining the obtained topography and the temporal evolution of the drainage network associated with constant tectonic forcing.
### 3.4. Stable versus Stick-Slip Behaviors

Over the last decade, many studies pointed out the diversity of slip modes along active faults, including stable slip [24], slow slip events [43], and earthquakes [44]. Distinguishing these different slip modes is of primary importance to improve the seismic hazard assessment. Here, we compared the previous results assuming a stable sliding with the topography obtained for a stick-slip fault. The cumulated fault slip remained the same, and we considered a return period $T_{ref}$ of 180 s and 1800 yr for the analog and numerical modeling, respectively (Table 1).

Our results suggest a specific morphological signature associated with the assumed fault slip mode (Figure 10). However, this signature is tenuous and depends on the modeling approach. To analyze these results, we defined the normalized differences in elevation as

$$ \Delta z = \frac{Elevation_{stick-slip} - Elevation_{Stable-slip}}{L} \qquad (9) $$

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Stick-slip versus stable fault slip mode. Positive differences in elevation indicate larger vertical displacements for stick-slip compared to stable slip mode. (a) Result from experimental modeling associated with 20 seismic events with a return period of 180 s and a vertical displacement of 0.34 mm. (b) Result from numerical modeling for 20 seismic events with a return period of 1800 yr and a vertical displacement of 3.4 m. (c) Repeatability of analog models. Blue and orange lines give the cumulative and reversed cumulative normalized error, respectively. At each point, this error is defined as the difference in elevation between analog models with the same properties and undergoing similar tectonic forcing. A perfectly reproducible model will have a zero error.</figcaption></figure>

The absolute value of $\Delta z$ reached up to $10^{-3}$ in analog models. This value is larger than that obtained in the numerical simulation. While $\Delta z$ was positive in analog models, $\Delta z$ was negative in numerical simulations in the hanging wall. Analog models showed a more complex pattern in the footwall with $\Delta z > 0$ at the center of the neoformed fan and $\Delta z < 0$ at its boundaries (Figure 10a). One can also notice that the highest values of numerical $\Delta z$ were located near the fault scarp (Figure 10b).

One explanation of these differences in $\Delta z$ between both approaches comes from the repeatability of analog models (Figure 10c). Based on the comparison of four laboratory experiments, our repeatability test highlights that nearly 40% of the analog model has a normalized error of more than $5 \times 10^{-4}$, which is of the same magnitude as the elevation difference $\Delta z$. This error shows the limitations of analog modeling, in which both the water and sedimentary input flux needs to be better controlled. The analog results are therefore not robust enough to be conclusive. In the following section, we discuss in more detail the morphological signature of the fault slip mode using only the numerical approach.

To sum up, in this results section, we showed the interest of a joint approach combining analog and numerical models to study the morpho-tectonic evolution of an alluvial fan. The consistency of results obtained during the fan building and tectonic stages makes our
results more robust. Analog models are advantageous to calibrate erosion law parameters used in numerical simulations. The numerical approach can overcome some of the analog model limitations, such as the control of sediment input. Our results suggest that the tested analog models correspond to an intermediate case between detachment-limited and transport-limited models, with a low river bed-rock erodibility compensated by a high slope exponent. The mean channel slope appears as a critical parameter controlling the temporal evolution of our modeling. The obtained fault slip mode signature is tenuous and seems challenging to observe in the field.

# 4. Discussion

The tectonic forcing depends on several parameters, such as the slip rate *v*, the slip return period *T*, the fault dip angle *α*, and the fault location *y*<sub>fault</sub>. As previously mentioned, the elevation change due to horizontal tectonic advection is not considered in the numerical simulations. Changing *α* is thus equivalent to changing *v*. Here, we only studied the role played by *v*, *T*, and *y*<sub>fault</sub> on the morphotectonic evolution of the alluvial fan. Therefore, we changed these tectonic parameters one by one and compared the results to the reference model used until now (Table 1).

## 4.1. Effect of Slip Rate on Surface Morphology

For the reference model, we used a homogeneous slip rate $v_{ref} = 2.8 \text{ mm.yr}^{-1}$ along the fault. This value is relatively low compared to field observations in active fault zones (e.g., Marechal et al., 2016). Three slip rates were tested: one is the reference model $v = v_{ref}$ (Table 1), and the other two are associated with a faster slip rate $v = 1.5 \times v_{ref}$ and $v = 3 \times v_{ref}$.

Due to tectonic-erosion feedback, whatever the slip mode, the channel incision was more intense when the slip rate increased (Figure 11). In contrast with the reference model, the obtained topography exhibited a well-marked central valley, even for a slight slip rate increase. This topography is in better agreement with the analog model, which showed a valley in the hanging wall widening from its source to its outlet at the fault zone (Figure 6). This result underlines the too low fault slip rate used in the numerical reference model. It should be remembered that only the order of magnitude of the scaling factor is assessed in Section 2.1.3. The coefficient of 1.5 required to explain a large central valley does not question the previously obtained results.

Regardless of the slip rate, finding a difference between a stable or unstable slip pattern is not straightforward. The elevation deviation was calculated from Equation (8). For the three models, we obtained a rather similar pattern for $\Delta z$. The main difference was in the two zones of the fault scarps with high $|\Delta z|$, suggesting a more localized deformation with stable sliding (Figure 11). These two zones were parallel to the fault and had a width and an amplitude dependent on the fault slip rate. Ten additional models with slip rates ranging from $1.2 \times v_{ref}$ to $3.6 \times v_{ref}$ were tested to confirm this finding. Our result suggests linear relationships between the slip rate, the zone width, and the deviation amplitude: the higher the fault slip velocity, the higher the width, and the higher the amplitude (Figure 12).

## 4.2. Return Period

The reference return period $T_{ref}$ is equal to 1800 years (Table 1). Except in particular regions with a very low shortening rate (e.g., [43]), the return period of strong earthquakes ($M > 7$) is generally a few hundred years (e.g., [23,45]). Here, we ran seven models with return periods between $0.05 \times T_{ref}$ and $1.1 \times T_{ref}$. A limitation of numerical modeling is that long return periods are associated with a large displacement on nodes in a single time step, leading to numerical instabilities. All these models had the same cumulative slip and were associated with a slip rate of $3 \times v_{ref}$.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Morphological signature of fault slip rate. (a) Elevation associated with stick-slip behavior assuming the fault slip rate $v_{ref}$ used in the numerical reference model (see Table 1). (b,c) same as (a), assuming a fault slip rate of $1.5 \times v_{ref}$ and $3 \times v_{ref}$, respectively. (d-f) same as (a-c), assuming a stable sliding. (g-i) elevation deviation between slick-slip and stable sliding for a fault slip rate of $v_{ref}$, $1.5 \times v_{ref}$, and $3 \times v_{ref}$, respectively. W is the width of the zone with a high elevation deviation.</figcaption></figure>

The return period did not influence the main topographic features of the models. Following the approach described in the previous paragraph, the width of the fault scarp's zones of high $|\Delta z|$ remained constant with a variation of the return period (Figure 13a). Unsurprisingly, the return period influenced the amplitude of the area with higher differences in $|\Delta z|$. The stick-slip fault behavior tended towards a stable sliding when the return period became shorter. Our results suggest a linear relationship between the return period and the deviation amplitude: the longer the return period, the higher the amplitude (Figure 13b).

### 4.3. Distance between Fault and Sediment Source

The accommodation space of deformation along major thrust faults is often complex at the surface. Strain can be distributed along several branches in a fault zone extending for several kilometers [27,46–50]. In the reference model, the distance $d = d_{ref}$ between the sediment source and the fault is equal to 1000 m (Figure 2). Three other distances were tested: $d_{ref}/2$, $d_{ref} \times 2$, and $d_{ref} \times 3$. For all numerical models, we assumed a fault slip rate of $3 \times v_{ref}$.

Our simulations underline the significant impact of $d$ on morphology. For large distances, the radial geometry of the drainage network promoted incision of the scarp along its entire length (Figure 14a). This geometry led to a homogeneous drainage density. No secondary alluvial fan was then formed on the foothill. As this distance decreased, the flow of water and sediment became more channeled, forming (1) a central valley that drains most of the flow and (2) a secondary alluvial fan on the foothill (Figure 14b,c). Finally, all water flow and sediment transport occurred in a central valley for the shortest distances (Figure 14c,d). Note that the distance did not control the width of the central valley.
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Morphological signature of slip rate along the fault scarp outside the area affected by the alluvial fan (|X/L| > 0.2). (a) Mean elevation profiles associated with stick-slip (blue line) and stable sliding (orange line). The green line shows the deviation between these two fault slip behaviors. Note the coefficient of ten used to highlight this deviation. W and A are the deviation width and amplitude, respectively. (b) Red circles show the variation of W obtained from numerical modeling with respect to fault slip rate. The obtained linear trend is given by the red line. $v_{ref}$ is the fault slip rate used in the reference model (see Table 1). (c) Relationship between the deviation amplitude A and the fault slip rate obtained from numerical modeling (green circles). The green line shows the linear trend between these two parameters.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Role of the earthquake returns period in the deviation between slick slip and stable sliding morphology. A slip rate of $3 \times v_{ref}$ and modeling duration time of $t_{end}$ including a fan building period and tectonic forcing were assumed (see Table 1). (a) Red circles show the variation of W obtained from numerical modeling with respect to the earthquake return period. The obtained linear trend is given by the red line. $T_{ref}$ is the return period used in the reference model (see Table 1). (b) Relationship between the deviation amplitude A and the fault slip rate obtained from numerical modeling (green circles). The green line shows the linear trend between these two parameters.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Role of fault location on alluvial fan morphology using a fault slip rate of 3 × *v*<sub>ref</sub>. Color scale gives the normalized elevation. The red dashed line indicates the fault trace location. *d*<sub>ref</sub> = L/6 is the distance in the reference model between the top edge of the model and the fault trace. The effect of this distance is tested for (a) *d*<sub>ref</sub> × 3, (b) *d*<sub>ref</sub> × 2, (c) *d*<sub>ref</sub> and (d) *d*<sub>ref</sub> / 2.</figcaption></figure>

Together, our results demonstrate that *d* is a key parameter in the morpho-tectonic evolution of an alluvial fan. The closer the fault is to the sediment source, the higher the drainage density in the center of the hanging wall and the higher the surface of the neo-formed alluvial fan.

**5. Conclusions**

This paper investigated the role of surface processes and tectonics on the evolution
of the alluvial fan morphology. We used a joint approach based on analog and numeri-
cal models sharing identical geometry and boundary conditions. We demonstrated the
excellent repeatability of the analog models by showing that the main geometrical charac-
teristics of the modeled fans can be reproduced with a deviation of ca. 10%. Based on a
numerical inversion, we use these experimental results to calibrate erosion law parameters.
While our approach is not relevant to assess hillslope erodibility, we found a relationship
between bedrock river erodibility and the slope exponent of the incision power law. To
our knowledge, this type of combined approach is rarely developed. Our results illus-
trate its efficiency in better quantifying surface processes that would otherwise remain
poorly constrained.
The analog and the numerical results indicated a good consistency on both the final topography and the temporal evolution of the alluvial fan morphology. We obtained similar results concerning the genesis of specific landforms such as (1) an alluvial fan formed during the first modeling stage with no tectonics, (2) a wide central valley in the hanging wall, and (3) a neo-formed alluvial fan in the footwall. A detailed analysis of the model evolution reveals a threshold process controlling the temporal changes in the drainage network geometry. As long as the critical slope of ca. 0.081 was not reached, there was no channelization. Beyond this threshold, the channelization of the central valley took place.

The morphologic signature of fault slip mode is tenuous. We only spotted a minor deviation along the fault scarp. Although this deviation may increase with the fault slip rate and the return time of slip events, it remained low compared to fault scarp cumulated displacement. In nature, such a deviation is challenging to identify. High-resolution topographic data cannot be used alone but appears valuable to complement stratigraphic and geodetic approaches to study the slip mode of active faults.

This study assumed a constant external forcing associated with no time variation in the input water and sediment flux. However, such a constant regime may be seldomly found in nature, and the relative impact of extreme events in shaping landscape remains an open question [51,52]. The role of climate variations must be better considered in future works, including extreme processes such as heavy rainfalls or flash floods, and variations at regional base levels.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10.3390/geosciences11100412/s1, Supplement 1- Flow transfer monitoring, Figure S1: time evolution of the mass of sediment during the analog experiments, Supplement 2- Material Physical Properties, Table S1: Material and components characteristics, Supplement 3- Movie of the analog experiments, Movie S3: Creep-2-b and Movie S4: Earthquake-2-b.

**Author Contributions:** Funding acquisition, R.C.; Investigation, C.G.-E.; Methodology, C.G.-E., Y.C., R.C., S.D., and R.S.; Project administration, R.C.; Software, R.C.; Supervision, R.C.; Writing—original draft, C.G.-E., and R.C.; Writing—review & editing, C.G.-E., Y.C., R.C., S.D., and R.S.. All authors have read and agreed to the published version of the manuscript.

**Funding:** Clément Garcia-Estève's Ph.D. is supported by a fellowship from the French Ministry for Higher Education and Research. This study is supported by the French Agence National de la Recherche (project ANR-18-CE01-0017 (Topo-Extreme)).

**Acknowledgments:** We thank Christian Romano for the technical support.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Bishop, P. Long-term landscape evolution: Linking tectonics and surface processes. *Earth Surf. Process. Landf. J. Br. Geomorphol. Res. Group* **2007**, *32*, 329–365. [CrossRef]

2. Burbank, D.W.; Pinter, N. Landscape evolution: The interactions of tectonics and surface processes. *Basin Res.* **1999**, *11*, 1–6. [CrossRef]

3. Naylor, L.A.; Spencer, T.; Lane, S.N.; Darby, S.E.; Magilligan, F.J.; Macklin, M.G.; Möller, I. Stormy geomorphology: Geomorphic contributions in an age of climate extremes. *Earth Surf. Process. Landf.* **2017**, *42*, 166–190. [CrossRef]

4. Turowski, J.M.; Cook, K.L. Field techniques for measuring bedrock erosion and denudation. *Earth Surf. Process. Landf.* **2017**, *42*, 109–127. [CrossRef]

5. Giano, S.I.; Schiattarella, M. Age constraints and denudation rate of a multistage fault line scarp: An example from Southern Italy. *Geochronometria* **2014**, *41*, 245–255. [CrossRef]

6. Giano, S.I.; Pescatore, E.; Agosta, F.; Prosser, G. Geomorphic evidence of Quaternary tectonics within an underlap fault zone of southern Apennines, Italy. *Geomorphology* **2018**, *303*, 172–190. [CrossRef]

7. Migon, P. Book Review: Morphotectonics. *Prog. Phys. Geogr. Earth Environ.* **2006**, *30*, 138–140. [CrossRef]

8. Weldon, R.; Scharer, K.; Fumal, T.; and Biasi, G. Wrightwood and the earthquake cycle: What a long recurrence record tells us about how faults work. *GSA Today* **2004**, *14*, 4. [CrossRef]

9. Harel, M.A.; Mudd, S.M.; Attal, M. Global analysis of the stream power law parameters based on worldwide ¹⁰Be denudation rates. *Geomorphology* **2016**, *268*, 184–196. [CrossRef]
10. Richardson, P.W.; Perron, J.T.; Schurr, N.D. Influences of climate and life on hillslope sediment transport. *Geology* **2019**, *47*, 423–426.
[CrossRef]

11. Barnhart, K.R.; Hutton, E.W.; Tucker, G.E.; Gasparini, N.M.; Istanbulluoglu, E.; Hobley, D.E.; Lyons, N.J.; Mouchene, M.; Nudurupati, S.S.; Adams, J.M.; et al. Landlab v2. 0: A software package for Earth surface dynamics. *Earth Surf. Dyn.* **2020**, *8*, 379–397. [CrossRef]

12. Graveleau, F.; Dominguez, S. Analogue modelling of the interaction between tectonics, erosion and sedimentation in foreland thrust belts. *Comptes Rendus Geosci.* **2008**, *340*, 324–333. [CrossRef]

13. Graveleau, F.; Hurtrez, J.E.; Dominguez, S.; Malavieille, J. A new experimental material for modeling relief dynamics and interactions between tectonics and surface processes. *Tectonophysics* **2011**, *513*, 68–87. [CrossRef]

14. Graveleau, F.; Strak, V.; Dominguez, S.; Malavieille, J.; Chatton, M.; Manighetti, I.; Petit, C. Experimental modelling of tectonics–erosion–sedimentation interactions in compressional, extensional, and strike–slip settings. *Geomorphology* **2015**, *244*, 146–168.
[CrossRef]

15. Hobley, D.E.; Adams, J.M.; Nudurupati, S.S.; Hutton, E.W.; Gasparini, N.M.; Istanbulluoglu, E.; Tucker, G.E. Creative computing with Landlab: An open-source toolkit for building, coupling, and exploring two-dimensional numerical models of Earth-surface dynamics. *Earth Surf. Dyn.* **2017**, *5*, 21–46. [CrossRef]

16. Paola, C.; Straub, K.; Mohrig, D.; Reinhardt, L. The ‘unreasonable effectiveness’ of stratigraphic and geomorphic experiments. *Earth-Sci. Rev.* **2009**, *97*, 1–43. [CrossRef]

17. Strak, V.; Dominguez, S.; Petit, C.; Meyer, B.; Loget, N. Interaction between normal fault slip and erosion on relief evolution: Insights from experimental modelling. *Tectonophysics* **2011**, *513*, 1–19. [CrossRef]

18. Buiter, S.J.; Babeyko, A.Y.; Ellis, S.; Gerya, T.V.; Kaus, B.J.; Kellner, A.; Schreurs, G.; Yamada, Y. The numerical sandbox: Comparison of model results for a shortening and an extension experiment. *Geol. Soc. Lond. Spec. Publ.* **2006**, *253*, 29–64.
[CrossRef]

19. Reitano, R.; Faccenna, C.; Funiciello, F.; Corbi, F.; Willett, S.D. Erosional response of granular material in landscape models. *Earth Surf. Dyn.* **2020**, *8*, 973–993. [CrossRef]

20. Mao, Y.; Li, Y.; Yan, B.; Wang, X.; Jia, D.; Chen, Y. Response of Surface Erosion to Crustal Shortening and its Influence on Tectonic Evolution in Fold-and-Thrust Belts: Implications from Sandbox Modeling on Tectonic Geomorphology. *Tectonics* **2021**, *40*, e2020TC006515. [CrossRef]

21. Dal Zilio, L.; Jolivet, R.; van Dinther, Y. Segmentation of the Main Himalayan Thrust illuminated by Bayesian inference of interseismic coupling. *Geophys. Res. Lett.* **2020**, *47*, e2019GL086424. [CrossRef]

22. Langridge, R.M.; Villamor, P.; Almond, P.; Basili, R.; Hemphill-Haley, M.; Ries, W. Late Holocene paleoseismicity of the Australia-Pacific plate boundary in central South Island: The Alpine to Hope Fault transition. *Geol. Soc. Am. Abstr. Programs* **2009**, *41*, 691.

23. Le Roux-Mallouf, R.; Ferry, M.; Cattin, R.; Ritz, J.F.; Drukpa, D.; & Pelgay, P. A 2600-year-long paleoseismic record for the Himalayan Main Frontal Thrust (western Bhutan). *Solid Earth* **2020**, *11*, 2359–2375. [CrossRef]

24. Marechal, A.; Mazzotti, S.; Cattin, R.; Cazes, G.; Vernant, P.; Drukpa, D.; Thinley, K.; Tarayoun, A.; Le Roux-Mallouf, R.; Thapa, B.B.; et al. Evidence of interseismic coupling variations along the Bhutan Himalayan arc from new GPS data. *Geophys. Res. Lett.* **2016**, *43*, 12–399. [CrossRef]

25. Hubbert, M.K. Theory of scale models as applied to the study of geologic structures. *Bull. Geol. Soc. Am.* **1937**, *48*, 1459–1520.
[CrossRef]

26. Cobbold, P.R.; Jackson, M.P.A. Gum rosin (colophony): A suitable material for thermomechanical modelling of the lithosphere. *Tectonophysics* **1992**, *210*, 255–271. [CrossRef]

27. Davy, P.; Cobbold, P.R. Experiments on shortening of a 4-layer model of the continental lithosphere. *Tectonophysics* **1991**, *188*, 1–25.
[CrossRef]

28. Horsfield, K. Morphology of branching trees related to entropy. *Respir. Physiol.* **1977**, *29*, 179–184. [CrossRef]

29. Malverti, L.; Lajeunesse, E.; Métivier, F. Small is beautiful: Upscaling from microscale laminar to natural turbulent rivers. *J. Geophys. Res. Earth Surf.* **2008**, *113*. [CrossRef]

30. Niemann, J.D.; Hasbargen, L.E. A comparison of experimental and natural drainage basin morphology across a range of scales. *J. Geophys. Res. Earth Surf.* **2005**, *110*. [CrossRef]

31. Peakall, J.; Ashworth, P.; Best, J. Physical modelling in fluvial geomorphology: Principles, applications and unresolved issues. *Sci. Nat. Geomorphol.* **1996**, *221*, 253–263.

32. Liu, Y.; Métivier, F.; Lajeunesse, É.; Lancien, P.; Narteau, C.; Ye, B.; & Meunier, P. Measuring bedload in gravel-bed mountain rivers: Averaging methods and sampling strategies. *Geodin. Acta* **2008**, *21*, 81–92. [CrossRef]

33. Jolivet, M.; Barrier, L.; Dominguez, S.; Guerit, L.; Heilbronn, G.; Fu, B. Unbalanced sediment budgets in the catchment–alluvial fan system of the Kuitun River (northern Tian Shan, China): Implications for mass-balance estimates, denudation and sedimentation rates in orogenic systems. *Geomorphology* **2014**, *214*, 168–182. [CrossRef]

34. Girod, L.; Nuth, C.; Kääb, A.; McNabb, R.; Galland, O. MMASTER: Improved ASTER DEMs for elevation change monitoring. *Remote Sens.* **2017**, *9*, 704. [CrossRef]
35. James, M.R.; Robson, S.; Smith, M.W. 3-D uncertainty-based topographic change detection with structure-from-motion photogrammetry: Precision maps for ground control and directly georeferenced surveys. *Earth Surf. Process. Landf.* **2017**, *42*, 1769–1788.
[CrossRef]

36. Carretier, S.; Martinod, P.; Reich, M.; Godderis, Y. Modelling sediment clasts transport during landscape evolution. *Earth Surf. Dyn.* **2016**, *4*, 237–251. [CrossRef]

37. Davy, P.; and Lague, D. Fluvial erosion/transport equation of landscape evolution models revisited. *J. Geophys. Res. Earth Surf.* **2009**, *114*. [CrossRef]

38. Barnhart, K.R.; Glade, R.C.; Shobe, C.M.; Tucker, G.E. Terrainbento 1.0: A Python package for multi-model analysis in long-term drainage basin evolution. *Geosci. Model Dev.* **2019**, *12*, 1267–1297. [CrossRef]

39. Howard, A.D.; Dietrich, W.E.; Seidl, M.A. Modeling fluvial erosion on regional to continental scales. *J. Geophys. Res. Solid Earth* **1994**, *99*, 13971–13986. [CrossRef]

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

41. Stock, J.D.; Montgomery, D.R. Geologic constraints on bedrock river incision using the stream power law. *J. Geophys. Res. Solid Earth* **1999**, *104*, 4983–4993. [CrossRef]

42. Field, J. Channel avulsion on alluvial fans in southern Arizona. *Geomorphology* **2001**, *37*, 93–104. [CrossRef]

43. Radiguet, M.; Cotton, F.; Vergnolle, M.; Campillo, M.; Walpersdorf, A.; Cotte, N.; Kostoglodov, V. Slow slip events and strain accumulation in the Guerrero gap, Mexico. *J. Geophys. Res. Solid Earth* **2012**, *117*. [CrossRef]

44. Dal Zilio, L.; Hetényi, G.; Hubbard, J.; Bollinger, L. Building the Himalaya from tectonic to earthquake scales. *Nat. Rev. Earth Environ.* **2021**, *2*, 251–268. [CrossRef]

45. Ren, J.; Zhang, S. Estimation of recurrence interval of large earthquakes on the Central Longmen Shan fault zone based on seismic moment accumulation/release model. *Sci. World J.* **2013**, *2013*, 458341. [CrossRef] [PubMed]

46. Butler, R.; Bond, C. Chapter 9—Thrust systems and contractional tectonics. In *Regional Geology and Tectonics*, 2nd ed.; Scarselli, N., Adam, J., Chiarella, D., Roberts, D.G., Bally, A.W., Eds.; Elsevier: Amsterdam, The Netherlands, 2020; pp. 149–167.

47. Dean, S.L.; Morgan, J.K.; & Fournier, T. Geometries of frontal fold and thrust belts: Insights from discrete element simulations. *J. Struct. Geol.* **2013**, *53*, 43–53. [CrossRef]

48. Faulkner, D.; Jackson, C.; Lunn, R.; Schlische, R.; Shipton, Z.; Wibberley, C.; and Withjack, M. A review of recent developments concerning the structure, mechanics and fluid flow properties of fault zones. *J. Struct. Geol.* **2010**, *32*, 1557–1575. [CrossRef]

49. Fossen, H. *Structural Geology*; Cambridge University Press: Cambridge, UK, 2010.

50. Yang, H.; Quigley, M.; King, T. Surface slip distributions and geometric complexity of intraplate reverse-faulting earthquakes. *GSA Bull.* **2021**. [CrossRef]

51. Adams, B.A.; Whipple, K.X.; Forte, A.M.; Heimsath, A.M.; Hodges, K.V. Climate controls on erosion in tectonically active landscapes. *Sci. Adv.* **2020**, *6*, eaaz3166. [CrossRef]

52. Wolman, M.G.; Miller, J.P. Magnitude and frequency of forces in geomorphic processes. *J. Geol.* **1960**, *68*, 54–74. [CrossRef]