Erosional response of granular material in landscape models

. Tectonics and erosion/sedimentation are the main processes responsible for shaping the earth surface. The link between these processes has strong influence on the evolution of landscapes. One of the tools we have for investigating coupled process models is analogue modeling. Here we contribute to the utility of this tool by presenting laboratory-scaled analogue 10 models of erosion. We explore the erosional response of different materials to imposed boundary conditions, trying to find the composite material that best mimics the behaviour of the natural prototype. The models recreate conditions in which tectonic uplift is no longer active, but there is an imposed fixed slope. On this slope the erosion is triggered by precipitation and gravity, with the formation of channels in valleys and diffusion on hillslope that are function of the analogue material. Using Digital Elevation Models (DEMs) and laser-scan correlation technique, we show model evolution and measure sediment discharge 15 rates. We propose three main components of our analogue material (silica powder, glass microbeads and PVC powder) and we investigate how different proportions of these components affect the model evolution and the development of landscapes. We find that silica powder is the main responsible for creating a realistic landscape in laboratory. Furthermore, we find that varying the concentration of silica powder between 40 wt.% and 50 wt.% (with glass microbeads and PVC powder in the range 35-40 wt.% and 15-20 wt.%, respectively) results in metrics and morphologies that are comparable with those from natural 20 prototypes.


Introduction
Whenever tectonics creates topography, erosion and surface processes act in response to the imposed gradient, tending to reduce topography, and with time, remove all relief.During the last decades a strong theoretical background has been built based on field-and analytical observations (e.g.Howard, 1994;Kirby andWhipple, 2001, 2012;Tucker and Whipple, 2002;Whipple et al., 1999;Whipple andTucker, 1999, 2002) but since natural observations provide only a snapshot of processes acting at different timescales (e.g.Castillo et al., 2014;Cyr et al., 2014;Pederson and Tressler, 2012;Sembroni et al., 2016;Vanacker et al., 2015) a quantitative framing of the existing feedbacks between surface processes and tectonics in modifying the topography remain a difficult task.Analytical, numerical and analogue models are often used by tectonic geomorphologists to improve the understanding of the feedbacks between tectonics and surface processes.Numerical models have the advantage

Experimental Approach
In this study we analyzed four different brittle granular materials to be used as rock analogues for the upper crust: Silica Powder (SP), Glass Microbeads (GM), Crushed Quartz (CQ), PolyVinyl Chloride powder (PVC).These pure materials are suited in five different mixes, in different proportions, while one single material (SP) is tested on its own (Tab.1).The selection has fallen on granular materials because: a) They have the proper physical properties to simulate downscaled crustal rock behavior under laboratory conditions in a natural gravity field (e.g.Davis et al., 1983;Lallemand et al., 1994;Mulugeta and Koyi, 1987;Schreurs et al., 2001Schreurs et al., , 2006Schreurs et al., , 2016;;Storti et al., 2000;Storti and McClay, 1995).As a matter of fact, they obey to Mohr-Coulomb failure criterion, showing strain-hardening prior to failure at peak strength and strain-softening until a stable value is reached (stable friction) (Lohrmann et al., 2003;Schreurs et al., 2006); b) they reproduce reasonable geomorphic features due to developing of erosion/sedimentation processes like incision, mass wasting and diffusive erosion, transport and sedimentation, although important differences in behavior and characteristics have emerged (Graveleau et al., 2011(Graveleau et al., , 2015;;Graveleau and Dominguez, 2008b;Viaplana-Muzas et al., 2019).
In the following, we describe a) scaling to the natural prototype; b) geotechnical characterization of the materials (including geometrical and physical/chemical properties, frictional properties and permeability); c) erosional characterization.We also define which conditions a proper analogue material should satisfy to be used in landscape evolution models.

Mechanical properties
Here we describe the mechanical properties of four granular materials mixed in different proportions.The aims of this analysis are to study how mechanical properties affect erosion style and to define properties that will be used in future works, where the materials will be involved in experiments where active tectonics and erosion act simultaneously.In five experiments the analogue material are mixtures of the previous materials at various percentage (CM1-5, Tab. 1).In Graveleau et al. (2011), the authors display how these different pure materials (excepting CQ, which was not considered in their work) show advantages and disadvantages in responding to erosion and sedimentation, in terms of morphological features developed, and brittle behavior.For example, GM and PVC produce high scarps with almost no channelization when erosion is applied, but realistically reproduce the brittle rheology of the upper crust (Graveleau et al., 2011).SP morphologies scale well with natural landscapes, but the higher strength of the powder induces unreliable structures under deformation.To overcome the above limitations of materials used as single "ingredient", the authors suggested that a mixture of these three granular materials can be the most appropriate choice.Following Graveleau et al. (2011), here we focus on these mixes rather than pure materials.
The latter are still analyzed, highlighting their role in the mix.
We measured geometrical, physical and frictional properties such as grain size and shape, density, porosity and permeability, internal friction angle and cohesion.We measured frictional properties of experimental granular material (internal friction coefficient μ and cohesion C) with a Casagrande shear box.We performed tests for peak and stable friction at variable normal stresses.Density has been measured with a helium pycnometer.The grain size has been estimated via a series of sieves of decreasing openings dimension (from 250 µm to 45 µm).The material passing the 45 µm sieve has been analyzed using sedimentation in a distilled-water tank, with a hydraulic pump for recirculation of water and a thermal control for estimation of water density.We also used a laser diffractometer for checking the reliability of the previous measurements.A qualitative analysis has been carried out using a Scanning Electron Micrograph (SEM) for the shape of grains and composition.

