\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.

\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
%the 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. 

The excellence of the fit is
not surprising, since the original slip model was chosen
by~\cite{chlieh07} to fit the motion and seismic data well. 
Consequently, the replication of the generation data has limited meaning for 
the model structure presented in Section~\ref{sec:models}. But for 
uncalibrated source models or source models calibrated on other data 
this test of generation would be more meaningful.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth,keepaspectratio=true]{figures/surface_deformation.jpg}
\end{center}

\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 compared with observed deformation 
  (arrows). 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 of where there
  were no changes in the amount of reef or land exposure before and after
the earthquake (\cite{meltzner06}). All the
  observational data are from the dataset collated
  by~\cite{chlieh07}.}
\label{fig:surface_deformation}
\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$^\circ$ to 100$^\circ$ East and $0^\circ$ to 15$^\circ$ 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 and three arc seconds in the third grid. The computational domain is
shown in Figure~\ref{fig:computational_domain}.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.8\textwidth,keepaspectratio=true]{figures/extent_of_ANUGA_model.jpg}
\end{center}

\caption{Computational domain of the \textsc{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{figure}


Figure \ref{fig:jasonComparison} provides a comparison of the
\textsc{ursga}-predicted sea surface elevation with the \textsc{jason}
satellite altimetry data. The \textsc{ursga} model replicates the 
amplitude and timing of the the wave observed at $2.5^\circ$ South,
but underestimates the amplitude of the wave further to the south at
$4^\circ$ 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 \textsc{jason} 
data at $3^\circ$ to $7^\circ$ 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 which could be improved by nesting
grids near Aceh.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=\textwidth,keepaspectratio=true]{figures/jasonComparison.jpg}
\end{center}

\caption{Comparison of the \textsc{ursga}-predicted surface elevation
  with the \textsc{jason} satellite altimetry data. The \textsc{ursga} wave
  heights have been corrected for the time the satellite passed
  overhead compared to \textsc{jason} sea level anomaly.}
\label{fig:jasonComparison}
\end{figure}

After propagating the tsunami in the open ocean using \textsc{ursga},
the approximated ocean and surface elevation and horizontal flow
speeds 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. Data from the
three~second grid which is approximately 30~m apart was decimated to
match the resolution chosen in \textsc{Anuga}.  The computational
domain is shown in Figure~\ref{fig:computational_domain}.

The domain was discretised into 386,338 triangles. The resolution of
the grid was increased in regions inside the bay and on-shore to 
efficiently increase the simulation accuracy for the impact area. 
The grid resolution ranged between a
maximum triangle area of $1\times 10^5$~m$^2$ 
(corresponding to approximately 440~m between mesh points) 
near the western ocean
boundary (roughly following the 100~m depth contour) 
to $20$~m$^2$ (corresponding to 
approximately 6 m between mesh points) 
in the small regions surrounding the inundation
region in Patong Bay. The coarse resolution was chosen to balance accuracy 
with computational costs while the fine resolution was chosen to match the available resolution of topographic 
data and building data in Patong City. Figure~\ref{fig:mesh} shows a 
section of the mesh covering the southern part of the City.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.8\textwidth, keepaspectratio=true]{figures/mesh_section.jpg}
\end{center}

\caption{Section of the mesh used by ANUGA to simulate the tsunami inundation.
The finest mesh resolution is approximately 6 m between nodes which 
is sufficient to resolve individual buildings affecting the flows.}
\label{fig:mesh}
\end{figure}


Due to a lack of available roughness 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 information about the incident wave or 
its velocity field is available from the \textsc{Ursga} model 
was chosen as a transmissive
boundary condition, effectively replicating the time dependent wave
height present just inside the computational domain. 
The velocity field on these boundaries was kept at
to zero during the simulation. Other choices include applying the mean tide value as a
Dirichlet boundary condition. 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 in the offshore region 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 three
hours before it reached Patong Bay, the period of time during which the
wave propagated through the \textsc{anuga} domain is much
smaller of the order of 2 hours. Consequently the assumption of constant tide height is
reasonable. The initial water level for the river was set to 0.

\subsection{Inundation}\label{sec:inundation results}
The \textsc{anuga} simulation described in the previous section and used to
 model shallow water propagation also predicts
inundation. Maximum onshore inundation depth was computed from the inundation model
and used to generate a measure of the inundated area. 
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. The precision of the inundation boundary 
generated by the on-site survey is most likely less than this 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 and Figure~\ref{fig:anuga screenshot} shows a screenshot from the inundation model.
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=\textwidth,keepaspectratio=true]{figures/threshold.jpg}
\end{center}

