%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is just an example/guide for you to refer to when submitting manuscripts to Frontiers, it is not mandatory to use Frontiers .cls files nor frontiers.tex  %
% This will only generate the Manuscript, the final article will be typeset by Frontiers after acceptance.   
%                                              %
%                                                                                                                                                         %
% When submitting your files, remember to upload this *tex file, the pdf generated with it, the *bib file (if bibliography is not within the *tex) and all the figures.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Version 3.4 Generated 2018/06/15 %%%
%%% You will need to have the following packages installed: datetime, fmtcount, etoolbox, fcprefix, which are normally inlcuded in WinEdt. %%%
%%% In http://www.ctan.org/ you can find the packages and how to install them, if necessary. %%%
%%%  NB logo1.jpg is required in the path in order to correctly compile front page header %%%

\documentclass[utf8]{frontiersSCNS} % for Science, Engineering and Humanities and Social Sciences articles
%\documentclass[utf8]{frontiersHLTH} % for Health articles
%\documentclass[utf8]{frontiersFPHY} % for Physics and Applied Mathematics and Statistics articles

%\setcitestyle{square} % for Physics and Applied Mathematics and Statistics articles

\usepackage{url,hyperref,lineno,microtype,subcaption}
\hypersetup{
    colorlinks=true,
    linkcolor=blue,
    filecolor=magenta,      
    urlcolor=cyan,
    citecolor=blue,
    pdftitle={LEM Model},
    pdfpagemode=FullScreen,
}
\usepackage[onehalfspacing]{setspace}

\linenumbers

\usepackage[utf8]{inputenc}
\usepackage{amsmath}
\usepackage[T1]{fontenc}
\graphicspath{ {./images/} }
\usepackage{tikz}
\usetikzlibrary{shapes.geometric}
\usetikzlibrary{shapes.arrows}
\usetikzlibrary{arrows}
\usepackage{array}
\usepackage{tikz}
%\usetikzlibrary{arrows.meta}
\usepackage{rotating}
\usepackage{setspace}
\usepackage{pdflscape,array,booktabs}
\usepackage{placeins}
\usepackage[english]{babel}


% Leave a blank line between paragraphs instead of using \\


\def\keyFont{\fontsize{8}{11}\helveticabold }
\def\firstAuthorLast{Cheraghi {et~al.}} %use et al only if is more than 1 author
\def\Authors{Mohsen Cheraghi\,$^{1,*}$, Andrea Rinaldo\,$^{2,3}$, Graham C. Sander\,$^{4}$, Paolo Perona\,$^{5,6}$, Andrea Cimatoribus\,$^{1}$, Seifeddine Jomaa\,$^{7}$ and D. A. Barry\,$^{1}$,}
% Affiliations should be keyed to the author's name with superscript numbers and be listed as follows: Laboratory, Institute, Department, Organization, City, State abbreviation (USA, Canada, Australia), and Country (without detailed address information such as city zip codes or street names).
% If one of the authors has a change of address, list the new address below the correspondence details using a superscript symbol and use the same symbol to indicate the author in the author list.
\def\Address{$^{1}$Ecological Engineering Laboratory (ECOL), Institute of Environmental Engineering (IIE), School of Architecture, Civil and Environmental Engineering (ENAC), École Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland \\
$^{2}$Ecohydrology Laboratory (ECHO), Institute of Environmental Engineering (IIE), School of Architecture, Civil and Environmental Engineering (ENAC), École Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland\\
$^{3}$Dipartimento di Ingegneria Civile Edile e Ambientale, Universitá di Padova, Padua, Italy\\
$^{4}$School of Architecture, Building and Civil Engineering, Loughborough University, Loughborough, UK\\
$^{5}$Hydraulic Platform-LCH, Civil Engineering Institute (IIC), School of Architecture, Civil and Environmental Engineering (ENAC), École Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland\\
$^{6}$School of Engineering, Institute for Infrastructure and Environment, The University of Edinburgh, UK\\
$^{7}$Department of Aquatic Ecosystem Analysis and Management, Helmholtz Centre for Environmental Research—UFZ, Magdeburg, Germany}