Geometrical and physical/chemical properties
The materials physical properties like grain size, density, porosity and permeability are listed in Table 2.The SP is a very fine powder (D50 = 20 µm), with clasts of different shape and size (Fig. 1): the smallest ones, are elongated and may lay on bigger clast, with a very high roughness.These characteristics require a careful use of this powder, due to the danger for the respiratory system.The density is the highest among the studied components (2661 ± 1 kg m -3 ).We obtain composition of ~95% SiO2, ~3% Al2O3, ~1% K2O and < 1% of Na2O, MgO and CaO (Fig. 2).The CQ has bigger dimensions with respect to SP (D50 = 87 µm), with medium sphericity of the clasts and high roughness.The composition is very similar to silica powder (with ~0.5% of FeO) and so is the density, slightly lower (2588 ± 1 kg m -3 ).GM (D50 = 98 µm) have a very high sphericity and a very low roughness, where density (2452 ± 1 kg m -3 ) is lower than the CQ.The GM qualitative composition is made of ~69% of SiO2, ~15% of Na2O, ~10% of CaO, ~4% of MgO, ~1% of Al2O3 and < 1% of K2O and fluorine.Finally, the PVC (D50 = 181 µm) has a similar shape respect to glass microbeads but less uniform.The density is the lowest among the components (1402 ± 1 kg m -3 ), and this has a strong effect on the erosive properties of the material, as we will show afterwards.We did not perform SEM measurements of the PVC, due to the complexity of its chemical composition ((C2H3Cl)n).The grain size, density and shape of grains in the mixes are function of the percentage of every single material that form it (Tab.2).

