\documentclass{deltares_report}
%\usepackage{booktabs}
%\usepackage{epigraph} %Epigrafes 
%\usepackage{eurosym} %euro symbol
%\usepackage[hyphens]{url}
%W%\usepackage{pdfpages} %para incluir pdf como anejos
\usepackage{natbib} %references
%\usepackage{fancyhdr}
%\setlength{\headheight}{14.5pt}
%\usepackage{authblk} %affiliations
%\usepackage{pdflscape} %\begin{landscape}
%\usepackage{xspace}
%
%---------version
%\usepackage{mVersion}
%\setVersion{3.1} %x.y y=1: V
%\increaseBuild
%
%
%---------figures
%\usepackage{graphics}
%\usepackage{graphicx}
%\usepackage{wrapfig} % Allows in-line images ?
%\usepackage{color}
%\usepackage{subfigure}
%\usepackage[section]{placeins}
%\usepackage[section,subsection,subsubsection]{extraplaceins}

%---------math
%\usepackage{cancel} %cancel term in equation
%\usepackage{amsmath}
%\usepackage{amssymb}
%\usepackage{bbm} %bold math 0 
%\usepackage{bbold} %bold math characters
%\usepackage{array} %math arrays
%\usepackage{chngpage} % allows for temporary adjustment of side margins
%\makeatletter
%\renewcommand*\env@matrix[1][*\c@MaxMatrixCols c]{%
%  \hskip -\arraycolsep
%  \let\@ifnextchar\new@ifnextchar
%  \array{#1}}
%\makeatother
%\usepackage{siunitx}

\hyphenation{ver-ge-lij-king-en}
\hyphenation{af-zon-der-lij-ke}

\begin{document}
%\title{Mathematical analysis of the well-posedness of Hirano \mbox{active} layer concept}
\title{Mathematical analysis of the well-posedness of the Hirano \mbox{active} layer concept in 2D \mbox{models}}
\subtitle{}
\date{\today}
\client{Rijkswaterstaat (WVL)}
\keywords{Hirano, graded sediment, morphodynamics, eigenvalues, well-posedness, Delft3D}
\references{}
\summary{Het voorliggende rapport beschrijft een wiskundige analyse van het actieve laag concept van Hirano, dat wordt gebruikt voor het berekenen van morfologische ontwikkeling met meerdere sedimentfracties in combinatie met de ondiep-water-vergelijkingen. 
	
Het criterium voor goedgesteldheid, zonder het effect van secundaire stroming of hellingseffekten, is geformaliseerd op basis van de bewegingsvergelijkingen. Het wordt aanbevolen om het effect van deze ontbrekende termen nader te bestuderen. 

De bewegingsvergelijkingen zijn omgeschreven tot een matrix vector vermenigvuldiging, waarbij de matrix in $x$- en $y$-richting worden beschreven door de Jacobiaan van het systeem in de respectievelijke richtingen. Deze matrices zijn bepaald met een symbolische toolbox in MATLAB. Hiermee zijn verschillende Monge kegels af te leiden om het effect van de modelparameters te bestuderen. Bovendien kan de analytische oplossing met behulp van MATLAB implementatie worden gebruikt om de toets voor goedsteldheid in Delft3D te valideren.

Een routine is toegevoegd aan de code van Delft3D, waarmee de Jacobiaan van de afzonderlijke asrichtingen van het systeem kan worden bepaald. Voor de hydrodynamica worden de analytische uitdrukkingen gebruikt, terwijl voor de morphologische vergelijkingen de afgeleiden worden bepaald via een numerieke differentiatie. Deze Jacobiaanen kunnen worden gebruikt om het systeem van de vergelijkingen te karakteriseren (d.w.z. goed- of slechtgesteld). De met Delft3D berekende Jacobiaanen geven hetzelfde resultaat als die berekend met behulp van MATLAB. Enkele keuzes met betrekking tot de de locatie van de gebruikte variabelen, heeft nog wat aandacht nodig. 

Daarnaast is een simulatie opgezet van het Pannerdensche Kop splitsingspunt. Deze kan worden gebruikt om de goed- of slechtgesteldheid van echte morfologische berekeningen te bepalen. Als de eigenwaarde berekening is toegevoegd aan Delft3D, kan deze simulatie worden gebruikt om te toetsen of de berekening goed- of slechtgesteld is, en wat er gebeurt in de Delft3D berekening als het probleem slechtgesteld is. 
}
\version{1.0}
\projectnumber{1230044}
\classification{}

\versioni{1.0 $Revision$}
\datei{9-12-2016} 
\authori{V\'ictor Chavarr\'ias (TU~Delft)} 
\authorii{Willem Ottevanger} 
\revieweri{Kees Sloff} 
\reviewerii{Astrid~Blom~(TU~Delft)} 
\approvali{Johan Boon}

\status{definitief}
\disclaimer{}
\reference{1230044-000-ZWS-0035}

\deltarestitle
%------------------------------------------------------------------------------

%----------------------------
%---------INTRODUCTION
%----------------------------
%
\chapter{Introduction}
\label{sec:introduction}
%
The current model to predict mixed-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 in 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 layer 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. Nevertheless in some circumstances the model may be ill-posed \citep{Ribberink87_PhD,Stecca14} and in these conditions can lead to unphysical and unrealistic results.

There are alternatives to the active layer model as the one by \citet{Ribberink87_PhD} where a second active layer is taken into account. This model was derived for dune-dominated cases and it does not solve the problem of ill-posedness \citep{Sieben97_PhD}. The continuous model by \citet{Blom03,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{Viparelli16} which is applicable at large scale but it can also become ill-posed \citep{ChavarriasXX}. Thus the active layer model is still the main model for mixed sediment river morphodynamics.

All previous analyses of the active layer model that we are aware of have been done assuming one-dimensional flow. Yet, to model the effects of flow curvature it is necessary to (at least) consider the bi-dimensionality of flow, which introduces a degree of complexity. While the equations to model the evolution of the bed are intrinsically the same only adding one dimension to the divergence terms, the flow equations substantially change between one and two dimensions. In two dimensions the water flow is modelled using the Shallow Water Equations (which in one dimension reduce into the \citet{SaintVenant71} equations). The Shallow Water Equations do not resolve for the vertical distribution of the velocity and are thus inadequate when, for instance, large curvature in a river bend induces a three dimensional flow (secondary flow). Yet, the secondary flow can be parametrized to avoid the expense of a three dimensional model \citep{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 this variable which is found in extra terms in the momentum equations. Thus, the consideration of secondary flow not only modifies the momentum equations but also adds an extra equation to the system. 

Another degree of complexity is added when taking into consideration the effect of the bed slope (in streamwise and lateral direction) in the sediment transport. This effect corrects for the fact that sediment is more easily transported downslope than upslope. A fundamental difference between this second degree of complexity when compared to the parametrization of secondary flow is that this does not introduce a dependent variable as it is modelled as a closure relation. Nevertheless, this closure relation also introduces a diffusive character just as the secondary flow.

Here we aim at building a check function to assess the well-posedness of 2D simulations run in Delft3D. To this end in Section \ref{chap:mmd} we describe the equations that model shallow water flow coupled to the active layer model. In Section \ref{chap:mc} we describe how we characterize a 2D system of partial differential equations and we obtain the matrix-vector formulation of the model. In Section \ref{chap:pr} we present the (preliminary) results of the implementation of the check function.
%
%----------------------------
%

%
%----------------------------
%
\chapter{Mathematical Model Description}
%
%----------------------------
%
\label{chap:mmd}
%
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 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. 
%
%----------------------------
\section{Equations}
%----------------------------
%
%
\label{subsec:sec_eq}
%
The equations we use 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}
		\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}
\item Water momentum conservation in $y$ direction:
	\begin{equation}
	\label{eq:wmomy_sf}
		\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}
\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 (Exner):
	\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 (Hirano):
	\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$ 					= x coordinate [$m$]
	\item $y$ 					= 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 $I$ = secondary flow intensity [$m/s$]
	\item $D_H$ = diffusion coefficient [$m^2/s$]
	\item $S_s$ = source term [$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*}

%----------------------------
\section{Simplifications}
%----------------------------
%
%
\label{subsec:sec_sim}
%
The active layer thickness is assumed 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}
%
%
%----------------------------
\section{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.

\subsection{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}

\subsection{Sediment Transport Rate}

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_{\tau}},\sin{\varphi_{\tau}}) \quad k\in\{1,N\}
\end{equation}
where:
\begin{itemize}
	\item $\varphi_{\tau}$ = direction of the sediment transport rate [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_{\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}

The constant $\alpha_I$ [$-$] is:
\begin{equation}
\label{eq:alphaI}
	\alpha_I=\frac{2}{\kappa^2}E_s\left(1-\frac{\sqrt{C_f}}{2\kappa}\right)
\end{equation}
where:
\begin{itemize}
	\item $\kappa$ = Von K\'arm\'an constant [$-$] 
	\item $E_s$ = calibration parameter [$-$]
\end{itemize}

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)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 $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$]
	\item $d_k$ = characteristic grain size of size fraction $k$ [$m$]
	\item $A$ = nondimensional parameter [$-$]
	\item $B$ = nondimensional parameter [$-$]
	\item $\theta_k$ = module of the Shields stress of size fraction $k$ [-] 
	\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}
	\item $D_m$ = characteristic mean grain size of the mixture [$m$]
