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

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

Added info from conversations with Richard.

File size: 16.9 KB
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
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
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.
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}.
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.
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.
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}.}
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}.
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).}
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.
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.}
127FIXME (Jane): This graph does not look nice.
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}.
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?)
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.
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.
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{}) 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
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
211Figure~\ref{fig:inundationcomparison1cm} (right) shows the simulated
212inundation using a larger threshold of 10 cm.
215The datasets necessary for reproducing the results
216of the inundation stage are available on Sourceforge under the \textsc{anuga}
217project (\url{}).
218At the time of
219writing the direct link is \url{}.
221The scripts required are part of the \textsc{anuga} distribution also
222available from Sourceforge \url{} under
223the validation section.
225An animation of this simulation is available on the \textsc{anuga} website at \url{} or directly from \url{}.
233\caption{Simulated inundation versus observed inundation using an
234inundation threshold of 1cm (left) and 10cm (right).}
239To quantify the agreement between the observed and simulated inundation we
240introduce the measure
242\rho_{in}=\frac{A(I_m\cap I_o)}{A(I_o)}
244representing the ratio of the area of the observed
245inundation region $I_o$ and the area of the observed inundation region
246captured by the model $I_m$. Another useful
247measure is the fraction of the modelled inundation area that falls
248outside the observed inundation area given by the formula
250\rho_{out}=\frac{A(I_m\setminus (I_m\cap I_o))}{A(I_o)}
252These values for the two aforementioned simulations are given in
253Table~\ref{table:inundationAreas}. High value of both $\rho_{in}$ and $\rho_{out}$ indicate
254that the model overestimates the impact whereas low values of both quantities would indicate
255an underestimation. A high value of $\rho_{in}$ combined with a low value of $\rho_{out}$ 
256indicates a good model prediction of the survey.
258Discrepancies between the survey data and the modelled inundation
259include: unknown distribution of surface roughness, inappropriate
260parameterisation of the source model, effect of humans structures on
261flow, as well as uncertainties in the elevation data, effects of
262erosion and deposition by the tsunami event,
263measurement errors in the GPS survey recordings, and
264missing data in the field survey data itself. The impact of some of
265these sources of uncertainties are is investigated in
268\subsection{Eye-witness accounts}
269Figure \ref{fig:gauge_locations} shows four locations where time
270series have been extracted from the model. The two offshore time series
271are shown in Figure \ref{fig:offshore_timeseries} and the two onshore
272timeseries are shown in Figure \ref{fig:onshore_timeseries}. The
273latter coincide with locations where video footage from the event is
274available as described in Section \ref{sec:eyewitness data}.
280\caption{Time series obtained from the two offshore gauge locations,
2817C and 10C, shown in Figure \protect \ref{fig:gauge_locations}.}
290\caption{Time series obtained from the two onshore locations, North and South,
291shown in Figure \protect \ref{fig:gauge_locations}.}
297The estimated depths and flow rates given in Section
298\ref{sec:eyewitness data} are shown together with the modelled depths
299and flow rates obtained from the model in
300Table \ref{tab:depth and flow comparisons}.
301The predicted maximum depths and speeds are all of the same order
302of what was observed. However, unlike the real event,
303the model estimates complete withdrawal of the water between waves at the
304chosen locations and shows that the model must be used with caution at this
305level of detail.
306Nonetheless, this comparison serves to check that depths and speeds
307predicted are within the range of what is expected.
312  \begin{array}{|l|cc|cc|}
313  \hline
314                 & \multicolumn{2}{|c|}{\mbox{Depth [m]}}
315                 & \multicolumn{2}{c|}{\mbox{Flow [m/s]}} \\ 
316                 & \mbox{Observed} & \mbox{Modelled}
317                 & \mbox{Observed} & \mbox{Modelled} \\ \cline{2-5}                 
318    \mbox{North} & 1.5-2 & 1.4 & 5-7 & 0.1 - 3.3 \\
319    \mbox{South} & 1.5-2 & 1.5 & 0.5-2 & 0.2 - 2.6 \\ \hline
320  \end{array}
322\label{tab:depth and flow comparisons}
324FIXME (Jane): We should perhaps look at average data in area surrounding these points
326%can be estimated with landmarks found in
327%satellite imagery and the use of a GIS and were found to be in the
328%range of 5 to 7 metres per second (+/- 2 m/s) in the north and 0.5 to
329%2 metres per second (+/- 1 m/s) in the south.
331Given the uncertainties in both model and observations, there is agreement
332between the values obtained from the videos and the simulations.
334% Our modelled flow rates show
335%maximum values in the order of 0.2 to 2.6 m/s in the south and 0.1 to
336%3.3 m/s for the north as shown in the figures. Water depths could also
337%be estimated from the videos by the level at which water rose up the
338%sides of buildings such as shops. Our estimates are in the order of
339%1.5 to 2.0 metres (+/- 0.5 m). This is in the same range as our
340%modelled maximum depths of 1.4 m in the north and 1.5 m in the south
341%as seen in the figure.
Note: See TracBrowser for help on using the repository browser.