Frictional properties
A good crustal analog material must fail following the Mohr-Coulomb failure criterion (e.g.Davis et al., 1983;Davy and Cobbold, 1991;Krantz, 1991): where τ is the shear stress corresponding to the normal stress σ on the failure plane, and μ is coefficient of internal friction defined as  = tan(), with φ angle of internal friction.For geomorphic experiments, where water is added to the system, parameters like μ and C strongly control the evolution of the experiment, because they change with the amount of water.The brittle granular materials typically used in laboratory models have low elasticity and undergo plastic deformation when their yield strength is reached, sliding along discrete fault-analog planes (e.g.Panien et al., 2006;Ritter et al., 2018;Rossi and Storti, 2003;Schreurs et al., 2016).Deforming granular materials satisfy Mohr-Coulomb criterion (Eq.1), that highlights the relationship between shear stress τ and normal stress σ on the failure plane.The criterion typically shows a linear trend for σ in the order of kPa and MPa, but a convex-outward envelope for normal stresses in the order of hundreds of Pa or lower (Schellart, 2000).In this work the normal stress applied are in the 25-200 kPa range.Estimation of normal stress at the base of our models is about 1 kPa.The experimental tests are then only an underestimate of the real mechanical behavior within the models.Nevertheless, the tests provide useful insight into the frictional properties of the analogue materials.For every test (four tests per material) we defined peak friction (the first peak in the shear curve (Fig. 3) reflecting hardening-weakening during strain localization and then fault initiation) and stable friction (the plateau after peak friction representing friction during sliding (e.g.Montanari et al., 2017).A shear box has been used.It consists of a steel box (Rossi and Storti, 2003) split across its middle in two small blocks with an area of 6×6 cm 2 each.The bottom block is fixed, while the top block moves at constant velocity of 0.165 mm min -1 .Two dynamometers record horizontal and vertical displacement.The tests were made in water saturated condition.From every measurement we defined the material internal friction coefficient (μ, slope) and cohesion (C, intercept) for peak and stable friction.Results for every material and mix are listed in Table 2. Between pure materials SP shows the highest values for φ and C (33-40° and 0-8.5 kPa respectively).For GM internal friction angle is about 23-25° and the cohesion is close to 0 or negative, so is not considered in this analysis.PVC shows the same conditions for C, with an internal friction angle between 25° and 32°.The frictional properties for the CQ are similar to SP, with φ between 30-33° and C between 4.5-6.6 kPa.The mixes show a strong variability for φ and C, with average values about 31° and 6 kPa, respectively.
The values of φ and C do not highlight the increasing of SP concentration from CM2 to CM5.

Porosity and permeability
Porosity is defined as the ratio between the volume of voids or pore spaces (VV) and the total volume (VT): (2) The porosity has been computed measuring the volumetric change of a weighted amount of material respect to an ideal condition where no pores are in the samples.We used a vibrating plate looking for the best compaction of material and measured the variation in volume.We used this technique to draw close to the experimental conditions.This procedure has been repeated several times for each composite material.Unfortunately, porosity is dependent on the handling technique, it is thus impossible to precisely control the porosity of the materials during preparation.The values of porosity for single materials and mixes are shown in Table 2.We also measured porosity by weighing the same volume of material in water-saturated condition and after drying it in an air oven.The results are directly comparable between the two different techniques of measurement.GM shows the lowest porosity (0.26) between pure materials while SP and CQ the highest one (0.36-0.37, respectively).As far as the mixes is concerned, only CM1 shows values higher than 0.40, while porosity is around 0.30 from CM2 to CM5.
Permeability represents the ability of a material to transmit fluids.This property has been tested using an oedometer and measuring the velocity of water flowing through the sample.This parameter is essential in controlling the evolution of our models, as it will be explained later in the text.SP has the lowest permeability (3.56•10 -14 m 2 ) while GM has the highest permeability (2.87•10 -12 m 2 ).CM1 shows the highest permeability between the various mixes (7.42•10 -13 m 2 ).From CM2 to CM5 (from 40 and 70 wt.% of SP) the permeability slightly decreases.However, CM3 and CM5 do not strictly follow this trend: the former has the lowest permeability observed in mixes (9.25•10 -14 m 2 ), while the latter show a value comparable with CM2 and CM4 but slightly higher (4.06•10 -13 m 2 ).The permeability values for mixes are then in the order of 10 -13 m 2 .

Erosion laws and erosive properties
In a stream channel, the relationship between channel slope S and contributing area A is often expressed through Flint's Law (Flint, 1974), and takes the form where ks and θ are the steepness and concavity index respectively.The most common erosion law, consistent with slope-area scaling, for channelized processes is a power law function of the contributing area A and surface gradient S, defined as the "stream power" law (e.g.Howard, 1994;Tucker and Whipple, 2002;Whipple and Tucker, 1999) where z is the elevation of the stream channel (i.e.dz/dt elevation trough time), K is the erosional constant (bound up to the erosional efficiency) that contains information about lithology, climate and channel geometry (Howard et al., 1994;Whipple and Tucker, 1999), while m and n are two positive dimensionless exponents, with the ratio m/n (i.e. the concavity index θ) that typically falls between 0.4 and 0.7 (Tucker and Whipple, 2002).This model is better known as a "detachment-limited" stream power model, because in tectonically active regions or in condition of steep topography, the channel erosive power is high and limited by its capacity to detach particles from the bedrock (Tucker and Whipple, 2002;Whipple and Tucker, 2002).It is possible to rewrite Eq. (4) in terms of distance x along the stream using Hack's Law (Hack, 1957) where ka is a scaling coefficient and H is the reciprocal of the Hack's exponent.Combining Eq. (4) and Eq. ( 5) we obtain where κ = Kka m .In our experiments, K, κ, m, n, H should be constant (for the same experiment), due to the homogeneous lithology and constant precipitation rate.
A proper analogue material for landscape evolution models should erode via localized, area-dependent processes (i.e., advection in valleys), diffusion (i.e., on hillslopes) and mass wasting.Primarily, it must develop channelization in response to accumulated flow.This requires that precipitation collects in surface drainage networks, with branching channels in order to be consistent with Hack's Law and slope-area scaling.For the erosional behavior of the composite material, the ratio between precipitation rate and infiltration capacity appears to be the main factor controlling the geomorphological response.If the precipitation rate is higher than the infiltration capacity, the model can develop surface runoff (Graveleau, 2008).Otherwise, the water flows through and inside the model, inducing fast erosion through discrete and rapid events.Fine sand or powders have been typically used for geomorphic experiments (e.g.Babault et al., 2005;Hasbargen and Paola, 2000), so that runoff could develop (i.e., low infiltration capacity due to the grain size).Nevertheless, different materials exhibit different emergent morphological characteristics when precipitation is imposed.Among the pure materials presented above, only SP (or a mix with SP) successfully reproduces linear incision (e.g.Bonnet and Crave, 2006;Graveleau et al., 2011;Schumm and Parker, 1973;Tejedor et al., 2017), while GM, PVC (Graveleau et al., 2011) and CQ erode mostly by diffusion or mass wasting.

Experimental Setup for erosional characterization
For studying the composite materials response to applied boundary conditions (precipitation rate and slope), we developed a new experimental setup depositing the material into a box on an inclined plane under rainfall (precipitation).Both the initial slope for the apparatus and the precipitation rate are kept constant.No kinematic conditions of sidewalls are applied: no active tectonic is thus reproduced in our model.
The experimental setup is made of three different devices (Fig. 4): the box, the rainfall-and the monitoring systems.The material box controls the imposed initial slope, while the rainfall system triggers surface erosion.The evolution of the model is recorded with digital images and a laser scanner.The only forcing applied to the models is due to the gravity acceleration, that allows for the erosion triggered by slope and rainfall.

Box
The box is a plexiglass tank 0.35×0.3×0.05m 3 , filled with the experimental material, water saturated (about 25 wt.%).After pouring the material into the box and leveling, it is left flat at least 12 hours, to avoid prior deformations.The slope of the box is then fixed at 15°, in analogy to what has been done in Graveleau et al. (2011).

Rainfall system
Three nozzles fixed to an aluminum frame produce a high-density fog in which the droplet size is small enough (≤ 100 µm) to avoid rainsplash erosion (e.g.Bonnet and Crave, 2006;Graveleau et al., 2012;Lague et al., 2003;Viaplana-Muzas et al., 2015).The precipitation rate is controlled by both water pressure and the number of sprinklers.In our models the precipitation rate is fixed to 25-30 mm h -1 .The configuration allows for a homogeneous droplets distribution with a spatial variation of about 20%.The precipitation rate induces surface runoff, channel incision and gravity-driven processes that are responsible for the erosion of the model.

Monitoring system
Each experiment is recorded using one camera and a laser scanner.The camera records the model evolution in oblique view.
The laser horizontal and vertical resolution are 0.05 mm and 0.07 mm, respectively.The scans are converted in digital elevation models (DEMs) using MATLAB.DEMs are analyzed with TopoToolbox (Schwanghart and Scherler, 2014) for geomorphological quantifications.Erosion and sediment discharge are computed with ad-hoc MATLAB algorithms (see Data availability for accessing the codes).Stopping the rainfall and letting the surface dry is required to avoid distortions during the laser scan.For the first hour pictures and scans are taken every 15 min, then every 30 min and 60 min, depending on model evolution rate.

Scaling analysis
An analogue model should be scaled by its geometry, kinematic, dynamic and rheological properties (e.g.Hubbert, 1937;Ramberg, 1981).The length scaling factor (h* = hmodel/hnature) commonly used in sandbox experiments for deformation in the upper crust with granular material ranges between 10 -5 and 10 -6 (e.g.Cruz et al., 2008;Davy and Cobbold, 1991;Konstantinovskaya and Malavieille, 2011;Persson et al., 2004;Storti et al., 2000).Hence, 1 cm in the model may be in the range of 1 -10 km in nature.The gravitational acceleration model-to-nature ratio is set to 1 (g* = gmodel/gnature) working in the natural gravitational field.The density of the dry quartz sand or corundum sand employed in literature defines the model-tonature ratio for density (ρ* = ρmodel/ ρnature) to be around 0.5 -0.6.Since the dimensionless coefficient of internal friction is very similar between the analogue material and the natural crustal rocks, the cohesion/body forces scaling factor can be expressed as Typical values of this scaling factor are in the order of 10 -6 , so that 1 Pa in the model would correspond to about 1 MPa in nature (e.g.Buiter, 2012;Graveleau et al., 2012).
The classical approach of analog modelling of convergent wedges is time independent (i.e., the evolution of the model is independent from convergence rate).Here, together with compressional tectonics we also model erosion with the perspective of implementing the tested materials in analog models where tectonics and erosion are coupled.Therefore, following Willett (1999), we introduce a time scaling factor  * which is the ratio between the mass flux added to the system Fin over the mass flux removed Fout.Fin is defined as follow: where vc and h are the convergence velocity and initial thickness of the layers considered, respectively.Fout, is defined as: where K is a constant proportional to bedrock incision efficiency and precipitation rate (Eq. 4 and 6) and L is the wedge width.
Assuming m = 0.5 and n = 1 in Eq. ( 4) (e.g.Whipple and Tucker, 1999), the K parameter has the dimension of time -1 .Therefore, the mass balance can be expressed as the ratio between these two fluxes: Mb is a dimensionless number, so keeping it the same in the experiment as in nature (in the same gravity field) it is possible to derive the scaling factor for time rewriting Eq. ( 4) as and considering that k* has the dimension t -1 , so that where * marks the model over nature dimensionless ratio for every quantity, as defined before.With this scaling factor, considering L* = 10 -5 -10 -6 and vc* = 8 • 10 4 , 1 min in the model corresponds to 3.8 -38 kyr in nature.
Experimentalists are always in search of a perfect dynamic scaling for their models.Scaling all the aspects of geological processes is very difficult to achieve, if not impossible (Reber et al., 2020).For example, using granular materials leads to a length scaling inherently not perfect: grains in the order of 0.1 -1 mm in laboratory, assuming a length scaling factor of 10 -5 , would correspond to 10 to 100 m in nature, which is obviously overestimated.For landscape evolution models even more issues are linked to fluid flow or sediment transport (Paola et al., 2009).Nevertheless, these experiments provide an "unreasonable effectiveness" (Paola et al., 2009) that allows for their interpretation in terms of scaling by similarity (Reber et al., 2020).When the models and their natural prototype behave in a similar way, it is indeed possible to infer information about the prototype studying the processes acting on the model (Reber et al., 2020), because of the scale independency of the processes.

Results
We present six different models carried out with the same boundary conditions (i.e.imposed slope, precipitation rate and experimental time) and different mixed materials (Tab.1).The same models have been conducted multiple times to ensure the reproducibility of the experimental results.The results presented in this paper are just a selection from amongst these repetitions.We tested mixes with different concentration of CQ, SP, GM and PVC, accounting for the differences in erosional responses, starting from what is already known in literature.These materials (except CQ) are also used in Graveleau et al. (2011) for the authors selected mix, named "MatIV", composed of 40% SP, 40% GM, 18% PVC and 2% graphite powder.
The same mix is here represented by CM2, where PVC is 20 wt.% due to the absence of graphite.Therefore, this mix has been set as our reference model for further analysis.From this starting point, we increased the SP concentration (from 40 to 70 wt.%,from CM2 to CM5) lowering the GM and PVC concentration.Two experiments were carried out using only SP (SM1) and a mix with the same proportion of CM2 but with SP replaced by CQ (CM1).
The analysis illustrates the erosional properties showing the influence of different composition on the morphology of the landscape, the river longitudinal profile, the sediment discharge and erosion.All eroded material leaves the system; therefore, sedimentation is not modelled.

Morphology of erosion
In the reference experiment CM2 the rainfall system induces channel incision and triggers mass wasting processes of portion of the analogue materials adjacent to stream channels.Both advection in channels and diffusion processes on hillslopes are present.A well-developed river network evolves 5 minutes after initiation.Single channels coalesce in basins with the increase of erosion and they are separated by sharp ridges.Three main basins are located at the upper part of the model (Fig. 5).The planar surfaces developing close to the lowermost side of the experiment have a slope of about 9°, six degrees lower than the initial imposed slope.The bottom of the valleys and the peak are generally separated by 1-2 cm relief (Fig. 6).CM1 evolution differs substantially from the reference model CM2.No channel incision is observed (Fig. 5), while diffusion and mass wasting are the dominant processes.Two different planar surfaces are formed, separated by a vertical scarp convex upward.An elongated elevated body stands close to the left boundary, related to the boundary effect itself.The planar surfaces have a slope of 12° and 10° for the lower and upper surface, respectively.The 12° slope is reached after 5 min from the beginning of the experiment, where a proto scarp is already formed.Subsequently, the scarp moves backward for about 60 min at a rate > 0.4 cm min -1 , following which the scarp continues retreating but at a slower rate (< 0.1 cm min -1 ).The difference in elevation between the two planar surfaces is 2-3 cm.In the experiment SM1 channel incision is strong and affects all the model.Almost no mass wasting processes are observed.The landscape evolution for this model is similar to CM2.Four main basins are observed (Fig. 5), with a series of smaller basins linked to the major ones.Two of these basins stand on the leftmost part of the model and are separated by two main ridges (Fig. 6).The rightmost basins have a small ridge separating them.The ridges between different basins can attain a slope close or even higher than 90°.The planar surfaces that form at the end of the experiment have a slope between 9° and 10°, 6° or 5° lower than the imposed initial slope.On the slopes bordering the basins several small channels form.Further increasing the SP concentration changes the erosional response of the model (Fig. 5, 6).
Channel incision becomes the main process acting on the model with the SP concentration from 40 wt.% to 50 wt.%.Further increase in the amount of SP produces more and narrower channels (Fig. 5).An anastomose system develops in CM5.In CM3, CM4 and CM5 the morphologies develop after around 10 minutes from the beginning of the experiment and are almost constant through the evolution of the models.No proper basin develops in these models and there is no evidence of diffusive processes on hillslopes.As a matter of fact, swath profiles transverse to the rivers show strong variation in elevation with a small wavelength (Fig. 6).The valleys are sharp and are very close each other and are not very incised.

River longitudinal profiles
The river longitudinal profile represents the variation of stream elevation relatively to the distance from the outlet.In Fig. 7 we show a river profile for each experiment at four different timesteps.The river evolution in the reference model CM2 follows a well-known path, starting from the undisturbed initial slope and arriving at the final profile with a concave-upward shape.
We can also observe how the propagation of the perturbation, from the initial condition, migrates from the outlet to the headwater as a knickpoint separating the transient from the equilibrium channel profile.In the upper parts of the model, the erosion removes up to 3-3.5 cm of material.In CM1 no proper river develops and the main knickpoint defines the entire topography, with two planar surfaces separated by a sharp scarp (Fig. 5).The experiments of CM3, CM4 and CM5 show a common behavior.Channels do not show a concave-upward shape, or maybe only in the uppermost part of the model, while generally straight rivers develop.Nevertheless, we can observe the propagation of the erosion wave from the bottom to the top.As in the other models, the incision is strong, but it does not produce very deep valleys (1.5-2 cm).Finally, in SM1 it is difficult to observe a proper concave upward river profile.Some knickpoints in the earlier stage of river development are later obliterated (green profile in Fig. 7).The incision removes almost 3 cm in the northernmost portions of the model.

Sediment discharge and erosion
Sediment (mass) discharge can be characterized by the amount of material that leaves the model.Keeping the boundary conditions constant in all the experiments, its evolution is only function of the analogue material composition.Sediment discharge plotted over time shows always two main phases (Fig. 8): phase I, fast removal of material from the model; phase II, slower removal of material with a lower discharge rate that is kept constant until the end of the experiment.After an initial period in which the material quickly responds to the boundary condition with a high discharge rate that varies through time (the slope of the solid line in Fig. 8, phase I), the system reaches an equilibrium with an almost constant discharge rate (the slope of the dashed line in Fig. 8, phase II).In the reference model CM2 and in SM1 this occurs when basins reach the dimension of 40-80 cm 2 .Different behaviors are shown by CM3, CM4 and CM5, where the phase I is extremely short (Fig. 8).In the reference model CM2, phase I last at least 80-90 min with a discharge rate around 15 g min -1 , while for phase II is 6 g min -1 (both values are comparable with SM1).In CM1 the discharge rate decreases from phase I to phase II, from ca. 31 g min -1 in the first 60 min to ca. 7 g min -1 respectively.The loss in discharge rate between phase I and phase II is around 76% and corresponds, in time, when the morphological evolution of the experiment significantly decreases, reaching an almost stationary condition.SM1 shows a similar trend, but a late and smoother transition from phase I to phase II happening after 140 min from the beginning of the experiment.The discharge rate is around 17 g min -1 during the phase I, and 6 g min -1 during the phase II.A strong decrease in sediment discharge over time is observed between CM2 and CM3, while from CM3 to CM5 the difference is smaller.In CM3, CM4 and CM5 phase I is very short in time (< 20 min) with a discharge rate that decreases with the increasing of SP concentration (from 19 g min -1 to 13 g min -1 ).Phase II then lasts for the rest of the experiments, with a discharge rate of about 3 g min -1 .
In Fig. 9 we show the evolution of the erosion for the selected mix.In CM2 and SM1 river channels and basins are observed.
In CM2 they are initially wider than in SM1.CM2 develops less basins, and they are less elongated.In CM2 the erosion appears to be more efficient, with a removal of material up to 3 cm depth in the uppermost portion of the model, close to the middle part of it.The erosion in CM1 follows the retreat of the scarps and it is mainly focused from the scarp to 5-6 cm over the outlet (Fig. 9).It erodes at least 3 cm of the model, although no channels form.The erosion is likely homogeneous on the lower planar surface.The channel incision in CM3 is evident and it reaches 2 to 2.5 cm depth.Here the erosion is focused in the channels, and they do not coalesce into basins.In CM4 and CM5 this affects even more the models.The erosion is extremely low, carving the models with incision lower than 2/2.5 cm (even lower in CM5).

Discussion
The experimental setting presented in this paper allows for the investigation of the materials and composite materials response to the applied boundary conditions (15° slope and precipitation rate of 25-30 mm h -1 ).Despite simplifications (i.e., lack of tectonics, vegetation, storms, chemical weathering, seasoning, presence of infrastructures), these models highlight how the composition of the experimental material controls its erosional response.Respect to other works on the same topic (e.g.Bonnet and Crave, 2006;Graveleau et al., 2008Graveleau et al., , 2011)), we focus on how varying the concentration of materials in a mix affect the models.Increasing the concentration of SP respect to GM and PVC results in straighter channels and lower incision.Using CQ instead of SP results in an almost uniform morphology where no proper basins or channels form, and where the erosion is mainly due to fast discrete events.Three main aspects arise from our results: a) SP is a key ingredient to properly model erosion, b) the physical properties of the experimental mix influence the sediment discharge rate and c) the experimental results can be used to better understand how surface processes act in nature.These considerations lay the foundations for choosing the proper analogue material for landscape evolution models.This indeed must satisfy conditions like morphology of the river channels, geomorphic indexes and erosional behavior, that should fall in the range of natural observations.We find that the best analogue material for landscape models settled in this work is represented by the composite materials used in models CM2 and CM3 (40 wt.% SP, 40 wt.%GM, 20 wt.% PVC and 50 wt.%SP, 35 wt.% GM, 15 wt.%PVC, respectively).Even if our compositions are comparable with the ones used in other laboratories (e.g.Graveleau et al., 2011), it is important to define how different combinations of the proposed materials may result in extremely different model evolution.

