source: anuga_work/publications/boxing_day_validation_2008/patong_validation.tex @ 7303

\documentclass[a4paper]{article}
% to be added when submitted to ocean dynamics
%\documentclass[smallcondensed,draft]{svjour3}
%\usepackage{mathptmx}
%\journalname{Ocean Dynamics}

\usepackage{graphicx}
\usepackage{hyperref}
\usepackage{amsfonts}
\usepackage{url}      % for URLs and DOIs
\newcommand{\doi}[1]{\url{http://dx.doi.org/#1}}


%----------title-------------%
\title{Benchmarking Tsunami Models using the December 2004 Indian
  Ocean Tsunami and its Impact at Patong Beach}

%-------authors-----------
\author{J.~D. Jakeman \and O. Nielsen \and K. VanPutten \and
  D. Burbidge \and R. Mleczko \and N. Horspool}

% to be added when submitted to ocean dynamics
%\institute{J.~D. Jakeman \at 
%       The Australian National University, Canberra, \textsc{Australia}\
%       \email{john.jakeman@anu.edu.au}
%       \and
%       O. Nielsen \and R. Mleczko \and D. Burbidge \and K. VanPutten \and N. Horspool \at 
%       Geoscience Australia, Canberra, \textsc{Australia}
%}


%================Start of Document================
\begin{document}
\maketitle
%------Abstract--------------
\begin{abstract}
In this paper a new benchmark for tsunami model validation is
proposed. The benchmark is based upon the 2004 Indian Ocean tsunami,
which affords a uniquely large amount of observational data for model
comparison. Unlike the small number of existing benchmarks, the
proposed test validates all three stages of tsunami evolution -
generation, propagation and inundation. Specifically we use geodetic
measurements of the Sumatra--Andaman earthquake to validate the
tsunami source, altimetry data from the \textsc{jason} satellite to
test open ocean propagation, eye-witness accounts to assess near shore
propagation and a detailed inundation survey of Patong Bay, Thailand
to compare model and observed inundation. Furthermore we utilise this
benchmark to further validate the hydrodynamic modelling tool
\textsc{anuga} which is used to simulate the tsunami
inundation. Important buildings and other structures were incorporated
into the underlying computational mesh and shown to have a large
influence on inundation extent. Sensitivity analysis also showed that
the model predictions are comparatively insensitive to large changes
in friction and small perturbations in wave weight at the 100 m depth
contour.
% to be added when submitted to ocean dynamics
%\keywords{Tsunami \and modelling \and validation and verification \and benchmark}
\end{abstract}

\tableofcontents
%================Section===========================

\section{Introduction}
Tsunami are a potential hazard to coastal communities all over the
world. A number of recent large events have increased community and
scientific awareness of the need for effective detection, forecasting,
and emergency preparedness. Probabilistic, geological, hydrodynamic,
and economic models are required to predict the location and
likelihood of an event, the initial sea floor deformation and
subsequent propagation and inundation of the tsunami, the
effectiveness of hazard mitigation procedures and the economic impact
of such measures and of the event itself. Here we focus on modelling of
the physical processes. 
%OLE: I commented this out 23 June 2009 as there was no reference.
%For discussion on economic and decision based
%models refer to~\cite{} and the references therein.

Various approaches are currently used to assess the potential impact
of tsunami. These methods differ in both the formulation used to
describe the evolution of the tsunami and the numerical methods used
to solve the governing equations. However any legitimate model must
address each of the three distinct stages of tsunami evolution---
generation, propagation and inundation. Models combining observed seismic, 
geodetic and sometimes tsunami data must be used
to provide estimates of initial sea floor and ocean surface
deformation. The complexity of these models ranges from empirical to
non-linear three-dimensional mechanical models. The shallow water wave
equations, linearised shallow water wave equations, and
Boussinesq-type equations are frequently used to simulate tsunami
propagation. These models are typically used to predict quantities
such as arrival times, wave speeds and heights, and inundation extents
which are used to develop efficient hazard mitigation plans.

Inaccuracies in model prediction can result in inappropriate
evacuation plans and town zoning, which may result in loss of life and
large financial losses. Consequently tsunami models must undergo
sufficient end-to-end testing to increase scientific and community
confidence in the model predictions.

Complete confidence in a model of a physical system cannot be
established.  One can only hope to state under what conditions and to what extent the
model hypothesis holds true. Specifically the utility of a model can
be assessed through a process of verification and
validation. Verification assesses the accuracy of the numerical method
used to solve the governing equations and validation is used to
investigate whether the model adequately represents the physical
system~\cite{bates01}. Together these processes can be used to
establish the likelihood that a model represents a legitimate
hypothesis.

The sources of data used to validate and verify a model can be
separated into three main categories: analytical solutions, scale
experiments and field measurements. Analytical solutions of the
governing equations of a model, if available, provide the best means
of verifying any numerical model. However, analytical solutions are
frequently limited to a small set of idealised examples that do not
completely capture the more complex behaviour of `real' events. Scale
experiments, typically in the form of wave-tank experiments, provide a
much more realistic source of data that better captures the complex
dynamics of flows such as those generated by a tsunami, whilst allowing
control of the event and much easier and accurate measurement of the
tsunami properties. Comparison of numerical predictions with field
data provides the most stringent test. The use of field data increases
the generality and significance of conclusions made regarding model
utility. On the other hand, it must be noted that the use of field
data also significantly increases the uncertainty of the validation
experiment that may constrain the ability to make unequivocal
statements~\cite{bates01}.

Currently, the extent of tsunami-related field data is limited. The
cost of tsunami monitoring programs, bathymetry and topography surveys
prohibits the collection of data in many of the regions in which
tsunamis pose greatest threat. The resulting lack of data has limited
the number of field data sets available to validate tsunami
models. Synolakis et al~\cite{synolakis07} have developed a set of
standards, criteria and procedures for evaluating numerical models of
tsunami. They propose three analytical solutions to help identify the
validity of a model, and five scale comparisons (wave-tank benchmarks)
and two field events to assess model veracity.

The first field data benchmark introduced in \cite{synolakis07} compares model
results against observed data from the Hokkaido-Nansei-Oki tsunami
that occurred around Okushiri Island, Japan on the 12 July
1993. This tsunami provides an example of extreme run-up generated from
reflections and constructive interference resulting from local
topography and bathymetry. The benchmark consists of two tide gauge
records and numerous spatially-distributed point sites at which
modelled maximum run-up elevations can be compared. The second
benchmark is based upon the Rat Islands tsunami that occurred off the
coast of Alaska on the 17 November 2003. The Rat Island tsunami
provides a good test for real-time forecasting models since the tsunami
was recorded at three tsunameters. The test requires matching the
propagation model data with the DART recording to constrain the
tsunami source model, and then using it to reproduce the tide gauge
record at Hilo.

In this paper we develop a field data benchmark to be used in
conjunction with the other tests proposed by Synolakis et
al~\cite{synolakis07} to validate and verify tsunami models. 
The benchmark proposed here allows evaluation of
model structure during all three distinct stages tsunami evolution.
It consists of geodetic measurements of the
Sumatra--Andaman earthquake that are used to validate the description
of the tsunami source, altimetry data from the JASON satellite to test
open ocean propagation, eye-witness accounts to assess near shore
propagation, and a detailed inundation survey of Patong Bay, Thailand
to compare model and observed inundation. A description of the data
required to construct the benchmark is given in
Section~\ref{sec:data}.

An associated aim of this paper is to illustrate the use of this new
benchmark to validate an operational tsunami inundation model called
\textsc{anuga} used by Geoscience Australia. A description of
\textsc{anuga} is given in Section~\ref{sec:models} and the validation
results are given in Section~\ref{sec:results}.

The numerical models used to model tsunami are extremely
computationally intensive. Full resolution models of the entire
evolution process will often take a number of days to
run. Consequently the uncertainty in model predictions is difficult to
quantify. However model uncertainty should not be ignored. Section
~\ref{sec:sensitivity} provides a simple sensitivity analysis that can
be used to investigate the sensitivity of model predictions to model
parameters.

%================Section===========================
\section{Data}\label{sec:data}
The sheer magnitude of the 2004 Sumatra-Andaman earthquake and the
devastation caused by the subsequent tsunami have generated much
scientific interest. As a result an unusually large amount of post
seismic data has been collected and documented. Data sets from
seismometers, tide gauges, \textsc{gps} surveys, satellite overpasses,
subsequent coastal field surveys of run-up and flooding, and
measurements of coseismic displacements and bathymetry from ship-based
expeditions, have now been made
available. %~\cite{vigny05,amnon05,kawata05,liu05}.  
In this section we present the corresponding data necessary to implement 
the proposed benchmark for each of the three stages of the tsunami's evolution.

\subsection{Generation}\label{sec:gen_data}
All tsunami are generated from an initial disturbance of the ocean
which develops into a low frequency wave that propagates outwards from
the source. The initial deformation of the water surface is most
commonly caused by coseismic displacement of the sea floor, but
submarine mass failures, landslides, volcanoes or asteroids can also
cause tsunami. In this section we detail the information used in
this study to validate models of the sea floor deformation generated
by the 2004 Sumatra--Andaman earthquake.

The 2004 Sumatra--Andaman tsunami was generated by a severe coseismic
displacement of the sea floor as a result of one of the largest
earthquakes on record. The mega-thrust earthquake started on the 26
December 2004 at 0h58'53'' UTC (or just before 8 am local time)
approximately 70 km offshore of North Sumatra
(\url{http://earthquake.usgs.gov/eqcenter/eqinthenews/2004/usslav}). The
rupture propagated 1000-1300 km along the Sumatra-Andaman trench to
the north at a rate of 2.5-3 km.s$^{-1}$ and lasted approximately 8-10
minutes~\cite{ammon05}. Estimates of the moment magnitude of this
event range from about 9.1 to 9.3 $M_w$~\cite{chlieh07,stein07}.

The unusually large surface deformation caused by this earthquake
means that there were a range of different geodetic measurements of
the surface deformation available. These include field measurements of
uplifted or subsided coral heads, continuous or campaign \textsc{GPS}
measurements and remote sensing measurements of uplift or subsidence
(see~\cite{chlieh07} and references therein). Here we use the the near-field
estimates of vertical deformation in northwestern Sumatra and
the Nicobar-Andaman islands collated by~\cite{chlieh07} to validate
that our crustal deformation model of the 2004 Sumatra--Andaman
earthquake is producing reasonable results. Note that the geodetic
data used here is a combination of the vertical deformation that
happened in the $\sim$10 minutes of the earthquake plus the
deformation that followed in the days following the earthquake before
each particular measurement was actually made (typically of order
days). Therefore some of the observations may not contain the purely
co-seismic deformation but could include some post-seismic deformation
as well~\cite{chlieh07}.

%DAVID: I commented out the figure since we can combine it with the model result without obscuring it. That will keep the number of figures down.

%\begin{figure}[ht]
%\begin{center}
%\includegraphics[width=8.0cm,keepaspectratio=true]{geodeticMeasurements.jpg}
%\caption{Near field geodetic measurements used to validate tsunami generation. FIXME: Insert appropriate figure here}
%\label{fig:geodeticMeasurements}
%\end{center}
%\end{figure}

\subsection{Propagation}
Once generated, a tsunami will propagate outwards from the source until
it encounters the shallow water bordering coastal regions. This period
of the tsunami evolution is referred to as the propagation stage. The
height and velocity of the tsunami is dependent on the local
bathymetry in the regions through which the wave travels and the size
of the initial wave. This section details the bathymetry data needed
to model the tsunami propagation and the satellite altimetry transects
used here to validate open ocean tsunami models.

\subsubsection{Bathymetry Data}
A number of raw data sets were obtained, analysed and checked for
quality and subsequently gridded for easier visualisation and input
into the tsunami models. The resulting grid data is relatively coarse
in the deeper water and becomes progressively finer as the distance to
Patong Bay decreases.

The nested bathymetry grid was generated from: 
\begin{itemize}
\item a two arc minute grid data set covering the Bay of Bengal,
  DBDB2, obtained from US Naval Research Labs;
\item a 3 second arc grid covering the whole of the Andaman Sea based
  on Thai Navy charts no. 45 and no. 362; and
\item a one second grid created from the digitised Thai Navy
  bathymetry chart, no. 358, which covers Patong Bay and the
  immediately adjacent regions.
\end{itemize}

The final bathymetry data set consists of four nested grids obtained
via interpolation and resampling of the aforementioned data sets. The
four grids are shown in Figure~\ref{fig:nested_grids}.  The coarsest
bathymetry was obtained by interpolating the DBDB2 grid to a 27 second
arc grid. A subsection of this region was then replaced by nine second
data which was generated by sub-sampling the three second of arc grid from
NOAA. A subset of the nine second grid was replaced by the three second
data. Finally, the one second grid was used to approximate the
bathymetry in Patong Bay and the immediately adjacent regions. Any
points that deviated from the general trend near the boundary were
deleted.

The sub-sampling of larger grids was performed by using {\bf resample},
a Generic Mapping Tools (\textsc{GMT}) program (\cite{wessel98}). The
gridding of data was performed using {\bf Intrepid}, a commercial
geophysical processing package developed by Intrepid Geophysics. The
gridding scheme employed the nearest neighbour algorithm followed by
an application of minimum curvature akima spline smoothing. 
See \url{http://www.intrepid-geophysics.com/ig/manuals/english/gridding.pdf} 
for details on the Intrepid model.


\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.75\textwidth,keepaspectratio=true]{nested_grids}
\caption{Nested grids of the elevation data.}
\label{fig:nested_grids}
\end{center}
\end{figure}

\subsubsection{JASON Satellite Altimetry}\label{sec:data_jason}
During the 26 December 2004 event, the Jason satellite tracked from
north to south and over the equator at 02:55 UTC nearly two hours
after the earthquake \cite{gower05}. The satellite recorded the sea
level anomaly compared to the average sea level from its previous five
passes over the same region in the 20-30 days prior.  This data was
used to validate the propagation stage in Section
\ref{sec:resultsPropagation}.
%DB I suggest we combine with model data to reduce the number of figures. The satellite track is shown in Figure~\ref{fig:satelliteTrack}.

%\begin{figure}[ht]
%\begin{center}
%\includegraphics[width=8.0cm,keepaspectratio=true]{sateliteTrack.jpg}
%\caption{URS wave heights 120 minutes after the initial earthquake with the JASON satellite track and its observed sea level anomalies overlaid. Note the URS data has not been corrected for the flight path time. FIXME: should we just have track and not URS heights.}
%\label{fig:satelliteTrack}
%\end{center}
%\end{figure}

%\begin{figure}[ht]
%\begin{center}
%\includegraphics[width=8.0cm,keepaspectratio=true]{jasonAltimetry.jpg}
%\caption{JASON satellite altimetry seal level anomaly. FIXME: should we include figure here with just JASON altimetry.}
%\label{fig:jasonAltimetry}
%\end{center}
%\end{figure}

%FIXME: Can we compare the urs model against the TOPEX-poseidon satellite as well? DB No (we don't have the data currently).

\subsection{Inundation}
Inundation refers to the final stages of the evolution of a tsunami and
covers the propagation of the tsunami in shallow coastal water and the
subsequent run-up onto land. This process is typically the most
difficult of the three stages to model due to thin layers of water
flowing rapidly over dry land.  Aside from requiring robust solvers
which can simulate such complex flow patterns, this part of the
modelling process also requires high resolution and quality elevation
data which is often not available. In the case of model validation
high quality field measurements are also required. For the proposed
benchmark the authors have obtained a high resolution bathymetry and
topography data set and a high quality inundation survey map from the
CCOP in Thailand (\cite{szczucinski06}) to validate model inundation.

The datasets necessary for reproducing the results
of the inundation stage are available on Sourceforge under the \textsc{anuga}
project (\url{http://sourceforge.net/projects/anuga}). At the time of
writing the direct link is \url{http://tinyurl.com/patong2004-data}.
%
%\url{http://sourceforge.net/project/showfiles.php?group_id=172848&package_id=319323&release_id=677531}.

In this section we also present eye-witness accounts which can be used to qualitatively validate tsunami inundation.

\subsubsection{Topography Data}
A one second grid was used to approximate the topography in Patong
Bay. This elevation data was again created from the digitised Thai
Navy bathymetry chart, no 358. A visualisation of the elevation data
set used in Patong Bay is shown in
Figure~\ref{fig:patong_bathymetry}. The continuous topography is an
interpolation of known elevation measured at the coloured dots.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=8.0cm,keepaspectratio=true]{patong_bay_data.jpg}
\caption{Visualisation of the elevation data set used in Patong Bay showing data points, contours, rivers and roads draped over the final model.}
\label{fig:patong_bathymetry}
\end{center}
\end{figure}

\subsubsection{Buildings and Other Structures}
Human-made build and structures can significantly effect tsunami
inundation. The location and size and number of floors of the
buildings in Patong Bay were extracted from a GIS data set provided by
the CCOP in Thailand (see acknowledgements at the end of this paper). 
The heights of these
buildings were estimated assuming that each floor has a height of 3 m.

\subsubsection{Inundation Survey}
Tsunami run-up is often the cause of the largest financial and human
losses, yet run-up data that can be used to validate model run-up
predictions is scarce. Of the two field benchmarks proposed in~\cite{synolakis07},
only the Okushiri benchmark facilitates comparison between
modelled and observed run-up. One of the major strengths of the
benchmark proposed here is that modelled run-up can be compared to an
inundation survey which maps the maximum run-up along an entire coastline
rather than at a series of discrete sites. The survey map is
shown in Figure~\ref{fig:patongescapemap} and plots the maximum run-up
of the 2004 tsunami in Patong Bay. Refer to Szczucinski et
al~\cite{szczucinski06} for further details.

\subsubsection{Eyewitness Accounts}\label{sec:eyewitness data}
Eyewitness accounts detailed in~\cite{papadopoulos06}
report that most people at Patong Beach observed an initial retreat of
the shoreline of more than 100 m followed a few minutes later, by a
strong wave (crest). Another less powerful wave arrived another five
or ten minutes later. Eyewitness statements place the arrival time of
the strong wave between 2 hours and 55 minutes to 3 hours and 5
minutes after the source rupture (09:55am to 10:05am local time).

Two videos were sourced from the internet\footnote{The footage is 
widely available and can for example be obtained from 
\url{http://www.archive.org/download/patong_bavarian/patong_bavaria.wmv}
(Comfort Hotel) and
\url{http://www.archive.org/download/tsunami_patong_beach/tsunami_patong_beach.wmv}
(Novotel)}
%http://wizbangblog.com/content/2005/01/01/wizbang-tsunami.php
which include footage of the tsunami in Patong Bay on the day
of the Indian Ocean Tsunami. Both videos show an already inundated
group of buildings. They also show what is to be assumed as the second
and third waves approaching and further flooding of the buildings and
street.  The first video is in the very north, filmed from what is
believed to be the roof of the Novotel Hotel marked ``north'' in Figure
\ref{fig:gauge_locations}. The second video is in the very south,
filmed from the second story of a building next door to the Comfort
Resort near the corner of Ruam Chai St and Thaweewong Road.  This
location is marked ``south'' in Figure \ref{fig:gauge_locations}. 
Figure~\ref{fig:video_flow} shows stills from this video. Both videos
were used to estimate flow speeds and inundation depths over time.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=5.0cm,keepaspectratio=true]{flow_rate_south_0_00sec.jpg}
\includegraphics[width=5.0cm,keepaspectratio=true]{flow_rate_south_5_04sec.jpg}
\includegraphics[width=5.0cm,keepaspectratio=true]{flow_rate_south_7_12sec.jpg}
\includegraphics[width=5.0cm,keepaspectratio=true]{flow_rate_south_7_60sec.jpg}
\caption{Four frames from a video where flow rate could be estimated,
  circle indicates tracked debris, from top left: 0.0 sec, 5.0 s, 7.1
  s, 7.6 s.}
\label{fig:video_flow}
\end{center}
\end{figure}

Flow rates were estimated using landmarks found in both videos and
were found to be in the range of 5 to 7 metres per second (+/- 2 m/s)
in the north and 0.5 to 2 metres per second (+/- 1 m/s) in the south.
Water depths could also
be estimated from the videos by the level at which water rose up the
sides of buildings such as shops. Our estimates are in the order of
1.5 to 2.0 metres (+/- 0.5 m). 
Fritz ~\cite{fritz06} performed a detailed
analysis of video frames taken around Banda Aceh and arrived at flow
speeds in the range of 2 to 5 m/s.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=8.0cm,keepaspectratio=true]{patongescapemap.jpg}
\caption{Tsunami survey mapping the maximum observed inundation at
  Patong beach courtesy of the Thai Department of Mineral Resources
  \protect \cite{szczucinski06}.}
\label{fig:patongescapemap}
\end{center}
\end{figure}

\subsection{Validation Check-List}
\label{sec:checkList}
The data described in this section can be used to construct a
benchmark to validate all three stages of the evolution of a
tsunami. In particular we propose that a legitimate tsunami model
should reproduce the following behaviour:
\begin{itemize}
 \item reproduce the vertical deformation observed in north-western
   Sumatra and along the Nicobar--Andaman islands (see
   Section~\ref{sec:gen_data}).
 \item reproduce the \textsc{jason} satellite altimetry sea surface
   anomalies (see Section~\ref{sec:data_jason}).
 \item reproduce the inundation survey map in Patong bay
   (Figure~\ref{fig:patongescapemap}).
 \item simulate a leading depression followed by two distinct crests
   of decreasing magnitude.
 \item predict the water depths and flow speeds, at the locations of
   the eye-witness videos, that fall within the bounds obtained from
   the videos.
\end{itemize}

Ideally, the model should also be compared to measured timeseries of
waveheights and velocities but the authors are not aware of the
availability of such data.



%================Section===========================
\section{Modelling the Event}\label{sec:models}
Numerous models are currently used to model and predict tsunami
generation, propagation and run-up~\cite{titov97a,satake95}. Here we
introduce the three part modelling methodology employed by Geoscience Australia
to illustrate the utility of the proposed benchmark.

\subsection{Generation}\label{sec:modelGeneration}

There are various approaches to modelling the expected crustal
deformation from an earthquake at depth. Most approaches model the
earthquake as a dislocation in a linear elastic medium. Here we use
the method of Wang et al~\cite{wang03}. One of the main advantages
of their method is that it allows the dislocation to be located in a
stratified linear elastic half-space with an arbitrary number of
layers. Other methods (such as those based on Okada's equations) can
only model the dislocation in a homogeneous elastic half space, or can
only include a limited number of layers, and thus cannot model the
effect of the depth dependence of the elasticity of the
Earth~\cite{wang03}. The original versions of the codes described here
are available from \url{http://www.iamg.org/CGEditor/index.htm}. The
first program, \textsc{edgrn}, calculates elastic Green's function for
a set of point sources at a regular set of depths out to a specified
distance. The equations controlling the deformation are solved by
using a combination of Hankel's transform and Wang et al's
implementation of the Thomson-Haskell propagator
algorithm~\cite{wang03}. Once the Green's functions are calculated we
use a slightly modified version of \textsc{edcmp} to calculate the sea
floor deformation for a specific subfault. This second code
discretises the subfault into a set of unit sources and sums the
elastic Green's functions calculated from \textsc{edgrn} for all the
unit sources on the fault plane in order to calculate the final static
deformation caused by a two dimensional dislocation along the
subfault. This step is possible because of the linearity of the
governing equations. For this study, we have made minor modifications
to \textsc{edcmp} in order for it to provide output in a file format
compatible with the propagation code in the following section. Otherwise it
is similar to the original code.

In order to calculate the crustal deformation using these codes we
need a model that describes the variation in elastic
properties with depth and a slip model of the earthquake to describe
the dislocation. The elastic parameters used for this study are the
same as those in Table 2 of Burbidge~\cite{burbidge08}. For the slip
model, there are many possible models for the 2004 Andaman--Sumatran
earthquake to select from
~\cite{chlieh07,asavanant08,arcas06,grilli07,ioualalen07}. Some are
determined from various geological surveys of the site. Others solve
an inverse problem which calibrates the source based upon the tsunami
wave signal, the seismic signal and/or the run-up. The source
parameters used here to simulate the 2004 Indian Ocean tsunami were
taken from the slip model G-M9.15 of Chlieh
et al~\cite{chlieh07}. This model was created by inversion of wide
range of geodetic and seismic data. The slip model consists of 686
20km x 20km subsegments each with a different slip, strike and dip
angle. The dip subfaults go from $17.5^0$ in the north and $12^0$ in
the south. Refer to Chlieh et al~\cite{chlieh07} for a detailed
discussion of this model and its derivation. Note that the geodetic
data used in the validation was also included by~\cite{chlieh07} in
the inversion used to find G-M9.15. Thus the validation is not
completely independent. However, a reasonable validation would still
show that the crustal deformation and elastic properties model used
here is at least as valid as the one used by Chlieh
et al~\cite{chlieh07} and can reproduce the observations just as
accurately.

\subsection{Propagation}\label{sec:modelPropagation}
We use the \textsc{ursga} model described below to simulate the
propagation of the 2004 tsunami in the deep ocean ocean, based on a
discrete representation of the initial deformation of the sea floor, as
described in Section~\ref{sec:modelGeneration}. For the models shown
here, we assume that the uplift is instantaneous and creates a wave of
the same size and amplitude as the co-seismic sea floor deformation.

\subsubsection{URSGA}
\textsc{ursga} is a hydrodynamic code that models the propagation of
the tsunami in deep water using a finite difference method on a staggered grid.
It solves the depth integrated linear or nonlinear shallow water equations in
spherical co-ordinates with friction and Coriolis terms. The code is
based on Satake~\cite{satake95} with significant modifications made by
the \textsc{urs} corporation~\cite{thio08} and Geoscience
Australia~\cite{burbidge08}. The tsunami is propagated via the nested
grid system. Coarse grids are used in the open ocean and the finest
resolution grid is employed in the region of most
interest. \textsc{Ursga} is not publicly available.

\subsection{Inundation}\label{sec:modelInundation}
The utility of the \textsc{ursga} model decreases with water depth
unless an intricate sequence of nested grids is employed. In
comparison \textsc{anuga}, described below, is designed to produce
robust and accurate predictions of on-shore inundation, but is less
suitable for earthquake source modelling and large study areas because
it is based on projected spatial coordinates. Consequently, the
Geoscience Australia tsunami modelling methodology is based on a
hybrid approach using models like \textsc{ursga} for tsunami
propagation up to a 100 m depth contour.
%Specifically we use the \textsc{ursga} model to simulate the
%propagation of the 2004 Indian Ocean tsunami in the deep ocean, based
%on a discrete representation of the initial deformation of the sea
%floor, described in Section~\ref{sec:modelGeneration}.
The wave signal is then used as a time varying boundary condition for
the \textsc{anuga} inundation simulation.
% A description of \textsc{anuga} is the following section.

\subsubsection{ANUGA}
\textsc{Anuga} is an Open Source hydrodynamic inundation tool that
solves the conserved form of the depth-integrated nonlinear shallow
water wave equations. The scheme used by \textsc{anuga}, first
presented by Zoppou and Roberts~\cite{zoppou99}, is a high-resolution
Godunov-type method that uses the rotational invariance property of
the shallow water equations to transform the two-dimensional problem
into local one-dimensional problems. These local Riemann problems are
then solved using the semi-discrete central-upwind scheme of Kurganov
et al~\cite{kurganov01} for solving one-dimensional conservation
equations. The numerical scheme is presented in detail in 
Roberts and Zoppou~\cite{zoppou99,roberts00} and
Nielsen et al~\cite{nielsen05}. An important capability of the
software is that it can model the process of wetting and drying as
water enters and leaves an area. This means that it is suitable for
simulating water flow onto a beach or dry land and around structures
such as buildings. It is also capable of adequately resolving
hydraulic jumps due to the ability of the finite-volume method to
handle discontinuities. The numerical scheme can also handle
transitions between sub-critical and super-critical flow regimes
seamlessly. \textsc{Anuga} has been validated against a number of
analytical solutions and the wave tank simulation of the 1993 Okushiri
Island tsunami~\cite{nielsen05,roberts06}.

%================Section===========================
\section{Results}\label{sec:results}
This section presents a validation of the modelling practice of Geoscience
Australia against the new proposed benchmarks. The criteria outlined 
in Section~\ref{sec:checkList} are addressed for each of the three stages
of tsunami evolution.

\subsection{Generation}\label{modelGeneration}
The location and magnitude of the sea floor displacement associated
with the 2004 Sumatra--Andaman tsunami calculated from the G-M9.15
model of~\cite{chlieh07} is shown in
Figure~\ref{fig:surface_deformation}. The magnitude of the sea floor
displacement ranges from about $-3.0$ to $5.0$ metres. The region near
the fault is predicted to uplift, while that further away from the
fault subsides. Also shown in Figure~\ref{fig:surface_deformation} are
the areas that were observed to uplift (arrows pointing up) or subside
(arrows point down) during and immediately after the earthquake. Most
of this data comes from uplifted or subsided coral heads. The length of
vector increases with the magnitude of the displacement; the length
corresponding to 1 m of observed motion is shown in the top right
corner of the figure. As can be seen, the source model detailed in
Section~\ref{sec:modelGeneration} produces a crustal deformation that
matches the vertical displacements in the Nicobar-Andaman islands and
Sumatra very well. Uplifted regions are close to the fault and
subsided regions are further away. The crosses on
Figure~\ref{fig:surface_deformation} show estimates of the pivot line
from the remote sensing data~\cite{chlieh07} and they follow the
predicted pivot line quite accurately. The average difference between
the observed motion and the predicted motion (including the pivot line
points) is only 0.06 m, well below the typical error of the
observations of between 0.25 and 1.0 m. However, the occasional point
has quite a large error (over 1 m); for example a couple
uplifted/subsided points appear to be on a wrong side of the predicted
pivot line~\ref{fig:surface_deformation}. The excellence of the fit is
not surprising, since the original slip model was chosen
by~\cite{chlieh07} to fit this (and the seismic data) well. 
This does demonstrate, however, that \textsc{edgrn} and our modified version of
\textsc{edstat} can reproduce the correct pattern of vertical
deformation very well when the slip distribution is well constrained
and when reasonable values for the elastic properties are used.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=5cm,keepaspectratio=true]{surface_deformation.jpg}
%\includegraphics[totalheight=0.3\textheight,width=0.8\textwidth]{surface_deformation.jpg}
\caption{Location and magnitude of the vertical component of the sea
  floor displacement associated with the 2004 Indian Ocean tsunami
  based on the slip model, G-M9.15. The black arrows which point up
  show areas observed to uplift during and immediately after the
  earthquake; those pointing down are locations which subsided. The
  length of the arrow increases with the magnitude of the deformation. The arrow
  length corresponding to 1 m of deformation is shown in the top right
  hand corner of the figure. The cross marks show the location of
  the pivot line (the region between the uplift and subsided region
  where the uplift is zero) derived from remote sensing. All the
  observational data are from the dataset collated
  by~\cite{chlieh07}.}
\label{fig:surface_deformation}
\end{center}
\end{figure}

\subsection{Propagation}\label{sec:resultsPropagation}
The deformation results described in Section~\ref{sec:modelGeneration}
were used to provide a profile of the initial ocean surface
displacement. This wave was used as an initial condition for
\textsc{ursga} and was propagated throughout the Bay of Bengal. The
rectangular computational domain of the largest grid extended from
90$^0$ to 100$^0$ East and 0 to 15$^0$ North and contained
1335$\times$1996 finite difference points. Inside this grid, a nested
sequence of grids was used. The grid resolution of the nested grids
went from 27 arc seconds in the coarsest grid, down to nine arc seconds
in the second grid, three arc seconds in the third grid and finally one arc
second in the finest grid near Patong. The computational domain is
shown in Figure~\ref{fig:computational_domain}.

Figure \ref{fig:jasonComparison} provides a comparison of the
\textsc{ursga}-predicted sea surface elevation with the JASON
satellite altimetry data. The \textsc{ursga} model replicates the 
amplitude and timing of the the wave observed at $2.5^0$ South,
but underestimates the amplitude of the wave further to the south at
$4^0$ South. In the model, the southern most of these two waves
appears only as a small bump in the cross section of the model (shown
in Figure~\ref{fig:jasonComparison}) instead of being a distinct peak
as can be seen in the satellite data. Also note
that the \textsc{ursga} model prediction of the ocean surface
elevation becomes out of phase with the JASON data at $3^0$ to $7^0$ North
latitude. Chlieh et al~\cite{chlieh07} also observed these misfits and
suggest it is caused by a reflected wave from the Aceh Peninsula that
is not resolved in the model due to insufficient resolution of the
computational mesh and bathymetry data. This is also a limitation of
the model presented here, but probably could be improved by nesting
grids near Aceh.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=12.0cm,keepaspectratio=true]{jasonComparison.jpg}
\caption{Comparison of the \textsc{ursga}-predicted surface elevation
  with the JASON satellite altimetry data. The \textsc{ursga} wave
  heights have been corrected for the time the satellite passed
  overhead compared to JASON sea level anomaly.}
\label{fig:jasonComparison}
\end{center}
\end{figure}

\subsection{Inundation}
After propagating the tsunami in the open ocean using \textsc{ursga},
the approximated ocean and surface elevation and horisontal flow
velocities were extracted and used to construct a boundary condition
for the \textsc{anuga} model. The interface between the \textsc{ursga}
and \textsc{anuga} models was chosen to roughly follow the 100 m depth
contour along the west coast of Phuket Island. The computational
domain is shown in Figure~\ref{fig:computational_domain}.
\begin{figure}[ht]
\begin{center}
%\includegraphics[width=5.0cm,keepaspectratio=true]{extent_of_ursga_model.jpg}
%\includegraphics[width=5.0cm,keepaspectratio=true]{ursgaDomain.jpg}
\includegraphics[width=5.0cm,keepaspectratio=true]{extent_of_ANUGA_model.jpg}
\caption{Computational domain of the URSGA simulation (inset: white and black squares and main: black square) and the \textsc{anuga} simulation (main and inset: red polygon).}
\label{fig:computational_domain}
\end{center}
\end{figure}

The domain was discretised into 386,338 triangles. The resolution of
the grid was increased in certain regions to efficiently increase the
accuracy of the simulation. The grid resolution ranged between a
maximum triangle area of $1\times 10^5$ m$^2$ near the Western ocean
boundary to $20$ m$^2$ in the small regions surrounding the inundation
region in Patong Bay. Due to a lack of available data, friction was
set to a constant throughout the computational domain. For the
reference simulation a Manning's coefficient of 0.01 was chosen to
represent a small resistance to the water flow. See Section
\ref{sec:friction sensitivity} for details on model sensitivity to
this parameter.


The boundary condition at each side of the domain towards the south
and the north where no data was available was chosen as a transmissive
boundary condition, effectively replicating the time dependent wave
height present just inside the computational domain. Momentum was set
to zero. Other choices include applying the mean tide value as a
Dirichlet type boundary condition. But experiments as well as the
result of the verification reported here showed that this approach
tends to underestimate the tsunami impact due to the tempering of the
wave near the side boundaries, whereas the transmissive boundary
condition robustly preserves the wave.

During the \textsc{anuga} simulation the tide was kept constant at
$0.80$ m. This value was chosen to correspond to the tidal height
specified by the Thai Navy tide charts
(\url{http://www.navy.mi.th/hydro/}) at the time the tsunami arrived
at Patong Bay. Although the tsunami propagated for approximately 3
hours before it reach Patong Bay, the period of time during which the
wave propagated through the \textsc{anuga} domain is much
smaller. Consequently the assumption of constant tide height is
reasonable.

Maximum onshore inundation elevation was computed from the model
throughout the entire Patong Bay region. 
Figure~\ref{fig:inundationcomparison1cm} (left) shows very good
agreement between the measured and simulated inundation. However
these results are dependent on the classification used to determine
whether a region in the numerical simulation was inundated. In
Figure~\ref{fig:inundationcomparison1cm} (left) a point in the computational
domain was deemed inundated if at some point in time it was covered by
at least 1 cm of water. However, the precision of the inundation boundary 
generated by the on-site survey is most likely less than that as it 
was determined by observing water marks and other signs
left by the receding waters. Consequently the measurement error along
the inundation boundary of the survey is likely to vary significantly 
and somewhat unpredictably. 
An inundation threshold of 10 cm therefore was selected for inundation 
extents reported in this paper to reflect
the more likely accuracy of the survey, and subsequently facilitate a more
appropriate comparison between the modelled and observed inundation
area.
Figure~\ref{fig:inundationcomparison1cm} (right) shows the simulated
inundation using a larger threshold of 10 cm. 


An animation of this simulation is available on the \textsc{anuga} website at \url{https://datamining.anu.edu.au/anuga} or directly from \url{http://tinyurl.com/patong2004}.

%\url{https://datamining.anu.edu.au/anuga/attachment/wiki/AnugaPublications/patong_2004_indian_ocean_tsunami_ANUGA_animation.mov}.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=6.0cm,keepaspectratio=true]{final_1cm.jpg}
\includegraphics[width=6.0cm,keepaspectratio=true]{final_10cm.jpg}
\caption{Simulated inundation versus observed inundation using an 
inundation threshold of 1cm (left) and 10cm (right).}
\label{fig:inundationcomparison1cm}
\end{center}
\end{figure}

To quantify the agreement between the observed and simulated inundation we
introduce the measure
\begin{equation}
\rho_{in}=\frac{A(I_m\cap I_o)}{A(I_o)}
\end{equation}
representing the ratio $\rho_{in}$ of the observed
inundation region $I_o$ captured by the model $I_m$. Another useful
measure is the fraction of the modelled inundation area that falls
outside the observed inundation area given by the formula
\begin{equation}
\rho_{out}=\frac{A(I_m\setminus (I_m\cap I_o))}{A(I_o)}
\end{equation}
These values for the two aforementioned simulations are given in
Table~\ref{table:inundationAreas}. High value of both $\rho_{in}$ and $\rho_{out}$ indicate 
that the model overestimates the impact whereas low values of both quantities would indicate 
an underestimation. A high value of $\rho_{in}$ combined with a low value of $\rho_{out}$ 
indicates a good model prediction of the survey. 

Discrepancies between the survey data and the modelled inundation
include: unknown distribution of surface roughness, inappropriate
parameterisation of the source model, effect of humans structures on
flow, as well as uncertainties in the elevation data, effects of
erosion and deposition by the tsunami event, measurement errors, and
missing data in the field survey data itself. The impact of some of
these sources of uncertainties are is investigated in
Section~\ref{sec:sensitivity}

\subsection{Eye-witness accounts}
Figure \ref{fig:gauge_locations} shows four locations where time
series have been extracted from the model. The two offshore time series
are shown in Figure \ref{fig:offshore_timeseries} and the two onshore
timeseries are shown in Figure \ref{fig:onshore_timeseries}. The
latter coincide with locations where video footage from the event is
available as described in Section \ref{sec:eyewitness data}.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=8.0cm,keepaspectratio=true]{gauge_locations.jpg}
\caption{Location of timeseries extracted from the model output.}
\label{fig:gauge_locations}
\end{center}
\end{figure}


\begin{figure}[ht]
\begin{center}
\includegraphics[width=10.0cm,keepaspectratio=true]{gauge_bay_depth.jpg}
\includegraphics[width=10.0cm,keepaspectratio=true]{gauge_bay_speed.jpg}
\caption{Time series obtained from the two offshore locations shown in Figure \protect \ref{fig:gauge_locations}.}
\end{center}
\label{fig:offshore_timeseries}
\end{figure}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=10.0cm,keepaspectratio=true]{gauges_hotels_depths.jpg}
\includegraphics[width=10.0cm,keepaspectratio=true]{gauges_hotels_speed.jpg}
\caption{Time series obtained from the two onshore locations shown in Figure \protect \ref{fig:gauge_locations}.}
\end{center}
\label{fig:onshore_timeseries}
\end{figure}


The estimated max depths and flow rates given in Section
\ref{sec:eyewitness data} are shown together with the modelled depths
and flow rates obtained from the model in Table \ref{tab:depth and
  flow comparisons}. The minimum depths shown in the model are clearly
lower than expected and an indication that the tsunami model does not
predict flow dynamics accurately at this level of detail. However,
this comparison serves to check that depths and speeds predicted are
within the range of what is expected.


\begin{table}
\[
  \begin{array}{|l|cc|cc|} 
  \hline
                 & \multicolumn{2}{|c|}{\mbox{Depth [m]}} 
                 & \multicolumn{2}{c|}{\mbox{Flow [m/s]}} \\ 
                 & \mbox{Observed} & \mbox{Modelled}
                 & \mbox{Observed} & \mbox{Modelled} \\ \cline{2-5}                 
    \mbox{North} & 1.5-2 & 1.4 & 5-7 & 0.1 - 3.3 \\
    \mbox{South} & 1.5-2 & 1.5 & 0.5-2 & 0.2 - 2.6 \\ \hline
  \end{array} 
\]
\label{tab:depth and flow comparisons}
\end{table} 

%can be estimated with landmarks found in
%satellite imagery and the use of a GIS and were found to be in the
%range of 5 to 7 metres per second (+/- 2 m/s) in the north and 0.5 to
%2 metres per second (+/- 1 m/s) in the south. 

Given the uncertainties in both model and observations, there is agreement
between the values obtained from the videos and the simulations.

% Our modelled flow rates show
%maximum values in the order of 0.2 to 2.6 m/s in the south and 0.1 to
%3.3 m/s for the north as shown in the figures. Water depths could also
%be estimated from the videos by the level at which water rose up the
%sides of buildings such as shops. Our estimates are in the order of
%1.5 to 2.0 metres (+/- 0.5 m). This is in the same range as our
%modelled maximum depths of 1.4 m in the north and 1.5 m in the south
%as seen in the figure. 





%================Section===========================
\section{Sensitivity Analysis}
\label{sec:sensitivity}
This section investigates the effect of different values of Manning's
friction coefficient, changing waveheight at the 100 m depth contour,
and the presence and absence of buildings in the elevation dataset on
model maximum inundation. The reference model is the one reported in 
Figure~\ref{fig:inundationcomparison1cm} (right) with a friction coefficient of 0.01, 
buildings included and the boundary condition produced by the URSGA model.

%========================Friction==========================%
\subsection{Friction}
\label{sec:friction sensitivity}
The first study investigated the impact of surface roughness on the
predicted run-up. According to Schoettle~\cite{schoettle2007}
appropriate values of Manning's coefficient range from 0.007 to 0.03
for tsunami propagation over a sandy sea floor and the reference model
uses a value of 0.01.  To investigate sensitivity to this parameter,
we simulated the maximum onshore inundation using a Manning's
coefficient of 0.0003 and 0.03. The resulting inundation maps are
shown in Figure~\ref{fig:sensitivity_friction} and the maximum flow
speeds in Figure~\ref{fig:sensitivity_friction_speed}. These figures
show that the on-shore inundation extent decreases with increasing
friction and that small perturbations in the friction cause bounded
changes in the output. This is consistent with the conclusions of
Synolakis~\cite{synolakis05} et al, who state that the long wavelength of
tsunami tends to mean that friction is less important in
comparison to the motion of the wave.

%========================Wave-Height==========================%
\subsection{Input Wave Height}\label{sec:waveheightSA}
The effect of the wave height used as input to the inundation model
\textsc{anuga} was also investigated.
Figure~\ref{fig:sensitivity_boundary} indicates that the inundation
severity is directly proportional to the boundary waveheight but small
perturbations in the input wave height of 10 cm appear to have little
effect on the final on-shore run-up. Obviously larger perturbations
will have greater impact. However, this value is generally well
predicted by the generation and propagation models such as
\textsc{ursga}. See e.g.\ \cite{thomas2009}.



%========================Buildings==========================%
\subsection{Buildings and Other Structures}
The presence of buildings has the greatest influence on the maximum
on-shore inundation extent. Figure~\ref{fig:sensitivity_nobuildings}
shows the maximum run-up and associated flow speeds in the presence and absence of buildings. It
is apparent that the inundation is much more severe when the presence
of human made structures and buildings are ignored.

\begin{table}
\begin{center}
\label{table:inundationAreas}
\caption{$\rho_{in}$ and $\rho_{out}$ of the reference simulation and all sensitivity studies.}
\begin{tabular}{|l|c|c|}
\hline
 & $\rho_{in}$ & $\rho_{out}$ \\ 
\hline\hline
Reference model & 0.79 & 0.20\\ 
Friction = 0.0003 & 0.83 & 0.26 \\ 
Friction = 0.03 & 0.67 & 0.09\\ 
Boundary wave hight minus 10 cm & 0.77 & 0.17 \\
Boundary wave hight plus 10 cm & 0.82 & 0.22 \\
No Buildings & 0.94 & 0.44 \\
\hline 
\end{tabular}
\end{center}
\end{table}

%================Section===========================

\section{Conclusion}
This paper proposes an additional field data benchmark for the
verification of tsunami inundation models. Currently, there is a
scarcity of appropriate validation datasets due to a lack of well-documented
historical tsunami impacts. The benchmark proposed here
utilises the uniquely large amount of observational data for model
comparison obtained during, and immediately following, the
Sumatra--Andaman tsunami of 26 December 2004. Unlike the small
number of existing benchmarks, the proposed test validates all three
stages of tsunami evolution - generation, propagation and
inundation. In an attempt to provide higher visibility and easier
accessibility for tsunami benchmark problems, the data used to
construct the proposed benchmark is documented and freely available at
\url{http://tinyurl.com/patong2004-data}.

This study also shows that the tsunami impact modelling methodology
adopted is sane and able to predict inundation extents with reasonable
accuracy.  An associated aim of this paper was to further validate the
hydrodynamic modelling tool \textsc{anuga} which is used to simulate
the tsunami inundation and run rain-induced floods. Model predictions
matched well geodetic measurements of the Sumatra--Andaman earthquake,
altimetry data from the \textsc{jason}, eye-witness accounts of wave
front arrival times and flow speeds and a detailed inundation survey
of Patong Bay, Thailand.

A simple sensitivity analysis was performed to assess the influence of
small changes in friction, wave height at the 100 m depth contour and
the presence of buildings and other structures on the model
predictions. The presence of buildings has the greatest influence on
the simulated inundation extent. The value of friction and small
perturbations in the waveheight at the \textsc{anuga} boundary have
comparatively little effect on the model results.

%================Acknowledgement===================
\section*{Acknowledgements}
This project was undertaken at Geoscience Australia and the Department
of Mathematics, The Australian National University. The authors would
like to thank Niran Chaimanee from the CCOP, Thailand for providing
the post 2004 tsunami survey data, building footprints, aerial
photography and the elevation data for Patong beach, Prapasri Asawakun
from the Suranaree University of Technology and Parida Kuneepong for
supporting this work; and Drew Whitehouse from the Australian National
University for preparing the animation of the inundation model.

\clearpage
\section{Appendix}

This appendix present the images used to assess the model sensitivities described in 
Section~\ref{sec:sensitivity}.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=6cm,keepaspectratio=true]{sensitivity_reference_depth}
\includegraphics[width=6cm,keepaspectratio=true]{sensitivity_reference_speed}
\caption{Results from reference model as reported in Section \protect \ref{sec:results}, 
  i.e.\ including buildings and a friction value of 0.01. The seaward boundary condition is as
  provided by the URSGA model. The left image shows the maximum
  modelled depth while the right hand image shows the maximum modelled
  flow velocities.}
\label{fig:reference_model}
\end{center}
\end{figure}



\begin{figure}[ht]
\begin{center}
\includegraphics[width=6cm,keepaspectratio=true]{sensitivity_minus10cm_depth}
\includegraphics[width=6cm,keepaspectratio=true]{sensitivity_plus10cm_depth}
\caption{Model results with wave height at \textsc{anuga} boundary artificially
  modified to assess sensitivities. The reference inundation extent is shown in Figure
  \protect \ref{fig:reference_model} (left).  The left and right images 
  show the inundation results if the wave at the \textsc{anuga} boundary
  is reduced or increased by 10cm respectively. The inundation
  severity varies in proportion to the boundary waveheight, but the
  model results are only slightly sensitive to this parameter for the
  range of values tested.}
\label{fig:sensitivity_boundary}
\end{center}
\end{figure}


\begin{figure}[ht]
\begin{center}
\includegraphics[width=6cm,keepaspectratio=true]{sensitivity_minus10cm_speed}
\includegraphics[width=6cm,keepaspectratio=true]{sensitivity_plus10cm_speed}
\caption{The maximal flow speeds for the same model parameterisations
  found in Figure \protect \ref{fig:sensitivity_boundary}. The
  reference flow speeds are shown in Figure \protect
  \ref{fig:reference_model} (right).}
\label{fig:sensitivity_boundary_speed}
\end{center}
\end{figure}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=6cm,keepaspectratio=true]{sensitivity_nobuildings_depth}
\includegraphics[width=6cm,keepaspectratio=true]{sensitivity_nobuildings_speed}
\caption{This figure shows the effect of having buildings as part of
  the elevation data set. 
  The left hand image shows the inundation depth results for
  a model entirely without buildings.  As expected, the absence of
  buildings will increase the inundation extent beyond what was
  surveyed. The right hand image shows the corresponding flow speeds in the absence of buildings.  
  The reference results are as shown in Figure
  \protect \ref{fig:reference_model}.}
\label{fig:sensitivity_nobuildings}
\end{center}
\end{figure}


\begin{figure}[ht]
\begin{center}
\includegraphics[width=6cm,keepaspectratio=true]{sensitivity_f0_0003_depth}
\includegraphics[width=6cm,keepaspectratio=true]{sensitivity_f0_03_depth}
\caption{Model results for different values of Manning's friction
  coefficient shown to assess sensitivities. The reference inundation extent for a 
  friction value of 0.01 is shown in Figure
  \protect \ref{fig:reference_model} (left).  The left and right images 
  show the inundation results for friction values of 0.0003 and
  0.03 respectively. The inundation extent increases for the lower
  friction value while the higher slows the flow and decreases the
  inundation extent. Ideally, friction should vary across the entire
  domain depending on terrain and vegetation, but this is beyond the
  scope of this study.}
\label{fig:sensitivity_friction}
\end{center}
\end{figure}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=6cm,keepaspectratio=true]{sensitivity_f0_0003_speed}
\includegraphics[width=6cm,keepaspectratio=true]{sensitivity_f0_03_speed}
\caption{The maximal flow speeds for the same model parameterisations
  found in Figure \protect \ref{fig:sensitivity_friction}. The
  reference flow speeds are shown in Figure \protect
  \ref{fig:reference_model} (right).}
\label{fig:sensitivity_friction_speed}
\end{center}
\end{figure}

\clearpage

%====================Bibliography==================
\bibliographystyle{spmpsci}
\bibliography{tsunami07}
\end{document}