\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}

\subsection{Secondary Flow Terms}

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 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}


\subsection{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}
%  \end{linenomath*}
%
%----------------------------
%
%
\chapter{Model Characterization}
%
%----------------------------
%
\label{chap:mc}
%
In this section we describe the assessment of the model comprised by the Shallow Water Equations coupled to the active layer model \citet{Hirano71}. In Section \ref{sec:hyp2D} we introduce the necessary mathematical concepts. To asses the model we decompose the equations' terms into spatial gradients of the dependent variables (Section \ref{subsec:sec_ed}). In Section \ref{subsec:sec_mv} we present the matrix-vector formulation of the model. 

\section{Hyperbolicity of a 2D System of Equations}
\label{sec:hyp2D}

We first assume a system of equations without diffusive terms. In a 1D problem the equations can be recast in a matrix form:
\begin{equation}
\label{eq:matrixf_1D}
\frac{\partial \mathbf{Q}}{\partial t}+\mathbf{A}\frac{\partial \mathbf{Q}}{\partial x}=\mathbf{S}\;.
\end{equation}
%  \end{linenomath*}
This is the one-dimensional quasi-linear non-conservative form of the advection equation. $\mathbf{Q}$ is the vector of dependent variables, $\mathbf{A}$ is the system matrix, and $\mathbf{S}$ is the vector of source terms. $t$ [$s$] denotes the time coordinate, $x$ [$m$] the streamwise coordinate.

The problem in Equation \ref{eq:matrixf_1D} will be well-posed if wave-like solutions exist. This is, if the $n$ system equations give rise to $n$ waves propagating along the domain. This property is called hyperbolicity. A system is hyperbolic or parabolic at a point $(x,t)$ if all the eigenvalues of matrix $\mathbf{A}$ are real. Physical propagation problems are modelled with hyperbolic or parabolic systems of equations. If all the eigenvalues are complex, the system is said to be elliptic. Elliptic systems model equilibrium physical problems. If matrix $\mathbf{A}$ has both real and complex eigenvalues it is a mixed-type system. 

A 2D problem can be recast in a matrix form:
\begin{equation}
\label{eq:matrixf_2D}
\frac{\partial \mathbf{Q}}{\partial t}+\mathbf{A_x}\frac{\partial \mathbf{Q}}{\partial x}+\mathbf{A_y}\frac{\partial \mathbf{Q}}{\partial y}=\mathbf{S}
\end{equation}
where $\mathbf{A_x}$ and $\mathbf{A_y}$ are the system matrices in $x$ and $y$ direction.