% The Corresponding Author should be marked with an asterisk
% Provide the exact contact address (this time including street name and city zip code) and email of the corresponding author
\def\corrAuthor{Mohsen Cheraghi}

\def\corrEmail {\href{mailto:msn.cheraghi@gmail.com}{msn.cheraghi@gmail.com}}

\begin{document}
\onecolumn
\firstpage{1}

\title[LEM in absence of rills]{Applicability of the landscape evolution model in the absence of rills} 

\author[\firstAuthorLast ]{\Authors} %This field will be automatically populated
\address{} %This field will be automatically populated
\correspondance{} %This field will be automatically populated

\extraAuth{}% If there are more than 1 corresponding author, comment this line and uncomment the next one.
%\extraAuth{corresponding Author2 \\ Laboratory X2, Institute X2, Department X2, Organization X2, Street X2, City X2 , State XX2 (only USA, Canada and Australia), Zip Code2, X2 Country X2, email2@uni2.edu}


\maketitle

\begin{abstract}
\section{}
Despite numerous applications of physically-based models for incised landscapes, their applicability for overland flow on unchanneled surfaces is not known. This work challenges a widely used landscape evolution model for the case of non-uniform rainfall and absence of rills using laboratory flume experiment. Rainfall with an average intensity of 85\,$\text{mm\,h}^{-1}$ was applied for 16\,h during which high resolution laser scans of the morphology were captured. The overland flow was modeled as a network that preserves the water flux for each cell in the discretized domain. This network represented the gravity-driven surface flow and determined the evolution direction. The model was calibrated using the first 8\,h of the experiment and was then used to predict the second 8\,h. The calibrated model predicted, as expected, a smoother surface morphology (and less detailed overland flow network) than that measured. This difference resulted from quenched randomness (e.g., small pebbles) within the experimental soil that emerged during erosion and that were captured by the laser scans. To investigate the quality of the prediction, a low-pass filter was applied to remove the small-scale variability of the surface morphology. This step confirmed that the model simulations captured the main characteristics of the measured morphology. The experimental results were found to satisfy a scaling relation for the exceedance probability of discharge even in absence of rills. However, the model did not reproduce the experimental scaling relation as the detailed surface micro-roughness was not accounted for by the model. A lower cutoff on the scale of applicability of the general landscape evolution equation is thus suggested, complementing other work on the upper cutoff underpinned by runoff-producing areas. 
\tiny
 \keyFont{ \section{Keywords:} Landscape evolution model, flume experiment, absence of rills, non-uniform rainfall, high resolution morphology, spectral analysis, drainage network, scaling laws} %All article types: you may provide up to 8 keywords; at least 5 are mandatory.
\end{abstract}

\section{Introduction}
\label{sec:intro} 
The complexity of natural landscapes reflects the numerous factors involved in their formation such as climate \citep{Han2015,Francipane2015,Hooshyar2019,vanderMeij2020}, chemical and physical properties of the sediments  \citep{Massong2000, Sklar2001, Park2015}, episodic gully erosion \citep{Pazzaglia2015}, different vegetation types \citep{Istanbulluoglu2005,Jeffery2014, Corenblit2015}, tectonic effects \citep{Pedrazzini2016} and the rate of sediment production \citep{Rodriguez-Lloveras2015,  Reusser2015, forte2016}. Landscape formation and evolution are not directly measurable in the field due to the time scales involved. Consequently, suitably designed laboratory experiments are of considerable value in that they provide empirical data that are amenable to testing hypotheses of mechanisms underlying observed geomorphological changes \citep{Paola2009}.

Physically-based landscape evolution models (LEMs) \citep{Willgoose1989,Willgoose1991,Willgoose1992,Howard1994,Whipple1999,Whipple2002} are useful tools to explain the surface geometry in landscapes and laboratory experiments \citep{Mudd2016,Whipple2016,Whipple2017,Sinclair2017a}. In LEMs, the complex fluid-particle interactions within landscapes are described by a governing equation for the surface elevation, with an additional model for surface flow \citep{Chen2014a}. For catchments, the focus of surface flow is the stream/river drainage network, rather than the overland flow. Typically, flow is modeled in a simplified manner that conserves the volume flux at each cell in the landscape \citep{OCallaghan1984,Freeman1991,Quinn1991,Costa-Cabral1994,Tarboton1997}. Broadly speaking, the relative importance of the advective and diffusive processes described by the LEM controls the landscape geometry produced. That is, considering an initially non-incised morphology, localized (channel-forming) landscape incision is favored when advection (or surface shear stress) dominates, whereas more gradually varying landscapes will result when (effective) diffusion dominates.