Comparison with previous works
The study by Graveleau et al. (2011) represents the recent foundation for modelling landscape evolution.Therefore, we start discussing differences and similarities of our results with respect to their study.In Graveleau et al. (2011) the authors tested four pure materials (silica powder, glass microbeads, PVC powder and graphite) and a single composite material (named "MatIV").The composite material is composed by 40% silica powder, 40% glass microbeads, 18% PVC powder and 2% graphite.In our study we did not analyze graphite.We tested the effect of crushed quartz in composite materials, instead.
Concerning pure materials (silica powder, glass microbeads and PVC powder), our estimations and measurements for sphericity, grain size, density and permeability match the ones made by Graveleau et al. (2011) with unavoidable minor differences.We measured higher values of porosity for PVC, while permeability is in the same order of magnitude.Internal friction angles at peak and stable friction here presented are consistently lower for all our materials in comparison with Graveleau et al. (2011).The authors settled the tests at lower normal stresses than our measurements (< 5 kPa and < 250 kPa, respectively).The Mohr-Coulomb failure criterion shows that when low normal stresses are applied to the sample, the failure envelope tends to steepen, inducing values of internal friction angle higher than if measured at higher normal stresses (Schellart, 2000).This could explain the differences in results.
The composite granular material presented by the author ("MatIV") has been proposed in this work as analogue material used in CM2, except for graphite powder replaced by a slightly higher amount of PVC powder (see Results for more details).
Density, porosity and permeability are comparable with what has been measured by Graveleau et al. (2011) for MatIV.The values for peak friction and stable friction measured in this work are comparable to what has been measured by Graveleau et al. (2011).
Erosion of the models also show similar evolution.The same rate for precipitation has been adopted in both works.We can observe strong similarity in the landscape reorganization between "MatIV" and CM2, looking at the frames in the evolution of the models at 15° slope (Fig. 9a-e, 10c in Graveleau et al. (2011) and Fig. 5 of this manuscript).The mass discharges over time (Fig. 9f in Graveleau et al. (2011) and Fig. 8 of this manuscript) for SilPwd (SM1 in this work) and "MatIV" coincide with our curves for SM1 and CM2 (at least for the first 90 min).During our phase II, the mass discharge observed in CM2 grows faster than "MatIV" (Fig. 9f in Graveleau et al. (2011)).The mass discharge rate during phase II in CM2 is 6 g min -1 , while in Graveleau et al. (2011) the discharge rate for "MatIV" is 2.8 g min -1 .
The analytical approach regarding the temporal scaling here proposed (1 min = 3.8 -38 kyr) can be compared with what has been proposed by Graveleau et al. (2011) and Strak et al. (2011).These authors proposed different approaches for time scaling with respect to this manuscript.However, they obtained the same order of magnitude of the model-to-nature scaling factors here proposed.In Graveleau et al. (2011), in a context tectonically quiescent, the authors compared erosion rate in models and in natural tectonically inactive areas, starting from the geometric scaling.They estimated that 1 min in the models corresponds to 4.1 -16.8 kyr in nature.In Strak et al. (2011) the authors compared the models' denudation rate with the one computed for the Weber and the Salt Lake City segments of the Wasatch fault, obtaining that 1 min = 3.9 -22.5 kyr.Both estimates fall in the range for temporal scaling proposed in this manuscript.The convergence of results grants the reliability of the experimental method, carried out independently at two different laboratories and by different working groups.

