
\title{GRASS GIS and Modelling of Natural Hazards}
\subtitle{An Integrated Approach for Debris Flow Simulation --- First results of an application in the Central Andes}
\author{Martin Mergili and Wolfgang Fellin}

\maketitle

\section{Background}
\label{sec:sbackground}
Debris flows are rapid mass movements of water and debris, constituting a considerable hazard when interfering with people, buildings, or infrastructure.
They are often triggered by heavy or prolonged rainfall or by extreme snow melt. Mobilization of the material occurs due to translational or rotational failure of saturated or undercut slopes, or by detachment due to surface runoff or the debris flow itself. Various models do exist for simulating sub-processes included into debris flows, for example for detachment (\textit{r.sim.sediment} within the GRASS GIS environment), for soil hydrology and slope stability \citep{xieetal:2004}, or for debris flow runout \citep{savagehutter:1989, pudasainihutter:2007}. More integrated GIS-based approaches as attempted for example by \cite{burtonbathurst:1998} or \cite{wichmann:2006} are scarce. Such approaches would be valuable for a quick assessment of hydrological thresholds for potential debris flow hazard regarding specified features at risk. This paper describes and discusses the development of such a model as GRASS GIS raster module. The model is designed for small catchments (few square kms) and is tested at the moment with seven study areas along the international road corridor from Mendoza (Western Argentina) to Central Chile, crossing the highest section of the Andes (figure~\ref{fig:mergili_fig01}). The preliminary results for the study area \textit{Guido~A} are presented.

\begin{figure*}[htb]
	\centering
	\includegraphics[width=0.8\textwidth]{./proceedings/mergili/mergili_fig01.jpg}
	\caption{Study areas. The preliminary results for Guido~A are presented.}
	\label{fig:mergili_fig01}
\end{figure*}

\section{Model}
\label{sec:smodel}

\subsection{Implementation and model design}
\label{subsec:ssimpdesign}
The simulation model is implemented as a GRASS GIS raster module called \textit{r.debrisflow}, based on the C programming language. Data management is facilitated using shell scripts. The model is in an intermediate stage of development right now, with major technical and methodical enhancements prospected for the near future. Additionally, a GUI for data management
shall be created. By now, the latest development version can be downloaded from the homepage of the first author. \textit{r.debrisflow} constitutes of a framework of a number of sub-modules described in more detail below, the general model design
is illustrated in figure~\ref{fig:mergili_fig02}. The sub-modules can be combined in two different ways, depending on the availability of input information:

\begin{description}
   \item \textbf{Simulation mode 1:} The entire hydrological, stability, detachment and runout modelling is executed for a defined number of time steps during a rainfall or snowmelt event, requiring an extensive set of information as input, including meteorological data, an elevation model, soil mechanical and hydrological parameters and surface hydrological characteristics (including land cover).
	\item \textbf{Simulation mode 2:} The zones of debris flow initiation are defined manually (e.g. from mapping in the field or from orthophotos), and only runout is computed. The advantage of this mode is that it requires much less input than the others, but, on the other hand, it is not suitable for predicting future events.
\end{description}

\begin{figure}[htbp]
	\centering
	\includegraphics[scale=0.19]{./proceedings/mergili/mergili_fig02.jpg}
	\caption{General model design.}
	\label{fig:mergili_fig02}
\end{figure}

\subsection{Water input}
\label{subsec:sswaterinput}
Rainfall is read from the input file and added to the system by increasing the surface water table of each cell, reduced by interception. If a snow cover exists, snowmelt is computed for each cell with a user defined degree-day-factor and added to the surface water table.

\subsection{Soil hydrology and slope stability}
\label{subsec:sssoilhydroslopestab}
For this sub-module, a three-dimensional raster approach is used, down to the depth of bedrock (if known), or to a user-defined maximum soil depth. The soil is assumed to be homogeneous over its entire depth regarding its physical, hydrological and mechanical properties. Vertical flow between cells is computed with the Darcy-Equation. If the water content of a cell exceeds 90 \% of the maximum content, groundwater flow is assumed to be parallel to the slope and it is tested whether the cell is stable or not, using an infinite slope stability approach \citep{xieetal:2004}. For each pixel, the bottom of the deepest cell with a factor of safety lower than 1 is considered as failure plane (figure~\ref{fig:mergili_fig03}).
It has to be pointed out that this approach constitutes a rough approximation to the reality with the character of a worst-case assumption: the stabilizing role of vertical water movement is neglected, and the destabilizing role of the assumed slope-parallel component is fully included in case of saturation. In the real world, both components are combined, resulting in more stable conditions than in the model.

