source: anuga_work/publications/boxing_day_validation_2008/results.tex @ 7451

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

\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. 
%
%However, the occasional point
%has quite a large error (over 1 m); for example a couple of
%uplifted/subsided points appear to be on a wrong 
%(FIXME (Jane): This is incorrect) 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 the motion and seismic data well. 
Consequently the replication of the generation data has little 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.
%
%This does demonstrate, however, that \textsc{edgrn} and our modified version of
%\textsc{edstat} (FIXME(Jane): This has never been mentioned before) 
%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=0.8\textwidth,keepaspectratio=true]{figures/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 
  (FIXME(Jane): How was that possible?). 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}.

FIXME (Ole): I know that a nested ursga model was trialled for the 
end-to-end modelling. However, for the study done here, where models
were coupled, I didn't think nested grids were used with URSGA - 
and certainly not down to 1 arc second. Can someone shed some light 
on this please? RICHARD

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.8\textwidth,keepaspectratio=true]{figures/extent_of_ANUGA_model.jpg}
\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{center}
\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^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 \textsc{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 which could be improved by nesting
grids near Aceh.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=\textwidth,keepaspectratio=true]{figures/jasonComparison.jpg}
\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{center}
\end{figure}
FIXME (Jane): This graph does not look nice.

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

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 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 be 
commensurate with the model output from the \textsc{ursga} model
while the latter was chosen to match the available resolution of topographic 
data and building data in Patong city.
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 
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 set
to zero. Other choices include applying the mean tide value as a
Dirichlet 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 three
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.

\subsection{Inundation}
The \textsc{anuga} simulation described in the previous section and used to
 model shallow water propgation also predicts
inundation. Maximum onshore inundation depth was computed from the model
throughout the entire Patong Bay region 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. 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. 


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}.
The scripts required are part of the \textsc{anuga} distribution also 
available from Sourceforge \url{http://sourceforge.net/projects/anuga} under 
the validation section.

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}
\includegraphics[width=\textwidth,keepaspectratio=true]{figures/threshold.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 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}. 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 in the GPS survey recordings, 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=0.8\textwidth,keepaspectratio=true]{figures/gauge_bay_depth.jpg}
\includegraphics[width=0.8\textwidth,keepaspectratio=true]{figures/gauge_bay_speed.jpg}
\caption{Time series obtained from the two offshore gauge locations, 
7C and 10C, shown in Figure \protect \ref{fig:gauge_locations}.}
\end{center}
\label{fig:offshore_timeseries}
\end{figure}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=\textwidth,keepaspectratio=true]{figures/gauges_hotels_depths.jpg}
\includegraphics[width=\textwidth,keepaspectratio=true]{figures/gauges_hotels_speed.jpg}
\caption{Time series obtained from the two onshore locations, North and South, 
shown in Figure \protect \ref{fig:gauge_locations}.}
\end{center}
\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. 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 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} 
FIXME (Jane): We should perhaps look at average data in area surrounding these points

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