Silica powder for erosion models
SP is widely used in geomorphic experiments showing a good qualitative response to erosion/sedimentation and developing geomorphic markers that morphologically approximate the natural prototype (e.g.Bonnet and Crave, 2006;Graveleau et al., 2011;Schumm and Parker, 1973;Tejedor et al., 2017).As already stated by Graveleau et al. (2011), pure granular materials such as SP, GM or PVC are not able to fully satisfy the requirements for our analogue models.Pure GM and PVC show a deformational style characterized by a few localized thrusts and backthrusts, with a low surface slope around 10° coherently with the one of convergent margins (Graveleau et al., 2011).However, these materials fail to reproduce a realistic landscape morphology.SP shows the opposite behavior (Graveleau et al., 2011 and references therein).This material allows the models to develop streams, channels, basins and other geomorphic features, whereas deformation produces a high number of thrusts and backthrusts closely spaced, and a taper slope higher than 14°.A weighted mixture of these three components is then needed to fulfill the requirements for a scaled analogue model, in terms of deformation and erosional style.We managed to pin down silica as main component in our composite materials, and we tested two different siliceous materials: CQ and SP.They are

Sediment discharge rate as function of physical properties of analogue materials
Phase I and phase II differ in terms of sediment discharge rate (sdr).Independently from the material, phase I displays higher sdr than phase II.In phase I the models equilibrate with the boundary conditions imposed.The amount of potential energy triggers a fast reorganization of the system.In the first time steps the materials quickly leave the models, until the energy decreases and a new equilibrium is reached.Lowering the model slope toward an equilibrium shape slows down the erosional response of the model, entering in phase II (Fig. 8).In this latter phase the system has reached a balance with the boundary condition, and sdr shows a lower variability for a given model (Fig. 10).Despite these common points, features like the onset of the phases, their duration or the amount of discharged material differ among the models.During phase I all the models show a high variability in sdr (Fig. 10).Nevertheless, CM1 still has the highest mean sdr (red line in Fig. 10), that later decreases from CM2 to CM5.CM1 erodes through fast and discrete processes (e.g.mass wasting), while all the other models show also (or only) channel incision (Fig. 5).From CM2 to CM5 the sdr decreases.Due to the absence of chemical reaction between components and water, the differences in sdr are linked to the physical properties of the materials.In CM1 the subsurface water flow induces collapse of material in catastrophic events.In SM1 this does not happen.The very low permeability of the material inhibits significant subsurface flows.Here the erosion of the material is mainly linked to the ability of water to detach particles from the riverbed and carrying them outside the model.Initially, particles are detached when shear stresses exerted by water overcome the threshold for detachment of the grains in the analogue materials (Howard, 1994).The strength of SP thus controls the rate of incision.Similar considerations for CM2, but here the higher permeability can trigger both channel incision and gravitational processes.Surprisingly, CM3, CM4 and CM5 response to precipitation rate is very different from SM1 or CM2.From CM3 to CM5 the sdr strongly decreases in both phases (Fig. 10), even if components like GM and PVC are added to the composite materials.The properties of these two pure materials would produce an erosional response similar to CM1 (Graveleau et al., 2011).However, when SP is in a proportion higher than 50 wt.% the water capacity of detaching particles strongly decrease, so that even the incision is very shallow (Fig. 9) and so the sdr.No mass wasting processes act on these models, as suggested by the low permeability.We propose that the higher grain size distribution allows SP to fill the voids between the GM and PVC particles, lowering the permeability (Tab.2) and increasing the material resistance with capillarity and electrostatic forces.In phase II the sdr variability is smaller, and the mean values are more representative of the whole sdr.The previous consideration on the role of grain size applies also in this phase, even if sdr is significantly lower.Up to now, we have talked about the balance between shear stress exerted by water and strength of the riverbed, as responsible for incision in the models.In the works from Lague et al. (2003) and Graveleau et al. (2011), the authors acknowledge the presence of an erosion threshold that must be overcome before significant erosion and transport occurs.Graveleau et al. (2011) proposed that the tilted downstream zone observed in the models may be related to the presence of this erosion threshold.We also observe this tilted downstream zone, accounting likely for the same erosional threshold.In the experiments performed by Lague et al. (2003), all the experiments reach a final height about 1 cm.This was independent of the initial condition.They related this limit elevation to an intrinsic threshold to erosion.They also approach the problem analytically, defining erosion laws in which a threshold term is present.Here, we recognized that the erosion threshold is mainly controlled by the mechanical strength of the materials used in the models, together with topographical parameters (e.g.gradient, rivers organization).From CM2 to CM5, the mechanical strength of the material appears to increase with the increasing of the concentration of SP in the mix, producing less incised landscapes, as highlighted by the morphologies (Fig. 5), the swath profiles (Fig. 6), the incision maps (Fig. 9) and the sediment discharge charts (Fig. 8 and 10).