\begin{figure*}[htb]
	\centering
	\includegraphics[width=0.7\textwidth]{./proceedings/mergili/mergili_fig03.jpg}
	\caption{Subsurface hydrology and slope stability model.}
	\label{fig:mergili_fig03}
\end{figure*}

\subsection{Surface runoff and detachment}
\label{subsec:sssurfdet}
All the water which is not infiltrated into the soil is added to the surface water table of the corresponding cell. Flow velocity from one cell to the other is computed using the Manning equation. If no flow channel is defined for any downward cell, flow takes place to all downward cells, weighted for slope angle. If a channel is defined, the entire flow moves there. Transport capacity is computed using the \cite{yang:1973} equation, which it is supposed to suit best for the conditions in the study areas \citep{obrien:2003}. Rates of detachment are derived from the transport capacity. The model design is illustrated in figure~\ref{fig:mergili_fig04}.

\begin{figure*}
	\centering
	\includegraphics[width=0.8\textwidth]{./proceedings/mergili/mergili_fig04.jpg}
	\caption{Surface runoff and detachment model.}
	\label{fig:mergili_fig04}
\end{figure*}

\subsection{Debris flow runout and deposition}
\label{subsec:ssrunoutdep}
Debris flows are here understood in a strict sense, with a non-hydraulic flow regime and excluding heavily sediment-loaded water discharge. For this
reason, the sediment concentration is computed for the mobilized soil volume of each cell:

\begin{enumerate}
	\item If the sediment concentration exceeds a threshold value \citep{obrien:2003}, the entire mobilized volume is considered to develop into a debris flow. If mobilization occurred due to slope failure, the cell is marked and runout of all the failed cells is computed at the end of the event. If mobilization occurred due to detachment, runout is calculated immediately at
the end of the time step. This is a first approximation as the empirical runout models used (compare below) neglect the time required for runout.
	\item If the sediment concentration is too low for the development of a debris flow, the material is removed with surface runoff and either deposited downslope or removed from the system, following the \cite{yang:1973} equation (compare figure~\ref{fig:mergili_fig04}). Though deposition from surface runoff does not fall into the concept of debris flow, it can provide valuable complementary information and is therefore regarded, too.
\end{enumerate}

Debris flow runout itself can be simulated using physically-based models \citep{savagehutter:1989, pudasainihutter:2007}, but they are complex and hard to be integrated into a GIS environment. Therefore it was decided to use a combination of empirical approaches first for estimating the runout distance and the distribution of the deposited volume (figure~\ref{fig:mergili_fig05}). \cite{vandre:1985} developed an approach enabling the distinction between scouring and deposition areas, using threshold slope angles and the ratio between vertical distance of scouring and horizontal distance of deposition. \cite{corominasetal:2003} and \cite{rickenmann:1999} used mobilized volume, angles, and runout distances for estimating the reach of the debris flow. These approaches, however, have the disadvantage that they don't allow distinguishing between areas of scouring and deposition, and therefore the distribution of the deposited material. The approach was included into the model as follows:

\begin{enumerate}
	\item The Vandre approach was applied with user-defined parameters for estimating distributed scoured and deposited volumes.
	\item The Corominas et al. and the Rickenmann approach
were applied independently and then combined to an index.
\end{enumerate}

The debris flow is routed downwards separately for each unstable cell, following a random walk \citep{hughes:1995} weighted for slope angle and the existence of a defined channel, until the stop criteria for all of the three approaches is fulfilled. Though each cell is treated separately, the mobilized volume required for runout distance is calculated for each patch of debris flow initiation. In the area of scouring, the entire saturated soil column is considered to be removed, but never exceeding the depth of initiation. The initiated and scoured volumes are considered to distribute over the area of deposition as wedge shape rising towards the front.