For a 2D problem to be hyperbolic, wavelike solutions must exist for any arbitrary direction specified by the unit vector $\vec{n}=(n_x,n_y)$. Thus, the problem specified by Equation \ref{eq:matrixf_2D} is said to be hyperbolic if the matrix:
\begin{equation}
\label{eq:A2D}
\mathbf{A_{2D}}=\mathbf{A_x}n_x+\mathbf{A_y}n_y
\end{equation}
diagonalizes with real eigenvalues. In this case, the derivative of the eigenvalues with respect to $n_x$ and $n_y$ yield the eigenvalues in the $x$-$y$ domain. These eigenvalues form the intersection of the Monge cone(s) with a unit time (see e.g.\ \citet{Sloff92,Sieben94}). Thus the Monge cone is a visualization of the propagation celerity of the information in the domain. If at least one cone is not real the system is elliptic. 
%
%
%----------------------------
\section{Decomposition of the System Equations}
%----------------------------
%
%
\label{subsec:sec_ed}
%
In this section we decompose the divergence terms of the model equations.

\subsection{Decomposition of the Sediment Transport Rate}

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}) \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}  \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}) \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} \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 model 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.

\subsection{Secondary Flow Terms}

The secondary flow t, 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)$.

\subsection{Substitution and Expansion}

In this section we substitute the decomposed divergence terms in the model equations and we expand the terms to obtain products of the gradient of the dependent variables. 

All the terms in the water mass equation, Equation (\ref{eq:wmass_sf}), are already products of the gradient of the dependent variables. 

Expanding the terms of the equation for the conservation of water momentum in the $x$ direction, Equation (\ref{eq:wmomx_sf}), we obtain: 
\begin{equation}
\label{eq:wmomx_sf_s}
\begin{split}
	\frac{\partial q_x}{\partial t}+\left(gh-\left(\frac{q_x}{h}\right)^2\right)\frac{\partial h}{\partial x}-\frac{q_xq_y}{h^2}\frac{\partial h}{\partial y}\\
	+\left(2\frac{q_x}{h}-\frac{\partial T^\prime_{xx}}{\partial q_x}\right)\frac{\partial q_x}{\partial x}+\left(\frac{q_y}{h}-\frac{\partial T^\prime_{xy}}{\partial q_x}\right)\frac{\partial q_x}{\partial y}\\
	-\frac{\partial T^\prime_{xx}}{\partial q_y}\frac{\partial q_y}{\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_{xx}}{\partial I}\frac{\partial I}{\partial x}-\frac{\partial T^\prime_{xy}}{\partial I}\frac{\partial I}{\partial y}+gh\frac{\partial \eta}{\partial x}=-ghS_{fx}
\end{split}
\end{equation}

Expanding the terms of the equation for the conservation of water momentum in the $y$ direction, Equation (\ref{eq:wmomy_sf}), we obtain: 
\begin{equation}
\label{eq:wmomy_sf_s}
\begin{split}
	\frac{\partial q_y}{\partial t}-\frac{q_xq_y}{h^2}\frac{\partial h}{\partial x}+\left(gh-\left(\frac{q_y}{h}\right)^2\right)\frac{\partial h}{\partial y}\\
	+\left(\frac{q_y}{h}-\frac{\partial T^\prime_{yx}}{\partial q_x}\right)\frac{\partial q_x}{\partial x}-\frac{\partial T^\prime_{yy}}{\partial q_x}\frac{\partial q_x}{\partial y}\\
	+\left(\frac{q_x}{h}-\frac{\partial T^\prime_{yx}}{\partial q_y}\right)\frac{\partial q_y}{\partial x}+\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_{yx}}{\partial I}\frac{\partial I}{\partial x}-\frac{\partial T^\prime_{yy}}{\partial I}\frac{\partial I}{\partial y}+gh\frac{\partial \eta}{\partial y}=-ghS_{fy}
\end{split}
\end{equation}

All the terms in the constitutive equation of the secondary flow intensity, Equation (\ref{eq:consI}), are already products of the gradient of the dependent variables. 

Substitution of equations \ref{eq:qsfuncx_sf}, and \ref{eq:qsfuncy_sf} in \ref{eq:exner_sf} yields the Exner equation:
\begin{equation}
\label{eq:exner_sf_s}
\begin{split}
	\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}=0
\end{split}
\end{equation}

Substitution of equations \ref{eq:qsfuncx_sf}, and \ref{eq:qsfuncy_sf} in \ref{eq:hirano} yields the Hirano equation:
\begin{equation}
\label{eq:hirano_sf_s}
\begin{split}
	\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_{bky}}{\partial h}-f^{I}_{k}\frac{\partial q_{by}}{\partial h}\right)\frac{\partial h}{\partial y}+\left(\frac{\partial q_{bky}}{\partial q_x}-f^{I}_{k}\frac{\partial q_{by}}{\partial q_{x}}\right)\frac{\partial q_x}{\partial x} \\
	+\left(\frac{\partial q_{bky}}{\partial q_{y}}-f^{I}_{k}\frac{\partial q_{by}}{\partial q_y}\right)\frac{\partial q_y}{\partial y}+\left(\frac{\partial q_{bky}}{\partial I}-f^{I}_{k}\frac{\partial q_{by}}{\partial I}\right)\frac{\partial I}{\partial y}\\
	+\sum_{l=1}^{N-1}\left(\frac{\partial q_{bky}}{\partial M_{al}}-f^{I}_{k}\frac{\partial q_{by}}{\partial M_{al}}\right)\frac{\partial M_{al}}{\partial y}=0 \quad k\in\{1,N-1\} 
\end{split}
\end{equation}

%----------------------------
\section{Matrix-Vector Formulation}
%----------------------------
%
%
\label{subsec:sec_mv}
%
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}$. 

