%layout with cover
%
%
\documentclass[a4paper, 10pt]{article}
%
%----------------------------
%---------PACKAGES-----------
%----------------------------
%
%---------linenumbers and double spacing
%\usepackage{lineno}
%	\linenumbers
%\usepackage{setspace}
%	\doublespacingb
%
%---------spanish
%\usepackage[latin1]{inputenc}
%\usepackage[spanish, es-tabla]{babel}
%
%---------geometry
\usepackage{geometry}
	\geometry{
		a4paper,
		left=20mm,
		right=20mm,
		top=20mm,
		bottom=20mm,
		}
	\setlength{\columnsep}{5mm}
%
%---------version
\usepackage{mVersion}
\setVersion{1.1}
\increaseBuild
%
%---------math
\usepackage{array} %math arrays
\usepackage{amsmath} %math matrices
\usepackage{amssymb}
\usepackage{cancel} %cancel term in equation
\usepackage{bbm} %bold math 0 
\usepackage{bbold} %bold math characters
\usepackage{multirow} %necessary for bracket in matrix
\usepackage{bigdelim} %necessary for bracket in matrix
\usepackage{arydshln} %\hdashline[2pt/2pt] {c;{2pt/2pt}c;{2pt/2pt}c;{2pt/2pt}c}
%
%---------first page stuff
\usepackage[yyyymmdd,hhmmss,24hr]{datetime}
\usepackage{abstract} %abstract (and \saythanks)
\usepackage{authblk} %affiliations
\usepackage{afterpage} %add blank page \blankpage
	\newcommand\blankpage{
		\null
		\thispagestyle{empty}
		\addtocounter{page}{-1}
		\newpage}
