The same experimental data set was used previously (Cheraghi et al., 2018), and only essential details are presented here. We applied a nonuniform rainfall with an average intensity of 85 mm h⁻¹ on a 2-m × 1-m erosion flume with 5% slope (Figure 1). There was no infiltration at the flume’s bottom. Sandy sediment was added to the flume to a depth of 15 cm (Table 1). Digital Elevation Models (DEMs) at 0.25, 0.5, 1, 2, 4, 8 and 16 h were captured by a 3D laser scanner with about 4-mm resolution. The rainfall stopped while the morphology was measured and then restarted.

We critically assess the applicability of Fokker–Planck form of landscape evolution model at a scale which is less than a computational cell of catchments. Due to large spatial scales, the model is derived based on simplification of turbulent flows and erosion mechanisms. In this work, the idea is to consider the large-scale model at the plot scale and assess how well the simple LEM can represent the complex fluid-particle interactions during heterogeneous rainfall and in absence of incision. The LEM used here is a modification of Howard’s model (Howard, 1994) and is derived similarly to the approach of Perron et al. (2008) with two differences: 1) it is based on the critical stream power rather than critical shear stress, and 2) instead of the drainage area at large scales, the discharge (Q) is used in the nonlinear part of the advection term. In absence of tectonic effects the model is: ∂ z ∂ t = D ∇ 2 z − K f Q m S − Ω c r (1) where Cartesian coordinates are used. In this equation, z(x, y, t), t, S, are elevation from the horizontal, time and slope, respectively, and f(ζ) = ζ H(ζ) where H(ζ) is the Heaviside step function. The parameter Ω _cr is proportional to the critical stream power by which fluvial sediment transport is initiated. The two parameters D and K are functions of rainfall intensity, droplet size, sediment properties (e.g., density, particle size, cohesion) and surface roughness (Furbish et al., 2007; Dunne et al., 2010; Mahmoodabadi and Sajjadi, 2016; Sadeghi et al., 2017). The exponent m takes values in the range 0.41–0.857 (Willgoose, 1989; Rodriguez-Iturbe and Rinaldo, 2001). We took m = 1/2 as this value was also used in studies on optimal channel networks based on the theory of minimum energy expenditure for river network evolution (Rinaldo et al., 1992, 1993). At catchment scales, the first term on the right side of Eq. 1 is dominant for hillslopes and includes different processes such as weathering (Perron, 2017), soil creep and rain splash (Culling, 1960, 1963, 1965). The second term on the right is usually assumed to model sediment transport within the river network (since elsewhere this term vanishes). However, in our flume-scale experiment, there were no surface incisions, and the overland flow is at all times continuous across the entire flume. Thus, the overland flow in the experiment is modeled as a network, although there is no river network as found in previous applications of the LEM (Tucker and Hancock, 2010; Benaïchouche et al., 2016; Whipple et al., 2016; Hancock et al., 2017; Perron, 2017).

During the initial stage of the experiment, the surface morphology at the flume exit changed rapidly due to the location of the flume drain. Such behavior is not captured by the model. Therefore, the scanned morphology at t = 15 min was used as the initial condition for the numerical simulations. In applying the LEM, at each time step, the pit points were removed using the algorithm of Planchon and Darboux (2002). After determining the flow direction via the d8 algorithm (O’Callaghan and Mark, 1984), the discharge at the ith cell, Q _i , was calculated using: Q i = ∑ j = 1 8 w j i Q j + R i Δ x Δ y (2) where the summation over j refers to the eight cells surrounding the ith cell. The value of the element of the connectivity matrix, w _ji (−) is unit if the cell j flows into cell i, and zero otherwise. R _i (mm h⁻¹) is rainfall intensity on the cell i (Figure 2) and Δx (mm) and Δy (mm) are the grid sizes in x and y directions, respectively. Numerical results presented here used cell sizes of 8 mm × 8 mm.

