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