%% Copernicus Publications Manuscript Preparation Template for LaTeX Submissions
%% ---------------------------------
%% This template should be used for copernicus.cls
%% The class file and some style files are bundled in the Copernicus Latex Package, which can be downloaded from the different journal webpages.
%% For further assistance please contact Copernicus Publications at: production@copernicus.org
%% https://publications.copernicus.org/for_authors/manuscript_preparation.html

%% Please use the following documentclass and journal abbreviations for discussion papers and final revised papers.

%% 2-column papers and discussion papers
\documentclass[hess, manuscript]{copernicus}

%% Journal abbreviations (please use the same for discussion papers and final revised papers)

% Advances in Geosciences (adgeo)
% Advances in Radio Science (ars)
% Advances in Science and Research (asr)
% Advances in Statistical Climatology, Meteorology and Oceanography (ascmo)
% Annales Geophysicae (angeo)
% Archives Animal Breeding (aab)
% ASTRA Proceedings (ap)
% Atmospheric Chemistry and Physics (acp)
% Atmospheric Measurement Techniques (amt)
% Biogeosciences (bg)
% Climate of the Past (cp)
% DEUQUA Special Publications (deuquasp)
% Drinking Water Engineering and Science (dwes)
% Earth Surface Dynamics (esurf)
% Earth System Dynamics (esd)
% Earth System Science Data (essd)
% E&G Quaternary Science Journal (egqsj)
% Fossil Record (fr)
% Geochronology (gchron)
% Geographica Helvetica (gh)
% Geoscience Communication (gc)
% Geoscientific Instrumentation, Methods and Data Systems (gi)
% Geoscientific Model Development (gmd)
% History of Geo- and Space Sciences (hgss)
% Hydrology and Earth System Sciences (hess)
% Journal of Micropalaeontology (jm)
% Journal of Sensors and Sensor Systems (jsss)
% Mechanical Sciences (ms)
% Natural Hazards and Earth System Sciences (nhess)
% Nonlinear Processes in Geophysics (npg)
% Ocean Science (os)
% Primate Biology (pb)
% Proceedings of the International Association of Hydrological Sciences (piahs)
% Scientific Drilling (sd)
% SOIL (soil)
% Solid Earth (se)
% The Cryosphere (tc)
% Web Ecology (we)
% Wind Energy Science (wes)


% \usepackage commands included in the copernicus.cls:
\usepackage[german, english]{babel}
\usepackage{tabularx}
\usepackage{cancel}
\usepackage{multirow}
\usepackage{supertabular}
\usepackage{algorithmic}
\usepackage{algorithm}
\usepackage{amsthm}
\usepackage{float}
\usepackage{subfig}
\usepackage{rotating}
\renewcommand{\thetable}{S\arabic{table}}
\renewcommand{\thefigure}{S\arabic{figure}}
\usepackage{multirow}
\usepackage[table,xcdraw]{xcolor}
\usepackage{hyperref}

\begin{document}

\title{Physically-based model for gully simulation: application to the Brazilian Semiarid Region - Supplement}


% \Author[affil]{given_name}{surname}

\Author[1,2]{Pedro Henrique Lima}{Alencar}
\Author[2]{José Carlos}{de Araújo}
\Author[2]{Adunias dos Santos}{Teixeira}

\affil[1]{TU Berlin, Institut für Ökologie, 10587 Berlin, Germany}
\affil[2]{Federal University of Ceará, Departamento de Engenharia Agrícola, Fortaleza, Brazil}

%% The [] brackets identify the author with the corresponding affiliation. 1, 2, 3, etc. should be inserted.



\runningtitle{Physically-based model for gully simulation}

\runningauthor{P.H. L. Alencar et al.}

\correspondence{Pedro Alencar (pedro.alencar@campus.tu-berlin.de); José Carlos de Araújo (jcaraujo@ufc.br)}

\received{}
\pubdiscuss{} %% only important for two-stage journals
\revised{}
\accepted{}
\published{}

%% These dates will be inserted by Copernicus Publications during the typesetting process.

\graphicspath{./figures/}
\firstpage{1}
\maketitle

\section{Rain data}
Daily rainfall data from the Foundation of Meteorology and Water Resources of Ceará (FUNCEME) were used in the study. The interval of the study was nearly sixty years and to cover the whole period, data from three stations were used:

\begin{enumerate}
    \item From 1958 to 1973: Station Coroatá (5.03º S; 39.33º W)
    \item From 1974 to 1987: Station Boa Viagem (5.12º S; 39.73º W)
    \item From 1988 onward: Station Madalena (4.85º S; 39.57º W)
    \item Validation: Station Uruquê (5.14º S; 39.18º W)
    \item Validation: Station Paus Ferro (4.99º S; 39.48º W)
\end{enumerate}

All stations were equipped with Pluviometer Ville de Paris with a standard opening of 400 cm\textsuperscript{2}. The final rain series is presented in Figure \ref{fig:rain}. To assess the quality of the data series, double mass analysis with other stations in the vicinity was performed, as ilustrated in Figure \ref{fig:dmass}, from which we can see that there is no relevant bias in the data.

\begin{figure} [H]
    \centering
    \includegraphics[width = 8cm]{figures/S1.png}
    \caption{Rain series with monthly (left vertical axes) and yearly (right vertical axes) rainfall. The red line indicates the average anual rainfall (over 600 mm yr\textsuperscript{-1})}
    \label{fig:rain}
\end{figure}

\begin{figure} [H]
    \centering
    \includegraphics[width = 7cm]{figures/S2.png}
    \caption{Double mass diagram. Neighbour stations in a radius of 50 km.} 
    \label{fig:dmass}
\end{figure}

\section{Gully models from literature}

\subsection{Foster and Lane Model (FLM)}