The numerical solution was obtained by two fractional time steps (Press et al., 2007). The first step used the second-order Runge-Kutta scheme to solve the advection term (Eq. 1), while the second used the Alternating-Direction Implicit (ADI) method for the diffusion part. This approach was used previously by Perron et al. (2008), who explain it in detail. At the side wall, a conditional boundary condition was used as shown in Figure 3. At each time step, after the matrix operations (Runge-Kutta + ADI), the gradient normal to the boundary ( ∂ z ∂ n | w t , where w refers to the wall) was calculated using the experimental data ( ∂ z e ∂ n | w t ) , which were linearly interpolated between scans. Afterwards, the sediment elevations adjacent to the walls ( z w t ) were calculated based on the gradients. In the decision block, the value of ( z w t ) was checked; the condition was that the elevation had to decrease at each time step ( z w t < z w t − Δ t ) . If the condition was not met, the slope remained unchanged (equal to the last time step (t − Δt)) and, based on that, the new elevation at the boundary was calculated.

An evolutionary algorithm, Borg MOEA (Hadka and Reed, 2013), was used to find the optimal parameters. The initial condition of the numerical modeling was the experimental data at t = 0.25 h (results in Figure 4). The model was calibrated using the root-mean-square error of the morphology (z) at t = 8: f = ∑ i = 1 N z i , m o d e l − z i , exp 2 N | t = 8 h (3) where z _i,model and z _i, exp are the modeled and measured elevations for cell i.

The calibrated parameters are presented in Table 2. The average diffusion coefficient (D) is also higher than those reported for the field scale (0.16–222.4 mm² h⁻¹) (Martin, 2000; Benaïchouche et al., 2016). This is consistent with the lack of rilling in our experiments, as rill formation is favored for lower diffusivities (Sweeney et al., 2015). In contrast to the present experimental setup, previous laboratory experiments analyzed with the LEM used fine particles, low rainfall intensities and small droplets, i.e., they were designed to produce rills (Hancock and Willgoose, 2002; Bonnet and Crave, 2003; Hasbargen and Paola, 2003; Lague et al., 2003; Bonnet and Crave, 2006; Bonnet, 2009; Paola et al., 2009; Graveleau et al., 2012; Rohais et al., 2012; Reinhardt and Ellis, 2015; Singh et al., 2015; Sweeney et al., 2015).

The numerical results are compared with the experimental data in Figure 4. The model is able to capture the main characteristics of the morphology, i.e., the downstream area (z < 60 mm) has a symmetric shape while the upstream (z ≥ 60 mm) area is being eroded. In agreement with the experiment, the model shows a noticeable erosion and growth of downstream region (z < 60 mm) at t = 4 and 8 h. Recall that the model was calibrated using measurements at 8 h, and then used to predict the morphology at t = 16 h, where the agreement is satisfactory. The evolution of the surface elevation change in both experiment and model is more evident in the high precipitation area (right hand side of the flume).

Despite simulating the main characteristics of the morphology evolution, some differences are seen between the model and the experiment. These differences are likely due to the local scale fluid-particle and particle-particle interactions that are not accounted for in the model. One of these processes is the armoring effect (Polyakov and Nearing, 2003; Wang et al., 2014; Cheraghi et al., 2016; Lisle et al., 2017). Due to shorter erosion time scales, the fine sediment particles are rapidly removed while the larger particles are deposited on the surface or are not moved at all, resulting in a surface covered by pebbles (Figure 5). Nonetheless, the calibrated LEM can reproduce the main features of the morphology.

In the LEM (Eq. 1), the advective term is involved in the computation when the value of Q ^1/2 S exceeds the critical value of Ω _cr . At the catchment scale, the advection-dominated area are observed as rivers networks (Howard, 1994; Perron et al., 2008). However, in this experiment, there is no visual geomorphological criterion to identify the advective regions as there is no incision.