Drainage network morphology
Our models are not meant to simulate any specific landscape, but to explore how material properties influence landscape development.However, despite the unavoidable limitations and simplifications of the model, it is tempting to compare the experimental and natural data.To do so, one can make a comparison to Hack's Law.Hack's Law (Eq.5) can be also written as where L is the length of the channel in a basin and A its drainage area, c is a scaling coefficient and h is the scaling exponent referred to as Hack's exponent.The scaling coefficient c and the scaling exponent h in Eq. 13 are related to ka and H in Eq. 5 by c = ka -1/H and h = 1/H, respectively.Hack's Law represents the relationship between channel length and drainage area, and allows to analyze the geometry of the drainage network.Dodds and Rothman (2000) show that in nature h is in the range 0.44-0.56,while c is between 1.3 and 6.6 (for individual basins compared at their outlets).In our models, SM1 and CM2 show the lowest values of h (< 0.8).The scaling coefficient c for SM1 and CM2 is in the range 1-4 and 0-4, respectively (between the 25 th and 75 th percentile) (Fig. 11).CM4 shows values for c and h close to 1 and in the range 0.8-1, respectively, while CM5 has a slightly larger distribution.CM1 and CM3 show the lowest values for the scaling coefficient c.Comparing the lengtharea scaling of our analog models with observations (Fig. 5) we notice that the models are characterized by a very low degree of branching of the drainage network.Moreover, our calculations of h are systematically higher than 0.5 (Fig. 11).Values of h greater than 0.5 are typically interpreted as indicating basin elongation with increasing size (Rigon et al., 1996).In fact, the drainage basins in our models are typically elongated, especially for CM3, CM4 and CM5.SM1 and CM2 still have high values for h, but lower with respect to the other models.For SM1 and CM2 the basins are morphologically better defined (Fig. 5).