The \citet{Foster1983} ephemeral-gully model aims at explaining “erosion by concentrated flow in farm fields” for single runoff intensive events. The gullies are considered ephemeral as productive farmlands usually provide periodic tillage to diminish or remove the gullies generated by previous events. The model is physically based and uses the Manning equation, mass balance, and shear stress mobilisation; it assumes an equilibrium channel width and the gully evolution in two steps. The first step is the vertical incision when the concentrated overflow starts digging the channel with a constant width. The second step starts after the bottom of the channel reaches a non-erodible layer. Then, the section starts a sideward erosive process, widening until the end of the effective runoff, i.e., with a shear stress below the critical stress. Detachment ratio ($D_r$) and shear stress ($\tau$) are given by the Equations (\ref{eq:dr}) and (\ref{eq:tau}).

\begin{equation}
     D_r=\ K_r\ (\tau-\tau_c\ )
     \label{eq:dr}
\end{equation}
\begin{equation}
    \tau(X)=1.35\gamma R_hS\left[1-\left(1-2\frac{X}{WP}\right)^{2.9}\right]\
    \label{eq:tau}
\end{equation}

In Equation \ref{eq:tau}, $X$ is the position of a point on the channel bed, varying from zero to $WP$ (wet perimeter). $S$ is the longitudinal slope of the channel; $r_h$ the hydraulic radius; and $\gamma$ is the specific gravity of water (assumed 9.81 kN m\textsuperscript{-3}).  

In order to model long-term gullies using \citet{Foster1983} equations, the following assumptions were made: First, all mobilised sediment is carried away by the discharges, i.e., there is no sediment deposition on the channel bed. This assumption was confirmed by field surveys in many sections where very little loose sediment was identified. Secondly, in the long run, the effect of each intense runoff can be piled in a cumulative model of widths/depths layers. This implies that each erosive event does not suffer significant influence of the previous, and the total eroded soil is related only with the energy of each event. The piling process considered all events with runoff. The following figures display the flow charts for the original Foster and Lane Model and how it allows to model multiple events by piling area available in the supplementary material (Fig. \ref{fig:flowflm1} and \ref{fig:flowflm2}). 

\begin{figure} [H]
    \centering
    \includegraphics[width = 15cm]{figures/S3.png}
    \caption{Flowchart of the original Foster and Lane Model (1983)}
    \label{fig:flowflm1}
\end{figure}

\begin{figure} [H]
    \centering
    \includegraphics[width = 15cm]{figures/S4.png}
    \caption{Flowchart of the modified Foster and Lane Model, in order to allow multiple successive events.}
    \label{fig:flowflm2}
\end{figure}

\subsection{Sidorchuk Model (SM)}

The Sidorchuk model \citep{Sidorchuk1999} is both physically and empirically based. It considers mass balance of sediment, shear stress (in terms of critical velocity), soil cohesion and the Manning equation to estimate the cross-section geometry and channel slope. It also uses empirical equations based on field measurement to estimate the flow depth and width. The model gives special attention to the processes involving gully wall transformation, as shown in Equations (\ref{eq:dvrc}) and (\ref{eq:ch}). 

\begin{equation}
    D_{vcr}=\frac{2C_h}{g\ \rho_s}\cos{\left(\varphi\right)}\sin^{-2}{\left[\frac{1}{2}\left(\varphi+\frac{\pi}{2}\right)\right]} 
    \label{eq:dvrc}
\end{equation}

\begin{equation}
    \frac{C_h}{g\ \rho_s\ D_v}=\  \frac{\rho_s-w\rho}{\rho}\tan{\varphi\cos^2{\phi}}-\frac{\sin{2}\phi}{2}
    \label{eq:ch}
\end{equation}
 
In Equations (\ref{eq:dvrc}) and (\ref{eq:ch}), $C_h$  is soil cohesion; $D_v$ the depth incision; $D_{v\!c\!r}$ the critical value of depth for wall failure; $w$ the volumetric soil water content; $\rho_s$ the bulk density of the soil; $\rho$\ the density of water; $g$ gravity's acceleration; $\varphi$\ the soil internal friction angle; and $\phi$\ the wall slope, in degrees. A flow chart of the model is available in the supplementary material (Fig. \ref{fig:flowsm}).

\begin{figure} [h]
    \centering
    \includegraphics[width = 15cm]{figures/S5.png}
    \caption{Flowchart of the Sidorchuk Model (1999).}
    \label{fig:flowsm}
\end{figure}

\section{Model evaluators}

To allow us to assess a more realistic, concerning uncertainty, evaluation of the model’s quality, we implemented the routine proposed by \citet{Ritter2013}. The routine, called FITEVAL, allows the modeller to obtain a confidence interval for the Nash-Sutcliffe Efficiency (NSE). The common values of NSE is, therefore, interpreted as the centre of this interval and of the probabilistic distribution of its values depending on how the data are resampled. The method classifies the model as Acceptable to Very good – NSE $\in$ [0.66, 0.95]; (p-value = 0.05). A conservative interpretation of these results is to understand the lowest values as the minimum state of information, i.e., the one that contains (almost) no unproved hypothesis and does not rely on luck. Therefore, the model can be classified, in the worst case scenario, as Acceptable. The output of the FITEVAL is presented in (Figure \ref{fig:ritter}). On the left we see a scatter of points Computed (modelled) and Observed (measured) points, corresponding to the cross-section areas. Below, we have a list of the percentiles of each category. On the right we have the cumulative probabilistic distribution of values of efficiency and below that, a plot showing the observed (diamonds) and computed (continuous line) in decreasing order.

\begin{figure} [H]
    \centering
    \includegraphics[width = 14cm]{figures/S6.png}
    \caption{Output of FITEVAL \citep{Ritter2013}, indicating the model proposed as, in the worst case scenario, acceptable. On the left we have a scatter plot of modelled and measured data, on the right  the cumulative probability distribution of NSE (top) and the values of measured (black line) and modelled (diamond) in decreasing order (bottom).}
    \label{fig:ritter}
\end{figure}

\section{General data and code}

Following we present basic data and results used in the model and the code used. For more, please access the following link for the permanent repository:

\noindent \textcolor{blue}{\href{https://github.com/PedroAlencarTUB/GullyModel-FLSM}{github.com/PedroAlencarTUB/GullyModel-FLSM}}

% Please add the following required packages to your document preamble:

% If you use beamer only pass "xcolor=table" option, i.e. \documentclass[xcolor=table]{beamer}


\begin{sidewaystable}
\centering
\begin{tabular}{
>{\columncolor[HTML]{EFEFEF}}c 
>{\columncolor[HTML]{EFEFEF}}c cccccc
>{\columncolor[HTML]{EFEFEF}}c 
>{\columncolor[HTML]{EFEFEF}}c 
>{\columncolor[HTML]{EFEFEF}}c 
>{\columncolor[HTML]{EFEFEF}}c }
\hline
\cellcolor[HTML]{EFEFEF} & \cellcolor[HTML]{EFEFEF} &  &  &  &  &  &  & \multicolumn{4}{c}{\cellcolor[HTML]{EFEFEF}Measured} \\ \cline{9-12} 
\multirow{-2}{*}{\cellcolor[HTML]{EFEFEF}Gully} & \multirow{-2}{*}{\cellcolor[HTML]{EFEFEF}Section} & \multirow{-2}{*}{\begin{tabular}[c]{@{}c@{}}Basin Area \\ (m²)\end{tabular}} & \multirow{-2}{*}{\begin{tabular}[c]{@{}c@{}}Perimeter \\ (m)\end{tabular}} & \multirow{-2}{*}{Length (m)} & \multirow{-2}{*}{Kc} & \multirow{-2}{*}{Kf} & \multirow{-2}{*}{\begin{tabular}[c]{@{}c@{}}Conc. Time\\  (s)\end{tabular}} & Area (m²) & Perimeter (m) & Width (m) & Prof \\ \hline
\cellcolor[HTML]{EFEFEF} & S1 & 3200 & 432.5 & 195 & 2.141 & 0.084 & 25 & 0.248 & 2.322 & 2.264 & 0.215 \\
\cellcolor[HTML]{EFEFEF} & S2 & 3080 & 415.5 & 187.8 & 2.096 & 0.087 & 23 & 0.742 & 3.451 & 3.188 & 0.582 \\
\cellcolor[HTML]{EFEFEF} & S3 & 1183 & 337.4 & 179.3 & 2.747 & 0.037 & 22 & 0.684 & 3.734 & 3.558 & 0.498 \\
\cellcolor[HTML]{EFEFEF} & S4 & 491 & 107.7 & 35.9 & 1.361 & 0.381 & 6 & 0.505 & 3.826 & 3.763 & 0.288 \\
\cellcolor[HTML]{EFEFEF} & S5 & 78 & 49.8 & 22.9 & 1.579 & 0.149 & 5 & 0.084 & 1.520 & 1.498 & 0.103 \\
\multirow{-6}{*}{\cellcolor[HTML]{EFEFEF}Gully 1} & S6 & 3060 & 408 & 184.5 & 2.065 & 0.090 & 23 & 0.770 & 2.707 & 2.180 & 0.709 \\ \hline
\cellcolor[HTML]{EFEFEF} & S1 & 2000 & 383 & 123 & 2.398 & 0.132 & 31 & 0.721 & 2.235 & 1.853 & 0.479 \\
\cellcolor[HTML]{EFEFEF} & S2 & 1870 & 379 & 114.3 & 2.454 & 0.143 & 29 & 2.428 & 7.559 & 7.322 & 0.694 \\
\cellcolor[HTML]{EFEFEF} & S3 & 1470 & 336 & 99.7 & 2.454 & 0.148 & 28 & 1.105 & 6.950 & 6.891 & 0.307 \\
\cellcolor[HTML]{EFEFEF} & S4 & 252 & 79.9 & 23.5 & 1.409 & 0.456 & 7 & 0.644 & 3.085 & 2.951 & 0.356 \\
\cellcolor[HTML]{EFEFEF} & S5 & 40 & 43.7 & 12.4 & 1.935 & 0.260 & 4 & 0.008 & 0.488 & 0.479 & 0.033 \\
\cellcolor[HTML]{EFEFEF} & S6 & 1930 & 381 & 119 & 2.428 & 0.136 & 30 & 2.403 & 7.633 & 7.438 & 0.653 \\
\multirow{-7}{*}{\cellcolor[HTML]{EFEFEF}Gully 2} & S7 & 1425 & 320 & 92.2 & 2.374 & 0.168 & 25 & 1.470 & 6.263 & 6.163 & 0.451 \\ \hline
\cellcolor[HTML]{EFEFEF} & S1 & 3500 & 554 & 204 & 2.622 & 0.084 & 40 & 1.126 & 6.459 & 6.321 & 0.521 \\
\cellcolor[HTML]{EFEFEF} & S2 & 3400 & 540 & 193 & 2.593 & 0.091 & 39 & 0.962 & 5.596 & 5.543 & 0.306 \\
\cellcolor[HTML]{EFEFEF} & S3 & 1000 & 248 & 105 & 2.196 & 0.091 & 25 & 0.200 & 2.239 & 2.181 & 0.198 \\
\cellcolor[HTML]{EFEFEF} & S4 & 56 & 49.1 & 13.9 & 1.837 & 0.290 & 4 & 0.103 & 1.351 & 1.303 & 0.117 \\
\cellcolor[HTML]{EFEFEF} & S5 & 2230 & 473 & 185.7 & 2.805 & 0.065 & 38 & 0.511 & 3.098 & 3.016 & 0.306 \\
\cellcolor[HTML]{EFEFEF} & S6 & 4800 & 620 & 213 & 2.506 & 0.106 & 43 & 3.751 & 6.485 & 6.057 & 0.897 \\
\cellcolor[HTML]{EFEFEF} & S7 & 4964 & 639 & 218 & 2.539 & 0.104 & 45 & 2.151 & 5.744 & 4.679 & 0.785 \\
\multirow{-8}{*}{\cellcolor[HTML]{EFEFEF}Gully 3} & S8 & 4712 & 609 & 205 & 2.484 & 0.112 & 41 & 1.543 & 5.935 & 5.327 & 0.844 \\ \hline
\end{tabular}
\label{tab: data1}
\caption{Paper data compilation part 1 - Topography results}
\end{sidewaystable}


\begin{sidewaystable}
\centering
\centering
\begin{tabular}{
>{\columncolor[HTML]{EFEFEF}}c cc
>{\columncolor[HTML]{EFEFEF}}c 
>{\columncolor[HTML]{EFEFEF}}c ccccc
>{\columncolor[HTML]{EFEFEF}}c c}
\hline
\cellcolor[HTML]{EFEFEF} & \multicolumn{2}{c}{FLM-Iav} & \multicolumn{2}{c}{\cellcolor[HTML]{EFEFEF}FLM-I60} & \multicolumn{2}{c}{FLM-I30} & \multicolumn{2}{c}{\cellcolor[HTML]{EFEFEF}FLM-I15} & FLM-λ & SM & FL-SM \\ \cline{2-12} 
\multirow{-2}{*}{\cellcolor[HTML]{EFEFEF}Gully} & Area (m²) & Width (m) & Area (m²) & Width (m) & Area (m²) & Width (m) & \cellcolor[HTML]{EFEFEF}Area (m²) & \cellcolor[HTML]{EFEFEF}Width (m) & Area (m²) & Area (m²) & Area (m²) \\ \hline
\cellcolor[HTML]{EFEFEF} & 0.071 & 0.299 & 0.208 & 0.970 & 0.301 & 1.193 & \cellcolor[HTML]{EFEFEF}0.288 & \cellcolor[HTML]{EFEFEF}1.309 & 0.306 & 0.380 & 0.301 \\
\cellcolor[HTML]{EFEFEF} & 0.283 & 0.476 & 0.862 & 1.423 & 1.047 & 1.720 & \cellcolor[HTML]{EFEFEF}1.106 & \cellcolor[HTML]{EFEFEF}1.877 & 1.094 & 1.685 & 1.047 \\
\cellcolor[HTML]{EFEFEF} & 0.165 & 0.331 & 0.500 & 0.992 & 0.871 & 1.215 & \cellcolor[HTML]{EFEFEF}0.866 & \cellcolor[HTML]{EFEFEF}1.325 & 0.494 & 1.118 & 0.871 \\
\cellcolor[HTML]{EFEFEF} & 0.050 & 0.173 & 0.172 & 0.542 & 0.216 & 0.663 & \cellcolor[HTML]{EFEFEF}0.210 & \cellcolor[HTML]{EFEFEF}0.726 & 0.641 & 0.388 & 0.216 \\
\cellcolor[HTML]{EFEFEF} & 0.006 & 0.061 & 0.023 & 0.222 & 0.028 & 0.274 & \cellcolor[HTML]{EFEFEF}0.039 & \cellcolor[HTML]{EFEFEF}0.298 & 0.043 & 0.082 & 0.028 \\
\multirow{-6}{*}{\cellcolor[HTML]{EFEFEF}Gully 1} & - & - & - & - & 1.170 & 1.669 & \cellcolor[HTML]{EFEFEF}- & \cellcolor[HTML]{EFEFEF}- & 1.248 & 2.173 & 1.170 \\ \hline
\cellcolor[HTML]{EFEFEF} & 0.207 & 0.431 & 0.624 & 1.303 & 1.359 & 1.592 & \cellcolor[HTML]{EFEFEF}1.142 & \cellcolor[HTML]{EFEFEF}1.738 & 2.324 & 2.276 & 2.276 \\
\cellcolor[HTML]{EFEFEF} & 0.278 & 0.401 & 0.855 & 1.199 & 1.369 & 1.467 & \cellcolor[HTML]{EFEFEF}1.474 & \cellcolor[HTML]{EFEFEF}1.556 & 2.029 & 1.968 & 1.968 \\
\cellcolor[HTML]{EFEFEF} & 0.100 & 0.326 & 0.315 & 0.984 & 0.567 & 1.206 & \cellcolor[HTML]{EFEFEF}0.505 & \cellcolor[HTML]{EFEFEF}1.318 & 0.860 & 0.590 & 0.567 \\
\cellcolor[HTML]{EFEFEF} & 0.049 & 0.137 & 0.153 & 0.431 & 0.279 & 0.532 & \cellcolor[HTML]{EFEFEF}0.211 & \cellcolor[HTML]{EFEFEF}0.581 & 0.940 & 0.473 & 0.279 \\
\cellcolor[HTML]{EFEFEF} & 0.002 & 0.060 & 0.005 & 0.154 & 0.006 & 0.187 & 0.007 & 0.203 & 0.014 & 0.006 & 0.006 \\
\cellcolor[HTML]{EFEFEF} & - & - & - & - & 1.202 & 1.287 & - & - & 1.721 & 1.690 & 1.690 \\
\multirow{-7}{*}{\cellcolor[HTML]{EFEFEF}Gully 2} & - & - & - & - & 0.477 & 1.081 & - & - & 0.791 & 0.917 & 0.477 \\ \hline
\cellcolor[HTML]{EFEFEF} & 0.328 & 0.641 & 0.961 & 1.877 & 1.171 & 2.288 & 2.003 & 2.493 & 1.192 & 1.752 & 1.171 \\
\cellcolor[HTML]{EFEFEF} & 0.160 & 0.522 & 0.478 & 1.562 & 1.021 & 1.912 & 0.914 & 2.086 & 1.101 & 0.803 & 1.021 \\
\cellcolor[HTML]{EFEFEF} & 0.048 & 0.242 & 0.152 & 0.740 & 0.254 & 0.909 & 0.234 & 0.991 & 0.273 & 0.289 & 0.254 \\
\cellcolor[HTML]{EFEFEF} & 0.008 & 0.067 & 0.022 & 0.189 & 0.028 & 0.239 & 0.031 & 0.268 & 0.068 & 0.086 & 0.028 \\
\cellcolor[HTML]{EFEFEF} & 0.120 & 0.392 & 0.379 & 1.183 & 0.660 & 1.450 & 0.582 & 1.584 & 0.558 & 0.661 & 0.660 \\
\cellcolor[HTML]{EFEFEF} & - & - & - & - & 2.193 & 2.275 & - & - & 2.626 & 3.596 & 3.596 \\
\cellcolor[HTML]{EFEFEF} & - & - & - & - & 1.832 & 2.318 & - & - & 2.173 & 3.025 & 1.832 \\
\multirow{-8}{*}{\cellcolor[HTML]{EFEFEF}Gully 3} & - & - & - & - & 1.891 & 2.250 & - & - & 2.359 & 3.279 & 1.891 \\ \hline
\end{tabular}
\label{tab: data2}
\caption{Paper data compilation part 2 - Model results}
\end{sidewaystable}

% Please add the following required packages to your document preamble:
% \usepackage{multirow}
% \usepackage[table,xcdraw]{xcolor}
% If you use beamer only pass "xcolor=table" option, i.e. \documentclass[xcolor=table]{beamer}
\begin{table}[h]
\centering
\begin{tabular}{|
>{\columncolor[HTML]{EFEFEF}}c |c|
>{\columncolor[HTML]{EFEFEF}}c |c|
>{\columncolor[HTML]{EFEFEF}}c |c|
>{\columncolor[HTML]{EFEFEF}}c |c|
>{\columncolor[HTML]{EFEFEF}}c |
>{\columncolor[HTML]{EFEFEF}}c |}
\hline
\cellcolor[HTML]{EFEFEF} & \textbf{\begin{tabular}[c]{@{}c@{}}Volume \\ (Total Station)\end{tabular}} & \textbf{\begin{tabular}[c]{@{}c@{}}Volume\\ (UAV)\end{tabular}} & \textbf{\begin{tabular}[c]{@{}c@{}}Error\\ volume\end{tabular}} & \textbf{Length} & \textbf{\begin{tabular}[c]{@{}c@{}}Max. \\ Width\end{tabular}} & \textbf{\begin{tabular}[c]{@{}c@{}}Max. \\ Depth\end{tabular}} & \textbf{\begin{tabular}[c]{@{}c@{}}Channel \\ area\end{tabular}} & \multicolumn{2}{c|}{\cellcolor[HTML]{EFEFEF}\textbf{\begin{tabular}[c]{@{}c@{}}Coordenates UTM\\ (zone: 23S)\end{tabular}}} \\ \cline{2-10} 
\multirow{-2}{*}{\cellcolor[HTML]{EFEFEF}\textbf{Gully}} & (m³) & (m³) & (\%) & (m) & (m) & (m) & (m²) & X (m) & Y (m) \\ \hline
1 & 20 & 22.5 & 11\% & 45.5 & 3.7 & 0.709 & 150 & 445297 & 9449332 \\ \hline
2 & 38 & 33.5 & 13\% & 33.1 & 5.8 & 0.694 & 320 & 444614 & 9447240 \\ \hline
3 & 69 & - & - & 37.4 & 8.4 & 0.844 & 480 & 444883 & 9447479 \\ \hline
\end{tabular}
\label{tab:vols}
\caption{Comparison of gullies dimensions.}
\end{table}

\begin{table}[h]
\centering
\begin{tabular}{cccccc}
\rowcolor[HTML]{C0C0C0} 
\cellcolor[HTML]{C0C0C0} & \cellcolor[HTML]{C0C0C0} & \multicolumn{4}{c}{\cellcolor[HTML]{C0C0C0}\textbf{Soil Class}} \\
\rowcolor[HTML]{C0C0C0} 
\multirow{-2}{*}{\cellcolor[HTML]{C0C0C0}\textbf{Land-use description}} & \multirow{-2}{*}{\cellcolor[HTML]{C0C0C0}\textbf{Condition}} & \textbf{A} & \textbf{B} & \textbf{C} & \textbf{D} \\
Cultuivated land & with conservation treatment & 62 & 71 & 78 & 81 \\
 & \cellcolor[HTML]{EFEFEF}without conservation treatment & \cellcolor[HTML]{EFEFEF}72 & \cellcolor[HTML]{EFEFEF}81 & \cellcolor[HTML]{EFEFEF}88 & \cellcolor[HTML]{EFEFEF}91 \\
Pasture & good & 39 & 61 & 74 & 80 \\
 & \cellcolor[HTML]{EFEFEF}bad & \cellcolor[HTML]{EFEFEF}68 & \cellcolor[HTML]{EFEFEF}79 & \cellcolor[HTML]{EFEFEF}(86) & \cellcolor[HTML]{EFEFEF}89
\end{tabular}
\label{tab:cn}
\caption{CN values by land-use \citep{Chow1988}. In parentesis the CN adopted in this study.}

\end{table}

\newpage
\subsection{Fortran Code}
\footnotesize
\begin{verbatim}
program FLSM

IMPLICIT NONE
CHARACTER arquivo*30,arquivo1*30,arquivo2*30
integer No,i
REAL Q(1500),n,S,Kr,Ds,Pa,Tc,Xnc,Gc,Rn,WP,Rh,Ta,Wn,Weq,Le,Delta,Param,Tmax,tm(1500),t(1500),tn
REAL NEL,Tb,Gcf,Xncf,Xnc1,Xnc2,Gc1,Gc2,Rn1,Rn2,WP1,Lim,Ch,Phi,Ang,Dvcr,g,pi,es,Ws
REAL WP2,Rh1,Rh2,Ta1,Ta2,Dg1,Dg2,Tnc1,Tnc2,Xncf1,Xncf2,gcf1,gcf2,tnf1,tnf2,dgf1,dgf2
REAL wi(1500),wf,dp(1500),depth,DeltaW,PerimF,AreaI,AreaF,Rh_S,Max_wi,Max_W,dW
COMMON /grand1/Gcf,Xncf1,Xncf2
COMMON /grand2/Ch,g,Ds,es,Pa,Phi,depth,Ang,pi

!Phi = internal friction angle of the soil (degrees) > the data can be obtained in the NAVFAC 7.02, 
!based in texture; in Ortiz et al 1989 based in texture, Plasticity intex and porosity or 
!via laboratory experiments
!Ch = Cohesion (Pa) > the data can be obtained in the NAVFAC 7.02, based in texture; 
!in Ortiz et al 1989 based in texture, Plasticity intex and porosity or via laboratory experiments
!Lim = a threshold for aplication of sidorchuk routine
!Q = peak discharge (m3/s)
!tm  = discharge's duration (min)
!n = Manning's number
!S = Slope (abs) - (not in degrees or percentage)
!Kr = rill erodibility factor (s/m)
!Tc = Critical shear stress (Pa)
!Xnc = Normalized distance in the WP where T=Tc (abs)
!Gc = Conveyance function at Xnc (abs)
!Rn = Normalized Hydraulic Radius (abs)
!WP = Wet Perimeter (m)
!Rh = Hydraulic Radius (m)
!Ta = Average shear stress (Pa)
!Wn = Normalized width (abs)
!Weq = equilibrium width (m)
!Er = Erosion rate (kg/m.s)
!Ve = Velocity of movement down (m/s)
!Pse = Weigth of soil eroded in the event (kg)
!Le = Sheet thickness of eroded soil(m)
!Param = A parameter that is repeated throughout the calculations
!t = time (s)
!NEL = Depth of nonerodible layer(m)
!Tb = Tension when the erosion reach the NEL (Pa)
!Gcf = conveyance funcion final (abs)
!Wf = final width of the channel (m)
!tne = time to reach  the NEL (s)

        WRITE(*,*)'******************************************'
        WRITE(*,*)'           Foster_Lane program            '
        WRITE(*,*)'                                          '
        WRITE(*,*)'    Calculates the final total of soil    '
        WRITE(*,*)'       eroded in an ephemeral gully       '
        WRITE(*,*)'        for a time-serie rainfall         '
        WRITE(*,*)'                                          '
        WRITE(*,*)'      Universidade Federal do Ceara       '
        WRITE(*,*)'    Departamento de Engenharia Agricola   '
        WRITE(*,*)'      Doutorado em Engenharia Agricola     '
        WRITE(*,*)'                                          '
        WRITE(*,*)'            Pedro Alencar, 2019           '
        WRITE(*,*)'                                          '
        WRITE(*,*)'******************************************'

No = 947 !Number of events

!1. READS INPUT DATA
!write(*,*)'Insert the number of observed events'
!Read(*,*)No
write(*,'(a)')'Insert the name of the file containing the runoff data (Discharge e Duration): 
' !in the absense of measured data we sugest use the 30-min intensity
Read(*,'(a30)')arquivo1
write(*,'(a)')'Insert the name of the file containing the hillslope and soil data: '
Read(*,'(a30)')arquivo2
WRITE(*,'(a)',advance='no') 'Insert the name of the output file: '
READ(*,'(a30)')arquivo
open(50,file=arquivo1)
    do i=1, No
        read(50,*) Q(i),tm(i)
        write(*,*) Q(i),tm(i)
        t(i)=tm(i)*60.
    end do
close(50)
open(60,file=arquivo2)
    read(60,*)n,S,Tc,kr,Ds,es,NEL,Ch,Phi,Lim
    write(*,*)'The number of Manning of the channel is....',n
    write(*,*)'The declivity of the hillslope is..........',S
    write(*,*)'The critical shear stress is...............',Tc
    write(*,*)'The rill erodibility coefficient is........',Kr
    write(*,*)'The soil Bulk density is...................',Ds
    write(*,*)'The porosity of the soil is................',es
    write(*,*)'The depth of the nonerodible layer is......',NEL
    write(*,*)'The soil cohesion is.......................',Ch
    write(*,*)'The internal friction angle is.............',Phi
    write(*,*)'The threshold for wall erosion is..........',Lim
close(60)
Pa = 9803. !N/m³
g = 9.803 !gravity
i = 1
depth = 0
AreaI = 0
Max_wi = 0
Max_W  = 0
pi = 3.1415926536
Phi = Phi*pi/180.
Ws = 0.
i=1

!2. FOSTER AND LANE EQUATIONS FOR INCISION BEFORE REACHING THE NEL
open(40,file=arquivo)
do while(i.le.No)
    Param = (n*Q(i)/(S**(0.5)))**0.375

    if (depth.lt.NEL) then
        Gc = Pa*S*Param/Tc
        if (Gc.lt.1.79) then
            i=i+1
            else

                Delta = sqrt(3.9429**2 - 4*6.9594*(1/Gc))
                Xnc1 = (3.9429-Delta)/(2*6.9594)
                Xnc2 = (3.9429+Delta)/(2*6.9594)
                Rn1 = -0.8834*Xnc1+0.1395*Xnc1+0.151
                Rn2 = -0.8834*Xnc2+0.1395*Xnc2+0.151
        end if
        if (Rn1.gt.0) then
            WP1 = Param/Rn1**(0.625)
            Rh1 = Rn1*WP1
            Ta1 = Pa*Rh1*S
            Tnc1 = Tc/Ta1
            Gc1 = 1./(Tnc1*(Rn1**0.375))
            Dg1 = abs(Gc-Gc1)
            else
                Gc1 = 0
                Dg1 = abs(Gc-Gc1)
        end if
        if (Rn2.gt.0) then
            WP2 = Param/Rn2**(0.625)
            Rh2 = Rn2*WP2
            Ta2 = Pa*Rh2*S
            Tnc2 = Tc/Ta2
            Gc2 = 1/(Tnc2*(Rn2**(0.375)))
            Dg2 = abs(Gc-Gc2)
            else
                Gc2 = 0
                Dg2 = abs(Gc-Gc2)
        end if
        if (Dg1.gt.Dg2) then
            Xnc = Xnc2
            else
                Xnc = Xnc1
        end if
        Write(*,*)'Xnc ',Xnc
        Wn = -1.4873*Xnc+0.7436
        Rn = -0.8834*Xnc**2+0.1395*Xnc+0.151
        WP = Param/Rn**(0.625)
        Weq = WP*Wn
        Rh = Rn*WP
        Ta = Pa*Rh*S
        Tmax = 1.35*Ta
        if (Tmax.lt.Tc) then
            Write(*,*)'The event didn´t cause erosion.'
            wi(i) = 0.
            dp(i) = 0.
            else
                wi(i) = Weq
                Le = Kr*(Tmax-Tc)*t(i)/Ds
                if (Le.gt.NEL) then
                    dp(i) = NEL
                    else
                        dp(i) = Le
                end if
                depth = depth + dp(i)
                if (Max_wi.lt.Weq) then
                    Max_wi = wi(i)
                    Max_W  = wi(i)
                end if
                AreaI = AreaI+wi(i)*dp(i)
            Write(40,*)'Fase1',i,wi(i),dp(i),Max_wi,AreaI
        end if
    end if

!3. FOSTER AND LANE EQUATIONS FOR INCISION AFTER REACHING THE NEL
    if (depth.ge.NEL) then
        Gcf = Pa*S*Param/Tc
        if (Gcf.lt.1.78) then
            i=i+1
            else
                CALL Newton_NEL
        end if
        tnf1 = 1.35*(1-(1-2.*xncf1)**2.9)
        gcf1 = 1/(tnf1*(xncf1*(1-2*xncf1))**(0.375))
        dgf1 = abs(gcf-gcf1)
        tnf2 = 1.35*(1-(1-2.*xncf2)**2.9)
        gcf2 = 1/(tnf2*(xncf2*(1-2.*xncf2))**(0.375))
        dgf2 = abs(gcf-gcf2)
            if (dgf1.gt.dgf2) then
                xncf = xncf2
                else
                    xncf = xncf1
            end if
        Tb = 1.35*Ta*(1-(1-2*xncf)**2.9)*Pa*S*Param*(xncf*(1-2*xncf))**0.375
        if (Tb.le.Tc)then
            dW = 0
            else
            dW = Kr*(Tb-Tc)/Ds
        end if

        wf = Param*((1-2.*xncf)/(xncf**1.667))**0.375
        if (Max_W.ge.Wf) then
            tn = 0.
            else
                tn = t(i)*dW/(Wf-Max_W)
        end if
        Max_W = (1-exp(-tn))*(Wf-Max_W)+Max_W
        write(40,*)'Fase2',i, Max_W, NEL, dw
    end if
    i = i+1
end do
DeltaW = Max_W-Max_Wi
AreaF = AreaI + DeltaW*depth
PerimF = Max_W + 2.*depth
Rh_S = AreaF/PerimF

!4. WALL TRANSFORMATION
if (AreaF.gt.Lim) then
    Dvcr = (2.*Ch*cos(Phi)/(g*Ds))/(sin(0.5*(Phi+pi/2.)))**2.
    if  (Dvcr.le.NEL) then
        CALL Newton_PHI
        Ws = Max_W +2.*NEL/tan(Ang)
        AreaF = NEL*(Max_W+Ws)/2.
    end if
end if
write(40,*)'The cross section area is: ',AreaF,' m2'
if (depth.lt.NEL) then
    write(40,*)'The depth is: ',depth,'m'
    write(40,*)' '
    write(40,*)'The erosion didn´t reach the non erodible layer.'
    else
        write(40,*)'The depth is: ',NEL,'m'
        write(40,*)' '
        write(40,*)'The erosion reached the non erodible layer.'
    end if
if (Ws.gt.0.) then
    write(40,*)'There is erosion of the walls.'
    write(40,*)'The top width is',Ws,'m and the wall slope is',100.*tan(Ang),'%'
    write(40,*)Ang*180/pi
end if
write(*,*)' '
write(*,*)'Please, see the output file'

end program

Subroutine Newton_NEL
!Subrotina para calcular as raizes da equação de Xncf
implicit none
integer i,j
real x,xa,xb,f,df,erro,gcf,xncf1,xncf2
COMMON /grand1/Gcf,Xncf1,Xncf2
Xa = 0.1
Xb = 0.4
i=0
j=0
erro=1000000.
X=Xa
do while (erro.gt.0.000001)
    f = -321.29*x**6+249.02*x**5+9.059*x**4-58.033*x**3+13.012*x**2+1.7199*x-0.004883*COS(x) -1/Gcf
    df = -321.29*6*x**5+249.02*5*x**4+9.059*4*x**3-58.033*3*x**2+13.012*2*x+1.7199+0.004883*sin(x)

    x = xa - f/df
    erro = abs(x-xa)/xa
    xa=x
    i=i+1
end do
erro=1000000.
x=Xb
do while (erro.gt.0.000001)
    f = -321.29*x**6+249.02*x**5+9.059*x**4-58.033*x**3+13.012*x**2+1.7199*x-0.004883*COS(x) -1/Gcf
    df = -321.29*6*x**5+249.02*5*x**4+9.059*4*x**3-58.033*3*x**2+13.012*2*x+1.7199+0.004883*sin(x)

    x = xb - f/df
    erro = abs(x-xb)/xb
    xb=x
    j=j+1
end do
Xncf1 = xa
Xncf2 = xb
end subroutine

Subroutine Newton_PHI
!Subroutine to calculate wall stable angle
    implicit none

    integer i
    real k1,k2,f1,df1,x,xa,erro,Ch,g,Ds,es,Pa,Phi,depth,Ang,pi
    COMMON /grand2/Ch,g,Ds,es,Pa,Phi,depth,Ang,pi

    i=0
    k1 = Ch/(g*Ds*depth)
    k2 = tan(Phi)*(Ds - es*1000.)/1000.

    erro=1000000.
    x = pi/4.
        do while (erro.gt.0.00000001)

        f1 = k2*(cos(x))**2 - sin(2*x)/2. - k1
        df1 = -2*k2*cos(x) - cos(2*x)

        xa = x - f1/df1
        erro = abs(x-xa)
        x=xa
        i=i+1
    end do
Ang = x
end subroutine

\end{verbatim}

\codedataavailability{Code and data available in the link: https://github.com/PedroAlencarTUB/GullyModel-FLSM}

\begin{acknowledgements}
This study was partly financed by the \textit{Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil} (CAPES) - Finance Code 001, Capes/Print Grant 88881.311770/2018-01 and by the \textit{Edital Universal} CNPq - number 407999/2016-7. Pedro Alencar is funded by the DAAD -- Award Number: 91693642. Many thanks to Prof. Eva Paton for the support and advising at the Technische Universität Berlin.
\end{acknowledgements}


%% REFERENCES

%% The reference list is compiled as follows:

%\begin{thebibliography}{}

\bibliographystyle{copernicus.bst}
\bibliography{references.bib}

%\bibitem[AUTHOR(YEAR)]{LABEL1}
%REFERENCE 1

%\bibitem[AUTHOR(YEAR)]{LABEL2}
%REFERENCE 2

%\end{thebibliography}

%% Since the Copernicus LaTeX package includes the BibTeX style file copernicus.bst,
%% authors experienced with BibTeX only have to include the following two lines:
%%
%% \bibliographystyle{copernicus}
%% \bibliography{example.bib}
%%
%% URLs and DOIs can be entered in your BibTeX file as:
%%
%% URL = {http://www.xyz.org/~jones/idx_g.htm}
%% DOI = {10.5194/xyz}


%% LITERATURE CITATIONS
%%
%% command                        & example result
%% \citet{jones90}|               & Jones et al. (1990)
%% \citep{jones90}|               & (Jones et al., 1990)
%% \citep{jones90,jones93}|       & (Jones et al., 1990, 1993)
%% \citep[p.~32]{jones90}|        & (Jones et al., 1990, p.~32)
%% \citep[e.g.,][]{jones90}|      & (e.g., Jones et al., 1990)
%% \citep[e.g.,][p.~32]{jones90}| & (e.g., Jones et al., 1990, p.~32)
%% \citeauthor{jones90}|          & Jones et al.
%% \citeyear{jones90}|            & 1990



%% FIGURES

%% When figures and tables are placed at the end of the MS (article in one-column style), please add \clearpage
%% between bibliography and first table and/or figure as well as between each table and/or figure.


%% ONE-COLUMN FIGURES

%%f
%\begin{figure}[t]
%\includegraphics[width=8.3cm]{FILE NAME}
%\caption{TEXT}
%\end{figure}
%
%%% TWO-COLUMN FIGURES
%
%%f
%\begin{figure*}[t]
%\includegraphics[width=12cm]{FILE NAME}
%\caption{TEXT}
%\end{figure*}
%
%
%%% TABLES
%%%
%%% The different columns must be seperated with a & command and should
%%% end with \\ to identify the column brake.
%
%%% ONE-COLUMN TABLE
%
%%t
%\begin{table}[t]
%\caption{TEXT}
%\begin{tabular}{column = lcr}
%\tophline
%
%\middlehline
%
%\bottomhline
%\end{tabular}
%\belowtable{} % Table Footnotes
%\end{table}
%
%%% TWO-COLUMN TABLE
%
%%t
%\begin{table*}[t]
%\caption{TEXT}
%\begin{tabular}{column = lcr}
%\tophline
%
%\middlehline
%
%\bottomhline
%\end{tabular}
%\belowtable{} % Table Footnotes
%\end{table*}
%
%%% LANDSCAPE TABLE
%
%%t
%\begin{sidewaystable*}[t]
%\caption{TEXT}
%\begin{tabular}{column = lcr}
%\tophline
%
%\middlehline
%
%\bottomhline
%\end{tabular}
%\belowtable{} % Table Footnotes
%\end{sidewaystable*}
%
%
%%% MATHEMATICAL EXPRESSIONS
%
%%% All papers typeset by Copernicus Publications follow the math typesetting regulations
%%% given by the IUPAC Green Book (IUPAC: Quantities, Units and Symbols in Physical Chemistry,
%%% 2nd Edn., Blackwell Science, available at: http://old.iupac.org/publications/books/gbook/green_book_2ed.pdf, 1993).
%%%
%%% Physical quantities/variables are typeset in italic font (t for time, T for Temperature)
%%% Indices which are not defined are typeset in italic font (x, y, z, a, b, c)
%%% Items/objects which are defined are typeset in roman font (Car A, Car B)
%%% Descriptions/specifications which are defined by itself are typeset in roman font (abs, rel, ref, tot, net, ice)
%%% Abbreviations from 2 letters are typeset in roman font (RH, LAI)
%%% Vectors are identified in bold italic font using \vec{x}
%%% Matrices are identified in bold roman font
%%% Multiplication signs are typeset using the LaTeX commands \times (for vector products, grids, and exponential notations) or \cdot
%%% The character * should not be applied as mutliplication sign
%
%
%%% EQUATIONS
%
%%% Single-row equation
%
%\begin{equation}
%
%\end{equation}
%
%%% Multiline equation
%
%\begin{align}
%& 3 + 5 = 8\\
%& 3 + 5 = 8\\
%& 3 + 5 = 8
%\end{align}
%
%
%%% MATRICES
%
%\begin{matrix}
%x & y & z\\
%x & y & z\\
%x & y & z\\
%\end{matrix}
%
%
%%% ALGORITHM
%
%\begin{algorithm}
%\caption{...}
%\label{a1}
%\begin{algorithmic}
%...
%\end{algorithmic}
%\end{algorithm}
%
%
%%% CHEMICAL FORMULAS AND REACTIONS
%
%%% For formulas embedded in the text, please use \chem{}
%
%%% The reaction environment creates labels including the letter R, i.e. (R1), (R2), etc.
%
%\begin{reaction}
%%% \rightarrow should be used for normal (one-way) chemical reactions
%%% \rightleftharpoons should be used for equilibria
%%% \leftrightarrow should be used for resonance structures
%\end{reaction}
%
%
%%% PHYSICAL UNITS
%%%
%%% Please use \unit{} and apply the exponential notation


\end{document}