%
%
%---------references
\usepackage{natbib} %references
	\DeclareRobustCommand{\van}[3]{#2} %add \DeclareRobustCommand{\van}[3]{#3} before the bibliography
%
%---------symbols and text
\usepackage{siunitx} %SI units
\sisetup{list-final-separator = {, and }}
\usepackage[hyphens]{url}
%\usepackage{eurosym} %euro symbol
\usepackage{color}
\usepackage{enumitem} %\begin{enumerate}[resume]
%
%---------others
\usepackage{booktabs}
%\usepackage{epigraph} %Epigrafes 
%\usepackage{pdfpages} %para incluir pdf como anejos
%
%---------pdf view
\usepackage{hyperref}
\hypersetup{pdfpagemode=UseNone}
%
%---------header and footer
\usepackage{fancyhdr}
	\pagestyle{fancy}
	\fancyhead[L]{\small{Ill-posedness in 2D River Morphodynamics}}
	\fancyhead[R]{\small{V. Chavarr\'ias et al.}}
%	\fancyfoot[L]{\ttfamily{\small{DRAFT COMPILED ON \date{\today} AT \currenttime}}}
%	\fancyfoot[R]{\thepage}
%
%---------figures
\usepackage{graphics}
\usepackage{graphicx}
%	\graphicspath{{C:/Users/Victor/Dropbox/}} %TOSH
%	\graphicspath{{D:/victorchavarri/Dropbox/}} %TUD
%\usepackage{wrapfig} % Allows in-line images ?
%\usepackage{subfigure}
%\usepackage[section]{placeins}
\usepackage[section,subsection]{extraplaceins}
%
%---------flow chart
\usepackage{tikz}
\usetikzlibrary{shapes.geometric,arrows}
%
%---------external references
\usepackage{xr}
%\externaldocument{chapterI}
%
%---------tables
%\usepackage{multirow} %columns spaning multiple rows
%\usepackage{multicol} %rows spaning multiple columns
%
%---------minitoc for table of content in every chapter
%\usepackage{minitoc} 
%\ProvidesFile{spanish.mld}[1999/03/16]
%% Spanish titles for minitoc.sty
%\def\mtctitle{\'Indice General}%	
%\def\mlftitle{\'Indice de Figuras}%
%\def\mlttitle{\'Indice de Tablas}%
%\def\ptctitle{\'Indice General}%
%\def\plftitle{\'Indice de Figuras}%
%\def\plttitle{\'Indice de Tablas}%
%\def\stctitle{\'Indice General}%
%\def\slftitle{\'Indice de Figuras}%
%\def\slttitle{\'Indice de Tablas}%
%
%----------------------------
%---------NEW COMMANDS
%----------------------------
%
%\newcommand{\HRule}{\rule{\linewidth}{0.5mm}}
%\addto\captionsspanish{\renewcommand{\chaptername}{}}
%\def\?#1{}
%
%--------------------------------------------------------------------------
%
%
\begin{document}
%
%----------------------------
%---------TITLE
%----------------------------
%
\begin{titlepage}
\null
\vfill
\begin{flushright}
\huge{Ill-posedness in 2D Mixed Sediment River Morphodynamics}\\
\vspace{5mm}
\Large{Report}\\
\vspace{10mm}
\end{flushright}
\end{titlepage}

\blankpage
%
\pagenumbering{roman}
%\date{\today}
\author{}
\date{}
%\author[1]{V\'ictor Chavarr\'ias} %\thanks{Corresponding author: \ttfamily{v.chavarriasborras@tudelft.nl}}
%\author[2]{Willem Ottevanger}
%\author[3]{Ralph Schielen}
%\author[1]{Robert Jan Labeur}
%\author[1]{Astrid Blom}
%
%
%\affil[1]{\small{Faculty of Civil Engineering and Geosciences, Delft University of Technology, Delft, The Netherlands.}}
%\affil[2]{\small{Deltares, Delft, The Netherlands.}}
%\affil[3]{\small{Rijkswaterstaat, Lelystad, The Netherlands.}}

\title{
\begin{flushright}
Authors: \\
\vspace{10mm}
\large{
V\'ictor Chavarr\'ias $^1$ \\
Willem Ottevanger $^2$ \\
Ralph Schielen $^3$ \\
%Robert Jan Labeur $^1$ \\
Astrid Blom $^1$ \\
}
\vspace{10mm}
\small{$^1$ Faculty of Civil Engineering and Geosciences, Delft University of Technology, Delft, The Netherlands.} \\
\small{$^2$ Deltares, Delft, The Netherlands.} \\
\small{$^3$ Rijkswaterstaat, Lelystad, The Netherlands.} \\
\end{flushright}
%\vspace{10mm}
%Version: \version \quad \small{DRAFT COMPILED ON  AT \currenttime}\\
}
\maketitle
\begin{center}
Version: \version \quad \small{Compiled on \date{\today} at \currenttime}
\end{center}
\pagebreak
\blankpage
%--------------------------------------------------------------------------
%

%
%
%----------------------------
%---------COVER
%----------------------------
%
%\clearpage %for not having number page in the first page
%\maketitle
%\thispagestyle{empty} %for not having number page in the first page
%\pagebreak
%
%
%----------------------------
%---------PREFACE
%----------------------------
%
\section*{Preface}
This document is the final report of the project ``Ellipticity implementation in D3D'' (RWS \emph{bestelnummer} 4500268550, TUD \emph{kenmerk} 17363). 
\pagebreak
\clearpage 
\blankpage
% 
%
%----------------------------
%---------ABSTRACT
%----------------------------
%
\begin{abstract}
The currently applied model equations for predicting mixed-sediment river morphodynamics at large spatial scales and long time scales may be ill-posed under certain circumstances. Under these conditions it produces unphysical and unrealistic results. These conditions, in which the model loses its predictive capabilities, have been thoroughly studied under the assumption of one-dimensional flow. In this project we extend the analysis to two-dimensional conditions. We find that the consideration of secondary flow, a necessary mechanism to account for the intrinsic three-dimensionality of the flow in bends, may cause the model to be ill-posed if the diffusion coefficient of the transport equation of the secondary flow is below a certain threshold. Modeling a laboratory experiment we find that the diffusion coefficient necessary to obtain a well-posed model may be larger than a physically realistic value. The sediment transport rate depends primarily on the bed shear stress. In two-dimensional simulations it is crucial to take into consideration the effect of the bed slope (i.e., the downward gravitational pull) on the direction of the sediment transport rate. This effect is accounted for using a closure relation. We find that the consideration of the effect of the bed slope in the direction of the sediment transport rate is not only a physical mechanism necessary to reproduce natural situations but also a key process to obtain a well-posed model and that not all closure relations yield a well-posed model under mixed-size sediment conditions. We implement a routine in the software package Delft3D to check whether simulations suffers from ill-posedness caused by the mixed-size sediment model. This routine is used in a field scale simulation (DVR) and we find that it is ill-posed. In this case it is up to the modeler to decide whether the consequences are significant and whether the solution is reasonable and usable to address the problem for which the model was set up.
\end{abstract}
\pagebreak
\blankpage
%
%
%----------------------------
%---------TABLE OF CONTENTS
%----------------------------
%
\tableofcontents
\pagebreak
\pagenumbering{arabic}
\listoffigures
\pagebreak
\listoftables
\pagebreak
%
%----------------------------
%
\section{Introduction}
%
%----------------------------
%
\label{sec:int}
%
The mixed-size character of the sediment in alluvial rivers is a property necessary to explain physical phenomena such as downstream fining \citep{Sternberg75,Blom16}, the gravel sand transition zone \citep{Yatsu55,Blom17}, the formation of bedload sheets \citep{Seminara96}, the evolution of bars \citep{Lanzoni99}, and bed surface armoring \citep{Parker82}. 

In modeling the morphodynamic change of a river, estuary or the coast, most often a set of equations representing flow is solved in combination with a set of equations accounting for mass conservation of bed sediment. The currently applied model for predicting mixed-size sediment river morphodynamics at large spatial scales and long time scales is the active layer model \citep{Hirano71}. This is a deterministic model where the bed is discretized into layers and the sediment mixture in a finite number of size fractions. The topmost layer (the active layer) interacts with the flow (i.e., the sediment transport and friction depend on the conditions of this layer) while the layers below (the substrate) only act as bookkeeping system. There is a sediment flux from and to the active layer if there is a time variation of the elevation of the interface between the active layer and the substrate. This model has proven its value over a wide range of different situations for decades. 

Yet, under some circumstances the model may be ill-posed \citep{Ribberink87_PhD,Stecca14,Chavarrias18}. An ill-posed model produces unphysical and unrealistic results \citep{Joseph90}. The numerical solution of an ill-posed problem is characterized by the exponential growth of spurious oscillations dependent on the grid size. Thus, the solution does not converge with the grid, meaning that a refining of the grid causes a larger divergence of the solution. In these conditions the model loses its predictive capabilities.

There are alternatives to the active layer model such as the one by \citet{Ribberink87_PhD} where a layer below the active layer accounts for the mixing effects of large bedforms. This model was specifically derived for dune-dominated cases and it does not solve the problem of ill-posedness \citep{Sieben97_PhD}. The continuous model by \citet{Blom06_2,Blom08_2} was also derived for dune dominated cases and it is not applicable at large space scales or long time scales \citep{Blom08}. A simplified continuous model derived by \citet{Viparelli17} is applicable at large scales but it can also become ill-posed \citep{Chavarrias18}. Thus, the active layer model is still the main model for predicting mixed sediment river morphodynamics.

In a one-dimensional approach, the flow is represented by the \citet{SaintVenant71} equations, or in its steady form, the backwater equation. The conditions under which the one-dimensional system of equations can be ill-posed have been studied by \citet{Ribberink87_PhD}, \citet{Sieben97_PhD}, \citet{Stecca14}, and \citet{Chavarrias18}. When flow curvature plays a significant role in the morphodynamics the one-dimensional approach is no longer valid. This fact is clearly seen in river bends where the curvature of the flow creates a shallower inner bend and deeper outer bend as was first observed and explained by \citet{Thomson76}. In these cases the flow is intrinsically three-dimensional and is characterized by a mean flow and a secondary or spiral flow \citep{vanBendegom47}. The secondary flow can be of the same order of magnitude as the mean flow \citep{Flokstra77}. For large temporal and spatial scales as typically occurs in engineering applications, it is not feasible to solve a three-dimensional system of equations. In that case, the secondary flow can be parametrized and included in a two-dimensional system of equations (i.e., the Shallow Water Equations), considerably reducing the computational time \citep{Kalkwijk80,deVriend81_PhD,Kalkwijk86}. The parametrization is based on the intensity of the secondary flow $I$, which is a measure of the magnitude of the velocity component normal to the depth-averaged velocity. An advection-diffusion equation models the transport of this magnitude \citep{Jagers03_PhD}. Moreover, secondary flow introduces extra terms in the momentum equations \citep{Rodi82}.

A mechanism that counteracts secondary flow is the effect of the bed slope on the sediment transport direction. The direction of the sediment transport rate is mainly dependent on the direction of the bed shear stress. Yet, the bed slope has an effect, tending to move the sediment downslope \citep{vanBendegom47}. This mechanism is crucial to reproduce phenomena such as bar formation, braided rivers, and meandering \citep{Engelund73}.

\citet{Chavarrias16} studied the mathematical character of the system of equations formed by the Shallow Water Equations in combination with the active layer model. However, they did not study the effect of considering secondary flow and bed slope effects. Due to their simplifications the analysis was conducted studying the existence of the so-called Monge cones \citep{Sloff92,Sloff93,sieben94, Sieben97_PhD}. These cones represent the wave front of the propagation of disturbance. If they exist, the model is well-posed and vice versa. \citet{Chavarrias16} implemented a routine in the software package Delft3D to obtain the parameters needed to assess the well-posedness of a simulation. Yet, the actual check of the parameters was not implemented. Moreover, they did not study the implications of ill-posedness in 2D numerical simulations. 

The objective of this project is to extend the previous analysis by including secondary flow and bed slope effects. The report is organized as follows. In Section \ref{sec:rq} we summarize the research questions. The model equations are shown in Section \ref{sec:me}. In Section \ref{sec:pa} we linearize the system of equations. In Section \ref{sec:res_lin} we present the results of the linear analysis. In Section \ref{sec:ch_d3d} we discuss the implementation of the check in Delft3D. We apply the results of the linear analysis and the check in Delft3D to numerical simulations in Section \ref{sec:apl}. In Section \ref{sec:disc} we discuss our results. In Section \ref{sec:conc} we present the main conclusions, and in Section \ref{sec:rec} we summarize our recommendations.
%
%----------------------------
%
\section{Research Questions and Methodology}
%
%----------------------------
%
\label{sec:rq}
%
Our objective is to study the mathematical character of the system of equations used to model mixed-sediment morphodynamics in two dimensions. We will focus on the following research questions:
\begin{enumerate}
\item What is the role of secondary flow and bed slope effects as regards to the mathematical character of the system of equations?
\item What are the consequences of ill-posedness in 2D numerical simulations?
\item What are the implications of ill-posedness in field cases?
\end{enumerate}

To answer Research Question 1 we will conduct a mathematical analysis of the system of equations composed of the Shallow Water Equations in combination with the active layer model and a model to account for secondary flow and bed slope effects. We will assess the role of diffusive processes, as well as the relative importance of different physical mechanisms.

To answer Research Question 2 we will finish the implementation of the check tool in Delft3D. Then, the tool will be extended accounting for secondary flow. Once the tool is validated, it will be applied to idealized simulations to get insight on the consequences of ill-posedness. 

To answer Research Question 3 we will run a simulation which, rather than idealized, resembles the conditions frequently found in field cases. 
%
%----------------------------
%
\section{Model Equations}
%
%----------------------------
%
\label{sec:me}
%
In this section we describe the equations that model shallow water flow with a parametrized secondary flow correction together with the active layer model of \citet{Hirano71}. These equations represent two-dimensional hydrostatic flow over a mobile bed composed of several ($N$) non-cohesive sediment fractions. In Section \ref{subsec:sec_eq} we write the model equations. In Section \ref{subsec:sec_sim} we simplify the equations. In Section \ref{subsec:sec_cr} we describe the closure relations of the model. The equations are expanded in Section \ref{subsec:sec_ed} to write them in matrix form (Section \ref{subsec:sec_mf}). 
%
%----------------------------
%
\subsection{Balance Equations}
%
%----------------------------
%
\label{subsec:sec_eq}
%
In this section we present the model equations. 

The starting point of the flow equations are the Reynolds-averaged Navier-Stokes equations. We assume that the density variations are small (i.e, Boussinesq approximation) and we that the Reynolds stresses are dependent on the gradient of the mean flow velocity (i.e., eddy viscosity concept, Boussinesq closure relation). We assume shallow water flow, that is, that the vertical length scale and velocity is negligible respect to the horizontal length scale and velocity. 

Then, the equations are:
\begin{itemize}
\item Water mass conservation:
	\begin{equation}
	\label{eq:wmass_sf}
		\frac{\partial h}{\partial t}+\frac{\partial q_x}{\partial x}+\frac{\partial q_y}{\partial y}=0
	\end{equation}
\item Water momentum conservation in $x$ direction:
	\begin{equation}
	\label{eq:wmomx_sf}
	\begin{split}
		\frac{\partial q_x}{\partial t}+\frac{\partial (q_x^2/h+gh^2/2)}{\partial x}+\frac{\partial \left(\frac{q_xq_y}{h}\right)}{\partial y}+gh\frac{\partial \eta}{\partial x}-F^\prime_{sx} = \\
		=2\frac{\partial}{\partial x}\left(\nu_{t,H} h\frac{\partial (\frac{q_x}{h})}{\partial x}\right)+\frac{\partial}{\partial y}\left(\nu_{t,H} h \left(\frac{\partial (\frac{q_x}{h})}{\partial y} + \frac{\partial (\frac{q_y}{h})}{\partial x}\right)\right) -gh S_{fx}
		\end{split}
	\end{equation}
\item Water momentum conservation in $y$ direction:
	\begin{equation}
	\label{eq:wmomy_sf}
	\begin{split}
		\frac{\partial q_y}{\partial t}+\frac{\partial (q_y^2/h+gh^2/2)}{\partial y}+\frac{\partial \left(\frac{q_xq_y}{h}\right)}{\partial x}+gh\frac{\partial \eta}{\partial y}-F^\prime_{sy} =  \\
		2\frac{\partial}{\partial y}\left(\nu_{t,H} h\frac{\partial (\frac{q_y}{h})}{\partial y}\right)+\frac{\partial}{\partial x}\left(\nu_{t,H} h \left(\frac{\partial (\frac{q_y}{h})}{\partial x} + \frac{\partial (\frac{q_x}{h})}{\partial y} \right)\right) -gh S_{fy}
		\end{split}
	\end{equation}
\item Constitutive equation for the secondary flow intensity:
	\begin{equation}
	\label{eq:consI}
		\frac{\partial I}{\partial t}+\frac{q_x}{h}\frac{\partial I}{\partial x}+\frac{q_y}{h}\frac{\partial I}{\partial y}-D_H\frac{\partial^2 I}{\partial x^2}-D_H\frac{\partial^2 I}{\partial y^2}=S_s
	\end{equation}
\item Sediment mass conservation for the entire mixture (\citeauthor{Exner20}):
	\begin{equation}
	\label{eq:exner_sf}
		\frac{\partial \eta}{\partial t}+\frac{\partial q_{bx}}{\partial x}+\frac{\partial q_{by}}{\partial y}=0
	\end{equation}
\item Sediment mass conservation per grain size in the active layer (\citeauthor{Hirano71}):
	\begin{equation}
	\label{eq:hirano_full_sf}
		\frac{\partial M_{ak}}{\partial t}+f^{I}_{k}\frac{\partial (\eta-L_{a})}{\partial t}+\frac{\partial q_{bkx}}{\partial x}+\frac{\partial q_{bky}}{\partial y}=0 \quad k\in\{1,N-1\} 
	\end{equation}
\item Sediment mass conservation per grain size in the substrate:
	\begin{equation}
	\label{eq:substrate_full_sf}
		\frac{\partial M_{sk}}{\partial t}-f^{I}_{k}\frac{\partial (\eta-L_{a})}{\partial t}=0 \quad k\in\{1,N-1\} 
	\end{equation}
\end{itemize}
%
where:
\begin{itemize}
	\item $x$ 					= Cartesian x coordinate [$m$]
	\item $y$ 					= Cartesian y coordinate [$m$]
	\item $t$ 					= time coordinate [$s$]
	\item $h$ 					= flow depth [$m$]
	%\item $u$ 					= flow velocity in $x$ direction [$m/s$]
	%\item $v$ 					= flow velocity in $y$ direction [$m/s$]
	\item $q_x$ 					= specific water discharge in $x$ direction [$m^2/s$]
	\item $q_y$ 					= specific water discharge in $y$ direction [$m^2/s$]
	\item $\eta$  			= bed elevation [$m$]
	\item $f^{I}_{k}$ 			= volume fraction content of size fraction $k$ at the interface between the active layer and the substrate [$-$]
	\item $L_{a}$				= active layer thickness [$m$]
	%\item $\hat{q_{bkx}}$= sediment transport of size fraction $k$ in $x$ direction [${m^2}/s$]
	%\item $\hat{q_{bky}}$= sediment transport of size fraction $k$ in $y$ direction [${m^2}/s$]
	%\item $\hat{Q_{bkx}}$= sediment transport capacity of size fraction $k$ in $x$ direction [${m^2}/s$]
	%\item $\hat{Q_{bky}}$= sediment transport capacity of size fraction $k$ in $y$ direction [${m^2}/s$]
	\item $q_{bkx}$			= sediment transport of size fraction $k$ (including pores) in $x$ direction [$m^2/s$]
	\item $q_{bky}$			= sediment transport of size fraction $k$ (including pores) in $y$ direction [$m^2/s$]
	\item $q_{bx}$			= total sediment transport (including pores) in $x$ direction [$m^2/s$]
	\item $q_{by}$			= total sediment transport (including pores) in $y$ direction [$m^2/s$]
	%\item $Q_{bkx}$			= sediment transport capacity of size fraction $k$ (including pores) in $x$ direction [${m^2}/s$]
	%\item $Q_{bky}$			= sediment transport capacity of size fraction $k$ (including pores) in $y$ direction [${m^2}/s$]
	%\item $d_{sk}$  		= characteristic grain size of size fraction $k$ [$m$]
	%\item $D_{m}$	  		= mean grain size [$m$]
	\item $S_{fx}$   		= friction slope in the $x$ direction [$-$]
	\item $S_{fy}$   		= friction slope in the $y$ direction [$-$]
	%\item $p$ 					= porosity [$-$]
	\item $g$ 					= acceleration due to gravity [$m/s^2$]
	%\item $Fr$ 					= Froude number [$-$]
	%\item $C_f$					= dimensionless friction coefficient [$-$]
	\item $M_{ak}=F_{ak}L_{a}$ 			= volume of sediment of size fraction $k$ in the active layer per unit of surface area [$m$]
	\item $F_{ak}$ 			= volume fraction content of size fraction $k$ in the active layer [$-$]
	\item $M_{sk}=\int_{\eta_0}^{\eta-L_{a}}f_{sk}(z)\mathrm{d}z$ = volume of sediment of size fraction $k$ in the substrate per unit of surface area [$m$]
	\item $\eta_{0}$		= reference bed elevation [$m$]
	\item $f_{sk}$ 			= volume fraction of size fraction $k$ in the substrate [$-$]
	\item $N$ = number of size fractions [$-$]
  %	
	\item $F^\prime_{sx}$ = force per unit mass along the flow depth due to the secondary flow in the $x$ direction [$m^2/s^2$]
	\item $F^\prime_{sy}$ = force per unit mass along the flow depth due to the secondary flow in the $y$ direction [$m^2/s^2$]
	\item $\nu_{t,H}$ = horizontal eddy viscosity [$m^2/s$]
	\item $I$ = secondary flow intensity [$m/s$]
	\item $D_H$ = horizontal diffusion coefficient [$m^2/s$]
	\item $S_s$ = source term due to secondary flow [$m/s$]
\end{itemize}

%We use a notation with prime to avoid confusion with the Deltares notation which is based on the non-conservative Shallow Water Equations. 
% Delft3D manual F [m/s^2] force per unit mass

The volume fraction content of sediment at the active layer and the substrate are constrained by the equations:
%  \begin{linenomath*}
\begin{equation}
\label{eq:Fcons}
	 \sum_{k=1}^N F_{ak}=1,\quad \sum_{k=1}^N f_{sk}(z)=1 \;,
\end{equation}
%  \end{linenomath*}
thus, the volume of sediment per unit area are constrained by the equations:
\begin{equation}
\label{eq:Mcons}
	 \sum_{k=1}^N M_{ak}=L_a, \quad \sum_{k=1}^N M_{sk}=\eta -L_a-\eta_0 \;.
\end{equation}
%  \end{linenomath*}
%
%----------------------------
\subsection{Simplifications}
%----------------------------
%
%
\label{subsec:sec_sim}
%
In this section we simplify the system of equations to ease the mathematical analysis. 

As regards to the flow we assume the eddy viscosity to be zero. This implies that there are no diffusive processes in the flow equations. The momentum equations are in this case:
	\begin{equation}
	\label{eq:wmomx_sf_s}
		\frac{\partial q_x}{\partial t}+\frac{\partial (q_x^2/h+gh^2/2)}{\partial x}+\frac{\partial \left(\frac{q_xq_y}{h}\right)}{\partial y}+gh\frac{\partial \eta}{\partial x}-F^\prime_{sx} =-gh S_{fx} \;,
	\end{equation}
and 
	\begin{equation}
	\label{eq:wmomy_sf_s}
		\frac{\partial q_y}{\partial t}+\frac{\partial (q_y^2/h+gh^2/2)}{\partial y}+\frac{\partial \left(\frac{q_xq_y}{h}\right)}{\partial x}+gh\frac{\partial \eta}{\partial y}-F^\prime_{sy} = -gh S_{fy} \;.
	\end{equation}
	
We assume the active layer thickness to be constant:
\begin{equation}
\label{eq:cntLa}
	\frac{\partial L_a}{\partial t}=0 \;.
\end{equation}
Substitution of Equation \ref{eq:cntLa} in \ref{eq:hirano_full_sf} and \ref{eq:substrate_full_sf} yields:
\begin{equation}
\label{eq:hirano}
	\frac{\partial M_{ak}}{\partial t}-f^{I}_{k}\frac{\partial q_{bx}}{\partial x}-f^{I}_{k}\frac{\partial q_{by}}{\partial y}+\frac{\partial q_{bkx}}{\partial x}+\frac{\partial q_{bky}}{\partial y}=0 \quad k\in\{1,N-1\}
\end{equation}
\begin{equation}
\label{eq:substrate}
	\frac{\partial M_{sk}}{\partial t}+f^{I}_{k}\frac{\partial q_{bx}}{\partial x}+f^{I}_{k}\frac{\partial q_{by}}{\partial y}=0 \quad k\in\{1,N-1\}
\end{equation}
%
%
%----------------------------
\subsection{Closure Relations}
%----------------------------
%
%
\label{subsec:sec_cr}
%
The governing equations still need closure relations for the friction slope, the sediment transport rate, the secondary flow terms, and the fractions at the interface between the substrate and the active layer to form a complete set of equations. In this section we describe those closure relations.

\subsubsection{Friction Slope}

The friction slope is:
\begin{equation}
\label{eq:Sf}
	S_{fx}=\frac{C_fq_xQ}{gh^3} \quad S_{fy}=\frac{C_fq_yQ}{gh^3}
\end{equation}
%
where:
\begin{itemize}
	\item $C_f$					= dimensionless friction coefficient [$-$]
	\item $Q=|\vec{q}|$ = module of the specific water discharge [$m/s^2$]
\end{itemize}

\subsubsection{Sediment Transport Rate}
\label{subsubsec:sedtrans}

The sediment transport rate per size fraction (including pores) $\vec{q_{bk}}$ [$m^2/s$] can be expressed as:
\begin{equation}
\label{eq:qbk_sf}
	\vec{q_{bk}}=(q_{bkx},q_{bky})=q_{bk}(\cos{\varphi_{sk}},\sin{\varphi_{sk}}) \quad k\in\{1,N\}
\end{equation}
where:
\begin{itemize}
	\item $\varphi_{sk}$ = direction of the sediment transport rate (correcting for bed slope and secondary flow) [rad]
	\item $q_{bk}$					= absolute value of the sediment transport rate including pores [$m^2/s$]
\end{itemize}

The direction of the sediment transport rate $\varphi_{sk}$ [rad] is:
\begin{equation}
\label{eq:phis}
	\tan{\varphi_{sk}}=\frac{\sin{\varphi_{\tau}}-\frac{1}{g_{sk}}\frac{\partial \eta}{\partial y}}{\cos{\varphi_{\tau}}-\frac{1}{g_{sk}}\frac{\partial \eta}{\partial x}} \quad k\in\{1,N\}
\end{equation}
where:
\begin{itemize}
	\item $\varphi_{\tau}$ = direction of the sediment transport rate only correcting for secondary flow [$rad$]
	\item $g_{sk}$ = bed slope function [$-$]
\end{itemize}

The direction of the sediment transport rate $\varphi_{\tau}$ [rad] is:
\begin{equation}
\label{eq:phitau0}
	\tan{\varphi_{\tau}}=\frac{q_y-h\alpha_I\frac{q_x}{Q}I}{q_x-h\alpha_I\frac{q_y}{Q}I}
\end{equation}
%
where:
\begin{itemize}
	\item $\alpha_I$ = constant [$-$]
\end{itemize}

To obtain the bed shear stress deviation due to secondary flow it is necessary to estimate the radial velocity. To this end under uniform flow conditions are considered, the advection term in the momentum equations is neglected, a mixing length is assumed, the radial and vertical velocities are assumed to be small relative to the streamwise velocity, and a logarithmic velocity profile is considered \citep{Rozovskii57}. In these conditions the constant $\alpha_I$ [$-$] is:
\begin{equation}
\label{eq:alphaI}
	\alpha_I=\frac{2}{\kappa^2}E_s\left(1-\frac{\sqrt{C_f}}{\kappa}\right)
\end{equation}
where:
\begin{itemize}
	\item $\kappa$ = Von K\'arm\'an constant [$-$] 
	\item $E_s$ = calibration parameter [$-$]
\end{itemize}

The bed slope function $g_{sk}$ is:
\begin{equation}
\label{eq:g_sk}
	g_{sk}=A_s\theta_k^{B_s}\left(\frac{d_k}{h}\right)^{C_s}\left(\frac{d_k}{D_m}\right)^{D_s}
\end{equation}
where:
\begin{itemize}
	\item $A_s$ = calibration constant for gravitational pull on sediment transport [$-$]
	\item $B_s$ = calibration constant for gravitational pull on sediment transport [$-$]
	\item $C_s$ = calibration constant for gravitational pull on sediment transport [$-$]
	\item $D_s$ = calibration constant for gravitational pull on sediment transport [$-$]
	\item $\theta_k$ = module of the Shields stress of size fraction $k$ [$-$]
	\item $d_k$ = characteristic grain size of size fraction $k$ [$m$]
	\item $D_m$ = characteristic mean grain size of the mixture [$m$]
\end{itemize}
Values for the constants are summarized in Table \ref{tab:gsval}.
\begin{table}[ht]
	\begin{center}
	\begin{tabular}{rrrrr}		
			\toprule
			$A_s$ & $B_s$ & $C_s$ & $D_s$ & Note\\
			\midrule
	1.70	& 0.5	& 0.0 & 0.0 & \citet{Talmon95} (their experiments, straight flume)  \\ 
	0.85	& 0.5 & 0.0 & 0.0 & \citet{Talmon95} (experiments by \citet{Zimmerman78}, circular flume)  \\ 
	9.00	& 0.5 & 0.3 & 0.0 & \citet{Talmon95} (2D runs of field and lab)\\ 
			\bottomrule
			\end{tabular}
      \caption{Parameters of the function for the bed slope.}
      \label{tab:gsval}
	\end{center}
\end{table}

The absolute value of the sediment transport rate is:
%  \begin{linenomath*}
\begin{equation}
\label{eq:qbk}
	q_{bk}=F_{ak}\sqrt{g R d_k^3}\left(1-p\right)q^*_{bk} \quad k\in\{1,N\} \;,
\end{equation}
%  \end{linenomath*}
where:
\begin{itemize}
	\item $q^*_{bk}$ = nondimensional sediment transport rate [$-$]
	\item $p$ 					= porosity [$-$]
	\item $R=\rho_s/\rho_w-1$ =  submerged sediment density [$-$]
	\item $\rho_s=2650$  = sediment density [$kg/m^3$]
	\item $\rho_w=1000$  = water density [$kg/m^3$]
\end{itemize}

The nondimensional sediment transport rate is computed using, for instance, a generalized form of the relation developed by \citet{MeyerPeter48}
%  \begin{linenomath*}
\begin{equation}
\label{eq:qbk_st}
	q^*_{bk}=A\max\left(\theta_k-\xi_k \theta_c,0\right)^B \quad k\in\{1,N\} \;,
\end{equation}
%  \end{linenomath*}
where:
\begin{itemize}
	\item $A$ = nondimensional parameter [$-$]
	\item $B$ = nondimensional parameter [$-$]
	\item $\theta_c$ = nondimensional critical bed shear stress [$-$]
	\item $\xi_k$ = hiding coefficient [$-$]
\end{itemize}

The module of the Shields stress is:
\begin{equation}
\label{eq:thetak}
	\theta_k=\frac{C_f\left(\frac{Q}{h}\right)^2}{gRd_k} \quad k\in\{1,N\}
\end{equation}

A common hiding functions is the one due to \citet{Egiazaroff65}:
%  \begin{linenomath*}
\begin{equation}
\label{eq:egiazaroff}
	\xi_{k}=\left(\frac{\log_{10}\left(19\right)}{\log_{10}\left(19\frac{d_k}{D_m}\right)}\right)^2 \quad k\in\{1,N\}
\end{equation}
%  \end{linenomath*}
%
%where:
%\begin{itemize}
	%
%\end{itemize}

A simpler expression was developed by \citet{Parker82_3}:
%  \begin{linenomath*}
\begin{equation}
\label{eq:xi_parker}
	\xi_{k}=\left(\frac{D_m}{d_k}\right)^b \quad k\in\{1,N\}
\end{equation}
%  \end{linenomath*}
%
where:
\begin{itemize}
	\item $b$ = nondimensional parameter [$-$]
\end{itemize}

We define the sediment transport capacity $Q_{bk}$ as:
\begin{equation}
\label{eq:Qbk}
	Q_{bk}=\sqrt{g R d_k^3}\left(1-p\right)q^*_{bk} \quad k\in\{1,N\} \;,
\end{equation}
which is the sediment transport rate as if all the sediment at the bed surface was composed by the same size fraction, yet retaining hiding effects. This term is of mathematical convenience but has no true physical interpretation. In general terms, $Q_{bk}>Q_{bk+1}$. Reverse mobility conditions occur when hiding so strong that $b>1$ implying $Q_{bk}<Q_{bk+1}$.

%-------------
\subsubsection{Secondary Flow Terms}
\label{subsubsec:sf_cr}

The secondary flow terms in the momentum equations are:
\begin{equation}
\label{eq:Fsx}
	F^\prime_{sx}=\frac{\partial T^\prime_{xx}}{\partial x}+\frac{\partial T^\prime_{xy}}{\partial y}
\end{equation}
\begin{equation}
\label{eq:Fsy}
	F^\prime_{sy}=\frac{\partial T^\prime_{yx}}{\partial x}+\frac{\partial T^\prime_{yy}}{\partial y}
\end{equation}
%
where:
\begin{itemize}
	\item $T^\prime_{lm}$ = shear stress per unit mass and volume along the flow depth in the direction $l$-$m$ [$m^3/s^2$] 
\end{itemize}
% Delft3D manual T [m^2/s^2] shear stress per unit density

The closure relation for the secondary flow force terms are:
\begin{equation}
\label{eq:Txx}
	T^\prime_{xx}=-2\frac{\beta^*I}{Q}q_x q_y
\end{equation}
\begin{equation}
\label{eq:Txy}
	T^\prime_{xy}=T^\prime_{yx}=\frac{\beta^*I}{Q}\left(q_x^2-q_y^2\right)
\end{equation}
\begin{equation}
\label{eq:Tyy}
	T^\prime_{yy}=T^\prime_{yy}=2\frac{\beta^*I}{Q}q_x q_y
\end{equation}
%
where:
\begin{itemize}
	\item $\beta^*=\beta_c\left(5\alpha-15.6\alpha^2+37.5\alpha^3\right)$ = constant [$-$]
	\item $\beta_c$ $\in[0,1]$ = calibration parameter [$-$]
	\item $\alpha=\frac{\sqrt{C_f}}{\kappa}<0.5$ = constant [$-$]
\end{itemize}

The closure relation for the source term of the secondary flow equation (\ref{eq:consI}) is:
\begin{equation}
\label{eq:Ss}
	Ss=-\frac{I-I_e}{T_I}
\end{equation}
%
where:
\begin{itemize}
	\item $I_e=Q/R_s$ = equilibrium secondary flow intensity [$m/s$]
	\item $T_I$  = adaptation time scale [$s$]
	\item $R_s$  = radius of curvature of the streamlines [$m$]
\end{itemize}

The radius of curvature of the streamlines is defined as \citep[e.g.][]{Legleiter06}:
\begin{equation}
\label{eq:Rs_1}
	\frac{1}{R_s}=\frac{\frac{\mathrm{d}x}{\mathrm{d}t}\frac{\mathrm{d}^2y}{\mathrm{d}t^2}-\frac{\mathrm{d}y}{\mathrm{d}t}\frac{\mathrm{d}^2x}{\mathrm{d}t^2}}{\left(\left(\frac{\mathrm{d}x}{\mathrm{d}t}\right)^2+\left(\frac{\mathrm{d}y}{\mathrm{d}t}\right)^2\right)^{3/2}} \;,
\end{equation}
substituting $u=\frac{\mathrm{d}x}{\mathrm{d}t}$ and $v=\frac{\mathrm{d}y}{\mathrm{d}t}$ we obtain:
\begin{equation}
\label{eq:Rs_2}
	\frac{1}{R_s}=\frac{u\frac{\mathrm{d}v}{\mathrm{d}t}-v\frac{\mathrm{d}u}{\mathrm{d}t}}{\left(u^2+v^2\right)^{3/2}} \;,
\end{equation}
expanding the material derivatives and assuming steady flow we obtain:
\begin{equation}
\label{eq:Rs_3}
	\frac{1}{R_s}=\frac{u^2\frac{\partial v}{\partial x}+uv\frac{\partial v}{\partial y}-uv\frac{\partial u}{\partial x}-v^2\frac{\partial u}{\partial y}}{\left(u^2+v^2\right)^{3/2}} \;,
\end{equation}
in terms of water discharge we obtain:
\begin{equation}
\label{eq:Rs_4}
	\frac{1}{R_s}=\frac{-q_xq_y\frac{\partial q_x}{\partial x}+q_x^2\frac{\partial q_y}{\partial x}-q_y^2\frac{\partial q_x}{\partial y}+q_xq_y\frac{\partial q_y}{\partial y}}{\left(q_x^2+q_y^2\right)^{3/2}} \;.
\end{equation}

The adaptation time scale $T_I$ is:
\begin{equation}
\label{eq:TI}
	T_I=\frac{L_I h}{Q}
\end{equation}
%
where:
\begin{itemize}
	\item $L_I=L^*_Ih$ = adaptation length scale [$m$]
	\item $L^*_I=\frac{1-2\alpha}{2\kappa^2\alpha}$ = nondimensional adaptation length scale [$-$]
\end{itemize}

The diffusion coefficient $D_H$ accounts for the filtered processes, i.e., those physical processes not resolved by the equations. In the derivation of the Shallow Water Equations, we have filtered \citep[e.g.][]{Madsen88}: (1) the random molecular motion, (2) the turbulent motion below a certain scale, and (3) the vertical velocity profile. In general terms, the third filtering process is the most relevant one \citep{Fischer73,Struiksma85_2,Madsen88}. It is important to notice that, given the anisotropy of the flow field the diffusion tensor has in principle non-diagonal terms which are of importance close to corners \citep{Fischer73}. Far from corners the diagonal terms (i.e, streamwise direction $D_{Hs}$ and transverse direction $D_{Hn}$) dominate. \citet{Elder59} derives a diffusion coefficient assuming a logarithmic profile for the vertical velocity:
\begin{equation}
\label{eq:DHs}
	D_{Hs}=\left(\frac{0.4041}{\kappa^3}+\frac{1}{6}\kappa\right)hu^* \;,
\end{equation}
%
\begin{equation}
\label{eq:DHn}
	D_{Hn}=\frac{1}{6}\kappa hu^* \;.
\end{equation}
These values of the horizontal diffusion coefficient apply if the spatial grid is of the order of magnitude of the flow depth. \citeauthor{Elder59} neglects the effect of the viscous sublayer, which causes his analytical expression to be a lower limit of the diffusion coefficient \citet{Fischer67}. Yet, it provides a reasonable order of magnitude. 

Transverse mixing due to secondary flow contribute to a reduction of the longitudinal diffusion \citet{Erdogan67}. \citet{Fischer69} derives a diffusion coefficient which accounts for secondary flow:
%
\begin{equation}
\label{eq:DHn_sf}
	D_{Hn}=0.25\kappa \frac{u^2h^3}{\kappa^5R_s^2u^*} \;.
\end{equation}

%However, in general terms it is used as a calibration coefficient and uniform value is used. 

\subsubsection{Volume Fraction Content at the Interface}

The volume fraction content at the interface between the active layer and the substrate under degradational conditions is assumed to be equal to the volume fraction content at the top part of the substrate. Under aggradational conditions \citet{Hirano71} proposed the flux to the substrate to have the same grain size distribution as the active layer. \citet{Parker91_1} introduced the concept that the aggradational flux to the substrate is also influenced by the grain size distribution of the bed load. \citet{Hoey94} combined both concepts in a parameter that sets the contribution of the bed load relative to the active layer. Currently, only the initial concept of Hirano is implemented in Delft3D which in mathematical terms can be written as:
%  \begin{linenomath*}
\begin{equation}
\label{eq:hoey}
f^{I}_{k}=
  \begin{cases}
    f_{sk}(z=\eta-L_a) & \text{if} \quad \frac{\partial \left(\eta-L_a\right)}{\partial t}<0 \\
    F_{ak}						 & \text{if} \quad \frac{\partial \left(\eta-L_a\right)}{\partial t}>0 \\
  \end{cases}
\end{equation}

%
%
%----------------------------
\subsection{Expansion of the System Equations}
%----------------------------
%
%
\label{subsec:sec_ed}
%
All the terms in the model equations need to be expressed as a function of the second derivative, first derivative, linear relation, or source of the dependent variables of the system (i.e., $h,q_x,q_y,I,\eta,M_{ak} \forall k\in \{1,N-1\}$). The terms that need to be decomposed are the secondary flow terms in the momentum equations (Section \ref{subsubsec:sec_f_mom}), the sediment transport rate (Section \ref{subsubsec:sed_tra}), and the source term in the constitutive equation for secondary flow (Section \ref{subsubsec:stsf}). The system equations are written in expanded form in Sections \ref{subsubsec:exp_wmass}-\ref{subsubsec:exp_hirano_sf}.

\subsubsection{Expansion of the Secondary Flow Terms in the Momentum Equations}
\label{subsubsec:sec_f_mom}

The secondary flow terms in the momentum equations, Equations (\ref{eq:Txx}), (\ref{eq:Txy}), and (\ref{eq:Txy}) are a function of the specific water discharge and secondary flow intensity only:
\begin{equation}
\label{eq:Tiifuncx_sf}
\begin{split}
	T^\prime_{ij}=f(q_x,q_y,I) \Rightarrow \\ 
	\Rightarrow \frac{\partial T^\prime_{ij}}{\partial x}=\frac{\partial T^\prime_{ij}}{\partial q_x}\frac{\partial q_x}{\partial x}+\frac{\partial T^\prime_{ij}}{\partial q_y}\frac{\partial q_y}{\partial x}+\frac{\partial T^\prime_{ij}}{\partial I}\frac{\partial I}{\partial x}\\
	\\
	\text{and}\\
	\\
	\Rightarrow \frac{\partial T^\prime_{ij}}{\partial y}=\frac{\partial T^\prime_{ij}}{\partial q_x}\frac{\partial q_x}{\partial y}+\frac{\partial T^\prime_{ij}}{\partial q_y}\frac{\partial q_y}{\partial y}+\frac{\partial T^\prime_{ij}}{\partial I}\frac{\partial I}{\partial y}\\
\end{split}
\end{equation} 
where $i=(x,y)$ and $j=(x,y)$.

\subsubsection{Expansion of the Sediment Transport Rate}
\label{subsubsec:sed_tra}

The directional sediment transport rates $q_{bkx}$ and $q_{bky}$ are not only a function of the flow velocity in their respective directions but also of the flow velocity in the other direction and the secondary flow intensity: 
\begin{equation}
\label{eq:qsfuncx_sf}
\begin{split}
	q_{bkx}=f(h,q_x,q_y,I,M_{ak},\frac{\partial \eta}{\partial x},\frac{\partial \eta}{\partial y}) \Rightarrow \\ 
	\Rightarrow \frac{\partial q_{bkx}}{\partial x}=\frac{\partial q_{bkx}}{\partial h}\frac{\partial h}{\partial x}+\frac{\partial q_{bkx}}{\partial q_x}\frac{\partial q_x}{\partial x}+\frac{\partial q_{bkx}}{\partial q_y}\frac{\partial q_y}{\partial x}+\frac{\partial q_{bkx}}{\partial I}\frac{\partial I}{\partial x}\\
	+\sum_{l=1}^{N-1}\frac{\partial q_{bkx}}{\partial M_{al}}\frac{\partial M_{al}}{\partial x}+\frac{\partial q_{bkx}}{\partial \frac{\partial \eta}{\partial x}}\frac{\partial^2 \eta}{\partial x^2}+\frac{\partial q_{bkx}}{\partial \frac{\partial \eta}{\partial y}}\frac{\partial^2 \eta}{\partial x \partial y}  \quad k\in\{1,N\} 
\end{split}
\end{equation}
and, 
\begin{equation}
\label{eq:qsfuncy_sf}
\begin{split}
	q_{bky}=f(h,q_x,q_v,I,M_{ak},\frac{\partial \eta}{\partial x},\frac{\partial \eta}{\partial y}) \Rightarrow \\ 
	\Rightarrow \frac{\partial q_{bky}}{\partial y}=\frac{\partial q_{bky}}{\partial h}\frac{\partial h}{\partial y}+\frac{\partial q_{bky}}{\partial q_x}\frac{\partial q_x}{\partial y}+\frac{\partial q_{bky}}{\partial q_y}\frac{\partial q_y}{\partial y}+\frac{\partial q_{bky}}{\partial I}\frac{\partial I}{\partial y}\\
	+\sum_{l=1}^{N-1}\frac{\partial q_{bky}}{\partial M_{al}}\frac{\partial M_{al}}{\partial y}+\frac{\partial q_{bky}}{\partial \frac{\partial \eta}{\partial x}}\frac{\partial^2 \eta}{\partial x \partial y}+\frac{\partial q_{bky}}{\partial \frac{\partial \eta}{\partial y}}\frac{\partial^2 \eta}{\partial y^2}  \quad k\in\{1,N\} 
\end{split}
\end{equation}

The same holds for the total bed load in $x$ and $y$ direction ($q_{bx}$ and $q_{by}$). Note that $M_{aN}$ is not an independent variable. 

%The derivatives of the directional sediment transport rates with respect to the dependent variables have the same expression as when secondary flow is not considered. However, in this case the direction of the sediment transport depends on all the variables but the volume of sediment in the active layer. Thus, some of the derivatives of the sediment transport direction that are zero when secondary flow is not considered are not 0 when this is considered.

\subsubsection{Expansion of the Source Term in the Constitutive Equation for the Secondary Flow}
\label{subsubsec:stsf}

Here we expand the term $S_s$ of Equation (\ref{eq:consI}).

\begin{equation}
\label{eq:Ssfunc_sf}
\begin{split}
	S_s=-\frac{I-I_e}{T_I}=\frac{-I}{T_I}+\frac{I_e}{T_I}=\frac{-QI}{h^2L^*_I}+S_{sxx}\frac{\partial q_x}{\partial x}+S_{syx}\frac{\partial q_y}{\partial x}+S_{sxy}\frac{\partial q_x}{\partial y}+S_{syy}\frac{\partial q_y}{\partial y}
\end{split}
\end{equation}
where:
\begin{equation}
\begin{split}
	S_{sxx}=\frac{1}{Qh^2L^*_I}(-q_xq_y) \\
	S_{syx}=\frac{1}{Qh^2L^*_I}(q_x^2) \\
	S_{sxy}=\frac{1}{Qh^2L^*_I}(-q_y^2) \\
	S_{syy}=\frac{1}{Qh^2L^*_I}(q_xq_y) \\
\end{split}
\end{equation}

\subsubsection{Water Mass Conservation}
\label{subsubsec:exp_wmass}

All the terms of Equation (\ref{eq:wmass_sf}) are already expanded:
\begin{equation}
\label{eq:wmass_sf_exp}
	\frac{\partial h}{\partial t}+\frac{\partial q_x}{\partial x}+\frac{\partial q_y}{\partial y}=0
\end{equation}

\subsubsection{Water Momentum Conservation in \emph{x} Direction}
\label{subsubsec:exp_wmomx}

\begin{equation}
\label{eq:wmomx_sf_exp}
\begin{split}
	\frac{\partial q_x}{\partial t}+\frac{\partial (q_x^2/h+gh^2/2)}{\partial x}+\frac{\partial \left(\frac{q_xq_y}{h}\right)}{\partial y}+gh\frac{\partial \eta}{\partial x}-F^\prime_{sx}=-gh S_{fx}\implies\\
	\frac{\partial q_x}{\partial t}+\left(gh-\left(\frac{q_x}{h}\right)^2\right)\frac{\partial h}{\partial x}+\left(2\frac{q_x}{h}-\frac{\partial T^\prime_{xx}}{\partial q_x}\right)\frac{\partial q_x}{\partial x}-\frac{\partial T^\prime_{xx}}{\partial q_y}\frac{\partial q_y}{\partial x}-\frac{\partial T^\prime_{xx}}{\partial I}\frac{\partial I}{\partial x}+gh\frac{\partial \eta}{\partial x}+\\
	+\frac{-q_x q_y}{h^2}\frac{\partial h}{\partial y}+\left(\frac{q_y}{h}-\frac{\partial T^\prime_{xy}}{\partial q_x}\right)\frac{\partial q_x}{\partial x}-\left(\frac{q_x}{h}-\frac{\partial T^\prime_{xy}}{\partial q_y}\right)\frac{\partial q_y}{\partial y}-\frac{\partial T^\prime_{xy}}{\partial I}\frac{\partial I}{\partial y}+\\
	+\frac{C_fq_xQ}{h^2}=0
\end{split}
\end{equation}

\subsubsection{Water Momentum Conservation in \emph{y} Direction}
\label{subsubsec:exp_wmomy}

\begin{equation}
\label{eq:wmomy_sf_exp}
\begin{split}
		\frac{\partial q_y}{\partial t}+\frac{\partial (q_y^2/h+gh^2/2)}{\partial y}+\frac{\partial \left(\frac{q_xq_y}{h}\right)}{\partial x}+gh\frac{\partial \eta}{\partial y}-F^\prime_{sy}=-gh S_{fy} \implies\\
	\frac{\partial q_y}{\partial t}+\frac{-q_{x}q_{y}}{h^2}\frac{\partial h}{\partial x}+\left(\frac{q_y}{h}-\frac{\partial T^\prime_{yx}}{\partial q_x}\right)\frac{\partial q_x}{\partial x}+\left(\frac{q_x}{h}-\frac{\partial T^\prime_{yx}}{\partial q_y}\right)\frac{\partial q_y}{\partial x}-\frac{\partial T^\prime_{yx}}{\partial I}\frac{\partial I}{\partial x}+\\
	+\left(gh-\left(\frac{q_y}{h}\right)^2\right)\frac{\partial h}{\partial y}-\frac{\partial T^\prime_{yy}}{\partial q_x}\frac{\partial q_x}{\partial y}+\left(2\frac{q_y}{h}-\frac{\partial T^\prime_{yy}}{\partial q_y}\right)\frac{\partial q_y}{\partial y}-\frac{\partial T^\prime_{yy}}{\partial I}\frac{\partial I}{\partial y}+gh\frac{\partial \eta}{\partial y}+\\
	+\frac{C_fq_yQ}{h^2}=0
\end{split}
\end{equation}

\subsubsection{Constitutive Equation for the Secondary Flow}
\label{subsubsec:exp_consI}

\begin{equation}
\label{eq:consI_exp}
\begin{split}
\frac{\partial I}{\partial t}+\frac{q_x}{h}\frac{\partial I}{\partial x}+\frac{q_y}{h}\frac{\partial I}{\partial y}-D_H\frac{\partial^2 I}{\partial x^2}-D_H\frac{\partial^2 I}{\partial y^2}=S_s
 \implies\\
\frac{\partial I}{\partial t}-D_H\frac{\partial^2 I}{\partial x^2}-D_H\frac{\partial^2 I}{\partial y^2}-S_{sxx}\frac{\partial q_x}{\partial x}-S_{syx}\frac{\partial q_y}{\partial x}+\frac{q_x}{h}\frac{\partial I}{\partial x}-S_{sxy}\frac{\partial q_x}{\partial y}-S_{syy}\frac{\partial q_y}{\partial y}+\frac{q_y}{h}\frac{\partial I}{\partial y}+\frac{QI}{h^2L^*_I}=0
\end{split}
\end{equation}

\subsubsection{Sediment Mass Conservation for the Entire Mixture}
\label{subsubsec:exp_exner_sf}

\begin{equation}
\label{eq:exner_sf_exp}
\begin{split}
		\frac{\partial \eta}{\partial t}+\frac{\partial q_{bx}}{\partial x}+\frac{\partial q_{by}}{\partial y}=0
 \implies\\
\frac{\partial \eta}{\partial t}+\frac{\partial q_{bx}}{\partial h}\frac{\partial h}{\partial x}+\frac{\partial q_{bx}}{\partial q_x}\frac{\partial q_x}{\partial x}+\frac{\partial q_{bx}}{\partial q_y}\frac{\partial q_y}{\partial x}+\frac{\partial q_{bx}}{\partial I}\frac{\partial I}{\partial x}+\sum_{l=1}^{N-1}\frac{\partial q_{bx}}{\partial M_{al}}\frac{\partial M_{al}}{\partial x}+\\
+\frac{\partial q_{by}}{\partial h}\frac{\partial h}{\partial y}\frac{\partial q_{by}}{\partial q_x}\frac{\partial q_x}{\partial y}+\frac{\partial q_{by}}{\partial q_y}\frac{\partial q_y}{\partial y}+\frac{\partial q_{by}}{\partial I}\frac{\partial I}{\partial y}+\sum_{l=1}^{N-1}\frac{\partial q_{by}}{\partial M_{al}}\frac{\partial M_{al}}{\partial y}
\end{split}
\end{equation}

\subsubsection{Sediment Mass Conservation per Grain Size in the Active Layer}
\label{subsubsec:exp_hirano_sf}

\begin{equation}
\label{eq:hirano_sf_exp}
\begin{split}
	\frac{\partial M_{ak}}{\partial t}-f^{I}_{k}\frac{\partial q_{bx}}{\partial x}-f^{I}_{k}\frac{\partial q_{by}}{\partial y}+\frac{\partial q_{bkx}}{\partial x}+\frac{\partial q_{bky}}{\partial y}=0 \quad k\in\{1,N-1\}
 \implies\\
\frac{\partial M_{ak}}{\partial t}+\left[\frac{\partial q_{bkx}}{\partial h}-f^{I}_{k}\frac{\partial q_{bx}}{\partial h}\right]\frac{\partial h}{\partial x}+\left[\frac{\partial q_{bkx}}{\partial q_x}-f^{I}_{k}\frac{\partial q_{bx}}{\partial q_x}\right]\frac{\partial q_x}{\partial x}+\left[\frac{\partial q_{bkx}}{\partial q_y}-f^{I}_{k}\frac{\partial q_{bx}}{\partial q_y}\right]\frac{\partial q_y}{\partial x}+\\
+\left[\frac{\partial q_{bkx}}{\partial I}-f^{I}_{k}\frac{\partial q_{bx}}{\partial I}\right]\frac{\partial I}{\partial x}+\sum_{l=1}^{N-1}\left[\frac{\partial q_{bkx}}{\partial M_{al}}-f^{I}_{k}\frac{\partial q_{bx}}{\partial M_{al}}\right]\frac{\partial M_{al}}{\partial x}+\\
\left[\frac{\partial q_{bkx}}{\partial h}-f^{I}_{k}\frac{\partial q_{bx}}{\partial h}\right]\frac{\partial h}{\partial y}+\left[\frac{\partial q_{bkx}}{\partial q_x}-f^{I}_{k}\frac{\partial q_{bx}}{\partial q_x}\right]\frac{\partial q_x}{\partial y}+\left[\frac{\partial q_{bkx}}{\partial q_y}-f^{I}_{k}\frac{\partial q_{bx}}{\partial q_y}\right]\frac{\partial q_y}{\partial y}+\\
+\left[\frac{\partial q_{bkx}}{\partial I}-f^{I}_{k}\frac{\partial q_{bx}}{\partial I}\right]\frac{\partial I}{\partial y}+\sum_{l=1}^{N-1}\left[\frac{\partial q_{bkx}}{\partial M_{al}}-f^{I}_{k}\frac{\partial q_{bx}}{\partial M_{al}}\right]\frac{\partial M_{al}}{\partial y}+\\
\end{split}
\end{equation}
%
%
%----------------------------
\subsection{Matrix Formulation}
%----------------------------
%
%
\label{subsec:sec_mf}
%
In this section we write the system of equations in matrix formulation. Equation (\ref{eq:substrate}) is a linear combination of equations (\ref{eq:hirano}) and (\ref{eq:exner_sf}). The rest of the equations do not depend on $M_{sk}$. Thus, the substrate equations provide a zero eigenvalue with multiplicity $N-1$. To simplify the writing we omit the substrate equations.

We recast the model Equations, (\ref{eq:wmass_sf_exp}), (\ref{eq:wmomx_sf_exp}), (\ref{eq:wmomy_sf_exp}), (\ref{eq:consI_exp}), (\ref{eq:exner_sf_exp}), and (\ref{eq:hirano_sf_exp}) in matrix form:
\begin{equation}
\label{eq:matrixf_sf}
	\frac{\partial \mathbf{Q}}{\partial t}+\mathbf{D_x}\frac{\partial^2 \mathbf{Q}}{\partial x^2}+\mathbf{D_y}\frac{\partial^2 \mathbf{Q}}{\partial y^2}+\mathbf{C}\frac{\partial^2 \mathbf{Q}}{\partial x \partial y}+\mathbf{A_x}\frac{\partial \mathbf{Q}}{\partial x}+\mathbf{A_y}\frac{\partial \mathbf{Q}}{\partial y}+\mathbf{S}=0
\end{equation}

The dependent variables are $h$, $q_x$, $q_y$, $I$, $\eta$, and $M_{ak}$ for $1\leq k \leq N-1$:
\begin{equation}
\label{eq:Q_sf}
	\mathbf{Q}=
 \begin{bmatrix}
	h \\
	q_x \\
	q_y \\
	I \\
  \eta \\
	\hdashline[2pt/2pt]
  [M_{ak}]  
 \end{bmatrix}
\end{equation}

The diffusive matrix in $x$ direction is:
%\begin{adjustwidth}{-4cm}{-3cm}
\begin{equation}
\label{eq:D_x}
\makebox[\textwidth][c]{$
		\mathbf{D_x}=\left[
		\begin{array}{ccccc;{2pt/2pt}c}
  0 & 0 & 0 & 0 & 0 & \left[\mathbb{0}\right] \\
	0 & 0 & 0 & 0 & 0 & \left[\mathbb{0}\right] \\
	0 & 0 & 0 & 0 & 0 & \left[\mathbb{0}\right] \\
	0 & 0 & 0 & -D_H & 0 & \left[\mathbb{0}\right] \\
	0 & 0 & 0 & 0 & \frac{\partial q_{bx}}{\partial \frac{\partial \eta}{\partial x}} & \left[\mathbb{0}\right] \\
	\hdashline[2pt/2pt]
	\left[\mathbb{0}\right] & \left[\mathbb{0}\right] & \left[\mathbb{0}\right] & \left[\mathbb{0}\right] & \left[\frac{\partial q_{bkx}}{\partial \frac{\partial \eta}{\partial x}}-f^{I}_{k}\frac{\partial q_{bx}}{\partial \frac{\partial \eta}{\partial x}}\right] & \left[\mathbb{0}\right] \\
%
 \end{array}\right]
$}
\end{equation}
%\end{landscape}

The diffusive matrix in $y$ direction is:
%\begin{adjustwidth}{-4cm}{-3cm}
\begin{equation}
\label{eq:D_y}
\makebox[\textwidth][c]{$
		\mathbf{D_y}=\left[
		\begin{array}{ccccc;{2pt/2pt}c}
  0 & 0 & 0 & 0 & 0 & \left[\mathbb{0}\right] \\
	0 & 0 & 0 & 0 & 0 & \left[\mathbb{0}\right] \\
	0 & 0 & 0 & 0 & 0 & \left[\mathbb{0}\right] \\
	0 & 0 & 0 & -D_H & 0 & \left[\mathbb{0}\right] \\
	0 & 0 & 0 & 0 & \frac{\partial q_{by}}{\partial \frac{\partial \eta}{\partial y}} & \left[\mathbb{0}\right] \\
	\hdashline[2pt/2pt]
	\left[\mathbb{0}\right] & \left[\mathbb{0}\right] & \left[\mathbb{0}\right] & \left[\mathbb{0}\right] & \left[\frac{\partial q_{bky}}{\partial \frac{\partial \eta}{\partial y}}-f^{I}_{k}\frac{\partial q_{by}}{\partial \frac{\partial \eta}{\partial y}}\right] & \left[\mathbb{0}\right] \\
%
 \end{array}\right]
$}
\end{equation}
%\end{landscape}

The matrix of cross derivatives is:
\begin{equation}
\label{eq:C}
\makebox[\textwidth][c]{$
		\mathbf{C}=\left[
		\begin{array}{ccccc;{2pt/2pt}c}
  0 & 0 & 0 & 0 & 0 & \left[\mathbb{0}\right] \\
	0 & 0 & 0 & 0 & 0 & \left[\mathbb{0}\right] \\
	0 & 0 & 0 & 0 & 0 & \left[\mathbb{0}\right] \\
	0 & 0 & 0 & 0 & 0 & \left[\mathbb{0}\right] \\
	0 & 0 & 0 & 0 & \frac{\partial q_{bx}}{\partial \frac{\partial \eta}{\partial y}}+\frac{\partial q_{by}}{\partial \frac{\partial \eta}{\partial x}} & \left[\mathbb{0}\right] \\
	\hdashline[2pt/2pt]
	\left[\mathbb{0}\right] & \left[\mathbb{0}\right] & \left[\mathbb{0}\right] & \left[\mathbb{0}\right] & \left[\frac{\partial q_{bkx}}{\partial \frac{\partial \eta}{\partial y}}-f^{I}_{k}\frac{\partial q_{bx}}{\partial \frac{\partial \eta}{\partial y}}+\frac{\partial q_{bky}}{\partial \frac{\partial \eta}{\partial x}}-f^{I}_{k}\frac{\partial q_{by}}{\partial \frac{\partial \eta}{\partial x}}\right] & \left[\mathbb{0}\right] \\
%
 \end{array}\right]
$}
\end{equation}
%\end{landscape}

%\begin{landscape}
%\thispagestyle{empty}
The system matrix in $x$ direction is:
\begin{equation}
\label{eq:Ax_sf}
\makebox[\textwidth][c]{$
		\mathbf{A}_x=\left[
 \begin{array}{ccccc;{2pt/2pt}c}
%wmass 
  0 & 1 & 0 & 0 & 0 & \left[\mathbb{0}\right] \\
%water mom x
	gh-\left(\frac{q_x}{h}\right)^2 & 2\frac{q_x}{h}-\frac{\partial T^\prime_{xx}}{\partial q_x} & -\frac{\partial T^\prime_{xx}}{\partial q_y} & -\frac{\partial T^\prime_{xx}}{\partial I} & gh & \left[\mathbb{0}\right] \\
%water mom y
	\frac{-q_{x}q_{y}}{h^2} & \frac{q_y}{h}-\frac{\partial T^\prime_{yx}}{\partial q_x} & \frac{q_x}{h}-\frac{\partial T^\prime_{yx}}{\partial q_y} & -\frac{\partial T^\prime_{yx}}{\partial I} & 0 & \left[\mathbb{0}\right] \\
%const I
0 & -S_{sxx} & -S_{syx} & \frac{q_x}{h} & 0 & \left[\mathbb{0}\right] \\
%exner
	\frac{\partial q_{bx}}{\partial h} & \frac{\partial q_{bx}}{\partial q_x} & \frac{\partial q_{bx}}{\partial q_y} & \frac{\partial q_{bx}}{\partial I} & 0 & \left[\frac{\partial q_{bx}}{\partial M_{al}}\right] \\
%hirano
	\hdashline[2pt/2pt]
	\left[\frac{\partial q_{bkx}}{\partial h}-f^{I}_{k}\frac{\partial q_{bx}}{\partial h}\right] & \left[\frac{\partial q_{bkx}}{\partial q_x}-f^{I}_{k}\frac{\partial q_{bx}}{\partial q_x}\right] & \left[\frac{\partial q_{bkx}}{\partial q_y}-f^{I}_{k}\frac{\partial q_{bx}}{\partial q_y}\right] & \left[\frac{\partial q_{bkx}}{\partial I}-f^{I}_{k}\frac{\partial q_{bx}}{\partial I}\right] & \left[\mathbb{0}\right] & \left[\frac{\partial q_{bkx}}{\partial M_{al}}-f^{I}_{k}\frac{\partial q_{bx}}{\partial M_{al}}\right]
%end
 \end{array}\right]
$}
\end{equation}
%\end{landscape}

\vspace{1cm}
The system matrix in $y$ direction is:

%\begin{landscape}
%\thispagestyle{empty}
\begin{equation}
\label{eq:Ay_sf}
\makebox[\textwidth][c]{$
		\mathbf{A}_y=\left[
 \begin{array}{ccccc;{2pt/2pt}c}
%water mass
  0 & 0 & 1 & 0 & 0 & \left[0\right] \\
%water mom x
	\frac{-q_{x}q_{y}}{h^2} & \frac{q_y}{h}-\frac{\partial T^\prime_{xy}}{\partial q_x} & \frac{q_x}{h}-\frac{\partial T^\prime_{xy}}{\partial q_y} & -\frac{\partial T^\prime_{xy}}{\partial I} & 0 & \left[\mathbb{0}\right] \\
%water mom y	
	gh-\left(\frac{q_y}{h}\right)^2 & -\frac{\partial T^\prime_{yy}}{\partial q_x} & 2\frac{q_y}{h}-\frac{\partial T^\prime_{yy}}{\partial q_y} & -\frac{\partial T^\prime_{yy}}{\partial I} & gh & \left[\mathbb{0}\right] \\
%const I
0 & -S_{sxy} & -S_{syy} & \frac{q_y}{h} & 0 & \left[\mathbb{0}\right] \\
%exner
	\frac{\partial q_{by}}{\partial h} & \frac{\partial q_{by}}{\partial q_x} & \frac{\partial q_{by}}{\partial q_y} & \frac{\partial q_{by}}{\partial I} & 0 & \left[\frac{\partial q_{by}}{\partial M_{al}}\right] \\
%hirano
	\hdashline[2pt/2pt]
	\left[\frac{\partial q_{bky}}{\partial h}-f^{I}_{k}\frac{\partial q_{by}}{\partial h}\right] & \left[\frac{\partial q_{bky}}{\partial q_x}-f^{I}_{k}\frac{\partial q_{by}}{\partial q_x}\right] & \left[\frac{\partial q_{bky}}{\partial q_y}-f^{I}_{k}\frac{\partial q_{by}}{\partial q_y}\right] & \left[\frac{\partial q_{bky}}{\partial I}-f^{I}_{k}\frac{\partial q_{by}}{\partial I}\right] & \left[\mathbb{0}\right] & \left[\frac{\partial q_{bky}}{\partial M_{al}}-f^{I}_{k}\frac{\partial q_{by}}{\partial M_{al}}\right]
%
 \end{array}\right]
$}
\end{equation}
%\end{landscape}

The vector of source terms is:
\begin{equation}
\label{eq:S_sf}
		\mathbf{S}=
 \begin{bmatrix}
  0 \\
	ghS_{fx} \\
	ghS_{fy} \\
	\frac{QI}{h^2L^*_I} \\
	0 \\
	\hdashline[2pt/2pt]
	\left[\mathbb{0}\right] \\
 \end{bmatrix}
\end{equation}
%
%----------------------------
%
\section{Perturbation Analysis}
%
%----------------------------
%
\label{sec:pa}
%
In this Section we conduct a perturbation analysis of the model equations (Section \ref{subsec:lin}). We test the linear analysis by means of numerical simulations (Section \ref{subsec:nosec}).
%
%
%----------------------------
\subsection{Linearization}
%----------------------------
%
%
\label{subsec:lin}
%
We consider a reference state of dependent variables $\mathbf{Q_0}$ which is a solution of the equations and a small perturbation to the state $\mathbf{Q'}$ so that $\mathbf{Q}=\mathbf{Q_0}+\mathbf{Q'}$. The reference state is that of steady uniform straight flow in an arbitrary direction over a flat sloping bed composed of an arbitrary uniform grain size distribution. 

Mathematically, $h_0=\text{c1}$, $q_{x0}=\text{c2}$, $q_{y0}=\text{c3}$, $I_0=0$, $\frac{\partial \eta}{\partial x}=\text{c4}$, $\frac{\partial \eta}{\partial y}=\text{c5}$, $M_{ak0}=\text{c6}$ $\forall k \in \{1,N-1\}$, where c$X$ stands for a constant such that:
\begin{equation}
\label{eq:matrixf_sf_0}
\begin{split}
	\frac{\partial \mathbf{Q_0}}{\partial t}+\mathbf{D_x}\left(\mathbf{Q_0}\right)\frac{\partial^2 \mathbf{Q_0}}{\partial x^2}+\mathbf{D_y}\left(\mathbf{Q_0}\right)\frac{\partial^2 \mathbf{Q_0}}{\partial y^2}+\mathbf{C}\left(\mathbf{Q_0}\right)\frac{\partial^2 \mathbf{Q_0}}{\partial x \partial y}+\mathbf{A_x}\left(\mathbf{Q_0}\right)\frac{\partial \mathbf{Q_0}}{\partial x}+\mathbf{A_y}\left(\mathbf{Q_0}\right)\frac{\partial \mathbf{Q_0}}{\partial y}+\mathbf{S}\left(\mathbf{Q_0}\right)=0 \implies \\
		\mathbf{A_x}\left(\mathbf{Q_0}\right)\frac{\partial \mathbf{Q_0}}{\partial x}+\mathbf{A_y}\left(\mathbf{Q_0}\right)\frac{\partial \mathbf{Q_0}}{\partial y}+\mathbf{S}\left(\mathbf{Q_0}\right)=0 
\end{split}
\end{equation}

We substitute the perturbed solution into the system of equations (\ref{eq:matrixf_sf}):
\begin{equation}
\label{eq:matrixf_sf_s1}
\begin{split}
	\frac{\partial \mathbf{Q}}{\partial t}+\mathbf{D_x}\left(\mathbf{Q}\right)\frac{\partial^2 \mathbf{Q}}{\partial x^2}+\mathbf{D_y}\left(\mathbf{Q}\right)\frac{\partial^2 \mathbf{Q}}{\partial y^2}+\mathbf{C}\left(\mathbf{Q}\right)\frac{\partial^2 \mathbf{Q}}{\partial x \partial y}+\mathbf{A_x}\left(\mathbf{Q}\right)\frac{\partial \mathbf{Q}}{\partial x}+\mathbf{A_y}\left(\mathbf{Q}\right)\frac{\partial \mathbf{Q}}{\partial y}+\mathbf{S}\left(\mathbf{Q}\right)=0\implies\\
	\frac{\partial \left(\mathbf{Q_0}+\mathbf{Q'}\right)}{\partial t}+\mathbf{D_x}\left(\mathbf{Q_0}+\mathbf{Q'}\right)\frac{\partial^2 \left(\mathbf{Q_0}+\mathbf{Q'}\right)}{\partial x^2}+\mathbf{D_y}\left(\mathbf{Q_0}+\mathbf{Q'}\right)\frac{\partial^2 \left(\mathbf{Q_0}+\mathbf{Q'}\right)}{\partial y^2}+\mathbf{C}\left(\mathbf{Q_0}+\mathbf{Q'}\right)\frac{\partial^2 \left(\mathbf{Q_0}+\mathbf{Q'}\right)}{\partial x \partial y}+\\
	+\mathbf{A_x}\left(\mathbf{Q_0}+\mathbf{Q'}\right)\frac{\partial \left(\mathbf{Q_0}+\mathbf{Q'}\right)}{\partial x}+\mathbf{A_y}\left(\mathbf{Q_0}+\mathbf{Q'}\right)\frac{\partial \left(\mathbf{Q_0}+\mathbf{Q'}\right)}{\partial y}+\mathbf{S}\left(\mathbf{Q_0}+\mathbf{Q'}\right)=0
\end{split}
\end{equation}

The terms of the matrices $\mathbf{D_x}$, $\mathbf{D_y}$, $\mathbf{C}$, $\mathbf{A_x}$, and $\mathbf{A_y}$ evaluated at $\left(\mathbf{Q_0}+\mathbf{Q'}\right)$ are linearized such that:
\begin{equation}
\label{eq:l_gen}
\begin{split}
	\mathbf{G}\left(\mathbf{Q_0}+\mathbf{Q'}\right)=\mathbf{G}\left(\mathbf{Q_0}\right)+\mathbf{G_{J}}\left(\mathbf{Q_0}\right)\mathbf{Q'}+\mathcal{O}\left(\mathbf{Q'}^2\right)
\end{split}
\end{equation}
where $\mathbf{G}$ is one of the matrices and $\mathbf{G_{J}}$ is a matrix that contains the first order terms. For instance, for the first term of the momentum equation in the $x$ direction we obtain:
\begin{equation}
\label{eq:lin1}
\begin{split}
	gh-\left(\frac{q_x}{h}\right)^2=g\left(h_0+h'\right)-\left(\frac{\left(q_{x0}+q'_x\right)}{\left(h_0+h'\right)}\right)^2=\\
	=gh_0+gh'-\left(q_{x0}^2+q^{'2}_{x}+2q_{x0}q'_{x}\right)\left(h_0^{-2}-2h_0^3h'\right)=\\
	=gh_0-\left(\frac{q_{x0}}{h_0}\right)^2+\left(g+\frac{q_{x0}^2}{h_0^3}\right)h'+2\frac{q_{x0}}{h_0^2}q'_x-h_0^{-2}q^{'2}_{x}+2h_0^{-3}h'q^{'2}_{x}-4h_0^{-3}q_{x0}q'_{x}h'
\end{split}
\end{equation}

Neglecting secondary terms we obtain:
\begin{equation}
\label{eq:D_lin}
\begin{split}
	\mathbf{D_j}\left(\mathbf{Q_0}+\mathbf{Q'}\right)\frac{\partial^2 \left(\mathbf{Q_0}+\mathbf{Q'}\right)}{\partial j^2}=\left[\mathbf{D_j}\left(\mathbf{Q_0}\right)+\mathbf{D_{Jj}}\left(\mathbf{Q_0}\right)\mathbf{Q'}\right]\left[\frac{\partial^2 \mathbf{Q_0}}{\partial j^2}+\frac{\partial^2 \mathbf{Q'}}{\partial j^2}\right]=\mathbf{D_j}\left(\mathbf{Q_0}\right)\frac{\partial^2 \mathbf{Q'}}{\partial j^2}
\end{split}
\end{equation}
%
\begin{equation}
\label{eq:C_lin}
\begin{split}
	\mathbf{C}\left(\mathbf{Q_0}+\mathbf{Q'}\right)\frac{\partial^2 \left(\mathbf{Q_0}+\mathbf{Q'}\right)}{\partial xy}=\left[\mathbf{C}\left(\mathbf{Q_0}\right)+\mathbf{C_{J}}\left(\mathbf{Q_0}\right)\mathbf{Q'}\right]\left[\frac{\partial^2 \mathbf{Q_0}}{\partial xy}+\frac{\partial^2 \mathbf{Q'}}{\partial xy}\right]=\mathbf{C}\left(\mathbf{Q_0}\right)\frac{\partial^2 \mathbf{Q'}}{\partial xy}
\end{split}
\end{equation}
%
\begin{equation}
\label{eq:A_lin}
\begin{split}
	\mathbf{A_j}\left(\mathbf{Q_0}+\mathbf{Q'}\right)\frac{\partial \left(\mathbf{Q_0}+\mathbf{Q'}\right)}{\partial j}=\left[\mathbf{A_j}\left(\mathbf{Q_0}\right)+\mathbf{A_{Jj}}\left(\mathbf{Q_0}\right)\mathbf{Q'}\right]\left[\frac{\partial \mathbf{Q_0}}{\partial j}+\frac{\partial \mathbf{Q'}}{\partial j}\right]=\mathbf{A_{Jj}}\left(\mathbf{Q_0}\right)\mathbf{Q'}\frac{\partial \mathbf{Q_0}}{\partial j}+\mathbf{A_j}\left(\mathbf{Q_0}\right)\frac{\partial \mathbf{Q'}}{\partial j}
\end{split}
\end{equation}
where subindex $j=x,y$. We have used that $\frac{\partial^2 \mathbf{Q_0}}{\partial j^2}=0$ and that $\mathbf{G_{J}}\left(\mathbf{Q_0}\right)\mathbf{Q'}\frac{\partial^2 \mathbf{Q'}}{\partial j^2}$ is negligible. Note that $\frac{\partial \mathbf{Q_0}}{\partial j}\neq0$ because the bed slope is not zero. In this case:
%
\begin{equation}
\label{eq:eta_lin2}
\begin{split}
 gh\frac{\partial \eta}{\partial j}=g\left(h_0+h'\right)\frac{\partial \left(\eta_0+\eta'\right)}{\partial j}=gh_0\frac{\partial \eta_0}{\partial j}+gh_0\frac{\partial \eta'}{\partial j}+gh'\frac{\partial \eta_0}{\partial j}
\end{split}
\end{equation}
where we have neglected secondary terms. The first right hand side term is part of $\mathbf{A_j}\left(\mathbf{Q_0}\right)\frac{\partial \mathbf{Q_0}}{\partial j}$, the second term is part of $\mathbf{A_j}\left(\mathbf{Q_0}\right)\frac{\partial \mathbf{Q'}}{\partial j}$ and the third term is the only term of $\mathbf{A_{Jj}}\left(\mathbf{Q_0}\right)\mathbf{Q'}\frac{\partial \mathbf{Q_0}}{\partial j}$. Thus, we write:
%
\begin{equation}
\label{eq:AJt}
\makebox[\textwidth][c]{$
		\mathbf{A_{Jt}}=\left(\mathbf{A_{Jx}}\left(\mathbf{Q_0}\right)\frac{\partial \mathbf{Q_0}}{\partial x}+\mathbf{A_{Jy}}\left(\mathbf{Q_0}\right)\frac{\partial \mathbf{Q_0}}{\partial y}\right)\mathbf{Q'}=\left[
		\begin{array}{ccccc;{2pt/2pt}c}
  0 & 0 & 0 & 0 & 0 & \left[\mathbb{0}\right] \\
	g\frac{\partial \eta_0}{\partial x} & 0 & 0 & 0 & 0 & \left[\mathbb{0}\right] \\
	g\frac{\partial \eta_0}{\partial y} & 0 & 0 & 0 & 0 & \left[\mathbb{0}\right] \\
	0 & 0 & 0 & 0 & 0 & \left[\mathbb{0}\right] \\
	0 & 0 & 0 & 0 & 0 & \left[\mathbb{0}\right] \\
	\hdashline[2pt/2pt]
	\left[\mathbb{0}\right] & \left[\mathbb{0}\right] & \left[\mathbb{0}\right] & \left[\mathbb{0}\right] & \left[\mathbb{0}\right] & \left[\mathbb{0}\right] \\
%
 \end{array}\right]
$}
\end{equation}
where $\frac{\partial \eta_0}{\partial j}=-\frac{C_fq_{j0}Q_0}{gh_0^3}$ which we obtain from the fact that $\mathbf{Q_0}$ is a uniform steady solution.

The source term is linearized:
\begin{equation}
\label{eq:S_lin}
\begin{split}
	\mathbf{S}\left(\mathbf{Q_0}+\mathbf{Q'}\right)=\mathbf{S}\left(\mathbf{Q_0}\right)+\mathbf{S_{J}}\left(\mathbf{Q_0}\right)\mathbf{Q'}+\mathcal{O}\left(\mathbf{Q'}^2\right)
\end{split}
\end{equation}
where $\mathbf{S_{J}}$ is the jacobian of $S$:
\begin{equation}
\label{eq:SJ}
\makebox[\textwidth][c]{$
		\mathbf{S_J}=\left[
		\begin{array}{ccccc;{2pt/2pt}c}
  0 & 0 & 0 & 0 & 0 & \left[\mathbb{0}\right] \\
	\frac{-2QC_fq_x}{h^3} & \frac{C_f\left(Q^2+q_x^2\right)}{h^2Q} & \frac{C_fq_xq_y}{Qh^2} & 0 & 0 & \left[\mathbb{0}\right] \\
	\frac{-2QC_fq_y}{h^3} & \frac{C_fq_xq_y}{Qh^2} & \frac{C_f\left(Q^2+q_y^2\right)}{h^2Q} & 0 & 0 & \left[\mathbb{0}\right] \\
	\frac{-2QI}{L_I^*h^3} & \frac{q_xI}{L_I^*h^2Q} & \frac{q_yI}{L_I^*h^2Q} & \frac{Q}{L_I^* h^2} & 0 & \left[\mathbb{0}\right] \\
	0 & 0 & 0 & 0 & 0 & \left[\mathbb{0}\right] \\
	\hdashline[2pt/2pt]
	\left[\mathbb{0}\right] & \left[\mathbb{0}\right] & \left[\mathbb{0}\right] & \left[\mathbb{0}\right] & \left[\mathbb{0}\right] & \left[\mathbb{0}\right] \\
%
 \end{array}\right]
$}
\end{equation}

	%[                 0,                              0,                              0,          0]
	%[  -(2*Q*cf*qx)/h^3, (Q*cf)/h^2 + (cf*qx^2)/(Q*h^2),             (cf*qx*qy)/(Q*h^2),          0]
	%[  -(2*Q*cf*qy)/h^3,             (cf*qx*qy)/(Q*h^2), (Q*cf)/h^2 + (cf*qy^2)/(Q*h^2),          0]
	%[ -(2*I*Q)/(LI*h^3),              (I*qx)/(LI*Q*h^2),              (I*qy)/(LI*Q*h^2), Q/(LI*h^2)]

Substituting Equations (\ref{eq:D_lin}), (\ref{eq:C_lin}), (\ref{eq:A_lin}), (\ref{eq:AJt}), and (\ref{eq:S_lin}) in Equation (\ref{eq:matrixf_sf_s1}) and using the fact $\mathbf{Q_0}$ satisfies Equation (\ref{eq:matrixf_sf} we obtain:
\begin{equation}
	\label{eq:adv_p}
	\frac{\partial \mathbf{Q'}}{\partial t}+\mathbf{D_{x0}}\frac{\partial^2 \mathbf{Q'}}{\partial x^2}+\mathbf{D_{y0}}\frac{\partial^2 \mathbf{Q'}}{\partial y^2}+\mathbf{C_0}\frac{\partial^2 \mathbf{Q'}}{\partial x \partial y}+\mathbf{A_{x0}}\frac{\partial \mathbf{Q'}}{\partial x}+\mathbf{A_{y0}}\frac{\partial \mathbf{Q'}}{\partial y}+\mathbf{B_{0}}\mathbf{Q'}=0
\end{equation}
where the subindex $\mathbf{_0}$ indicates that the matrix is evaluated at the unperturbed state and $\mathbf{B}=\mathbf{S_{J}}+\mathbf{A_{Jt}}$. 

Assuming a wave type perturbation $\mathbf{Q'}=\mathbf{\hat{Q}}\left(t\right)e^{ik_{wx} x+ik_{wy} y}$ where $i$ is the imaginary unit and $k_{wx}$ and $k_{wy}$ are wave numbers in $x$ and $y$ direction respectively and $\mathbf{\hat{Q}}\left(t\right)$ is a time dependent amplitude we obtain the system:
\begin{equation}
\begin{split}
	\label{eq:adv_pw}
	\frac{\partial \mathbf{\hat{Q}}\left(t\right)}{\partial t}-\mathbf{D_{x0}}k_{wx}^2\mathbf{\hat{Q}}\left(t\right)-\mathbf{D_{y0}}k_{wy}^2\mathbf{\hat{Q}}\left(t\right)-\mathbf{C_0}k_{wx} k_{wy}\mathbf{\hat{Q}}\left(t\right)+\mathbf{A_{x0}}ik_{wx}\mathbf{\hat{Q}}\left(t\right)+\mathbf{A_{y0}}ik_{wy}\mathbf{\hat{Q}}\left(t\right)+\mathbf{B_{0}}\mathbf{\hat{Q}}\left(t\right)=0\implies \\
		\frac{\partial \mathbf{\hat{Q}}\left(t\right)}{\partial t}=\left[\mathbf{D_{x0}}k_{wx}^2+\mathbf{D_{y0}}k_{wy}^2+\mathbf{C_0}k_{wx} k_{wy}-\mathbf{A_{x0}}ik_{wx}-\mathbf{A_{y0}}ik_{wy}-\mathbf{B_{0}}\right]\mathbf{\hat{Q}}\left(t\right) \;.
\end{split}
\end{equation}
Note that we have assumed that the perturbation can be represented as a Fourier series which implies that it is bounded for $x=\pm \inf$. Using the properties of the eigenvalues we obtain the ordinary differential equation:
\begin{equation}
	\label{eq:ode}
	\frac{\mathrm{d} \mathbf{\Phi}}{\mathrm{d} t}=\lambda\mathbf{\Phi} \;,
\end{equation}
where $\mathbf{\Phi}$ $[-]$ and $\lambda$ [\si{1/s}] are an eigenvector and eigenvalue of matrix:
\begin{equation}
	\label{eq:R}
	\mathbf{R}=\mathbf{D_{x0}}k_{wx}^2+\mathbf{D_{y0}}k_{wy}^2+\mathbf{C_0}k_{wx} k_{wy}-\mathbf{A_{x0}}ik_{wx}-\mathbf{A_{y0}}ik_{wy}-\mathbf{B_{0}} \;,
\end{equation}
respectively. The real part of $\lambda_m$ $\forall m \in\{1,N+4\}$ is the growth/decay rate (if positive/negative) of the perturb solution. The imaginary part yields the angular frequency. 

%
%
%----------------------------
\subsection{Test Using Numerical Simulations}
%----------------------------
%
%
\label{subsec:nosec}
%
Before we use the linear solution of the model equations obtained in Section \ref{subsec:lin}, in this Section we test the solution. 

To test the linear solution we study the initial development of free alternate bars in straight channels which are recognized to be the result of an instability mechanism \citep{Engelund73,Parker76,Fredsoe78}. The main result of the liner analysis of free alternate bars is that there exists a width to depth ratio threshold below which alternate bars decay, and grow if the ratio is larger than the critical value \citep{Olesen83,Kuroki84,Colombini87,Tubino99}. In other words, the linear analysis predicts the growth or decay of a combinations of wavenumbers in streamwise and transverse direction for a given a reference state. We first compare the growth and decay domain predicted by our solution to the domain predicted by \citet{Olesen83}, \citet{Colombini87}, and \citet{Schielen93} and we obtain similar domains. The agreement is not perfect since each model assumes certain simplifications. For instance, all the previous models assume steady flow in order to obtain analytical expressions of the growth domain while we include the effect of the flow unsteadiness. 

To test our solution thoroughly we compare the linear theory to the results of numerical simulations that solve the same equations that we have analyzed following a similar strategy as \citet{Siviglia13}. We use standard Delft3D (Version 6.02.13.7545M compiled for Windows 32bit). A double-periodic perturbation (characterized by a certain wave number in $x$ and $y$ direction) is added to a uniform flow and if the real part of all the eigenvalues of the resulting system, Equation (\ref{eq:R}), are negative, the perturbations will decay. If at least one is positive, they will grow. 

We consider a uniform state with the parameters shown in Table \ref{tab:L_ph}. We consider unisize sediment. The sediment transport rate is computed using the relation by \citet{Engelund67}. The bed slope effects are taken into account using the relation by \citet{Sekine92_1}. This is equivalent to use the relation by \citet{Koch80}, Equation (\ref{eq:g_sk}), with parameters $A_s=1$, $B_s=0$, $C_s=0$, and $D_s=0$. The instability domain of the reference state is plotted in Figure \ref{fig:L_re}. The curves are the separatrix between the linearly stable and unstable domains (i.e., growth rate, real part of the eigenvalues, equal to 0). The continuous line is obtained assuming no secondary flow while the discontinuous line considers secondary flow with a diffusion coefficient $D_H=\SI{5}{m^2/s}$. The left plot presents the domain in terms of wave number ($k_{wi}$ [\si{m^{-1}}]). The right plot present the same domain in terms of wavelength ($l_{wi}$ [\si{m}]). These are related by the equation $k_{wi}=\frac{2*\pi}{l_{wi}}$ for $i=x,y$. Due to the impermeable boundary conditions at the closed domains (i.e., the river banks), only natural multiples of $\pi$ are valid transverse wave numbers. Since the first mode (i.e., transverse bars) is the most unstable one \citep{Colombini87,Schielen93}, the river width $B$ [\si{m}] is obtained as $B=\pi/k_{wy}$. Note that secondary flow reduces the unstable domain thus reducing the conditions in which bars grow but it is not a necessary mechanism to take into account to predict bar formation \citep{Schielen93,Iwasaki17}. Provided that we account for the advection and diffusion of secondary flow intensity, the result that the growth domain decreases when considering secondary flow needs to be treated with caution. This result can be strongly dominated by the diffusion coefficient and it is possible that if secondary flow is considered in equilibrium (i.e., there is no advection diffusion equation of secondary flow intensity) the result is opposite. 
%
\begin{table}[ht]
%\tiny
	\begin{center}
	\begin{tabular}{ccccccc}		
			\toprule
$u$ [\si{m/s}] & $v$ [\si{m/s}] & $h$ [\si{m}] & $C_f$ [\si{-}] & $d_k$ [\si{m}] & $s$ [$-$] & $q_{bk}$ [\si{m^2/s}]\\
			\midrule
$1$ & $0$ & $1$ & $0.007$ & $0.001$ & \num{7.13541e-4} & \num{1.118e-4} \\ 
			\bottomrule
			\end{tabular}
      \caption{Physical parameters of the reference state to test linear analysis.}
      \label{tab:L_ph}
	\end{center}
\end{table}
%
\begin{figure}[ht]
    \centering
    \includegraphics{L_re.pdf}
    \caption{Domain of linear stability of the reference situation (Table \ref{tab:L_ph}) in terms of: (\textbf{a}) wavenumber, and (\textbf{b}) wave length. The continuous line is the separatrix for the case without secondary flow while the dashed line is the separatrix considering secondary flow and a diffusion coefficient $D_H=\SI{5}{m^2/s}$. The circle and square mark the conditions of the numerical simulations (Table \ref{tab:L_num_diff}).}
    \label{fig:L_re}
\end{figure}

We conduct three numerical simulations varying the width and the consideration of secondary flow (Figure \ref{fig:L_re}). The simulations start under the equilibrium conditions of the uniform reference state in which a double-periodic perturbation in bed elevation that satisfies the boundary conditions is added:
\begin{equation}
	\label{eq:eta_p}
 \eta'=A\sin\left(\frac{\pi y}{B}\right)\cos\left(\frac{2 \pi x}{L_b}-\frac{\pi}{2}\right) \;,
\end{equation}
where $A=\SI{0.005}{m}$ is the amplitude, $B$ [\si{m}] is the width, and $L_b=\SI{100}{m}$ is the streamwise wavelength. The channel is aligned with the $x$ axis with the origin at the center of the channel. Table \ref{tab:L_num_eq} presents the parameters which are equal for all simulations and Table \ref{tab:L_num_diff} presents the ones which are different. The simulations have a spin up time of \SI{1800}{s} for the flow to adapt to the perturbed initial condition. We have visually checked that after the spin up time the flow is steady. 
%
\begin{table}[ht]
%\tiny
	\begin{center}
	\begin{tabular}{ccccc}		
			\toprule
$L$ [\si{m}] & $\Delta t$ [\si{s}] & $\Delta x$ [\si{m}] & $\Delta y$ [\si{m}] & $T$ [\si{s}] \\
			\midrule
$1000$ & $1$ & $1$ & $1$ & $(1800)+3600$ \\ 
			\bottomrule
			\end{tabular}
      \caption{Parameters which are equal for the simulations to test the linear analysis. The symbols not defined previously are the length of the domain $L$ [\si{m}], the time step $\Delta t$ [\si{s}], the space step in the $x$ direction $\Delta x$ [\si{m}], the space step in $y$ direction [\si{m}], and the simulation time $T$ [\si{s}]. The time in parenthesis is the spin-up time.}
      \label{tab:L_num_eq}
	\end{center}
\end{table}
%
\begin{table}[ht]
%\tiny
	\begin{center}
	\begin{tabular}{ccccccccc}		
			\toprule
  Simulation & $B$ [\si{m}] & sec. flow & linear analysis \\
			\midrule
L1 & 22 & N & growth \\ 
L2 & 17 & N & decay \\ 
L3 & 22 & Y & decay \\ 
			\bottomrule
			\end{tabular}
      \caption{Parameters which are different for the simulations to test the linear analysis.}
      \label{tab:L_num_diff}
	\end{center}
\end{table}

In Figure \ref{fig:L_t} we show the evolution of the flow depth with time at $y=B/2$ for a location in which initially we find the trough of a bar (i.e., maximum flow depth). For the three numerical simulations the initial perturbation grows or decays as predicted in the linear analysis. 
%
\begin{figure}[ht]
    \centering
    \includegraphics{L_t.pdf}
    \caption{Evolution of the flow depth with time at $y=B/2$ for a location in which initially we find the trough of a bar (i.e., maximum flow depth) for simulation (\textbf{a}) L1, (\textbf{b}) L2, and (\textbf{c}) L3 (Tables \ref{tab:L_num_eq} and \ref{tab:L_num_diff}).}
    \label{fig:L_t}
\end{figure}
%
%----------------------------
%
\section{Results of the Perturbation Analysis}
%
%----------------------------
%
\label{sec:res_lin}
%
In this section we apply the results of the linear analysis to study the mathematical character of the system of equations depending on the secondary flow (Section \ref{subsec:secflow_uni}) bed slope effects (Section \ref{subsec:bedslope_uni}), and mixed-size sediment (Section \ref{subsec:mixsed}). First, in Section \ref{subsec:def_ill} we introduce the topic of ill-posedness.
%
%
%----------------------------
\subsection{Definition of Ill-posedness}
%----------------------------
%
%
\label{subsec:def_ill}
%
In this section we explain the concept of ill-posedness and the reasons that make this concept important in modeling river morphodynamics.

The first step in modeling a physical process is its mathematization. That is, to interpret the phenomena and to represent the physical properties (e.g, flow velocity, water elevation, et cetera) as variables related to each other by means of equations subject to certain conditions. A first error is introduced, which is due to the difference between the actual physical processese occurring in nature and the system of equations that represents them. The second step concerns the solution of the system of equations. Since in general terms the problem is too complex to solve it analytically, numerical approximations are used to find a solution. A second error is introduced due to the difference between the ``physical'' and ``numerical'' set of equations.

A system of equations and its conditions (i.e., initial and boundary), needs to be classified as well-posed to be representable of the physical process under study \citep[e.g.][]{Joseph90}. A problem is well-posed if it fulfills three conditions: (1) a solution exists, (2) the solution is unique, and (3) the solution is continuous on the data \citep{Hadamard22}. The first two conditions are trivial. The third conditions means that an infinitesimal perturbation to the data (initial and boundary conditions) cannot yield a finite divergence in the solution in an arbitrarily short time. Worded differently, in a well-posed problem, as the size of the perturbation in the data tends to zero, the size of the perturbation in the solution also tends to zero at any fixed time. 

The remark on the arbitrarily short time is important since this is what makes chaotic systems not to be ill-posed. Deterministic chaos is also sensitive to initial conditions, which implies that models lose their predictive capabilities after a certain time \citep{Devaney89,Banks92}. Yet, there exists a finite time in which the dynamics are predictable. Moreover, chaotic problems are unpredictable after a certain time but continue to be predictable in statistical terms. 

To classify a problem we recall the linearized eigenvalue problem, Equation (\ref{eq:ode}). The above definition of ill-posedness is analogous to say that a problem is ill-posed if the growth rate of at least one perturbation (i.e., the real part of at least one eigenvalue of matrix $\mathbf{R}$, Equation (\ref{eq:R})) is positive for an infinite value of the wavenumber. Note that in the continuous limit a perturbation may be composed of an infinite amount of wavenumbers. Since the growth rate depends exponentially on the wavenumber, if it is positive for increasing wavenumbers, the perturbation is unbounded and there is no wavenumber cutoff after which we can neglect the remaining terms. A solution of this type is physically unrealistic. 

When the space domain is discretized the maximum wavenumber is limited by the grid size. Thus, in numerical simulation there exists a wavenumber cutoff. However, since in ill-posed problems larger wavenumbers give rise to faster growth of perturbations, the numerical solution of an ill-posed problem keeps changing for a decreasing grid size \citep{Joseph90,Kabanikhin08}. The solution does not converge. 

The numerical solution of a well-posed physical model may be ill-posed. This is the case, for instance, if the CFL number \citep{Courant28} of an explicit numerical scheme of an advection equation is larger than 1. The solution presents unphysical oscillations and it is straightforward to acknowledge that in this case the numerical discretization is not correct. Another possibility is that the numerical solution of an ill-posed model is well-posed. Contrary to the previous case, this one is of cumbersome identification since, if one does not know that the physical model is ill-posed, nothing in the solution gives the idea that there is an inherent problem. A situation in which this flaw occurs is, for instance, the inviscid two-fluid problem (i.e., the Kelvin-Helmholtz instability in vortex sheets \citep{Kelvin71}). This problem is ill-posed but diffusion introduce by the numerical scheme may provide a well-posed numerical model. Although the solution may seem realistic this is conceptually wrong since the solution is realistic for the wrong reason. The numerical diffusion solves the problems due to the lack of physical diffusion and a change of the grid and the solution does not converge with the grid since a change of grid changes the numerical diffusion. 

Let's assume that the model under study is formed out of one partial differential equation (one eigenvalue) with a second order derivative (diffusion), a first order derivative, and a source term (linear term). In this case, the mathematical character of the model depends exclusively on the diffusion coefficient since it is multiplied by the square of the wave number (Equation (\ref{eq:R})). If there is no diffusion, the first order term sets the character of the problem. In any case the source term does not play a role since it is negligible for an increasing wavenumber. The problem becomes more complicated if there is more than one equation and they are coupled. In a model with two equations, one of them may have a negative diffusion coefficient, which implies that a model with only this equation is well-posed, but the two-equation model may be ill-posed. 
%
%
%----------------------------
\subsection{Secondary Flow}
%----------------------------
%
%
\label{subsec:secflow_uni}
%
In this section we study the effect of secondary flow. That is, we study the consequences of modifying the momentum equations, (\ref{eq:wmomx_sf}) and (\ref{eq:wmomy_sf}) with stress terms related to the effect and secondary flow together with a constitutive equation to model the advection and diffusion of secondary flow as regards to the well-posedness of the model. 

We gain empirical insight studying a reference case with the parameters shown in Table \ref{tab:L_ph}. We first consider a diffusion coefficient equal to 0 (S1, Table \ref{tab:s_var}). We numerically compute the eigenvalues of matrix $\mathbf{R}$, Equation \ref{eq:R}, for varying wavenumber in the $x$ and $y$ direction. The growth domain is shown in Figure \ref{fig:dom_S1}. Note that this situation is ill-posed (Section \ref{subsec:def_ill}). 
%
\begin{table}[ht]
%\tiny
	\begin{center}
	\begin{tabular}{ccccc}		
			\toprule
Case & $D_{H}$ [\si{m^2/s}] & $A_s$ [$-$] & $\beta_c$ [$-$] & Math. charac. \\
			\midrule
S1 & 0.0 & 1          & 1 & I \\ 
S2 & 0.0 & \num{1e-4} & 1 & I \\ 
S3 & 0.5 & 1          & 1 & I \\ 
S4 & 1.0 & 1          & 1 & W \\ 
S5 & 0.0 & 1          & 0 & W \\ 
 		\bottomrule
			\end{tabular}
      \caption{Parameters that are different between secondary flow cases. The parameters that are equal are shown in Table \ref{tab:L_ph}. W=well-posed, I=ill-posed.}
      \label{tab:s_var}
	\end{center}
\end{table}
%
\begin{figure}[ht]
    \centering
    \includegraphics{dom_S1.png}
    \caption{Wave growth domain for the case S1 (Tables \ref{tab:L_ph} and \ref{tab:s_var}) as a function of the wavenumber (\textbf{a}) and wavelength (\textbf{b}). In the green area the growth rate is negative (dampening) while in the red area it is positive (growth). This situation is ill-posed (Section \ref{subsec:def_ill}).}
    \label{fig:dom_S1}
\end{figure}

A diffusion coefficient equal to 0 is physically unrealistic in our case (Section \ref{subsubsec:sf_cr}). Yet, the physically meaningful value may be below the threshold to obtain a well-posed model. We compute the necessary value of the diffusion coefficient that makes the reference case to be well-posed. This is done numerically by minimizing a function that, given matrix $\mathbf{R}(D_H)$, returns the maximum real part of the eigenvalues. Figure \ref{fig:Dh_min_1} presents the minimum diffusion coefficient as a function of the streamwise and crosswise wavenumber as well as the wavelengths. Note that for some values of wavelength, the linear analysis predicts growth but the necessary diffusion coefficient to obtain a well-posed simulation is 0. This is because the necessary value of diffusion is very small ($<\num1e-7$), below the accuracy of the minimization function. 
%
\begin{figure}[ht]
    \centering
    \includegraphics{Dh_min.png}
    \caption{Minimum diffusion coefficient for the case S1 (Tables \ref{tab:L_ph} and \ref{tab:s_var}) to obtain a well-posed model as a function the wavenumber (\textbf{a}) and wavelength (\textbf{b}).}
    \label{fig:Dh_min_1}
\end{figure}

We see that the minimum diffusion coefficient to obtain a well-posed simulation in the conditions of the reference case is slightly larger than \SI{0.8}{m^2/s}. However, this is above the theoretical values in streamwise direction, $D_{Hs}=\SI{0.49}{m^2/s}$, and transverse direction, $D_{Hn}=\SI{0.0057}{m^2/s}$ (Section \ref{subsubsec:sf_cr}). Figures \ref{fig:dom_S3} and \ref{fig:dom_S4} show the wave growth domain for diffusion coefficients equal to \SI{0.5}{m^2/s} and \SI{1}{m^2/s}, respectively (Table \ref{tab:s_var}). 
%
\begin{figure}[ht]
    \centering
    \includegraphics{dom_S3.png}
    \caption{Wave growth domain for the case S3 (Tables \ref{tab:L_ph} and \ref{tab:s_var}) as a function of the wavenumber (\textbf{a}) and wavelength (\textbf{b}). In the green area the growth rate is negative (dampening) while in the red area it is positive (growth). This situation is ill-posed (Section \ref{subsec:def_ill}).}
    \label{fig:dom_S3}
\end{figure}
%
\begin{figure}[ht]
    \centering
    \includegraphics{dom_S4.png}
    \caption{Wave growth domain for the case S4 (Tables \ref{tab:L_ph} and \ref{tab:s_var}) as a function of the wavenumber (\textbf{a}) and wavelength (\textbf{b}). In the green area the growth rate is negative (dampening) while in the red area it is positive (growth). This situation is well-posed (Section \ref{subsec:def_ill}).}
    \label{fig:dom_S4}
\end{figure}

We have studied a case with no sediment transport (fixed bed) and we obtain the same result, that the model is ill-posed if the diffusion coefficient is equal to 0. This shows that the origin of ill-posedness is on the secondary flow model itself and not on its coupling with morphology. When there is sediment transport the source of diffusion in the model is not only the horizontal diffusion coefficient of the secondary flow equation but also the bed slope effects on sediment transport. We empirically test whether bed slope effects have a significant impact on the ill-posedness due to secondary flow. To this end we compute the growth domain of a case (S2) in the same conditions as before but with a physically unrealistic value of parameter $A_s=0.0001$ in Equation (\ref{eq:g_sk}) such that diffusion due to bed slope effects is significantly increased (Figure \ref{fig:dom_S2}). The situation is ill-posed which indicates that the specific problem of ill-posedness when considering secondary flow cannot be tackled by artificially adding diffusion in the morphology equations. Note that bars cannot grow if bed slope effects have such a large effect.
%
\begin{figure}[ht]
    \centering
    \includegraphics{dom_S2.png}
    \caption{Wave growth domain for the case S2 (Tables \ref{tab:L_ph} and \ref{tab:s_var}) as a function of the wavenumber (\textbf{a}) and wavelength (\textbf{b}). In the green area the growth rate is negative (dampening) while in the red area it is positive (growth). This situation is ill-posed (Section \ref{subsec:def_ill}).}
    \label{fig:dom_S2}
\end{figure}

To account for the flow curvature not only an equation for the advection and diffusion of the secondary flow intensity is added to the shallow water equations but also the momentum equations are modified with extra stress terms (Section \ref{subsec:sec_eq}). Yet, due to added complexity in solving the extra terms, these are sometimes neglected \citep{Olesen82,Duan98,Jia99,Duan01,Kassem02,Dulal10}. \citet{Duan04} shows the consequences of neglecting this term in numerical simulations. We study the consequences of this approximation in Case S5. This case is equal to Case S1 but the extra terms in the momentum equation are not taken into consideration by setting the parameter $\beta_c$ in Equations (\ref{eq:Txx}), (\ref{eq:Txy}), and (\ref{eq:Tyy}) equal to 0. This case is well-posed (Figure \ref{fig:dom_S5}). This does not mean that the stress terms due to secondary flow in the momentum equations are the origin of ill-posedness, but that the coupling of the secondary flow model with the Shallow Water Equations may provide unrealistic results under certain conditions. 
%
\begin{figure}[ht]
    \centering
    \includegraphics{dom_S5.png}
    \caption{Wave growth domain for the case S5 (Tables \ref{tab:L_ph} and \ref{tab:s_var}) as a function of the wavenumber (\textbf{a}) and wavelength (\textbf{b}). In the green area the growth rate is negative (dampening) while in the red area it is positive (growth). This situation is well-posed (Section \ref{subsec:def_ill}).}
    \label{fig:dom_S5}
\end{figure}

We test the results of the linear analysis as regards to secondary flow running numerical simulations using standard Delft3D (Version 6.02.13.7545M compiled for Windows 32bit) in the conditions of cases S1 and S4. The domain is \SI{100}{m} long and \SI{10}{m} wide and it is discretized is square cells of length \SI{0.1}{m}. The initial condition is a flat unperturbed bed. The time step is \SI{0.1}{s} for case S1 and \SI{0.001}{s} for S4. Apparently, the larger diffusion coefficient requires a smaller time step.

In Figure \ref{fig:S1_xyd_1} we plot the flow depth after \SI{15}{s} for the simulation of Case S1. Truncation errors in the initial condition have grown significantly and continue to grow due to ill-posedness of the system of equations. After \SI{1}{min} the flow depth is physically unrealistic with flow depths as small as \SI{0.2}{m} (Figure \ref{fig:S1_xyd_2}). Figure \ref{fig:S4_xyd} presents the flow depth after \SI{1}{min} for the simulation of Case S4. This case is well-posed and no oscillations are present.
%
\begin{figure}[ht]
    \centering
    \includegraphics{S1_xyd_1.png}
    \caption{Flow depth after \SI{15}{s} in the simulation of Case S1.}
    \label{fig:S1_xyd_1}
\end{figure}
%
\begin{figure}[ht]
    \centering
    \includegraphics{S1_xyd_2.png}
    \caption{Flow depth after \SI{60}{s} in the simulation of Case S1.}
    \label{fig:S1_xyd_2}
\end{figure}
%
\begin{figure}[ht]
    \centering
    \includegraphics{S4_xyd.png}
    \caption{Flow depth after \SI{60}{s} in the simulation of Case S4.}
    \label{fig:S4_xyd}
\end{figure}

A closer look at the results shows that the whole domain present growth of oscillations since the start of the simulation. This is what we would expect as the numerical result of an ill-posed model. Yet, apparently, the upstream boundary introduces a perturbation larger than the ones initially in the domain and gives the effect that the solution is unperturbed until the upstream information has arrived to a certain location. Note also the fact that, if we see the consequence of ill-posedness is because the numerical scheme is introducing a minimum amount of diffusion (Section \ref{subsec:def_ill}). 
%
%We numerically test this result running a simulation in the conditions of the reference state considering secondary flow with no diffusion. The domain is \SI{200}{m} long and \SI{4}{m} wide and it is discretized is square cells of length \SI{0.2}{m}. The initial condition is a sinusoidal perturbation as in Simulations L (Section \ref{subsec:nosec}) with parameters $A=\SI{0.005}{m}$ and $L_b=\SI{20}{m}$. We have used a smaller grid size to gain accuracy and reduce numerical diffusion. The smaller grid size determines the domain length and initial condition wave length to run the simulations in a reasonable time. The time step is \SI{0.01}{s}.
%
%In Figure \ref{fig:sf_t} we show the flow depth along the streamwise direction at the side of the domain. The solution presents oscillations as we expect from the fact that it is ill-posed. 
%%
%\begin{figure}[ht]
    %\centering
    %\includegraphics{sf_t.pdf}
    %\caption{Flow depth along the streamwise direction at the side of the domain as a function of time.}
    %\label{fig:sf_t}
%\end{figure}
%
%
%----------------------------
\subsection{Bed Slope Effect}
%----------------------------
%
%
\label{subsec:bedslope_uni}
%
In this section we study the consequences of considering bed slope effects in the sediment transport rate. 

We gain empirical insight studying the reference state with the parameters of Table \ref{tab:L_ph} without considering secondary flow. The base case (B1, Table \ref{tab:B_var}) uses the bed slope function developed by \citet{Sekine92_1} (i.e., Equation (\ref{eq:g_sk}) with parameters $A_s=1$, $B_s=0$, $C_s=0$, and $D_s=0$). We numerically compute the eigenvalues of matrix $\mathbf{R}$, Equation \ref{eq:R}, for a varying wavenumber obtaining in this way the wave growth domain (Figure \ref{fig:dom_B1}). This result is in accordance with the linear theory of bar growth (Section \ref{subsec:nosec}). Note that it is well-posed since the growth rate decays for increasing wavenumber (Figure \ref{fig:dom_B1}a) and that there is a critical transversal wavelength below which no growth occurs (Figure \ref{fig:dom_B1}b). 
%
\begin{table}[ht]
%\tiny
	\begin{center}
	\begin{tabular}{ccc}		
			\toprule
Case & bed slope & Math. charac. \\
			\midrule
B1 & Y & W \\ 
B2 & N & I \\ 
 		\bottomrule
			\end{tabular}
      \caption{Parameters that are different between secondary flow cases. W=well-posed, I=ill-posed.}
      \label{tab:B_var}
	\end{center}
\end{table}
%
\begin{figure}[ht]
    \centering
    \includegraphics{dom_B1.png}
    \caption{Wave growth domain for the case B1 (Tables \ref{tab:L_ph} and \ref{tab:B_var}) as a function of the wavenumber (\textbf{a}) and wavelength (\textbf{b}). In the green area the growth rate is negative (dampening) while in the red area it is positive (growth). This situation is well-posed (Section \ref{subsec:def_ill}).}
    \label{fig:dom_B1}
\end{figure}

Now we analyze the same situation but without considering bed slope effects (Figure \ref{fig:dom_B2}). This case is ill-posed, which shows that bed slope effects are not only a refinement of the morphology module but a key element to obtain a physically realistic model.
%
\begin{figure}[ht]
    \centering
    \includegraphics{dom_B2.png}
    \caption{Wave growth domain for the case B2 (Tables \ref{tab:L_ph} and \ref{tab:B_var}) as a function of the wavenumber (\textbf{a}) and wavelength (\textbf{b}). In the green area the growth rate is negative (dampening) while in the red area it is positive (growth). This situation is ill-posed (Section \ref{subsec:def_ill}).}
    \label{fig:dom_B2}
\end{figure}

We test the results of the linear analysis as regards to bed slope effects running numerical simulations using Delft3D (Version 6.02.13.7545M compiled for Windows 32bit) in the conditions of cases B1 and B2. The domain is \SI{100}{m} long and \SI{10}{m} wide and it is discretized is square cells of length \SI{0.1}{m}. The initial condition is a flat unperturbed bed. The time step is \SI{0.1}{s}.

In Figures \ref{fig:B2_xyd} and \ref{fig:B2_xd} we plot the flow depth after \SI{8}{h} for the simulation of Case B2. Truncation errors in the initial condition have grown significantly and continue to grow due to ill-posedness of the system of equations. Figures \ref{fig:B1_xyd} and \ref{fig:B1_xd} present the flow depth after \SI{8}{h} for the simulation of Case B1. This case is well-posed and no oscillations are present. The initial condition has a slight deviation from the numerical uniform solution. This causes a small perturbation to travel, unamplified, in downstream direction (Figure \ref{fig:B1_xd}). 
%
\begin{figure}[ht]
    \centering
    \includegraphics{B2_xyd.png}
    \caption{Flow depth at the end of the simulation of Case B2.}
    \label{fig:B2_xyd}
\end{figure}
%
\begin{figure}[ht]
    \centering
    \includegraphics{B2_xd.eps}
    \caption{Longitudinal profile of the flow depth at the end of the simulation of Case B2.}
    \label{fig:B2_xd}
\end{figure}
%
\begin{figure}[ht]
    \centering
    \includegraphics{B1_xyd.png}
    \caption{Flow depth at the end of the simulation of Case B1.}
    \label{fig:B1_xyd}
\end{figure}
%
\begin{figure}[ht]
    \centering
    \includegraphics{B1_xd.eps}
    \caption{Longitudinal profile of the flow depth at the end of the simulation of Case B1.}
    \label{fig:B1_xd}
\end{figure}
%
%
%----------------------------
\subsection{Mixed-size Sediment}
%----------------------------
%
%
\label{subsec:mixsed}
%
In this section we study the effect of secondary flow and bed slope effects on the well-posedness of mixed-size sediment simulations. 

It is known that under certain conditions (mainly degradational conditions into a fine substrate) the 1D active layer model becomes ill-posed \citet{Ribberink87_PhD,Stecca14,Chavarrias18}. We consider a reference mixed-size sediment situation with the parameters shown in Table \ref{tab:ref_mixed}. We consider three situation varying the volume fraction content at the interface between the active layer and the substrate and the consideration of secondary flow (Table \ref{tab:mixed_var}). 
%
\begin{table}[ht]
%\tiny
	\begin{center}
	\begin{tabular}{cccccccc}		
			\toprule
$u$ [\si{m/s}] & $v$ [\si{m/s}] & $h$ [\si{m}] & $C_f$ [\si{-}] & $d_k$ [\si{m}] & $L_a$ [\si{m}] & $F_{ak}$ [$-$]\\
			\midrule
$1$ & $0$ & $1$ & $0.007$ & [0.001, 0.002] &  $0.1$ & $[0.50, 0.50]$ \\ 
 		\bottomrule
			\end{tabular}
      \caption{Physical parameters of the reference situation for mixed-size sediment conditions.}
      \label{tab:ref_mixed}
	\end{center}
\end{table}
%
\begin{table}[ht]
%\tiny
	\begin{center}
	\begin{tabular}{cccccccc}		
			\toprule
 Case & $f^I_{k}$ [$-$]  & Sec. flow & Math. charac. \\
			\midrule
M1 & $[0.50, 0.50]$ & NO  & W \\ 
M2 & $[1.00, 0.00]$ & NO  & I \\ 
M3 & $[1.00, 0.00]$ & YES & I \\ 
 		\bottomrule
			\end{tabular}
      \caption{Parameters that are different between mixed-size sediment cases. W=well-posed, I=ill-posed.}
      \label{tab:mixed_var}
	\end{center}
\end{table}

In case M1, the volume fraction content of the sediment at the interface between the active layer and the substrate is equal to the one in the active layer. This corresponds to aggradational conditions, or degradation into a substrate with the same grain size distribution as the active layer. This situation is well-posed in 1D \citep{Chavarrias18}. In Figure \ref{fig:dom_M1} we show the growth domain in 2D. Note that it is well-posed (Section \ref{subsec:def_ill}) and that we obtain a domain of linear bar growth under mixed-size sediment conditions. 
%
\begin{figure}[ht]
    \centering
    \includegraphics{dom_M1.png}
    \caption{Wave growth domain for the case M1 (Tables \ref{tab:ref_mixed} and \ref{tab:mixed_var}) as a function of the wavenumber (\textbf{a}) and wavelength (\textbf{b}). In the green area the growth rate is negative (dampening) while in the red area it is positive (growth). This situation is well-posed (Section \ref{subsec:def_ill}).}
    \label{fig:dom_M1}
\end{figure}
	
Case M2 differs from M1 in the volume fraction content of the sediment at the interface between the active layer and the substrate. In this case the substrate is finer than the active layer. In 1D this condition is ill-posed \citep{Chavarrias18}. In 2D we account for a diffusion mechanism (bed slope effects on the sediment transport direction) which is not present in the 1D analysis. Thus, it could be possible that in 2D the situation is well-posed. Figure \ref{fig:dom_M2} presents the growth domain of this situation. Note that the situation is ill-posed (Section \ref{subsec:def_ill}). We have tried several conditions with a coarser substrate (changing $f^I_k$) and more diffusion (changing $A_s$, Equation (\ref{eq:g_sk})) obtaining the same qualitative result. Thus, we conclude that the the diffusion due to bed slope effects in the sediment transport rate direction does not affect the conditions for ill-posedness due to the mixed-size sediment model. 
%
\begin{figure}[ht]
    \centering
    \includegraphics{dom_M2.png}
    \caption{Wave growth domain for the case M2 (Tables \ref{tab:ref_mixed} and \ref{tab:mixed_var}) as a function of the wavenumber (\textbf{a}) and wavelength (\textbf{b}). In the green area the growth rate is negative (dampening) while in the red area it is positive (growth). This situation is ill-posed (Section \ref{subsec:def_ill}).}
    \label{fig:dom_M2}
\end{figure}

This empirical result can be explained based on the diffusion matrices in Equations (\ref{eq:D_x}) and (\ref{eq:D_y}). In the conditions we study, the directional sediment transport can be expanded as:
\begin{equation}
\label{eq:qbkx_m1}
	q_{bkx}=q_{bk}\cos\left(\varphi_{sk}\right)\approx q_{bk} \;,
\end{equation}
and: 
\begin{equation}
\label{eq:qbky_m1}
	q_{bky}=q_{bk}\sin\left(\varphi_{sk}\right)\approx q_{bk} \varphi_{sk} \approx -q_{bk}\frac{1}{g_{sk}}\frac{\partial \eta}{\partial y} \;,
\end{equation}
where we have used that the angle between the sediment transport rate and the mean flow velocity is small such that:
\begin{equation}
\label{eq:varphisk_m1}
	\varphi_{sk}\approx \tan\left(\varphi_{sk}\right) = \frac{\sin{\varphi_{\tau}}-\frac{1}{g_{sk}}\frac{\partial \eta}{\partial y}}{\cos{\varphi_{\tau}}-\frac{1}{g_{sk}}\frac{\partial \eta}{\partial x}} \approx \frac{-\frac{1}{g_{sk}}\frac{\partial \eta}{\partial y}}{1-\frac{1}{g_{sk}}\frac{\partial \eta}{\partial x}} \approx -\frac{1}{g_{sk}}\frac{\partial \eta}{\partial y}\;,
\end{equation}
where we have neglected the effect of secondary flow and we have assumed a small bed slope in $x$ direction (i.e., $\frac{\partial \eta}{\partial x}<<1$). Then, the diffusive terms in the $x$ direction are equal to 0 and in $y$ direction are equal to:
\begin{equation}
\label{eq:dqbky_dsy_m1}
	\frac{\partial q_{bky}}{\partial \frac{\partial \eta}{\partial y}}=-\frac{q_{bk}}{g_{sk}}=-q_{bk} \;,
\end{equation}
where we have assumed the simplest formulation for the bed slope effects, i.e, $g_{sk}=1$ $\forall k$ \citep{Sekine92_1}. Then, the term $\frac{\partial q_{by}}{\partial \frac{\partial \eta}{\partial y}}$ in the \citeauthor{Exner20} equation is negative (diffusion) since it is the addition of two negative numbers. However, the term in the active layer equation:
\begin{equation}
\label{eq:dqbky_dsy_m2}
	\frac{\partial q_{bky}}{\partial \frac{\partial \eta}{\partial y}}-f^I_k\frac{\partial q_{by}}{\partial \frac{\partial \eta}{\partial y}}=-q_{bk}+f^I_k\sum_{l=1}^N q_{bl} \;,
\end{equation}
can be positive (antidiffusion) and it is positive in degradational conditions into a fine substrate. This is easy to see in the simplest case of two size fractions ($N=2$):
\begin{equation}
\label{eq:dqbky_dsy_m3}
	\frac{\partial q_{bk1}}{\partial \frac{\partial \eta}{\partial y}}-f^I_1\frac{\partial q_{by}}{\partial \frac{\partial \eta}{\partial y}}=-q_{b1}\left(1-f^I_1\right)+f^I_1q_{b2} \;.
\end{equation}
Under aggradational conditions ($f^I_k$=$F_{ak}$) we always find diffusion:
\begin{equation}
\label{eq:dqbky_dsy_m3_agg}
	\frac{\partial q_{bk1}}{\partial \frac{\partial \eta}{\partial y}}-F_{a1}\frac{\partial q_{by}}{\partial \frac{\partial \eta}{\partial y}}=F_{a1}F_{a2}\left(Q_{b2}-Q_{b1}\right)<0 \;,
\end{equation}
assuming that reverse mobility (Section \ref{subsubsec:sedtrans}) does not occur. 

The above derivation explains the results of the linear analysis but please note that it is not formal proof that the ill-posedness due to mixed-size sediment is not affected by bed slope effects in the sediment transport rate. 

We have seen that for a unisize case it is necessary to account for bed slope effects in order to obtain a well-posed model. Under mixed-size sediment conditions we find a more restrictive condition because not all closure relations for bed slope effects guarantee that the diffusion terms is negative for all active layer equations. 

In Case M3 we test the effect of the diffusive term in the constitutive equation for the secondary flow on the ill-posedness due to mixed-size sediment. Case M3 is equal to Case M2 but we consider secondary flow with a relatively large diffusive coefficient $D_H=\SI{10}{m^2/s}$ (Figure \ref{fig:dom_M3}). There is hardly any difference when compared to Case M2 which provides evidence that the ill-posedness due to the mixed-size sediment model is independent of the flow. 
%
\begin{figure}[ht]
    \centering
    \includegraphics{dom_M3.png}
    \caption{Wave growth domain for the case M3 (Tables \ref{tab:ref_mixed} and \ref{tab:mixed_var}) as a function of the wavenumber (\textbf{a}) and wavelength (\textbf{b}). In the green area the growth rate is negative (dampening) while in the red area it is positive (growth). This situation is ill-posed (Section \ref{subsec:def_ill}).}
    \label{fig:dom_M3}
\end{figure}
%
%----------------------------
%
\section{Ill-posedness Check in Delft3D}
%
%----------------------------
%
\label{sec:ch_d3d}
%
We aim at implementing a routine that determines whether a numerical simulation run in Delft3D is under condition in which the model is ill-posed. In this section we explain the implementation of the routine. In Section \ref{subsec:desc} we describe the routine. Sections \ref{subsec:mat_test}, \ref{subsec:eig_test}, and \ref{subsec:out_test} present the test of the matrices, eigenvalues, and output. In Section \ref{subsec:time_test} we test the computational time added in conducting the ill-posedness test.
%
%
%----------------------------
\subsection{Description}
%----------------------------
%
%
\label{subsec:desc}
%
In this section we describe the routine implemented in Delft3D to check for ill-posedness. The routine is implemented on a branch (\texttt{/branches/research/Technical University of Delft/20161020\_ellipticity\_check}) created from the trunk of Delft3D at revision 6118 (20 October 2016 13:17:24).

The mathematical character (i.e., well-posed or ill-posed) is a property of each point in the space-time domain. Thus, the routine should check every node at every time step. Since we expect this to be reasonably expensive, the routine checks all nodes at the time steps before saving map results. The user can decrease the interval between map results to increase the time resolution of the ill-posed check output. The output is a variable defined for every node which is equal to 0 if the model is well-posed an equal to 1 if it is ill-posed. 

The results of the linear analysis show that a model can be ill-posed due to a lack of diffusion in the secondary flow, the lack of a bed slope effects function in unisize sediment conditions, the use of a wrong closure relation for bed slope effects under mixed-sediment conditions, and the grain size distribution conditions (mainly degradation into a substrate finer than the active layer). To test whether the problem at a certain node and at a certain time is ill-posed, we should build the matrix $\mathbf{R}$, Equation (\ref{eq:R}), and check that for the largest wavenumbers the growth rate is negative. This is, fixing the wavenumber in $x$ direction ($k_{wx}$) checking that for all $k_{wy}$ the real part of all eigenvalues is negative. Idem fixing the wavenumber in $y$ direction. 

The above approach is the most general one but it is also relatively expensive since it involves an eigenvalue problem for several values of the wavenumber in $x$ and $y$ direction. Usual values of the diffusion coefficient used in numerical simulations are reasonably large \citep{Madsen88}. Moreover, to check that the problem is well-posed for bed slope effects under mixed-size sediment conditions suffices to check that the diffusion coefficient is negative. Thus, we first focus on the problem of ill-posedness due to the grain size distribution conditions. In this case, we have empirical results that show that ill-posedness is independent on the diffusion coefficient and bed slope effect. Neglecting the diffusive matrices we obtain that the result is only dependent on matrices $\mathbf{A_{x}}$ and $\mathbf{A_{y}}$. The results also show that if ill-posedness does not depend in this cases on the wavenumbers, meaning that, if it is ill-posed, the growth rate is positive for all wavenumbers. Using these results we implement a routine that builds matrices $\mathbf{A_{x}}$ and $\mathbf{A_{y}}$ for all nodes and computes the eigenvalues of these matrices. If all eigenvalues of both matrices are real, the model is well-posed. If one single eigenvalue of one of the two matrices is complex, the problem is ill-posed. 

\citet{Chavarrias16} checked that the matrices were correctly implemented. The result is written in Section \ref{subsec:mat_test} for the sake of completeness. We test that the eigenvalues are correctly computed (Section \ref{subsec:eig_test}) and that the output is correctly passed (Section \ref{subsec:out_test}). 

We use the same three test simulations as \citet{Chavarrias16}. The simulations are 2 dimensional. The domain is rectangular with two open boundaries on opposite ends. The initial condition is in equilibrium for a certain upstream water and sediment discharge and a certain downstream water level. The initial water and sediment discharge remain constant throughout the simulation time. A slow lowering (\SI{0.001}{m/min}) of the downstream water level is imposed that causes degradational conditions. The simulations differ in the number of size fractions and the initial stratigraphy. Table \ref{tab:I_ph_eq} summarizes the initial conditions and physical parameters that are the same for all simulations. Table \ref{tab:I_ph_dif} contains the initial conditions and physical parameters which are different between simulations. Table \ref{tab:I_num} summarizes the numerical parameters of the simulations. 

\begin{table}[ht]
%\tiny
	\begin{center}
	\begin{tabular}{cccccccc}		
			\toprule
       $L$ [\si{m}] & $B$ [\si{m}] & $u$ [\si{m/s}] & $v$ [\si{m/s}] & $h$ [\si{m}] & $C_f$ [\si{-}] & $L_a$ [\si{m}] & $s$ [$-$] \\
			\midrule
$1$ & $1$ & $1$ & $0$ & $1$ & $0.007$ & $0.1$ & $7.13541\cdot10^{-4}$ \\ 
			\bottomrule
			\end{tabular}
      \caption{Physical parameters of the simulations to test the implementation that are equal for all simulations}
      \label{tab:I_ph_eq}
	\end{center}
\end{table}

\begin{table}[ht]
\tiny
	\begin{center}
	\begin{tabular}{ccccccc}		
			\toprule
      Simulation & $d_k$ [\si{m}] & $F_{ak}$ [$-$] & $f^I_{k}$ [$-$] & $q_{bk}$ [\si{m^2/s}] & Sed. trans. & Math. charac. \\
			\midrule
I1  & $[0.001, 0.002]$					& $[0.50, 0.50]$ 				& $[0.50, 0.50]$ 			 & $[5.588, 2.794]\cdot10^{-5}$ 				& EH 				& W \\ 
I2  & $[0.001, 0.002]$ 					& $[0.00, 1.00]$ 				& $[1.00, 0.00]$ 			 & $[0.000, 5.588]\cdot10^{-5}$ 				& EH 				& I \\ 
I3  & $[0.001, 0.002, 0.004]$ 	& $[0.33, 0.33, 0.33]$ 	& $[0.33, 0.33, 0.33]$ & $[6.723, 6.319, 5.865]\cdot10^{-5}$ 	& MPM-P0.8  & W \\ 
			\bottomrule
			\end{tabular}
      \caption{Physical parameters of the simulations to test the implementation that are different for all simulations. EH stands for \citet{Engelund67} and MPM-P0.8 stands for \citet{MeyerPeter48} with the hiding correction by \citet{Parker82_3} with parameter $b$=0.8. W means well-posed and I ill-posed.}
      \label{tab:I_ph_dif}
	\end{center}
\end{table}

\begin{table}[ht]
%\tiny
	\begin{center}
	\begin{tabular}{ccccccccc}		
			\toprule
$\Delta t$ [\si{s}] & $\Delta x$ [\si{m}] & $\Delta y$ [\si{m}] & $\Delta z$ [\si{m}] & $T$ [\si{s}]  \\
			\midrule
0.1 & 0.5 & 0.5 & 0.5 & 60  \\ 
			\bottomrule
			\end{tabular}
      \caption{Numerical parameters of the simulations to test the implementation}
      \label{tab:I_num}
	\end{center}
\end{table}

%The symbols not defined previously are the length $L$ [\si{m}] and width $B$ [\si{m}] of the domain, the bed slope $s$ [$-$], the time step $\Delta t$ [\si{s}], the space step in the $x$ direction $\Delta x$ [\si{m}], the space step in $y$ direction [\si{m}], the thickness of the substrate layers $\Delta z$ [\si{m}], and the simulation time $T$ [\si{s}]. 
%
%
%----------------------------
\subsection{Matrices Test}
%----------------------------
%
%
\label{subsec:mat_test}
%
In this section we test the implementation of matrices $\mathbf{A_x}$ and $\mathbf{A_y}$. This was done by \citet{Chavarrias16} and it is presented ere for the sake of completeness. The subscript $_\text{Delft3D}$ stands for the matrices computed by the routine implemented in the Delft3D branch. The matrices without subscript are computed analytically using the Matlab \textregistered~ Symbolic Toolbox.

\begin{equation}
\label{eq:Ax_S1}
\makebox[\textwidth][c]{$
		\mathbf{A}_{xI1}=\left[
 \begin{array}{rrrrr}
 0 									&	 1										&  0	&  0				&  0 									\\
 8.8100 						&	 2										&  0	&  9.8100 	&	 0 									\\
 0									&  0										&  1	&  0				&  0 									\\
-6.9853\cdot10^{-4} &  6.9853\cdot10^{-4} 	&  0 	&  0 				&  9.3138\cdot10^{-4}	\\
-1.1642\cdot10^{-4} &  1.1642\cdot10^{-4} 	&  0 	&  0 				&  1.3970\cdot10^{-3}	\\	
 \end{array}\right]
$}
\end{equation}

\begin{equation}
\label{eq:Ay_S1}
\makebox[\textwidth][c]{$
		\mathbf{A}_{yI1}=\left[
 \begin{array}{rrrrr}
 0 			&	 0	&  1										&  0			& 0 \\
 0 			&	 0	&  1										&  0 			&	0 \\
 9.8100	&  0	&  0										&  9.8100	& 0 \\
 0 			&  0 	&  1.3970\cdot10^{-4} 	&  0 			& 0 \\
 0 			&  0 	&  2.3284\cdot10^{-5} 	&  0 			& 0 \\	
 \end{array}\right]
$}
\end{equation}

\begin{equation}
	\label{eq:Ax_S1_Delft3D}
	\makebox[\textwidth][c]{$
		\mathbf{A}_{xI1,\text{Delft3D}}=\left[
		\begin{array}{rrrrr}
 0  								&  1.0000  						&  0  								& 	0  			&  0  \\
 8.8097  						&  2.0001  						&  0  								& 	9.8098	&  0  \\
 0  								&  0  								&  1.0001  						& 	0  			&  0  \\
-6.9828\cdot10^{-4} &  6.9859\cdot10^{-4}	&  2.7937\cdot10^{-8}	& 	0  			&  9.3135\cdot10^{-4}  \\
-1.1639\cdot10^{-4} &  1.1644\cdot10^{-4}	&  4.6564\cdot10^{-9}	& 	0  			&  1.3970\cdot10^{-3} \\
         \end{array}\right]
		$}
\end{equation}

\begin{equation}
	\label{eq:Ay_S1_Delft3D}
	\makebox[\textwidth][c]{$
		\mathbf{A}_{yI1,\text{Delft3D}}=\left[
		\begin{array}{rrrrr}
 0   			&  0 &  1.0000  						&  0    		&  0  \\
 0   			&  0 &  1.0001  						&  0    		&  0  \\
 9.8098   &  0 &  0  									&  9.8098   &  0  \\
 0   			&  0 &  1.3969\cdot10^{-4}  &  0    		&  0  \\
 0   			&  0 &  2.3283\cdot10^{-5}  &  0    		&  0  \\
		\end{array}\right]
		$}
\end{equation}

\begin{equation}
\label{eq:Ax_S2}
\makebox[\textwidth][c]{$
		\mathbf{A}_{xI2}=\left[
 \begin{array}{rrrrr}
 0 									&	 1									& 0	& 0				& 0 \\
 8.8100 						&	 2									& 0	& 9.8100 	&	0 \\
 0									&  0									& 1	& 0				& 0 \\
-4.6569\cdot10^{-4} &  4.6569\cdot10^{-4} & 0 & 0 			& 9.3138\cdot10^{-4} \\
 4.6569\cdot10^{-4} & -4.6569\cdot10^{-4} & 0 & 0 			& 9.3138\cdot10^{-4} \\	
 \end{array}\right]
$}
\end{equation}

\begin{equation}
\label{eq:Ay_S2}
\makebox[\textwidth][c]{$
		\mathbf{A}_{yI2}=\left[
 \begin{array}{rrrrr}
0 			&	0	&  1									& 0				& 0 \\
0 			&	0	&  1									& 0 			&	0 \\
9.8100	& 0	&  0									& 9.8100	& 0 \\
0 			& 0 &  9.3138\cdot10^{-5} & 0 			& 0 \\
0 			& 0 & -9.3138\cdot10^{-5} & 0 			& 0 \\	
 \end{array}\right]
$}
\end{equation}

\begin{equation}
	\label{eq:Ax_S2_Delft3D}
	\makebox[\textwidth][c]{$
		\mathbf{A}_{xI2,\text{Delft3D}}=\left[
		\begin{array}{rrrrr}
   0 									&  1.0000  						&  0  								&  0  			& 0 \\
   8.8097 						&  2.0001  						&  0  								&  9.8098  	& 0 \\
   0 									&  0  								&  1.0001  						&  0  			& 0 \\
  -4.6426\cdot10^{-4} &  4.6655\cdot10^{-4}	&  1.8624\cdot10^{-7}	&  0  			& 9.3135\cdot10^{-4} \\
   4.6426\cdot10^{-4} & -4.6655\cdot10^{-4}	& -1.8624\cdot10^{-7}	&  0  			& 9.3135\cdot10^{-4}  \\
         \end{array}\right]
		$}
\end{equation}
					
\begin{equation}
	\label{eq:Ay_S2_Delft3D}
	\makebox[\textwidth][c]{$
		\mathbf{A}_{yI2,\text{Delft3D}}=\left[
		\begin{array}{rrrrr}
0     	&  0 &  1.0000 							&  0  			&  0   \\
0     	&  0 &  1.0001 							&  0  			&  0   \\
9.8098	&  0 &  0 									&  9.8098  	&  0   \\
0     	&  0 &  9.3125\cdot10^{-5} 	&  0  			&  0   \\
0     	&  0 & -9.3125\cdot10^{-5} 	&  0  			&  0   \\
		\end{array}\right]
		$}
\end{equation}

\begin{equation}
\label{eq:Ax_S3}
\makebox[\textwidth][c]{$
		\mathbf{A}_{xI3}=\left[
 \begin{array}{rrrrrr}
0 									&	 1									 &  0	&  0			& 0 								 &  0  								\\
8.8100 							&	 2									 &  0	&  9.8100 &	0 								 &  0									\\
0										&  0									 &  1	&  0			& 0 								 &  0									\\
-1.2553\cdot10^{-3} &  1.2553\cdot10^{-3}  &  0 &  0 			& 2.0233\cdot10^{-3} &  1.2898\cdot10^{-3} \\
-9.2609\cdot10^{-6} &  9.2609\cdot10^{-6}  &  0 &  0 			& 3.1578\cdot10^{-3} & -1.1604\cdot10^{-4} \\
-5.1076\cdot10^{-7} &  5.1076\cdot10^{-7}  &  0	&  0 			&-1.4466\cdot10^{-4} &  3.0825\cdot10^{-3} \\
 \end{array}\right]
$}
\end{equation}

\begin{equation}
\label{eq:Ay_S3}
\makebox[\textwidth][c]{$
		\mathbf{A}_{yI3}=\left[
 \begin{array}{rrrrrr}
0 			&	0	& 1									 & 0		  & 0 & 0\\
0 			&	0	& 1									 & 0 			&	0 & 0\\
9.8100	& 0	& 0									 & 9.8100	& 0 & 0\\
0 			& 0 & 3.1510\cdot10^{-4} & 0 			& 0 & 0\\
0 			& 0 & 7.0133\cdot10^{-6} & 0 		  & 0 & 0\\
0 			& 0 & 2.7599\cdot10^{-7} & 0 		  & 0 & 0\\	
 \end{array}\right]
$}
\end{equation}

\begin{equation}
	\label{eq:Ax_S3_Delft3D}
	\makebox[\textwidth][c]{$
		\mathbf{A}_{xI3,\text{Delft3D}}=\left[
 \begin{array}{rrrrrr}
0 									&  1.0000 						&  0 									&  0  			&  0 									&  0  \\
8.8097 							&  2.0001 						&  0 									&  9.8098  	&  0 									&  0  \\
0 									&  0 									&  1.0001 						&  0  			&  0 									&  0  \\
-1.2549\cdot10^{-3} &  1.2553\cdot10^{-3} &  4.7001\cdot10^{-8} &  0  			&  2.0232\cdot10^{-3} &  1.2898\cdot10^{-3}  \\
-9.2611\cdot10^{-6} &  9.2605\cdot10^{-6} &  1.1237\cdot10^{-10}&  0  			&  3.1590\cdot10^{-3} & -1.1600\cdot10^{-4}  \\
-5.1066\cdot10^{-7} &  5.1063\cdot10^{-7} &  1.1737\cdot10^{-11}&  0  			& -1.4453\cdot10^{-4} &  3.0833\cdot10^{-3}  \\
         \end{array}\right]
		$}
\end{equation}

\begin{equation}
	\label{eq:Ay_S3_Delft3D}
	\makebox[\textwidth][c]{$
		\mathbf{A}_{yI3,\text{Delft3D}}=\left[
		\begin{array}{rrrrrr}
0 			& 0 &  1.0000  						& 0   		& 0 &  0   \\
0 			& 0 &  1.0001  						& 0   		& 0 &  0   \\
9.8098 	& 0 &  0  								& 9.8098 	& 0 &  0   \\
0 			& 0 &  3.1505\cdot10^{-4}	& 0   		& 0 &  0   \\
0 			& 0 &  7.0130\cdot10^{-6}	& 0   		& 0 &  0   \\
0 			& 0 &  2.7588\cdot10^{-7}	& 0   		& 0 &  0   \\
		\end{array}\right]
		$}
\end{equation}

We conclude that the matrices are correctly implemented. 
%
%
%----------------------------
\subsection{Eigenvalues Test}
%----------------------------
%
%
\label{subsec:eig_test}
%
In this section we test that the eigenvalues are correctly computed. In Delft3D, we compute eigenvalues using the double-precision Fortran subroutine 'rg' of the EISPACK collection. We compare the outcome to the eigenvalues computed using the Matlab \textregistered~ Symbolic Toolbox.

In Tables \ref{tab:I_mat_x} and \ref{tab:I_mat_y} we show the eigenvalues of matrices $\mathbf{A_x}$ and $\mathbf{A_y}$, respectively, as computed for the reference situation using Matlab. In Tables \ref{tab:I_d3d_x} and \ref{tab:I_d3d_y} we show the eigenvalues given in Delft3D at the node $n=2$ and $m=2$ at $t=2$ $s$. The minimal difference is mainly due to the fact that while the eigenvalues computed in Matlab are those at the exact reference situation, the ones in Delft3D are computed with the flow conditions after \SI{2}{s}. 

\begin{table}[ht]
%\tiny
	\begin{center}
	\begin{tabular}{
	l
	S[retain-explicit-plus,table-align-uncertainty=false,table-align-exponent=false,table-format=1.3e1]
	S[retain-explicit-plus,table-align-uncertainty=false,table-align-exponent=false]
	S[retain-explicit-plus,table-align-uncertainty=false,table-align-exponent=false]
	}		
			\toprule
Eigenvalue & I1 & I2 & I3 \\
			\midrule
$\lambda_1$ & \num{+4.133 +0.000 i e0 } & \num{+4.133 +0.000 i e0 } & \num{+4.134 +0.000 i e0 } \\ 
$\lambda_2$ & \num{-2.134 +0.000 i e0 } & \num{-2.133 +0.000 i e0 } & \num{-2.135 +0.000 i e0 } \\ 
$\lambda_3$ & \num{+6.216 +0.000 i e-4} & \num{+7.251 +6.630 i e-4} & \num{+1.382 +0.000 i e-3} \\ 
$\lambda_4$ & \num{+1.552 +0.000 i e-3} & \num{+7.251 -6.630 i e-4} & \num{+3.258 +0.000 i e-3} \\ 
$\lambda_5$ & \num{+1.000 +0.000 i e0 } & \num{+1.000 +0.000 i e0 } & \num{+2.995 +0.000 i e-3} \\ 
$\lambda_6$ & 													& \num{+0.000 +0.000 i e0 } & \num{+1.000 +0.000 i e0 } \\ 
			\bottomrule
			\end{tabular}
      \caption{Eigenvalues in the $x$ direction for the reference situation computed with Matlab.}
      \label{tab:I_mat_x}
	\end{center}
\end{table}

\begin{table}[!ht]
%\tiny
	\begin{center}
	\begin{tabular}{
	l
	S[retain-explicit-plus,table-align-uncertainty=false,table-align-exponent=false,table-format=1.3e1]
	S[retain-explicit-plus,table-align-uncertainty=false,table-align-exponent=false]
	S[retain-explicit-plus,table-align-uncertainty=false,table-align-exponent=false]
	}		
			\toprule
Eigenvalue & I1 & I2 & I3 \\
			\midrule
$\lambda_1$ & \num{+3.132 +0.000 i e0} & \num{+3.132 +0.000 i e0} & \num{-3.133 +0.000 i e0} \\ 
$\lambda_2$ & \num{-3.132 +0.000 i e0} & \num{-3.132 +0.000 i e0} & \num{+3.133 +0.000 i e0} \\ 
$\lambda_3$ & \num{+0.000 +0.000 i e0} & \num{+0.000 +0.000 i e0} & \num{+0.000 +0.000 i e0} \\ 
$\lambda_4$ & \num{+0.000 +0.000 i e0} & \num{+0.000 +0.000 i e0} & \num{+0.000 +0.000 i e0} \\ 
$\lambda_5$ & \num{+0.000 +0.000 i e0} & \num{+0.000 +0.000 i e0} & \num{+0.000 +0.000 i e0} \\ 
$\lambda_6$ & \num{+0.000 +0.000 i e0} & \num{+0.000 +0.000 i e0} & \num{+0.000 +0.000 i e0} \\ 
			\bottomrule
			\end{tabular}
      \caption{Eigenvalues in the $y$ direction for the reference situation computed with Matlab.}
      \label{tab:I_mat_y}
	\end{center}
\end{table}

\begin{table}[!ht]
%\tiny
	\begin{center}
	\begin{tabular}{
	l
	S[retain-explicit-plus,table-align-uncertainty=false,table-align-exponent=false,table-format=1.3e1]
	S[retain-explicit-plus,table-align-uncertainty=false,table-align-exponent=false]
	S[retain-explicit-plus,table-align-uncertainty=false,table-align-exponent=false]
	}		
			\toprule
Eigenvalue & I1 & I2 & I3 \\
			\midrule
$\lambda_1$ & \num{+4.133 +0.000 i e0 } & \num{+4.133 +0.000 i e0 } & \num{+4.133 +0.000 i e0 } \\ 
$\lambda_2$ & \num{-2.134 +0.000 i e0 } & \num{-2.133 +0.000 i e0 } & \num{-2.135 +0.000 i e0 } \\ 
$\lambda_3$ & \num{+6.216 +0.000 i e-4} & \num{+7.249 +6.622 i e-4} & \num{+1.383 +0.000 i e-3} \\ 
$\lambda_4$ & \num{+1.552 +0.000 i e-3} & \num{+7.249 -6.622 i e-4} & \num{+3.261 +0.000 i e-3} \\ 
$\lambda_5$ & \num{+1.000 +0.000 i e0 } & \num{+1.000 +0.000 i e0 } & \num{+2.998 +0.000 i e-3} \\ 
$\lambda_6$ & 													& \num{+0.000 +0.000 i e0 } & \num{+1.000 +0.000 i e0 } \\ 
			\bottomrule
			\end{tabular}
      \caption{Eigenvalues in the $x$ direction computed in Delft3D.}
      \label{tab:I_d3d_x}
	\end{center}
\end{table} 

\begin{table}[!ht]
%\tiny
	\begin{center}
	\begin{tabular}{
	l
	S[retain-explicit-plus,table-align-uncertainty=false,table-align-exponent=false,table-format=1.3e1]
	S[retain-explicit-plus,table-align-uncertainty=false,table-align-exponent=false]
	S[retain-explicit-plus,table-align-uncertainty=false,table-align-exponent=false]
	}		
			\toprule
Eigenvalue & I1 & I2 & I3 \\
			\midrule
$\lambda_1$ & \num{+3.132 +0.000 i e0} & \num{+3.132 +0.000 i e0} & \num{-3.132 +0.000 i e0} \\ 
$\lambda_2$ & \num{-3.132 +0.000 i e0} & \num{-3.132 +0.000 i e0} & \num{+3.132 +0.000 i e0} \\ 
$\lambda_3$ & \num{+0.000 +0.000 i e0} & \num{+0.000 +0.000 i e0} & \num{+0.000 +0.000 i e0} \\ 
$\lambda_4$ & \num{+0.000 +0.000 i e0} & \num{+0.000 +0.000 i e0} & \num{+0.000 +0.000 i e0} \\ 
$\lambda_5$ & \num{+0.000 +0.000 i e0} & \num{+0.000 +0.000 i e0} & \num{+0.000 +0.000 i e0} \\ 
$\lambda_6$ & \num{+0.000 +0.000 i e0} & \num{+0.000 +0.000 i e0} & \num{+0.000 +0.000 i e0} \\ 
			\bottomrule
			\end{tabular}
      \caption{Eigenvalues in the $y$ direction computed in Delft3D.}
      \label{tab:I_d3d_y}
	\end{center}
\end{table}

We conclude that the eigenvalues are correctly computed. 
%
%
%----------------------------
\subsection{Output Test}
%----------------------------
%
%
\label{subsec:out_test}
%
In this section we test that the output of the ill-posedness check routine is correct. 

A critical point is the fact that the volume fraction content of sediment at the interface between the active layer and the substrate is input to the routine. This volume fraction content is taken as the one of the first substrate layer, i.e., the first layer below the active layer. Due to the manner in which the substrate is handled in Delft3D, a step before the check is to identify which is the first substrate layer, i.e., the first layer in the substrate with a thickness different than 0. The simple cases used before do not allow to test the implementation of this part. To this end, we run a fourth case (I4) with the same parameters as I1. However, the width is reduced to \SI{0.1}{m} and the domain is discretized in cells \SI{0.1}{m} long and \SI{0.01}{m} wide. The width reduction is necessary to preserve the one-dimensional character when the simulation becomes ill-posed. Otherwise degradation is not uniform along the cross-section which complicates the analysis. 

The substrate is formed out of \SI{1}{cm} thick layers alternating coarse and fine material (Figure \ref{fig:I4_fsk}). A \SI{5}{cm} lowering of the downstream water level in \SI{5}{h} causes the entrainment of the substrate sediment. The system is well-posed when degradation of the coarse layers occur and it is ill-posed when the layers composed out of fine sediment are entrained (Figure \ref{fig:I4_i}). 
%
\begin{figure}[ht]
    \centering
    \includegraphics{I4_fsk.png}
    \caption{Longitudinal profile at the center of the domain ($y$=\SI{5}{m}) at the initial state of the simulation of Case I4.}
    \label{fig:I4_fsk}
\end{figure}
%
\begin{figure}[ht]
    \centering
    \includegraphics{I4_i.png}
    \caption{Volume fraction content at the first substrate layer (\textbf{a}) and output of ill-posedness (\textbf{b}). The data concerns the center of the domain ($y$=\SI{5}{m}).}
    \label{fig:I4_i}
\end{figure}

At the downstream end of the domain there are four short periods of time for which the routine's output is ``well-posed'' while the first substrate layer is made out of fine sediment. This is because, since the volume fraction content of sediment at the interface between the active layer and the substrate is not an output of Delft3D, we post-process the Delft3D output to find the first substrate layer, but we do know whether the bed was aggrading or degrading. We have checked the bed elevation with time and we confirm that at the locations where the first substrate layer is made out of fine sediment but the output is ``well-posed'' aggradation occurs. The origin of the localized aggradational conditions is the fact that the model is ill-posed and generates oscillations (Section \ref{subsec:def_ill}). 

We conclude that the output is correct for all $M$ and $N$ nodes and for all time. 
%
%
%----------------------------
\subsection{Computational Time Test}
%----------------------------
%
%
\label{subsec:time_test}
%
In this section we study the added computational cost of the ill-posedness routine. 

We consider a case with the parameters in Table \ref{tab:T_ph1}. We use the sediment transport relation developed by \citet{Engelund67} with the bed slope effects closure relation developed by \citet{Sekine92_1}. The initial condition is in equilibrium. At the upstream end, the bed level and the composition are fixed. The downstream water level raises \SI{6}{mm} in \SI{600}{s}.

We run 10 simulations in which we test the cost of the ill-posedness check routine compared to the cost of computing secondary flow (Table \ref{tab:T_var}). Since the ill-posedness check routine is executed before results are saved, we also test the cost of increasing the results time. Since the size of the matrices increases with the power 2 of the number of sediment size fractions we test the consequence of increasing the number of fractions. We compare the latest tagged version (6.02.13.7545M) to the branch in which we implemented the ill-posedness check routine (6.02.07.7900M). 

For the case $N=2$ the characteristic grain sizes are \SI{1}{mm} and \SI{2}{mm} and the initial volume fraction content of sediment in the active layer and in the substrate is equal to 0.5. For the case $N=10$ the characteristic grain sizes are linearly interpolated between \SI{1}{mm} and \SI{2}{mm}. The volume fraction content is equal to 0.1 for all size fractions both at the initial active layer and the substrate. The numerical parameters are shown in Table \ref{tab:T_ph2}. 

Both Delft3D versions were compiled with Visual Fortran Compiler \textregistered 12.1 using Visual Studio \textregistered 2010 for Windows 64. Simulations were computed on an Intel \textregistered Xeon \textregistered CPU E5-1650 v3 at \SI{3.50}{GHz} processor with \SI{32.0}{GB} of RAM on a Windows 7 Enterprise operating system. Simulations run one at a time and no other expensive process (i.e., other than text editing) was in use while the simulations run. 
%
\begin{table}[ht]
%\tiny
	\begin{center}
	\begin{tabular}{cccccccc}		
			\toprule
$L$ [\si{m}] & $B$ [\si{m}] & $u$ [\si{m/s}] & $v$ [\si{m/s}] & $h$ [\si{m}] & $C_f$ [\si{-}] & $L_a$ [\si{m}] & $s$ [$-$] \\
			\midrule
$10$ & $10$ & $1$ & $0$ & $1$ & $0.007$ & $0.1$ & $7.13541\cdot10^{-4}$ \\ 
			\bottomrule
			\end{tabular}
      \caption{Phisical parameters of the simulations to test the computational cost.}
      \label{tab:T_ph1}
	\end{center}
\end{table}

\begin{table}[ht]
%\tiny
	\begin{center}
	\begin{tabular}{cccccccc}		
			\toprule
$\Delta t$ [\si{s}] & $\Delta x$ [\si{m}] & $\Delta y$ [\si{m}] & $\Delta z$ [\si{m}] & $z_{tot}$ [\si{m}] & $T$ [\si{s}]  \\
			\midrule
0.1 & 0.1 & 0.1 & 0.5 & 5 & 600  \\ 
			\bottomrule
			\end{tabular}
      \caption{Numerical parameters of the simulations to test the computational cost.}
      \label{tab:T_ph2}
	\end{center}
\end{table}

\begin{table}[ht]
%\tiny
	\begin{center}
	\begin{tabular}{cccccc}		
			\toprule
      Simulation & Version & Sec. flow & ill-p. check & $N$ & $T/\Delta t_{\text{r}}$ \\
			\midrule
T1   & tag	    & Y	& - & 2  & 10 \\ 
T2   & tag	    & N	& - & 2  & 10 \\ 
T3   & branch	  & Y	& N & 2  & 10 \\ 
T4   & branch	  & Y	& Y & 2  & 10 \\ 
T5   & branch	  & N	& N & 2  & 10 \\ 
T6   & branch	  & N	& Y & 2  & 10 \\ 
T7   & branch	  & N	& N & 2  & 100 \\ 
T8   & branch	  & N	& Y & 2  & 100 \\ 	
T9   & branch	  & N	& N & 10 & 10 \\ 
T10  & branch	  & N	& Y & 10 & 10 \\ 		
			\bottomrule
			\end{tabular}
      \caption{Simulations to test the computation cost of the ill-posedness check routine.}
      \label{tab:T_var}
	\end{center}
\end{table}

The results are shown in Table \ref{tab:T_times} and in Figure \ref{fig:T}. The time needed to check for ill-posedness is insignificant (compare T3-T4, T5-T6, T7-T8, T9-T10). It is specially insignificant when comparing to the added cost of computing the transport of a constituent such as secondary flow (compare T3-T6). The computation of eigenvalues does not significantly increase with an increase of the number of sediment size fractions (compare T9-T10). Ten times more results (and times we check for ill-posedness) does not significantly increase the overall computational time (compare T6-T8). Note that in this difference in time there is also the time of actually writing the results (compare T5-T7). Without secondary flow the branch is faster than the latest tagged version (compare T1-T3). The opposite happens when solving secondary flow (compare T2-T5).
\begin{table}[ht]
%\tiny
	\begin{center}
	\begin{tabular}{rrrrrr|r}		
			\toprule
Simulation & Mom. eq. [\si{s}] & Cont. eq. [\si{s}] & Trans. eq. [\si{s}] & Turb. [\si{s}] & 3D Morph. [\si{s}] & \textbf{TOTAL} [\si{s}] \\
			\midrule
T1  &   54.63  & 61.39  & 16.71 &  0.98 &  461.40 &\textbf{  595.11}  \\
T2  &   28.17  & 34.05  &  0.00 &  1.01 &  399.26 &\textbf{  462.49}  \\
T3  &   50.81  & 57.33  & 15.50 &  0.89 &  430.66 &\textbf{  555.19}  \\
T4  &   51.01  & 57.11  & 15.63 &  0.92 &  431.74 &\textbf{  556.41}  \\
T5  &   27.71  & 33.80  &  0.00 &  1.09 &  407.11 &\textbf{  469.71}  \\
T6  &   27.79  & 33.93  &  0.00 &  0.97 &  411.10 &\textbf{  473.79}  \\ 
T7  &   27.90  & 33.99  &  0.00 &  1.04 &  409.47 &\textbf{  472.30}  \\
T8  &   27.68  & 33.58  &  0.00 &  0.95 &  417.99 &\textbf{  480.10}  \\
T9  &   27.91  & 33.72  &  0.00 &  0.98 & 1599.59 &\textbf{ 1662.10}  \\
T10 &   27.66  & 33.81  &  0.00 &  0.84 & 1601.12 &\textbf{ 1663.43}  \\   
			\bottomrule
			\end{tabular}
      \caption{Time spent in each module.}
      \label{tab:T_times}
	\end{center}
\end{table} 
%
\begin{figure}[ht]
    \centering
    \includegraphics{T.pdf}
    \caption{Time spent in each module.}
    \label{fig:T}
\end{figure}
%%
%%
%%----------------------------
%\subsection{Implementation Test Considering Secondary Flow and Bed Slope Effects}
%%----------------------------
%%
%%
%\label{subsec:im_sec}
%%

%
%----------------------------
%
\section{Application to Real Scenarios}
%
%----------------------------
%
\label{sec:apl}
%
In this section we apply the results of the linear analysis and the check routine to real scenarios at a laboratory scale (Section \ref{subsec:ashida90}) and at a field scale (Section \ref{subsec:dvr}). In the first case we reproduce a set of laboratory experiments were we face the problem of ill-posedness, we show the consequences and the possibilities that a modeler has to obtain a well-posed model. We also show the limitations of the linear analysis when applied to a real scenario. In the second case we re-run a simulation used in the DVR project obtaining the output of the ill-posed routine. 
%
%
%----------------------------
\subsection{Laboratory Scale}
%----------------------------
%
%
\label{subsec:ashida90}
%
In this section we reproduce the laboratory experiments conducted by \citet{Ashida90}.

They conducted several laboratory experiments in a curved flume. The width of the flume was \SI{0.2}{m}. The planform was formed out of 4.5 wavelengths of a sine-generated curve with a maximum angle equal to $\theta_m=35^\circ$ and a wavelength equal to $\lambda_s=\SI{2.2}{m}$. The bed was covered with sediment. Initially, the bed was flat and sloping in downstream direction $s=0.009$ (Figure \ref{fig:ashida_dom}). A constant water discharge was set upstream and the water level was kept constant downstream. The experiments run until equilibrium conditions were found at the third wavelength starting from upstream. To our understanding, the was no sediment discharge input at the upstream end. Thus, in strict terms, mobile bed equilibrium conditions cannot be reached. Yet, since the measurement location is far downstream from the boundary conditions (in terms of propagation celerity of bed disturbances) there is a period of steady state before disturbances from upstream affect the measurement region. We find that equilibrium conditions are well obtained after \SI{10}{min}. For all cases we run \SI{2}{min} of warm up without morphodynamic changes which is sufficient for the flow to achieve a steady state.  
%
\begin{figure}[ht]
    \centering
    \includegraphics{ashida_dom.png}
    \caption{Initial bed elevation of the model of the experiments by \citet{Ashida90}.}
    \label{fig:ashida_dom}
\end{figure}

We focus on their cases 4 and 2 which we label A1 and A2 respectively. The experimental parameters are shown in Table \ref{tab:A_eq}. Case A1 was run with unisize sediment and A2 with mixed-size sediment (Figure \ref{fig:A_gsd}). 
%
\begin{figure}[ht]
    \centering
    \includegraphics{A_gsd.eps}
    \caption{Grain size distribution of the sediment in Case A2, original (blue) and discretized (orange).}
    \label{fig:A_gsd}
\end{figure}
%
\begin{table}[ht]
%\tiny
	\begin{center}
	\begin{tabular}{cccccc}		
			\toprule
Case & Q [\si{m^2/s}] & $h$ [\si{m}] & $q_b$ [\si{m^2/s}] & $C_f$ [$-$] & sed. \\
			\midrule
A1 & 0.0036 & 0.0411 & \num{1.008e-05} & 0.0189 & uniform \\ 
A2 & 0.0036 & 0.0426 & \num{1.258e-05} & 0.0211 & mixed  \\
 		\bottomrule
			\end{tabular}
      \caption{Parameters of the laboratory experiments conducted by \citet{Ashida90}. The flow depth $h$ is an average value on one wavelength. The sediment transport rate is a cross sectional average.}
      \label{tab:A_eq}
	\end{center}
\end{table}

We try to resemble the numerical simulations that \citet{Ashida90} conducted of their experiments. They used the sediment transport relation by \citet{Ashida72}. The critical shear stress was computed using the relation by \citet{Iwagaki56}. We use a constant value equal to 0.05 (the original value in the relation by \citet{Ashida72}). They did not model secondary flow but the direction of the sediment transport rate was corrected due to secondary flow assuming an equilibrium state \citep{Engelund74}. The bed slope effect on the sediment transport rate was accounted for using the relation developed by \citet{Hasegawa83}. This relation has a dependency on the bed shear stress as in the one by \citet{Talmon95}. We use the relation by \citet{Talmon95} with parameters $A_{s}=1.7$, $B_{s}=0.5$ and $C_s=D_s=0$. 

They account for diffusion in the momentum equations. Note that this is not equivalent to the diffusion in the constitutive equation for the secondary flow intensity and that we have not included this effect in our study. The streamwise diffusion of the streamwise momentum equation and the crosswise diffusion of the crosswise momentum equation are equal to $2D_{Hn}$ (Equation \ref{eq:DHn}). The crosswise diffusion of the streamwise momentum equation and the streamwise diffusion of the crosswise momentum equation are equal to $D_{Hn}$. We use the average values of flow depth and shear velocity that yields $D_{Hn}=\SI{1.6907e-04}{m^2/s}$.

First we analyze the modeling of secondary flow and its effect on the final bathymetry (Section \ref{subsubsec:ashida90_sf}). At continuation we study the role of bed slope effects under mixed-sediment conditions (Section \ref{subsubsec:ashida90_bs}). In Section \ref{subsubsec:ashida90_deg} we study a hypothetical experiment conducted under degradational conditions into a substrate finer than the active layer. 
%
%
%----------------------------
\subsubsection{Modeling of Secondary Flow}
%----------------------------
%
%
\label{subsubsec:ashida90_sf}
%
In this section we focus on the modeling of secondary flow and its effect on the final bathymetry of the experiment conducted by \citet{Ashida90}. To this end we reproduce Case A1. Figure \ref{fig:A1_me} presents the measured bed elevation at equilibrium. We analyze the resulting topography for different choices on the secondary flow modeling (Table \ref{tab:A1_var}). 

Since secondary flow may yield an ill-posed model (Section \ref{subsec:secflow_uni}), a first possibility may be to not account for secondary flow (Case A1.1, Table \ref{tab:A1_var}). The final bathymetry does not seem realistic (Figure \ref{fig:A11_etab}). We do not predict the order of magnitude of the transverse slope. The inner part is significantly lower than measured and the opposite occurs in the outer part of the bed. 

The second approach is to model secondary flow using a diffusion coefficient computed using a closure relation. The linear analysis predicts for this case that for all closure relations (i.e., streamwise value according to \citeauthor{Elder59}, Equations (\ref{eq:DHs}), crosswise value according to \citeauthor{Elder59}, Equations (\ref{eq:DHs}), or the modification including flow curvature developed by \citet{Fischer69}, Equation (\ref{eq:DHn_sf})) the system of equations is ill-posed. We run a simulation using the crosswise diffusion coefficient developed by \citet{Elder59} (Case A1.2, Table \ref{tab:A1_var}). The bathymetry is more realistic than without secondary flow (Figure \ref{fig:A12_etab}). The order of magnitude of the transverse bed slope is correct. At the bed entrance, on the left-hand-side (up), the maximum bed elevation is slightly larger than \SI{3}{cm}, as measured in the experiment. However, on the right-hand-side, they measured a deeper scour than predicted. The position of the bar after the bed is well-captured although the maximum bed elevation is lower than measured. Despite the reasonable agreement, we find that the bed elevation is not stable. The bed elevation oscillates with time at the crests of the bars (Figure \ref{fig:A12_xz}). The wavelength of the oscillations is related to the grid size. We confirm that throughout the simulation the Courant condition has been fulfilled and the diffusion solver has always converged. 

Given that the solution oscillates for the relatively low value of the diffusion coefficient used in Case A1.2, we run a third case in which the diffusion coefficient is set to $D_{H}=\SI{0.5}{m^2/s}$ which is large enough to provide a well-posed simulation according to the linear analysis (Case A1.3, Table \ref{tab:A1_var}). The bed elevation does not present the oscillations. However, the final bathymetry is less realistic than in Case A1.2 (Figure \ref{fig:A13_etab}). 

Provided that \citet{Ashida90} did not measure the flow velocity for this experiment, there is no data to which compare our numerical results. Yet, we study the development of the secondary flow intensity from a theoretical perspective. \citet{deVriend81_PhD} studied the decay of secondary flow intensity in a straight reach after a bend. Assuming shallow hydrostatic flow, a logarithmic  vertical profile of flow velocity, and the streamwise velocity to be a function of the depth averaged streamwise velocity (i.e., similarity approximation) a simplified system of equations is obtained \citep{deVriend81_PhD}. This system is solved numerically to obtain that secondary flow decreases approximately exponentially with a relaxation length $L_I=1.3h/\sqrt{C_F}$. In our case the adaptation length is equal to $L_I=\SI{0.3886}{m}$ which implies that after \SI{0.89}{m} secondary flow has decayed to a 10\% of the initial value. The equilibrium value of secondary flow is $I_e=Q/R_s=\SI{0.0314}{m/s}$ where we have used that the curvature at the apex of a sine-generated curve is $R_s=\lambda_s/(2\pi\theta_m)=\SI{0.5732}{m}$. 

The secondary flow intensity of Case A1.2 is to a certain extent in accordance with the relaxation theory (Figure \ref{fig:A12_I}). The maximum absolute value is slightly larger than the equilibrium value. Secondary flow decays to zero in approximately \SI{0.7}{m}. A large value of the diffusion coefficient provides a well-posed model a no instability is found, however, the secondary flow pattern of Case A1.3 is unrealistic (Figure \ref{fig:A13_I}). There is no variation in the crosswise direction and the equilibrium value is not reached.
%
\begin{table}[ht]
%\tiny
	\begin{center}
	\begin{tabular}{ccccc}		
			\toprule
Case & Sec. flow & $D_H$ [\si{m^2/s}] & Math. char. & Result \\
			\midrule
A1.1 & N & - & W & unrealistic \\ 
A1.2 & Y & \num{1.6907e-04} & I & oscillations \\ 
A1.3 & Y & 0.50 & W & stable \\ 
 		\bottomrule
			\end{tabular}
      \caption{Simulations of Case A1 to study the consequences of the modeling choices as regards to secondary flow.}
      \label{tab:A1_var}
	\end{center}
\end{table}
%
\begin{figure}[ht]
    \centering
    \includegraphics[width=\textwidth]{A1_me.png}
    \caption{Measured bed elevation of case A1 (from \citet{Ashida90}).}
    \label{fig:A1_me}
\end{figure}
%
\begin{figure}[ht]
    \centering
    \includegraphics{A11_etab.png}
    \caption{Predicted bed elevation of Case A1 without secondary flow (Simulation A1.1).}
    \label{fig:A11_etab}
\end{figure}
%
\begin{figure}[ht]
    \centering
    \includegraphics{A12_etab.png}
    \caption{Predicted bed elevation of Case A1 with secondary flow and $D_H=\SI{1.6907e-04}{m^2/s}$ (Simulation A1.2).}
    \label{fig:A12_etab}
\end{figure}
%
\begin{figure}[ht]
    \centering
    \includegraphics{A13_etab.png}
    \caption{Predicted bed elevation of Case A1 with secondary flow and $D_H=\SI{0.50}{m^2/s}$ (Simulation A1.3).}
    \label{fig:A13_etab}
\end{figure}
%
\begin{figure}[ht]
    \centering
    \includegraphics{A12_xz.eps}
    \caption{Longitudinal profile of bed elevation with time at $y=\SI{0.0466}{m}$ to the left of the centerline ($N=12$). The results correspond to the Simulation of Case A1 with secondary flow and $D_H=\SI{1.6907e-04}{m^2/s}$ (Simulation A1.2).}
    \label{fig:A12_xz}
\end{figure}
%
\begin{figure}[ht]
    \centering
    \includegraphics{A12_I.png}
    \caption{Secondary flow intensity of Case A1 with $D_H=\SI{1.6907e-04}{m^2/s}$ (Simulation A1.2).}
    \label{fig:A12_I}
\end{figure}
%
\begin{figure}[ht]
    \centering
    \includegraphics{A13_I.png}
    \caption{Secondary flow intensity of Case A1 with $D_H=\SI{0.50}{m^2/s}$ (Simulation A1.3).}
    \label{fig:A13_I}
\end{figure}

%\errmessage{2do}
%\textcolor{red}{
%\begin{itemize}
%\item 3D simulation L026 see flow distribution
%\item 2D simulation with better distribution of secondary flow
%\end{itemize}
%}
%
%
%----------------------------
\subsubsection{Modeling of the Bed Slope Effects}
%----------------------------
%
%
\label{subsubsec:ashida90_bs}
%
In this section we center our attention on the modeling of the bed slope effects. To this end we try to model the mixed-size sediment experiment conducted by \citet{Ashida90} (Case A2). 

We start from the simplest mixed-size sediment case in which we consider two sediment size fractions (mixture M1 in Figure \ref{fig:A_gsd}). In Section \ref{subsubsec:ashida90_sf} we have seen that secondary flow without a large diffusion coefficient yields an ill-posed model. For this reason we first neglect secondary flow. The substrate sediment is initially the same as the active layer one. We analyze this case at a linear level and we obtain that it is well-posed when the bed slope effects are taken into consideration using the closure relation developed by \citet{Sekine92_1} but it is ill-posed when we use the relation by \citet{Talmon95}. This result is in agreement with the results of the linear analysis (Section \ref{subsec:mixsed}) we have shown that under aggradational conditions the diffusive term in the active layer equation due to bed slope effects is negative (which implies that the model is well-posed) provided that the bed slope effects are modeled using the relation by \citet{Sekine92_1}.

We run two numerical simulations which differ in the closure relation for the bed slope effects (Table \ref{tab:A2_var}). We would expect that the results of Case A2.2 show the oscillatory pattern characteristic from ill-posed simulations. However, this does not occur and the results of both simulations are smooth and the solution converges to a steady state (Figures \ref{fig:A21_xz} and \ref{fig:A22_xz}). 
%
\begin{table}[ht]
%\tiny
	\begin{center}
	\begin{tabular}{cccccc}		
			\toprule
Case & Sec. flow & Bed slp. & $N$ & Math. char. & Result \\
			\midrule
A2.1 & N & \citeauthor{Talmon95}   & 2  & I & stable \\ 
A2.2 & N & \citeauthor{Sekine92_1} & 2  & W & stable \\ 
A2.3 & Y & \citeauthor{Talmon95}   & 10 & I & stable \\ 
 		\bottomrule
			\end{tabular}
      \caption{Simulations of Case A2 to study the consequences of the modeling choices as regards to bed slope effects.}
      \label{tab:A2_var}
	\end{center}
\end{table}
%
\begin{figure}[ht]
    \centering
    \includegraphics{A21_xz.eps}
    \caption{Longitudinal profile of bed elevation with time at $y=\SI{0.0466}{m}$ to the left of the centerline ($N=12$). The results correspond to the Simulation of Case A2 using the bed slope effects relation by \citet{Talmon95} (Simulation A2.1).}
    \label{fig:A21_xz}
\end{figure}
%
\begin{figure}[ht]
    \centering
    \includegraphics{A22_xz.eps}
    \caption{Longitudinal profile of bed elevation with time at $y=\SI{0.0466}{m}$ to the left of the centerline ($N=12$). The results correspond to the Simulation of Case A2 using the bed slope effects relation by \citet{Sekine92_1} (Simulation A2.2).}
    \label{fig:A22_xz}
\end{figure}

The fact that the linear analysis predicts that the model is ill-posed but the solution does not present any symptoms stresses the limitations of the linear analysis. The linear analysis assumes a reference situation without flow curvature (Section \ref{subsec:lin}). Yet, the flow in the initial situation of this case is curved which implies, for instance, that the transverse bed slope is not equal to 0 as we used in the linear analysis. It may be possible that the model is initially ill-posed. However, superimposed to the hypothetical spurious oscillations there is the natural and physical growth of bars. Eventually, after a very short time the bar growth may turn the model to be well-posed. Note that, if this is the case, a decrease of the grid size would significantly change the solution but the results show that using this grid, even if those spurious oscillations are present, they are negligible when compared to the growth of bars. 

Eventually, we try to reproduce the experimental results of Case A2 using the knowledge gained with the previous simulations. We have seen that we need to account for secondary flow in order to obtain a reasonable bathymetry. Yet, the we need to use a large value of the diffusion coefficient to obtain a well-posed model. In principle the bed slope function by \citet{Talmon95} yields an ill-posed model but we have observed that this does not seem to be the case. For this case (A2.3) we discretize the sediment mixture in 10 grain sizes (mixture M2 in Figure \ref{fig:A_gsd}). We compare the grain size distribution of the bed surface at the end of the run (Figure \ref{fig:A23_dm}) with the measured values (Figure \ref{fig:A2_me}).
%
\begin{figure}[ht]
    \centering
    \includegraphics[width=\textwidth]{A2_me.png}
    \caption{Measured mean grain size at the bed surface of case A2 (from \citet{Ashida90}).}
    \label{fig:A2_me}
\end{figure}
%
\begin{figure}[ht]
    \centering
    \includegraphics{A23_dm.png}
    \caption{Predicted mean grain size at the bed surface of case A2 (Simulation A2.3).}
    \label{fig:A23_dm}
\end{figure}

The result is unsatisfactory. The predicted data show a more homogeneous grain size distribution along the bed surface than measured. This is somehow expected since the large diffusion coefficient, necessary to obtain a well-posed model, flattens the bed preventing the creation of coarse and fine regions. 
%
%
%----------------------------
\subsubsection{Modeling of a Degradational Case}
%----------------------------
%
%
\label{subsubsec:ashida90_deg}
%
In the previous sections (\ref{subsubsec:ashida90_sf} and \ref{subsubsec:ashida90_bs}) we dealt with ill-posedness due to secondary flow and bed slope effects. In this section we deal with ill-posedness due to degradation into a substrate finer than the active layer. 

The substrate of the mixed-size sediment experiment conducted by \citet{Ashida90} had the same grain size distribution than the active layer. In these conditions the active layer model is well-posed given that reverse mobility does not occur and that the aggradational flux to the substrate has the same grain size distribution than the active layer \citep{Chavarrias18}. Thus, we conduct a though experiment to study ill-posedness due to a fine substrate (Case A3). The experiment is the same as Case A2 (Table \ref{tab:A_eq}), however, the substrate is composed of only the fine fraction. The downstream water level is kept constant for \SI{10}{min} until equilibrium conditions are reached. At continuation, a lowering at a rate of \SI{1.5}{cm/h} is imposed. The sediment is discretized in two sediment size fractions (mixture M1 in Figure \ref{fig:A_gsd}). Secondary flow is taken into consideration using a diffusion coefficient equal to \SI{0.50}{m^2/s}. 

During the first \SI{10}{min} bars form what causes a relatively small amount of aggradation in some areas and degradation in others. Thus, in the equilibrium configuration the grain size distribution of the top part of the substrate is not equal to the initial one. After \SI{2}{h} of lowering of the downstream water level the substrate sediment is entirely composed of fine sediment since the coarse substrate sediment initially created during the bar formation period has been entrained. In this conditions the model is ill-posed according to the linear analysis and the bed elevation with time presents oscillations (Figure \ref{fig:A31_xz}). 
%
\begin{figure}[ht]
    \centering
    \includegraphics{A31_xz.eps}
    \caption{Longitudinal profile of bed elevation with time at $y=\SI{0.0466}{m}$ to the left of the centerline ($N=12$). The results correspond to the Simulation of Case A3.}
    \label{fig:A31_xz}
\end{figure}
	
%%   case 4: A1 (unisize)
%%   h=0.0411 m
%%   slope=0.009
%%   q_m=3.6/1000=0.0036 m^3/s
%%   q=q_m/B=0.0036/0.2=0.018 m^2/s
%%   u=q/h=0.0180/0.0411=0.4380 m/s
%%   Fr=u/sqrt(gh)=0.4380/sqrt(9.81*0.0411)=0.69
%%   Cf(1)=u_st^2/u^2=(6.02/100)^2/0.4380^2=0.0189
%%   Cf(2)=h^3*g*s/q^2=0.0411^3*9.81*0.009/0.018^2=0.0189
%%   qb_m=5.34 gr/s=320.4 gr/min=0.00534 kg/s=1.007547169811321e-05 m^2/s (without pores)=1.679245283018868e-05 m^2/s (with pores)=3.358490566037736e-06 m^3/s (with pores);
%%   \mu_2=0.6, \mu_k=0.48 -> A_sk=0.53, B=?? 1/2 -> A_sk=1.7 B_sk=0.5
%%   dk=0.00174
%
%%   L004 sf=N 							-> unrealistic pattern
%%   L006 sf=Y Dh crossmwise -> unstable
%%   L005 sf=Y Dh streamwise -> stable and realistic but strange bar after apex.
%%   
%%   Dh streamwise = \left(\frac{0.4041}{\kappa^3}+\frac{1}{6}\kappa\right)hu^* = (0.4041/0.41^3+1/6*0.41)*0.0411*6.02/100=0.0147 m^2/s 
%%   Dh crosswise = 1/6*\kappa*h^u_st=1/6*0.41*0.0411*6.02/100=1.6907e-04 m^2/s -> dampening of largest wave numbers = L006 unstable 
%%   Dh crosswise with curvature = 0.25*0.4380^2*0.0411^3/0.41^5/0.0602/0.5732^2 = 0.0145 m^2/s
%%   radious of curvature at apex = L/2/pi/maxangle= 2.2/2/pi/(35/180*pi) = 0.5732 m





%%   case 2: A2 (mixed-size)
%%   h=0.0426 m
%%   slope=0.009
%%   q_m=3.6/1000=0.0036 m^3/s
%%   q=q_m/B=0.0036/0.2=0.018 m^2/s
%%   u=q/h=0.0180/0.0426=0.4225 m/s
%%   Fr=u/sqrt(gh)=0.4225/sqrt(9.81*0.0426)=0.65
%%   Cf(1)=u_st^2/u^2=(6.13/100)^2/0.4225^2=0.0211
%%   Cf(2)=h^3*g*s/q^2=0.0426^3*9.81*0.009/0.018^2=0.0211
%%   qb_m=6.67 gr/s=400.2 gr/min=0.00667 kg/s=1.258490566037736e-05 m^2/s (without pores)=2.097484276729560e-05 m^2/s (with pores)=4.194968553459119e-06 m^3/s (with pores);
%%   \mu_2=0.6, \mu_k=0.48 -> A_sk=0.53, B=?? 1/2 -> A_sk=1.7 B_sk=0.5
%%   dk=0.001 0.002 -> Fak=0.56
%%   Dh streamwise = \left(\frac{0.4041}{\kappa^3}+\frac{1}{6}\kappa\right)hu^* = (0.4041/0.41^3+1/6*0.41)*0.0426*6.13/100=0.0155 m^2/s 
%%   Dh crosswise = 1/6*\kappa*h^u_st=1/6*0.41*0.0426*6.13/100=1.7844e-04 m^2/s
%%   Dh crosswise with curvature = 0.25*0.4225^2*0.0426^3/0.41^5/0.0613/0.5732^2 = 0.0148 m^2/s
%%   radious of curvature at apex = L/2/pi/maxangle= 2.2/2/pi/(35/180*pi) = 0.5732 m

%
%
%----------------------------
\subsection{DVR Simulation}
%----------------------------
%
%
\label{subsec:dvr}
%
In this section we apply the tool to a real scenario at a field scale. 

We choose a subdomain of the DVR model developed to study the effects of the sediment nourishments in the Rhine River at its entrance to The Netherlands \citep{Ottevanger15}. In particular, we select the upper part of the domain, the Bovenrijn, between the bridge at Emmerich (Germany) and the Pannerdenschekop (Figure \ref{fig:br2_dom}). In the DVR model the natural hydrograph is discretized into 8 stages that are combined to create a syntethic hydrograph with intrayear and interyear variations. We run the subdomain for one week of real time (the morphodynamic acceleration factor \citep{Latteux95, Roelvink06,Ranasinghe11} is set to 1) and we select the initial and boundary conditions for a water discharge at Lobith equal to \SI{2251}{m^3/s}. The DVR model runs in domain decomposition. Since we only run one subdomain, we use the initial condition to fix the water level at the downstream end. The sediment is discretized into 10 size fractions ranging between \SI{0.063}{mm} and \SI{64}{mm}. 
%
\begin{figure}[ht]
    \centering
    \includegraphics[width=\textwidth]{br2_dom2.jpg}
    \caption{Rhine at its entrance to The Netherlands. The red mesh is the computational grid (Bovenrijn). (Adapted from Google Earth)}
    \label{fig:br2_dom}
\end{figure}

The morphological simulation starts after three days of spin-up. Figure \ref{fig:D_h0} shows the initial flow depth where we can appreciate the deep outer bend and shallow inner side. Figures \ref{fig:D_Fak0} and \ref{fig:D_fIk0} present the mean grain size of the sediment at the bed surface and at the top part of the substrate, respectively. Note that the substrate is not defined for the whole domain. A pattern of fine stripes is clearly seen in the initial condition. The interpretation of this data is outside the scope of this project. 
%
\begin{figure}[ht]
    \centering
    \includegraphics{D_h0.png}
    \caption{Initial flow depth.}
    \label{fig:D_h0}
\end{figure}
%
\begin{figure}[ht]
    \centering
    \includegraphics{D_Fak0.png}
    \caption{Initial mean grain size at the bed surface.}
    \label{fig:D_Fak0}
\end{figure}
%
\begin{figure}[ht]
    \centering
    \includegraphics{D_fIk0.png}
    \caption{Initial mean grain size at the top substrate.}
    \label{fig:D_fIk0}
\end{figure}

In Figure \ref{fig:D_hir} we show the locations in which the routine implemented in Delft3D predicts the model to be ill-posed after 10 days. Note that at scattered locations, isolated nodes appear to be ill-posed. During the simulation some of the locations in which the model is ill-posed change. Figure \ref{fig:D_hir_c} present the locations that at some point have been ill-posed. We note that the majority of the locations that appear to have been ill-posed at some point were ill-posed during the initial part of the simulation. Figure \ref{fig:D_xzt2} shows the bed elevation with time for a longitudinal section in which none of the nodes are ill-posed throughout the simulation time. Small and gradual changes in bed elevation occur which is reasonable provided that the simulation time is only 7 days. Figure \ref{fig:D_xzt1} present the bed elevation for a longitudinal cross section along nodes which are almost continuously ill-posed during the simulation time. In this case we observe a substantial change in bed elevation and the formation of waves with wavelength related to grid size, as we would expect from an ill-posed simulation.
%
\begin{figure}[ht]
    \centering
    \includegraphics{D_hir.png}
    \caption{Ill-posed locations at the end of the simulation.}
    \label{fig:D_hir}
\end{figure}
%
\begin{figure}[ht]
    \centering
    \includegraphics{D_hir_c.png}
    \caption{Locations that at some point have been ill-posed.}
    \label{fig:D_hir_c}
\end{figure}
%
\begin{figure}[ht]
    \centering
    \includegraphics{D_xzt2.eps}
    \caption{Bed elevation with time along a cross section in which the nodes are mainly well-posed (M=50; N=100:130).}
    \label{fig:D_xzt2}
\end{figure}
%
\begin{figure}[ht]
    \centering
    \includegraphics{D_xzt1.eps}
    \caption{Bed elevation with time along a cross section in which the nodes are mainly ill-posed (M=60; N=80:110).}
    \label{fig:D_xzt1}
\end{figure}

The fact that the system of equations is in some nodes and at some point in time ill-posed implies that the solution is subjective and needs to be carefully assessed. First, the fact that the solution has been always well-posed in a certain subdomain does not imply that it does not suffer from the consequences of ill-posedness. This is because perturbations that are generated in the ill-posed part of the domain may travel to parts of the domain that are well-posed. Second, the fact that the solution is at some point ill-posed does not mean that the consequences are appreciable. For instance, if a limited number of nodes become ill-posed for a short time, oscillations may not have time to significantly grow. Third, the solution is grid dependent. A refinement of the grid will not cause convergence of the solution but the opposite (Section \ref{subsec:def_ill}). The difficulty arises in the fact that it is up to the modeler to decide whether the consequences are significant or not and whether the solution is still reasonable and usable to answer the questions for which the model was made. 
%
%----------------------------
%
\section{Discussion}
%
%----------------------------
%
\label{sec:disc}
%
We have studied the mathematical character (ill-posedness or well-posedness) of the set of equations used to model two-dimensional river morphodynamics. In this section we discuss the results, limitations, and we address further research on the topic.

To study the mathematical character we have linearized the system of equations around a reference state which is a solution to the system of equations. We have used the simplest solution which is that of uniform steady flow on a flat sloping bed with no transverse slope. \emph{Sensu stricto} our results are only valid in conditions close to this reference state. A valid question that arises is what would be the mathematical character of the same model equations but using as a reference state, for instance, a solution with constant secondary flow intensity (i.e., an infinite bend). We do not expect the core of the results to be different. This is because we have shown that there are essential problems in the model equations and it does not seem plausible that those problems do not exist if the reference state is different. Worded differently, ill-posedness highlights a lack of physical mechanisms in the model equations and this is independent of the reference solution we consider in the analysis. For instance, the main problem, that of ill-posedness caused by the active layer model under degradational conditions, is present even in 1D, so it does not seem logical that it will not be present, for instance, in the case of flow in a bend. Moreover, we have shown that the domain of ill-posedness obtained considering as a reference solution the simplest scenario yields a good estimation of the domain of ill-posedness in conditions with strong curvature (i.e., the experiment conducted by \citet{Ashida90}, Section \ref{subsec:ashida90}) and a complex flow pattern (i.e., the DVR simulation, Section \ref{subsec:dvr}). Although we do not expect the domain of ill-posedness to be significantly affected by the reference solution used to linearize the system of equations, insight can be gained from such an analysis. For instance, the propagation direction of disturbances, the absolute or convective character of instabilities, and the celerity of perturbations. This knowledge is essential to properly set the boundary conditions of a numerical model, predict the formation of bars, or to obtain analytical solutions to test numerical implementations \citep[e.g.][]{Lanzoni06}.

In the analysis we have neglected the diffusive character of the eddy viscosity in the momentum equations. This simplifies the analysis since the only diffusive component in the model is brought by the transport of the secondary flow intensity and bed slope effects. Since a diffusive component tends to turn the model to be well-posed, it would be of interest to study the effect of this variable. Note that the eddy viscosity does not only introduce diffusive (second derivative) terms but also first order (first derivative) terms. This mechanism most probably has an impact on the results considering ill-posedness due to the secondary flow model because both mechanisms affect the flow. However, we do expect that eddy viscosity has an impact on ill-posedness due to the active layer model since the model equations are only mildly coupled. Numerical simulations run with eddy viscosity point in this direction. 

An important limitation of our results as regards to the domains of ill-posedness is that it is purely based on empirical insight. Worded differently, our conclusions are based on observations of the ill-posed domain of selected cases. The complexity of the system hinders obtaining analytical expressions to prove the existence or absence of ill-posedness. Our results, however, prove that the system may be ill-posed. We recommend a thorough analysis in which the equations are non-dimensionalized to study the characteristic polynomial of the system in order to prove that a system is ill-posed or well-posed without the necessity of solving an eigenvalue problem. The analytical study would require a simplification of the system (e.g., steady flow, or simple closure relations). Thus for practical cases in numerical simulation we prefer to obtain the eigenvalues numerically including all possible complexity. Moreover, in numerical simulations there is a minimum possible wavenumber related to the domain size and a maximum related to the grid size. This means that when numerically obtaining the domain of ill-posedness one can use the grid size as a cut-off for the maximum wavenumber and be sure that no spurious oscillations can appear above that wavenumber if the results show that the model is well-posed. We have observed that diffusion plays a significant role and that there appears to be a minimum amount of diffusion that yields a well-posed model as regards to secondary flow. Once the role of eddy viscosity is introduced in the system of equations we recommend the study of the characteristic polynomial to obtain an analytical expression for the minimum amount of diffusion necessary to obtain a well-posed model. It is of interest to compare the minimum value with physically realistic estimations of the horizontal diffusion coefficient. 

Our study has focused on the eigenvalues of the system. This yields the growth rate of oscillations. If the growth rate is negative oscillations decay exponentially with time. Thus, in an infinite time oscillations disappear. This analysis, however, provides no insight in the temporal evolution of the solution. It is possible that, although in the long term oscillations decay, in the short term they grow due to non-orthogonality of the eigenvectors \citep{Darmofal96}. The growth may be sufficient to limit the assumption of linearity. In this case, the eigenvalue analysis predicts decay of oscillations and well-posedness but the model shows growth which may be catastrophic. We recommend further research in this regard.

We have shown that the secondary flow model has severe limitations. In particular, for the experiment conducted by \citet{Ashida90}, the necessary diffusion to obtain a well-posed model is unrealistically large, which causes an unrealistic flow field. The ill-posedness indicates that key physical parameters are not taken into consideration and severe simplifications have been conducted. For instance, the model does not consider the anisotropy of the diffusion coefficient. We recommend to implement this possibility in Delft3D, which is useful for any transported constituent. Yet, we do not expect that this will yield a well-posed model since even a large diffusion coefficient (when compared to the one derived by \citet{Elder59}) in all directions may yield an ill-posed model. We recommend to further develop the model for secondary flow for instance, parametrizing its intensity in the crosswise direction. 

We have implemented a routine in Delft3D that allows us to check whether a simulation turns to be ill-posed. The routine has been thought to efficiently identify the problem of ill-posedness due to the mixed-size sediment model. The diffusive terms are not included at this moment so it is no possible to check whether the horizontal diffusion coefficient in the transport equation of the secondary flow intensity is sufficient to yield a well-posed model. We suggest to first study the effect of eddy viscosity in the momentum equations and then include all terms in the ill-posedness routine. The implementation will depend on the results of the study of the model equations including diffusion. Yet, if the consequences of eddy viscosity are similar to those of the diffusion coefficient in the secondary flow intensity equation, the check routine will have the added complexity of a solution depending on the wavenumber. We estimate that the procedure should be: (1) to determine the maximum wavenumber in $x$ and $y$ direction based on the grid size, and a minimum wavenumber based on the domain length, (2) to discretize the range of possible wavenumbers, and (3) to loop and check that for the maximum possible value in one direction (e.g., the $x$ direction), all wavenumbers in $y$ direction have a negative growth rate. We expect that the computational cost will significantly increase. For this reason we recommend to set this test independently from the test for the mixed-size sediment model, which we consider to be the most important cause of ill-posedness. If found, simple analytical expressions could be implemented in the routine to reduce the computational time. 

This project has focused exclusively on the identification of the problem of ill-posedness and its consequences. Eventually we aim at solving the problem. A solution can only come from an improved set of equation that take into consideration key physical processes currently neglected \citep{Joseph90}. We are currently working on a possible solution for the problem of ill-posedness caused by the mixed-size sediment model. The solution consist on a modification of the time scale of the mixing processes to recover the well-posed character of the system of equations. This technique is applied at those spatio-temporal locations in which the model is ill-posed. Thus, the current routine in Delft3D to identify the locations in which the model is ill-posed can be seen as a first step towards the solution of the problem of ill-posedness. Eventually, we desire to develop a model for mixed-size sediment morphodynamics that includes the necessary physical mechanisms to obtain an unconditionally well-posed model.
%
%----------------------------
%
\section{Conclusions}
%
%----------------------------
%
\label{sec:conc}
%
%
We have studied the mathematical character (well-posedness or ill-posedness) of the system of equations used to model river morphodynamics in two dimensions including the effects of flow curvature (i.e., secondary flow and bed slope effects). To this end we have conducted a perturbation analysis of the linearized system of equations. Studying the growth rate of perturbations as a function of the wavenumber we conclude the following:
\begin{itemize}
\item The current secondary flow model may yield an ill-posed model. This happens even when morphodynamic changes not are taken into consideration. 
\item Above a certain threshold of horizontal diffusion in the transport equation for the secondary flow intensity the model is well-posed. The threshold sometimes reaches values that are physically unrealistic. 
\item The effect of bed slope on the sediment transport direction is a necessary mechanism to obtain a well-posed model. 
\item Without considering the effects of secondary flow, a unisize model is well-posed. This is independent on the closure relation to account for the effects of bed slope in the sediment transport direction.
\item Without considering the effects of secondary flow, a two-dimensional mixed-size sediment model may be ill-posed in the conditions in which a one-dimensional model is well-posed. This depends on the closure relation to account for the effects of bed slope in the sediment transport direction. For the simplest case \citep[i.e.][]{Sekine92_1} a two-dimensional model is ill-posed only in the conditions in which a one-dimensional model is ill-posed too.
\end{itemize}

These results are summarized in a flow chart that may be used as a general guideline for modelers to be aware of ill-posedness (Figure \ref{fig:ellipticity}). 
%
\begin{figure}[ht]
    \centering
    \includegraphics[width=\textwidth]{ellipticity.pdf}
    \caption{Ill-posedness flow chart.}
    \label{fig:ellipticity}
\end{figure}

%We have run idealized 2D numerical simulations under conditions in which the system of equations is ill-posed. We conclude that:
%\begin{itemize}
%\item The consequences of ill-posedness in 2D are similar to those in 1D, infinitesimal perturbations grow...
%\end{itemize}

We have modeled the experiments conducted by \citet{Ashida90} in a curved flume to study the effects of the different model components (i.e., secondary flow, mixed-sediment, bed slope effects) on the well-posedness of a numerical simulation. The numerical results are compared to the laboratory data. We conclude that:
\begin{itemize}
\item The value of the diffusion coefficient in the transport equation of the secondary flow intensity required to achieve a well-posed model is in this case physically unrealistic. 
\end{itemize}

We have implemented a routine in Delft3D to check if a numerical simulation is ill-posed. The routine focuses on the problem of ill-posedness due to the mixed-size sediment module. We have tested the implementation of the matrices, the computation of the eigenvalues, and whether the output is correct.

We have used the routine in Delft3D to check a simulation of the DVR model used to assess dredging and dumping operation in the Rhine river. We conclude that:
\begin{itemize}
\item The DVR simulation suffers from ill-posedness. This leads to the growth of spurious oscillations in the solution.
\item The solution is subjective and needs to be carefully assessed. 
\item It is up to the modeler to decide whether the consequences are significant and whether the solution is reasonable and usable to address the problem for which the model was set up.
\end{itemize}
%
%----------------------------
%
\section{Recommendations}
%
%----------------------------
%
\label{sec:rec}
%
%
In this section we collect and summarize our recommendations. The justification of the recommendations is in the Discussion (Section \ref{sec:disc}). \\

We recommend to further study the mathematical system of equations as regards to ill-posedness focusing on the effects of:
\begin{enumerate}
\item the reference state used in the linearization,
\item the diffusive processes in the momentum equations. \\
\end{enumerate}

To ease the assessment of ill-posedness we recommend to study the possibility of obtaining:
\begin{enumerate}[resume]
\item an analytical discriminant that discerns between ill-posedness and well-posedness of the system of equations in 2D,
\item an expression to find the minimum amount of diffusion to obtain a well-posed model. \\
\end{enumerate}

To improve the modeling of secondary flow we recommend to implement in Delft3D:
\begin{enumerate}[resume]
\item an anisotropic diffusion coefficient,
\item a cross-wise parametrization of the secondary flow intensity. \\
\end{enumerate}

We recommend to extend the ill-posedness routine in Delft3D:
\begin{enumerate}[resume]
\item accounting for all diffusive processes,
\item including an analytical discriminant or bounds to reduce the computational time.\\
\end{enumerate}

As regards to mixed-sediment morphodynamics we recommend to:
\begin{enumerate}[resume]
\item once it is further studied and validated, implement in Delft3D the regularization strategy we are currently researching.\\
\end{enumerate}

The second recommendation is important as it is a necessary step before Recommendations 3, 4, 7, and 8 can be considered. Recommendation 1 may provide valuable insight but it is secondary in terms of importance. Our recommendations as regards to the improvement of the modeling of secondary flow are justified from a physical perspective. Moreover, we have shown the limitations of the current modeling approach. However, our assessment is based on a limited study and the task of improving the modeling of secondary flow would require deep insight from other cases and studies. We consider important to extend the routine in Delft3D (Recommendations 7 and 8). Yet, we think that at this point the most urgent task is to solve the problem of ill-posedness due to mixed-size sediment (Recommendation 9).

%  \end{linenomath*}
%\begin{table}[h!]
	%\begin{center}
	%\begin{tabular}{ >{\raggedright}p{2cm} | >{\centering\arraybackslash}p{5cm} | >{\centering\arraybackslash}p{5cm}  }		
			%\toprule
       %& Plane bed & Dune-dominated 
			%\\ 
			%\cmidrule(r){1-1}\cmidrule(lr){2-2}\cmidrule(l){3-3}
			%Fully mobile
			%&
			%\raisebox{-\totalheight}{\includegraphics[scale=0.4]{p1.jpg}}
      %& 
			%\raisebox{-\totalheight}{\includegraphics[scale=0.4]{d1.jpg}}
      %\\ 
			%\cmidrule(r){1-1}\cmidrule(lr){2-2}\cmidrule(l){3-3}
			%Selective transport
			%&
			%\raisebox{-\totalheight}{\includegraphics[scale=0.4]{p2.jpg}}
      %& 
			%\raisebox{-\totalheight}{\includegraphics[scale=0.4]{d2.jpg}}
      %\\ 
			%\cmidrule(r){1-1}\cmidrule(lr){2-2}\cmidrule(l){3-3}
			%Partial transport
			%&
			%\raisebox{-\totalheight}{\includegraphics[scale=0.4]{p3.jpg}}
      %& 
			%\raisebox{-\totalheight}{\includegraphics[scale=0.4]{d3.jpg}}
      %\\ 			
			%\bottomrule
      %\end{tabular}
      %\caption{Armor regimes (A. Blom)}
      %\label{tab:armor_regimes}
	%\end{center}
%\end{table}
 
%\begin{figure}[h!]
    %\centering
    %\includegraphics[width=0.8\textwidth]{sketch_model_4.pdf}
    %\caption{Model sketch. Arrows indicate fluxes between layers.}
    %\label{fig:sketch_model}
%\end{figure}

%\begin{table}[h!]
	%\begin{center}
	%\begin{tabular}{ l | l }		
			%\toprule
      %Task & GS credits \\
			%\midrule
			%Speak at NCR 2014 conference 			& 0.5 \\
			%Poster presentation at AGU 2014 	& 2 \\
			%Organization of Sorting Workshop 2015 & 4 \\
			%Writing RCEM 2015 conference paper & 1 \\
			%Correction of River Engineering 2014-15 exams & 3 \\
			%\bottomrule
			%\end{tabular}
      %\caption{GS credits on Research Skills already acquired.}
      %\label{tab:research_skills_done}
	%\end{center}
%\end{table}

%	
%
%----------------------------
%---------SCHEDULE
%----------------------------
%
%\section{Schedule}
%\includepdf[angle=90]{schedule.pdf}
%
%
%----------------------------
%---------REFERENCES
%----------------------------
%
\DeclareRobustCommand{\van}[3]{#3}
\bibliographystyle{agufull08_mod}
\bibliography{references}
%
%
%----------------------------
%---------END OF DOCUMENT
%----------------------------
%
%Include pdf
%\includepdf[angle=90]{mapa_suelos_almacelles_2.pdf}
%
\end{document}