1 | \section{Results}\label{sec:results} |
---|
2 | This section presents a validation of the modelling practice of Geoscience |
---|
3 | Australia against the new proposed benchmarks. The criteria outlined |
---|
4 | in Section~\ref{sec:checkList} are addressed.S |
---|
5 | |
---|
6 | \subsection{Generation}\label{modelGeneration} |
---|
7 | The location and magnitude of the sea floor displacement associated |
---|
8 | with the 2004 Sumatra--Andaman tsunami calculated from the G-M9.15 |
---|
9 | model of~\cite{chlieh07} is shown in |
---|
10 | Figure~\ref{fig:surface_deformation}. |
---|
11 | %The magnitude of the sea floor |
---|
12 | %displacement ranges from about $-3.0$ to $5.0$ metres. The region near |
---|
13 | %the fault is predicted to uplift, while that further away from the |
---|
14 | %fault subsides. Also shown in Figure~\ref{fig:surface_deformation} are |
---|
15 | %the areas that were observed to uplift (arrows pointing up) or subside |
---|
16 | %(arrows point down) during and immediately after the earthquake. Most |
---|
17 | %of this data comes from uplifted or subsided coral heads. The length of |
---|
18 | %the vector increases with the magnitude of the displacement; the length |
---|
19 | %corresponding to 1 m of observed motion is shown in the top right |
---|
20 | %corner of the figure. |
---|
21 | As can be seen, the source model detailed in |
---|
22 | Section~\ref{sec:modelGeneration} produces a crustal deformation that |
---|
23 | matches the vertical displacements in the Nicobar-Andaman islands and |
---|
24 | Sumatra very well. Uplifted regions are close to the fault and |
---|
25 | subsided regions are further away. The crosses on |
---|
26 | Figure~\ref{fig:surface_deformation} show estimates of the pivot line |
---|
27 | from the remote sensing data~\cite{chlieh07} and they follow the |
---|
28 | predicted pivot line quite accurately. The average difference between |
---|
29 | the observed motion and the predicted motion (including the pivot line |
---|
30 | points) is only 0.06 m, well below the typical error of the |
---|
31 | observations of between 0.25 and 1.0 m. |
---|
32 | % |
---|
33 | %However, the occasional point |
---|
34 | %has quite a large error (over 1 m); for example a couple of |
---|
35 | %uplifted/subsided points appear to be on a wrong |
---|
36 | %(FIXME (Jane): This is incorrect) side of the predicted |
---|
37 | %pivot line~\ref{fig:surface_deformation}. |
---|
38 | % |
---|
39 | The excellence of the fit is |
---|
40 | not surprising, since the original slip model was chosen |
---|
41 | by~\cite{chlieh07} to fit the motion and seismic data well. |
---|
42 | Consequently the replication of the generation data has little meaning for |
---|
43 | the model structure presented in Section~\ref{sec:models}. But for |
---|
44 | uncalibrated source models or source models calibrated on other data |
---|
45 | this test of generation would be more meaningful. |
---|
46 | % |
---|
47 | %This does demonstrate, however, that \textsc{edgrn} and our modified version of |
---|
48 | %\textsc{edstat} (FIXME(Jane): This has never been mentioned before) |
---|
49 | %can reproduce the correct pattern of vertical |
---|
50 | %deformation very well when the slip distribution is well constrained |
---|
51 | %and when reasonable values for the elastic properties are used. |
---|
52 | |
---|
53 | |
---|
54 | \begin{figure}[ht] |
---|
55 | \begin{center} |
---|
56 | \includegraphics[width=0.8\textwidth,keepaspectratio=true]{figures/surface_deformation.jpg} |
---|
57 | %\includegraphics[totalheight=0.3\textheight,width=0.8\textwidth]{surface_deformation.jpg} |
---|
58 | \caption{Location and magnitude of the vertical component of the sea |
---|
59 | floor displacement associated with the 2004 Indian Ocean tsunami |
---|
60 | based on the slip model, G-M9.15. The black arrows which point up |
---|
61 | show areas observed to uplift during and immediately after the |
---|
62 | earthquake; those pointing down are locations which subsided. The |
---|
63 | length of the arrow increases with the magnitude of the deformation. The arrow |
---|
64 | length corresponding to 1 m of deformation is shown in the top right |
---|
65 | hand corner of the figure. The cross marks show the location of |
---|
66 | the pivot line (the region between the uplift and subsided region |
---|
67 | where the uplift is zero) derived from remote sensing |
---|
68 | (FIXME(Jane): How was that possible?). All the |
---|
69 | observational data are from the dataset collated |
---|
70 | by~\cite{chlieh07}.} |
---|
71 | \label{fig:surface_deformation} |
---|
72 | \end{center} |
---|
73 | \end{figure} |
---|
74 | |
---|
75 | \subsection{Propagation}\label{sec:resultsPropagation} |
---|
76 | The deformation results described in Section~\ref{sec:modelGeneration} |
---|
77 | were used to provide a profile of the initial ocean surface |
---|
78 | displacement. This wave was used as an initial condition for |
---|
79 | \textsc{ursga} and was propagated throughout the Bay of Bengal. The |
---|
80 | rectangular computational domain of the largest grid extended from |
---|
81 | 90$^0$ to 100$^0$ East and 0 to 15$^0$ North and contained |
---|
82 | 1335$\times$1996 finite difference points. Inside this grid, a nested |
---|
83 | sequence of grids was used. The grid resolution of the nested grids |
---|
84 | went from 27 arc seconds in the coarsest grid, down to nine arc seconds |
---|
85 | in the second grid, three arc seconds in the third grid and finally one arc |
---|
86 | second in the finest grid near Patong. The computational domain is |
---|
87 | shown in Figure~\ref{fig:computational_domain}. |
---|
88 | |
---|
89 | \begin{figure}[ht] |
---|
90 | \begin{center} |
---|
91 | \includegraphics[width=0.8\textwidth,keepaspectratio=true]{figures/extent_of_ANUGA_model.jpg} |
---|
92 | \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).} |
---|
93 | \label{fig:computational_domain} |
---|
94 | \end{center} |
---|
95 | \end{figure} |
---|
96 | |
---|
97 | |
---|
98 | Figure \ref{fig:jasonComparison} provides a comparison of the |
---|
99 | \textsc{ursga}-predicted sea surface elevation with the \textsc{jason} |
---|
100 | satellite altimetry data. The \textsc{ursga} model replicates the |
---|
101 | amplitude and timing of the the wave observed at $2.5^0$ South, |
---|
102 | but underestimates the amplitude of the wave further to the south at |
---|
103 | $4^0$ South. In the model, the southern most of these two waves |
---|
104 | appears only as a small bump in the cross section of the model (shown |
---|
105 | in Figure~\ref{fig:jasonComparison}) instead of being a distinct peak |
---|
106 | as can be seen in the satellite data. Also note |
---|
107 | that the \textsc{ursga} model prediction of the ocean surface |
---|
108 | elevation becomes out of phase with the \textsc{jason} |
---|
109 | data at $3^0$ to $7^0$ North |
---|
110 | latitude. Chlieh et al~\cite{chlieh07} also observed these misfits and |
---|
111 | suggest it is caused by a reflected wave from the Aceh Peninsula that |
---|
112 | is not resolved in the model due to insufficient resolution of the |
---|
113 | computational mesh and bathymetry data. This is also a limitation of |
---|
114 | the model presented here which could be improved by nesting |
---|
115 | grids near Aceh. |
---|
116 | |
---|
117 | \begin{figure}[ht] |
---|
118 | \begin{center} |
---|
119 | \includegraphics[width=\textwidth,keepaspectratio=true]{figures/jasonComparison.jpg} |
---|
120 | \caption{Comparison of the \textsc{ursga}-predicted surface elevation |
---|
121 | with the \textsc{jason} satellite altimetry data. The \textsc{ursga} wave |
---|
122 | heights have been corrected for the time the satellite passed |
---|
123 | overhead compared to \textsc{jason} sea level anomaly.} |
---|
124 | \label{fig:jasonComparison} |
---|
125 | \end{center} |
---|
126 | \end{figure} |
---|
127 | FIXME (Jane): This graph does not look nice. |
---|
128 | |
---|
129 | After propagating the tsunami in the open ocean using \textsc{ursga}, |
---|
130 | the approximated ocean and surface elevation and horisontal flow |
---|
131 | velocities were extracted and used to construct a boundary condition |
---|
132 | for the \textsc{anuga} model. The interface between the \textsc{ursga} |
---|
133 | and \textsc{anuga} models was chosen to roughly follow the 100~m depth |
---|
134 | contour along the west coast of Phuket Island. The computational |
---|
135 | domain is shown in Figure~\ref{fig:computational_domain}. |
---|
136 | |
---|
137 | The domain was discretised into 386,338 triangles. The resolution of |
---|
138 | the grid was increased in regions inside the bay and on-shore to |
---|
139 | efficiently increase the simulation accuracy for the impact area. |
---|
140 | The grid resolution ranged between a |
---|
141 | maximum triangle area of $1\times 10^5$ m$^2$ |
---|
142 | (corresponding to approximately 440 m between mesh points) |
---|
143 | near the western ocean |
---|
144 | boundary to $20$ m$^2$ (corresponding to |
---|
145 | approximately 6 m between mesh points) |
---|
146 | in the small regions surrounding the inundation |
---|
147 | region in Patong Bay. The coarse resolution was chosen to be |
---|
148 | commensurate with the model output from the \textsc{ursga} model (FIXME (Ole): Richard says that the ursga model used all four grids which would mean that the |
---|
149 | resolution at the ANUGA boundary was 1 second or about 30m. |
---|
150 | This is not consistent with my memory and certainly not with us choosing a |
---|
151 | resolution of 440m. John, do you remember what the spacing was between the |
---|
152 | URSGA points? Did we weed them out or did we take them as they were?) |
---|
153 | |
---|
154 | while the latter was chosen to match the available resolution of topographic |
---|
155 | data and building data in Patong city. |
---|
156 | Due to a lack of available roughness data, friction was |
---|
157 | set to a constant throughout the computational domain. For the |
---|
158 | reference simulation, a Manning's coefficient of 0.01 was chosen to |
---|
159 | represent a small resistance to the water flow. See Section |
---|
160 | \ref{sec:friction sensitivity} for details on model sensitivity to |
---|
161 | this parameter. |
---|
162 | |
---|
163 | |
---|
164 | The boundary condition at each side of the domain towards the south |
---|
165 | and the north where no information about the incident wave or |
---|
166 | its velocity field is available |
---|
167 | was chosen as a transmissive |
---|
168 | boundary condition, effectively replicating the time dependent wave |
---|
169 | height present just inside the computational domain. |
---|
170 | The velocity field on these boundaries was set |
---|
171 | to zero. Other choices include applying the mean tide value as a |
---|
172 | Dirichlet boundary condition. But experiments as well as the |
---|
173 | result of the verification reported here showed that this approach |
---|
174 | tends to underestimate the tsunami impact due to the tempering of the |
---|
175 | wave near the side boundaries, whereas the transmissive boundary |
---|
176 | condition robustly preserves the wave. |
---|
177 | |
---|
178 | During the \textsc{anuga} simulation the tide was kept constant at |
---|
179 | $0.80$ m. This value was chosen to correspond to the tidal height |
---|
180 | specified by the Thai Navy tide charts |
---|
181 | (\url{http://www.navy.mi.th/hydro/}) at the time the tsunami arrived |
---|
182 | at Patong Bay. Although the tsunami propagated for approximately three |
---|
183 | hours before it reach Patong Bay, the period of time during which the |
---|
184 | wave propagated through the \textsc{anuga} domain is much |
---|
185 | smaller. Consequently the assumption of constant tide height is |
---|
186 | reasonable. |
---|
187 | |
---|
188 | \subsection{Inundation}\label{sec:inundation results} |
---|
189 | The \textsc{anuga} simulation described in the previous section and used to |
---|
190 | model shallow water propgation also predicts |
---|
191 | inundation. Maximum onshore inundation depth was computed from the model |
---|
192 | throughout the entire Patong Bay region and used to generate |
---|
193 | a measure of the inundated area. |
---|
194 | Figure~\ref{fig:inundationcomparison1cm} (left) shows very good |
---|
195 | agreement between the measured and simulated inundation. However |
---|
196 | these results are dependent on the classification used to determine |
---|
197 | whether a region in the numerical simulation was inundated. In |
---|
198 | Figure~\ref{fig:inundationcomparison1cm} (left) a point in the computational |
---|
199 | domain was deemed inundated if at some point in time it was covered by |
---|
200 | at least 1 cm of water. However, the precision of the inundation boundary |
---|
201 | generated by the on-site survey is most likely less than that as it |
---|
202 | was determined by observing water marks and other signs |
---|
203 | left by the receding waters. Consequently the measurement error along |
---|
204 | the inundation boundary of the survey is likely to vary significantly |
---|
205 | and somewhat unpredictably. |
---|
206 | An inundation threshold of 10 cm therefore was selected for inundation |
---|
207 | extents reported in this paper to reflect |
---|
208 | the more likely accuracy of the survey, and subsequently facilitate a more |
---|
209 | appropriate comparison between the modelled and observed inundation |
---|
210 | area. |
---|
211 | Figure~\ref{fig:inundationcomparison1cm} (right) shows the simulated |
---|
212 | inundation using a larger threshold of 10 cm. |
---|
213 | |
---|
214 | |
---|
215 | The datasets necessary for reproducing the results |
---|
216 | of the inundation stage are available on Sourceforge under the \textsc{anuga} |
---|
217 | project (\url{http://sourceforge.net/projects/anuga}). |
---|
218 | At the time of |
---|
219 | writing the direct link is \url{http://tinyurl.com/patong2004-data}. |
---|
220 | %%\url{http://sourceforge.net/project/showfiles.php?group_id=172848&package_id=319323&release_id=677531}. |
---|
221 | The scripts required are part of the \textsc{anuga} distribution also |
---|
222 | available from Sourceforge \url{http://sourceforge.net/projects/anuga} under |
---|
223 | the validation section. |
---|
224 | |
---|
225 | 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}. |
---|
226 | %\url{https://datamining.anu.edu.au/anuga/attachment/wiki/AnugaPublications/patong_2004_indian_ocean_tsunami_ANUGA_animation.mov}. |
---|
227 | |
---|
228 | \begin{figure}[ht] |
---|
229 | \begin{center} |
---|
230 | \includegraphics[width=\textwidth,keepaspectratio=true]{figures/threshold.jpg} |
---|
231 | \caption{Simulated inundation versus observed inundation using an |
---|
232 | inundation threshold of 1 cm (left) and 10 cm (right).} |
---|
233 | \label{fig:inundationcomparison1cm} |
---|
234 | \end{center} |
---|
235 | \end{figure} |
---|
236 | |
---|
237 | To quantify the agreement between the observed and simulated inundation we |
---|
238 | introduce the measure |
---|
239 | \begin{equation} |
---|
240 | \rho_{in}=\frac{A(I_m\cap I_o)}{A(I_o)} |
---|
241 | \end{equation} |
---|
242 | representing the ratio of the area of the observed |
---|
243 | inundation region $I_o$ and the area of the observed inundation region |
---|
244 | captured by the model $I_m$. Another useful |
---|
245 | measure is the fraction of the modelled inundation area that falls |
---|
246 | outside the observed inundation area given by the formula |
---|
247 | \begin{equation} |
---|
248 | \rho_{out}=\frac{A(I_m\setminus (I_m\cap I_o))}{A(I_o)} |
---|
249 | \end{equation} |
---|
250 | These values for the two aforementioned simulations are given in |
---|
251 | Table~\ref{table:inundationAreas}. High value of both $\rho_{in}$ and $\rho_{out}$ indicate |
---|
252 | that the model overestimates the impact whereas low values of both quantities would indicate |
---|
253 | an underestimation. A high value of $\rho_{in}$ combined with a low value of $\rho_{out}$ |
---|
254 | indicates a good model prediction of the survey. |
---|
255 | |
---|
256 | Discrepancies between the survey data and the modelled inundation |
---|
257 | include: unknown distribution of surface roughness, inappropriate |
---|
258 | parameterisation of the source model, effect of humans structures on |
---|
259 | flow, as well as uncertainties in the elevation data, effects of |
---|
260 | erosion and deposition by the tsunami event, |
---|
261 | measurement errors in the GPS survey recordings, and |
---|
262 | missing data in the field survey data itself. The impact of some of |
---|
263 | these sources of uncertainties are is investigated in |
---|
264 | Section~\ref{sec:sensitivity}. |
---|
265 | |
---|
266 | \subsection{Eye-witness accounts} |
---|
267 | |
---|
268 | \subsubsection{Arrival time} |
---|
269 | The arrival time of the first wave took place between 9:55 and 10:55 as described in |
---|
270 | Section~\ref{sec:eyewitness data}. The modelled arrival time at the beach is 10:01 |
---|
271 | as can be verified from the animation provided in \label{sec:inundation results}. |
---|
272 | Subsequent waves of variable magnitude appear over the next two hours |
---|
273 | approximately 20-30 minutes apart. |
---|
274 | % 10:01, 10:19, 10:46, 11:13, 11:43 |
---|
275 | The first arrival and overall dynamic behaviour is therefor reasonably consistent with the |
---|
276 | eye-witness accounts. |
---|
277 | |
---|
278 | \subsubsection{Observed wave dynamics} |
---|
279 | Figure \ref{fig:gauge_locations} shows four locations where time |
---|
280 | series have been extracted from the model. The two offshore time series |
---|
281 | are shown in Figure \ref{fig:offshore_timeseries} and the two onshore |
---|
282 | timeseries are shown in Figure \ref{fig:onshore_timeseries}. The |
---|
283 | latter coincide with locations where video footage from the event is |
---|
284 | available as described in Section \ref{sec:eyewitness data}. |
---|
285 | |
---|
286 | \begin{figure}[ht] |
---|
287 | \begin{center} |
---|
288 | \includegraphics[width=0.8\textwidth,keepaspectratio=true]{figures/gauge_bay_depth.jpg} |
---|
289 | \includegraphics[width=0.8\textwidth,keepaspectratio=true]{figures/gauge_bay_speed.jpg} |
---|
290 | \caption{Time series obtained from the two offshore gauge locations, |
---|
291 | 7C and 10C, shown in Figure \protect \ref{fig:gauge_locations}.} |
---|
292 | \end{center} |
---|
293 | \label{fig:offshore_timeseries} |
---|
294 | \end{figure} |
---|
295 | |
---|
296 | \begin{figure}[ht] |
---|
297 | \begin{center} |
---|
298 | \includegraphics[width=\textwidth,keepaspectratio=true]{figures/gauges_hotels_depths.jpg} |
---|
299 | \includegraphics[width=\textwidth,keepaspectratio=true]{figures/gauges_hotels_speed.jpg} |
---|
300 | \caption{Time series obtained from the two onshore locations, North and South, |
---|
301 | shown in Figure \protect \ref{fig:gauge_locations}.} |
---|
302 | \end{center} |
---|
303 | \label{fig:onshore_timeseries} |
---|
304 | \end{figure} |
---|
305 | |
---|
306 | |
---|
307 | The estimated depths and flow rates given in Section |
---|
308 | \ref{sec:eyewitness data} are shown together with the modelled depths |
---|
309 | and flow rates obtained from the model in |
---|
310 | Table \ref{tab:depth and flow comparisons}. |
---|
311 | The predicted maximum depths and speeds are all of the same order |
---|
312 | of what was observed. However, unlike the real event, |
---|
313 | the model estimates complete withdrawal of the water between waves at the |
---|
314 | chosen locations and shows that the model must be used with caution at this |
---|
315 | level of detail. |
---|
316 | Nonetheless, this comparison serves to check that depths and speeds |
---|
317 | predicted are within the range of what is expected. |
---|
318 | |
---|
319 | |
---|
320 | \begin{table} |
---|
321 | \[ |
---|
322 | \begin{array}{|l|cc|cc|} |
---|
323 | \hline |
---|
324 | & \multicolumn{2}{|c|}{\mbox{Depth [m]}} |
---|
325 | & \multicolumn{2}{c|}{\mbox{Flow [m/s]}} \\ |
---|
326 | & \mbox{Observed} & \mbox{Modelled} |
---|
327 | & \mbox{Observed} & \mbox{Modelled} \\ \cline{2-5} |
---|
328 | \mbox{North} & 1.5-2 & 1.4 & 5-7 & 0.1 - 3.3 \\ |
---|
329 | \mbox{South} & 1.5-2 & 1.5 & 0.5-2 & 0.2 - 2.6 \\ \hline |
---|
330 | \end{array} |
---|
331 | \] |
---|
332 | \label{tab:depth and flow comparisons} |
---|
333 | \end{table} |
---|
334 | |
---|
335 | %can be estimated with landmarks found in |
---|
336 | %satellite imagery and the use of a GIS and were found to be in the |
---|
337 | %range of 5 to 7 metres per second (+/- 2 m/s) in the north and 0.5 to |
---|
338 | %2 metres per second (+/- 1 m/s) in the south. |
---|
339 | |
---|
340 | Given the uncertainties in both model and observations, there is agreement |
---|
341 | between the values obtained from the videos and the simulations. |
---|
342 | |
---|
343 | % Our modelled flow rates show |
---|
344 | %maximum values in the order of 0.2 to 2.6 m/s in the south and 0.1 to |
---|
345 | %3.3 m/s for the north as shown in the figures. Water depths could also |
---|
346 | %be estimated from the videos by the level at which water rose up the |
---|
347 | %sides of buildings such as shops. Our estimates are in the order of |
---|
348 | %1.5 to 2.0 metres (+/- 0.5 m). This is in the same range as our |
---|
349 | %modelled maximum depths of 1.4 m in the north and 1.5 m in the south |
---|
350 | %as seen in the figure. |
---|
351 | |
---|