There are numerous applications of LEMs to understand different features of the surface morphology of natural landscapes \citep{Yang2015,Mudd2016,Whipple2016,Whipple2017,Sinclair2017a,bonetti2020,Hu2021,Kwang2021,Shelef2021,Hooshyar2021,Litwin2022}. For example, \citet{Perron2008a} derived an expression for the distance between first-order valleys and validated the formula via measurements from five different natural landscapes \citep{Perron2009a}.  \citet{Willett2014} showed that in disequilibrium conditions, the drainage divides between basins and migrates towards a steady-state condition and thereby reorganize their structure. The criteria defining the reorganization direction were found by using the steady-state solution of an LEM \citep{Perron2013}. It is noteworthy that the LEM equation in its original formulation (where the local landscape-forming fluxes are  surrogated by total contributing area thus postulating uniform precipitation) corresponds to the leading term of the small-gradient approximation of the general mass-balance equation under general reprametrization invariance \citep{Banavar2001a}. This has consequences, because at steady state optimal channel networks (OCNs) are exact constructions and each local, dynamically accessible minimum proves to be a tree \citep{Rodriguez-Iturbe1997,Rinaldo2014}.   

Similar to their application to natural landscapes, LEMs were also used in the analyses of laboratory experiments. Such experiments permit exploration of different initial conditions, e.g., a land surface composed of uniform and non-cohesive sediment grains, with surface morphology changes induced by rainfall. Advection-dominated setups focus on the evolution of the network structure (i.e., surface incisions) \citep{Hancock2002, Bonnet2003,Hasbargen2003,Lague2003,Bonnet2006, Paola2009, Bonnet2009,Graveleau2012,Rohais2012, Reinhardt2015, Singh2015,Cheraghi2018}. Small raindrop sizes minimize the kinetic energy of raindrop impact, leading to surface morphologies that evolve almost exclusively through shear stress of surface flows. For instance, \citet{Hasbargen2000a} set up an experiment in an elliptical plot (98\,cm $\times$ 87\,cm) subjected to a constant uplift rate. They pointed out that the oscillation of the erosion rate at steady-state was the result of knickpoint migration in the domain. The spatial patterns of landslides and knickpoints in a steady-state landscape were measured by \citet{Bigi2006}, who found a power-law relation between the number of landslides and drainage area. The power-law relation was consistent with the gamma-distribution probability density function of landslides in the field \citep{Malmud2004}. \citet{Sweeney2015} analyzed landscapes subjected to rainfall events with different droplet sizes considering the same initial condition. For the different landscape patterns so-created, they showed that the drainage density decreases with increasing droplet size (i.e., relative increase in diffusion).

There are numerous similar studies that examined rill formation on hillslopes \citep{Parker1977,Bryan1989, Gomez1992,Brunton2000, Romkens2001,Pelletier2003,Raff2004, Rieke-ZappD.H.Nearing2005,Tatard2008,Yao2008,Oliveto2010, Gordon2011,Gordon2012, Stefanon2012,Shit2013, He2014,Bennett2015,Bennett2016,Wu2020,Zongping2021}. In these experiments, droplet sizes are large enough to induce splash-impact erosion, which is manifested in changes in surface morphology. For instance, \citet{Gomez2003} and \citet{Berger2010} tested the minimal energy expenditure theory \citep{Rinaldo1992a} for rill networks considering different slopes, rainfall intensities and initial conditions. They found that, when the land was effectively incised, the total energy dissipation decreased as the rills evolved. \citet{McGuire2013} calibrated experimental data with the shallow water flow equations coupled with a process-based erosion model \citep{Hairsine1991,Hairsine1992a, Hairsine1992b,Simpson2003,Simpson2006}. They found that that the entrainment term differences between the models have a greater effect on morphology predictions of the models than the diffusive term. 

