Tsunami simulation

From AMCGMedia

Sumatra-Andaman Tsunami of 26th December 2004

The Sumatra–Andaman Tsunami of 26th December 2004 is an ideal validation case for ICOM due to recent studies focused on accurately constraining the initial sea surface perturbation and modelling the resultant wave. Pietrzak et al. (2007) reconstructed a series of possible initial free surface displacements derived from coseismic analysis of Global Positioning System (GPS) and uplift data. Validating the resultant waves against altimetric satellite track data and tsunami arrival times as recorded by tide gauges showed this to be an accurate way of determining the initial conditions. The same free surface initial conditions can be incorporated directly into ICOM for a tsunami validation simulation.

To simulate the Sumatra–Andaman Tsunami with ICOM, an unstructured mesh was created for the Indian Ocean domain with boundaries of 60°E to 110°E and−20°S and 10°N. Using GMT (Generic Mapping Tools, 2008), the 1 arcminute resolution GEBCO (General Bathymetric Chart of the Oceans, 2003) bathymetry datasetwas sub-sampled to 6 arcminutes, and this was subsequently processed with Terreno (Gorman et al., 2006, 2007, 2008) to generate a mesh of 103,864 nodes and 304,275 tetrahedral elements. The same source conditions as used in 'Model 1' of Pietrzak et al. (2007) were applied to generate the tsunami (Fig. 1). Detectors were placed at the same locations as used by Pietrzak et al. (2007) so as to enable a comparison of arrival times (Fig. 1). The first arrival time was taken as the first maxima in the free surface elevation time series, and is compared with tide gauge data.

Fig. 1: (A) Bathymetry for the Indian Ocean domain and (B) initial conditions for the Sumatra–Andaman Tsunami of the 26th December 2004 as used by Pietrzak et al. (2007) (their ‘Model 1’). Locations of tide gauges are also annotated.

ICOM accurately simulates the anisotropic nature of the tsunami as it propagates away from the source region (Fig. 2). There is a very good correlation between ICOM's simulated arrival times and the recorded tide gauge data with the exception of Port Blair (Table 1). The result for Port Blair, with a difference of 47 min, is highly anomalous and is probably a result of the complex, shallow bathymetry near this location that could not be captured in the mesh without unreasonable computational expense. Furthermore, the absence of a wetting and drying implementation in ICOM and the consequent requirement for a fixed depth of −20 m at the shoreline might act to compound any errors affecting the calculated arrival time. When Port Blair is omitted, both the mean of the absolute error (4.6 min) and RMS difference (6.7 min) are within a 5–10minute error margin that can be attributed to grid resolution or bathymetric uncertainty in coastal regions (Pietrzak et al., 2007). Such errors are inherent in the GEBCO dataset, which is poorly resolved and heavily dependent on interpolation in many areas (Ham, 2006). Many of the predicted arrival times are also within the 3–6 minute error margin of the tide gauges themselves (Pietrzak et al., 2007).

Fig. 2: Propagation of the Sumatra–Andaman Tsunami of the 26th December 2004. Free surface elevation contours are shown for plus and minus 5, 10, 20 and 40 cm. ICOM accurately simulates the anisotropic nature of the tsunami as it evolves in time.

Table 1: Comparison of recorded arrival times of the Sumatra–Andaman Tsunami at 9 tide gauge locations with simulated results from ICOM. Port Blair generates significant error. With Port Blair omitted, the mean absolute difference and Root Mean Square (RMS) of the differences lie within an error margin (5–10min) that could be attributed to the accuracy of the bathymetry.

Grand Banks Tsunami of 1929

ICOM has been used to simulate the propagation of 1929 landslide generated 'Grand Banks Tsunami' across the North Atlantic. The model uses a fixed unstructured tetrahedral mesh composed of ~210,000 nodes set to a spherical projection and incorporates a source generated with a Delta function as approximated by a Gaussian of a 12km arbitrarily determined width.

The model accurately simulates the anisotropic nature of the wave as it propagates away from the source region (see images below).

(a) 1 hour

(b) 2 hours

(c) 4 hours

(d) 8 hours

Palaeo-tsunami applications

ICOM has been used to simulate palaeo-tsunamis in palaeobathymetric domains to address the debate concerning the scarcity of tsunami deposits in the geological record. Despite tsunamis being a common phenomena in modern day seas and oceans there is relatively little sedimentological evidence for their occurrence in the expansive epicontinental sea deposits that dominate much of the marine stratigraphic record. Several reasons have been suggested for this including mis-identification of facies and poor preservation potential. It is also possible however that their infrequency is due to poor propagation potential for tsunamis in shallow water. Much work has been done addressing the way in which tidal waves are dissipated by frictional processes in epicontinental seas (e.g. see Wells, 2008). It therefore seems possible that tsunamis, also longwaves that interact with the bathymetry, are subject to similar hydrodynamic effects.

Fig. 4: Palaeobathymetry for the Hettangian Stage of the Laurasian Seaway. Main source references are Ziegler (1990); Pienkowski (1991); Dercourt et al. (2000) and Golonka (2004, 2007)

ICOM has been used to investigate this problem via assessing the propagation potential of idealised tsunamis in a typical epicontinental sea; the Lower Jurassic Laurasian Seaway (present day Europe). A palaeobathymetric domain was constructed for the Hettangian Stage from a compilation of published palaeogeographic maps. Palaeodepths were estimated at various locations in the domain from published descriptions and interpretations of sedimentary deposits. When combined with the mapped facies boundaries in the palaeogeographic maps, the palaeodepths could be interpolated over the whole domain (Fig. 4). A one layer deep unstructured mesh of tetrahedral elements was constructed from the palaeobathymetric data using Terreno (Gorman et al., 2006, 2007, 2008).