We recast the model Equations, (\ref{eq:wmass_sf}), (\ref{eq:wmomx_sf_s}), (\ref{eq:wmomy_sf_s}), (\ref{eq:consI}), (\ref{eq:exner_sf_s}), and (\ref{eq:hirano_sf_s}) in matrix form:
\begin{equation}
\label{eq:matrixf_sf}
	\frac{\partial \mathbf{Q}}{\partial t}+\mathbf{A_x}\frac{\partial \mathbf{Q}}{\partial x}+\mathbf{A_y}\frac{\partial \mathbf{Q}}{\partial y}+\mathbf{D_x}\frac{\partial^2 \mathbf{Q}}{\partial x^2}+\mathbf{D_y}\frac{\partial^2 \mathbf{Q}}{\partial y^2}+\mathbf{B}=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 \\
	\hline
  [M_{ak}]  
 \end{bmatrix}
\end{equation}


%\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|c}
%wmass 
  0 & 1 & 0 & 0 & 0 & \left[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[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[0\right] \\
%const I
0 & 0 & 0 & \frac{q_x}{h} & 0 & \left[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
	\hline
	\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[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|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[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[0\right] \\
%const I
0 & 0 & 0 & \frac{q_y}{h} & 0 & \left[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
	\hline
	\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[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 diffusive matrix in $x$ and $y$ direction are:
%\begin{adjustwidth}{-4cm}{-3cm}
\begin{equation}
\label{eq:D_x}
		\mathbf{D_x}=\mathbf{D_y}=
 \begin{bmatrix}
  0 & 0 & 0 & 0 & 0 & \left[0\right] \\
	0 & 0 & 0 & 0 & 0 & \left[0\right] \\
	0 & 0 & 0 & 0 & 0 & \left[0\right] \\
	0 & 0 & 0 & D_H & 0 & \left[0\right] \\
	0 & 0 & 0 & 0 & 0 & \left[0\right] \\
	\hline
	\left[0\right] & \left[0\right] & \left[0\right] & \left[0\right] & \left[0\right] & \left[0\right] \\
\end{bmatrix}
\end{equation}
%\end{adjustwidth}

The vector of source terms is:
\begin{equation}
\label{eq:B_sf}
		\mathbf{B}=
 \begin{bmatrix}
  0 \\
	-ghS_{fx} \\
	-ghS_{fy} \\
	S_s \\
	0 \\
	\hline
	\left[0\right] \\
 \end{bmatrix}
\end{equation}

\section{Implementation in Delft3D}
In this section a short description of the implementation of the Hirano check in Delft3D is given. For details about the files and the generated output see Appendix \ref{app:changeDelft}.

In the current implementation in Delft3D the entries of the Jacobian matrices, Equations (\ref{eq:Ax_sf}) and (\ref{eq:Ay_sf}), have been computed numerically. However, the effect of secondary flow and gravitational pull on the sediment in transport have been neglected. 

The entries of the Jacobian in top left block (3x3 related to $h$, $q_x$ and $q_y$) are directly computed based on the analytical expressions presented in Equations \ref{eq:Ax_sf} and \ref{eq:Ay_sf}. The entries of the Jacobian of the other blocks (related to $\eta$ and $M_{ak}$) contain derivatives of the sediment transport rate with respect to $h$, $q_x$, $q_y$, and $M_{ak}$. These derivatives are obtained using using a finite difference approach:
\begin{equation}
\label{eq:diff}
	\frac{\partial q_{bkx}}{\partial f}=\frac{q_{bkx}\left(f+\mathrm{d}f\right)-q_{bkx}\left(f\right)}{\mathrm{d}f} \;,
\end{equation}
where $f$ is one of the dependent variables. To compute these values many calls are made to the sediment transport computation without bed or composition update. For each variable a perturbation ($\mathrm{d}f$) is introduced. This leads to the following sequence of calls (only the resulting $x$-direction is shown here, but the $y$-direction is also computed in the same loop):

\begin{itemize}
\item Call 1: perturbation to the first dependent variable:\\
$h+\mathrm{d}h$, $q_x$, $q_y$, $\eta$, $M_{a,1}$, $\ldots$, $M_{a,N-1}$ \\
\indent $\rightarrow$ $q_{bx}(h+\mathrm{d}h)$, $q_{bx,1}(h+\mathrm{d}h)$ $\ldots $ $q_{bx,N-1}(h+\mathrm{d}h)$\\
\item Call 2: perturbation to the second dependent variable:\\ 
$h$, $q_x+\mathrm{d}q_x$, $q_y$, $\eta$, $M_{a,1}$, \ldots, $M_{a,N-1}$ \\
\indent $\rightarrow$ $q_{bx}(q_x+\mathrm{d}q_x)$, $q_{bx,1}(q_x+\mathrm{d}q_x)$ $\ldots $ $q_{bx,N-1}(q_x+\mathrm{d}q_x)$\\
$\vdots$\\
\\
\item Call $N+3$: perturbation to the $N+3$ (and last) dependent variable:\\ 
$h$, $q_x$, $q_y$, $\eta$, $M_{a,1}$, \ldots, $M_{a,N-1}+\mathrm{d}M_{a,N-1}$ \\
\indent $\rightarrow$ $q_{bx}(M_{a,N-1}+\mathrm{d}M_{a,N-1})$, $q_{bx,1}(M_{a,N-1}+\mathrm{d}M_{a,N-1})$ $\ldots $ $q_{bx,N-1}(M_{a,N-1}+\mathrm{d}M_{a,N-1})$\\
\item Call $N+4$: no perturbation, update of bed elevation and composition:\\
$h$, $q_x$, $q_y$, $\eta$, $M_{a,1}$, \ldots, $M_{a,N-1}$ \\
\indent $\rightarrow$ $q_{bx}$, $q_{bx,1}$ $\ldots $ $q_{bx,N-1}$\\
\end{itemize}

To compute the partial derivatives the differentials need to modify only one of the dependent variables at a time. Note that we have chosen to use as dependent variables $h$, $q_x$, $q_y$, $\eta$, and $M_{ak}$ for $k=1,N$. Thus, when a differential is added to, for instance, the flow depth $h$, we need to impose that there is no change in the specific water discharge in $x$ and $y$ direction by adjusting the mean flow velocity. 

In the case of the differentials added to the volume of each sediment fraction in the active layer the constrain that the addition of the volume fraction content of all the sediment in the active layer is equal to 1 (Equations (\ref{eq:Fcons}) and (\ref{eq:Mcons})) needs to hold. To this end the volume fraction of sediment in the active layer after one size fraction is perturbed (superscript $'$) is added to compute the volume of sediment of the last fraction $N$ as: 
\begin{equation}
\label{eq:MaN}
	 M_{aN}^{'}=L_a-\sum_{k=1}^{N-1} M_{ak}^{'} \;.
\end{equation}

A series of tests were consequently run during the development after every change to the code, to ensure that the implemented code does not interfere with the normal computation loop. During the development no changes to the outcome of the standard Delft3D simulation were made.


\chapter{Results}
\label{chap:pr}
%
\section{2D Eigenvalue Analysis}
\label{sec:2dea}
%
In this section we describe the eigenvalue analysis performed to be able to check the implementation of the ellipticity check tool in Delft3D. Once the function is implemented in Delft3D we will compare the eigenvalues obtained by that function to the ones obtained by means of a symbolic solver. The use of the symbolic solver allows us to obtain the eigenvalues without introducing the complexities of the interpolation and the numerical computation of the necessary derivatives in Delft3D. Moreover, we validate the symbolic results with previous analysis and analytical solutions.

Using the symbolic solver of Matlab we obtain a symbolic expression of the system matrix $\mathbf{A}_{2D}$, Equation \ref{eq:A2D}, for 4 simplified cases:
\begin{enumerate}
\item The Shallow Water Equations without secondary flow correction
\item The Shallow Water Equations without secondary flow correction coupled to the \citet{Exner20} model (unisize)
\item The Shallow Water Equations without secondary flow correction coupled to the \citet{Hirano71} model (2 size fractions)
\item The Shallow Water Equations without secondary flow correction coupled to the \citet{Hirano71} model (3 size fractions)
\end{enumerate}

We obtain the symbolic expressions of eigenvalues of $\mathbf{A}_{2D}$, Equation \ref{eq:A2D}, for each case. These expressions are symbolically derived with respect to $n_x$ and $n_y$ to obtain the eigenvalues in $x$ and $y$ coordinates. The evaluation of these last eigenvalues for certain values of unitary vectors covering all the space domain yield the Monge cone associated with each system. 

To validate our method we compare the Monge cones obtained using our method to analytical solutions and previous results. 

We first consider a case in which $u=1$ $m/s$, $v=0$ $m/s$, and $h=1$ $m$. In Figure \ref{fig:monge}a we plot the intersection of the Monge cone of the Shallow Water Equations at $t=1$ $s$ for this specific case. Note that because we do not consider the non-linearities of the system the time at which we plot the Monge cone is irrelevant. This solution is the same as obtained analytically \citep{Sloff92}. 
%
\begin{figure}
	\centering
	\includegraphics[width=1.0\textwidth]{validation_pictures/monge.eps}
%	\makebox[\textwidth][c]{\includegraphics[width=0.8\textwidth]{testcase_pictures/PannerdenscheKop.png}}%
	\caption{Intersection of the Monge cone at $t=1$ $s$ for the Shallow Water equations in combination with: (\textbf{a}) a fixed bed, (\textbf{b}-\textbf{c}) the \citet{Exner20} equation, (\textbf{d}-\textbf{e}) the \citet{Hirano71} equation considering 2 sediment size fractions in aggradational conditions, (\textbf{f}-\textbf{g}) the \citet{Hirano71} equation considering 2 sediment size fractions in degradational conditions, and (\textbf{h}-\textbf{i}) the \citet{Hirano71} equation considering 3 sediment size fractions in aggradational conditions.}
	\label{fig:monge}
\end{figure}

We consider the same situation with a mobile bed composed of one sediment size fraction with characteristic grain size equal to 0.001 $m$ and a dimensionless friction coefficient equal to 0.007. The sediment transport is computed using the relation by \citet{Engelund67}. In Figure \ref{fig:monge}b we plot the intersection of the Monge cones of the Shallow Water Equations coupled to the \citet{Exner20} equation at $t=1$ $s$ for this specific case. Note that a new cone appears close to the origin. A zoom in this cone is plotted in Figure \ref{fig:monge}c. This result is the same as obtained in \citep{Sloff92}. 

We consider the previous situation in which, rather than one single size fraction, we have two size fractions with characteristic grain sizes equal to 0.001 and 0.002 $m$. The volume fraction content of both size fractions in the active layer is equal to 0.5. The volume fraction content of both size fractions at the interface between the active layer and the substrate is equal to 0.5. Thus, this situation represents an aggradational case. The active layer thickness is equal to 0.1 $m$. In Figure \ref{fig:monge}d we plot the intersection of the Monge cones of the Shallow Water Equations coupled to the \citet{Hirano71} equation at $t=1$ $s$ for this specific case. A new cone appears close to the origin which can only be appreciated if we zoom in (Figure \ref{fig:monge}e). This result is the same as obtained in \citet{Sieben94}. 

We consider a situation with the same parameters as before, yet the active layer is composed of the coarse size fraction only and the volume fraction content of fine sediment at the interface between the active layer and the substrate is equal to 1. Thus, this situation represents a degradational case into a fine substrate. In Figure \ref{fig:monge}g we plot a zoom of the intersection of the Monge cone of the Shallow Water Equations coupled to the Hirano equations at $t=1$ $s$ for this specific case. In this case the cone does not pertain to the real numbers. Thus, this situation is elliptic and we can only plot the real part of the cone.  

Eventually, we consider the previous situation in which, rather than two size fractions, we have three size fractions with characteristic grain sizes equal to 0.001, 0.002, and 0.004 $m$. The volume fraction content all size fractions in the active layer is equal to 0.33. The volume fraction content of all size fractions at the interface between the active layer and the substrate is equal to 0.33. Thus, this situation represents an aggradational case. The sediment transport rate is computed using the relation by \citet{MeyerPeter48} with the hiding correction by \citet{Parker82_3} with parameter $b$=0.8. In Figure \ref{fig:monge}h we plot the intersection of the Monge cones of the Shallow Water Equations coupled to the \citet{Hirano71} equation at $t=1$ $s$ for this specific case. A new cone appears close to the origin which can only be appreciated if we zoom in (Figure \ref{fig:monge}i).

\section{Implementation Test}
\label{sec:impl_test}

In this section we test the implementation of the routine in Delft3D. To this end we compare the matrices given by Deft3D to the ones obtained in the external Matlab algorithm. We test the implementation in three short normal flow simulations which reproduce the last three cases described in Section \ref{sec:2dea}.

The domain is 1 $m$ long and 1 $m$ wide, discretized in square cells of side 0.5 $m$. The boundaries at $x=0$ $m$ and $x=1$ $m$ are open. The boundaries at $y=0$ $m$ and $y=1$ $m$ are closed. To obtain a mean flow velocity in the $x$ direction equal to 1 $m/s$ and 0 $m/s$ in the $y$ direction in combination with a mean flow depth equal to 1 $m$, the water discharge is equal to 1 $m^3/s$. The water level at the downstream end is such that the flow depth is equal to 1 $m$. A tiny lowering of 0.001 $m/min$ is imposed to guarantee that there is a minimal degradation so that the volume fraction content at the interface between the active layer and the substrate is equal to the volume fraction content at the substrate. The nondimensional friction coefficient is equal to 0.007 (Ch\'ezy equal to 37.44 $m^{1/2}/s$). The active layer thickness is equal to 0.1 $m$. The substrate is 1 $m$ deep discretized in two equal layers. Secondary flow is not taken into account. The time step is equal to 0.1 $s$ (which guarantees a small enough CFL number). The simulations last for 60 $s$.  

The matrices in Delft3D are obtained at the node $n=2$ and $m=2$ at $t=2$ $s$. 

In Simulation 1, the sediment mixture is discretized in two sediment size fractions with characteristic grain sizes equal to 0.001 and 0.002 $m$. The sediment transport rate is computed using the \citet{Engelund67} sediment transport relation. The sediment input at the upstream end of the fine and coarse size fractions is equal to $5.588\cdot10^{-5}$ and $2.794\cdot10^{-5}$ $m^3/s$. This boundary condition is in equilibrium with a bed slope equal to $7.13541\cdot10^{-4}$ and a volume fraction content of the fine sediment in the active layer equal to 0.5. The volume fraction content of fine sediment in the substrate is equal to 0.5. 

Under these conditions, matrix $\mathbf{A}_x$, Equation \ref{eq:Ax_sf}, and $\mathbf{A}_y$, Equation \ref{eq:Ay_sf}, (assuming no secondary flow) computed in Matlab are:

\begin{equation}
\label{eq:Ax_S1}
\makebox[\textwidth][c]{$
		\mathbf{A}_{x1}=\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}_{y1}=\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}

The same matrices computed in Delft3D are:

\begin{equation}
	\label{eq:Ax_S1_Delft3D}
	\makebox[\textwidth][c]{$
		\mathbf{A}_{x1,\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}_{y1,\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}

The absolute value of the relative error in the norm of matrices $\mathbf{A}_x$ and $\mathbf{A}_y$ in Simulation 1 is equal to $2.3617\cdot10^{-5}$ and $1.8428\cdot10^{-5}$ respectively. 

In Simulation 2, the sediment mixture is discretized in two sediment size fractions with characteristic grain sizes equal to 0.001 and 0.002 $m$. The sediment transport rate is computed using the \citet{Engelund67} sediment transport relation. The sediment input at the upstream end of the fine and coarse size fractions is equal to $0$ and $5.588\cdot10^{-5}$ $m^3/s$. This boundary condition is in equilibrium with a bed slope equal to $7.13541\cdot10^{-4}$ and a volume fraction content of the fine sediment in the active layer equal to 0. The volume fraction content of fine sediment in the substrate is equal to 1. 

Under these conditions, matrix $\mathbf{A}_x$, Equation \ref{eq:Ax_sf}, and $\mathbf{A}_y$, Equation \ref{eq:Ay_sf}, (assuming no secondary flow) computed in Matlab are:

\begin{equation}
\label{eq:Ax_S2}
\makebox[\textwidth][c]{$
		\mathbf{A}_{x2}=\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}_{y2}=\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}

The same matrices computed in Delft3D are:

\begin{equation}
	\label{eq:Ax_S2_Delft3D}
	\makebox[\textwidth][c]{$
		\mathbf{A}_{x2,\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}_{y2,\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}

The absolute value of the relative error in the norm of matrices $\mathbf{A}_x$ and $\mathbf{A}_y$ in Simulation 2 is equal to $2.3568\cdot10^{-5}$ and $1.8331\cdot10^{-5}$ respectively. 
						
In Simulation 3, the sediment mixture is discretized in three sediment size fractions with characteristic grain sizes equal to 0.001, 0.002, and 0.004 $m$. The sediment transport rate is computed using the \citet{MeyerPeter48} sediment transport relation with the hiding correction by \citet{Parker82_3} with parameter $b$ equal to 0.8 (Equation (\ref{eq:xi_parker})). The sediment input at the upstream end of the fine, medium, and coarse size fractions is equal to $6.723\cdot10^{-5}$, $6.319\cdot10^{-5}$, and $5.865\cdot10^{-5}$ $m^3/s$. This boundary condition is in equilibrium with a bed slope equal to $7.13541\cdot10^{-4}$ and a volume fraction content of the fine and medium fractions in the active layer equal to $1/3$. The volume fraction content of fine and medium fractions in the substrate is equal to $1/3$. 

Under these conditions, matrix $\mathbf{A}_x$, Equation \ref{eq:Ax_sf}, and $\mathbf{A}_y$, Equation \ref{eq:Ay_sf}, (assuming no secondary flow) computed in Matlab are:

\begin{equation}
\label{eq:Ax_S3}
\makebox[\textwidth][c]{$
		\mathbf{A}_{x3}=\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}_{y3}=\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}

The same matrices computed in Delft3D are:

\begin{equation}
	\label{eq:Ax_S3_Delft3D}
	\makebox[\textwidth][c]{$
		\mathbf{A}_{x3,\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}_{y3,\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}

The absolute value of the relative error in the norm of matrices $\mathbf{A}_x$ and $\mathbf{A}_y$ in Simulation 3 is equal to $2.3762\cdot10^{-5}$ and $1.8525\cdot10^{-5}$ respectively. 

\section{2D test case}
For the rapid development of the check of the well-posedness in the Hirano system, a real world test case near the Pannerdensche Kop is used. At this location the Bovenrijn bifurcates into the Waal heading west and the Pannerdensch Kanaal towards the north (see Figure \ref{fig:pannerdensch_kop}). 

\begin{figure}
	\centering
	\includegraphics[width=0.8\textwidth]{testcase_pictures/PannerdenscheKop.png}
%	\makebox[\textwidth][c]{\includegraphics[width=0.8\textwidth]{testcase_pictures/PannerdenscheKop.png}}%
	\caption{Pannerdensch Kop bifurcation (\copyright Google Earth).}
	\label{fig:pannerdensch_kop}
\end{figure}

The amount of sand transported into the Waal is about 8 times the amount of sand transported into the Pannerdensch Kanaal and the amount of gravel transport into the Waal branch equals 2 times the gravel transport into the Pannerdensch Kanaal \citep{Frings14_2}. This situation is therefore a suitable test for testing the well-posedness of the Hirano active layer equations.

Due to the fact that the implementation is not yet finished we cannot check this simulation. 

\chapter{Discussion and Further Research}
\label{chap:disc}

The Jacobians computed in Delft3D and in MATLAB are similar enough. A source of error is due not only to the fact that the derivatives in Delft3D are obtained with finite differences, but also due to the fact that the data used is not exactly the same. The values of the dependent variables in Delft3D differ from the imposed ones in the analytical computation because the simulation is not exactly a normal flow simulation. Thus, part of the difference between the analytical and the Delft3D solution is not error. 

We have only tested the implementation in 3 cases which do not include all the generalities. For instance, we have to check with flow in both $x$ and $y$ or the effect of a curvilinear grid.

At this moment the only output are the Jacobians. It is necessary to implement the eigenvalues computation of the directional Jacobian, Equation \ref{eq:A2D}, to obtain as an output the character of the system (i.e.\ hyperbolic or elliptic).
 
At this moment we use the sediment transport at the edges rather than at the cell centers. This could be a source of errors for more complicated domains. Furthermore, a more consistent implementation would use the sediment transport at the cell centers (or the other variables at the cell edges).

In the current implementation the volume fraction content at the interface between the active layer and the substrate is obtained from the second layer of the Delft3D data structure. This second layer correctly corresponds to the top substrate layer if the bed degradation is smaller than the initial thickness of the second layer. A loop that identifies the top substrate layer in a general case needs to be implemented.

Secondary flow has not been studied. The matrix-vector formulation has not yet been checked and the implementation is not yet performed. 

The role of the gravitational pull effects on the sediment transport has also not been studied. A thorough analytical study needs to be conducted before anything is implemented.

\chapter{Conclusion}
\label{chap:conclusions}

In the present report we present a mathematical analysis of the active layer model \citep{Hirano71} used in mixed-sediment river morphodynamics in combination with the Shallow Water Equations. This is relevant to assess the well-posedness of 2D simulations.

The criteria for well-posedness without considering the effects of secondary flow and slope effects have been formalized based on the continuous equations of motion. It is necessary to further study the effect of the neglected terms.

We obtain the system matrices using a symbolic solver. This allows us to plot the Monge cones to study the effect of the different model parameters. Moreover, this analytical solution allow us to validate the implementation of a well-posedness routine in Delft3D.

We have implemented a routine in Delft3D to characterize the system of equations (i.e.\ ill- or well-posed). The implementation appears to be correct. Yet, more tests need to be conducted. 

A test case near the Pannerdensch Kanaal is created, which can be used to assess the occurrence of the elliptic behavior. When the eigenvalue computation has been validated it is recommended to research this graded sediment model for the occurrence of elliptic behavior, and what happens in the Delft3D solution when such a circumstance is encountered. 
%
%	
%
%----------------------------
%---------SCHEDULE
%----------------------------
%
%\chapter{Schedule}
%\includepdf[angle=90]{schedule.pdf}
%
%
%----------------------------
%---------REFERENCES
%----------------------------
%
%\DeclareRobustCommand{\van}[3]{#3}
%\bibliographystyle{apalike}
\chapter{References}
\bibliography{references}
%
%
%----------------------------
%---------APPENDIX
%----------------------------
%
%
%\appendix
%%%
%%%%
%\chapter{Notation}
%%
%\begin{itemize}
%\item $x$ 					= streamwise 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 $\eta_{0}$		= reference bed elevation [$m$]
%\item $F_{ak}$ 			= volume fraction of size fraction $k$ in the active layer [$-$]
%\item $f^{I}_{k}$ 			= volume fraction of size fraction $k$ at the interface    [$-$]
%\item $f_{sk}$ 			= volume fraction of size fraction $k$ in 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_{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_{f}$   		= friction slope [$-$]
%\item $p$ 					= porosity [$-$]
%\item $g$ 					= gravity constant [$m/s^2$]
%\item $Fr$ 					= Froude number [$-$]
%\item $C_f$					= dimensionless friction coefficient [$-$]
%\end{itemize}
	

%\begin{itemize}
%\item $M_{ak}$ 			= $F_{ak}L_{a}$
%\item $M_{sk}$ 			= $\int_{\eta_0}^{\eta-L_{a}}f_{sk}(z)dz$
%\item $Q_{bk}$			= $\frac{1}{1-p}\hat{Q_{bk}}$
%\item $q_{bk}$			= $F_{ak}Q_{bk}$
%\item $\psi$  			= $\frac{\partial {q_b}}{\partial q}$
%\item $\gamma_{1}$ 	= $c_{1}-f^{I}_{1}$
%\item $c_{1}$ 			= $\frac{1}{\psi}\frac{\partial q_{b1}}{\partial q_{b}}$
%\item $\vec{u}$			= $\sqrt{u^2+v^2}$
%%\item $\kappa_{1}$	= $\gamma{1}+c_{1}$
%\end{itemize}
%
%\chapter{Consistency of the Regularization Strategy}
%\label{app:crs}
%%
%We write the \citet{Exner20} equation with the modifying factor $\beta$:
%\begin{equation}
%\label{eq:exnerm}
	%\beta\frac{\partial \eta}{\partial t}+\frac{\partial q_b}{\partial x}=0
%\end{equation}
%and the active layer equation with the $\alpha_k$ and $\beta$ parameters:
%\begin{equation}
%\label{eq:hiranom}
	%\beta\alpha_k\frac{\partial M_{ak}}{\partial t}+f^I_k\beta\frac{\partial (\eta-L_a)}{\partial t}+\frac{\partial q_{bk}}{\partial x}=0 %\quad \text{for } 1\leq k \leq N-1
%\end{equation}
%where $q_{bk}$ is the sediment transport of size fraction $k$ per unit width including pores [$m^2/s$] and $f^I_k$ the volume of size fraction $k$ at the interface between the active layer and the substrate [$-$]. Note the factor $\beta$ in the the time derivative of the interface elevation in the active layer equation. This is consistent with our definition in Equation \ref{eq:matrixm} where the spatial derivatives are not modified. We add all ($N$) the active layer equations and substitute the \citet{Exner20} equation obtaining:
%\begin{equation}
%\label{eq:condition}
	%\beta\left[\frac{\partial L_a}{\partial t}\left(\sum_{k=1}^{N}\alpha_kF_{ak}-1\right)+L_a\frac{\partial \sum_{k=1}^{N}\alpha_kF_{ak}}{\partial t}\right]=0
%\end{equation}
%

%----------------------------
%---------END OF DOCUMENT
%----------------------------
\appendix
\chapter{Changes to Delft3D}
\label{app:changeDelft}

A branch of the Delft3D modelling system was created for the development of the approximate eigenvalue analysis \citep{url:Delft3Dbranch}. Using a branch allows the implementation of new concepts, whilst ensuring the stability of the main development line. The following files have been updated for the implementation of the check of the well-posedness of the active layer description of \citet{Hirano71}, and subsequently determine the possibilities to circumvent the occurrence of elliptic behavior. 

\begin{verbatim}
engines_gpl/flow2d3d/packages/data/include/inttim.igs
engines_gpl/flow2d3d/packages/data/src/basics/globaldata.f90
engines_gpl/flow2d3d/packages/io/src/input/rdsedmortra.f90
engines_gpl/flow2d3d/packages/io/src/output/wrsedm.f90
engines_gpl/flow2d3d/packages/kernel/kernel.vfproj
engines_gpl/flow2d3d/packages/kernel/src/compute_sediment/erosed.f90
engines_gpl/flow2d3d/packages/kernel/src/compute_sediment/hirano_check.f90
engines_gpl/flow2d3d/packages/kernel/src/compute_sediment/hirano_check_init.f90
engines_gpl/flow2d3d/packages/kernel/src/compute_sediment/hirano_check_perturbation.f90
engines_gpl/flow2d3d/packages/kernel/src/compute_sediment/Makefile.am
engines_gpl/flow2d3d/packages/kernel/src/compute_sediment/rg_dependencies.f90
engines_gpl/flow2d3d/packages/kernel/src/main/trisol.f90
engines_gpl/flow2d3d/packages/manager/src/tricom_step.F90
utils_gpl/morphology/packages/morphology_data/src/morphology_data_module.f90
utils_gpl/morphology/packages/morphology_io/src/rdmor.f90
utils_gpl/morphology/packages/morphology_io/src/rdstm.f90
utils_gpl/morphology/packages/morphology_kernel/src/bedcomposition_module.f90
\end{verbatim}

In addition output fields 'HIRCHK', 'HIRJCU', 'HIRJCV' are added to 'map-sed-series' in the trim-file. 
The field 'HIRCHK' is used to indicate whether the approximate eigenvalues at the time step previous to writing the map file are elliptic or not. This is indicated using a 1 for elliptic and 0 for not elliptic. 
The fields 'HIRJCU', 'HIRJCV' represent the approximate Jacobian in $u$ and $v$ direction.
This check can be activated using the 'HiranoCheck' option in the morphology input file (with file extension .mor). For ease of use, the output has been added in such a way that it can be viewed and loaded into Matlab memory using Delft3D-QuickPlot. Furthermore the perturbation can be set using the keyword 'HiranoCheckPertubation'
\begin{verbatim}
[Numerics]
    HiranoCheck = 1     [ - ] Flag for well-posedness of Hirano (1971) check
                                0 (default): Off
                                1          : On
    HiranoCheckPertubation = 0.0001     # (default 0.0001)                            
\end{verbatim}





%----------------------------
%---------END OF DOCUMENT
%----------------------------

%------------------------------------------------------------------------------
\end{document}