Steepness and concavity index
We use the following metrics for quantifying erosion in both the laboratory and nature: ksn and θ.Both ksn and θ represent a 1:1 metric for lab-to nature comparison.θ is dimensionless while ksn, whose dimensions are a function of a reference concavity θref) has been computed considering the length scaling factor h* = 10 -5 (related to the strength of the granular materials), so that the values for A in analog models in Eq. ( 9) are in m 2 .For calculating ksn in analog models (ksn_MOD) we assumed θref = 0.45, similar to studies on natural landscapes.We analyzed river profiles of phase II of the experiments because this phase is linked to equilibrium of the system.In general, θ_MOD tends to be lower than 0.5 with the exception of CM2 and CM3, due to the straightening of river longitudinal profiles during the model run (Duvall, 2004;Whipple and Tucker, 1999).Despite the scattering of values for θ_MOD, SM1, CM2 and CM3 show average values higher than the other models, from 0.2 to 0.5 (for data between 25 th and 75 th percentile).For ksn_MOD we found that values computed during phase II range generally between 10 and 140 m 0.9 (Fig. 12).The values for ksn_MOD and θ_MOD do not allow for a unique discrimination between the types of erosion affecting the models, in terms of detachment-limited erosion and transport-limited erosion (Tucker and Slingerland, 1997;Tucker and Whipple, 2002;Whipple and Tucker, 2002), but we consider it likely that experiments are often transport-limited rather than detachment limited.The concavity index for detachment-limited streams are typically higher than for transportlimited streams (Brocard and Van der Beek, 2006;Whipple and Tucker, 2002), even if there is some evidence which suggest that this might not always be true (Gasparini, 1998;Massong and Montgomery, 2000;Tarboton et al., 1991).Of our models, CM2 and CM3 show the highest values for θ_MOD, while CM2 and SM1 show the highest values for ksn_MOD.Both ksn_MOD and θ_MOD (Fig. 12) are generally comparable with data coming from natural compilations (e.g.Kirby and Whipple, 2012) or studies on natural rivers in specific mountainous areas (e.g. .The matching of ksn and θ between models and nature supports future development and application of the analog materials tested in this study for modeling landscape evolution.

Conclusions
We used mixes of water-saturated granular materials as analogs for the upper brittle crust analyzing the role played by geometrical and physical properties in landscape evolution models.Our experimental results illustrate how small variations in the composition of an analogue material can strongly affect the evolution of the geomorphological features and the mechanical response of the materials.According with previous works, we find in SP the main component of analogue materials for landscape evolution models, better if mixed with GM and PVC.We can now conclude that: a) granular materials and mixes of them deform following Mohr-Coulomb criterion; b) composite materials with smaller grain size distribution and higher grain size (order of 100-200 μm) do not allow for advection in valleys, due to higher permeability.The sediment (mass) discharge rate is high, and the erosion happens quickly in time; c) composite materials with higher grain size distribution with particles in the order of 10s of μm allow for both channel incision in valleys and diffusion on hillslopes; d) composite materials where the percentage of SP higher or equal to 50 wt.%show high number of channels but with a very low incision.The discharge rate is extremely low, and erosion and incision affect less the model.
e) respect to the other models, SM1 and CM2 show more branching and well-defined basins, while CM2 and CM3 show higher values for concavity index.
The geomorphological observations carried out on the models here presented highlights how SM1, CM2, CM3, show features most similar to natural prototypes.Increasing the SP concentration from 40 to 50 wt.%(CM2 and CM3, respectively) leads to straighter channel better defined.For models coupling tectonics and surface processes, the material used in SM1 is not likely to be adequate, due to its deformational behavior (section 4.2).The Hack's exponent in all models was higher than observed in nature, but SM1 and CM2 exhibited the lowest values.Concavity index for all models tended towards values lower than in nature, except for CM2 and CM3, which showed good agreement with nature.All these considerations suggest that the materials used in models CM2 and CM3 should be implemented for reproducing analogue landscapes.Our findings are in agreement with previous works, and here we also quantified the differences between geomorphological indexes as function of composition of analogue materials, giving a further constrain on the choice of the materials.These mixes will be adopted in contexts of active tectonics in future works.
and physical properties of pure granular materials and mixes.Here we show values for grain size (D50), particle density (ρpart), porosity (γ), permeability (k), cohesion (C) and internal friction angle for peak (ϕ p wet) and stable (ϕ s wet) friction.Sphericity and roughness have been estimated after the acquisition of SEM imaging.The subscripts dry and wet indicate whatever 780 the tests were made with dry materials or water-saturated materials.

