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

: 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 bedrock 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 signiﬁcant 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.


Introduction
Geomorphological approaches consider the well-established principle that landscape morphology is controlled by the interactions between tectonics, sedimentation, and erosion [1][2][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][6][7][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][12][13][14][15][16][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, placements 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.  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.

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

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 microbeads. 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].
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.

Scaling
Since the pioneer paper of Hubbert (1937) [25], many studies have focused on analog models scaling (e.g., [19,20,[26][27][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.   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.

Scaling
Since the pioneer paper of Hubbert (1937) [25], many studies have focused on analog models scaling (e.g., [19,20,[26][27][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 −1 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]: where u is the flow velocity, l is a characteristic length as the flow thickness, ν 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, ν, 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 −1 -cm.s −1 , 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 −5 , 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: 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 −1 [33], leading to a V * of ca. 10 3 -10 5 . These values and the previously calculated L * give a time ratio T * of 6 × 10 −10 -12 × 10 −8 . 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 −1 , we set the device for an analog fault slip rate of 1 cm.h −1 , equivalent to 2.8 mm.yr −1 , in our reference model (Table 1).

Model Monitoring and Data Processing
Two camera acquisition systems monitored the model: 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 timelapse 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.

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) ∂z ∂t where ∂z/∂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]: where κ hill is an erodibility coefficient, S = ∂z ∂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 κ hill . A critical slope S c of 0.3 is set for all the numerical models to be consistent with the internal friction coefficient ϕ = 20 • ± 10 • of analog material (see Supplementary Materials S2). Estimating κ 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 = κ hill × 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 2 /yr. Here, we assumed a model node spacing dx of 50 m. This value leads to a range for κ hill between 2 × 10 −7 and 2 × 10 −3 m.yr −1 .
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] ∂z 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 ξ = q d * ×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 ξ → +∞ (i.e., small v s ) and ξ 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 ξ − q relationships tested by Davy and Lague (2009) [37].
The erosion ω was calculated from the classical power-law type equation 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 κ 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.

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.

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 • 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 × 6000 m model with a node spacing of 50 m. This discretization led to a mesh size of 120 × 120, which was a good compromise between computation time and model resolution of the calculated surface.

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 × 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 × 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 3 m 3 .yr −1 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 re f equaled 1 cm.h −1 , leading to v re f = 2.8 mm.yr −1 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 z = 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 z was then defined as the product v z × T. In the reference analog model, we assumed a return period T re f = 180 s, leading to T re f = 1800 yr for the numerical reference model.

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.

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.

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, L x , and in the y-direction, L y . We also used the fan DEM obtained from Micmac processing to calculate the elevation of the fan apex h a .
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 x = 30.7 ± 3.7 cm, L y = 36.7 ± 2.2 cm, h a = 18.5 ± 2 mm. These values were used in the following section to calibrate the erosion law parameters of the numerical models.
in the x-direction, , and in the y-direction, . We also used the fan DEM obtained from Micmac processing to calculate the elevation of the fan apex ℎ .
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 = 30.7 ± 3.7 cm, = 36.7 ± 2.2 cm, ℎ = 18.5 ± 2mm. These values were used in the following section to calibrate the erosion law parameters of the numerical models.

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 κ 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 κ hill , v s , n, and κ 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, The two vectors obs and σobs 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 κ river = 10 −5 , v s = 1 m.yr −1 , κ 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 bestfitting 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 ± 0.06, L y /L = 0.61 ± 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 ± 0.03, respectively. Next, we defined the likelihood function as  , and = 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 , which is equal to 60 cm and 6 km in analog and numerical models, respectively. In the analog experiments, the fan apex height ℎ is 1.85 cm, which corresponds to an elevation/L ratio of 0.031. 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 values. Indeed, alluvial fan building was associated with slopes <5°, which prevent hillslope processes. Our calibration approach was thus not relevant to constrain .  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.
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 κ hill values. Indeed, alluvial fan building was associated with slopes <5 • , which prevent hillslope processes. Our calibration approach was thus not relevant to constrain κ hill .
Geosciences 2021, 11, x FOR PEER REVIEW 12 of 24 Figure 5. Normalized likelihood distribution of numerical modeling parameters calculated from alluvial fan geometries obtained in eight experimental models (see figure 3). is the river bedrock erodibility. Its unit depends on the water discharge exponent , which is defined as = 0.5 in this study. is the slope exponent. is the net settling velocity.
is the erodibility coefficient associated with the hillslope process. A red cross is associated with the best-fitting model with = 10 , = 1 . , = 8. 10 . and = 2.5.
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 of ca. 1 m.yr -1 ( 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 between 0.5 and 3 and a river bedrock erodibility between 10 and 10 . In more detail, the obtained likelihood distribution showed a semi-logarithmic relationship between and : a low is associated with a high , 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 can compensate for a logarithmic increase in . 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 mor- Figure 5. Normalized likelihood distribution of numerical modeling parameters calculated from alluvial fan geometries obtained in eight experimental models (see Figure 3). κ 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. κ hill is the erodibility coefficient associated with the hillslope process. A red cross is associated with the best-fitting model with κ river = 10 −5 , v s = 1 m.yr −1 , κ hill = 8.10 −4 m.yr −1 and n = 2.5.
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 −1 ( 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 κ river : a low κ 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.

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 −1 and 1.9 mm.yr −1 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)  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 -1 and 1.9 mm.yr -1 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)   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 -1 , and a slope exponent of 2. Hillslope erodibility remained the least constrained parameter. A range 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 −4 , a settling velocity of 1 m.yr −1 , and a slope exponent of 2. Hillslope erodibility remained the least constrained parameter. A range between 5 × 10 −7 m.yr −1 and 2 × 10 −3 m.yr 1 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). between 5 × 10 m.yr -1 and 2 × 10 m.yr 1 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 nonstatic river network with lateral movements of central channels during the simulation (Figure 7b).  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. 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.  To further compare the two modeling approaches, we estimated the temporal evolution 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.
mentary 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 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.
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.  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.

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 re f 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 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 The absolute value of reached up to 10 in analog models. This value is larger than that obtained in the numerical simulation. While was positive in analog models, was negative in numerical simulations in the hanging wall. Analog models showed a more complex pattern in the footwall with > 0 at the center of the neoformed fan and < 0 at its boundaries (Figure 10a). One can also notice that the highest values of numerical were located near the fault scarp (Figure 10b). One explanation of these differences in 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 × 10 , which is of the same magnitude as the elevation difference . 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.  The absolute value of ∆z reached up to 10 −3 in analog models. This value is larger than that obtained in the numerical simulation. While ∆z was positive in analog models, ∆z was negative in numerical simulations in the hanging wall. Analog models showed a more complex pattern in the footwall with ∆z > 0 at the center of the neoformed fan and ∆z < 0 at its boundaries (Figure 10a). One can also notice that the highest values of numerical ∆z were located near the fault scarp (Figure 10b).
One explanation of these differences in ∆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 × 10 −4 , which is of the same magnitude as the elevation difference ∆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.

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 f ault . 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 f ault 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).

