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

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

Jane's (and some of mine) comments and suggestions.
See log file and FIXME's

File size: 59.7 KB
2% to be added when submitted to ocean dynamics
5%\journalname{Ocean Dynamics}
10\usepackage{url}      % for URLs and DOIs
15\title{Benchmarking Tsunami Models using the December 2004 Indian
16  Ocean Tsunami and its Impact at Patong Bay}
19\author{J.~D. Jakeman \and O. Nielsen \and K. VanPutten \and
20  D. Burbidge \and R. Mleczko \and N. Horspool}
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{}
26%       \and
27%       O. Nielsen \and R. Mleczko \and D. Burbidge \and K. VanPutten \and N. Horspool \at
28%       Geoscience Australia, Canberra, \textsc{Australia}
32%================Start of Document================
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
56% to be added when submitted to ocean dynamics
57%\keywords{Tsunami \and modelling \and validation and verification \and benchmark}
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.
75Various approaches are currently used to assess the potential impact
76of tsunami. These methods differ in both the formulation used to
77describe the evolution of the tsunami and the numerical methods used
78to solve the governing equations. However any legitimate model must
79address each of the three distinct stages of tsunami evolution---
80generation, propagation and inundation. Models combining observed seismic,
81geodetic and sometimes tsunami data must be used
82to provide estimates of initial sea floor and ocean surface
83deformation. The complexity of these models ranges from empirical to
84non-linear three-dimensional mechanical models. The shallow water wave
85equations, linearised shallow water wave equations, and
86Boussinesq-type equations are frequently used to simulate tsunami
87propagation. These models are typically used to predict quantities
88such as arrival times, wave speeds and heights, and inundation extents
89which are used to develop efficient hazard mitigation plans.
91Inaccuracies in model prediction can result in inappropriate
92evacuation plans and town zoning, which may result in loss of life and
93large financial losses. Consequently tsunami models must undergo
94sufficient end-to-end testing to increase scientific and community
95confidence in the model predictions.
97Complete confidence in a model of a physical system cannot be
98established.  One can only hope to state under what conditions and to what extent the
99model hypothesis holds true. Specifically the utility of a model can
100be assessed through a process of verification and
101validation. Verification assesses the accuracy of the numerical method
102used to solve the governing equations and validation is used to
103investigate whether the model adequately represents the physical
104system~\cite{bates01}. Together these processes can be used to
105establish the likelihood that a model represents a legitimate
108The sources of data used to validate and verify a model can be
109separated into three main categories: analytical solutions, scale
110experiments and field measurements. Analytical solutions of the
111governing equations of a model, if available, provide the best means
112of verifying any numerical model. However, analytical solutions are
113frequently limited to a small set of idealised examples that do not
114completely capture the more complex behaviour of `real' events. Scale
115experiments, typically in the form of wave-tank experiments, provide a
116much more realistic source of data that better captures the complex
117dynamics of flows such as those generated by a tsunami, whilst allowing
118control of the event and much easier and accurate measurement of the
119tsunami properties. Comparison of numerical predictions with field
120data provides the most stringent test. The use of field data increases
121the generality and significance of conclusions made regarding model
122utility. On the other hand, it must be noted that the use of field
123data also significantly increases the uncertainty of the validation
124experiment that may constrain the ability to make unequivocal
125statements~\cite{bates01}. FIXME (Jane): Because?
126FIXME (Phil): references to all of the paragraph above, please
128Currently, the extent of tsunami-related field data is limited. The
129cost of tsunami monitoring programs, bathymetry and topography surveys
130prohibits the collection of data in many of the regions in which
131tsunamis pose greatest threat. The resulting lack of data has limited
132the number of field data sets available to validate tsunami
135Synolakis et al~\cite{synolakis07} have developed a set of
136standards, criteria and procedures for evaluating numerical models of
137tsunami. They propose three analytical solutions to help identify the
138validity of a model, and five scale comparisons (wave-tank benchmarks)
139and two field events to assess model veracity.
141The first field data benchmark introduced in \cite{synolakis07} compares model
142results against observed data from the Hokkaido-Nansei-Oki tsunami
143that occurred around Okushiri Island, Japan on the 12 July
1441993. This tsunami provides an example of extreme run-up generated from
145reflections and constructive interference resulting from local
146topography and bathymetry. The benchmark consists of two tide gauge
147records and numerous spatially-distributed point sites at which
148modelled maximum run-up elevations can be compared. The second
149benchmark is based upon the Rat Islands tsunami that occurred off the
150coast of Alaska on the 17 November 2003. The Rat Island tsunami
151provides a good test for real-time forecasting models since the tsunami
152was recorded at three tsunameters. The test requires matching the
153tsunami propagation model output with the DART recording to constrain the
154tsunami source model, and then using it to reproduce the tide gauge
155record at Hilo, Hawaii.
156FIXME (Jane): Are the tsunameters and the DART recordings the same thing?
158In this paper we develop a field data benchmark to be used in
159conjunction with the other tests proposed by Synolakis et
160al~\cite{synolakis07} to validate and verify tsunami models.
161The benchmark proposed here allows evaluation of
162model structure during all three distinct stages tsunami evolution.
163It consists of geodetic measurements of the
164Sumatra--Andaman earthquake that are used to validate the description
165of the tsunami source, altimetry data from the JASON satellite to test
166open ocean propagation, eye-witness accounts to assess near shore
167propagation, and a detailed inundation survey of Patong city, Thailand
168to compare model and observed inundation. A description of the data
169required to construct the benchmark is given in
172An associated aim of this paper is to illustrate the use of this new
173benchmark to validate a dedicated inundation model called
174\textsc{anuga} used by Geoscience Australia. A description of
175\textsc{anuga} is given in Section~\ref{sec:models} and the validation
176results are given in Section~\ref{sec:results}.
178The numerical models used to simulate tsunami impact
179are computationally intensive and high resolution models of the entire
180evolution process will often take a number of days to
181run. Consequently, the uncertainty in model predictions is difficult to
182quantify as it would require a very large number of runs.
183However, model uncertainty should not be ignored. Section
184~\ref{sec:sensitivity} provides a simple analysis that can
185be used to investigate the sensitivity of model predictions to model
190The sheer magnitude of the 2004 Sumatra-Andaman earthquake and the
191devastation caused by the subsequent tsunami have generated much
192scientific interest. As a result an unusually large amount of post
193seismic data has been collected and documented. Data sets from
194seismometers, tide gauges, \textsc{gps} surveys, satellite overpasses,
195subsequent coastal field surveys of run-up and flooding, and
196measurements of coseismic displacements as well as bathymetry from ship-based
197expeditions, have now been made
198available. %~\cite{vigny05,amnon05,kawata05,liu05}. FIXME (Ole): Refs? 
199In this section we present the corresponding data necessary to implement
200the proposed benchmark for each of the three stages of the tsunami's evolution.
203All tsunami are generated from an initial disturbance of the ocean
204which develops into a low frequency wave that propagates outwards from
205the source. The initial deformation of the water surface is most
206commonly caused by coseismic displacement of the sea floor, but
207submarine mass failures, landslides, volcanoes or asteroids can also
208cause tsunami. In this section we detail the information used in
209this study to validate models of the sea floor deformation generated
210by the 2004 Sumatra--Andaman earthquake.
212The 2004 Sumatra--Andaman tsunami was generated by a coseismic
213displacement of the sea floor resulting from one of the largest
214earthquakes on record. The mega-thrust earthquake started on the 26
215December 2004 at 0h58'53'' UTC (or just before 8 am local time)
216approximately 70 km offshore of North Sumatra
217(\url{}). The
218rupture propagated 1000-1300 km along the Sumatra-Andaman trench to
219the north at a rate of 2.5-3 km.s$^{-1}$ and lasted approximately 8-10
220minutes~\cite{ammon05}. Estimates of the moment magnitude of this
221event range from about 9.1 to 9.3 $M_w$~\cite{chlieh07,stein07}.
223The unusually large surface deformation caused by this earthquake
224means that there were a range of different geodetic measurements of
225the surface deformation available. These include field measurements of
226uplifted or subsided coral heads, continuous or campaign \textsc{GPS}
227measurements and remote sensing measurements of uplift or subsidence
228(see~\cite{chlieh07} and references therein). Here we use the the near-field
229estimates of vertical deformation in northwestern Sumatra and
230the Nicobar-Andaman islands collated by~\cite{chlieh07} to validate
231that our crustal deformation model of the 2004 Sumatra--Andaman
232earthquake is producing reasonable results. Note that the geodetic
233data used here is a combination of the vertical deformation that
234happened in the $\sim$10 minutes of the earthquake plus the
235deformation that followed in the days following the earthquake before
236each particular measurement was actually made (typically of order
237days). Therefore some of the observations may not contain the purely
238co-seismic deformation but could include some post-seismic deformation
239as well~\cite{chlieh07}.
241%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.
246%\caption{Near field geodetic measurements used to validate tsunami generation. FIXME: Insert appropriate figure here}
252\label{sec:propagation data}
253Once generated, a tsunami will propagate outwards from the source until
254it encounters the shallow water bordering coastal regions. This period
255of the tsunami evolution is referred to as the propagation stage. The
256height and velocity of the tsunami is dependent on the local
257bathymetry in the regions through which the wave travels and the size
258of the initial wave. This section details the bathymetry data needed
259to model the tsunami propagation and the satellite altimetry transects
260used here to validate open ocean tsunami models.
262\subsubsection{Bathymetry Data}
263The bathymetry data used in this study was derived from the following
266\item a two arc minute grid data set covering the Bay of Bengal,
267  DBDB2, obtained from US Naval Research Labs;
268\item a 3 second arc grid covering the whole of the Andaman Sea based
269  on Thai Navy charts no. 45 and no. 362; and
270\item a one second grid created from the digitised Thai Navy
271  bathymetry chart, no. 358, which covers Patong Bay and the
272  immediately adjacent regions. (FIXME (Ole): How was the grid created from these digitised points?)
274FIXME (Jane): Refs for all these.
276%A number of raw data sets were obtained, analysed and checked for
277%quality and subsequently gridded for easier visualisation and input
278%into the tsunami models.
280These sets were combined via
281interpolation and resampling to produce four nested grids
282which are relatively coarse in the deeper water and
283progressively finer as the distance to
284Patong Beach decreases as shown in Figure~\ref{fig:nested_grids}
286The coarsest
287bathymetry was obtained by interpolating the DBDB2 grid to a 27 second
288arc grid. A subsection of this region was then replaced by nine second
289data which was generated by sub-sampling the three second of arc grid from
290NOAA (FIXME (Jane): This was not mentioned in the dots above).
291A subset of the nine second grid was replaced by the three second
292data. Finally, the one second grid was used to approximate the
293bathymetry in Patong Bay and the immediately adjacent regions. Any
294points that deviated from the general trend near the boundary were
295deleted as a quality check.
297The sub-sampling of larger grids was performed by using {\bf resample},
298a Generic Mapping Tools (\textsc{GMT}) program (\cite{wessel98}). The
299gridding of data was performed using {\bf Intrepid}, a commercial
300geophysical processing package developed by Intrepid Geophysics. The
301gridding scheme employed the nearest neighbour algorithm followed by
302an application of minimum curvature akima spline smoothing.
303See \url{} 
304for details on the Intrepid model.
310\caption{Nested bathymetry grids.}
315\subsubsection{JASON Satellite Altimetry}\label{sec:data_jason}
316During the 26 December 2004 event, the \textsc{jason} satellite tracked from
317north to south and over the equator at 02:55 UTC nearly two hours
318after the earthquake \cite{gower05}. The satellite recorded the sea
319level anomaly compared to the average sea level from its previous five
320passes over the same region in the 20-30 days prior. This data was
321used to validate the propagation stage in Section
323%DB I suggest we combine with model data to reduce the number of figures. The satellite track is shown in Figure~\ref{fig:satelliteTrack}.
328%\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.}
336%\caption{JASON satellite altimetry seal level anomaly. FIXME: should we include figure here with just JASON altimetry.}
341%FIXME: Can we compare the urs model against the TOPEX-poseidon satellite as well? DB No (we don't have the data currently).
344\label{sec:inundation data}
345FIXME (Ole): Technically propagation covers everything up to
346the coastline and inundation everything on-shore.
347This means that ANUGA covers the final part of the propagation and the inundation part. Should we adopt this distiction throughout the paper?
349Inundation refers to the final stages of the evolution of a tsunami and
350covers the propagation of the tsunami in coastal waters and the
351subsequent run-up onto land. This process is typically the most
352difficult of the three stages to model due to thin layers of water
353flowing rapidly over dry land.  Aside from requiring robust solvers
354which can simulate such complex flow patterns, this part of the
355modelling process also requires high resolution and quality elevation
356data which is often not available. In the case of model validation
357high quality field measurements are also required. For the proposed
358benchmark a high resolution bathymetry (FIXME (Ole): Bathymetry ?) and
359topography data set and a high quality inundation survey map from the
360Coordinating Committee Co-ordinating Committee for Geoscience Programmes
361in East and Southeast Asia (CCOP) (\cite{szczucinski06}) was obtained
362to validate model inundation. See also acknowledgements at the end of this paper.
364In this section we also present eye-witness accounts which can be used
365to qualitatively validate tsunami inundation.
367\subsubsection{Topography Data}
368A one second grid was used to approximate the topography in Patong
369Bay. This elevation data was again created from the digitised Thai
370Navy bathymetry chart, no 358.
371FIXME (Ole): I don't think so. The Navy chart is only offshore.
373 A visualisation of the elevation data
374set used in Patong Bay is shown in
375Figure~\ref{fig:patong_bathymetry}. The continuous topography
376(FIXME(Jane): What is meant by this?) is an
377interpolation of known elevation measured at the coloured dots. FIXME ??
382\caption{3D visualisation of the elevation data set used in Patong Bay showing data points, contours, rivers and roads draped over the final model.}
386FIXME (Jane): legend? Were the contours derived from the final dataset?
387This is not the entire mode, only the bay and the beach.
389\subsubsection{Buildings and Other Structures}
390Human-made buildings and structures can significantly affect tsunami
391inundation. The footprint and number of floors of the
392buildings 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).
393The heights of these
394buildings were estimated assuming that each floor has a height of 3 m and they
395were added to the topographic dataset.
397\subsubsection{Inundation Survey}
398Tsunami run-up is the cause of the largest financial and human
399losses, yet run-up data that can be used to validate model run-up
400predictions is scarce. Of the two field benchmarks proposed
402only the Okushiri benchmark facilitates comparison between
403modelled and observed run-up. One of the major strengths of the
404benchmark proposed here is that modelled run-up can be compared to an
405inundation survey which maps the maximum run-up along an entire coastline
406rather than at a series of discrete sites. The survey map is
407shown in Figure~\ref{fig:patongescapemap} and plots the maximum run-up
408of the 2004 Indian Ocean tsunami in Patong city. Refer to Szczucinski et
409al~\cite{szczucinski06} for further details.
414\caption{Tsunami survey mapping the maximum observed inundation at
415  Patong beach courtesy of the CCOP \protect \cite{szczucinski06}.}
421\subsubsection{Eyewitness Accounts}\label{sec:eyewitness data}
422Eyewitness accounts detailed in~\cite{papadopoulos06}
423report that most people at Patong Beach observed an initial retreat of
424the shoreline of more than 100 m followed a few minutes later, by a
425strong wave (crest). Another less powerful wave arrived another five
426or ten minutes later. Eyewitness statements place the arrival time of
427the strong wave between 2 hours and 55 minutes to 3 hours and 5
428minutes after the source rupture (09:55am to 10:05am local time).
430Two videos were sourced\footnote{The footage is
431widely available and can for example be obtained from
433(Comfort Hotel) and
437which include footage of the tsunami in Patong Bay on the day
438of the 2004 Indian Ocean Tsunami. Both videos show an already inundated
439group of buildings. They also show what is to be assumed as the second
440and third waves approaching and further flooding of the buildings and
441street.  The first video is in the very north, filmed from what is
442believed to be the roof of the Novotel Hotel marked ``north'' in Figure
443\ref{fig:gauge_locations}. The second video is in the very south,
444filmed from the second story of a building next door to the Comfort
445Resort near the corner of Ruam Chai St and Thaweewong Road.  This
446location is marked ``south'' in Figure \ref{fig:gauge_locations}.
447Figure~\ref{fig:video_flow} shows stills from this video. Both videos
448were used to estimate flow speeds and inundation depths over time.
456\caption{Four frames from a video where flow rate could be estimated,
457  circle indicates tracked debris, from top left: 0.0 sec, 5.0 s, 7.1
458  s, 7.6 s.}
463Flow rates were estimated using landmarks found in both videos and
464were found to be in the range of 5 to 7 metres per second (+/- 2 m/s)
465in the north and 0.5 to 2 metres per second (+/- 1 m/s) in the south.
466FIXME (Jane): How were these error bounds derived?
467Water depths could also
468be estimated from the videos by the level at which water rose up the
469sides of buildings such as shops. Our estimates are in the order of
4701.5 to 2.0 metres (+/- 0.5 m).
471Fritz ~\cite{fritz06} performed a detailed
472analysis of video frames taken around Banda Aceh and arrived at flow
473speeds in the range of 2 to 5 m/s.
476\subsection{Validation Check-List}
478The data described in this section can be used to construct a
479benchmark to validate all three stages of the evolution of a
480tsunami. In particular we propose that a legitimate tsunami model
481should reproduce the following behaviour:
483 \item reproduce the vertical deformation observed in north-western
484   Sumatra and along the Nicobar--Andaman islands (see
485   Section~\ref{sec:gen_data}),
486 \item reproduce the \textsc{jason} satellite altimetry sea surface
487   anomalies (see Section~\ref{sec:data_jason}),
488 \item reproduce the inundation survey map in Patong city
489   (Figure~\ref{fig:patongescapemap}),
490 \item simulate a leading depression followed by two distinct crests
491   of decreasing magnitude at the beach, and
492 \item predict the water depths and flow speeds, at the locations of
493   the eye-witness videos, that fall within the bounds obtained from
494   the videos.
497Ideally, the model should also be compared to measured timeseries of
498waveheights and velocities but the authors are not aware of the
499availability of such data near Patong Bay.
504\section{Modelling the Event}\label{sec:models}
505Numerous models are currently used to model and predict tsunami
506generation, propagation and run-up~\cite{titov97a,satake95}. Here we
507introduce the three part modelling methodology employed by Geoscience Australia
508to illustrate the utility of the proposed benchmark.
512There are various approaches to modelling the expected crustal
513deformation from an earthquake. Most approaches model the
514earthquake as a dislocation in a linear elastic medium. Here we use
515the method of Wang et al~\cite{wang03}. One of the main advantages
516of their method is that it allows the dislocation to be located in a
517stratified linear elastic half-space with an arbitrary number of
518layers. Other methods (such as those based on Okada's equations) can
519only model the dislocation in a homogeneous elastic half space, or can
520only include a limited number of layers, and thus cannot model the
521effect of the depth dependence of the elasticity of the
522Earth~\cite{wang03}. The original versions of the codes described here
523are available from \url{}. The
524first program, \textsc{edgrn}, calculates elastic Green's function for
525a set of point sources at a regular set of depths out to a specified
526distance. The equations controlling the deformation are solved by
527using a combination of Hankel's transform and Wang et al's
528implementation of the Thomson-Haskell propagator
529algorithm~\cite{wang03}. Once the Green's functions are calculated
530a slightly modified version of \textsc{edcmp}\footnote{For this study,
531we have made minor modifications
532to \textsc{edcmp} in order for it to provide output in a file format
533compatible with the propagation code in the following section. Otherwise it
534is similar to the original code.} is used to calculate the sea
535floor deformation for a specific subfault. This second code
536discretises the subfault into a set of unit sources and sums the
537elastic Green's functions calculated from \textsc{edgrn} for all the
538unit sources on the fault plane in order to calculate the final static
539deformation caused by a two dimensional dislocation along the
540subfault. This step is possible because of the linearity of the
541governing equations.
543In order to calculate the crustal deformation using these codes
544a model that describes the variation in elastic
545properties with depth and a slip model of the earthquake to describe
546the dislocation is required.
547The elastic parameters used for this study are the
548same as those in Table 2 of Burbidge et al~\cite{burbidge08}. For the slip
549model, there are many possible models for the 2004 Andaman--Sumatran
550earthquake to select from
551~\cite{chlieh07,asavanant08,arcas06,grilli07,ioualalen07}. Some are
552determined from various geological surveys of the site. Others solve
553an inverse problem which calibrates the source based upon the tsunami
554wave signal, the seismic signal and/or even the run-up.
555The source
556parameters used here to simulate the 2004 Indian Ocean tsunami were
557taken from the slip model G-M9.15 of Chlieh
558et al~\cite{chlieh07}. This model was created by inversion of wide
559range of geodetic and seismic data. The slip model consists of 686
56020km x 20km subsegments each with a different slip, strike and dip
561angle. The dip subfaults go from $17.5^0$ in the north and $12^0$ in
562the south. Refer to Chlieh et al~\cite{chlieh07} for a detailed
563discussion of this model and its derivation. Note that the geodetic
564data used in the validation was also included by~\cite{chlieh07} in
565the inversion used to find G-M9.15. Thus the validation is not
566completely independent. However, a reasonable validation would still
567show that the crustal deformation and elastic properties model used
568here is at least as valid as the one used by Chlieh
569et al~\cite{chlieh07} and can reproduce the observations just as
573The \textsc{ursga} model described below was used to simulate the
574propagation of the 2004 Indian Ocean tsunami across the open ocean, based on a
575discrete representation of the initial deformation of the sea floor, as
576described in Section~\ref{sec:modelGeneration}. For the models shown
577here, the uplift is assumed to be instantaneous and creates a wave of
578the same size and amplitude as the co-seismic sea floor deformation.
581\textsc{ursga} is a hydrodynamic code that models the propagation of
582the tsunami in deep water using a finite difference method on a staggered grid.
583It solves the depth integrated linear or nonlinear shallow water equations in
584spherical co-ordinates with friction and Coriolis terms. The code is
585based on Satake~\cite{satake95} with significant modifications made by
586the \textsc{urs} corporation, Thio et al~\cite{thio08} and Geoscience
587Australia, Burbidge et al~\cite{burbidge08}.
588The tsunami was propagated via the nested
589grid system described in Section \ref{sec:propagation data} where
590the coarse grids were used in the open ocean and the finest
591resolution grid was employed in the region closest to Patong bay.
592\textsc{Ursga} is not publicly available.
595The utility of the \textsc{ursga} model decreases with water depth
596unless an intricate sequence of nested grids is employed. In
597comparison \textsc{anuga}, described below, is designed to produce
598robust and accurate predictions of on-shore inundation, but is less
599suitable for earthquake source modelling and large study areas because
600it is based on projected spatial coordinates. Consequently, the
601Geoscience Australia tsunami modelling methodology is based on a
602hybrid approach using models like \textsc{ursga} for tsunami
603propagation up to an offshore depth contour, typically 100 m.
604%Specifically we use the \textsc{ursga} model to simulate the
605%propagation of the 2004 Indian Ocean tsunami in the deep ocean, based
606%on a discrete representation of the initial deformation of the sea
607%floor, described in Section~\ref{sec:modelGeneration}.
608The wave signal and the velocity field is then used as a
609time varying boundary condition for
610the \textsc{anuga} inundation simulation.
611% A description of \textsc{anuga} is the following section.
614\textsc{Anuga} is a Free and Open Source hydrodynamic inundation tool that
615solves the conserved form of the depth-integrated nonlinear shallow
616water wave equations using a Finite-Volume scheme on an
617unstructured triangular mesh.
618The scheme, first
619presented by Zoppou and Roberts~\cite{zoppou99}, is a high-resolution
620Godunov-type method that uses the rotational invariance property of
621the shallow water equations to transform the two-dimensional problem
622into local one-dimensional problems. These local Riemann problems are
623then solved using the semi-discrete central-upwind scheme of Kurganov
624et al~\cite{kurganov01} for solving one-dimensional conservation
625equations. The numerical scheme is presented in detail in
626Roberts and Zoppou~\cite{zoppou00,roberts00} and
627Nielsen et al~\cite{nielsen05}. An important capability of the
628finite-volume scheme is that discontinuities in all conserved quantities
629are allowed at every edge in the mesh. This means that the tool is
630well suited to adequately resolving hydraulic jumps, transcritical flows and
631the process of wetting and drying. This means that \textsc{Anuga} 
632is suitable for
633simulating water flow onto a beach or dry land and around structures
634such as buildings. \textsc{Anuga} has been validated against
635%a number of analytical solutions and  FIXME: These have not been published
636the wave tank simulation of the 1993 Okushiri
637Island tsunami~\cite{nielsen05,roberts06}.
638FIXME (Ole): Add reference to Tom Baldock's Dam Break valiadation of ANUGA.
643This section presents a validation of the modelling practice of Geoscience
644Australia against the new proposed benchmarks. The criteria outlined
645in Section~\ref{sec:checkList} are addressed for each of the three stages
646of tsunami evolution.
649The location and magnitude of the sea floor displacement associated
650with the 2004 Sumatra--Andaman tsunami calculated from the G-M9.15
651model of~\cite{chlieh07} is shown in
652Figure~\ref{fig:surface_deformation}. The magnitude of the sea floor
653displacement ranges from about $-3.0$ to $5.0$ metres. The region near
654the fault is predicted to uplift, while that further away from the
655fault subsides. Also shown in Figure~\ref{fig:surface_deformation} are
656the areas that were observed to uplift (arrows pointing up) or subside
657(arrows point down) during and immediately after the earthquake. Most
658of this data comes from uplifted or subsided coral heads. The length of
659the vector increases with the magnitude of the displacement; the length
660corresponding to 1 m of observed motion is shown in the top right
661corner of the figure. As can be seen, the source model detailed in
662Section~\ref{sec:modelGeneration} produces a crustal deformation that
663matches the vertical displacements in the Nicobar-Andaman islands and
664Sumatra very well. Uplifted regions are close to the fault and
665subsided regions are further away. The crosses on
666Figure~\ref{fig:surface_deformation} show estimates of the pivot line
667from the remote sensing data~\cite{chlieh07} and they follow the
668predicted pivot line quite accurately. The average difference between
669the observed motion and the predicted motion (including the pivot line
670points) is only 0.06 m, well below the typical error of the
671observations of between 0.25 and 1.0 m. However, the occasional point
672has quite a large error (over 1 m); for example a couple of
673uplifted/subsided points appear to be on a wrong
674(FIXME (Jane): This is incorrect) side of the predicted
675pivot line~\ref{fig:surface_deformation}. The excellence of the fit is
676not surprising, since the original slip model was chosen
677by~\cite{chlieh07} to fit this (and the seismic data) well.
678This does demonstrate, however, that \textsc{edgrn} and our modified version of
679\textsc{edstat} (FIXME(Jane): This has never been mentioned before)
680can reproduce the correct pattern of vertical
681deformation very well when the slip distribution is well constrained
682and when reasonable values for the elastic properties are used.
688\caption{Location and magnitude of the vertical component of the sea
689  floor displacement associated with the 2004 Indian Ocean tsunami
690  based on the slip model, G-M9.15. The black arrows which point up
691  show areas observed to uplift during and immediately after the
692  earthquake; those pointing down are locations which subsided. The
693  length of the arrow increases with the magnitude of the deformation. The arrow
694  length corresponding to 1 m of deformation is shown in the top right
695  hand corner of the figure. The cross marks show the location of
696  the pivot line (the region between the uplift and subsided region
697  where the uplift is zero) derived from remote sensing
698  (FIXME(Jane): How was that possible?). All the
699  observational data are from the dataset collated
700  by~\cite{chlieh07}.}
706The deformation results described in Section~\ref{sec:modelGeneration}
707were used to provide a profile of the initial ocean surface
708displacement. This wave was used as an initial condition for
709\textsc{ursga} and was propagated throughout the Bay of Bengal. The
710rectangular computational domain of the largest grid extended from
71190$^0$ to 100$^0$ East and 0 to 15$^0$ North and contained
7121335$\times$1996 finite difference points. Inside this grid, a nested
713sequence of grids was used. The grid resolution of the nested grids
714went from 27 arc seconds in the coarsest grid, down to nine arc seconds
715in the second grid, three arc seconds in the third grid and finally one arc
716second in the finest grid near Patong. The computational domain is
717shown in Figure~\ref{fig:computational_domain}.
724\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).}
730Figure \ref{fig:jasonComparison} provides a comparison of the
731\textsc{ursga}-predicted sea surface elevation with the \textsc{jason}
732satellite altimetry data. The \textsc{ursga} model replicates the
733amplitude and timing of the the wave observed at $2.5^0$ South,
734but underestimates the amplitude of the wave further to the south at
735$4^0$ South. In the model, the southern most of these two waves
736appears only as a small bump in the cross section of the model (shown
737in Figure~\ref{fig:jasonComparison}) instead of being a distinct peak
738as can be seen in the satellite data. Also note
739that the \textsc{ursga} model prediction of the ocean surface
740elevation becomes out of phase with the \textsc{jason} 
741data at $3^0$ to $7^0$ North
742latitude. Chlieh et al~\cite{chlieh07} also observed these misfits and
743suggest it is caused by a reflected wave from the Aceh Peninsula that
744is not resolved in the model due to insufficient resolution of the
745computational mesh and bathymetry data. This is also a limitation of
746the model presented here which could be improved by nesting
747grids near Aceh.
752\caption{Comparison of the \textsc{ursga}-predicted surface elevation
753  with the \textsc{jason} satellite altimetry data. The \textsc{ursga} wave
754  heights have been corrected for the time the satellite passed
755  overhead compared to \textsc{jason} sea level anomaly.}
759FIXME (Jane): This graph does not look nice. The legend URS Model should
760be URSGA model.
763After propagating the tsunami in the open ocean using \textsc{ursga},
764the approximated ocean and surface elevation and horisontal flow
765velocities were extracted and used to construct a boundary condition
766for the \textsc{anuga} model. The interface between the \textsc{ursga}
767and \textsc{anuga} models was chosen to roughly follow the 100~m depth
768contour along the west coast of Phuket Island. The computational
769domain is shown in Figure~\ref{fig:computational_domain}.
771The domain was discretised into 386,338 triangles. The resolution of
772the grid was increased in regions inside the bay and on-shore to
773efficiently increase the simulation accuracy for the impact area.
774The grid resolution ranged between a
775maximum triangle area of $1\times 10^5$ m$^2$ near the western ocean
776boundary to $20$ m$^2$ in the small regions surrounding the inundation
777region in Patong Bay. Due to a lack of available data, friction was
778set to a constant throughout the computational domain. For the
779reference simulation, a Manning's coefficient of 0.01 was chosen to
780represent a small resistance to the water flow. See Section
781\ref{sec:friction sensitivity} for details on model sensitivity to
782this parameter.
785The boundary condition at each side of the domain towards the south
786and the north where no information about the incident wave or
787its velocity field is available
788was chosen as a transmissive
789boundary condition, effectively replicating the time dependent wave
790height present just inside the computational domain.
791The velocity field on these boundaries was set
792to zero. Other choices include applying the mean tide value as a
793Dirichlet boundary condition. But experiments as well as the
794result of the verification reported here showed that this approach
795tends to underestimate the tsunami impact due to the tempering of the
796wave near the side boundaries, whereas the transmissive boundary
797condition robustly preserves the wave.
799During the \textsc{anuga} simulation the tide was kept constant at
800$0.80$ m. This value was chosen to correspond to the tidal height
801specified by the Thai Navy tide charts
802(\url{}) at the time the tsunami arrived
803at Patong Bay. Although the tsunami propagated for approximately three
804hours before it reach Patong Bay, the period of time during which the
805wave propagated through the \textsc{anuga} domain is much
806smaller. Consequently the assumption of constant tide height is
809Maximum onshore inundation depth was computed from the model
810throughout the entire Patong Bay region.
811Figure~\ref{fig:inundationcomparison1cm} (left) shows very good
812agreement between the measured and simulated inundation. However
813these results are dependent on the classification used to determine
814whether a region in the numerical simulation was inundated. In
815Figure~\ref{fig:inundationcomparison1cm} (left) a point in the computational
816domain was deemed inundated if at some point in time it was covered by
817at least 1 cm of water. However, the precision of the inundation boundary
818generated by the on-site survey is most likely less than that as it
819was determined by observing water marks and other signs
820left by the receding waters. Consequently the measurement error along
821the inundation boundary of the survey is likely to vary significantly
822and somewhat unpredictably.
823An inundation threshold of 10 cm therefore was selected for inundation
824extents reported in this paper to reflect
825the more likely accuracy of the survey, and subsequently facilitate a more
826appropriate comparison between the modelled and observed inundation
828Figure~\ref{fig:inundationcomparison1cm} (right) shows the simulated
829inundation using a larger threshold of 10 cm.
832The datasets necessary for reproducing the results
833of the inundation stage are available on Sourceforge under the \textsc{anuga}
834project (\url{}).
835At the time of
836writing the direct link is \url{}.
838The scripts required are part of the \textsc{anuga} distribution also
839available from Sourceforge \url{} under
840the validation section.
842An animation of this simulation is available on the \textsc{anuga} website at \url{} or directly from \url{}.
849\caption{Simulated inundation versus observed inundation using an
850inundation threshold of 1cm (left) and 10cm (right).}
855To quantify the agreement between the observed and simulated inundation we
856introduce the measure
858\rho_{in}=\frac{A(I_m\cap I_o)}{A(I_o)}
860representing the ratio $\rho_{in}$ of the observed
861inundation region $I_o$ captured by the model $I_m$. Another useful
862measure is the fraction of the modelled inundation area that falls
863outside the observed inundation area given by the formula
865\rho_{out}=\frac{A(I_m\setminus (I_m\cap I_o))}{A(I_o)}
867These values for the two aforementioned simulations are given in
868Table~\ref{table:inundationAreas}. High value of both $\rho_{in}$ and $\rho_{out}$ indicate
869that the model overestimates the impact whereas low values of both quantities would indicate
870an underestimation. A high value of $\rho_{in}$ combined with a low value of $\rho_{out}$ 
871indicates a good model prediction of the survey.
873Discrepancies between the survey data and the modelled inundation
874include: unknown distribution of surface roughness, inappropriate
875parameterisation of the source model, effect of humans structures on
876flow, as well as uncertainties in the elevation data, effects of
877erosion and deposition by the tsunami event,
878measurement errors in the GPS survey recordings, and
879missing data in the field survey data itself. The impact of some of
880these sources of uncertainties are is investigated in
883\subsection{Eye-witness accounts}
884Figure \ref{fig:gauge_locations} shows four locations where time
885series have been extracted from the model. The two offshore time series
886are shown in Figure \ref{fig:offshore_timeseries} and the two onshore
887timeseries are shown in Figure \ref{fig:onshore_timeseries}. The
888latter coincide with locations where video footage from the event is
889available as described in Section \ref{sec:eyewitness data}.
894\caption{Location of timeseries extracted from the model output.}
904\caption{Time series obtained from the two offshore gauge locations,
9057C and 10C, shown in Figure \protect \ref{fig:gauge_locations}.}
914\caption{Time series obtained from the two onshore locations, North and South,
915shown in Figure \protect \ref{fig:gauge_locations}.}
921The estimated depths and flow rates given in Section
922\ref{sec:eyewitness data} are shown together with the modelled depths
923and flow rates obtained from the model in Table \ref{tab:depth and
924  flow comparisons}. The minimum depths shown in the model are clearly
925lower than expected and an indication that the tsunami model does not
926predict flow dynamics accurately at this level of detail. However,
927this comparison serves to check that depths and speeds predicted are
928within the range of what is expected.
933  \begin{array}{|l|cc|cc|}
934  \hline
935                 & \multicolumn{2}{|c|}{\mbox{Depth [m]}}
936                 & \multicolumn{2}{c|}{\mbox{Flow [m/s]}} \\ 
937                 & \mbox{Observed} & \mbox{Modelled}
938                 & \mbox{Observed} & \mbox{Modelled} \\ \cline{2-5}                 
939    \mbox{North} & 1.5-2 & 1.4 & 5-7 & 0.1 - 3.3 \\
940    \mbox{South} & 1.5-2 & 1.5 & 0.5-2 & 0.2 - 2.6 \\ \hline
941  \end{array}
943\label{tab:depth and flow comparisons}
945FIXME (Jane): We should perhaps look at average data in area surrounding these points
947%can be estimated with landmarks found in
948%satellite imagery and the use of a GIS and were found to be in the
949%range of 5 to 7 metres per second (+/- 2 m/s) in the north and 0.5 to
950%2 metres per second (+/- 1 m/s) in the south.
952Given the uncertainties in both model and observations, there is agreement
953between the values obtained from the videos and the simulations.
955% Our modelled flow rates show
956%maximum values in the order of 0.2 to 2.6 m/s in the south and 0.1 to
957%3.3 m/s for the north as shown in the figures. Water depths could also
958%be estimated from the videos by the level at which water rose up the
959%sides of buildings such as shops. Our estimates are in the order of
960%1.5 to 2.0 metres (+/- 0.5 m). This is in the same range as our
961%modelled maximum depths of 1.4 m in the north and 1.5 m in the south
962%as seen in the figure.
969\section{Sensitivity Analysis}
971This section investigates the effect of different values of Manning's
972friction coefficient, changing waveheight at the 100 m depth contour,
973and the presence and absence of buildings in the elevation dataset on
974model maximum inundation. The reference model is the one reported in
975Figure~\ref{fig:inundationcomparison1cm} (right) with a friction coefficient of 0.01,
976buildings included and the boundary condition produced by the
977\textsc{ursga} model.
981\label{sec:friction sensitivity}
982The first sensitivity study investigated the impact of surface roughness on the
983predicted run-up. According to Schoettle~\cite{schoettle2007}
984appropriate values of Manning's coefficient range from 0.007 to 0.03
985for tsunami propagation over a sandy sea floor and the reference model
986uses a value of 0.01.  To investigate sensitivity to this parameter,
987we simulated the maximum onshore inundation using a Manning's
988coefficient of 0.0003 and 0.03. The resulting inundation maps are
989shown in Figure~\ref{fig:sensitivity_friction} and the maximum flow
990speeds in Figure~\ref{fig:sensitivity_friction_speed}. These figures
991show that the on-shore inundation extent decreases with increasing
992friction and that small perturbations in the friction cause bounded
993changes in the output. This is consistent with the conclusions of
994Synolakis~\cite{synolakis05} et al, who state that the long wavelength of
995tsunami tends to mean that friction is less important in
996comparison to the motion of the wave.
999\subsection{Input Wave Height}\label{sec:waveheightSA}
1000The effect of the wave height used as input to the inundation model
1001\textsc{anuga} was also investigated.
1002Figure~\ref{fig:sensitivity_boundary} indicates that the inundation
1003severity is directly proportional to the boundary waveheight but small
1004perturbations in the input wave height of 10 cm appear to have little
1005effect on the final inundated area. Obviously larger perturbations
1006will have greater impact. However, wave heights in the open ocean are
1007generally well
1008predicted by the generation and propagation models such as
1009\textsc{ursga} as demonstrated in Section \ref{sec:resultsPropagation} 
1010and also in \cite{thomas2009}.
1015\subsection{Buildings and Other Structures}
1016The presence or absence of physical buildings in the elevation model was also
1019shows the inundated area and the associated maximum flow speeds
1020in the presence and absence of buildings. It
1021is apparent that densely built-up areas act as
1022dissipators greatly reducing the inundated area. However, flow speeds
1023tend to increase in passages between buildings.
1029\caption{$\rho_{in}$ and $\rho_{out}$ of the reference simulation and all sensitivity studies.}
1032 & $\rho_{in}$ & $\rho_{out}$ \\ 
1034Reference model & 0.79 & 0.20\\ 
1035Friction = 0.0003 & 0.83 & 0.26 \\ 
1036Friction = 0.03 & 0.67 & 0.09\\ 
1037Boundary wave hight minus 10 cm & 0.77 & 0.17 \\
1038Boundary wave hight plus 10 cm & 0.82 & 0.22 \\
1039No Buildings & 0.94 & 0.44 \\
1048This paper proposes an additional field data benchmark for the
1049verification of tsunami inundation models. Currently, there is a
1050scarcity of appropriate validation datasets due to a lack of well-documented
1051historical tsunami impacts. The benchmark proposed here
1052utilises the uniquely large amount of observational data for model
1053comparison obtained during, and immediately following, the
1054Sumatra--Andaman tsunami of 26 December 2004. Unlike the small
1055number of existing benchmarks, the proposed test validates all three
1056stages of tsunami evolution - generation, propagation and
1057inundation. In an attempt to provide higher visibility and easier
1058accessibility for tsunami benchmark problems, the data used to
1059construct the proposed benchmark is documented and freely available at
1062This study also shows that the tsunami impact modelling methodology
1063adopted is credible and able to predict inundation extents with reasonable
1064accuracy.  An associated aim of this paper was to further validate the
1065hydrodynamic modelling tool \textsc{anuga} which is used to simulate
1066the tsunami inundation. Model predictions
1067matched well the geodetic measurements of the Sumatra--Andaman earthquake,
1068altimetry data from the \textsc{jason}, eye-witness accounts of wave
1069front arrival times and flow speeds and a detailed inundation survey
1070of Patong Bay, Thailand.
1072A simple sensitivity analysis was performed to assess the influence of
1073small changes in friction, wave height at the 100 m depth contour and
1074the presence of buildings and other structures on the model
1075predictions. Of these three, the presence of buildings was shown to
1076have the greatest influence on
1077the simulated inundation extent. The value of friction and small
1078perturbations in the waveheight at the \textsc{anuga} boundary have
1079comparatively little effect on the model results.
1083This project was undertaken at Geoscience Australia and the Department
1084of Mathematics, The Australian National University. The authors would
1085like to thank Niran Chaimanee from the CCOP for providing
1086the post 2004 tsunami survey data, building footprints, aerial
1087photography and the elevation data for Patong city, Prapasri Asawakun
1088from the Suranaree University of Technology and Parida Kuneepong for
1089supporting this work; and Drew Whitehouse from the Australian National
1090University for preparing the animation of the simulated impact.
1095This appendix present the images used to assess the model sensitivities described in
1102\caption{Results from reference model as reported in Section \protect \ref{sec:results},
1103  i.e.\ including buildings and a friction value of 0.01. The seaward boundary condition is as
1104  provided by the \textsc{ursga} model. The left image shows the maximum
1105  modelled depth while the right hand image shows the maximum modelled
1106  flow velocities.}
1117\caption{Model results with wave height at \textsc{anuga} boundary artificially
1118  modified to assess sensitivities. The reference inundation extent is shown in Figure
1119  \protect \ref{fig:reference_model} (left).  The left and right images
1120  show the inundation results if the wave at the \textsc{anuga} boundary
1121  is reduced or increased by 10 cm respectively. The inundation
1122  severity varies in proportion to the boundary waveheight, but the
1123  model results are only slightly sensitive to this parameter for the
1124  range of values tested.}
1128FIXME (Jane): How and why was the +/- 10 cm chosen?
1135\caption{The maximal flow speeds for the same model parameterisations
1136  found in Figure \protect \ref{fig:sensitivity_boundary}. The
1137  reference flow speeds are shown in Figure \protect
1138  \ref{fig:reference_model} (right).}
1147\caption{Model results show the effect of buildings in
1148  the elevation data set.
1149  The left hand image shows the maximum inundation depth results for
1150  a model entirely without buildings.  As expected, the absence of
1151  buildings will increase the inundation extent beyond what was
1152  surveyed. The right hand image shows the corresponding flow speeds in the absence of buildings. 
1153  The reference results are as shown in Figure
1154  \protect \ref{fig:reference_model}.}
1164\caption{Model results for different values of Manning's friction
1165  coefficient shown to assess sensitivities. The reference inundation extent for a
1166  friction value of 0.01 is shown in Figure
1167  \protect \ref{fig:reference_model} (left).  The left and right images
1168  show the inundation results for friction values of 0.0003 and
1169  0.03 respectively. The inundation extent increases for the lower
1170  friction value while the higher slows the flow and decreases the
1171  inundation extent. Ideally, friction should vary across the entire
1172  domain depending on terrain and vegetation, but this is beyond the
1173  scope of this study.}
1182\caption{The maximal flow speeds for the same model parameterisations
1183  found in Figure \protect \ref{fig:sensitivity_friction}. The
1184  reference flow speeds are shown in Figure \protect
1185  \ref{fig:reference_model} (right).}
Note: See TracBrowser for help on using the repository browser.