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

Last change on this file since 7404 was 7404, checked in by ole, 15 years ago

Corrected reference to eye witness accounts

File size: 61.2 KB
Line 
1\documentclass[a4paper]{article}
2% to be added when submitted to ocean dynamics
3%\documentclass[smallcondensed,draft]{svjour3}
4%\usepackage{mathptmx}
5%\journalname{Ocean Dynamics}
6
7\usepackage{graphicx}
8\usepackage{hyperref}
9\usepackage{amsfonts}
10\usepackage{url}      % for URLs and DOIs
11\newcommand{\doi}[1]{\url{http://dx.doi.org/#1}}
12
13
14%----------title-------------%
15\title{Benchmarking Tsunami Models using the December 2004 Indian
16  Ocean Tsunami and its Impact at Patong Bay}
17
18%-------authors-----------
19\author{J.~D. Jakeman \and O. Nielsen \and K. VanPutten \and
20  D. Burbidge \and R. Mleczko \and N. Horspool}
21
22% to be added when submitted to ocean dynamics
23%\institute{J.~D. Jakeman \at
24%       The Australian National University, Canberra, \textsc{Australia}\
25%       \email{john.jakeman@anu.edu.au}
26%       \and
27%       O. Nielsen \and R. Mleczko \and D. Burbidge \and K. VanPutten \and N. Horspool \at
28%       Geoscience Australia, Canberra, \textsc{Australia}
29%}
30
31
32%================Start of Document================
33\begin{document}
34\maketitle
35%------Abstract--------------
36\begin{abstract}
37This paper proposes a new benchmark for tsunami model validation.
38The benchmark is based upon the 2004 Indian Ocean tsunami,
39which affords a uniquely large amount of observational data for events of this kind.
40Unlike the small number of existing benchmarks, the
41proposed test validates all three stages of tsunami evolution -
42generation (FIXME (Jane): really?), propagation and inundation. Specifically we use geodetic
43measurements of the Sumatra--Andaman earthquake to validate the
44tsunami source, altimetry data from the \textsc{jason} satellite to
45test open ocean propagation, eye-witness accounts to assess near shore
46propagation and a detailed inundation survey of Patong city, Thailand
47to compare model and observed inundation. Furthermore we utilise this
48benchmark to further validate the hydrodynamic modelling tool
49\textsc{anuga} which is used to simulate the tsunami
50inundation. Important buildings and other structures were incorporated
51into the underlying computational mesh and shown to have a large
52influence on inundation extent. Sensitivity analysis also showed that
53the model predictions are comparatively insensitive to large changes
54in friction and small perturbations in wave weight at the 100 m depth
55contour.
56% to be added when submitted to ocean dynamics
57%\keywords{Tsunami \and modelling \and validation and verification \and benchmark}
58\end{abstract}
59
60\tableofcontents
61%================Section===========================
62
63\section{Introduction}
64Tsunami is a potential hazard to coastal communities all over the
65world. A number of recent large events have increased community and
66scientific awareness of the need for effective detection, forecasting,
67and emergency preparedness. Probabilistic, geophysical and hydrodynamic
68models are required to predict the location and
69likelihood of an event, the initial sea floor deformation and
70subsequent propagation and inundation of the tsunami. Engineering, economic and social vulnerability models can then be used to estimate the
71impact of the event as well as the effectiveness of hazard mitigation
72procedures. In this paper, we focus on modelling of
73the physical processes only.
74
75Various approaches are currently used to assess the potential tsunami
76inundation of coastal communities.
77These methods differ in both the formulation used to
78describe the evolution of the tsunami and the numerical methods used
79to solve the governing equations. However any legitimate model must
80address each of the three distinct stages of tsunami evolution---
81generation, propagation and inundation. Models combining observed seismic,
82geodetic and sometimes tsunami data must be used
83to provide estimates of initial sea floor and ocean surface
84deformation. The complexity of these models ranges from empirical to
85non-linear three-dimensional mechanical models. The shallow water wave
86equations, linearised shallow water wave equations, and
87Boussinesq-type equations are frequently used to simulate tsunami
88propagation. These models are typically used to predict quantities
89such as arrival times, wave speeds and heights, and inundation extents
90which are used to develop efficient hazard mitigation plans.
91
92Inaccuracies in model prediction can result in inappropriate
93evacuation plans and town zoning, which may result in loss of life and
94large financial losses. Consequently tsunami models must undergo
95sufficient end-to-end testing to increase scientific and community
96confidence in the model predictions.
97
98Complete confidence in a model of a physical system cannot be
99established.  One can only hope to state under what conditions and to what extent the
100model hypothesis holds true. Specifically the utility of a model can
101be assessed through a process of verification and
102validation. Verification assesses the accuracy of the numerical method
103used to solve the governing equations and validation is used to
104investigate whether the model adequately represents the physical
105system~\cite{bates01}. Together these processes can be used to
106establish the likelihood that a model represents a legitimate
107hypothesis.
108
109The sources of data used to validate and verify a model can be
110separated into three main categories: analytical solutions, scale
111experiments and field measurements. Analytical solutions of the
112governing equations of a model, if available, provide the best means
113of verifying any numerical model. However, analytical solutions are
114frequently limited to a small set of idealised examples that do not
115completely capture the more complex behaviour of `real' events. Scale
116experiments, typically in the form of wave-tank experiments, provide a
117much more realistic source of data that better captures the complex
118dynamics of flows such as those generated by a tsunami, whilst allowing
119control of the event and much easier and accurate measurement of the
120tsunami properties. Comparison of numerical predictions with field
121data provides the most stringent test. The use of field data increases
122the generality and significance of conclusions made regarding model
123utility. On the other hand, it must be noted that the use of field
124data also significantly increases the uncertainty of the validation
125experiment that may constrain the ability to make unequivocal
126statements~\cite{bates01}.
127FIXME (Jane): Why would that increase the uncertainty?
128FIXME (Phil): references to all of the paragraph above, please
129
130Currently, the extent of tsunami-related field data is limited. The
131cost of tsunami monitoring programs, bathymetry and topography surveys
132prohibits the collection of data in many of the regions in which
133tsunamis pose greatest threat. The resulting lack of data has limited
134the number of field data sets available to validate tsunami
135models.
136
137Synolakis et al~\cite{synolakis07} have developed a set of
138standards, criteria and procedures for evaluating numerical models of
139tsunami. They propose three analytical solutions to help identify the
140validity of a model, and five scale comparisons (wave-tank benchmarks)
141and two field events to assess model veracity.
142
143The first field data benchmark introduced in \cite{synolakis07} compares model
144results against observed data from the Hokkaido-Nansei-Oki tsunami
145that occurred around Okushiri Island, Japan on the 12 July
1461993. This tsunami provides an example of extreme run-up generated from
147reflections and constructive interference resulting from local
148topography and bathymetry. The benchmark consists of two tide gauge
149records and numerous spatially-distributed point sites at which
150modelled maximum run-up elevations can be compared. The second
151benchmark is based upon the Rat Islands tsunami that occurred off the
152coast of Alaska on the 17 November 2003. The Rat Island tsunami
153provides a good test for real-time forecasting models since the tsunami
154was recorded at three tsunameters. The test requires matching the
155tsunami propagation model output with the DART recording to constrain the
156tsunami source model, and then using it to reproduce the tide gauge
157record at Hilo, Hawaii.
158FIXME (Jane): Are the tsunameters and the DART recordings the same thing?
159
160In this paper we develop a field data benchmark to be used in
161conjunction with the other tests proposed by Synolakis et
162al~\cite{synolakis07} to validate and verify tsunami models.
163The benchmark proposed here allows evaluation of
164model structure during all three distinct stages tsunami evolution.
165It consists of geodetic measurements of the
166Sumatra--Andaman earthquake that are used to validate the description
167of the tsunami source, altimetry data from the \textsc{jason} satellite to test
168open ocean propagation, eye-witness accounts to assess near shore
169propagation, and a detailed inundation survey of Patong city, Thailand
170to compare model and observed inundation. A description of the data
171required to construct the benchmark is given in
172Section~\ref{sec:data}.
173
174An associated aim of this paper is to illustrate the use of this new
175benchmark to validate a dedicated inundation model called
176\textsc{anuga} used by Geoscience Australia. A description of
177\textsc{anuga} is given in Section~\ref{sec:models} and the validation
178results are given in Section~\ref{sec:results}.
179
180The numerical models used to simulate tsunami impact
181are computationally intensive and high resolution models of the entire
182evolution process will often take a number of days to
183run. Consequently, the uncertainty in model predictions is difficult to
184quantify as it would require a very large number of runs.
185However, model uncertainty should not be ignored. Section
186~\ref{sec:sensitivity} provides a simple analysis that can
187be used to investigate the sensitivity of model predictions to model
188parameters.
189
190%================Section===========================
191\section{Data}\label{sec:data}
192The sheer magnitude of the 2004 Sumatra-Andaman earthquake and the
193devastation caused by the subsequent tsunami have generated much
194scientific interest. As a result an unusually large amount of post
195seismic data has been collected and documented. Data sets from
196seismometers, tide gauges, \textsc{gps} surveys, satellite overpasses,
197subsequent coastal field surveys of run-up and flooding, and
198measurements of coseismic displacements as well as bathymetry from ship-based
199expeditions, have now been made
200available. %~\cite{vigny05,amnon05,kawata05,liu05}. FIXME (Ole): Refs? 
201
202In this section we present the corresponding data necessary to implement
203the proposed benchmark for each of the three stages of the tsunami's evolution.
204
205\subsection{Generation}\label{sec:gen_data}
206All tsunami are generated from an initial disturbance of the ocean
207which develops into a low frequency wave that propagates outwards from
208the source. The initial deformation of the water surface is most
209commonly caused by coseismic displacement of the sea floor, but
210submarine mass failures, landslides, volcanoes or asteroids can also
211cause tsunami. In this section we detail the information used in
212this study to validate models of the sea floor deformation generated
213by the 2004 Sumatra--Andaman earthquake.
214
215The 2004 Sumatra--Andaman tsunami was generated by a coseismic
216displacement of the sea floor resulting from one of the largest
217earthquakes on record. The mega-thrust earthquake started on the 26
218December 2004 at 0h58'53'' UTC (or just before 8 am local time)
219approximately 70 km offshore of North Sumatra
220(\url{http://earthquake.usgs.gov/eqcenter/eqinthenews/2004/usslav}). The
221rupture propagated 1000-1300 km along the Sumatra-Andaman trench to
222the north at a rate of 2.5-3 km.s$^{-1}$ and lasted approximately 8-10
223minutes~\cite{ammon05}. Estimates of the moment magnitude of this
224event range from about 9.1 to 9.3 $M_w$~\cite{chlieh07,stein07}.
225
226The unusually large surface deformation caused by this earthquake
227means that there were a range of different geodetic measurements of
228the surface deformation available. These include field measurements of
229uplifted or subsided coral heads, continuous or campaign \textsc{GPS}
230measurements and remote sensing measurements of uplift or subsidence
231(see~\cite{chlieh07} and references therein). Here we use the the near-field
232estimates of vertical deformation in northwestern Sumatra and
233the Nicobar-Andaman islands collated by~\cite{chlieh07} to validate
234that our crustal deformation model of the 2004 Sumatra--Andaman
235earthquake is producing reasonable results. Note that the geodetic
236data used here is a combination of the vertical deformation that
237happened in the $\sim$10 minutes of the earthquake plus the
238deformation that followed in the days following the earthquake before
239each particular measurement was actually made (typically of order
240days). Therefore some of the observations may not contain the purely
241co-seismic deformation but could include some post-seismic deformation
242as well~\cite{chlieh07}.
243
244%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.
245
246%\begin{figure}[ht]
247%\begin{center}
248%\includegraphics[width=8.0cm,keepaspectratio=true]{geodeticMeasurements.jpg}
249%\caption{Near field geodetic measurements used to validate tsunami generation. FIXME: Insert appropriate figure here}
250%\label{fig:geodeticMeasurements}
251%\end{center}
252%\end{figure}
253
254\subsection{Propagation}
255\label{sec:propagation data}
256Once generated, a tsunami will propagate outwards from the source until
257it encounters the shallow water bordering coastal regions.
258FIXME (Ole): Need to change this definition. I believe propagation takes place all the way to the shore line and not just up to shallow waters.
259
260This period
261of the tsunami evolution is referred to as the propagation stage. The
262height and velocity of the tsunami is dependent on the local
263bathymetry in the regions through which the wave travels and the size
264of the initial wave. This section details the bathymetry data needed
265to model the tsunami propagation and the satellite altimetry transects
266used here to validate open ocean tsunami models.
267
268\subsubsection{Bathymetry Data}
269The bathymetry data used in this study was derived from the following
270sources:
271\begin{itemize}
272\item a two arc minute grid data set covering the Bay of Bengal,
273  DBDB2, obtained from US Naval Research Labs;
274\item a 3 second arc grid covering the whole of the Andaman Sea based
275  on Thai Navy charts no. 45 and no. 362; and
276\item a one second grid created from the digitised Thai Navy
277  bathymetry chart, no. 358, which covers Patong Bay and the
278  immediately adjacent regions.
279  (FIXME (Ole): How was the grid created from these digitised points?)
280\end{itemize}
281FIXME (Jane): Refs for all these.
282
283%A number of raw data sets were obtained, analysed and checked for
284%quality and subsequently gridded for easier visualisation and input
285%into the tsunami models.
286
287These sets were combined via
288interpolation and resampling to produce four nested grids
289which are relatively coarse in the deeper water and
290progressively finer as the distance to
291Patong Beach decreases as shown in Figure~\ref{fig:nested_grids}
292
293The coarsest
294bathymetry was obtained by interpolating the DBDB2 grid to a 27 second
295arc grid. A subsection of this region was then replaced by nine second
296data which was generated by sub-sampling the three second of arc grid from
297NOAA (FIXME (Jane): This was not mentioned in the dots above).
298
299A subset of the nine second grid was replaced by the three second
300data. Finally, the one second grid was used to approximate the
301bathymetry in Patong Bay and the immediately adjacent regions. Any
302points that deviated from the general trend near the boundary were
303deleted as a quality check.
304
305The sub-sampling of larger grids was performed by using {\bf resample},
306a Generic Mapping Tools (\textsc{GMT}) program (\cite{wessel98}). The
307gridding of data was performed using {\bf Intrepid}, a commercial
308geophysical processing package developed by Intrepid Geophysics. The
309gridding scheme employed the nearest neighbour algorithm followed by
310an application of minimum curvature akima spline smoothing.
311See \url{http://www.intrepid-geophysics.com/ig/manuals/english/gridding.pdf} 
312for details on the Intrepid model.
313
314
315\begin{figure}[ht]
316\begin{center}
317\includegraphics[width=0.75\textwidth,keepaspectratio=true]{nested_grids}
318\caption{Nested bathymetry grids.}
319\label{fig:nested_grids}
320\end{center}
321\end{figure}
322
323\subsubsection{JASON Satellite Altimetry}\label{sec:data_jason}
324During the 26 December 2004 event, the \textsc{jason} satellite tracked from
325north to south and over the equator at 02:55 UTC nearly two hours
326after the earthquake \cite{gower05}. The satellite recorded the sea
327level anomaly compared to the average sea level from its previous five
328passes over the same region in the 20-30 days prior. This data was
329used to validate the propagation stage in Section
330\ref{sec:resultsPropagation}.
331FIXME (Ole): See Phil's second point and email with help from David
332
333
334%DB I suggest we combine with model data to reduce the number of figures. The satellite track is shown in Figure~\ref{fig:satelliteTrack}.
335
336%\begin{figure}[ht]
337%\begin{center}
338%\includegraphics[width=8.0cm,keepaspectratio=true]{sateliteTrack.jpg}
339%\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.}
340%\label{fig:satelliteTrack}
341%\end{center}
342%\end{figure}
343
344%\begin{figure}[ht]
345%\begin{center}
346%\includegraphics[width=8.0cm,keepaspectratio=true]{jasonAltimetry.jpg}
347%\caption{JASON satellite altimetry seal level anomaly. FIXME: should we include figure here with just JASON altimetry.}
348%\label{fig:jasonAltimetry}
349%\end{center}
350%\end{figure}
351
352%FIXME: Can we compare the urs model against the TOPEX-poseidon satellite as well? DB No (we don't have the data currently).
353
354\subsection{Inundation}
355\label{sec:inundation data}
356FIXME (Ole): Technically propagation covers everything up to
357the coastline and inundation everything on-shore.
358This means that ANUGA covers the final part of the propagation and the inundation part. Should we adopt this distiction throughout the paper?
359
360Inundation refers to the final stages of the evolution of a tsunami and
361covers the propagation of the tsunami in coastal waters and the
362subsequent run-up onto land. This process is typically the most
363difficult of the three stages to model due to thin layers of water
364flowing rapidly over dry land.  Aside from requiring robust solvers
365which can simulate such complex flow patterns, this part of the
366modelling process also requires high resolution and quality elevation
367data which is often not available. In the case of model validation
368high quality field measurements are also required. For the proposed
369benchmark a high resolution bathymetry (FIXME (Ole): Bathymetry ?) and
370topography data set and a high quality inundation survey map from the
371Coordinating Committee Co-ordinating Committee for Geoscience Programmes
372in East and Southeast Asia (CCOP) (\cite{szczucinski06}) was obtained
373to validate model inundation. See also acknowledgements at the end of this paper.
374
375In this section we also present eye-witness accounts which can be used
376to qualitatively validate tsunami inundation.
377
378\subsubsection{Topography Data}
379A one second grid was used to approximate the topography in Patong
380Bay. This elevation data was again created from the digitised Thai
381Navy bathymetry chart, no 358.
382FIXME (Ole): I don't think so. The Navy chart is only offshore.
383
384 A visualisation of the elevation data
385set used in Patong Bay is shown in
386Figure~\ref{fig:patong_bathymetry}. The continuous topography
387(FIXME(Jane): What is meant by this?) is an
388interpolation of known elevation measured at the coloured dots. FIXME ??
389
390\begin{figure}[ht]
391\begin{center}
392\includegraphics[width=8.0cm,keepaspectratio=true]{patong_bay_data.jpg}
393\caption{3D visualisation of the elevation data set used in Patong Bay showing data points, contours, rivers and roads draped over the final model.}
394\label{fig:patong_bathymetry}
395\end{center}
396\end{figure}
397FIXME (Jane): legend? Were the contours derived from the final dataset?
398This is not the entire model, only the bay and the beach.
399
400\subsubsection{Buildings and Other Structures}
401Human-made buildings and structures can significantly affect tsunami
402inundation. The footprint and number of floors of the
403buildings in Patong Bay were extracted from a GIS data set which was also provided by the CCOP (see Section \ref{sec:inundation data} for details).
404The heights of these
405buildings were estimated assuming that each floor has a height of 3 m and they
406were added to the topographic dataset.
407
408\subsubsection{Inundation Survey}
409Tsunami run-up in built-up areas can be the cause of large financial and human
410losses, yet run-up data that can be used to validate model run-up
411predictions is scarce because such events are relatively infrequent.
412Of the two field benchmarks proposed
413in~\cite{synolakis07},
414only the Okushiri benchmark facilitates comparison between
415modelled and observed run-up. One of the major strengths of the
416benchmark proposed here is that modelled run-up can be compared to an
417inundation survey which maps the maximum run-up along an entire coastline
418rather than at a series of discrete sites. The survey map is
419shown in Figure~\ref{fig:patongescapemap} and plots the maximum run-up
420of the 2004 Indian Ocean tsunami in Patong city. Refer to Szczucinski et
421al~\cite{szczucinski06} for further details.
422
423\begin{figure}[ht]
424\begin{center}
425\includegraphics[width=8.0cm,keepaspectratio=true]{patongescapemap.jpg}
426\caption{Tsunami survey mapping the maximum observed inundation at
427  Patong beach courtesy of the CCOP \protect \cite{szczucinski06}.}
428\label{fig:patongescapemap}
429\end{center}
430\end{figure}
431
432
433\subsubsection{Eyewitness Accounts}\label{sec:eyewitness data}
434Eyewitness accounts detailed in~\cite{papadopoulos06}
435report that many people at Patong Beach observed an initial
436retreat (trough or draw down) of
437the shoreline of more than 100 m followed a few minutes later by a
438strong wave (crest). Another less powerful wave arrived another five
439or ten minutes later. Eyewitness statements \cite{papadopoulos06} 
440place the arrival time of
441the first wave between 9:55 am and 10:05 am local time or about 2 hours
442after the source rupture.
443FIXME (Ole): We should add observed arrival time and later relate that to
444the modelled dynamics.
445
446
447\begin{figure}[ht]
448\begin{center}
449\includegraphics[width=8.0cm,keepaspectratio=true]{gauge_locations.jpg}
450\caption{Location of timeseries extracted from the model output.}
451\label{fig:gauge_locations}
452\end{center}
453\end{figure}
454
455
456
457
458
459Two videos were sourced\footnote{The footage is
460widely available and can for example be obtained from
461\url{http://www.archive.org/download/patong_bavarian/patong_bavaria.wmv}
462(Comfort Hotel) and
463\url{http://www.archive.org/download/tsunami_patong_beach/tsunami_patong_beach.wmv}
464(Novotel)}
465%http://wizbangblog.com/content/2005/01/01/wizbang-tsunami.php
466which include footage of the tsunami in Patong Bay on the day
467of the 2004 Indian Ocean Tsunami. Both videos show an already inundated
468group of buildings. They also show what is to be assumed as the second
469and third waves approaching and further flooding of the buildings and
470street.  The first video is in the very north, filmed from what is
471believed to be the roof of the Novotel Hotel marked ``north'' in Figure
472\ref{fig:gauge_locations}. The second video is in the very south,
473filmed from the second story of a building next door to the Comfort
474Resort near the corner of Ruam Chai St and Thaweewong Road.  This
475location is marked ``south'' in Figure \ref{fig:gauge_locations}.
476Figure~\ref{fig:video_flow} shows stills from this video. Both videos
477were used to estimate flow speeds and inundation depths over time.
478
479\begin{figure}[ht]
480\begin{center}
481\includegraphics[width=5.0cm,keepaspectratio=true]{flow_rate_south_0_00sec.jpg}
482\includegraphics[width=5.0cm,keepaspectratio=true]{flow_rate_south_5_04sec.jpg}
483\includegraphics[width=5.0cm,keepaspectratio=true]{flow_rate_south_7_12sec.jpg}
484\includegraphics[width=5.0cm,keepaspectratio=true]{flow_rate_south_7_60sec.jpg}
485\caption{Four frames from a video where flow rate could be estimated,
486  circle indicates tracked debris, from top left: 0.0 sec, 5.0 s, 7.1
487  s, 7.6 s.}
488\label{fig:video_flow}
489\end{center}
490\end{figure}
491
492Flow rates were estimated using landmarks found in both videos and
493were found to be in the range of 5 to 7 metres per second (+/- 2 m/s)
494in the north and 0.5 to 2 metres per second (+/- 1 m/s) in the south.
495FIXME (Jane): How were these error bounds derived?
496Water depths could also
497be estimated from the videos by the level at which water rose up the
498sides of buildings such as shops. Our estimates are in the order of
4991.5 to 2.0 metres (+/- 0.5 m).
500Fritz ~\cite{fritz06} performed a detailed
501analysis of video frames taken around Banda Aceh and arrived at flow
502speeds in the range of 2 to 5 m/s.
503
504
505\subsection{Validation Check-List}
506\label{sec:checkList}
507The data described in this section can be used to construct a
508benchmark to validate all three stages of the evolution of a
509tsunami. In particular we propose that a legitimate tsunami model
510should reproduce the following behaviour:
511\begin{itemize}
512 \item reproduce the vertical deformation observed in north-western
513   Sumatra and along the Nicobar--Andaman islands (see
514   Section~\ref{sec:gen_data}),
515 \item reproduce the \textsc{jason} satellite altimetry sea surface
516   anomalies (see Section~\ref{sec:data_jason}),
517 \item reproduce the inundation survey map in Patong city
518   (Figure~\ref{fig:patongescapemap}),
519 \item simulate a leading depression followed by two distinct crests
520   of decreasing magnitude at the beach, and
521 \item predict the water depths and flow speeds, at the locations of
522   the eye-witness videos, that fall within the bounds obtained from
523   the videos.
524\end{itemize}
525
526Ideally, the model should also be compared to measured timeseries of
527waveheights and velocities but the authors are not aware of the
528availability of such data near Patong Bay.
529
530
531
532%================Section===========================
533\section{Modelling the Event}\label{sec:models}
534Numerous models are currently used to model and predict tsunami
535generation, propagation and run-up~\cite{titov97a,satake95}. Here we
536introduce the three part modelling methodology employed by Geoscience Australia
537to illustrate the utility of the proposed benchmark.
538
539\subsection{Generation}\label{sec:modelGeneration}
540FIXME (Ole and Jane): Does this need to be so long?
541
542There are various approaches to modelling the expected crustal
543deformation from an earthquake. Most approaches model the
544earthquake as a dislocation in a linear elastic medium. Here we use
545the method of Wang et al~\cite{wang03}. One of the main advantages
546of their method is that it allows the dislocation to be located in a
547stratified linear elastic half-space with an arbitrary number of
548layers. Other methods (such as those based on Okada's equations) can
549only model the dislocation in a homogeneous elastic half space, or can
550only include a limited number of layers, and thus cannot model the
551effect of the depth dependence of the elasticity of the
552Earth~\cite{wang03}. The original versions of the codes described here
553are available from \url{http://www.iamg.org/CGEditor/index.htm}. The
554first program, \textsc{edgrn}, calculates elastic Green's function for
555a set of point sources at a regular set of depths out to a specified
556distance. The equations controlling the deformation are solved by
557using a combination of Hankel's transform and Wang et al's
558implementation of the Thomson-Haskell propagator
559algorithm~\cite{wang03}. Once the Green's functions are calculated
560a slightly modified version of \textsc{edcmp}\footnote{For this study,
561we have made minor modifications
562to \textsc{edcmp} in order for it to provide output in a file format
563compatible with the propagation code in the following section. Otherwise it
564is similar to the original code.} is used to calculate the sea
565floor deformation for a specific subfault. This second code
566discretises the subfault into a set of unit sources and sums the
567elastic Green's functions calculated from \textsc{edgrn} for all the
568unit sources on the fault plane in order to calculate the final static
569deformation caused by a two dimensional dislocation along the
570subfault. This step is possible because of the linearity of the
571governing equations.
572
573In order to calculate the crustal deformation using these codes
574a model that describes the variation in elastic
575properties with depth and a slip model of the earthquake to describe
576the dislocation is required.
577The elastic parameters used for this study are the
578same as those in Table 2 of Burbidge et al~\cite{burbidge08}. For the slip
579model, there are many possible models for the 2004 Andaman--Sumatran
580earthquake to select from
581~\cite{chlieh07,asavanant08,arcas06,grilli07,ioualalen07}. Some are
582determined from various geological surveys of the site. Others solve
583an inverse problem which calibrates the source based upon the tsunami
584wave signal, the seismic signal and/or even the run-up.
585The source
586parameters used here to simulate the 2004 Indian Ocean tsunami were
587taken from the slip model G-M9.15 of Chlieh
588et al~\cite{chlieh07}. This model was created by inversion of wide
589range of geodetic and seismic data. The slip model consists of 686
59020km x 20km subsegments each with a different slip, strike and dip
591angle. The dip subfaults go from $17.5^0$ in the north and $12^0$ in
592the south. Refer to Chlieh et al~\cite{chlieh07} for a detailed
593discussion of this model and its derivation. Note that the geodetic
594data used in the validation was also included by~\cite{chlieh07} in
595the inversion used to find G-M9.15. Thus the validation is not
596completely independent. However, a reasonable validation would still
597show that the crustal deformation and elastic properties model used
598here is at least as valid as the one used by Chlieh
599et al~\cite{chlieh07} and can reproduce the observations just as
600accurately.
601
602\subsection{Propagation}\label{sec:modelPropagation}
603The \textsc{ursga} model described below was used to simulate the
604propagation of the 2004 Indian Ocean tsunami across the open ocean, based on a
605discrete representation of the initial deformation of the sea floor, as
606described in Section~\ref{sec:modelGeneration}. For the models shown
607here, the uplift is assumed to be instantaneous and creates a wave of
608the same size and amplitude as the co-seismic sea floor deformation.
609
610\subsubsection{URSGA}
611\textsc{ursga} is a hydrodynamic code that models the propagation of
612the tsunami in deep water using a finite difference method on a staggered grid.
613It solves the depth integrated linear or nonlinear shallow water equations in
614spherical co-ordinates with friction and Coriolis terms. The code is
615based on Satake~\cite{satake95} with significant modifications made by
616the \textsc{urs} corporation, Thio et al~\cite{thio08} and Geoscience
617Australia, Burbidge et al~\cite{burbidge08}.
618The tsunami was propagated via the nested
619grid system described in Section \ref{sec:propagation data} where
620the coarse grids were used in the open ocean and the finest
621resolution grid was employed in the region closest to Patong bay.
622\textsc{Ursga} is not publicly available.
623
624\subsection{Inundation}\label{sec:modelInundation}
625The utility of the \textsc{ursga} model decreases with water depth
626unless an intricate sequence of nested grids is employed. In
627comparison \textsc{anuga}, described below, is designed to produce
628robust and accurate predictions of on-shore inundation, but is less
629suitable for earthquake source modelling and large study areas because
630it is based on projected spatial coordinates. Consequently, the
631Geoscience Australia tsunami modelling methodology is based on a
632hybrid approach using models like \textsc{ursga} for tsunami
633propagation up to an offshore depth contour, typically 100 m.
634%Specifically we use the \textsc{ursga} model to simulate the
635%propagation of the 2004 Indian Ocean tsunami in the deep ocean, based
636%on a discrete representation of the initial deformation of the sea
637%floor, described in Section~\ref{sec:modelGeneration}.
638The wave signal and the velocity field is then used as a
639time varying boundary condition for
640the \textsc{anuga} inundation simulation.
641% A description of \textsc{anuga} is the following section.
642
643\subsubsection{ANUGA}
644\textsc{Anuga} is a Free and Open Source hydrodynamic inundation tool that
645solves the conserved form of the depth-integrated nonlinear shallow
646water wave equations using a Finite-Volume scheme on an
647unstructured triangular mesh.
648The scheme, first
649presented by Zoppou and Roberts~\cite{zoppou99}, is a high-resolution
650Godunov-type method that uses the rotational invariance property of
651the shallow water equations to transform the two-dimensional problem
652into local one-dimensional problems. These local Riemann problems are
653then solved using the semi-discrete central-upwind scheme of Kurganov
654et al~\cite{kurganov01} for solving one-dimensional conservation
655equations. The numerical scheme is presented in detail in
656Roberts and Zoppou~\cite{zoppou00,roberts00} and
657Nielsen et al~\cite{nielsen05}. An important capability of the
658finite-volume scheme is that discontinuities in all conserved quantities
659are allowed at every edge in the mesh. This means that the tool is
660well suited to adequately resolving hydraulic jumps, transcritical flows and
661the process of wetting and drying. This means that \textsc{Anuga} 
662is suitable for
663simulating water flow onto a beach or dry land and around structures
664such as buildings. \textsc{Anuga} has been validated against
665%a number of analytical solutions and 
666%FIXME (Ole): Analytical solutions have not been published. Ask Steve.
667the wave tank simulation of the 1993 Okushiri
668Island tsunami~\cite{nielsen05,roberts06}.
669FIXME (Ole): Add reference to Tom Baldock's Dam Break valiadation of ANUGA.
670
671
672%================Section===========================
673\section{Results}\label{sec:results}
674This section presents a validation of the modelling practice of Geoscience
675Australia against the new proposed benchmarks. The criteria outlined
676in Section~\ref{sec:checkList} are addressed for each of the three stages
677of tsunami evolution.
678
679\subsection{Generation}\label{modelGeneration}
680The location and magnitude of the sea floor displacement associated
681with the 2004 Sumatra--Andaman tsunami calculated from the G-M9.15
682model of~\cite{chlieh07} is shown in
683Figure~\ref{fig:surface_deformation}. The magnitude of the sea floor
684displacement ranges from about $-3.0$ to $5.0$ metres. The region near
685the fault is predicted to uplift, while that further away from the
686fault subsides. Also shown in Figure~\ref{fig:surface_deformation} are
687the areas that were observed to uplift (arrows pointing up) or subside
688(arrows point down) during and immediately after the earthquake. Most
689of this data comes from uplifted or subsided coral heads. The length of
690the vector increases with the magnitude of the displacement; the length
691corresponding to 1 m of observed motion is shown in the top right
692corner of the figure. As can be seen, the source model detailed in
693Section~\ref{sec:modelGeneration} produces a crustal deformation that
694matches the vertical displacements in the Nicobar-Andaman islands and
695Sumatra very well. Uplifted regions are close to the fault and
696subsided regions are further away. The crosses on
697Figure~\ref{fig:surface_deformation} show estimates of the pivot line
698from the remote sensing data~\cite{chlieh07} and they follow the
699predicted pivot line quite accurately. The average difference between
700the observed motion and the predicted motion (including the pivot line
701points) is only 0.06 m, well below the typical error of the
702observations of between 0.25 and 1.0 m. However, the occasional point
703has quite a large error (over 1 m); for example a couple of
704uplifted/subsided points appear to be on a wrong
705(FIXME (Jane): This is incorrect) side of the predicted
706pivot line~\ref{fig:surface_deformation}. The excellence of the fit is
707not surprising, since the original slip model was chosen
708by~\cite{chlieh07} to fit this (and the seismic data) well.
709This does demonstrate, however, that \textsc{edgrn} and our modified version of
710\textsc{edstat} (FIXME(Jane): This has never been mentioned before)
711can reproduce the correct pattern of vertical
712deformation very well when the slip distribution is well constrained
713and when reasonable values for the elastic properties are used.
714
715\begin{figure}[ht]
716\begin{center}
717\includegraphics[width=5cm,keepaspectratio=true]{surface_deformation.jpg}
718%\includegraphics[totalheight=0.3\textheight,width=0.8\textwidth]{surface_deformation.jpg}
719\caption{Location and magnitude of the vertical component of the sea
720  floor displacement associated with the 2004 Indian Ocean tsunami
721  based on the slip model, G-M9.15. The black arrows which point up
722  show areas observed to uplift during and immediately after the
723  earthquake; those pointing down are locations which subsided. The
724  length of the arrow increases with the magnitude of the deformation. The arrow
725  length corresponding to 1 m of deformation is shown in the top right
726  hand corner of the figure. The cross marks show the location of
727  the pivot line (the region between the uplift and subsided region
728  where the uplift is zero) derived from remote sensing
729  (FIXME(Jane): How was that possible?). All the
730  observational data are from the dataset collated
731  by~\cite{chlieh07}.}
732\label{fig:surface_deformation}
733\end{center}
734\end{figure}
735
736\subsection{Propagation}\label{sec:resultsPropagation}
737The deformation results described in Section~\ref{sec:modelGeneration}
738were used to provide a profile of the initial ocean surface
739displacement. This wave was used as an initial condition for
740\textsc{ursga} and was propagated throughout the Bay of Bengal. The
741rectangular computational domain of the largest grid extended from
74290$^0$ to 100$^0$ East and 0 to 15$^0$ North and contained
7431335$\times$1996 finite difference points. Inside this grid, a nested
744sequence of grids was used. The grid resolution of the nested grids
745went from 27 arc seconds in the coarsest grid, down to nine arc seconds
746in the second grid, three arc seconds in the third grid and finally one arc
747second in the finest grid near Patong. The computational domain is
748shown in Figure~\ref{fig:computational_domain}.
749
750FIXME (Ole): I know that a nested ursga model was trialled for the
751end-to-end modelling. However, for the study done here, where models
752were coupled, I didn't think nested grids were used with URSGA -
753and certainly not down to 1 arc second. Can someone shed some light
754on this please?
755
756\begin{figure}[ht]
757\begin{center}
758%\includegraphics[width=5.0cm,keepaspectratio=true]{extent_of_ursga_model.jpg}
759%\includegraphics[width=5.0cm,keepaspectratio=true]{ursgaDomain.jpg}
760\includegraphics[width=5.0cm,keepaspectratio=true]{extent_of_ANUGA_model.jpg}
761\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).}
762\label{fig:computational_domain}
763\end{center}
764\end{figure}
765
766
767Figure \ref{fig:jasonComparison} provides a comparison of the
768\textsc{ursga}-predicted sea surface elevation with the \textsc{jason}
769satellite altimetry data. The \textsc{ursga} model replicates the
770amplitude and timing of the the wave observed at $2.5^0$ South,
771but underestimates the amplitude of the wave further to the south at
772$4^0$ South. In the model, the southern most of these two waves
773appears only as a small bump in the cross section of the model (shown
774in Figure~\ref{fig:jasonComparison}) instead of being a distinct peak
775as can be seen in the satellite data. Also note
776that the \textsc{ursga} model prediction of the ocean surface
777elevation becomes out of phase with the \textsc{jason} 
778data at $3^0$ to $7^0$ North
779latitude. Chlieh et al~\cite{chlieh07} also observed these misfits and
780suggest it is caused by a reflected wave from the Aceh Peninsula that
781is not resolved in the model due to insufficient resolution of the
782computational mesh and bathymetry data. This is also a limitation of
783the model presented here which could be improved by nesting
784grids near Aceh.
785
786\begin{figure}[ht]
787\begin{center}
788\includegraphics[width=12.0cm,keepaspectratio=true]{jasonComparison.jpg}
789\caption{Comparison of the \textsc{ursga}-predicted surface elevation
790  with the \textsc{jason} satellite altimetry data. The \textsc{ursga} wave
791  heights have been corrected for the time the satellite passed
792  overhead compared to \textsc{jason} sea level anomaly.}
793\label{fig:jasonComparison}
794\end{center}
795\end{figure}
796FIXME (Jane): This graph does not look nice. The legend URS Model should
797be URSGA model.
798
799\subsection{Inundation}
800After propagating the tsunami in the open ocean using \textsc{ursga},
801the approximated ocean and surface elevation and horisontal flow
802velocities were extracted and used to construct a boundary condition
803for the \textsc{anuga} model. The interface between the \textsc{ursga}
804and \textsc{anuga} models was chosen to roughly follow the 100~m depth
805contour along the west coast of Phuket Island. The computational
806domain is shown in Figure~\ref{fig:computational_domain}.
807
808The domain was discretised into 386,338 triangles. The resolution of
809the grid was increased in regions inside the bay and on-shore to
810efficiently increase the simulation accuracy for the impact area.
811The grid resolution ranged between a
812maximum triangle area of $1\times 10^5$ m$^2$ 
813(corresponding to approximately 440 m between mesh points)
814near the western ocean
815boundary to $20$ m$^2$ (corresponding to
816approximately 6 m between mesh points)
817in the small regions surrounding the inundation
818region in Patong Bay. The coarse resolution was chosen to be
819commensurate with the model output from the \textsc{ursga} model
820(FIXME - this has to be clearly stated in ursga section)
821while the latter was chosen to match the available resolution of topographic
822data and building data in Patong city.
823Due to a lack of available roughness data, friction was
824set to a constant throughout the computational domain. For the
825reference simulation, a Manning's coefficient of 0.01 was chosen to
826represent a small resistance to the water flow. See Section
827\ref{sec:friction sensitivity} for details on model sensitivity to
828this parameter.
829
830
831The boundary condition at each side of the domain towards the south
832and the north where no information about the incident wave or
833its velocity field is available
834was chosen as a transmissive
835boundary condition, effectively replicating the time dependent wave
836height present just inside the computational domain.
837The velocity field on these boundaries was set
838to zero. Other choices include applying the mean tide value as a
839Dirichlet boundary condition. But experiments as well as the
840result of the verification reported here showed that this approach
841tends to underestimate the tsunami impact due to the tempering of the
842wave near the side boundaries, whereas the transmissive boundary
843condition robustly preserves the wave.
844
845During the \textsc{anuga} simulation the tide was kept constant at
846$0.80$ m. This value was chosen to correspond to the tidal height
847specified by the Thai Navy tide charts
848(\url{http://www.navy.mi.th/hydro/}) at the time the tsunami arrived
849at Patong Bay. Although the tsunami propagated for approximately three
850hours before it reach Patong Bay, the period of time during which the
851wave propagated through the \textsc{anuga} domain is much
852smaller. Consequently the assumption of constant tide height is
853reasonable.
854
855Maximum onshore inundation depth was computed from the model
856throughout the entire Patong Bay region and used to generate
857a measure of the inundated area.
858Figure~\ref{fig:inundationcomparison1cm} (left) shows very good
859agreement between the measured and simulated inundation. However
860these results are dependent on the classification used to determine
861whether a region in the numerical simulation was inundated. In
862Figure~\ref{fig:inundationcomparison1cm} (left) a point in the computational
863domain was deemed inundated if at some point in time it was covered by
864at least 1 cm of water. However, the precision of the inundation boundary
865generated by the on-site survey is most likely less than that as it
866was determined by observing water marks and other signs
867left by the receding waters. Consequently the measurement error along
868the inundation boundary of the survey is likely to vary significantly
869and somewhat unpredictably.
870An inundation threshold of 10 cm therefore was selected for inundation
871extents reported in this paper to reflect
872the more likely accuracy of the survey, and subsequently facilitate a more
873appropriate comparison between the modelled and observed inundation
874area.
875Figure~\ref{fig:inundationcomparison1cm} (right) shows the simulated
876inundation using a larger threshold of 10 cm.
877
878
879The datasets necessary for reproducing the results
880of the inundation stage are available on Sourceforge under the \textsc{anuga}
881project (\url{http://sourceforge.net/projects/anuga}).
882At the time of
883writing the direct link is \url{http://tinyurl.com/patong2004-data}.
884%%\url{http://sourceforge.net/project/showfiles.php?group_id=172848&package_id=319323&release_id=677531}.
885The scripts required are part of the \textsc{anuga} distribution also
886available from Sourceforge \url{http://sourceforge.net/projects/anuga} under
887the validation section.
888
889An 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}.
890%\url{https://datamining.anu.edu.au/anuga/attachment/wiki/AnugaPublications/patong_2004_indian_ocean_tsunami_ANUGA_animation.mov}.
891
892\begin{figure}[ht]
893\begin{center}
894\includegraphics[width=6.0cm,keepaspectratio=true]{final_1cm.jpg}
895\includegraphics[width=6.0cm,keepaspectratio=true]{final_10cm.jpg}
896\caption{Simulated inundation versus observed inundation using an
897inundation threshold of 1cm (left) and 10cm (right).}
898\label{fig:inundationcomparison1cm}
899\end{center}
900\end{figure}
901
902To quantify the agreement between the observed and simulated inundation we
903introduce the measure
904\begin{equation}
905\rho_{in}=\frac{A(I_m\cap I_o)}{A(I_o)}
906\end{equation}
907representing the ratio $\rho_{in}$ of the observed
908inundation region $I_o$ captured by the model $I_m$. Another useful
909measure is the fraction of the modelled inundation area that falls
910outside the observed inundation area given by the formula
911\begin{equation}
912\rho_{out}=\frac{A(I_m\setminus (I_m\cap I_o))}{A(I_o)}
913\end{equation}
914These values for the two aforementioned simulations are given in
915Table~\ref{table:inundationAreas}. High value of both $\rho_{in}$ and $\rho_{out}$ indicate
916that the model overestimates the impact whereas low values of both quantities would indicate
917an underestimation. A high value of $\rho_{in}$ combined with a low value of $\rho_{out}$ 
918indicates a good model prediction of the survey.
919
920Discrepancies between the survey data and the modelled inundation
921include: unknown distribution of surface roughness, inappropriate
922parameterisation of the source model, effect of humans structures on
923flow, as well as uncertainties in the elevation data, effects of
924erosion and deposition by the tsunami event,
925measurement errors in the GPS survey recordings, and
926missing data in the field survey data itself. The impact of some of
927these sources of uncertainties are is investigated in
928Section~\ref{sec:sensitivity}
929
930\subsection{Eye-witness accounts}
931Figure \ref{fig:gauge_locations} shows four locations where time
932series have been extracted from the model. The two offshore time series
933are shown in Figure \ref{fig:offshore_timeseries} and the two onshore
934timeseries are shown in Figure \ref{fig:onshore_timeseries}. The
935latter coincide with locations where video footage from the event is
936available as described in Section \ref{sec:eyewitness data}.
937
938\begin{figure}[ht]
939\begin{center}
940\includegraphics[width=10.0cm,keepaspectratio=true]{gauge_bay_depth.jpg}
941\includegraphics[width=10.0cm,keepaspectratio=true]{gauge_bay_speed.jpg}
942\caption{Time series obtained from the two offshore gauge locations,
9437C and 10C, shown in Figure \protect \ref{fig:gauge_locations}.}
944\end{center}
945\label{fig:offshore_timeseries}
946\end{figure}
947
948\begin{figure}[ht]
949\begin{center}
950\includegraphics[width=10.0cm,keepaspectratio=true]{gauges_hotels_depths.jpg}
951\includegraphics[width=10.0cm,keepaspectratio=true]{gauges_hotels_speed.jpg}
952\caption{Time series obtained from the two onshore locations, North and South,
953shown in Figure \protect \ref{fig:gauge_locations}.}
954\end{center}
955\label{fig:onshore_timeseries}
956\end{figure}
957
958
959The estimated depths and flow rates given in Section
960\ref{sec:eyewitness data} are shown together with the modelled depths
961and flow rates obtained from the model in
962Table \ref{tab:depth and flow comparisons}.
963The predicted maximum depths and speeds are all of the same order
964of what was observed. However, unlike the real event,
965the model estimates complete withdrawal of the water between waves at the
966chosen locations and shows that the model must be used with caution at this
967level of detail.
968Nonetheless, this comparison serves to check that depths and speeds
969predicted are within the range of what is expected.
970
971
972\begin{table}
973\[
974  \begin{array}{|l|cc|cc|}
975  \hline
976                 & \multicolumn{2}{|c|}{\mbox{Depth [m]}}
977                 & \multicolumn{2}{c|}{\mbox{Flow [m/s]}} \\ 
978                 & \mbox{Observed} & \mbox{Modelled}
979                 & \mbox{Observed} & \mbox{Modelled} \\ \cline{2-5}                 
980    \mbox{North} & 1.5-2 & 1.4 & 5-7 & 0.1 - 3.3 \\
981    \mbox{South} & 1.5-2 & 1.5 & 0.5-2 & 0.2 - 2.6 \\ \hline
982  \end{array}
983\]
984\label{tab:depth and flow comparisons}
985\end{table} 
986FIXME (Jane): We should perhaps look at average data in area surrounding these points
987
988%can be estimated with landmarks found in
989%satellite imagery and the use of a GIS and were found to be in the
990%range of 5 to 7 metres per second (+/- 2 m/s) in the north and 0.5 to
991%2 metres per second (+/- 1 m/s) in the south.
992
993Given the uncertainties in both model and observations, there is agreement
994between the values obtained from the videos and the simulations.
995
996% Our modelled flow rates show
997%maximum values in the order of 0.2 to 2.6 m/s in the south and 0.1 to
998%3.3 m/s for the north as shown in the figures. Water depths could also
999%be estimated from the videos by the level at which water rose up the
1000%sides of buildings such as shops. Our estimates are in the order of
1001%1.5 to 2.0 metres (+/- 0.5 m). This is in the same range as our
1002%modelled maximum depths of 1.4 m in the north and 1.5 m in the south
1003%as seen in the figure.
1004
1005
1006
1007
1008
1009%================Section===========================
1010\section{Sensitivity Analysis}
1011\label{sec:sensitivity}
1012This section investigates the effect of different values of Manning's
1013friction coefficient, changing waveheight at the 100 m depth contour,
1014and the presence and absence of buildings in the elevation dataset on
1015model maximum inundation. The reference model is the one reported in
1016Figure~\ref{fig:inundationcomparison1cm} (right) with a friction coefficient of 0.01,
1017buildings included and the boundary condition produced by the
1018\textsc{ursga} model.
1019
1020%========================Friction==========================%
1021\subsection{Friction}
1022\label{sec:friction sensitivity}
1023The first sensitivity study investigated the impact of surface roughness on the
1024predicted run-up. According to Schoettle~\cite{schoettle2007}
1025appropriate values of Manning's coefficient range from 0.007 to 0.03
1026for tsunami propagation over a sandy sea floor and the reference model
1027uses a value of 0.01.  To investigate sensitivity to this parameter,
1028we simulated the maximum onshore inundation using a Manning's
1029coefficient of 0.0003 and 0.03. The resulting inundation maps are
1030shown in Figure~\ref{fig:sensitivity_friction} and the maximum flow
1031speeds in Figure~\ref{fig:sensitivity_friction_speed}. These figures
1032show that the on-shore inundation extent decreases with increasing
1033friction and that small perturbations in the friction cause bounded
1034changes in the output. This is consistent with the conclusions of
1035Synolakis~\cite{synolakis05} et al, who state that the long wavelength of
1036tsunami tends to mean that friction is less important in
1037comparison to the motion of the wave.
1038
1039%========================Wave-Height==========================%
1040\subsection{Input Wave Height}\label{sec:waveheightSA}
1041The effect of the wave height used as input to the inundation model
1042\textsc{anuga} was also investigated.
1043Figure~\ref{fig:sensitivity_boundary} indicates that the inundation
1044severity is directly proportional to the boundary waveheight but small
1045perturbations in the input wave height of 10 cm appear to have little
1046effect on the final inundated area. Obviously larger perturbations
1047will have greater impact. However, wave heights in the open ocean are
1048generally well
1049predicted by the generation and propagation models such as
1050\textsc{ursga} as demonstrated in Section \ref{sec:resultsPropagation} 
1051and also in \cite{thomas2009}.
1052
1053
1054
1055%========================Buildings==========================%
1056\subsection{Buildings and Other Structures}
1057The presence or absence of physical buildings in the elevation model was also
1058investigated.
1059Figure~\ref{fig:sensitivity_nobuildings}
1060shows the inundated area and the associated maximum flow speeds
1061in the presence and absence of buildings. It
1062is apparent that densely built-up areas act as
1063dissipators greatly reducing the inundated area. However, flow speeds
1064tend to increase in passages between buildings.
1065 
1066
1067\begin{table}
1068\begin{center}
1069\label{table:inundationAreas}
1070\caption{$\rho_{in}$ and $\rho_{out}$ of the reference simulation and all sensitivity studies.}
1071\begin{tabular}{|l|c|c|}
1072\hline
1073 & $\rho_{in}$ & $\rho_{out}$ \\ 
1074\hline\hline
1075Reference model & 0.79 & 0.20\\ 
1076Friction = 0.0003 & 0.83 & 0.26 \\ 
1077Friction = 0.03 & 0.67 & 0.09\\ 
1078Boundary wave hight minus 10 cm & 0.77 & 0.17 \\
1079Boundary wave hight plus 10 cm & 0.82 & 0.22 \\
1080No Buildings & 0.94 & 0.44 \\
1081\hline 
1082\end{tabular}
1083\end{center}
1084\end{table}
1085
1086%================Section===========================
1087
1088\section{Conclusion}
1089This paper proposes an additional field data benchmark for the
1090verification of tsunami inundation models. Currently, there is a
1091scarcity of appropriate validation datasets due to a lack of well-documented
1092historical tsunami impacts. The benchmark proposed here
1093utilises the uniquely large amount of observational data for model
1094comparison obtained during, and immediately following, the
1095Sumatra--Andaman tsunami of 26 December 2004. Unlike the small
1096number of existing benchmarks, the proposed test validates all three
1097stages of tsunami evolution - generation, propagation and
1098inundation. In an attempt to provide higher visibility and easier
1099accessibility for tsunami benchmark problems, the data used to
1100construct the proposed benchmark is documented and freely available at
1101\url{http://tinyurl.com/patong2004-data}.
1102
1103This study also shows that the tsunami impact modelling methodology
1104adopted is credible and able to predict inundation extents with reasonable
1105accuracy.  An associated aim of this paper was to further validate the
1106hydrodynamic modelling tool \textsc{anuga} which is used to simulate
1107the tsunami inundation. Model predictions
1108matched well the geodetic measurements of the Sumatra--Andaman earthquake,
1109altimetry data from the \textsc{jason}, eye-witness accounts of wave
1110front arrival times and flow speeds and a detailed inundation survey
1111of Patong Bay, Thailand.
1112
1113A simple sensitivity analysis was performed to assess the influence of
1114small changes in friction, wave height at the 100 m depth contour and
1115the presence of buildings and other structures on the model
1116predictions. Of these three, the presence of buildings was shown to
1117have the greatest influence on
1118the simulated inundation extent. The value of friction and small
1119perturbations in the waveheight at the \textsc{anuga} boundary have
1120comparatively little effect on the model results.
1121
1122%================Acknowledgement===================
1123\section*{Acknowledgements}
1124This project was undertaken at Geoscience Australia and the Department
1125of Mathematics, The Australian National University. The authors would
1126like to thank Niran Chaimanee from the CCOP for providing
1127the post 2004 tsunami survey data, building footprints, aerial
1128photography and the elevation data for Patong city, Prapasri Asawakun
1129from the Suranaree University of Technology and Parida Kuneepong for
1130supporting this work; and Drew Whitehouse from the Australian National
1131University for preparing the animation of the simulated impact.
1132
1133\clearpage
1134\section{Appendix}
1135
1136This appendix present the images used to assess the model sensitivities described in
1137Section~\ref{sec:sensitivity}.
1138
1139\begin{figure}[ht]
1140\begin{center}
1141\includegraphics[width=6cm,keepaspectratio=true]{sensitivity_reference_depth}
1142\includegraphics[width=6cm,keepaspectratio=true]{sensitivity_reference_speed}
1143\caption{Results from reference model as reported in Section \protect \ref{sec:results},
1144  i.e.\ including buildings and a friction value of 0.01. The seaward boundary condition is as
1145  provided by the \textsc{ursga} model. The left image shows the maximum
1146  modelled depth while the right hand image shows the maximum modelled
1147  flow velocities.}
1148\label{fig:reference_model}
1149\end{center}
1150\end{figure}
1151
1152
1153
1154\begin{figure}[ht]
1155\begin{center}
1156\includegraphics[width=6cm,keepaspectratio=true]{sensitivity_minus10cm_depth}
1157\includegraphics[width=6cm,keepaspectratio=true]{sensitivity_plus10cm_depth}
1158\caption{Model results with wave height at \textsc{anuga} boundary artificially
1159  modified to assess sensitivities. The reference inundation extent is shown in Figure
1160  \protect \ref{fig:reference_model} (left).  The left and right images
1161  show the inundation results if the wave at the \textsc{anuga} boundary
1162  is reduced or increased by 10 cm respectively. The inundation
1163  severity varies in proportion to the boundary waveheight, but the
1164  model results are only slightly sensitive to this parameter for the
1165  range of values tested.}
1166\label{fig:sensitivity_boundary}
1167\end{center}
1168\end{figure}
1169FIXME (Jane): How and why was the +/- 10 cm chosen?
1170
1171
1172\begin{figure}[ht]
1173\begin{center}
1174\includegraphics[width=6cm,keepaspectratio=true]{sensitivity_minus10cm_speed}
1175\includegraphics[width=6cm,keepaspectratio=true]{sensitivity_plus10cm_speed}
1176\caption{The maximal flow speeds for the same model parameterisations
1177  found in Figure \protect \ref{fig:sensitivity_boundary}. The
1178  reference flow speeds are shown in Figure \protect
1179  \ref{fig:reference_model} (right).}
1180\label{fig:sensitivity_boundary_speed}
1181\end{center}
1182\end{figure}
1183
1184\begin{figure}[ht]
1185\begin{center}
1186\includegraphics[width=6cm,keepaspectratio=true]{sensitivity_nobuildings_depth}
1187\includegraphics[width=6cm,keepaspectratio=true]{sensitivity_nobuildings_speed}
1188\caption{Model results show the effect of buildings in
1189  the elevation data set.
1190  The left hand image shows the maximum inundation depth results for
1191  a model entirely without buildings.  As expected, the absence of
1192  buildings will increase the inundation extent beyond what was
1193  surveyed. The right hand image shows the corresponding flow speeds in the absence of buildings. 
1194  The reference results are as shown in Figure
1195  \protect \ref{fig:reference_model}.}
1196\label{fig:sensitivity_nobuildings}
1197\end{center}
1198\end{figure}
1199
1200
1201\begin{figure}[ht]
1202\begin{center}
1203%\includegraphics[width=6cm,keepaspectratio=true]{sensitivity_f0_0003_depth}
1204%\includegraphics[width=6cm,keepaspectratio=true]{sensitivity_f0_03_depth}
1205\includegraphics[width=12cm,keepaspectratio=true]{friction_comparison_depth}
1206\caption{Model results for different values of Manning's friction
1207  coefficient shown to assess sensitivities. The reference inundation extent for a
1208  friction value of 0.01 is shown in Figure
1209  \protect \ref{fig:reference_model} (left).  The left and right images
1210  show the inundation results for friction values of 0.0003 and
1211  0.03 respectively. The inundation extent increases for the lower
1212  friction value while the higher slows the flow and decreases the
1213  inundation extent. Ideally, friction should vary across the entire
1214  domain depending on terrain and vegetation, but this is beyond the
1215  scope of this study.}
1216\label{fig:sensitivity_friction}
1217\end{center}
1218\end{figure}
1219
1220\begin{figure}[ht]
1221\begin{center}
1222\includegraphics[width=6cm,keepaspectratio=true]{sensitivity_f0_0003_speed}
1223\includegraphics[width=6cm,keepaspectratio=true]{sensitivity_f0_03_speed}
1224\caption{The maximal flow speeds for the same model parameterisations
1225  found in Figure \protect \ref{fig:sensitivity_friction}. The
1226  reference flow speeds are shown in Figure \protect
1227  \ref{fig:reference_model} (right).}
1228\label{fig:sensitivity_friction_speed}
1229\end{center}
1230\end{figure}
1231
1232\clearpage
1233
1234%====================Bibliography==================
1235\bibliographystyle{spmpsci}
1236\bibliography{tsunami07}
1237\end{document}
Note: See TracBrowser for help on using the repository browser.