Figure 2 :
Figure 2: SEM qualitative analysis of the material composition.On the left is the BackScattered-Electron (BSE) imaging, while on the right the Energy-Dispersive Detector (EDS) spectrum.A qualitative composition is presented for a) silica powder, b) crushed quartz and c) glass microbeads.PVC powder has not been analyzed, due to its complex composition.

Figure 3 :
Figure 3: Ideal output from a shear test.The peak friction corresponds to the high point of the curve, while stable friction is defined by the subsequent plateau.

Figure 4 :
Figure 4: Schematic representation of the experimental setup used for models with only erosion: block of experimental material (35×30×5 cm 3 ), rainfall system (commercial sprinklers) and reclining table.A single camera and a high-definition laser scan provide 795

Figure 5 :
Figure 5: Pictures of the models at 200 min from the beginning of the experiment.The red box is the trace of the swath profiles shown in Fig. 6.The blue line indicates the stream, for every experiment, analyzed in Fig. 7.

Figure 6 :
Figure 6: Swath profiles, transverse to the experiment slope.Profiles are plotted at the same experimental time, at which the system keeps its morphologies almost constant through time (ca.200 min).Location as shown in Fig. 5.

Figure 7 :
Figure 7: Longitudinal profiles of the streams highlighted in Fig. 5.This analysis is performed converting the laser scan into DEM and applying the TopoToolbox tool for MATLAB (Schwanghart and Scherler, 2014).The laser horizontal and vertical resolution 805

Figure 8 :
Figure 8: Cumulative mass (sediment) discharge over time for the six experiments.The solid lines correspond to phase I, while dashed lines to phase II.The black dot highlights the transition from phase I to phase II.

Figure 9 :Figure 10 :
Figure 9: Erosion evolution for the experiments, here represented by the cumulative difference in elevation (Δz) of the same point at consecutive times.Time is indicated in columns.Each row corresponds to a model, where the mix adopted is indicated in the first

Figure 11 :
Figure 11: Values distribution for h and c for Hack's Law, as given in Eq. (13), in the models.The black lines indicate the median, while the bottom and top edges of the green box indicate the 25 th and 75 th percentiles, respectively.The green whiskers outside the box cover the data point at < 25 th and > 75 th percentile that are not considered outliers, here indicated by green crosses.

Figure 12 :
Figure 12: Steepness (ksn) and concavity index (θ) in the experiments and in nature (field data).We use θref = 0.45 for computing ksn.The black dots indicate the median, while the bottom and top edges of the blue box indicate the 25 th and 75 th percentiles, respectively.The thin blue lines cover data <25 th and >75 th percentiles that are not considered outliers, and the outliers are indicated by the blue empty dots.Cyr10 = Cyr et al., 2010; DiB10 = DiBiase et al., 2010; TaW02 = Tucker and Whipple, 2002; Dik15 = Dikshit, 2015;