\begin{figure}
	\centering
	\includegraphics[width=1\columnwidth]{./proceedings/mergili/mergili_fig05.jpg}
	\caption{Runout models.}
	\label{fig:mergili_fig05}
\end{figure}

\section{First results}
\label{sec:sfirstresults}
The model was tested within the study area Guido~A (compare figure~\ref{fig:mergili_fig01}). Mainly consisting of granite residuals, the soils of the catchment (area: 2 km$^{2}$) are relatively homogeneous. Therefore it was decided to use one single set of soil parameters for the entire area:
\newline
\newline
\begin{tabular}{|c|c|c|c|c|c|}
\hline
texture & $\rho_d$ & $c_s$ & $\varphi $ & $\Theta_s$ & $k_f$ \\
& kg/m$^3$ & N/m$^2$ & deg. & & cm/h \\
\hline
Sand & 1850 & 0 & 40.0 & 0.43 & 29.7 \\
\hline
\end{tabular}
\newline
\newline
$\rho _d$ is the dried bulk density of the soil, $c_s$ is soil cohesion, $\varphi$ stands for the angle of internal friction, $\Theta _s$ is the maximum (saturated) water content, and $k_f$ is the saturated hydraulic conductivity.

The figures ~\ref{fig:mergili_fig06} and \ref{fig:mergili_fig07} illustrate the mapped areas of debris flow initiation in the study area Guido~A, and the patterns of surface change due to a debris flow event, based on the mapped areas of initiation and the computed patterns of scouring and deposition (simulation mode 2). The white line crossing the right part of the maps represents the international road, roughly coinciding with the distal part of the observed debris flow depositions.
%\newline
The figures ~\ref{fig:mergili_fig08} to \ref{fig:mergili_fig11} show some of the simulation results for a hypothetical 100~mm rainfall event, corresponding to the maximum daily sum ever recorded at the nearby meteorological station, and therefore constituting a worst-case assumption (simulation mode 1). All maps show plausible patterns when compared to field observations. The areas of debris flow initiation and deposition are located correctly, but are overestimated compared to the patterns observed in the field (what is not surprising for a worst-case assumption).
%\newline
The calculated sediment volumes deposited on the international road are within the same magnitude as those reported by the road authorities. When simulating the impacts of smaller rainfall events, the model results correspond well to the findings of \cite{moreiras:2005} that debris flows in the Mendoza valley usually occur at daily rainfall sums exceeding 6.6 to 12.9~mm.

%%% testing

\begin{figure}[htbp]
%	\begin{minipage}[t]{0.475\linewidth}
		\centering
		\includegraphics[width=1\textwidth]{./proceedings/mergili/mergili_fig06.jpg}
		\caption{Mapped areas of clearly identifiable previous debris flow initiation, depth of initiation assumed as 0.75~m according to field evidence.}
		\label{fig:mergili_fig06}	
%	\end{minipage}
%	\hspace{0.05\textwidth} 
%	\begin{minipage}[t]{0.475\linewidth}
\end{figure}
\begin{figure}[htbp]
		\centering
		\includegraphics[width=1\textwidth]{./proceedings/mergili/mergili_fig07.jpg}
		\caption{Simulated change of terrain height due to debris flow, using mapped areas of initiation.}
		\label{fig:mergili_fig07}
%	\end{minipage}
%\end{figure*}
\end{figure}

\begin{figure*}[htbp]
	\begin{minipage}[t]{0.475\linewidth}
	\centering
		\includegraphics[width=1\textwidth]{./proceedings/mergili/mergili_fig08.jpg}
		\caption{Maximum depth of saturation computed for 100~mm rainfall event.}
		\label{fig:mergili_fig08}
	\end{minipage}
	\hspace{0.05\textwidth} 
	\begin{minipage}[t]{0.475\linewidth}
		\centering
		\includegraphics[width=1\textwidth]{./proceedings/mergili/mergili_fig09.jpg}
		\caption{Simulated areas of potential debris flow initiation computed for 100~mm rainfall event.}
		\label{fig:mergili_fig09}
	\end{minipage}
\end{figure*}