In this work, we report a laboratory experiment with two significant differences to those described above. First, by design rainfall was highly non-uniform in space; and second, the surface morphology did not become incised at all in the experiment (i.e., rills did not form), but still exhibiting self-organized aggregation. \citet{Cheraghi2018} showed that the morphology measured in the laboratory flume exhibits the same scale-invariant power laws during its evolution as found for catchments. In other words, the statistical features of catchment channel networks were found in the unchanneled area when the morphology was captured and analyzed at high resolution. That finding motivated this investigation based on a widely-used catchment-scale LEM \citep{Howard1994,Perron2008a} as a tool to predict the evolution of its unchanneled morphology, i.e., for the situation where diffusive processes dominate such that channels do not form. Previous applications of this LEM considered an incised landscape in which discharge occurs either as overland flow (to channels/rivers), or as flow within the channels. In our experiment, only overland flow occurs, and the surface is never incised. Nevertheless, we model the overland flow as a discharge network, as done in other applications of the LEM. Additionally, we test the ability of the model to simulate directly the measured surface morphology, a test that is not possible at the landscape scale, and was hitherto not attempted at the laboratory experiment scale. 

\section{Materials and methods}
\subsection{Experiment}
The same experimental data set was used previously \citep{Cheraghi2018}, and only essential details are presented here. We applied a nonuniform rainfall with an average intensity of 85\,mm\,h\textsuperscript{-1} on a 2\textnormal{-}m $\times$ 1\textnormal{-}m erosion flume with 5\% slope (Figure\,\ref{fig:flume}). There was no infiltration at the flume's bottom. Sandy sediment was added to the flume to a depth of 15\,cm (Table\,\ref{tab:sand}). 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.

\begin{figure}
	\includegraphics[scale=6]{FiguresJPEG3}
	\centering
	\caption{Schematic of the flume experiment. The sediment elevation was measured from the outlet ($z = 0$) and with respect to the horizon. The relief height (the maximum elevation) was 0.19\,m. Details on the experimental setup are given in the Supplementary materials of \citet{Cheraghi2018}.}
	\label{fig:flume}
\end{figure}

\begin{figure}
	\includegraphics[scale=0.7]{FIG6_rainCountours}
	\centering
	\caption{Rainfall distribution on the flume. The average precipitation and the uniformity coefficient were 85\,mm\,h\textsuperscript{-1} and 26\%, respectively. }
	\label{fig:rain}
	
\end{figure}

\begin{table}[h]
	\caption{Characteristics of the sand used in the experiment.}
	\begin{tabular}{c c c c c}
		\hline
		Bulk density    &     Particle size ($d$) range    &    $d_{50}$                    &   $d\,<\,0.6\,mm$             &  $d\,>\,2\,mm		$\\
		\hline
		1584\,(m\textsuperscript{3}\,kg\textsuperscript{-1})  & 0-6\,mm & 0.53\,mm	 & 70\,(\% mass)       & 12\,(\% mass)      \\
		\hline
	\end{tabular}
	\label{tab:sand}
\end{table}

