
\documentclass[nhessd, online, hvmath]{copernicus}

\begin{document}\hack{\sloppy}

\title{TITAN2F: a pseudo-3-D  model of 2-phase debris flows}

\Author[1]{G.}{C\'ordoba}
\Author[2]{M.~F.}{Sheridan}
\Author[3]{E.~B.}{Pitman}

\affil[1]{Universidad de Nari\~no, Pasto, Colombia}
\affil[2]{Center for Geohazards Studies and Department of Geology, University of Buffalo, Buffalo, NY 14260, USA}
\affil[3]{Center for Geohazards Studies and Department of Mathematics, University of Buffalo, Buffalo, NY 14260, USA}

\runningtitle{A model for two phase debris flows}
\runningauthor{G.~C\'ordoba et~al.}

\correspondence{G.~C\'ordoba (gcordobaguerrero@gmail.com)}

\received{15~January 2015}
\accepted{22~April 2015}
\published{}

\firstpage{1}

\maketitle


\begin{abstract}
Debris flows, avalanches, landslides, and other geophysical mass flows can contain
 $O$(10$^6$--10$^{10}$)\,m$^3$ or more of material. These flows commonly consist
 of mixture of soil and rocks with a significant quantity of
interstitial fluid. They can be tens of meters deep, and their runouts can
extend many kilometers.
The complicated rheology of such a mixture challenges every constitutive model
that can reasonably be applied; the range of length and timescales involved in such mass flows
challenges the computational capabilities of existing models.
This paper extends recent efforts to develop a depth
averaged ``thin layer'' model for geophysical mass flows
that contain a mixture of solid material and fluid. Concepts from the engineering community are
integrated with phenomenological findings in geoscience, resulting in a theory that accounts
for the principal solid and fluid forces as well as interactions between the phases, across a wide range of
solid volume fractions. A principal contribution here is  to present drag and phase interaction terms
that comport with the literature in geosciences. The program predicts the evolution of the concentration and dynamic pressure.
The theory is validated with data from one dimensional dam break solutions
and it is verified  with data from artificial channel experiments.
\end{abstract}


\introduction

Globally there are about 50 volcanoes that erupt every year. During the past
century tens of thousands of people have been killed by volcanic flows and
hundreds of thousands forced from their homes
\citep{tilling96,nrc91,nrc94}. Two-phase mass flows containing water and solid
particles, called lahars, are common in volcanic regions. They can be
initiated by several mechanisms. A~volcanic explosion can be accompanied by
large plumes and pyroclastic flows consisting of rock and gas that race along
the surface of the mountain at speeds as high as 100\,\unit{m} per second
\citep{sheridan79}. The hot ash can melt snow, creating a~muddy mixture that
knock down trees and entrain rocks and boulders into the flow. Cotopaxi
Volcano in Ecuador is an example of a~volcano that has produced many large
lahars by this process in the past \citep{pistolesi13}. Crater lakes on
volcanoes are another source of mud flows, a~recent example being the 2007
lahar of Ruapehu in New Zealand \citep{procter10}.  A~third mechanism for
initiating lahars is intense rainfall on hillsides that are devoid of
vegetation and exposures of loose material like clay soils or volcanic ash. An
example of this type of lahar is the 1998 mudflow at Casita Volcano in
Nicaragua that occurred during Hurricane Mitch and caused hundreds of deaths
\citep{sheridan99}. Lahars can carry constituent particles that are typically
from clay to boulder size and can propagate tens of kilometers before coming
to rest \citep{procter10}. As solid particle sediment out of the flowing mass
the resulting deposits can be up to one hundred meters thick
\citep{longrunout02}. However, the typical deposits left after a~debris flow passes are on the scale of meters.


In order to develop a~complete mathematical model of mud and debris flows, two
principal challenges must be overcome: rheology, and scale.
First, constitutive relations must be developed to describe
granular material including clays, sands,
pebbles and rocks, with interstitial water.
Second, a~computational method must be developed that extends over six orders of magnitude, as clay diameters are $O(10^{-6}$ $~m$ and boulders are $O(1)$ $~m$.
Neither of these challenges can be fully met at this time. This paper tries to strike a~balance
between fidelity to the physics of mass flows and computational feasibility.
We describe
a~modeling effort that draws on the wisdom from engineering and geoscience,
to postulate constitutive theory and fluid-solid interaction effects, and, through
a~depth averaging process, results in a~system of equations that is computationally
tractable.

The modeling effort here has its origins the pioneering work of
\citet{savage89}.
They began with mass and momentum balance laws
based on a~Coulomb constitutive description of dry granular material.
By scaling and depth averaging, they
develop a ``thin layer'' model for granular flows down inclines.
Flow over general topography was addressed in \citet{gray99}, \citet{patra05}, and \citet{pudasaini03}.
Comparison of thin layer model results to
historic flows is presented in \citet{sheridan05}.
In \citet{hutter03},
the appropriateness of these thin layer models is considered, for several different types of
geophysical flows. Much of the modeling effort is summarized in \citet{pudasain07}.


Iverson and his co-workers \citep{iverson97, iverson01}
argue that the presence of
interstitial fluid fundamentally
alters the behavior of geophysical flows, and fluid effects
should be included in the constitutive behavior of
the flowing material.
Starting with equations of mixture theory \citep{bedford83}
and through a~careful examination of experiments,
these papers
developed a~system of mass and momentum balance laws for the mixture.
Unfortunately in this development an equation for the motion of pore fluid was not readily available.
Instead, based on experimental data, a~transport equation for the fluid was postulated.

A~different approach, based on a~fully three dimensional model of two phase flows, can be
found in \citet{meruane10} and \citet{dartevelle04}.

\citet{pitman05} rigorously developed a~two phase thin layer model of fluid and
granular material. They begin with a~fully three dimensional model of two phase flows,
based on model equations in engineering \citep{jackson00}. The model equations are scaled
and depth averaged. The resulting system of equations is not complete, and
closure assumptions are required. With these assumptions, the mathematical
system is shown to be hyperbolic
under common conditions, and thus well posed (see \citealp{pelanti08}).
The model of \citet{pitman05} includes a~drag term, which is the only term describing the interaction
of the two phases; the coefficient of drag must be fitted to experiments. That model assumes the fluid is
inviscid, and that there is no frictional dissipation in the fluid phase at the basal surface.  Both of
these features, which are reasonable in bench-scale fluidized bed experiments,
 are suspect for large mass flows. In order to address some of the shortcomings of the \citet{pitman05} model, this paper reconsiders the model equations proposed by them
and proposes
a~revision that better represents two phase geophysical flows, for example accounting for the friction at the wall of the fluid phase and no longer assuming a constant volumetric fraction of solids. Related work can be found
in \citet{valentine89b} and \citet{dobran91}.


A~different approach to modeling mud flows employ a~visco-plastic constitutive assumptions
\citep{mei01, coussot97}; see also \citet{ancey06} and \citet{balmforth99}. In these papers,  the
choice of a visco-plastic flow model drives the subsequent derivation,
as well as the required parameter-fitting necessary for the constitutive relations.
The process of depth averaging a~visco-plastic flow is always difficult. The
interface between yielding and
non-yielding material is itself a~free surface that must be determined. This attribute requires the use of
multiple layers in the model system, with all the resulting complexity.


\section{Model derivation}

This paper uses a~similar  framework to that  developed in
 \citet{pitman05}. However, a~complete set of model equations for a~granular phase and
for a~fluid phase are written. Phase interaction terms are modeled, and scaling of
all terms suggests simplifications that can be made.
Depth averaging and closure assumptions
completes the derivation.

A~note on sign convention: in soil mechanics it is common to consider compressive stresses
as positive; by contrast, in fluid mechanics, as an increase in pressure results in a reduction of the volume,  compression is negative. We caution the reader
to observe the sign convention in the equations below.


\subsection{Fundamental assumptions}

The fundamental theory of two phase flows used here can be found in \citet{dobran91} and \citet{jackson00}.
In two space dimensions we consider a~thin layer of
granular material ($\mathrm{s}$) and interstitial fluid ($\mathrm{f}$), each of constant
specific density $\rho^{\mathrm{s}}$ and $\rho^{\mathrm{f}}$,
respectively, flowing over a~smooth
basal surface, $b$. Erosion and deposition are neglected.
Along the basal surface, we define a~Cartesian coordinate system $Oxyz$, with origin $O$ defined so the
$Oxy$ is tangent to the basal surface, with the $x$~direction the downstream direction, and $Oz$ is in the normal direction.
Writing $\vec{v},\,\vec{u}$ for the velocities of the solid and fluid constituents, respectively, $\varphi$ for the solid volume fraction and $\varphi^f$ for the fluid volume fraction. We assume the mass is fully saturated, so
the sum of the solid and fluid volume fractions adds to one ($\varphi^f = 1-\varphi$). When writing equations in component
form, we use subscripts to denote the component of the vectors, and superscripts the phase of the flow (either solids or fluid).

Mass conservation for the two constituent phases
may be written as in \citet{anderson67}:
\begin{align}
&\partial_t \rho^{\mathrm{s}} \varphi + \nabla \cdot ( \rho^{\mathrm{s}} \varphi \vec{v}) = 0 \label{eq:mass-s-1} \\
&\partial_t  \rho^{\mathrm{f}} (1-\varphi) + \nabla \cdot  ( \rho^{\mathrm{f}} (1-\varphi) \vec{u}) = 0
\end{align}

The momentum equations are:
\begin{align}
&\partial_t (\rho^{\mathrm{s}} \varphi \vec{v} ) + \nabla \cdot (   \rho^{\mathrm{s}}
\varphi   \vec{v} \vec{v})
  = \nabla \cdot  T^{\mathrm{s}} + f^{\mathrm{s}} + \rho^{\mathrm{s}} \varphi \vec{g}  \\
 & \partial_t  (\rho^{\mathrm{f}} (1-\varphi) \vec{u} ) + \nabla \cdot
 (\rho^{\mathrm{f}}  (1-\varphi)      \vec{u} \vec{u}
 = \nabla \cdot  T^{\mathrm{f}} + f^{\mathrm{f}} + \rho^{\mathrm{f}} (1-\varphi) \vec{g}
\end{align}
Here $T^{\mathrm{s}}, \, T^{\mathrm{f}}$ are the stress tensors for the solid and fluid,
and $f^{\mathrm{s}}, \, f^{\mathrm{f}}$  the interaction force between the phases.
We  must postulate constitutive relations and an equation for the interphase force,
to close the system.
\citet{jackson00} presents an argument for separating buoyancy from other
interphase force terms (such as drag or virtual mass), and for properly accounting
for buoyancy in a~field with fluid pressure variations. Similar modeling can be found in
 \citet{neri03}, \citet{neri07}, \citet{dobran01}, \citet{dobran93b} and
  \citet{valentine89b}.
Neglecting capillarity, virtual mass and lift, we postulate
\begin{align}
&f^{\mathrm{s}} = - \varphi \nabla T^{\mathrm{f}} + D (\vec{v} - \vec{u}) \\ \nonumber
&f^{\mathrm{f}} = -f^{\mathrm{s}}
\end{align}
Here the total fluid stress is $T^{\mathrm{f}} = -P^{\mathrm{f}} + \tau^{\mathrm{f}}$,
where $P^{\mathrm{f}}$ is the fluid pressure and $\tau^{\mathrm{f}}$ is the viscous
contribution to the fluid stress.
The drag term exchanges momentum between the phases, with a~coefficient $D$ that is
phenomenological; \citet{wallis69}, \citet{ergun52}, \citet{gidaspow94}, \citet{fan98},
\citet{dobran01}, \citet{panneerselvan07}, and \citet{mazzei07}
among other sources, suggest values.
As $\varphi \rightarrow 1$ the effect of fluid forces becomes less important, relative to frictional forces. When
$\varphi \rightarrow 0$, the drag vanishes.
Following \citet{mazzei07}, we set
\begin{align}
&D ~=~ \frac{3}{4} C_{\mathrm{d}} \frac{\rho^{\mathrm{f}} \varphi |\vec{v} - \vec{u}| }{d} (1-\varphi)^{-\beta}
 \label{dragforce}
\end{align}
where $d$ is the mean particle diameter, $\beta$ is a~constant related to the constant $n$ in Ricardson--Zaki equation \citep{khan89}. According to  \citet{mazzei07}, this constant equals $2.80$ either when $R_{e_{p}} \to 0$ or $R_{e_{p}} \to \infty$, thus we use $\beta=2.80$ in Eq.~(\ref{dragforce}).  Finally, he drag coefficient is assumed to be constant $C_{\mathrm{d}} = 1$.

\subsection{Scaling}

The characteristic
thickness of the flowing granular material is $H$ and the characteristic length $L$. We scale
$x$ and $y$ by $L$, and $z$ by $H$, time by the free fall time $\sqrt{L/g}$, and
the $x, \, y$ and $z$~velocities by $\sqrt{Lg}$ and $\frac{H}{L}\sqrt{Lg}$, respectively. The
stresses are  scaled by $\rho^{\mathrm{s}} g H$ for the solids  phase and $\rho^{\mathrm{f}} g H$ for the fluid phase. After scaling the mass balance equations are unchanged.
Several terms in the momentum equations contain the factor $\epsilon = H/L$ which is small; $\epsilon \approx0.01 -  0.001$ is
not uncommon \citep{iverson01}.  Writing $x,\,y,\,z$ for $x_1,\,x_2,\,x_3$,
 the solid momentum balance equations become
\begin{align}
\label{scaledsolidmom}
&\partial_t ( \varphi v_{x}) + \partial_x ( \varphi v_x v_x)
+ \partial_y (  \varphi v_y v_x)
+ \partial_z ( \varphi v_z v_x)  \\ \nonumber
&  = \partial_x    \epsilon T_{xx}^{\mathrm{s}}   +  \partial_y   \epsilon T_{xy}^{\mathrm{s}}
+  \partial_z    T_{xz}^{\mathrm{s}}
  - \epsilon \frac{\rho^{\mathrm{f}}}{\rho^{\mathrm{s}}} \varphi \partial_x T_{xx}^{\mathrm{f}} \\ \nonumber   &
  - \epsilon \frac{\rho^{\mathrm{f}}}{\rho^{\mathrm{s}}} \varphi \partial_y T_{xy}^{\mathrm{f}}
 - \frac{\rho^{\mathrm{f}}}{\rho^{\mathrm{s}}} \varphi \partial_z T_{xz}^{\mathrm{f}}
   + \frac{D}{\rho^{\mathrm{s}}} (v_x - u_x)
 +  \varphi  g_x   \\
  & \partial_t ( \varphi v_{y}) + \partial_x( \varphi v_x v_y)
 +\partial_y( \varphi v_y v_y)
 +\partial_z (  \varphi v_z v_y) \\  \nonumber
& = \partial_x    \epsilon T_{xy}^{\mathrm{s}}   +  \partial_y    \epsilon T_{yy}^{\mathrm{s}}
+  \partial_z    T_{yz}^{\mathrm{s}}
  - \epsilon \frac{\rho^{\mathrm{f}}}{\rho^{\mathrm{s}}} \varphi \partial_x T_{xy}^{\mathrm{f}}  \\ \nonumber   &
   - \epsilon \frac{\rho^{\mathrm{f}}}{\rho^{\mathrm{s}}} \varphi \partial_y T_{yy}^{\mathrm{f}}
 - \frac{\rho^{\mathrm{f}}}{\rho^{\mathrm{s}}} \varphi \partial_z T_{yz}^{\mathrm{f}}
   +  \frac{D}{\rho^{\mathrm{s}}}  (v_y - u_y)
   +  \varphi  g_y   \\
  &     \epsilon (  \partial_t (\varphi v_z) + \partial_x (\varphi v_x v_z)
+ \partial_y(  \varphi v_y v_z)
+ \partial_z ( \varphi v_z v_z)  ) \\ \nonumber
&= \partial_x    \epsilon T_{xz}^{\mathrm{s}}   +  \partial_y   \epsilon T_{yz}^{\mathrm{s}}
+  \partial_z     T_{zz}^{\mathrm{s}}
- \epsilon \frac{\rho^{\mathrm{f}}}{\rho^{\mathrm{s}}} \varphi \partial_x T_{xz}^{\mathrm{f}}   \\ \nonumber
 &- \epsilon \frac{\rho^{\mathrm{f}}}{\rho^{\mathrm{s}}} \varphi \partial_y T_{yz}^{\mathrm{f}}
 - \frac{\rho^{\mathrm{f}}}{\rho^{\mathrm{s}}} \varphi \partial_z T_{zz}^{\mathrm{f}}
   +  \epsilon \frac{D}{\rho^{\mathrm{s}}}  (v_z - u_z) +  \varphi  g_z
\end{align}
Note that
components of gravity have been scaled by the magnitude $|g|$, so $(g_x,\, g_y,\,g_z)$ is  a~unit vector.


With the same scaling, the fluid momentum balance equations become
\begin{align}
\label{scaledfluidmom}
 & \partial_t ( (1-\varphi) u_{x}) + \partial_x (  (1- \varphi) u_x u_x)
 + \partial_y (   (1-\varphi) u_y u_x)
+  \partial_z (  (1-\varphi) u_z u_x) \nonumber \\
& =   \partial_x   \epsilon (1-\varphi)  T_{xx}^{\mathrm{f}}   +
 \partial_y   \epsilon (1-\varphi)   T_{xy}^{\mathrm{f}}
+  \partial_z    (1-\varphi)  T_{xz}^{\mathrm{f}}
- \frac{D}{\rho^{\mathrm{f}}}   (v_x - u_x)    + (1-\varphi)  g_x  \\
  &\partial_t ( (1-\varphi) u_{y})
   + \partial_x (  (1- \varphi) u_x u_y)
+ \partial_y (   (1-\varphi) u_y u_y)
+  \partial_z (  (1-\varphi) u_z u_y) \nonumber\\
& =   \partial_x   \epsilon (1-\varphi)   T_{xy}^{\mathrm{f}}   +  \partial_y
\epsilon (1-\varphi)    T_{yy}^{\mathrm{f}}
+  \partial_z  (1-\varphi)   T_{yz}^{\mathrm{f}}
- \frac{D}{\rho^{\mathrm{f}}}   (v_y - u_y)  + (1-\varphi)  g_y   \\
 & \epsilon (\partial_t ( (1-\varphi) u_z) + \partial_x(  (1-\varphi)u_x  u_z)
 \nonumber
+\partial_y ( (1-\varphi) u_y u_z) +\partial_z ( (1-\varphi) u_z u_z)  ) \\
& = \partial_x   \epsilon  (1-\varphi)    T_{xz}^{\mathrm{f}}   +  \partial_y
 \epsilon   (1-\varphi)    T_{yz}^{\mathrm{f}}
+  \partial_z   (1-\varphi)  T_{zz}^{\mathrm{f}}
- \epsilon \frac{D}{\rho^{\mathrm{f}}}  (v_z - u_z)  +  (1-\varphi)  g_z
\end{align}


In summary, then, the proposed equation system consists of the
solid volume fraction $\varphi$, the three solid velocities $\vec{v}$,  and  three fluid velocities $\vec{u}$.
These variables evolve according to the six momentum balance laws for the species, and the mass conservation relations
for each species.


\subsection{Constitutive assumptions and boundary conditions}

The upper surface of the flowing mass at $F_h ({x,y}, t) = 0$
is assumed to be a~material surface and stress free. At the base of the mass,
material is assumed to flow tangent to the basal surface $F_{\mathrm{b}} = 0$, and to satisfy a~sliding
friction law. For the solid constituent, this friction relation specifies that the shear traction and the normal stress
are proportional: $S|_{F_{\mathrm{b}}} = -{\text{sgn}}({v}) \, N|_{F_{\mathrm{b}}} \tan(\phi_{\text{bed}})$,
where $\phi_{\text{bed}}$ is the basal friction angle
and the $-{\text{sgn}}({v}) $ specifies that the shear traction
opposes motion.


We now discuss constitutive relations.
A~Coulomb constitutive relation is postulated for the material.  The Coulomb law
is a~nonlinear relation among the components of ${\mathbf{T}}^{\mathrm{s}}$, and stipulated that
material deforms when the total stress reaches yield, $||\text{dev}({\mathbf{T}}^{\mathrm{s}})|| = \kappa \, \text{tr}({\mathbf{T}}^{\mathrm{s}})$,
where $\text{dev}({\mathbf{T}}^{\mathrm{s}}) = {\mathbf{T}}^{\mathrm{s}} - {\frac{1}{2}} \text{tr}({\mathbf{T}}^{\mathrm{s}}) \mathbf{I}$ is the stress deviator,
$\text{tr}({\mathbf{T}}^{\mathrm{s}})$ is the trace of the stress (the sum of the diagonal components), $\mathbf{I}$ is the
identity tensor, and $\kappa$ is a~material parameter, and that as
deformation occurs, the stress and strain-rate tensors are aligned. That is, the alignment condition specifies
$\text{dev}({\mathbf{T}}^{\mathrm{s}}) = \lambda \, \text{dev}(\mathbf{V})$, where the strain-rate $\mathbf{V}= - {\frac{1}{2}} (\nabla {\vec{v}} +
\nabla {\vec{v}}^{\dagger})$ and $\dagger$ denotes the transpose.
To avoid a~switching between plastic and non-plastic behavior, we assume
 the solid material is everywhere in plastic yield.

The full Coulomb relations are too complex to be used here.
Two simplifications are proposed. First, at the basal surface the boundary condition ensures  proportionality and alignment of the tangential and normal forces; we
assume the same proportionality and alignment holds  throughout
the thin flowing layer of material. Written in components $ij$, this implies ${T_{\text{ij}}}^{\mathrm{s}} = \nu _{\text{ij}}{T_{zz}}^{\mathrm{s}}$,
where the proportionality constant $\nu$ is a~function of $\phi_{\text{bed}}$.
Second, following  \citet{rankine57}
and  \citet{terzaghi36},
an earth pressure relation is assumed for the diagonal
stress components,  $\nu _{\text{ii}}=k_{\text{ap}}$, or ${T_{xx}}^{\mathrm{s}} = k_{\text{ap}} {T_{zz}}^{\mathrm{s}}$, where
\begin{align*}
k_{\text{ap}} =  2 \frac{1 \pm [1 - \cos^2 (\phi_{\text{int}}) [1 + \tan^2 (\phi_{\text{bed}})]]^{{\frac{1}{2}}} }
  {\cos^2 (\phi_{\text{int}})}      -1 .
\end{align*}
\hack{\noindent}%
Here $\phi_{\text{int}}$ is the internal friction angle  and
the choice of the plus or minus sign depends on whether
flow is locally expanding (the active state, with $\nabla \cdot \vec{v} > 0$, and the $-$ sign) or contracting
(the passive state, with  $\nabla \cdot \vec{v} < 0$, and the $+$ sign).


For the fluid, the stress terms in \ref{scaledfluidmom} should be such that in case of $\varphi_s \to 0$ the equations agrees with the pure fluid depth averaged equations.

For pure water, the shear stress at wall can be approached by \citet{chow69}
\begin{align*}
\tau_w =  \rho^f g R_h S_f    -1 .
\end{align*}
\hack{\noindent}%

where $S_f$ is the slope friction and $R_h$ is the hydraulic ratio. Note that for shallow water problems $R_h \to h$. There are several approaches for the slope friction, thus for the shear stress at wall as well. For example, \citet{cun-hong06} use the empirical Manning approach, whereas \citet{liu78} and \citet{zhou99} use  the Chezy equation. Both Manning and Chezy appropaches pose numerical problmes when $h \to 0$. Thus we use use the Darcy equation \citet{guo95, xu06}: 

\begin{align*}
%\label{eq:wshear}
\tau_{xz} = -\rho^{\mathrm{f}} C_f u^2_x
\end{align*}
\hack{\noindent} 

where

\begin{align*}
\label{eq:cf}
C_f=\rho^f \frac{f_r}{8}
\end{align*}
\hack{\noindent} 
 
is the Chezy coefficient. The friction coefficient $f_r$   depends on the Reynolds number and the roughness of the channel ($k_s$). 

\section{Depth averaging}

The final step in the derivation is a~depth averaging of the mass and momentum balance equations. In this and the following sections we will use the same notation used by \citet{savage89}. If $h(x,y,t)$ is the unsteady surface of the flow and $b(x,y)$ is the terrain surface, for some function $f$, we compute
\begin{align*}
\frac{1}{|h-b|} \int_{b(x,y)}^{h(x,y,t)} \, f(x,y,z,t) \, {{\mathrm{d}}}z  = \overline f .
\end{align*}
where $h-b$ is the flow depth at a point $(x,y)$ and time $t$. 
Repeated use of the Leibnitz rule is made to
interchange integration and differentiation, and boundary conditions are employed to evaluate
terms at $b$ and $h$. In addition, several approximations must be made during the depth averaging process.
In what follows, we only briefly sketch the depth averaging process, noting as appropriate
those places where approximations are made. \citet{pitman05} provides an estimation of the errors typically made by these assumptions.

The terms of order $\epsilon$ are assumed small and we hope to drop all such terms from the model. However
 \citet{savage89} argues that
  diagonal contributions to the solid stress must be retained.    Because there is
no preferential downslope direction in our applications, and the flow direction
may change during a~flow, we retain the stress terms in both
the $x$ and $y$~directions,
dropping only
$O(\epsilon)$ terms in the $z$~direction; see the discussion in \citet{iverson01}.
Other $O(\epsilon)$ terms are dropped.


%\hack{\noindent}\textbf{Mass balance equation}
\subsection*{Mass balance equation}
As $\rho^s$ and $\rho^f$ are constants, equations \ref{eq:mass-s-1} can be reduced to
\begin{align*}
 \nabla \cdot ( \varphi {\vec{v}} + (1-\varphi) {\vec{u}} ) =0 ~.
\end{align*}
which says that the volume-weighted mixture flow is divergence free.
This equation is integrated from $z=b$ to $z=h$:
\begin{align}
&\int_b^h \nabla \cdot ( \varphi {\vec{v}} + (1-\varphi) {\vec{u}} ) \, \mathrm{d}z =0  .
\end{align}
\hack{\noindent}%
Use the Leibniz rule to interchange differentiation
and integration.  The upper free surface $F_h  = 0$ is
a~material surface for the mixture, so, following \citet{savage89} and \citet{pitman05},  at $z=h(x,y,t)$
\begin{align}
&\label{freesurf}
\partial_t ( \varphi h +  (1-\varphi) h) +
( \varphi v_{x} + (1-\varphi) u_x) \partial_x h \\ \nonumber
&+ ( \varphi v_{y} + (1-\varphi) u_y) \partial_y h
-   ( \varphi v_{z} + (1-\varphi) u_z)~ =~0 ~.
\end{align}
\hack{\noindent}%

At the basal surface $F_{\mathrm{b}} = 0$, the flow is tangent to the fixed bed and the bed is fixed in time, thus we can drop in \ref{freesurf} the terms in $\partial_t$, taking into account that at the surface $z=b(x,y)$ \citep{pitman05}
\begin{align}
\label{basalsurf}
&( \varphi v_{x} +
 (1-\varphi) u_x) \partial_x b  &
+  ( \varphi v_{y} +
  (1-\varphi) u_y) \partial_y b  -  ( \varphi v_{z} +
  (1-\varphi) u_z)=~0~.
\end{align}
\hack{\noindent}%
In arriving at these equations, we have ignored sedimentation, entrainment, and infiltration
of fluid into the bed.

Using these formulae and algebraic manipulation,the depth averaged equation for the total mass of the solid and fluid can be written
\begin{align}
&\partial_t \hat h + \partial_x (\hat h( {\overline{\varphi  v_{x}}} +  {\overline{(1-\varphi ) u_{x}}}))
+ \partial_y (\hat h( {\overline{\varphi  v_{y}}} +  {\overline{(1-\varphi ) u_{y}}}))  =  0  .
\label{totalmass}
\end{align}
\hack{\noindent}%
In writing this equation, the
depth averaged velocities are
$\hat h \overline {v_x} =  \int_b^h v_x \,{{\mathrm{d}}}z$, with a~similar expression for the volume fraction
$\overline {\varphi}$ and the other velocity components, and as in \citet{savage89},  $\hat h = h-b$.



\hack{\noindent}%
%z$\textbf{-Momentum}
\subsection*{$z$ momentum}
Observe that, upon setting $\epsilon \to 0$  in the
fluid $z$~momentum equation, we find the fluid to be hydrostatic:
\begin{align*}
\partial_z T_{zz}^{\mathrm{f}} =  g_{z} .
\end{align*}
\hack{\noindent}%
Integrating and imposing boundary conditions, we find
\begin{align}
&T_{zz}^{\mathrm{f}} (x,y,z) ~=~ -g_z [h - z] ,
\end{align}
\hack{\noindent}%
and the average
\begin{align}
&\label{fluidzstress}
\overline T_{zz}^{\mathrm{f}} (x,y) ~=~ - {\frac{1}{2}} g_z \hat h .
\end{align}


In the same manner, for the solid $z$~momentum
we find equation for an effective stress:
\begin{align*}
\partial_z T_{zz}^{\mathrm{s}} + {\frac{\rho^{\mathrm{f}}}{\rho^{\mathrm{s}}}}
\varphi \partial_z T_{zz}^{\mathrm{f}} = \varphi g_{z} .
\end{align*}
\hack{\noindent}%
Substituting,
\begin{align*}
\partial_z T_{zz}^{\mathrm{s}} ~=~  \left(1-{\frac{\rho^{\mathrm{f}}}{\rho^{\mathrm{s}}}}\right) g_z \varphi 
\end{align*}
\hack{\noindent}%
Thus the normal solid stress in the $z$~direction at any height
is equal to the (buoyancy) reduced weight of the solid material overburden.

In scaling these equations, the $z$~velocities have been dropped.
 Of course neglecting motion in the $z$~direction is
a~central component of a~thin layer theory.
Furthermore, any contribution to the $z$~momentum from fluid shearing  --  terms
such a~$T_{xz}^{\mathrm{f}},\, T_{yz}^{\mathrm{f}}$  --
are dropped due to scaling. This observation suggests that
only pressure contributions to
the fluid stress are important, an assumption we will make below, albeit with a~modification at the basal surface.


\hack{\noindent}%
%$x$\textbf{- and }$y$\textbf{ momentum}
\subsection*{$x$ and $y$ momentum}
We now must depth average the remaining momentum
equations.
The nonlinearity of these equations present difficulties in formulating
a~depth averaged theory,
complicate the derivation, and in several places,
it is necessary to take the depth
average of products of terms.
When necessary we approximate the required closure relation, for example as
$\overline {\varphi f}  \approx \overline \varphi  \overline f$ \citep{savage89, pitman05}.

Consider first the equation for the motion of the solid phase.
The left-hand side of the  $x$~momentum equation (\ref{scaledsolidmom})
can be written
\begin{align*}
\text{LHS} = \partial_t \varphi v_{x} + \partial_x \varphi {v_{x}}  +
\partial_y \varphi v_{x}v_{y} + \partial_z \varphi v_{x}v_{z} ~.
\end{align*}
Depth average and use boundary conditions to find
\begin{align}
&\int_{b}^{h}\text{LHS} \,{{\mathrm{d}}}z ~=~
 \partial_t\int_{b}^{h}\varphi v_{x} \,{{\mathrm{d}}}z +
\partial_x \int_{b}^{h}\varphi v_{x}^2 \,{{\mathrm{d}}}z +
\partial_y \int_{b}^{h}\varphi v_{x}v_{y} \,{{\mathrm{d}}}z
\label{LHS}
\end{align}

Now depth average the right hand side of (\ref{scaledsolidmom}):
\begin{align}
&\int_{b}^{h}\text{RHS} \,{{\mathrm{d}}}z = \underbrace{-\int_{b}^{h}(\epsilon \partial_x T^{\mathrm{s}}_{xx} + \epsilon \partial_y T^{\mathrm{s}}_{xy}
+ \partial_z T^{\mathrm{s}}_{xz})\, {{\mathrm{d}}}z}_{\text{(i)}}  \\ \nonumber
&\underbrace{- {\frac{\rho^{\mathrm{f}}}{\rho^{\mathrm{s}}}} \int_{b}^{h}\epsilon \varphi \partial_x T^{\mathrm{f}}_{xx}\,{{\mathrm{d}}}z}_{\text{(ii)}} +\underbrace{  \frac{D}{\rho^{\mathrm{s}}}  (v_x - u_x)}_{\text{(iii)}}+ \int_{b}^{h} \varphi g_x \,{\mathrm{d}}z  ~.
\end{align}

In order to proceed, several assumptions are made:
\begin{itemize}
\item{This equation governs the motion of the solid phase, and we assume the upper free surface for the mixture is
a~free surface for both of the individual phases.}
\item{The drag term $ \frac{D}{\rho^{\mathrm{s}}}  (v_x - u_x) $
is highly nonlinear and a~correct depth average is all but impossible to calculate. We postulate that
a~experiments could fit an averaged phenomenological drag
of a~similar form. As seen in equation \ref{dragforce} we assume a drag coeficient $C_d = 1$. In addition, we use the typical mean particle for lahars $d=1$ $~mm$ \citep{schmid81} }
\item{The earth pressure relation for the solid phase is employed.
That is, the basal shear stresses are assumed to be
proportional to the normal stress:
\begin{align*}
T^{\mathrm{s}}_{iz}=
- \frac{v_{i}}{||\vec{v}||} \tan(\phi_{\text{bed}}) T^{\mathrm{s}}_{zz} \equiv \alpha_{*z} T^{\mathrm{s}}_{zz} ,
\end{align*}
\hack{\noindent}%
where $i$ can be either $x$ or $y$,
and the velocity ratio enforces that friction opposes motion in the designated direction
 \citep{savage89,patra05}
The  $\alpha$ notation will provide a~convenient shorthand that we use in other places.
Likewise the diagonal stresses are taken to be proportional to the normal
solid stress
\begin{align*}
T^{\mathrm{s}}_{ii}= k_{\text{ap}} T^{\mathrm{s}}_{zz} \equiv \alpha_{ii} T^{\mathrm{s}}_{zz} ~.
\end{align*}
\hack{\noindent}%
Finally, following \citet{iverson01},
$xy$ shear stresses are determined by a~Coulomb relation
\begin{align*}
T^{\mathrm{s}}_{xy} = -{\text{sgn}}(\partial_y v_x) \sin(\phi_{\text{int}}) k_{\text{ap}} T^{\mathrm{s}}_{zz} \equiv \alpha_{xy} T^{\mathrm{s}}_{zz} ~,
\end{align*}
\hack{\noindent}%
where the sgn function ensures that friction
opposes straining in the $(x,y)$ plane. }
\item{For the fluid phase, the basal shear stresses  are assumed to be proportional to the square of the depth averaged velocities \citep{guo95,xu06}:
\begin{align*}
T^{\mathrm{f}}_{iz}=\rho^{\mathrm{f}} C_f u_i||\vec{u}||
~,
\end{align*}
\hack{\noindent}%
where $C_f$ is the Chezy coefficient, which depends on the friction coefficient $f_r$ (see equation \ref{eq:cf}. A~physical approach for this coefficient is the Colebrook--White equation \citep{colebrook37}, which for rough channels can be approximated as:
\begin{align*}
\label{eq:colebrook}
  \frac{1}{\sqrt{f_{\mathrm{r}}}}=-2 \log_{10}\left(\frac{k_{\mathrm{s}}}{14.8 R_{\mathrm{h}}} \right)
~,
%\hack{\noindent}%
\end{align*}
where $k_{\mathrm{s}}$ is the roughness of the channel and $ R_{\mathrm{h}}$ is the hydraulic radii, which for shallow water problems can be approached as the depth of the flow ($R_{\mathrm{h}} \to h$). Equation \ref{eq:colebrook} is logarithmic, thus large uncertainties in $k_{s}$ result only in small variations in $C_f$ \citep{swaffield83}. \citet{trrl75} proposes values of $k_s$ for different materials and channel types. We use $k_s = 1$ $~mm$ for channels in volcanic environments. Therefore, here $C_r$ will depend  on the flow depth $\hat h$, $\rho^{\mathrm{f}}$ and the fluid velocity $u_i$. Note that this is a physical approach for $T^f_{iz}$, which does not depend on empiric.

For pure fluids, the diagonal stresses and shear stress are zero.


\end{itemize}

From Leibniz's rule and the stress computations above,
we find
\begin{align}
&\text{(ii)} =  -\epsilon{\frac{\rho^{\mathrm{f}}}{\rho^{\mathrm{s}}}}
\int_{b}^{h} \varphi \partial_x T^{\mathrm{f}}_{xx} \,{{\mathrm{d}}}z =   \nonumber\\  &
-  \epsilon {\frac{\rho^{\mathrm{f}}}{\rho^{\mathrm{s}}}} \overline \varphi
[\partial_x \left(\frac{{\hat h}^2}{2} (-g^z)\right)
+  \hat h (-g_z) \partial_x b \,] .
\end{align}


Now using the fluid and solid stress relation
$\overline T^{\mathrm{s}}_{zz} = \left(1-{\frac{\rho^{\mathrm{f}}}{\rho^{\mathrm{s}}}}\right)\overline \varphi \overline T^{\mathrm{f}}_{zz}$,
term (i) is approximated as:
\begin{align}
&\text{(i)} = -\epsilon \int_{b}^{h} \partial_x \alpha_{xx} T^{\mathrm{s}}_{zz} \, {{\mathrm{d}}}z
- \epsilon  \int_{b}^{h} \partial_y \alpha_{xy} T^{\mathrm{s}_{zz}} \, {{\mathrm{d}}}z
-   \int_{b}^{h} \partial_z \alpha_{xz} T^{\mathrm{s}}_{zz} \, {{\mathrm{d}}}z \\ \nonumber
&= -\epsilon [\partial_x \int_{b}^{h} \alpha_{xx} T^{\mathrm{s}}_{zz} \, {{\mathrm{d}}}z
 - \alpha_{xx} T^{\mathrm{s}}_{zz} |_{z=h} {\partial_x h}  \nonumber
+ \alpha_{xx} T^{\mathrm{s}}_{zz} |_{z=b} {\partial_x b} ]
-\epsilon [\partial_y \int_{b}^{h} \alpha_{xy} T^{\mathrm{s}}_{zz} \, {{\mathrm{d}}}z \\ \nonumber
& - \alpha_{xy} T^{\mathrm{s}}_{zz} |_{z=h} {\partial_y h}
     + \alpha_{xy} T^{\mathrm{s}}_{zz} |_{z=b} {\partial_y b} ]  \nonumber
- \alpha_{xz} [T^{\mathrm{s}zz} |_{z=h} - T^{\mathrm{s}zz} |_{z=b} ] ~.
\end{align}
\hack{\noindent}%
The upper free surface is stress free, so all terms involving  $T^{\mathrm{s}}_{zz} |_{z=h}$
vanish.  The expression for (i) becomes
\begin{align}
&\text{(i)} =  -\epsilon \left(1-{\frac{\rho^{\mathrm{f}}}{\rho^{\mathrm{s}}}}\right)
\partial_x\left( \hat h  \alpha_{xx} \overline \varphi \overline T^{\mathrm{f}}_{zz} \right)  \nonumber
 - \epsilon \left(1-{\frac{\rho^{\mathrm{f}}}{\rho^{\mathrm{s}}}}\right) \partial_y
 \left( \hat h  \alpha_{xy} \overline \varphi \overline T^{\mathrm{f}}_{zz} \right) \\
& +\left(1-{\frac{\rho^{\mathrm{f}}}{\rho^{\mathrm{s}}}}\right)(-\epsilon\alpha_{xx}\partial_x
b-\epsilon\alpha_{xy}\partial_y b +
\alpha_{xz}) (-g_z) \hat h  \overline \varphi ~.
\end{align}
\hack{\noindent}%
Note that the factor $(-g_z)$ originates in the evaluation and depth averaging
of the fluid stress; in typical flows, this factor is positive.

Combining all terms
yields a~solids $x$~momentum equation:
\begin{align}
\label{solidxmom}
 & \partial_t \left( \hat h  \overline \varphi \overline {v_{x}} \right) +
 \partial_x \left( \hat h  \overline{\varphi} \overline{v_{x}} \overline {v_{x}} \right)
+ \partial_y \left( \hat h   \overline \varphi \overline{v_{x}} \overline{v_{y}} \right)
= - \frac{\epsilon}{2} \left(1-{\frac{\rho^{\mathrm{f}}}{\rho^{\mathrm{s}}}}\right)
\partial_x ( \alpha_{xx} {\hat h}^2 \overline \varphi (-g_z) ) \\ \nonumber
&  - \frac{\epsilon}{2}\left(1-{\frac{\rho^{\mathrm{f}}}{\rho^{\mathrm{s}}}}\right)
\partial_y ( \alpha_{xy} {\hat h}^2 \overline \varphi (-g_z) )  \nonumber
 +\left(1-{\frac{\rho^{\mathrm{f}}}{\rho^{\mathrm{s}}}}\right)(-\epsilon\alpha_{xx}
 \partial_x b -\epsilon\alpha_{xy}\partial_y b
+\alpha_{xz})\hat h \overline \varphi  (-g_z)  \\ \nonumber
& -  \frac{\epsilon}{2} {\frac{\rho^{\mathrm{f}}}{\rho^{\mathrm{s}}}} \overline
\varphi \partial_x ( {\hat h}^2 (-g_z) ) \nonumber
 -\epsilon{\frac{\rho^{\mathrm{f}}}{\rho^{\mathrm{s}}}} \hat h \overline
 \varphi (-g_z) \partial_x b
+ \left(\frac{\overline{D}}{\rho^{\mathrm{s}}} \right) \left(\overline {u_{x}}-
\overline {v_{x}} \right)
+ \hat h \overline \varphi g^{x}.
\end{align}

The $y$~solid momentum equation can be derived in a~similar fashion.

The equation for the fluid motion presents fewer difficulties.
 The depth averaged
$x$~momentum equation takes the form
\begin{align}
 &  \partial_t \left( \hat h  \overline {\varphi^{\mathrm{f}}}  \overline {u_{x}} \right)+  \partial_x
 \left( \hat h \overline {u_{x}} \overline {\varphi^{\mathrm{f}}} \overline {u_{x}} \right)
+  \partial_y \left( \hat h  \overline {\varphi^{\mathrm{f}}} \overline{u_{x}} \overline{u_{y}} \right) \\
&= - {\frac{1}{2}} \epsilon  \partial_x {\hat h}^2 \overline {\varphi^{\mathrm{f}}} (-g_z)
- \left(\frac{\overline{D}}{\rho^{\mathrm{s}}}  \right) \left(\overline {u_{x}}- \overline {v_{x}} \right)
+ \hat h \overline{\varphi^{\mathrm{f}}} g_{x}
 + \overline{\varphi^{\mathrm{f}}} C_f u_x ||\vec{u}|| \nonumber .
\end{align}
\hack{\noindent}%
Where $\overline {\varphi^{\mathrm{f}}}=1- \overline \varphi$. Again, the fluid $y$~momentum equation has a~similar form. Note that
if $\overline {\varphi^{\mathrm{f}}} \to 1$ Eq.~(25) becomes the typical shallow water approach of hydraulics \citep{chow69} (\citealp{kowalski08} describes how debris flows reduces to  a~shallow water flow as solid volume fraction vanishes).

 As noted above, we solve for volume fractions ($\varphi$), thus the bulk density can be calculated from:
\begin{align}
&
  \rho =\rho^{\mathrm{s}} \varphi +\rho^{\mathrm{f}} (1-\varphi)\label{eq:rho}
\end{align}
Then, we obtain the dynamic pressure $p$ from:
\begin{align}
&\label{eq:dpress}
   p = \frac{1}{2} \rho {\overline{\vec{v}}}^2
\end{align}
where $\overline{\vec{v}} $ is the mixture averaged velocity of the flow \citep{fan98}
\begin{align}
&\label{eq:dpress_v}
   \overline{\vec{v}}=\frac{\rho^s \overline{v} \varphi+\rho^f \overline{u} (1-\varphi)}{\rho}
\end{align}
where  $\overline{v}=\sqrt{v^2_x+v^2_y}$ and $\overline{u}=\sqrt{u^2_x+u^2_y}$ are the  speeds of each phase.
  
The use of the impacting dynamic pressure information
on structures and living beings  allows to  estimate levels of damage, as in
\citet{valentine98}, \citet{valentine11}, and \citet{jones12}, useful for vulnerability analysis.

\section*{Numerical solution}
This system of equations is then solved using the finite volume method, whose  solution  provides results of the velocity field, flow depth and the volumetric concentration of solids at the center of the each finite volume computational grid.

 To solve the  balance laws, we use the parallel, adaptive mesh, Godunov solver developed by \citet{davis88} already implemented in \citet{patra05} and \citet{pitman05}. The adaptive meshing is used as well, which allows to have very fine grids where indicators show high gradients, and coarser grids where low gradients are detected. The time step is adjusted from a CFL condition \citep{courant28}. The complexity of the equation system results in typical time steps of the order of $10^{-4}$.  However, the maximum time step allowed was reduced from $10^{-1}$\,s used within Titan2D to $10^{-3}$\,s to ensure stability.  As consequence of this small time step,  Titan2F become a~computationally expensive tool.


The numerical solution of the above set of equations presents strong numerical sensitivity to small changes and Digital Elevation Model (DEM) errors and the quality of that maps. The DEMS can have regions where elevations are not well defined, they can have crossed contour levels or even infinite holes. Abrupt terrain changes, both actual or DEM artifacts, cause computations of gradients and curvatures to become unstable. In order to avoid such a~numerical problems patching and intelligent smoothing of the DEMs was needed.

From the hyperbolicity analysis done by \citet{pitman05}, we try to ensure hyperbolicity imposing a minimum $\hat h_{min} =10^{-5}$ $~m$, a maximum $\varphi_{max}=\varphi_{pc}$ correspondent to a maximum packing concentration of $0.65$ and a minimum $\varphi_{min}=10^{-8}$. This constrains makes the program very stable if a good DEM is used.

The needed initial conditions are the location of the pile, its geometry, the volumetric solid concentration and initial velocity. There is no inflow condition implemented yet. The bed and internal friction are set internally as fixed values. Nevertheless, both of this parameters can affect  $k_{ap}$. The red line in Figure \ref{fig:kap_bed_phi} shows how $K_{ap}$ changes with $\phi_{bed}$ for {\em any} value of $\phi_{int}$ (Note that \citet{williams08} shows that Titan2D results are not strongly affected by  $\phi_{int}$). The black line in the same figure shows how $k_{ap}$ evolves when the fix value of $\phi_{bed}=40^o$, used in the program, is multiplied by $0<\varphi<0.65$. We set internally  $\phi_{int}=42^o$ as well. Then, we use in the evaluation of $k_{ap}$, the result of the multiplication $\varphi \times \phi_{bed} $ and  $\varphi \times \phi_{int} $ instead of $\phi_{bed}$ and  $\phi_{int}$ respectively.

\section*{Differences with \citet{pitman05} model}

There are five major differences between the present paper and
\citet{pitman05}:
\begin{enumerate}
\item{In \citet{pitman05}, mass and momentum conservation laws are derived
for the solid material and for the phase averaged mixture of solid and fluid, whereas here
the final model presents mass and momentum equations for both individual phases. }
\item{ Any two phase model system must postulate several phenomenological
functions, such as inter-phase drag coefficient.
In the present derivation these functions are better adapted to geophysical flows whose fluid phase corresponds to water and the solid phase are rounded solid particle. The drag is calculated from an expression valid for the whole range of Reynolds numbers. It only needs the mean particle diameter of the flow as a parameter. We use a typical mean particle diameter of lahars. }
\item{The volumetric particle concentration is no longer a~fixed parameter, which in our approach is calculated for every time step and grid point. This mean that instead of having constant $\varphi$ like in equations 3.2 in \citet{pitman05} or simply not accounted for, like in the depth averaged fluid mass and momentum balance equations 3.27 and 3.28 in \citet{pitman05}, we include the variable $\varphi$ within all derivatives and depth averaged equations. Thus, in regions when the particle concentration vanishes, the solid phase role in the equation system vanishes as well. In that way,  the equation system becomes the typical hydraulic shallow water approach, which does not happen in \citet{pitman05}.}
\item{We account for the fluid stresses at wall. We use a physical approach that only needs the roughness of the channel and the flow depth. The last term is calculated by the program at every time step, whereas that the former is set as fixed value.}
\item{The only input parameters needed by the program are the location of the pile of material, its volume and the volumetric solids concentration. The friction coefficients are no longer needed as input as they are automatically adjusted according to the evolution of the volumetric fraction of solids across the grid and time. The bed and internal  friction are set in such way that when the volumetric fraction of solids tends to an assumed maximum packing concentration   ($\varphi=0.65$), both internal and bed frictions tends to the values used in that cases in Titan2D (see \citet{sheridan05} and \citet{williams08} for ranges of values).}
\end{enumerate}

\section{Validation and verification}

In order to validate and verify the proposed approach, we did a~series of tests using a~one dimensional approach of the proposed system of equations. First, we test the consistency of the results verifying the expected symmetries. Then we used an exponentially decaying topography with and without obstacles (bumps). They show expected features like how the mathematical model can show a~reflected wave when the bump was big enough or the flow splitting in a~reflected dense part and bump overpass by a~very dilute flow. In addition we verify the predictions with two dimensional experimental measurements done by \citet{liu96} and \citet{iverson10}. 

Verification of the accuracy of the code was done with analytical solutions of the Dam Break problem and with several experimental results. Among them, we check the deposited pattern predicted by the program with the results shown by   \citet{liu96}. The prediction of the arrival time and the flow depth profile was compared with the experimental results shown by  \citet{iverson10} from his recent work done on his large channel facility.

Analytical solutions for shallow water problems are scarce. Only one dimensional analytical approaches are available in the literature, especially for the well known Dam Break problem (e.g. \citealp{dressler54,mangeney00,fernandez06,ancey08,wrachien09}).  Unfortunately, analytical solutions for geo-mass flows are almost impossible to find due to  the complexity of the  non-linear partial differential equations \citep{pudasain07}. Such solutions can be obtained only in special cases like the similarity based solutions proposed by Savage and Hutter \citep{savage89} for dry avalanches. In our test we use the solution proposed by  \citet{fernandez06} for the Dam Break problem on an incline for pure water. In our program we assume  $\varphi^{\mathrm{s}} \to 0$.  Figure~\ref{fig:dambreak} shows a~comparison of the Titan2F prediction with this analytical solution. The statistical {\em t}-test shows that there is no statistical difference between the analytical solutions and the 1D prediction. That test shows that about $70 \%$ of the predictions of Titan2F can explain the analytical solution. At least for the one dimensional case, the program successfully reproduces analytical solutions for different initial conditions down to  very low particle concentrations (less than $1\,\unit{\%}$). However, it can be noted in Figure~\ref{fig:dambreak} that the 1D version of the program tends to over-estimate the front advance of the flow, at least in the initial parts of its evolution. 

\citet{liu96} performed several experiments for geo-mass flows in an inclined  channel. He modified the initial volume, the channel slope and the particle concentration to find the final size of the debris flows measured by their resulting width and length. We reproduced experimental final width and length after the simulation reached the same time corresponding to the duration reported by  Liu. Figure~\ref{fig:liu_exp} shows the correspondence of the model with the experiments for (a) the width of the deposit and (b) the length. A~Pearson correlation shows that $90\,\unit{\%}$ of the experimental data for the deposit length can be explained with the predictions of Titan-Two-Phase, whereas  $80\,\unit{\%}$ of the data for the deposit width can be explained by the predictions of the program. This illustrates the high accuracy of the program in predicting the deposit characteristics for different initial volumes and high initial solid concentrations.

The experiments performed by Iverson and co-workers \citep{iverson10} done on a~$95\,\unit{m}$ long artificial channel was used to verify the accuracy of the predictions  of the flow front arrival time and the temporal evolution of the flow depth. These flows were unsteady and nonuniform. \citet{iverson10} show time series data for several measured properties:  flow thickness, pore pressure, basal normal stress and arrival time of the front. Raw data sent to us by Dr. Iverson were used to test the Two-Phase-Titan prediction concerning time evolution of the flow depth and arrival times at the check points located at $32$ and $66\,\unit{m}$ distance from the lock. As shown in Figures~\ref{fig:iverson1} and \ref{fig:iverson2}, the arrival time and the temporal evolution of the predicted depth fits very well within the confidence interval of the experiments. In both of the cases, Titan2F tends to over-estimate the flow depth just after the arrival of the front. This is probably an effect of the slight difference in the shape of the initial pile, as the free surface of the numerical pile follows the same slope of the channel, whereas the actual free surface within the lock is horizontal. 

The range of concentrations  that the program cope with, are from $\varphi_s=10^{-8}$ (almost pure water) to $\varphi_s =0.65$ (maximum packing concentration). Finally,  as expected, the program predicts high particle concentrations at the front of the flow and low particle concentrations at the tail of the flow (at some cases, even  near pure water concentrations, or $\varphi \to 0$), as can be seen in Fig.~\ref{fig:longitudinal_concentration} where a~longitudinal solids particle distribution predicted by Titan2F is shown. The predictions fits with \citet{iverson10} observation that the tail of the flow remains very watery.  Using predicted concentration of solids, the density is assessed (Eq.~\ref{eq:rho}) and together with the speed of the flow, the dynamic pressure distribution is calculated as well (Eq.~\ref{eq:dpress}). Figure~\ref{fig:dpress} shows  longitudinal and cross distributions of the dynamic pressures after $32\,\unit{s}$ simulation. As the flow velocity at lateral limits of the flow from the end of the channel are very low, the dynamic pressure shows to be low as well. Knowledge of the dynamic pressure information is of vital importance in risk analysis as structural damage and risk for human life can be assessed from it.

Verification with actual mud flows has been done as well, showing very good fit with field data. For example, \citet{sheridan11} shows that the Titan2F predictions are within $10\,\unit{\%}$  of the data shown by \citet{procter10} for the highly channeled mud flow at Ruapehu, New Zealand. In addition, the theory was tested against  field data assessed by \citet{williams08} for the 2006 Vazcun creek lahar at Tungurahua volcano, Ecuador, as shown by \citet{cordoba15}, where Titan2F predictions about velocity are within $10\,\unit{\%}$ and within $15\,\unit{\%}$ for the measured super-elevation. 


\conclusions

In this paper we present a~new computational two phase tool for lahar hazard assessment that has no constrains on initial volumetric particle concentration. The program computes space--time evolution of the particle concentration, flow depth, velocity field and dynamic pressure at each point of the computational grid.

The model is valid for two phase flows whose phases consist in solids and water. However, the phenomenological approach used for the interphase drag model assumes an average diameter of the solids, which mean individual boulders or particles cannot be tracked. In addition, the model is depth averaged, assuming thin layer and shadow water approaches.  Thus, our model correctly predicts the dynamic of gravity driven flows providing the depth averaged values for the particle concentration, flow and phases velocities and flow depth in a~three dimensional topography. In order to model other kind of geophysical mass flows, adjustments to the code must be done, for example pyroclastic flows can be modeled if the flow density of the fluid phase is appropriately addressed (eg. thermocoupled air density using ideal gas law and an additional equation for temperature).

The proposed mathematical approach allows to study the whole range of flow behavior. Regions with  almost pure fluid to regions of  friction dominated flows are correctly described by the algorithm. Using this information, dynamic pressure is deduced, which becomes a~very useful tool for risk assessment.

The highly non-linear coupled equation system makes the time step very small. The use of this new tool on natural terrains  or detailed Digital Elevation Models (DEMs)  requires higher computational power than the one provided for home PCs. The use of a~high velocity Work Station with multiple cores  is advised.

Important processes that are not addressed by this tool include the effect of turbulence, incorporation of solid material from the bed of the channel, and incorporation of water into the flow from existing water bodies. Nevertheless, this two-phase flow is an important step forward in forming an acceptable computational model for simulating a~hazardous natural phenomena. Currently, we are applying this tool for  real lahar hazard assessments like in \citet{cordoba15}.



\Supplementary{zip}

\authorcontribution{%% appendix
Gustavo C\'ordoba, Michael Sheridan and Bruce Pitman developed the mathematical approach of the model. Gustavo C\'ordoba did the coding.}


\begin{acknowledgements}
The authors wish to thank the economic support provided by the National Science  Foundation Grant  EAR 711497 aimed to the developing of Two-Phase-Titan.  We also want to thanks  Ramona Stefanescu,  Abani Patra,  Dinesh Kumar and  Keith Dalbey from University of Buffalo for the useful discussions about this development, and to Dr. Richard Iverson from USGS who sent us the raw data of his experiments.  The Universidad de  Nari\~no, Colombia, allowed the mobility of  Cordoba to assume the post-doctoral position at the University of Buffalo (2008--2010) where the mathematical approach and coding of the beta version was developed,  and later gave him the necessary dedication  time for this research. Finally, we thanks for the suggestions of the anonymous referees. who help us in the substantial improvement of the paper.
\end{acknowledgements}




\begin{thebibliography}{65}%


\bibitem[Ancey(2006)]{ancey06}
Ancey,~C.: Plasticity and geophysical flows: a~review, J. Non-Newton Fluid, 142, 4--35, 2006.


\bibitem[Ancey et~al.(2008)]{ancey08}
 Ancey,~C., Iverson,~R., Rentschler,~M., and Denlinger,~R.: An exact solution for ideal dam-break floods on steep slopes, Water Resour. Res., 44, W01430,
doi:\href{http://dx.doi.org/10.1029/2007WR006353}{10.1029/2007WR006353}, 2008.


\bibitem[Anderson and Jackson(1967)]{anderson67}
 Anderson,~T.~B. and Jackson,~R.: A~fluid mechanical description of fluidized beds: equations of motion, Ind. Eng. Chem. Fund., 6, 527--539, 1967.


\bibitem[Balmforth and Craster(1999)]{balmforth99}
 Balmforth,~N.~J. and Craster,~R.~V.: A~consistent thin-layer theory for Bingham plastics, J. Non-Newton Fluid, 84, 65--81, 1999.


\bibitem[Bedford and Drumheller(1983)]{bedford83}
 Bedford,~A. and Drumheller,~D.~S.: Recent advances: theories of immiscible and structured mixtures, Int.~J. Eng. Sci, 21, 863--960, 1983.


\bibitem[Xu(2004)]{xu06}
 Xu, K. : The gas-kinetic scheme for shallow water equations, Journal of Hydrodynamics, Proceedings of the Conference of
Global Chinese Scholars on Hydrodynamics,  73--76, Shangai, China, 11-14 July, 2006.


\bibitem[Chow(1969)]{chow69}
 Chow,~V.: Open Channel Hydraulics, McGraw Hill, New York, 1969


\bibitem[Colebrook and White(1937)]{colebrook37}
 Colebrook,~C. and White,~C.: Experiments with fluid friction in Roughened pipes, P.~R. Soc. London, 161, 367--381, 1937.


\bibitem[U.S.-Geodynamics-Committee(1994)]{nrc94}
 U.S.-Geodynamics-Committee: Mount Rainier: Active Cascade Volcano, National Academic Press, Washington, D.C.,  1994.


\bibitem[Cordoba et~al.(2015)]{cordoba15}
 C\'{o}rdoba,~G., Villarosa,~G., Sheridan,~M.~F., Viramonte,~J.~G., Beigt,~D., and Salmuni,~G.: Secondary lahar hazard assessment for Villa la Angostura, Argentina, using Two-Phase-Titan modelling code during 2011 Cord\'{o}n Caulle eruption, Nat. Hazards Earth Syst. Sci., 15, 757--766,
doi:\href{http://dx.doi.org/10.5194/nhess-15-757-2015}{10.5194/nhess-15-757-2015}, 2015.


\bibitem[Courant et~al.(1928)]{courant28}
 Courant,~R., Friedrichs,~K., and Lewy,~H.: On the partial difference equations of mathematical physics, Math. Ann., 100, 32--74, 1928.


\bibitem[Coussot(1997)]{coussot97}
 Coussot,~P.: Mudflow Rheology and Dynamics, IAHR Monograph Series, Balkema, 1997.


\bibitem[Dartevelle(2004)]{dartevelle04}
 Dartevelle,~S.: Numerical modeling of geophysical granular flows: 1. A~comprehensive approach to granular rheologies and geophysical multiphase flows, G-cubed-Geochemistry-Geophysics-Geosystems, 5, Q08003,
doi:\href{http://dx.doi.org/10.1029/2003GC000636}{10.1029/2003GC000636}, 2004.


\bibitem[Davis(1988)]{davis88}
 Davis,~S.~F.: Simplified second order Godunov-type methods, SIAM~J. Sci. Stat. Comp., 9, 445-473, 1988.


\bibitem[Dobran(1991)]{dobran91}
 Dobran,~F.: Theory of structured multiphase mixtures, Lect. Notes Phys., Volume 372, Springer, Berlin, Deutschland, 1991.


\bibitem[Dobran(2001)]{dobran01}
 Dobran,~F.: Volcanic Processes: Mechanisms in Material Transport, Kluwer Academic, NY,  2001.


\bibitem[Dobran et~al.(1993)]{dobran93b}
 Dobran,~F., Neri,~A., and Macedonio,~G.: Numerical Simulation of collapsing volcanic columns, J. Geophys. Res.,  98, 4231--4259, 1993.


\bibitem[Dressler(1954)]{dressler54}
 Dressler,~R.~F.: Comparison of theories and experiments for hydraulic dam
break waves, Proc. Int. Assoc. of Scientific Hydrology Assemblee Generale, 38, 319--328, 1954.


\bibitem[Ergun(1952)]{ergun52}
 Ergun,~S.: Fluid flow through packed columns, Chem. Eng. Prog., 49, 89-94, 1952.


\bibitem[Fan and Zhu(1998)]{fan98}
 Fan,~L. and Zhu,~C.: Principles of Gas-Solid Flows, Cambridge CB2 2RU, Cambridge, 1998.


\bibitem[Fernandez-Feira(2006)]{fernandez06}
 Fernandez-Feira,~R.: Dam-break flow for arbitrary slopes on the bottom, J. Eng. Math., 54, 319--331, 2006.


\bibitem[Gidaspow(1994)]{gidaspow94}
 Gidaspow,~D.: Multiphase Flow and Fluidization, Academic Press, New York, 1994.


\bibitem[Hutter et~al.(2003)]{hutter03}
 Hutter,~K., Wang,~Y., and Pudasaini,~S.~P.: The Savage-Hutter Avalanche Model: How far can it be pushed?, Preprint, NI03072-GPF, Isaac Newton Institute for Mathematical Sciences, Cambridge, {UK},  2003.


\bibitem[Gray et~al.(1999)]{gray99}
 Gray,~J.~M., Wieland,~M., and Hutter, K.: Gravity-driven free surface flow of granular avalanches over complex basal terrain, P. Roy. Soc. Lond. A, 455, 1841--1874, 1999.

\bibitem[Guo(1995)]{guo95}
 Guo,~Z.: Velocity-depth coupling in shallow-water flows,
J. Hydrol. Eng., 121, 717--724, 1995.

\bibitem[Iverson(1997)]{iverson97}
 Iverson,~R.~M.: The physics of debris flows, Rev. Geophys.,  35, 245--296, 1997.

\bibitem[Iverson and Denlinger(2001)]{iverson01}
 Iverson,~R.~M. and Denlinger,~R.~P.: Flow of variably fluidized granular material across three-dimensional terrain 1: Coulomb mixture theory,~J. Geophys. Res., 6, 537--552, 2001.

\bibitem[Iverson et~al.(2010)]{iverson10}
 Iverson,~R.~M., Logan,~M., La{H}usen,~R.~G., and Berti, M.: The perfect debris flow? Aggregated results from 28 large-scale experiments,~J. Geophys. Res., 115, F03005,
doi:\href{http://dx.doi.org/10,1029/2009FJ001514}{10,1029/2009FJ001514}, 2010.

\bibitem[Jackson(2000)]{jackson00}
 Jackson,~R.: The Dynamics of Fluidized Particles, Cambridge Monographs on Mechanics, Cambridge University Press, Cambridge, {UK}, 2000.


\bibitem[Jones(2012)]{jones12}
 Jones,~N.: Damage of plates due to impact, dynamic pressure and explosive loads, Lat. Am. J. Solids Stru., 10, 767--780, 2012.


\bibitem[Khan and Richardson(1989)]{khan89}
 Khan,~A.~R. and Richardson,~J.~F.: Fluid-particle interactions and flow characteristics of fluidized and settling beds ans settling suspensions of spherical particles, Chem. Eng. Commun., 78, 111--130, 1989.


\bibitem[Kowalski(2008)]{kowalski08}
 Kowalski,~J.: Two-{P}hase Modeling of Debris Flows, PhD thesis, Eidgenossische Technische Hochschule ({ETH}), Zurich, Switzerland, 135 pp., 2008.

\bibitem[Le(2006)]{le06}
Le, L.: New Models for geophysical mass flows. PhD Thesis. Department of Mathematics. State University of New York at Buffalo, 2006.


\bibitem[Legros(2002)]{longrunout02}
 Legros,~F.: The mobility of long runout landslides, Eng. Geol., 63, 301--331, 2002.


\bibitem[Liu(1996)]{liu96}
 Liu,~X.: Size of a~debris flow deposition: model experiment approach, Environ Geol., 28, 70--77, 1996.


\bibitem[Mangeney et~al.(2000)]{mangeney00}
 Mangeney,~A., Heinrich, P., and Roche,~R.: Analytical solution for testing debris avalanche numerical models, Pure Appl. Geophys.,  157, 1081--1096, 2000.


\bibitem[Mazzei and Lettieri(2007)]{mazzei07}
 Mazzei,~L. and Lettieri,~P.: A~drag force closure for uniformly dispersed fluidized suspensions, Chem. Eng. Sci., 62, 6129--6142, 2007.


\bibitem[Mei at al(2001)]{mei01}
 Mei,~K., Liu,~F., and Yuhi,~M.: Mud flow -- slow and fast, Geomorph. Fluid Mech., 22, 548--577, 2001.


\bibitem[Meruane et~al.(2010)]{meruane10}
 Meruane,~C., Tamburrino,~A., and Roche,~O.: On the role of the ambient fluid
on gravitational granular flow dynamics, Journal of Fluid Mechanics, 648, 381--404, 2010.


\bibitem[Neri et~al.(2003)]{neri03}
 Neri,~A., Ongaro,~T.~E., Macedonio,~G., and Gidaspow,~D.: Multiparticle simulation of collapsing volcanic columns and pyroclastic flows,~J. Geophys. Res., 108, 508--531, 2003.


\bibitem[Neri et~al.(2007)]{neri07}
 Neri,~A., Ongaro,~T.~E., Menconi,~G., Vitturu,~M., Cavazzoni,~C., Erbacci,~G., and Baxter,~P.: 4{D} simulation of explosive eruption dynamics at Vesuvius, Geophys. Res. Lett., 34, 1--7, 2007.



\bibitem[Panneerselvan et~al.(2007)]{panneerselvan07}
 Panneerselvan,~R., Savithri,~S., and Surender,~G.~D.: CFD based investigations on hydrodynamics and energy dissipation due to solid motion in liquid fluidized bed, Chem. Eng. J., 132, 159--171, 2007.


\bibitem[Patra et~al.(2005)]{patra05}
 Patra,~A., Bauer,~C., Nichita,~C., Pitman,~E., Sheridan,~M., Bursik,~M., Rupp,~B., Webber,~A., Stinton,~A., Namikawa,~L., and Renschler,~C.: Parallel adaptive numerical simulation of dry avalanches over natural terrain, J. Volcanol. Geoth. Res.,  139, 1--21, 2005.


\bibitem[Pelanti et~al.(2008)]{pelanti08}
 Pelanti,~M., Bouchut,~F., and Mangeney,~A.: A Roe-{t}ype scheme for
two-phase shallow granular flows over variable topography,
Math. Model. Numeric. Analy., 42, 851--885, 2008.


\bibitem[Pistolesi et~al.(2013)]{pistolesi13}
 Pistolesi,~M., Cioni,~R., Rosi,~M., Cashman,~K.~V., Rossotti,~A., and Aguilera,~E.: Evidence for lahar-triggering mechanisms in complex stratigraphic sequences: the post-twelfth century eruptive activity of Cotopaxi Volcano, Ecuador, B. Volcanol.,  75, 698, \href{http://dx.doi.org/10.1007/s00445-013-0698-1}, 2013.


\bibitem[Pitman and Le(2005)]{pitman05}
 Pitman,~E.~B. and Le,~L.: A~two-fluid model for avalanche and debris flows, Philos. T. R. Soc., 363, 1573--1601, 2005.


\bibitem[Procter et~al.(2010)]{procter10}
 Procter,~J.~N., Cronin,~S.~J., Fuller,~I.~C., Sheridan,~M., Neall,~V.~E., and Keys, H.: Lahar hazard assessment using Titan2D for an alluvial fan with rapidly changing geomorphology: Whangaehu River, Mt. Ruapehu, Geomorphology, 116, 162--174, 2010.


\bibitem[Pudasaini and Hutter(2003)]{pudasaini03}
 Pudasaini,~S.~P. and Hutter,~K.: Rapid shear flows of dry granular material down curved and twisted channels,~J. Fluid Mech., 495, 193--208, 2003.


\bibitem[Pudasaini and Hutter(2007)]{pudasain07}
 Pudasaini,~S.~P. and Hutter,~K.: Avalanche Dynamics, Springer-Verlag, Berlin, 2007.


\bibitem[Rankine(1957)]{rankine57}
 Rankine,~W.~J.~M.: On the stability of loose earth, Philos. T. R. Soc. Lond., 147, 9--27, 1957.


\bibitem[Savage and Hutter(1989)]{savage89}
 Savage,~S.~B. and Hutter,~K.: The motion of a~finite mass of granular material down a~rough incline,~J. Fluid Mech., 199, 177--215, 1989.

\bibitem[Schmid(1981)]{schmid81}
Schmid, R., Descriptive nomenclature and classification of pyroclastic deposits and fragments: recommendations of the {IUGS} subcommision on systematics of igneous rocks. Geology, 9, 41-43., 1981.

\bibitem[Sheridan(1979)]{sheridan79}
 Sheridan,~M.: Emplacement of pyroclastic flows: a~review, Geol. S. Am S., 180, 125--136, 1979.


\bibitem[Sheridan(1980)]{sheridan80}
 Sheridan,~M.: Pyroclastic block flow from the September 1976, eruption of {L}a
{S}oufri\`{e}re {V}olcano, {G}uadaloupe, Bull. Volcanol.,  43-2, 397--402,
1980.

\bibitem[Sheridan and Malin(1983)]{sheridan83}
 Sheridan,~M. and Malin,~M.: Application of computer-assisted mapping to volcanic hazard evaluation of surge eruptions: Vulcano, Lipari and Vesuvius, J. Volcanol. Geoth. Res.,  17, 187--202, 1983.


\bibitem[Sheridan et~al.(1999)]{sheridan99}
 Sheridan,~M.~F., Bonnard,~C., Carreno,~R., Siebe,~C., Strauch,~W.,
Navarro,~M., Calero,~J.~C., and Trujillo,~N.~B.: Rockfall/avalanche and
breakout flow of Casita Volcano, Nicaragua, Landslide News, 12, 2-4, 1999.

\bibitem[Sheridan et~al.(2005)]{sheridan05}
 Sheridan,~M.~F., Stinton,~A.~J., Patra,~A.~K., Bauer,~A.~C., Nichita,~C.~C., and Pitman,~E.~B.: Evaluating TITAN2D mass-flow model using the 1963 Little Tahoma Peak avalanches, Mount Rainier, Washington., J. Volcanol. Geoth. Res.,  139, 89--102, 2005.

\bibitem[Sheridan et~al.(2011)]{sheridan11}
 Sheridan,~M., C\'ordoba,~G., Pitman,~E., Cronin,~S., and Procter,~J.:
Application of a~wide-ranging two-phase Debris Flow Model to the 2007 Crater
Lake break-out lahar at Mt. Ruapehu, New Zealand, American Geophysical Union Fall Meeting, San Francisco ({CA}),  5--9~December 2011, abstract V53E-2691, 2011.

\bibitem[swaffield and Bridge(1983)]{swaffield83}
Swaffield, J.A., and Bridge, S., Applicability of the {C}olebrook-{W}hite formula to represent frictional losses in partially filled unsteady pipeflow. J. of research of the National Bureau of Standards, 88, 389-393, 1983.

\bibitem[Terzaghi(1936)]{terzaghi36}
 Terzaghi,~K.: The shearing resistance of saturated soils and the angle
between planes of shear, Proc. 1st Int. Conf. Soil Mech., June 22-26, 1936, Cambridge, {Mass.},  54--56, 1936.

\bibitem[Transport and Road Research Laboratory(1975)]{trrl75}
Transport and Road Research Laboratory, A guide for engineers to the design of storm sewer systems. Roadnote 35, 1975

\bibitem[The-Committee-on-Natural-Disasters(1991)]{nrc91}
 The-Committee-on-Natural-Disasters: The Eruption of Nevado del Ruiz Volcano,
Colombia, South America, November 13, 1985, Natural Disaster Studies, V4, National Academic Press, Washington, D.C., 1991.



\bibitem[Tilling(1996)]{tilling96}
 Tilling,~R.~I.: Volcanoes,~U.S. Geological Survey, U.S. Government Peinting Office, Washington, D.C., 1996.


\bibitem[Valentine(1998)]{valentine98}
 Valentine,~G.: Damage to structures by pyroclastic flows and surges, inferred from nuclear weapons effects, JVGR, 87, 117--140, 1998.

\bibitem[Valentine and Wohletz(1989)]{valentine89b}
 Valentine,~G. and Wohletz,~K.: Numerical models of plinian eruptions columns
and pyroclastic flows, J. Geophys. Res.,  94, 1867--1887, 1989.


\bibitem[Valentine et~al.(2011)]{valentine11}
 Valentine,~G., Doronzo,~D., Dellino,~P., and Tullio,~M.: Effects of volcano profile on dilute pyroclastic density currents: numerical simulations, Geology, 39, 947--950, 2011.

\bibitem[Wallis(1969)]{wallis69}
 Wallis,~G.~B.: One-Dimensional Two-Phase Flow, McGraw Hill,
NY, 1969.


\bibitem[Williams et~al.(2008)]{williams08}
 Williams,~R., Stinton,~A.~J., and Sheridan,~M.~F.: Evaluation of the Titan2{D} two-{p}hase flow model using an actual event: case study of the 2005 Vazc\'un Valley Lahar, J. Volcanol. Geoth. Res.,  177, 760--766, 2008.


\bibitem[Wrachien and Mambretti(2009)]{wrachien09}
 Wrachien,~D. and Mambretti,~S.: Dam-Break Problems, Solutions and Case Studies, WITpress, Southampton, UK, 2009.

\end{thebibliography}


%%
\begin{figure}%
\includegraphics[width=80mm]{nhess-2015-18-discussions-f01.png}
\caption{%
Effect on $k_{ap}$ due to changes in $\phi_{bed}$ and $42^{o} \times \varphi$ given a fix $\phi_{int}$. The black line shows the changes related with $\phi_{bed}$ for {\em any} fixed $\phi_{int}$. The red line shows the relation $k_{ap}$ with  $42^{o} \times \varphi$, used within the program.}
%
\label{fig:kap_bed_phi}
\end{figure}%

\begin{figure}%
\includegraphics[width=80mm]{nhess-2015-18-discussions-f02.png}
\caption{%
Dam-break 1-D problem. Evolution  of the flow front depth with
distance.  The figure shows the result of the theoretical solution and the results from the numerical model for two initial conditions ($h_o=3$ $m$ and $h_o=1.5$ $m$). The statistical {\em t}-test shows that there is no statistical significant difference between the analytical solution and the predictions of Titan2F.}
%
\label{fig:dambreak}
\end{figure}%


\begin{figure}%
\includegraphics[width=120mm]{nhess-2015-18-discussions-f03.png}
\caption{%
Graphic comparison of the predictions of Titan2F for the experiments performed by  \citet{liu96}. The results for several initial volumes compare \textbf{(a)} the width of the deposit and  \textbf{(b)} the length of the deposit. The circles show data from the experiments, and the asterisks the predictions of Titan2F.  }
%
\label{fig:liu_exp}
\end{figure}%


\begin{figure}%
\includegraphics[width=100mm]{nhess-2015-18-discussions-f04.png}
\caption{%
Verification of the accuracy of the Titan2F predictions for the time evolution of the  flow depth and arrival of the front at $32\,\unit{m}$ distance from the lock. The  black line represents the Titan2F results. The gray shadowed interval  represents the confidence interval of the experiments performed by Iverson and co-workers \citep{iverson10}.}
%
\label{fig:iverson1}
\end{figure}%


\begin{figure}%
\includegraphics[width=100mm]{nhess-2015-18-discussions-f05.png}
\caption{%
Verification of the accuracy of the Two-Phase-Titan predictions for the time evolution of the  flow depth and arrival of the front at $66\,\unit{m}$ distance from the lock.  The  black line represents the Titan2F results. The gray shadowed interval represents the confidence interval of the experiments performed by Iverson and co-workers \citep{iverson10}.}
%
\label{fig:iverson2}
\end{figure}%


\begin{figure}%
\includegraphics[width=80mm]{nhess-2015-18-discussions-f06.png}
\caption{%
Longitudinal distribution of the particle concentration after  $14\,\unit{s}$. The particle concentration at the upper part of the channel are very low, whereas at the end, the front of the  flow becomes very concentrated, as observed by \citet{iverson10}. }
%
\label{fig:longitudinal_concentration}
\end{figure}%


\begin{figure}%
\includegraphics[width=80mm]{nhess-2015-18-discussions-f07.png}
\caption{%
Longitudinal distribution of the predicted dynamic pressure at $14\,\unit{s}$. The peak over-pressure occurs just at the end of the inclined part of the channel.}
%
\label{fig:dpress}
\end{figure}


\end{document}
\endinput%%%%%%