\begin{figure*}[htbp]
	\begin{minipage}[t]{0.475\linewidth}
		\centering
		\includegraphics[width=1\textwidth]{./proceedings/mergili/mergili_fig10.jpg}
		\caption{Simulated change of terrain height due to debris flow caused by 100~mm rainfall event.}
		\label{fig:mergili_fig10}
	\end{minipage}
	\hspace{0.05\textwidth} 
	\begin{minipage}[t]{0.475\linewidth}
		\centering
		\includegraphics[width=1\textwidth]{./proceedings/mergili/mergili_fig11.jpg}
		\caption{Debris flow index for 100~mm rainfall event, denoting number of cells the mobilized material of which hits the pixel.}
		\label{fig:mergili_fig11}
	\end{minipage}
\end{figure*}

\section{Discussion and Preview}
\label{sec:sdiscprec}
Though the preliminary results for Guido A appear plausible, the simulation model still shows a number of shortcomings that have to be reworked.

\begin{enumerate}
   \item Infiltration of water into the soil is not yet modelled in a satisfactory way, so that the approach will have to be refined (Green-Ampt model). Slope-parallel groundwater flow will be included, too, for enabling a closer approximation to the reality of soil hydrology.
	\item The slope stability model as applied at the moment is only valid for plane, infinite, cohesionless slopes. For very shallow failures, this assumption is sufficiently close to reality, but for more deep-seated rotational failures, it is unsatisfactory and slope curvature has to be taken into account. \cite{wilkinsonetal:2002} and \cite{xieetal:2003} could serve as examples for such an approach.
	\item The empirical approaches for debris flow runout shall be complemented by the implementation of a physically-based runout model according to \cite{savagehutter:1989} and \cite{pudasainihutter:2007}, or at least of an interface to a non GIS-based runout model.
	\item An interesting extension would be to introduce some probabilistic elements into the slope stability model (regarding the SINMAP model as an example) and into the distinction rules between debris flow and other types of movements. 
	\item Finally, the model has to be applied to the remaining study areas (compare figure ~\ref{fig:mergili_fig01}), and the results have to be tested carefully against the field observations and the validation data (reports about volumes of material deposited on the international road).
\end{enumerate}
With the mentioned optimizations, r.debrisflow shall be a valuable tool for evaluating the potential magnitude of debris flows as a response to defined rainfall or snow melt events, including the possibility to determine meteorological thresholds for debris flow hazard. However, it has to be pointed out that all the results only denote potential occurrences with the character of worst case scenarios or probabilities - it will probably never be possible to predict the actual response of a slope to a meteorological event in the real world, as nature is too complex to be fully understood in all details.


\section{Acknowledgements}