Idealised point sources were generated using a Gaussian with an amplitude of 25 m and a FWHM of 70 km. These are of sufficient size to generate large tsunamis, the kind of which would have probably been infrequent even on geological timescales. Two simulations were run on a single element deep Terreno generated mesh; the first ('Simulation 1') was positioned in the deep water (approx. 5 km) of the Tethys Ocean and the second ('Simulation 2') in the shallow water of the shelf itself. Bed shear stress was calculated to establish the potential for the respective tsunamis to influence sedimentation/erosion over the shelf.

Results show that the waves, limited by the bathymetry, propagate slowly in both cases (Fig. 5). The maximum bed shear stress plots for Simulations 1 and 2 show that the potential of tsunamis to move sediment within the Laurasian Seaway is spatially limited (Fig. 6). The region where the threshold bed shear stress for the entrainment of sand in water (0.06 Pa) is exceeded is never greater in extent than ∼2000 km from the shelf break or source region in ‘Simulation 1’ or ‘Simulation 2’ (Fig. 6C and D). The ability for the tsunami to move gravel grade material of 2mm in diameter is more restricted and never extends more than ∼500 km from the source (Fig. 6E and F).

Fig. 5: Results for both ‘Simulation 1’ and ‘Simulation 2’ after, 1, 2, 4 and 8 hour simulated times respectively. Both are forced with the ‘base case’ idealised Gaussian source (a=25m and c=70 km). The shallow depths of the Laurasian Sea decrease the wave speed causing a very slow propagation across the shelf. Notably, ‘Simulation 2’ has a higher wave amplitude than ‘Simulation 1’. ‘Simulation 2’ is however more localised due to the shallow water preventing widespread dispersion.

Maximum bed shear stress values are significantly higher for ‘Simulation 2’ on the shallow water shelf region than for ‘Simulation 1’, with values of 30 Pa recorded for the former (Fig. 6B). This can be attributed to the greater magnitude of the source perturbation relative to the depth of water. Values for ‘Simulation 1’ on the other hand barely exceed more than 2.4 Pa and only values of up to ∼0.5 Pa can be mapped over any large area (Fig. 6A).

Fig. 6: Maximum bed shear stress for (A) ‘Simulation 1’ and (B) ‘Simulation 2’. Regions where the critical threshold bed shear stress for the entrainment of very fine sand (0.06 Pa) and gravel (1.0 Pa) is exceeded are highlighted in red for ‘Simulation 1’ (C and E) and ‘Simulation 2’ (D and F) respectively. For both source regions the tsunamis are rapidly damped down and only capable of moving sand over a relatively small portion of the epicontinental sea. Notably the maximum bed shear stresses generated for ‘Simulation 2’ are significantly greater than those of ‘Simulation 1’. This can be attributed to the size of the initial free surface perturbation being much greater relative to the depth of water.

Ultimately, Fig. 6C-F shows how even these large tsunamis are incapable of reworking sediment over an expansive area. This is probably due to a combination of frictional drag, energy loss on shoaling (for the deep water source entering the sea) and interference with islands. This suggests that tsunami deposits are only likely to be formed in close proximity to the source for large events and thus may help explain their apparent scarcity on the geological record

References

Dercourt, J., Gaetani, M., Vrielynck, B., Barrier, E., Biju-Duval, B., Brunet, M. F., Cadet, J. P., Crasquin, S., Sandulescu, M. (Eds.), 2000. Atlas Peri-Tethys, palaeogeographical maps. CCGM/CGMW, Paris: 24 maps and explanatory notes: I-XX.

Golonka, J., 2004. Plate tectonic evolution of the southern margin of Eurasia in the Mesozoic and Cenozoic. Tectonophysics 381, 235–273.

Golonka, J., 2007. Phanerozoic paleoenvironment and paleolithofacies maps. Mesozoic. Geologica 33, 211–264.

Gorman, G.J., Piggott, M.D., Pain, C.C., 2007. Shoreline approximation for unstructured mesh generation. Computers and Geosciences 33, 666–677.

Gorman, G.J., Piggott, M.D., Pain, C.C., de Oliveira, C.R.E., Umpleby, A.P., Goddard, A.J.H., 2006. Optimisation based bathymetry approximation through constrained unstructured mesh adaptivity. Ocean Modelling 12, 436–452.

Gorman, G.J., Piggott, M.D., Wells, M.R., Pain, C.C., Allison, P.A., 2008. A systematic approach to unstructured mesh generation for ocean modelling using GMT and Terreno. Computers and Geosciences 34, 1721–1731.

Ham, D., 2006. On techniques for modelling coastal and ocean flow with unstructured meshes. Ph.D. thesis, Technische Universiteit Delft.

Pienkowski, G., 1991. Eustatically-controlled sedimentation in the Hettangian–Sinemurian (Early Jurassic) of Poland and Sweden. Sedimentology 38, 503–518.

Pietrzak, J., Socquet, A., Ham, D., Simons, W., Labeur, C.V.R.J., Schrama, E., Stelling, G., Vatvani, D., 2007. Defining the source region of the Indian Ocean Tsunami from GPS, altimeters, tide gauges and tsunami models. Earth and Planetary Science Letters 261, 49–64.

Wells, M. R., 2008. Tidal modeling of modern and ancient seas and oceans. Ph.D. thesis, Imperial College, London.

AMCG