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

Last change on this file since 7465 was 7465, checked in by ole, 11 years ago

Minor cleanup

File size: 16.8 KB
Line 
1\section{Results}\label{sec:results}
2This section presents a validation of the modelling practice of Geoscience
3Australia against the new proposed benchmarks. The criteria outlined
4in Section~\ref{sec:checkList} are addressed.S
5
6\subsection{Generation}\label{modelGeneration}
7The location and magnitude of the sea floor displacement associated
8with the 2004 Sumatra--Andaman tsunami calculated from the G-M9.15
9model of~\cite{chlieh07} is shown in
10Figure~\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.
21As can be seen, the source model detailed in
22Section~\ref{sec:modelGeneration} produces a crustal deformation that
23matches the vertical displacements in the Nicobar-Andaman islands and
24Sumatra very well. Uplifted regions are close to the fault and
25subsided regions are further away. The crosses on
26Figure~\ref{fig:surface_deformation} show estimates of the pivot line
27from the remote sensing data~\cite{chlieh07} and they follow the
28predicted pivot line quite accurately. The average difference between
29the observed motion and the predicted motion (including the pivot line
30points) is only 0.06 m, well below the typical error of the
31observations 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%
39The excellence of the fit is
40not surprising, since the original slip model was chosen
41by~\cite{chlieh07} to fit the motion and seismic data well.
42Consequently the replication of the generation data has little meaning for
43the model structure presented in Section~\ref{sec:models}. But for
44uncalibrated source models or source models calibrated on other data
45this 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}
76The deformation results described in Section~\ref{sec:modelGeneration}
77were used to provide a profile of the initial ocean surface
78displacement. This wave was used as an initial condition for
79\textsc{ursga} and was propagated throughout the Bay of Bengal. The
80rectangular computational domain of the largest grid extended from
8190$^0$ to 100$^0$ East and 0 to 15$^0$ North and contained
821335$\times$1996 finite difference points. Inside this grid, a nested
83sequence of grids was used. The grid resolution of the nested grids
84went from 27 arc seconds in the coarsest grid, down to nine arc seconds
85in the second grid, three arc seconds in the third grid and finally one arc
86second in the finest grid near Patong. The computational domain is
87shown 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
98Figure \ref{fig:jasonComparison} provides a comparison of the
99\textsc{ursga}-predicted sea surface elevation with the \textsc{jason}
100satellite altimetry data. The \textsc{ursga} model replicates the
101amplitude and timing of the the wave observed at $2.5^0$ South,
102but 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
104appears only as a small bump in the cross section of the model (shown
105in Figure~\ref{fig:jasonComparison}) instead of being a distinct peak
106as can be seen in the satellite data. Also note
107that the \textsc{ursga} model prediction of the ocean surface
108elevation becomes out of phase with the \textsc{jason} 
109data at $3^0$ to $7^0$ North
110latitude. Chlieh et al~\cite{chlieh07} also observed these misfits and
111suggest it is caused by a reflected wave from the Aceh Peninsula that
112is not resolved in the model due to insufficient resolution of the
113computational mesh and bathymetry data. This is also a limitation of
114the model presented here which could be improved by nesting
115grids 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}
127FIXME (Jane): This graph does not look nice.
128
129After propagating the tsunami in the open ocean using \textsc{ursga},
130the approximated ocean and surface elevation and horisontal flow
131velocities were extracted and used to construct a boundary condition
132for the \textsc{anuga} model. The interface between the \textsc{ursga}
133and \textsc{anuga} models was chosen to roughly follow the 100~m depth
134contour along the west coast of Phuket Island. The computational
135domain is shown in Figure~\ref{fig:computational_domain}.
136
137The domain was discretised into 386,338 triangles. The resolution of
138the grid was increased in regions inside the bay and on-shore to
139efficiently increase the simulation accuracy for the impact area.
140The grid resolution ranged between a
141maximum triangle area of $1\times 10^5$ m$^2$ 
142(corresponding to approximately 440 m between mesh points)
143near the western ocean
144boundary to $20$ m$^2$ (corresponding to
145approximately 6 m between mesh points)
146in the small regions surrounding the inundation
147region in Patong Bay. The coarse resolution was chosen to be
148commensurate 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
149resolution at the ANUGA boundary was 1 second or about 30m.
150This is not consistent with my memory and certainly not with us choosing a
151resolution of 440m. John, do you remember what the spacing was between the
152URSGA points? Did we weed them out or did we take them as they were?)
153
154while the latter was chosen to match the available resolution of topographic
155data and building data in Patong city.
156Due to a lack of available roughness data, friction was
157set to a constant throughout the computational domain. For the
158reference simulation, a Manning's coefficient of 0.01 was chosen to
159represent a small resistance to the water flow. See Section
160\ref{sec:friction sensitivity} for details on model sensitivity to
161this parameter.
162
163
164The boundary condition at each side of the domain towards the south
165and the north where no information about the incident wave or
166its velocity field is available
167was chosen as a transmissive
168boundary condition, effectively replicating the time dependent wave
169height present just inside the computational domain.
170The velocity field on these boundaries was set
171to zero. Other choices include applying the mean tide value as a
172Dirichlet boundary condition. But experiments as well as the
173result of the verification reported here showed that this approach
174tends to underestimate the tsunami impact due to the tempering of the
175wave near the side boundaries, whereas the transmissive boundary
176condition robustly preserves the wave.
177
178During the \textsc{anuga} simulation the tide was kept constant at
179$0.80$ m. This value was chosen to correspond to the tidal height
180specified by the Thai Navy tide charts
181(\url{http://www.navy.mi.th/hydro/}) at the time the tsunami arrived
182at Patong Bay. Although the tsunami propagated for approximately three
183hours before it reach Patong Bay, the period of time during which the
184wave propagated through the \textsc{anuga} domain is much
185smaller. Consequently the assumption of constant tide height is
186reasonable.
187
188\subsection{Inundation}
189The \textsc{anuga} simulation described in the previous section and used to
190 model shallow water propgation also predicts
191inundation. Maximum onshore inundation depth was computed from the model
192throughout the entire Patong Bay region and used to generate
193a measure of the inundated area.
194Figure~\ref{fig:inundationcomparison1cm} (left) shows very good
195agreement between the measured and simulated inundation. However
196these results are dependent on the classification used to determine
197whether a region in the numerical simulation was inundated. In
198Figure~\ref{fig:inundationcomparison1cm} (left) a point in the computational
199domain was deemed inundated if at some point in time it was covered by
200at least 1 cm of water. However, the precision of the inundation boundary
201generated by the on-site survey is most likely less than that as it
202was determined by observing water marks and other signs
203left by the receding waters. Consequently the measurement error along
204the inundation boundary of the survey is likely to vary significantly
205and somewhat unpredictably.
206An inundation threshold of 10 cm therefore was selected for inundation
207extents reported in this paper to reflect
208the more likely accuracy of the survey, and subsequently facilitate a more
209appropriate comparison between the modelled and observed inundation
210area.
211Figure~\ref{fig:inundationcomparison1cm} (right) shows the simulated
212inundation using a larger threshold of 10 cm.
213
214
215The datasets necessary for reproducing the results
216of the inundation stage are available on Sourceforge under the \textsc{anuga}
217project (\url{http://sourceforge.net/projects/anuga}).
218At the time of
219writing 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}.
221The scripts required are part of the \textsc{anuga} distribution also
222available from Sourceforge \url{http://sourceforge.net/projects/anuga} under
223the validation section.
224
225An 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
232inundation threshold of 1cm (left) and 10cm (right).}
233\label{fig:inundationcomparison1cm}
234\end{center}
235\end{figure}
236
237To quantify the agreement between the observed and simulated inundation we
238introduce the measure
239\begin{equation}
240\rho_{in}=\frac{A(I_m\cap I_o)}{A(I_o)}
241\end{equation}
242representing the ratio of the area of the observed
243inundation region $I_o$ and the area of the observed inundation region
244captured by the model $I_m$. Another useful
245measure is the fraction of the modelled inundation area that falls
246outside 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}
250These values for the two aforementioned simulations are given in
251Table~\ref{table:inundationAreas}. High value of both $\rho_{in}$ and $\rho_{out}$ indicate
252that the model overestimates the impact whereas low values of both quantities would indicate
253an underestimation. A high value of $\rho_{in}$ combined with a low value of $\rho_{out}$ 
254indicates a good model prediction of the survey.
255
256Discrepancies between the survey data and the modelled inundation
257include: unknown distribution of surface roughness, inappropriate
258parameterisation of the source model, effect of humans structures on
259flow, as well as uncertainties in the elevation data, effects of
260erosion and deposition by the tsunami event,
261measurement errors in the GPS survey recordings, and
262missing data in the field survey data itself. The impact of some of
263these sources of uncertainties are is investigated in
264Section~\ref{sec:sensitivity}.
265
266\subsection{Eye-witness accounts}
267Figure \ref{fig:gauge_locations} shows four locations where time
268series have been extracted from the model. The two offshore time series
269are shown in Figure \ref{fig:offshore_timeseries} and the two onshore
270timeseries are shown in Figure \ref{fig:onshore_timeseries}. The
271latter coincide with locations where video footage from the event is
272available as described in Section \ref{sec:eyewitness data}.
273
274\begin{figure}[ht]
275\begin{center}
276\includegraphics[width=0.8\textwidth,keepaspectratio=true]{figures/gauge_bay_depth.jpg}
277\includegraphics[width=0.8\textwidth,keepaspectratio=true]{figures/gauge_bay_speed.jpg}
278\caption{Time series obtained from the two offshore gauge locations,
2797C and 10C, shown in Figure \protect \ref{fig:gauge_locations}.}
280\end{center}
281\label{fig:offshore_timeseries}
282\end{figure}
283
284\begin{figure}[ht]
285\begin{center}
286\includegraphics[width=\textwidth,keepaspectratio=true]{figures/gauges_hotels_depths.jpg}
287\includegraphics[width=\textwidth,keepaspectratio=true]{figures/gauges_hotels_speed.jpg}
288\caption{Time series obtained from the two onshore locations, North and South,
289shown in Figure \protect \ref{fig:gauge_locations}.}
290\end{center}
291\label{fig:onshore_timeseries}
292\end{figure}
293
294
295The estimated depths and flow rates given in Section
296\ref{sec:eyewitness data} are shown together with the modelled depths
297and flow rates obtained from the model in
298Table \ref{tab:depth and flow comparisons}.
299The predicted maximum depths and speeds are all of the same order
300of what was observed. However, unlike the real event,
301the model estimates complete withdrawal of the water between waves at the
302chosen locations and shows that the model must be used with caution at this
303level of detail.
304Nonetheless, this comparison serves to check that depths and speeds
305predicted are within the range of what is expected.
306
307
308\begin{table}
309\[
310  \begin{array}{|l|cc|cc|}
311  \hline
312                 & \multicolumn{2}{|c|}{\mbox{Depth [m]}}
313                 & \multicolumn{2}{c|}{\mbox{Flow [m/s]}} \\ 
314                 & \mbox{Observed} & \mbox{Modelled}
315                 & \mbox{Observed} & \mbox{Modelled} \\ \cline{2-5}                 
316    \mbox{North} & 1.5-2 & 1.4 & 5-7 & 0.1 - 3.3 \\
317    \mbox{South} & 1.5-2 & 1.5 & 0.5-2 & 0.2 - 2.6 \\ \hline
318  \end{array}
319\]
320\label{tab:depth and flow comparisons}
321\end{table} 
322FIXME (Jane): We should perhaps look at average data in area surrounding these points
323
324%can be estimated with landmarks found in
325%satellite imagery and the use of a GIS and were found to be in the
326%range of 5 to 7 metres per second (+/- 2 m/s) in the north and 0.5 to
327%2 metres per second (+/- 1 m/s) in the south.
328
329Given the uncertainties in both model and observations, there is agreement
330between the values obtained from the videos and the simulations.
331
332% Our modelled flow rates show
333%maximum values in the order of 0.2 to 2.6 m/s in the south and 0.1 to
334%3.3 m/s for the north as shown in the figures. Water depths could also
335%be estimated from the videos by the level at which water rose up the
336%sides of buildings such as shops. Our estimates are in the order of
337%1.5 to 2.0 metres (+/- 0.5 m). This is in the same range as our
338%modelled maximum depths of 1.4 m in the north and 1.5 m in the south
339%as seen in the figure.
340
Note: See TracBrowser for help on using the repository browser.