\begin{footnotesize}
The project is funded by the \textit{Doktoratsstipendium aus der Nachwuchsf\"orderung der LFU} and by grants from the Office for International Relations and the Faculties for Natural Sciences, University of Innsbruck, as well as from the Austrian Academy of Sciences and the Federal Government of Upper Austria. Thanks to Hans Stoetter, Axel Borsdorf, Clemens Geitner and Stella Moreiras for their support.
\end{footnotesize}

\begin{footnotesize}
\begin{thebibliography}{99}

\bibitem[Burton and Bathurst (1998), Burton, Bathurst]{burtonbathurst:1998}
A. Burton, J.C. Bathurst (1998)
\newblock Physically based modelling of shallow landslide sediment yield at a catchment scale.
\newblock {\em Environmental Geology} 35,2-3: 89-99.

\bibitem[Corominas et~al. (2003), Corominas, Copons, Vilaplana, Altamir, Amig\'o]{corominasetal:2003}
J. Corominas, R. Copons, J.M. Vilaplana, J. Altamir, J. Amig\'o (2003)
\newblock Integrated Landslide Susceptibility Analysis and Hazard Assessment in the Principality of Andorra.
\newblock {\em Natural Hazards} 30: 421-435.

\bibitem[Hughes (1995) Hughes]{hughes:1995}
B. D. Hughes (1995)
\newblock Random walks and random environments: Volume 1: Random Walks.
\newblock Oxford University Press. 652 pp.

\bibitem[Li (2007) Li]{li:2007}
X. Li (2007)
\newblock Finite element analysis of slope stability using a nonlinear
failure criterion.
\newblock {\em Computers and Geotechnics} 34: 127-136.

\bibitem[Moreiras (2005) Moreiras]{moreiras:2005}
S.M. Moreiras (2005)
\newblock Climatic effect of ENSO associated with landslide occurrence in the Central Andes, Mendoza Province, Argentina.
\newblock {\em Landslides} 2: 53-59.

\bibitem[O'Brien (2003) O'Brien]{obrien:2003}
J.S. O'Brien (2003)
\newblock FLO-2D Users' Manual Version 2003.06, July 2003.
\newblock FLO-2D Software Inc., Nutrioso, Arizona, USA. 232 pp.

\bibitem[Pudasaini and Hutter (2007) Pudasaini, Hutter]{pudasainihutter:2007}
S.P. Pudasaini, K. Hutter (2007)
\newblock Avalanche Dynamics: Dynamics of rapid flows of dense granular avalanches.
\newblock Springer, Berlin Heidelberg. 602 pp.

\bibitem[Rickenmann (1999) Rickenmann]{rickenmann:1999}
D. Rickenmann (1999)
\newblock Empirical Relationships for Debris Flows.
\newblock {\em Natural Hazards} 19: 47-77.

\bibitem[Savage and Hutter (1989) Savage, Hutter]{savagehutter:1989}
S.B. Savage, K. Hutter (1989)
\newblock The motion of a finite mass of granular material down a rough incline. 
\newblock {\em Journal of Fluid Mechanics} 199: 177-215.

\bibitem[Vandre (1985) Vandre]{vandre:1985}
B.C. Vandre (1985)
\newblock Rudd Creek debris flow.
\newblock In: D.S. Bowles (ed): Delineation of landslide, flash flood, and debris flow hazards in Utah. Utah Water Research Laboratory, Utah State University, Logan, Utah, 117-131.

\bibitem[Wichmann (2006) Wichmann]{wichmann:2006}
V. Wichmann (2006)
\newblock Modellierung geomorphologischer Prozesse in einem alpinen Einzugsgebiet - Abgrenzung und Klassifizierung der Wirkungsr\"aume von Sturzprozessen und Muren mit einem GIS.
\newblock Eichstaeter geographische Arbeiten 15: 231 pp. In German.

\bibitem[Wilkinson et~al. (2002)]{wilkinsonetal:2002}
P.L. Wilkinson, M.G. Anderson, D.M. Lloyd (2002)
\newblock An integrated hydrological model for rain-induced landslide prediction.
\newblock {\em Earth Surface Processes and Landforms} 27: 1285-1297.

\bibitem[Xie et~al. (2003) Xie, Esaki, Zhou, Mitani]{xieetal:2003}
M. Xie, T. Esaki, G. Zhou, Y. Mitani (2003)
\newblock Three-dimensional stability evaluation of landslides and a sliding process simulation using a new geographic information systems component.
\newblock {\em Environmental Geology} 43: 503-512.

\bibitem[Xie et~al. (2004) Xie, Esaki, Cai]{xieetal:2004}
M. Xie, T. Esaki, M. Cai (2004)
\newblock A time-space based approach for mapping rainfall-induced shallow landslide hazard.
\newblock {\em Environmental Geology} 46: 840-850.

\bibitem[Yang (1973) Yang]{yang:1973}
C.T. Yang (1973)
\newblock Incipient motion and sediment transport.
\newblock {\em Journal of the Hydraulics Division. Proceedings of the American Society of Civil Engineers} 99: 1679-1703.

\end{thebibliography}
\end{footnotesize}

\address{Martin Mergili\\
  University of Innsbruck, Institute of Geography\\
  and Austrian Academy of Sciences, Mountain Research\\
  \url{www.uibk.ac.at/geographie/personal/mergili}\\
  \email{martin.mergili AT uibk.ac.at}}

\address{Wolfgang Fellin\\
  University of Innsbruck, Institute of Infrastructure\\
  \url{www.uibk.ac.at/geotechnik/staff/fellin-en}\\
  \email{wolfgang.fellin AT uibk.ac.at}}