\caption{Simulated inundation versus observed inundation using an 
inundation threshold of 1 cm (left) and 10 cm (right).}
\label{fig:inundationcomparison1cm}
\end{figure}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.8\textwidth,keepaspectratio=true]{figures/anuga_screenshot.jpg}
\end{center}

\caption{Screenshot of the \textsc{anuga} simulation of the inundation
 at Patong City at local time 10:53 about one hour after the arrival of 
 the first wave.}
\label{fig:anuga screenshot}
\end{figure}


\subsubsection{Comparison to survey}

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 of the area of the observed 
inundation region $I_o$ and the area of the observed inundation region 
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} along with results from the sensitivity analysis in 
Section~\ref{sec:sensitivity}. High values 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 arise from errors and uncertainties in both the field surveys and the models. 
The former include measurement errors in the GPS survey recordings and
missing data in the field survey data itself.
The latter include 
unknown distribution of surface roughness, uncertainties in the
parameterisation of the source model, discretisation errors, 
effect of humans structures on
flow, as well as uncertainties in the elevation data including effects of
erosion and deposition by the tsunami event. 
The impacts of some of the model uncertainties are as investigated in
Section~\ref{sec:sensitivity}.


\subsubsection{Reproducibility of inundation results}
As one aim of this paper is to provide a new benchmark for tsunami 
inundation modelling we have made the datasets available
available on \textsc{Sourceforge} in the \textsc{anuga}
project (\url{http://sourceforge.net/projects/anuga}) under the directory 
\url{validation\_data/patong-1.0}.
At the time of 
writing a 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}.

To reproduce the inundation modelling results using this data, the reader
will need to run the validation scripts (\url{anuga\_validation/automated\_validation\_tests/patong\_beach\_validation}) which are part of the 
\textsc{anuga} distribution also available from 
\textsc{Sourceforge} \url{http://sourceforge.net/projects/anuga}.


\subsection{Eye-witness accounts}

\subsubsection{Arrival time}
The arrival time of the first wave took place between 9:55 and 10:05 as described in
Section~\ref{sec:eyewitness data}. The modelled arrival time at the beach is around 10:02 
as can be verified from the animation provided in 
Section \ref{sec:inundation results} or from Figure~\ref{fig:onshore_timeseries} below. 
Subsequent waves of variable magnitude appear over the next two hours 
approximately 20-30 minutes apart.
% 10:02, 10:20, 10:47, 11:14, 11:44
The first arrival and overall dynamic behaviour is therefore reasonably consistent with the 
eye-witness accounts.

\subsubsection{Observed wave dynamics}
Figure \ref{fig:gauge_locations} shows two on-shore locations where time
series have been extracted from the model. These are the locations where video footage from the event is
available as described in Section \ref{sec:eyewitness data}.
The corresponding time series are shown in Figure \ref{fig:onshore_timeseries}.


%\begin{figure}[ht]
%\begin{center}
%\includegraphics[width=0.7\textwidth,keepaspectratio=true]{figures/gauge_bay_depth.jpg}
%\includegraphics[width=0.7\textwidth,keepaspectratio=true]{figures/gauge_bay_speed.jpg}
%\end{center}
%
%\caption{Time series obtained from the two offshore gauge locations, 
%7C and 10C, shown in Figure \protect \ref{fig:gauge_locations}.}
%\label{fig:offshore_timeseries}
%\end{figure}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth,keepaspectratio=true]{figures/gauges_hotels_depths.jpg}
\includegraphics[width=0.7\textwidth,keepaspectratio=true]{figures/gauges_hotels_speed.jpg}
\end{center}

\caption{Time series obtained from the two onshore locations, North and South, 
shown in Figure \protect \ref{fig:gauge_locations}. Time is given in hours since the earthquake event (7:59am).}
\label{fig:onshore_timeseries}
\end{figure}


The estimated 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 predicted maximum depths and speeds are all of the same order 
of what was observed as is the approximate arrival time at the two locations. However, unlike the real event, 
the model estimates complete withdrawal of the water between waves at the 
chosen locations and shows that the model must be used with caution at this
level of detail.
Nonetheless, this comparison serves to check that the peak 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 (peak)}
                 & \mbox{Observed} & \mbox{Modelled (peak)} \\ \cline{2-5}                 
    \mbox{North} & 1.5-2 & 1.4 & 5-7 & 3.3 \\
    \mbox{South} & 1.5-2 & 1.5 & 0.5-2 & 2.6 \\ \hline
  \end{array} 
\]
\caption{Observed depth and flows from the video footage compared to values extracted from the inundation model.}
\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. 