The differences in the modeled and measured morphologies were quantified by taking 2D Fourier transforms of each. We present the resulting power spectral densities (PSD) for the morphologies at t = 8 h (plots for other times are similar) in Figure 6, where we show results as a function of the horizontal wave number (azimuthally averaged). The results show that measured and modeled PSD agree up to a wave number of 5–9 m⁻¹. To check the effect of the differences in the spectra, we applied a low-pass filter (Blackman window, Blackman and Tukey (1958)) to both measurements and modeled morphologies, with a cutoff wave number of about 7 m⁻¹. In Figure 7, the original morphologies of the model and experiment are compared with the filtered ones. Filtering removes the small scale roughness (the large deposited stones shown in Figure 5B) of the experiment while the model results are much less sensitive to the low-pass filtering. In Figure 8, the discharge was calculated via Eq. 2 using the original and filtered morphologies. For the experiment, the complex flow patterns of the measured surfaces (due to microtopography, small scale roughness and armoring) are removed by filtering. On the other hand, there is little change to the discharge network of the model. In short, after filtering, the flow networks for both the experimental measurements and model results are in good agreement. Note that the cutoff wave number and the difference between drainage network before and after filtering are a function of many factors such as sediment particle distribution, microtopography, grain size, cohesion (which does not exist in this experiment) and precipitation rate.

The discharge network of the filtered experimental results are compared with the model network in Figure 9. In both experiment and model, the network has a dendritic form at t = 0.5 h. Afterwards, the network becomes more concentrated and the flow from the upstream is directed to two main streams. The two streams migrate from the left to the right of the flume as the precipitation proceeds. This migration is responsible for the dynamic change in the erosion pattern of the upstream (Figure 4). The model is able to reproduce the migration of the concentrated flow to the right hand side, which induces the downstream morphological evolution.

An important factor to be investigated in future research would be the scanning resolution. At the catchment-scale, the effect of resolution on the hydrological features such as landslide susceptibility mapping (Meena and Gudiyangada Nachappa, 2019), stream network position (McMaster, 2002), and flood modeling (Muthusamy et al., 2021) is studied. Muthusamy et al. (2021) concluded that the resolution has to be finer than the considered river width for modeling. In the absence of rills, it is highly possible to have similar limitation for the LEM simulation based on the grain sizes and the inter-particle distance after armoring.

Here we examine the ability of the model to reproduce experimental scaling laws derived from the drainage area network. The drainage area at the ith cell, A _i , was calculated using: A i = ∑ j = 1 8 w j i A j + Δ x Δ y (4) where the summation over j and w _ji has the same definition as in Eq. 2.

In the analysis of Cheraghi et al. (2018), the experimental data followed a time-invariant power-law relation (Hack’s law, A = l ^h ) between upstream length (l) and drainage area (A) with the scaling exponent h in the range [0.54–0.6] during the network evolution. Furthermore, even though no rills were observed, the exceedance probability of drainage area (P(A > a) = a ^−β ) and the exceedance probability of length (P(L > l) = l ^−ψ ) had time-invariant exponent values of β = 0.47 and ψ = 0.75. These findings show that, at the flume-scale and in absence of incision, the surface morphology is statistically similar to catchments (Rinaldo et al., 1992; Rodriguez-Iturbe and Rinaldo, 2001) in that the same scalings are observed in each case.

The relation between the upstream length (l) and drainage area (A) and the exceedance probabilities for the LEM, original experiment and its low-pass filtered morphology are presented in Figure 10, Figure 11, Figure 12. The results show that the LEM deviates from the original experimental data for the metrics shown in each plot. However, the results each plot show that the model agrees with the low-pass filtered experimental data. These comparisons further confirm the ability of the model to capture the low-pass filtered experimental morphology, but not the micro roughness. On the other hand, the results for the unfiltered experimental data produce scaling laws that agree with previous results for catchments (Rinaldo et al., 1992; Rodriguez-Iturbe and Rinaldo, 2001). Higher order streams of the network (higher drainage area, A > 2 × 10⁴ mm²) are not affected by the low-pass filtering. As a result, both model and filtered experimental results are in agreement with the original experiment in this range.