\subsection{Landscape evolution model}
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 \citep{Howard1994} and is derived similarly to the approach of \citet{Perron2008a} with two differences: (i) it is based on the critical stream power rather than critical shear stress, and (ii) 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:  
\begin{equation}  \label{eq:4_LEM_10}
\frac{\partial z}{\partial t} = D \nabla^2 z - K\,f(Q^mS - \Omega_{cr}) 
\end{equation} 
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(\zeta) = \zeta\,H( \zeta)$ where $H( \zeta)$ is the Heaviside step function. The parameter $\Omega_{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 \citep{Furbish2007a,Dunne2010,Mahmoodabadi2016, Sadeghi2017}. The exponent $m$ takes values in the range 0.41-0.857 \citep{Willgoose1989,Rodriguez-Iturbe1997}. 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 \citep{Rinaldo1992a, Rinaldo1993}. At catchment scales, the first term on the right side of Eq. \ref{eq:4_LEM_10} is dominant for hillslopes and includes different processes such as weathering \citep{Perron2017}, soil creep and rain splash \citep{Culling1960,Culling1963,Culling1965}. 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 \citep{Tucker2010b,Benaichouche2016,Whipple2016,Perron2017,Hancock2017}.  \\

\subsection{Numerical simulation and calibration}
\label{sec:simulation}
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$ mins 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 \citet{Planchon2002}. After determining the flow direction via the d8 algorithm \citep{OCallaghan1984}, the discharge at the $i$th cell, $Q_i$, was calculated using:
\begin{equation}  \label{Eq:d8}
Q_i =\sum_{j=1}^{8} w_{ji}\,Q_j + R_i \,\Delta x\,\Delta y
\end{equation}
\begin{figure}
	\centering
	% Define block styles
	\tikzstyle{decision} = [diamond, draw, fill=blue!20, 
	text width=7em, text badly centered, node distance=3cm, inner sep=0pt]
	\tikzstyle{block} = [rectangle, draw, fill=blue!20, 
	text width=15em, text centered, rounded corners, minimum height=4em]
	\tikzstyle{line} = [draw, -latex']
	\tikzstyle{cloud} = [draw, ellipse,fill=red!20, node distance=3cm,
	minimum height=2em]
	\begin{tikzpicture}[node distance = 3cm, auto]
	% Place nodes
	\node [block ] (domain) {network extraction\\+\\ matrix operations};
	%  \node [cloud, left of=init] (expert) {expert};
	%\node [cloud, right of=init] (system) {system};
	\node [block, below of=domain] (BC1) {outlet: $z (0,0) = 0$ \\ walls: $\frac{\partial z}{\partial n}|_w^t = \frac{\partial z^{e}}{\partial n}|_w^{t}$\\$z_w^t = z_{w+dn}^{t} - \frac{\partial z}{\partial n}|_w^t dn $};
	\node [block, right of=BC1, node distance=8cm] (BC2) {outlet: $z (0,0) = 0$ \\ walls: $\frac{\partial z}{\partial n}|_w^t = \frac{\partial z}{\partial n}|_w^{t - \Delta t}$\\$z_w^t = z_{w+dn}^{t} - \frac{\partial z}{\partial n}|_w^t dn $};
		%  \node [block, below of=BC] (evaluate) {evaluate candidate models};
	\node [decision, below of=BC1] (decide) {$z_w^{t} < z_w^{t - \Delta t}? $};
	%  \node [block, left of=evaluate, node distance=3cm] (update) {update model};
	% Draw edges
	\path [line] (domain) -- (BC1);
	\path [line] (BC1) -- (decide);
	%   	 \path [line] (decide)  -| node (_) (BC) ;
	%    \path [line] (evaluate) -- (decide);
	\path [line] (decide) --++ (4,0) node [near start] {No} -| (BC2) ;
	\path [line] (decide) --++ (-4,0) node [near start] {Yes} |- (domain) ;
	\path [line] (BC2)  --++ (0,3) node{}  |- (domain) ;
	%    \path [line] (decide) -- node {no}(stop);
	%    \path [line,dashed] (expert) -- (init);
	%    \path [line,dashed] (system) -- (init);
	%    \path [line,dashed] (system) |- (evaluate);
	\end{tikzpicture}
	\caption{The wall boundary condition assuming $\vec n$ as the outward unit vector normal to the wall. The term $\frac{\partial z^{e}}{\partial n}|_w^t$ was calculated by linear interpolation between the two consecutive scans before and after the time $t$.}
	\label{fig:Boundary}
	\end{figure}
\noindent where the summation over $j$ refers to the eight cells surrounding the $i$th 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\textsuperscript{-1}) is rainfall intensity on the cell $i$ (Figure\,\ref{fig:rain}) and $\Delta x$ (mm) and $\Delta y$ (mm) are the grid sizes in $x$ and $y$ directions, respectively. Numerical results presented here used cell sizes of 8\,mm $\times$ 8\,mm. 

The numerical solution was obtained by two fractional time steps \citep{Press2007}. The first step used the second-order Runge-Kutta scheme to solve the advection term (Eq. \ref{eq:4_LEM_10}), while the second used the Alternating-Direction Implicit (ADI) method for the diffusion part. This approach was used previously by \citet{Perron2008a}, who explain it in detail. At the side wall, a conditional boundary condition was used as shown in Figure\,\ref{fig:Boundary}. At each time step, after the matrix operations (Runge-Kutta + ADI), the gradient normal to the boundary ($\frac{\partial z}{\partial n}|_w^t$, where $w$ refers to the wall) was calculated using the experimental data ($\frac{\partial z^{e}}{\partial 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 - \Delta t}$). If the condition was not met, the slope remained unchanged (equal to the last time step ($t - \Delta t$)) and, based on that, the new elevation at the boundary was calculated.

An evolutionary algorithm, Borg MOEA \citep{Hadka2013}, 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 \ref{fig:Zmodel}). The model was calibrated using the root-mean-square error of the morphology ($z$) at $t = 8$: 