Effect of Slip Rate on Surface Morphology
For the reference model, we used a homogeneous slip rate v re f = 2.8 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 re f (Table 1), and the other two are associated with a faster slip rate v = 1. 5 × v re f and v = 3 × v re f . 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 ∆z. The main difference was in the two zones of the fault scarps with high |∆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 × v re f to 3. 6 × v re f 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).

Return Period
The reference return period T re f 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 × T re f and 1.1 × T re f . 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 × v re f .  Table 1  is the fault slip rate used in the reference model (see Table 1). (c) Relationship between the deviation amplitude and the fault slip rate obtained from numerical modeling (green circles). The green line shows the linear trend between these two parameters.  Table 1). (b,c) same as (a), assuming a fault slip rate of 1.5 × v re f and 3 × v re f , 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 re f , 1.5 × v re f , and 3 × v re f , respectively. W is the width of the zone with a high elevation deviation.
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 |∆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 |∆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).

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][47][48][49][50]. In the reference model, the distance d = d re f between the sediment source and the fault is equal to 1000 m ( Figure 2). Three other distances were tested: d re f /2, d re f × 2, and d re f × 3. For all numerical models, we assumed a fault slip rate of 3 × v re f . 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 11. Morphological signature of fault slip rate. (a) Elevation associated with stick-slip behavior assuming the fault slip rate used in the numerical reference model (see Table 1). (b) and (c) same as (a), assuming a fault slip rate of 1.5 × and 3 × , respectively. (d), (e), and (f) same as (a), (b), and (c), assuming a stable sliding. (g), (h), and (i) elevation deviation between slick-slip and stable sliding for a fault slip rate of , 1.5 × , and 3 × , respectively. W is the width of the zone with a high elevation deviation. is the fault slip rate used in the reference model (see Table 1). (c) Relationship between the deviation amplitude and the fault slip rate obtained from numerical modeling (green circles). The green line shows the linear trend between these two parameters. (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 re f 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.

Return Period
The reference return period 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 × and 1.1 × . 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 × . 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 | | 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 | |. 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). Figure 13. Role of the earthquake returns period in the deviation between slick slip and stable sliding morphology. A slip rate of 3 × and modeling duration time of including a fan building period and tectonic forcing were assumed (see Table 1). (a) Red circles show the variation of obtained from numerical modeling with respect to the earthquake return period. The obtained linear trend is given by the red line.
is the return period used in the reference model (see Table 1). (b) Relationship between the deviation amplitude and the fault slip rate obtained from numerical modeling (green circles). The green line shows the linear trend between these two parameters.

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][47][48][49][50]. In the reference model, the distance = between Figure 13. Role of the earthquake returns period in the deviation between slick slip and stable sliding morphology. A slip rate of 3 × v re f 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 re f 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.
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 neoformed alluvial fan.  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.

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 numerical models sharing identical geometry and boundary conditions. We demonstrated the excellent repeatability of the analog models by showing that the main geometrical characteristics 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 illustrate 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. 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)).