\begin{equation}
f = \sqrt{\frac{\sum_{i=1}^{N}{(z_{i,model}-z_{i,exp})^2}}{N}}|_{t = 8h}
\end{equation} 


\noindent where $z_{i,model}$ and $z_{i,exp}$ are the modeled and measured elevations for cell $i$. 

\section{Results and discussion}

Here, we critically examine the LEM’s ability to reproduce the experimental observations of \citet{Cheraghi2018}. 

\subsection{Morphological evolution and the corresponding discharge network}
The calibrated parameters are presented in Table\,\ref{tab:parameters}. The average diffusion coefficient ($D$) is also higher than those reported for the field scale (0.16-222.4\,mm\textsuperscript{2}\,h\textsuperscript{-1}) \citep{Benaichouche2016, Martin2000}. This is consistent with the lack of rilling in our experiments, as rill formation is favored for lower diffusivities \citep{Sweeney2015}. 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 \citep{Hancock2002, Bonnet2003,Hasbargen2003,Lague2003,Bonnet2006, Paola2009, Bonnet2009,Graveleau2012,Rohais2012, Reinhardt2015, Singh2015,Sweeney2015}. 
\begin{table}[h]
	\centering
	\caption{Calibrated parameters ($\pm$ 95 \% confidence interval) for the LEM.}
	\begin{tabular}{c c c}
		\hline 
		\multicolumn{1}{r}{} 
		$D$\,(mm\textsuperscript{2}\,h\textsuperscript{-1})                    &   $K$\,(mm\textsuperscript{$-\frac{1}{2}$}\,h\textsuperscript{$-\frac{1}{2}$})               & $\Omega_{cr}$\,(mm\textsuperscript{$\frac{3}{2}$}\,h\textsuperscript{$-\frac{1}{2}$})\\
		\hline
		17571 $\pm$ 50 	 & 0.184997 $\pm$ 0.0024 & 13.0765 $\pm$ 0.60 \\
		\hline
	\end{tabular}
	\label{tab:parameters}
\end{table}

\begin{sidewaysfigure}
	\vbox{\vspace{-35em} \includegraphics[width=\textwidth]{mans2_ver15_1_Zmodelexp}}
	\centering
	\caption{Numerical simulation of the LEM and the experimental data. In both experiment and simulation, the morphology evolves toward the high precipitation area.}
	\label{fig:Zmodel}
\end{sidewaysfigure}

The numerical results are compared with the experimental data in Figure\,\ref{fig:Zmodel}. The model is able to capture the main characteristics of the morphology, i.e., the downstream area ($z<60\,\textrm{mm}$) has a symmetric shape while the upstream ($z\geq60\,\textrm{mm}$) area is being eroded. In agreement with the experiment, the model shows a noticeable erosion and growth of downstream region ($z<60\,\textrm{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 \citep{Polyakov2003, Wang2014, Cheraghi2016,Lisle2017}. 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\,\ref{fig:surface}). Nonetheless, the calibrated LEM can reproduce the main features of the morphology. 

In the LEM (Eq. \ref{eq:4_LEM_10}), the advective term is involved in the computation when the value of $Q^{1/2}S$ exceeds the critical value of $\Omega_{cr}$. At the catchment scale, the advection-dominated area are observed as rivers networks \citep{Howard1994,Perron2008a}. However, in this experiment, there is no visual geomorphological criterion to identify the advective regions as there is no incision.  

\begin{figure}[h]
	\includegraphics[scale=5,trim={0 0.45cm 0 0.3cm},clip]{surface}
	\centering
	\caption{Flume surface at $t$ = 0 (a) and $t$ = 16 h (b). This figure is from \citet{Cheraghi2018}.}
	\label{fig:surface}
\end{figure}
\bigskip
\subsection{Spectral analysis}

\begin{figure}[h]
	\includegraphics[scale=1.5]{E_K}
	\centering
	\caption{Power spectral density of the morphology at $t$ = 8\,h. The error bars of the experimental data are based on the standard deviation from the azimuthally averaged PSD in the wave number domain.}
	\label{fig:fft}
\end{figure}

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 \ref{fig:fft}, 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\textsuperscript{-1}. To check the effect of the differences in the spectra, we applied a low-pass filter (Blackman window, \citet{Blackman1958}) to both  measurements and modeled morphologies, with a cutoff wave number of about 7\,m\textsuperscript{-1}. In Figure \ref{fig:filteredZ}, 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 \ref{fig:surface}b) of the experiment while the model results are much less sensitive to the low-pass filtering. In Figure \ref{fig:filteredNet}, the discharge was calculated via Eq. (\ref{Eq:d8}) 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.

\begin{figure}
	\includegraphics[scale=1]{mans2_ver15_mmmZ}
	\centering
	\caption{The original morphology of and the low-pass filtered morphology at $t$ = 8\,h.}
	\label{fig:filteredZ}
\end{figure}

\begin{figure}
	\includegraphics[scale=1]{mans2_ver15_mmm}
	\centering
	\caption{The original discharge network and the network extracted from the low-pass filtered morphology at $t$ = 8\,h.}
	\label{fig:filteredNet}
\end{figure}

\begin{sidewaysfigure}
	\vbox{\vspace{40em}} \includegraphics[width=\textwidth]{mans2_ver15_1_Qmodelexp}
	\centering
	\caption{Discharge distribution at different times (logarithmic scale). The experimental network was calculated based on the filtered morphology. }
	\label{fig:modelQ}
\end{sidewaysfigure}

The discharge network of the filtered experimental results are compared with the model network in Figure\,\ref{fig:modelQ}. 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 \ref{fig:Zmodel}). 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 \citep{Meena2019}, stream network position \citep{McMaster2002}, and flood modeling \citep{Muthusamy2021} is studied. \citet{Muthusamy2021} 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. 


\bigskip
\subsection{Drainage area network statistical analysis}
Here we examine the ability of the model to reproduce experimental scaling laws derived from the drainage area network. The drainage area at the $i$th cell, $A_i$, was calculated using:
\begin{equation}  \label{Eq:d8_A}
A_i =\sum_{j=1}^{8} w_{ji}\,A_j + \Delta x\,\Delta y
\end{equation} 
\noindent where the summation over $j$ and  $ w_{ji}$ has the same definition as in Eq.\,\ref{Eq:d8}.

\begin{figure}
	\includegraphics[scale=0.9]{mans2_ver15_HACK_exp_Binn}
	\centering
	\caption{Relation between the upstream length and drainage area (Hack's law) for the model and experiment. The model results differ from the (unfiltered) experimental data, but agrees with the experimental data subjected to the low-pass filter.}
	\label{fig:5_Hacks_1_6_8}
\end{figure}
   
\begin{figure}
	\includegraphics[scale=0.9]{mans2_ver15_PDF_AmodelExp}
	\centering
	\caption{Exceedance probability of drainage area for the model, original experiment and its low-pass filtered morphology. The model results and filtered experimental results agree at all times.}
	\label{fig:5_P(A)_1_6_8}
\end{figure} 

\begin{figure}
	\includegraphics[scale=0.9]{mans2_ver15_PDF_LmodelExp}
	\centering
	\caption{Exceedance probability of upstream length the model, original experiment and its low-pass filtered morphology. There is good agreement between the filtered data and the model.}
	\label{fig:5_P(L)_1_6_8}
\end{figure}  

In the analysis of \citet{Cheraghi2018}, 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^{-\beta}$) and the exceedance probability of length ($P(L>l) = l^{-\psi}$) had time-invariant exponent values of $\beta = 0.47$ and $\psi = 0.75$. These findings show that, at the flume-scale and in absence of incision, the surface morphology is statistically similar to catchments \citep{Rinaldo1992a,Rodriguez-Iturbe1997} 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 Figures \,\ref{fig:5_Hacks_1_6_8}-\ref{fig:5_P(L)_1_6_8}. 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 \citep{Rinaldo1992a,Rodriguez-Iturbe1997}. Higher order streams of the network (higher drainage area, $A> 2 \times 10^4$ mm\textsuperscript{2}) 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. 
 
\FloatBarrier

\section{Conclusions}
This study introduced the large scale landscape evolution model (LEM), coupled with the d8 algorithm for shallow overland flow, as a robust simulation tool for describing land surface morphology changes in the absence of rills. By calibrating only three model parameters, the LEM was able to capture morphological changes that evolved under a heterogeneous rainfall. This is a manifestation of a simple and efficient model as highlighted by \citet{Paola2011}. The erosive processes captured by the LEM are broadly characterized as raindrop impact-induced diffusion and surface flow-induced shear stress. Unlike previous investigations, here the surface is continuously covered by overland flow, which is modeled as a flow network.

As described by the LEM, shear stress-driven erosion does not occur until a soil-specific flow rate is exceeded. The LEM-based analysis revealed that even for the considered situation where the morphological evolution is dominated by diffusive processes (and so the surface is not incised), the advective term in the LEM is still necessary to predict the surface evolution.  Similarly to natural river networks at catchment scale and incised surfaces in previously reported laboratory experiments, a power-law relation was observed in the discharge exceedance probability for the experimental data. The intricate surface flow details captured by the detailed surface scans of the experimental flume were, however, not present in the modeled surface. Low-pass filtering of the measured and modeled surfaces markedly improved the agreement between the calculated flow networks. However, in that case, the discharge exceedance scaling relationships were not maintained. In general, the LEM could reproduce the low-pass filtered experimental results, i.e., small scale variability due surface roughness that increased over the course of the experiment is not present in the LEM, and had to be removed to obtain reasonable agreement between the experimental measurements and the model results. 

% For Original Research articles, please note that the Material and Methods section can be placed in any of the following ways: before Results, before Discussion or after Discussion.

\section*{Author Contributions}
Conceptualization, writing the original draft, investigation and modeling were done by Mohsen Cheraghi. D.A. Barry, Andrea Rinaldo and Graham Sander supervised the work and edited the manuscript. Paolo Perona provided the laser scanner and supervised the scanning calibration. Andrea Cimatoribus contributed to the spectral analysis. Seifeddine Jomaa guided the initial phases of experiment and further contributed to analyzing the rainfall patterns.  
\section*{Funding}
Financial support was provided by the Swiss National Science Foundation (200021-144320).

\section*{Acknowledgments}
Jacques Roland Golay, Antoine Wiedmer, Pierre-Alain Hildenbrand, Htet Kyi Wynn and Julian Barry provided essential technical support for the execution of the experiment. We also appreciate the collaboration of the SAGRAVE company, Lausanne, Switzerland.

\section*{Data Availability Statement}
The surface morphology data and the LEM source code are available at \href{https://zenodo.org/record/6025295}{Zenodo}. 
% Please see the availability of data guidelines for more information, at https://www.frontiersin.org/about/author-guidelines#AvailabilityofData

%%% Make sure to upload the bib file along with the tex file and PDF
%%% Please see the test.bib file for some examples of references

\bibliographystyle{frontiersinSCNS_ENG_HUMS} % for Science, Engineering and Humanities and Social Sciences articles, for Humanities and Social Sciences articles please include page numbers in the in-text citations
%\bibliographystyle{frontiersinHLTH&FPHY} % for Health, Physics and Mathematics articles
\bibliography{man2.bib} 

%%% If you are submitting a figure with subfigures please combine these into one image file with part labels integrated.
%%% If you don't add the figures in the LaTeX files, please upload them when submitting the article.
%%% Frontiers will add the figures at the end of the provisional pdf automatically
%%% The use of LaTeX coding to draw Diagrams/Figures/Structures should be avoided. They should be external callouts including graphics.

\end{document}
