- Email: [email protected]
A synoptic picture of the impact of the 26th December 2004 Indian Ocean tsunami on the coast of Sri Lanka
Environmental Modelling & Software 25 (2010) 1874e1880
Contents lists available at ScienceDirect
Environmental Modelling & Software journal homepage: www.elsevier.com/locate/envsoft
A synoptic picture of the impact of the 26th December 2004 Indian Ocean tsunami on the coast of Sri Lanka M. Ioualalen a, *, W. Rentería b, K. Ilayaraja c, M. Chlieh a, P. Arreaga-Vargas b a
Institut de Recherche pour le Développement, IRD, GéoAzur, Observatoire Océanologique, 2 Quai de la Darse, BP 48, F-06235 Villefranche-sur-mer Cedex, France Instituto Oceanografico de la Armada de Ecuador, INOCAR, Av. 25 de Julio-Km 3½, via Pto. Maritimo Base Naval Sur Guayaquil, Ecuador c Department of Applied Geology, University of Madras, Guindy Campus, Chennai 600 025, India b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 17 July 2009 Received in revised form 31 January 2010 Accepted 16 April 2010 Available online 21 May 2010
A numerical simulation of the 26th December 2004 Indian Ocean tsunami for the entire coast of Sri Lanka is presented. The simulation approach is based on a fully nonlinear Boussinesq tsunami propagation model and a robust coseismic source. The simulation is first confronted to available measured wave height. The agreement between observations and the predicted wave heights allowed a reasonable validation of the simulation. As a result a synoptic picture of the tsunami impact is provided over the entire coast of Sri Lanka. It is found that amplification due to shoaling applies mainly in the Eastern and Southern coast because, here, the wave is propagating across the sea floor slope, while propagating along the slope for the Western coast. Spots of high waves are due to wave focusing in some coastal areas while local submarine canyons in other areas inhibit the wave amplification. Ó 2010 Elsevier Ltd. All rights reserved.
Keywords: 2004 tsunami Sumatra Runup Numerical simulation Sri Lanka
1. Introduction On Sunday, 26 December 2004 at 00:58:53 UTC, the Mw ¼ 9.1e9.3 SumatraeAndaman earthquake triggered a large tsunami that severely damaged coastal communities in countries along the Indian Ocean, including Indonesia, Thailand, Sri Lanka, India, Maldives, and Somalia located at more than 6000 km, which corresponds to a wave travel time up to 8 h. The epicenter and first aftershocks indicate that approximately 1300 km of subduction interface ruptured along the northern Sunda Trench. The tectonics and seismicity of the bay of Bengal can be found in Socquet et al. (2006), Stein and Okal (2005) and the recurrence of earthquaketriggered tsunamis is discussed in Altinok and Ersoy (2000). As far as the tsunami is concerned, Sri Lanka was the second area after Sumatra, Indonesia, that encountered most damages and casualties. A series of waves impacted more than two third of the coast. The tsunami left e31,000 human deaths, e4110 missing, 440,000 displaced. It provoked severe damages in harbours, infrastructures including housing, livelihoods. Communication and transportation systems were seriously affected. A great level of pollution was
* Corresponding author. E-mail address: [email protected] (M. Ioualalen). 1364-8152/$ e see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.envsoft.2010.04.010
observed in the flooded areas (e.g., in drinking water reservoirs). Naturally, the effects were increased or reduced by the shape of the coastline, geomorphologic features and the bathymetric variation, i.e., the physical processes involved. But also and importantly, the effective effects of the hitting were different depending on the level of exposure and vulnerability of the coastal communities to the impact of waves. Such parameters may also vary at a daily time scale, e.g., transport- traffic, rush hours. The measurements associated to the physical phenomena, the runup and inundation mainly, cannot be taken as a measure of the impact of tsunami because there are other facts that influences in the effects. However, there are cases where the impact magnitude is related strongest to these values. For example, two of the highest values for the height of wave and inundation were observed in Kalmunai and Hambantota (Liu et al., 2005; Fig. 1; Table 1); in the latter location, the number of victims reached 1000 and, as resulted of flooding, the lagoon near to the town was contaminated with debris swept for the tsunami inland. In Kalmunai, the wave flooded and caused damages as far 1500 m inland, within this zone all structures was severely damaged (Liu et al., 2005). In other places like Payagala, were the wave heights were weaker (Liu et al., 2005; Fig. 1; Table 1), the amount of fatalities was significantly high, 800 deaths. Here, a train with numerous passengers happened to circulate along the coast at the time the tsunami reached it. Also, another important feature is the level of protection of a coastal area
M. Ioualalen et al. / Environmental Modelling & Software 25 (2010) 1874e1880
1875
Fig. 1. Plotted measured and simulated wave height at locations mentioned in Table 1. GI and CH stand for the simulations based on the coseismic source of Grilli et al. (2007) and Chlieh et al. (2007) respectively. Isobathes 1000 m (brown curve), 500 m (green) and 300 m (blue) are reported. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
by natural barriers, like coral reef, mangroves and sand dunes, whatever is the wave height. Contrary to the southern and eastern coastline, in the eastern coast, where none of them exist, the damages were more important. In that area, poor housing situated over the coastline were converted by the tsunami in debris that yielded the re-enforcement of the damages. The numerical simulation that we present here for the Sri Lanka case study has two objectives: The first purpose is to provide
a synoptic picture of the event and, in particular, the runup or wave height maxima distribution along a coastline. It is suitable that the numerical simulation is not constrained by runup observations in order to assess a predicted runup distribution (this is the case in the present study). The other interest of numerical modeling is to identify physical processes that are responsible for local wave amplification or attenuation. Then the picture reveals vulnerable areas as well as sheltered ones. These issues are crucial in regions
1876
M. Ioualalen et al. / Environmental Modelling & Software 25 (2010) 1874e1880
Table 1 Estimates of wave height maxima relative to the mean sea level along the Southern coast Sri Lanka sampled by Liu et al. (2005) for the 26th December 2004 Indian Ocean tsunami along with simulated ones at the nearest locations of the computational domain (not necessarily the highest simulated wave) for 3 grid spacing: 0.2780 (in bold), 0.5560 and 1.1120 . SGI stands for Grilli et al. (2007) and Ioualalen et al. (2007a) tsunamiederived source and SCh is for the one proposed by Chlieh et al. (2007). Location
Long. ( E)
Lat. ( N)
Wave height max. (m) Obs.
SGI 0.2780
SGI 0.5560
SGI 1.1120
SCh 0.2780
SCh 0.5560
SCh 1.1120
(1) Trincomalee Hotel (2) Trincomalee town (3) Trincomalee military ChcPt. (4) Trincomalee HI Hostel (5) China Bay Stop 1 (6) China Bay Stop 2 (7) China Bay Stop 3 (8) N of Bathticaloa (9) Kattativu (10) Karativu (11) Ninto (12) Nalaveli Hostel Pottuvil (13) Ibral Nagar Nalaveli 1 (14) Ibral Nagar Nalaveli 3 (15) Ibral Nagar Nalaveli 4-5 (16) 2 Km SE of Kuchchaveli (17) Mankeri (18) Kalmurai Kuddi 1 (19) Kalmurai Kuddi 2 (20) Moratuwa (21) Koralawella (22) Wadduwa (23) Hambantota (24) Nonagama (25) Tangalla 1 (26) Kamburugama (27) Weligama (28) Galle (29) Dodanduwa (30) Hikkaduwa (31) Thiranagama (32) Galbokka (33) Seenigama (34) Dehiwala (35) Panadura (36) Pinnatara (37) Kalutara (38) Payagala (39) Yala (40) Boosa
81.21845 81.24202 81.23478 81.21380 81.19094 81.19206 81.26444 81.69267 81.74025 81.85432 81.86083 81.18847 81.21794 81.21918 81.21674 81.12103 81.48953 81.84164 81.83044 79.88353 79.88879 79.92110 81.12752 80.98835 80.79562 80.49195 80.44682 80.24915 80.14692 80.10413 80.12357 80.03091 80.08908 79.85640 79.90334 79.91306 79.94800 79.97835 81.25503 80.09014
8.618075 8.561729 8.563063 8.579356 8.503758 8.492397 8.463147 7.744343 7.686330 7.364405 7.343535 8.706572 8.660616 8.660272 8.661204 8.790385 8.013957 7.405348 7.422988 6.762450 6.749667 6.673167 6.128450 6.093750 6.029367 5.940050 5.968984 6.009734 6.083717 6.127517 6.110650 6.323317 6.166100 6.877667 6.715483 6.689667 6.608450 6.521217 6.166390 6.047760
4.4 Mean: 4.30 2.5 3.2 2.8 3.35 3.25 2.7 3.7 4.85 4.5 4.1-6.1 4.65 4.5 Mean: 4.13 3.35 5.5 6.2 7 3.56 4.55 3.61 6.1 8.71 3.24 2.4 2.73 5.24 3.6 4.2 4.55 4.26 5.05 3.48 4.24 4.15 3.82 5.04 4.65 Mean: 2.00
4.35 3.20 3.29 3.00 3.80 1.95 3.84 3.70 4.77 10.43 10.09 3.69 3.02 3.02 3.02 3.46 7.11 9.80 9.75 2.82 2.71 4.15 5.24 6.51 5.23 3.78 5.52 5.62 4.64 5.70 5.18 6.77 6.17 2.59 2.20 3.50 3.88 6.28 4.93 2.71
4.22 3.33 3.33 3.69 2.01 2.34 3.64 3.86 4.78 11.41 9.17 2.7 2.59 2.59 2.59 3.39 4.33 9.37 10.14 2.22 2.22 2.86 6.93 7.40 6.01 4.15 4.99 3.78 4.29 2.59 1.56 6.04 6.07 2.63 3.26 2.85 3.89 5.17 3.68 2.73
3.46 3.41 3.41 3.17 2.23 2.23 2.63 3.77 4.64 9.37 8.99 4.46 2.61 2.61 2.61 3.10 8.62 9.27 8.29 3.90 4.00 4.21 4.43 5.86 4.89 3.71 5.96 5.00 6.38 7.28 6.72 9.68 9.75 1.70 3.07 3.21 4.46 6.46 4.76 3.55
4.21 3.11 3.4 2.82 5.66 5.81 5.17 3.30 4.80 6.82 6.44 3.67 3.20 3.20 3.20 4.43 7.71 5.19 4.80 1.88 1.63 3.09 5.11 6.76 3.83 2.90 6.43 4.45 2.84 4.05 3.71 4.22 3.88 1.38 2.15 2.54 3.14 2.77 3.96 1.90
3.85 4.10 4.10 4.85 7.59 8.71 7.56 4.22 4.72 5.78 6.41 3.33 2.79 2.79 2.79 4.56 7.61 5.55 5.45 2.50 2.30 3.36 4.35 5.86 3.04 2.85 6.63 4.46 3.21 4.61 3.81 4.26 4.08 1.49 2.25 3.18 2.76 2.89 3.09 3.75
4.62 3.78 3.78 3.56 3.95 3.06 4.45 3.38 4.80 5.59 5.25 4.73 2.71 2.71 2.71 3.91 9.33 4.44 5.74 3.01 2.83 3.43 4.44 5.31 3.78 3.00 4.79 4.05 3.88 3.67 3.39 6.31 4.76 1.30 2.29 3.04 3.04 4.12 4.74 5.11
where tsunami events are frequent. Such methodology has been applied successfully by Ioualalen et al. (2007a) for the same tsunami event but on its impact on the Andaman coast of Thailand. The runup/wave height data set that is used here for validation of the numerical simulation is the one gathered by Liu et al. (2005) immediately after the event during the period 9e15 January 2005 (Table 1). Most of their sampling concerned wave height maxima (WHM) that have been derived from debris in trees and watermarks on buildings. Both the observed and our simulated WHMs are referenced to the mean sea level: the comparison between the two data sets is then direct because in both cases we have the freetide tsunami signal. 2. Numerical simulation of the tsunami for Sri Lanka 2.1. The numerical procedure We used here FUNWAVE tsunami propagation and runup/ inundation model. The model is fully nonlinear and dispersive, retaining information to leading order in frequency dispersion O [(kh)2] and to all orders in nonlinearity a/h (where k denotes an inverse wavelength scale, a denotes a wave amplitude, and h denotes a water depth) (Wei and Kirby, 1995; Wei et al., 1995). The
model treats the entire computational domain as an active fluid domain by employing an improved version of the slot or permeable-seabed technique, i.e., the moving shoreline algorithm proposed by Chen et al. (2000) and Kennedy et al. (2000) for simulation of runup. The model includes bottom friction, energy dissipation to account for the wave breaking and a subgrid turbulence scheme. We have constructed a 0.278 arc minute grid-spacing for the computational domain (e520 m) ranging from 2 N to 13 N and from 77.2 E to 96 E (Fig. 2). The domain is composed of 3958 2377 nodes and we chose an optimal time step of 0.5 s to avoid numerical instabilities and to limit truncation errors. A propagation time of 6 h has been chosen (43,200 time steps) to take into account the highest crest observed in Sri Lanka (Pattiaratchi and Wijeratne, 2009; Liu et al., 2005). Although, the grid spacing seems to be relatively large, Ioualalen et al. (2007a) used approximately the same one for the same event applied for the Thailand case study and they discussed the relevance of the parameter for this specific event. They obtained a satisfactory comparison with observations, e.g., cross-correlation between observed and simulated runup were 0.94 and the RMS difference was only about 17%. Three sets of data have been used: Offshore the Islands 2arc minute ETOPO-2 bathymetry has been
M. Ioualalen et al. / Environmental Modelling & Software 25 (2010) 1874e1880
1877
Fig. 2. Computational domain including the vertical sea floor deformation at 1 m contour intervals predicted from the coseismic SumatraeAndaman model G-M9.15 of Chlieh et al. (2007). Dashed lines indicate subsided area and the continuous lines the uplifted regions. The background bathymetry is plotted in grey at 500 m contour intervals.
used until depth 200 m. Then for shallow water digitalized marine charts data have been used. Inland, the SRTM DEM, with a horizontal resolution of 90 m, have been used.
2.2. The initial wave solution As input in our numerical modeling, we used two already validated seismic sources to compute the vertical seafloor deformation associated with the SumatraeAndaman earthquake: The first one, based on GPS network, has been proposed by Chlieh et al. (2007) and the second, developed by Grilli et al. (2007) and Ioualalen et al. (2007a), has been calibrated with available tide gauges records in the Indian Ocean and with Jason sea level anomalies data. The source G-M9.15 developped by Chlieh et al. (2007) integrates observations from several geodetic observations (Figs. 2 and 3): a near-field survey-mode GPS sites on northern Sumatra, several reoccupied GPS sites in the Andaman and the Nicobar islands, a set of more distant continuous GPS stations in Thailand and Malaysia, measured patterns of coral reefs vertical motion and the hinge-line deduced from satellite imagery. The joint inversion of all these geodetic data provided good constraints on the coseismic source. This source has been shown to be in good agreement with many others geophysical and hydrographic data, i.e, seismic observations, normal modes up to 300 s T-wave recorded from hydroacoustic arrays and sea level variations monitored by altimeters. Grilli et al. (2007) and Ioualalen et al. (2007a) calibrated a more energetic Mw ¼ 9.25 five-segment Okada (1985) dislocation source (Fig. 3; Table 2) using available hydrodynamic data, in particular the sea level anomalies recorded by the JASON-1 altimeter, which happened to transit over the area spawned by the spreading tsunami about 2 h after the earthquake initiation, during cycle 109 of its pass 129 (Smith et al., 2005), and several digital tide gauge records installed along the Gulf of Bengal coasts (Merrifield et al., 2005). This source has been tested and validated by Ioualalen et al. (2007a) for the Thailand case study. They were able to reproduce most of the runup features over the Andaman coast of Thailand. The runups obtained in a predictive mode (the simulations were not constrained by runup observations) fitted observations at a relatively high degree of accuracy. The two coseismic sources differ significantly in terms of seismic moment amplitude and also in terms of vertical deformation
distribution along the subduction zone. The gaps are mainly due to the respective modes of calibration: The source of Chlieh et al. (2007) has been constrained more objectively because it relies on effective GPS sea floor observation. The one of Grilli et al. (2007) has been built through hydrographic observations and thus only integrates the wave mode of propagation only. Besides it is built within Okada (1985)’s idealized assumptions. Ideally the two sources should converge to each other but it is not the case because (1) the functions of transfer between the floor deformation and the initial wave are not represented in the simulation (the transient movement within the water column) and (2) the hydrographic source highly depends on the quality of the ensemble of numerical simulations that have been performed to derive it.
2.3. The numerical results and discussion We perform a numerical simulation to provide a synoptic picture of the event along the entire coast of Sri Lanka and to analyze the features that are responsible for their coastal vulnerability or sheltering. To reach this objective we validate first our wave height results with available observations (Table 3, Fig. 1). Elementary statistics on simulations skills for the 0.278 arc minute grid-spacing simulation indicate that the numerical simulations based on Grilli et al. (2007) and Ioualalen et al. (2007a) and Chlieh et al. (2007) co-sesimic sources (named SGI and SCh herafter respectively) seem to give complementary results: SCh shows a better-fit for wave heights (mean, std and the rmse, i.e., the model skill) while GI exhibits a better cross-correlation with observations (Table 3). The amplitudes obtained with simulation SCh (6% of mean under-prediction and 18% of std over-prediction) are quite satisfactory which indicates that the seismic momentum and spatial distribution of the rate of subsidence/uplift of the Chlieh et al. (2007)’s source fits the tsunami wave heights requirements. The relatively poor correlation with observations (0.38) might be explained by the fact that the initial wave is triggered at once with no lag along the ruptured fault. As a reference, Ilayaraja et al. (2009) using the same source of Chlieh et al. (2007) for the same tsunami event and the case study of Andaman and Nicobar Islands found a cross-correlation of 0.88. For that particular case the Islands are located within the ruptured area and thus the risk of tsunami phase-lag with observation is considerably reduced (not the case for the remote Sri Lanka). We may then expect that a proper triggering process of the initial wave may improve substantially our
1878
M. Ioualalen et al. / Environmental Modelling & Software 25 (2010) 1874e1880
Fig. 3. Vertical sea floor deformation at 1 m contour intervals predicted from (left panel) SCh, the coseismic SumatraeAndaman model G-M9.15 of Chlieh et al. (2007) and (right panel) SGI, the solution of Grilli et al. (2007) and Ioualalen et al. (2007a) which Okada (1985)’s dislocation parameters are provided in Table 2. Dashed lines indicate subsided area and the continuous lines the uplifted regions. The background bathymetry is plotted in grey at 500 m contour intervals.
results. Such parameter has been tested by Wang and Liu (2006) who studied both effects of and impulsive and a transient fault models for the 2004 Indian Ocean tsunami. This is visible when we Table 2 Pure thrust reverse fault solution dipping eastward (Grilli et al., 2007 and Ioualalen et al., 2007a). Input parameters for Okada (1985) dislocation method (first 5 lines) and outputs (last 4 lines) for 5 tsunami source segments (Fig. 3): time delay of segment rupture from earthquake time s (a 60 s rising time is added); longitude and latitude of segment centroid (xo, yo); the centroid depth is d ¼ 25 km for all segments, the fault strike angle 4 (clockwise from North); the fault rake angle is l ¼ 90 for all segments (counterclock-wise from strike); the fault dip angle is d ¼ 12 everywhere (positive from the horizontal plane); the maximum fault slip D; the segment length along and width across (L, W); and the medium shear modulous taken m ¼ 4 1010 Pa for all segments; the seismic moment Mo; the characteristic initial tsunami wavelength lo and period so; and the characteristic tsunami trough and peak amplitudes ho. Note that, in the simulation, slip is highest at the segments’ centroids and drops by 50% at a radius of L from it. The total seismic moment of all 5 segments is Mo ¼ 7.55 1022 or Mw ¼ 9.25. Parameters S1
S2
S3
S4
S5
s (s)
272 93.90 E, 5.22 N 348 23 150, 130 1.58 1022 130 17.46 3.84, þ8.59
588 93.21 E, 7.41 N 338 12 390, 120 2.05 1022 120 23.30 2.33, þ4.72
913 92.60 E, 9.70 N 356 12 150, 95 0.61 1022 95 18.72 2.08, þ4.49
1273 92.87 E, 11.70 N 10 12 35 , 95 1.46 1022 95 18.72 2.31, þ4.60
xo, yo
4 D (m) L, W (km) Mo (J) lo (km) so (min.) ho (m)
60 94.57 E, 3.83 N 323 18 220, 130 1.85 1022 130 24.77 3.27, þ7.02
consider the simulation SGI: With relatively poor amplitudes results we obtain a better correlation better than for simulation SCh probably because we took into account some mode of initial wave triggering including time lags between the successive components of the segmented ruptured area (Table 2). The wave height overprediction of the simulation SGI is probably due to the fact that the authors calibrated their source with available tide gauges to the east with a 1 arc-minute-only grid-spacing for their computational grid. It is fair to say that their source would have been much improved if they used a finer gridding: Ioualalen et al. (2007a) showed for the Thailand case study a significant improvement of the simulated runup when from a 1arc minute grid spacing (mean over-prediction of 34%) to a 0.25arc minute (like here) (mean overprediction of 17% only) for the same GRI triggering source. However, since the FUNWAVE version we use does not allow multiple-grids we encounter difficulties to compute a finer gridspacing because of computing limitations. Since our validation conclusions are reasonable but less satisfactory than those of Ioualalen et al. (2007a) for the Thailand case study we need to perform sensitivity tests on the grid-spacing in order to estimate the convergence of our simulation and at least say if the simulations are in a coherent trend. For that purpose we have performed two simulations for each tsunamis source SGI and SCh with gridspacings 0.558 (e1040 m) and 1.112 (e2080 m) arc minute. Thanks to these grid tests, the basic statistics indicate that the SGI simulation is in a better convergence trend (Table 3): as far as the amplitude is concerned (std and rmse), the SGI simulations indicate a fair improvement when improving the grid spacing while no
M. Ioualalen et al. / Environmental Modelling & Software 25 (2010) 1874e1880
1879
Table 3 Comparison between the observed (OBS) and simulated runup or wave height maxima for coseismic sources SGI and SCh used for grid spacings 0.2780 (in bold), 0.5560 and 1.1120 . The basic statistics are the mean, the standard deviation, the root mean squared error (RMSE) (and their deviation in % from the observations), the cross-correlation qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P 2 coefficient and the norm L2. The usual RMSE representing the model skill is defined as i ðroi rsi Þ =n, where roi is the observed runup for each (i)-location of Table 1, rsi is for qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P 2 P 2 2 the associated simulated one and n is 40. The Norm L is defined as i ðroi rsi Þ = i ðroi Þ
Mean (m) STD (m) RMSE (m) Cross-correlation Norm L2
OBS
SGI 0.2780
SGI 0.5560
SGI 1.1120
Sch 0.2780
Sch 0.5580
Sch 1.1120
4.23 1.28
4.79 (D13%) 2.17 (D69%) 1.89 (45%) 0.55 0.43
4.37 (þ3%) 2.32 (þ81%) 1.98 (47%) 0.52 0.45
4.95 (þ17%) 2.32 (þ81%) 2.17 (51%) 0.47 0.49
3.99 (-6%) 1.51 (D18%) 1.58 (37%) 0.38 0.36
4.29 (þ1%) 1.66 (þ30%) 1.99 (47%) 0.10 0.45
4.07 (4%) 1.33 (þ4%) 1.48 (35%) 0.37 0.33
convergence trend appears for SCh ones. We obtain the same conclusion for the special distribution of the simulated runup: we improve the cross-correlation with finer grid spacing while no trend appears for SCh. We also tend to minimize the Norm L2 with SGI while there is no such feature for SCh. As a consequence, although the 0.278 arc minute grid, for SGI do not provide excellent results, we may assert that we place in a reasonable convergence trend and we may select this simulation for further process studies. This sensitivity test allowed us to select in a robust way the more coherent simulation between SGI and SCh while it is difficult to do so with only the absolute values of the 0.278 arc minute grid simulations. It is also fair to say that despite the relatively scarce grid spacing, the inherent uncertainties on the data used for the gridding (for example SRTM topography data may introduce errors on the slope of the beach and thus on the wave height), the possible gap between the positioning for observations and simulation locations (we compared heights at same locations or nearest point), we obtained a relatively coherent simulation using SGI triggering source. Based on that reasonable validation of the SGI simulation, we propose now a synoptic picture of the event (wave height distribution along the entire coast of Sri Lanka) and identify the physical processes that are responsible for the disparities. We propose here only the vulnerability picture regards to the physical event. Obviously, we do not discuss other vulnerability parameters suggested in the introduction (natural barriers, socio-economic features etc). Fig. 4 shows the predicted runup along entire coast of Sri Lanka. We use here the term “prediction” because the simulation has not been constrained by any observed runup or wave height values so that we place in a prognostic mode. In our case study, the SouthEastern coast of India exhibits a relatively pronounced coastal variability on tsunami impacting with clear risk areas and more sheltered ones for any tsunami event. This even can be considered as the worst case scenario for Sri Lanka because the threatening recurrent earthquakes happen in the considered subduction area and also because the Mw ¼ 9.1e9.3 event is extreme (based on an extreme size of the ruptured area). For that simple reason, the study of this event can be considered as a relevant for tsunami threat mapping. In the following we try to identify the features that are responsible for the runup distribution along the coast of Sri Lanka. It is obviously not the objective of this work to provide a GIS mapping of the tsunami risk. The wave amplification factors are relatively easy to identify. First the shoaling applies everywhere in the Eastern and Southern coasts where the wave reaches the coast at a pronounced incidence (Fig. 4) and propagates across a slope. This is not the case for the Western coast where the wave propagates along the slope, especially North of 7 N. Here, the flux of wave energy is conserved along a quasi-constant slope and thus does not yield a slope (or Green’s) effect.
Fig. 4 exhibits areas of wave refraction/focusing at the Eastern coast where spots of very high waves appear (noted FC in the Figure with red arrows). This is the case at around 9.2 N, 8.0 N and 7.3 N where we observe convex-type isobaths. At around 8.0 N, although Google Earth indicates the presence of deserted beaches, this non-occupied area is likely to be vulnerable. Because of its poor density in population, the 2004 event did not generate significant damages in that area. However the period of recurrence of such event (centenaries time scale), future development plan of the area should take into account the potential threat.
Fig. 4. Simulated maxima of wave height (in meters) recorded over the entire time integration for simulation SCh. Isobaths 30, 100 and 500 m are represented in while curves. FC stands for “focusing” and UC for underwater canyons. Blue arrows indicate wave divergence and red ones correspond to wave convergence (focusing). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
1880
M. Ioualalen et al. / Environmental Modelling & Software 25 (2010) 1874e1880
Fig. 4 also indicates areas where underwater canyons are present (noted UC in the figure with blue arrows). At their termination (the coast line), the wave is diverging (concave-type isobaths) and escapes from its sides (a defocusing). The effects of these canyons are to inhibit the wave amplification and eventually balance shoaling effects. This has been shown already for the Bangladesh coast and the same event where a multitude of underwater canyons have defocused the wave which encountered highest waves only offshore at the continental shelf (Ioualalen et al., 2007b). This explains why the coast of Bangladesh has been so little affected by the event. Acknowledgements The authors would like to thank the French Agence Nationale Pour la Recherche, ANR, for funding this work through the grant TSUMOD ANR-05-CATT-016-02. M. Ioualalen addresses his thanks to the Department Milieu et Environement of the Institut de Recherche pour le Développement, IRD, for having granted him a long stay at INOCAR, Ecuador, especially J. Boulegue and P. Soler and P. Charvis and J.-Y. Collot (Géosciences Azur) for their constant support. References Altinok, Y., Ersoy, S., 2000. Tsunamis observed on and near the Turkish coast. Nat. Hazards 21, 185e205. Chen, Q., Kirby, J.T., Dalrymple, R.A., Kennedy, A.B., Chawla, A., 2000. Boussinesq modeling of wave transformation, breaking, and run-up. II: 2D. J. Wtrwy, Port, Coast, Oc. Engrg., ASCE 126 (1), 48e56. Chlieh, M., Avouac, J.-P., Hjorleifsdottir, V., Song, T.-R.A., Ji, C., Sieh, K., Sladen, A., Hebert, H., Prawirodirdjo, L., Bock, Y., Galetska, J., 2007. Coseismic slip and afterslip of the great Mw 9.15 SumatraeAndaman earthquake of 2004. Bull. Seism. Soc. Amer. 97 (1A), S152eS173. Grilli, S.T., Ioualalen, M., Asavanant, J., Shi, F., Kirby, J.T., Watts, P., 2007. Source constraints and model simulation of the December 26, 2004 Indian Ocean tsunami. J. Wtrwy, Port, Coast, Oc. Engrg., ASCE 133 (6), 414e428.
Ilayaraja, K., Ioualalen, M., Chlieh, M., Krishnamurthy, M.V., Ramana Murthy, M.V., 2009. Numerical modeling of the 26th December 2004 India Ocean tsunami at Andaman and Nicobar Islands. J. Coast. Conserv.. doi:10.1007/s11853-008-00338. Ioualalen, M., Asavanant, J., Kaewbanjak, N., Grilli, S.T., Kirby, J.T., Watts, P., 2007a. Modeling of the 26th December 2004 Indian Ocean tsunami: case study of impact in Thailand. J. Geophys. Res. Oceans 112, C07024. doi:10.1029/ 2006JC003850. Ioualalen, M., Pelinovsky, E., Asavanant, J., Lipikorn, R., Deschamps, A., 2007b. On the weak impact of the 26 December Indian Ocean tsunami on the Bangladesh coast. Nat. Hazards Earth Syst. Sci. 7, 141e147. Kennedy, A.B., Chen, Q., Kirby, J.T., Dalrymple, R.A., 2000. Boussinesq modeling of wave transformation, breaking, and run-up. I: 1D. J. Wtrwy, Port, Coast, Oc. Engrg., ASCE 126 (1), 39e47. Liu, P.L.-F., Lynett, P., Fernando, H., Jaffe, B.E., Fritz, H., Higman, B., Morton, R., Goff, J., Synolakis, C., 2005. Observations by the international tsunami survey team in Sri Lanka. Science 308 (5728), 1595. doi:10.1126/science.1110730. Merrifield, M.A., Firing, Y.L., Aarup, T., Agricole, W., Brundrit, G., Chang-Seng, D., Farre, R., Kilonsky, B., Knight, W., Kong, L., Magori, C., Manurung, P., McCreery, C., Mitchell, W., Pillay, S., Schindele, F., Shillington, F., Testut, L., Wijeratne, E.-M., Caldwell, P., Jardin, J., Nakahara, S., Porter, F.-Y., Turetsky, N., 2005. Tide gauges observations of the Indian Ocean tsunami, December 26, 2004. Geophys. Res. Lett. 32, L09603. doi:10.1029/2005GL022610. Okada, Y., 1985. Surface deformation due to shear and tensile faults in a half-space. Bull. Seis. Soc. Am. 75 (4), 1135e1154. Pattiaratchi, C.B., Wijeratne, E.M.S., 2009. Tide gauge observations of 2004e2007 Indian Ocean tsunamis from Sri Lanka and Western Australia. Pure Appl. Geophys. 166, 233e258. Smith, W.H.F., Scharroo, R., Titov, V.V., Arcas, D., Arbic, B.K., 2005. Satellite altimeters measure tsunami: early model estimates confirmed. Oceanography 18 (2), 10. Socquet, A., Vigny, C., Simons, W., Chamot-Rooke, N., Rangin, C., Ambrosius, B., 2006. GPS determination of the relative motion between India and Sunda, and its accommodation in Myanmar. J. Geophys. Res. 111, B05406. doi:10.1029/ 2005JB003877. Stein, S., Okal, E., 2005. Speed and size of the Sumatra earthquake. Nature 434, 581e582. Wang, X., Liu, P.L.-F., 2006. An analysis of 2004 Sumatra earthquake fault plane mechanisms and Indian Ocean tsunami. J. Hydraulic Res. 44 (2), 147e154. Wei, G., Kirby, J.T., 1995. Time-dependent numerical code for extended Boussinesq equations. J. Wtrwy, Port, Coast, Oc. Engrg., ASCE 121 (5), 251e261. Wei, G., Kirby, J.T., Grilli, S.T., Subramanya, R., 1995. A fully nonlinear Boussinesq model for free surface waves. Part 1: highly nonlinear unsteady waves. J. Fluid Mech. 294, 71e92.
Contents lists available at ScienceDirect
Environmental Modelling & Software journal homepage: www.elsevier.com/locate/envsoft
A synoptic picture of the impact of the 26th December 2004 Indian Ocean tsunami on the coast of Sri Lanka M. Ioualalen a, *, W. Rentería b, K. Ilayaraja c, M. Chlieh a, P. Arreaga-Vargas b a
Institut de Recherche pour le Développement, IRD, GéoAzur, Observatoire Océanologique, 2 Quai de la Darse, BP 48, F-06235 Villefranche-sur-mer Cedex, France Instituto Oceanografico de la Armada de Ecuador, INOCAR, Av. 25 de Julio-Km 3½, via Pto. Maritimo Base Naval Sur Guayaquil, Ecuador c Department of Applied Geology, University of Madras, Guindy Campus, Chennai 600 025, India b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 17 July 2009 Received in revised form 31 January 2010 Accepted 16 April 2010 Available online 21 May 2010
A numerical simulation of the 26th December 2004 Indian Ocean tsunami for the entire coast of Sri Lanka is presented. The simulation approach is based on a fully nonlinear Boussinesq tsunami propagation model and a robust coseismic source. The simulation is first confronted to available measured wave height. The agreement between observations and the predicted wave heights allowed a reasonable validation of the simulation. As a result a synoptic picture of the tsunami impact is provided over the entire coast of Sri Lanka. It is found that amplification due to shoaling applies mainly in the Eastern and Southern coast because, here, the wave is propagating across the sea floor slope, while propagating along the slope for the Western coast. Spots of high waves are due to wave focusing in some coastal areas while local submarine canyons in other areas inhibit the wave amplification. Ó 2010 Elsevier Ltd. All rights reserved.
Keywords: 2004 tsunami Sumatra Runup Numerical simulation Sri Lanka
1. Introduction On Sunday, 26 December 2004 at 00:58:53 UTC, the Mw ¼ 9.1e9.3 SumatraeAndaman earthquake triggered a large tsunami that severely damaged coastal communities in countries along the Indian Ocean, including Indonesia, Thailand, Sri Lanka, India, Maldives, and Somalia located at more than 6000 km, which corresponds to a wave travel time up to 8 h. The epicenter and first aftershocks indicate that approximately 1300 km of subduction interface ruptured along the northern Sunda Trench. The tectonics and seismicity of the bay of Bengal can be found in Socquet et al. (2006), Stein and Okal (2005) and the recurrence of earthquaketriggered tsunamis is discussed in Altinok and Ersoy (2000). As far as the tsunami is concerned, Sri Lanka was the second area after Sumatra, Indonesia, that encountered most damages and casualties. A series of waves impacted more than two third of the coast. The tsunami left e31,000 human deaths, e4110 missing, 440,000 displaced. It provoked severe damages in harbours, infrastructures including housing, livelihoods. Communication and transportation systems were seriously affected. A great level of pollution was
* Corresponding author. E-mail address: [email protected] (M. Ioualalen). 1364-8152/$ e see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.envsoft.2010.04.010
observed in the flooded areas (e.g., in drinking water reservoirs). Naturally, the effects were increased or reduced by the shape of the coastline, geomorphologic features and the bathymetric variation, i.e., the physical processes involved. But also and importantly, the effective effects of the hitting were different depending on the level of exposure and vulnerability of the coastal communities to the impact of waves. Such parameters may also vary at a daily time scale, e.g., transport- traffic, rush hours. The measurements associated to the physical phenomena, the runup and inundation mainly, cannot be taken as a measure of the impact of tsunami because there are other facts that influences in the effects. However, there are cases where the impact magnitude is related strongest to these values. For example, two of the highest values for the height of wave and inundation were observed in Kalmunai and Hambantota (Liu et al., 2005; Fig. 1; Table 1); in the latter location, the number of victims reached 1000 and, as resulted of flooding, the lagoon near to the town was contaminated with debris swept for the tsunami inland. In Kalmunai, the wave flooded and caused damages as far 1500 m inland, within this zone all structures was severely damaged (Liu et al., 2005). In other places like Payagala, were the wave heights were weaker (Liu et al., 2005; Fig. 1; Table 1), the amount of fatalities was significantly high, 800 deaths. Here, a train with numerous passengers happened to circulate along the coast at the time the tsunami reached it. Also, another important feature is the level of protection of a coastal area
M. Ioualalen et al. / Environmental Modelling & Software 25 (2010) 1874e1880
1875
Fig. 1. Plotted measured and simulated wave height at locations mentioned in Table 1. GI and CH stand for the simulations based on the coseismic source of Grilli et al. (2007) and Chlieh et al. (2007) respectively. Isobathes 1000 m (brown curve), 500 m (green) and 300 m (blue) are reported. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
by natural barriers, like coral reef, mangroves and sand dunes, whatever is the wave height. Contrary to the southern and eastern coastline, in the eastern coast, where none of them exist, the damages were more important. In that area, poor housing situated over the coastline were converted by the tsunami in debris that yielded the re-enforcement of the damages. The numerical simulation that we present here for the Sri Lanka case study has two objectives: The first purpose is to provide
a synoptic picture of the event and, in particular, the runup or wave height maxima distribution along a coastline. It is suitable that the numerical simulation is not constrained by runup observations in order to assess a predicted runup distribution (this is the case in the present study). The other interest of numerical modeling is to identify physical processes that are responsible for local wave amplification or attenuation. Then the picture reveals vulnerable areas as well as sheltered ones. These issues are crucial in regions
1876
M. Ioualalen et al. / Environmental Modelling & Software 25 (2010) 1874e1880
Table 1 Estimates of wave height maxima relative to the mean sea level along the Southern coast Sri Lanka sampled by Liu et al. (2005) for the 26th December 2004 Indian Ocean tsunami along with simulated ones at the nearest locations of the computational domain (not necessarily the highest simulated wave) for 3 grid spacing: 0.2780 (in bold), 0.5560 and 1.1120 . SGI stands for Grilli et al. (2007) and Ioualalen et al. (2007a) tsunamiederived source and SCh is for the one proposed by Chlieh et al. (2007). Location
Long. ( E)
Lat. ( N)
Wave height max. (m) Obs.
SGI 0.2780
SGI 0.5560
SGI 1.1120
SCh 0.2780
SCh 0.5560
SCh 1.1120
(1) Trincomalee Hotel (2) Trincomalee town (3) Trincomalee military ChcPt. (4) Trincomalee HI Hostel (5) China Bay Stop 1 (6) China Bay Stop 2 (7) China Bay Stop 3 (8) N of Bathticaloa (9) Kattativu (10) Karativu (11) Ninto (12) Nalaveli Hostel Pottuvil (13) Ibral Nagar Nalaveli 1 (14) Ibral Nagar Nalaveli 3 (15) Ibral Nagar Nalaveli 4-5 (16) 2 Km SE of Kuchchaveli (17) Mankeri (18) Kalmurai Kuddi 1 (19) Kalmurai Kuddi 2 (20) Moratuwa (21) Koralawella (22) Wadduwa (23) Hambantota (24) Nonagama (25) Tangalla 1 (26) Kamburugama (27) Weligama (28) Galle (29) Dodanduwa (30) Hikkaduwa (31) Thiranagama (32) Galbokka (33) Seenigama (34) Dehiwala (35) Panadura (36) Pinnatara (37) Kalutara (38) Payagala (39) Yala (40) Boosa
81.21845 81.24202 81.23478 81.21380 81.19094 81.19206 81.26444 81.69267 81.74025 81.85432 81.86083 81.18847 81.21794 81.21918 81.21674 81.12103 81.48953 81.84164 81.83044 79.88353 79.88879 79.92110 81.12752 80.98835 80.79562 80.49195 80.44682 80.24915 80.14692 80.10413 80.12357 80.03091 80.08908 79.85640 79.90334 79.91306 79.94800 79.97835 81.25503 80.09014
8.618075 8.561729 8.563063 8.579356 8.503758 8.492397 8.463147 7.744343 7.686330 7.364405 7.343535 8.706572 8.660616 8.660272 8.661204 8.790385 8.013957 7.405348 7.422988 6.762450 6.749667 6.673167 6.128450 6.093750 6.029367 5.940050 5.968984 6.009734 6.083717 6.127517 6.110650 6.323317 6.166100 6.877667 6.715483 6.689667 6.608450 6.521217 6.166390 6.047760
4.4 Mean: 4.30 2.5 3.2 2.8 3.35 3.25 2.7 3.7 4.85 4.5 4.1-6.1 4.65 4.5 Mean: 4.13 3.35 5.5 6.2 7 3.56 4.55 3.61 6.1 8.71 3.24 2.4 2.73 5.24 3.6 4.2 4.55 4.26 5.05 3.48 4.24 4.15 3.82 5.04 4.65 Mean: 2.00
4.35 3.20 3.29 3.00 3.80 1.95 3.84 3.70 4.77 10.43 10.09 3.69 3.02 3.02 3.02 3.46 7.11 9.80 9.75 2.82 2.71 4.15 5.24 6.51 5.23 3.78 5.52 5.62 4.64 5.70 5.18 6.77 6.17 2.59 2.20 3.50 3.88 6.28 4.93 2.71
4.22 3.33 3.33 3.69 2.01 2.34 3.64 3.86 4.78 11.41 9.17 2.7 2.59 2.59 2.59 3.39 4.33 9.37 10.14 2.22 2.22 2.86 6.93 7.40 6.01 4.15 4.99 3.78 4.29 2.59 1.56 6.04 6.07 2.63 3.26 2.85 3.89 5.17 3.68 2.73
3.46 3.41 3.41 3.17 2.23 2.23 2.63 3.77 4.64 9.37 8.99 4.46 2.61 2.61 2.61 3.10 8.62 9.27 8.29 3.90 4.00 4.21 4.43 5.86 4.89 3.71 5.96 5.00 6.38 7.28 6.72 9.68 9.75 1.70 3.07 3.21 4.46 6.46 4.76 3.55
4.21 3.11 3.4 2.82 5.66 5.81 5.17 3.30 4.80 6.82 6.44 3.67 3.20 3.20 3.20 4.43 7.71 5.19 4.80 1.88 1.63 3.09 5.11 6.76 3.83 2.90 6.43 4.45 2.84 4.05 3.71 4.22 3.88 1.38 2.15 2.54 3.14 2.77 3.96 1.90
3.85 4.10 4.10 4.85 7.59 8.71 7.56 4.22 4.72 5.78 6.41 3.33 2.79 2.79 2.79 4.56 7.61 5.55 5.45 2.50 2.30 3.36 4.35 5.86 3.04 2.85 6.63 4.46 3.21 4.61 3.81 4.26 4.08 1.49 2.25 3.18 2.76 2.89 3.09 3.75
4.62 3.78 3.78 3.56 3.95 3.06 4.45 3.38 4.80 5.59 5.25 4.73 2.71 2.71 2.71 3.91 9.33 4.44 5.74 3.01 2.83 3.43 4.44 5.31 3.78 3.00 4.79 4.05 3.88 3.67 3.39 6.31 4.76 1.30 2.29 3.04 3.04 4.12 4.74 5.11
where tsunami events are frequent. Such methodology has been applied successfully by Ioualalen et al. (2007a) for the same tsunami event but on its impact on the Andaman coast of Thailand. The runup/wave height data set that is used here for validation of the numerical simulation is the one gathered by Liu et al. (2005) immediately after the event during the period 9e15 January 2005 (Table 1). Most of their sampling concerned wave height maxima (WHM) that have been derived from debris in trees and watermarks on buildings. Both the observed and our simulated WHMs are referenced to the mean sea level: the comparison between the two data sets is then direct because in both cases we have the freetide tsunami signal. 2. Numerical simulation of the tsunami for Sri Lanka 2.1. The numerical procedure We used here FUNWAVE tsunami propagation and runup/ inundation model. The model is fully nonlinear and dispersive, retaining information to leading order in frequency dispersion O [(kh)2] and to all orders in nonlinearity a/h (where k denotes an inverse wavelength scale, a denotes a wave amplitude, and h denotes a water depth) (Wei and Kirby, 1995; Wei et al., 1995). The
model treats the entire computational domain as an active fluid domain by employing an improved version of the slot or permeable-seabed technique, i.e., the moving shoreline algorithm proposed by Chen et al. (2000) and Kennedy et al. (2000) for simulation of runup. The model includes bottom friction, energy dissipation to account for the wave breaking and a subgrid turbulence scheme. We have constructed a 0.278 arc minute grid-spacing for the computational domain (e520 m) ranging from 2 N to 13 N and from 77.2 E to 96 E (Fig. 2). The domain is composed of 3958 2377 nodes and we chose an optimal time step of 0.5 s to avoid numerical instabilities and to limit truncation errors. A propagation time of 6 h has been chosen (43,200 time steps) to take into account the highest crest observed in Sri Lanka (Pattiaratchi and Wijeratne, 2009; Liu et al., 2005). Although, the grid spacing seems to be relatively large, Ioualalen et al. (2007a) used approximately the same one for the same event applied for the Thailand case study and they discussed the relevance of the parameter for this specific event. They obtained a satisfactory comparison with observations, e.g., cross-correlation between observed and simulated runup were 0.94 and the RMS difference was only about 17%. Three sets of data have been used: Offshore the Islands 2arc minute ETOPO-2 bathymetry has been
M. Ioualalen et al. / Environmental Modelling & Software 25 (2010) 1874e1880
1877
Fig. 2. Computational domain including the vertical sea floor deformation at 1 m contour intervals predicted from the coseismic SumatraeAndaman model G-M9.15 of Chlieh et al. (2007). Dashed lines indicate subsided area and the continuous lines the uplifted regions. The background bathymetry is plotted in grey at 500 m contour intervals.
used until depth 200 m. Then for shallow water digitalized marine charts data have been used. Inland, the SRTM DEM, with a horizontal resolution of 90 m, have been used.
2.2. The initial wave solution As input in our numerical modeling, we used two already validated seismic sources to compute the vertical seafloor deformation associated with the SumatraeAndaman earthquake: The first one, based on GPS network, has been proposed by Chlieh et al. (2007) and the second, developed by Grilli et al. (2007) and Ioualalen et al. (2007a), has been calibrated with available tide gauges records in the Indian Ocean and with Jason sea level anomalies data. The source G-M9.15 developped by Chlieh et al. (2007) integrates observations from several geodetic observations (Figs. 2 and 3): a near-field survey-mode GPS sites on northern Sumatra, several reoccupied GPS sites in the Andaman and the Nicobar islands, a set of more distant continuous GPS stations in Thailand and Malaysia, measured patterns of coral reefs vertical motion and the hinge-line deduced from satellite imagery. The joint inversion of all these geodetic data provided good constraints on the coseismic source. This source has been shown to be in good agreement with many others geophysical and hydrographic data, i.e, seismic observations, normal modes up to 300 s T-wave recorded from hydroacoustic arrays and sea level variations monitored by altimeters. Grilli et al. (2007) and Ioualalen et al. (2007a) calibrated a more energetic Mw ¼ 9.25 five-segment Okada (1985) dislocation source (Fig. 3; Table 2) using available hydrodynamic data, in particular the sea level anomalies recorded by the JASON-1 altimeter, which happened to transit over the area spawned by the spreading tsunami about 2 h after the earthquake initiation, during cycle 109 of its pass 129 (Smith et al., 2005), and several digital tide gauge records installed along the Gulf of Bengal coasts (Merrifield et al., 2005). This source has been tested and validated by Ioualalen et al. (2007a) for the Thailand case study. They were able to reproduce most of the runup features over the Andaman coast of Thailand. The runups obtained in a predictive mode (the simulations were not constrained by runup observations) fitted observations at a relatively high degree of accuracy. The two coseismic sources differ significantly in terms of seismic moment amplitude and also in terms of vertical deformation
distribution along the subduction zone. The gaps are mainly due to the respective modes of calibration: The source of Chlieh et al. (2007) has been constrained more objectively because it relies on effective GPS sea floor observation. The one of Grilli et al. (2007) has been built through hydrographic observations and thus only integrates the wave mode of propagation only. Besides it is built within Okada (1985)’s idealized assumptions. Ideally the two sources should converge to each other but it is not the case because (1) the functions of transfer between the floor deformation and the initial wave are not represented in the simulation (the transient movement within the water column) and (2) the hydrographic source highly depends on the quality of the ensemble of numerical simulations that have been performed to derive it.
2.3. The numerical results and discussion We perform a numerical simulation to provide a synoptic picture of the event along the entire coast of Sri Lanka and to analyze the features that are responsible for their coastal vulnerability or sheltering. To reach this objective we validate first our wave height results with available observations (Table 3, Fig. 1). Elementary statistics on simulations skills for the 0.278 arc minute grid-spacing simulation indicate that the numerical simulations based on Grilli et al. (2007) and Ioualalen et al. (2007a) and Chlieh et al. (2007) co-sesimic sources (named SGI and SCh herafter respectively) seem to give complementary results: SCh shows a better-fit for wave heights (mean, std and the rmse, i.e., the model skill) while GI exhibits a better cross-correlation with observations (Table 3). The amplitudes obtained with simulation SCh (6% of mean under-prediction and 18% of std over-prediction) are quite satisfactory which indicates that the seismic momentum and spatial distribution of the rate of subsidence/uplift of the Chlieh et al. (2007)’s source fits the tsunami wave heights requirements. The relatively poor correlation with observations (0.38) might be explained by the fact that the initial wave is triggered at once with no lag along the ruptured fault. As a reference, Ilayaraja et al. (2009) using the same source of Chlieh et al. (2007) for the same tsunami event and the case study of Andaman and Nicobar Islands found a cross-correlation of 0.88. For that particular case the Islands are located within the ruptured area and thus the risk of tsunami phase-lag with observation is considerably reduced (not the case for the remote Sri Lanka). We may then expect that a proper triggering process of the initial wave may improve substantially our
1878
M. Ioualalen et al. / Environmental Modelling & Software 25 (2010) 1874e1880
Fig. 3. Vertical sea floor deformation at 1 m contour intervals predicted from (left panel) SCh, the coseismic SumatraeAndaman model G-M9.15 of Chlieh et al. (2007) and (right panel) SGI, the solution of Grilli et al. (2007) and Ioualalen et al. (2007a) which Okada (1985)’s dislocation parameters are provided in Table 2. Dashed lines indicate subsided area and the continuous lines the uplifted regions. The background bathymetry is plotted in grey at 500 m contour intervals.
results. Such parameter has been tested by Wang and Liu (2006) who studied both effects of and impulsive and a transient fault models for the 2004 Indian Ocean tsunami. This is visible when we Table 2 Pure thrust reverse fault solution dipping eastward (Grilli et al., 2007 and Ioualalen et al., 2007a). Input parameters for Okada (1985) dislocation method (first 5 lines) and outputs (last 4 lines) for 5 tsunami source segments (Fig. 3): time delay of segment rupture from earthquake time s (a 60 s rising time is added); longitude and latitude of segment centroid (xo, yo); the centroid depth is d ¼ 25 km for all segments, the fault strike angle 4 (clockwise from North); the fault rake angle is l ¼ 90 for all segments (counterclock-wise from strike); the fault dip angle is d ¼ 12 everywhere (positive from the horizontal plane); the maximum fault slip D; the segment length along and width across (L, W); and the medium shear modulous taken m ¼ 4 1010 Pa for all segments; the seismic moment Mo; the characteristic initial tsunami wavelength lo and period so; and the characteristic tsunami trough and peak amplitudes ho. Note that, in the simulation, slip is highest at the segments’ centroids and drops by 50% at a radius of L from it. The total seismic moment of all 5 segments is Mo ¼ 7.55 1022 or Mw ¼ 9.25. Parameters S1
S2
S3
S4
S5
s (s)
272 93.90 E, 5.22 N 348 23 150, 130 1.58 1022 130 17.46 3.84, þ8.59
588 93.21 E, 7.41 N 338 12 390, 120 2.05 1022 120 23.30 2.33, þ4.72
913 92.60 E, 9.70 N 356 12 150, 95 0.61 1022 95 18.72 2.08, þ4.49
1273 92.87 E, 11.70 N 10 12 35 , 95 1.46 1022 95 18.72 2.31, þ4.60
xo, yo
4 D (m) L, W (km) Mo (J) lo (km) so (min.) ho (m)
60 94.57 E, 3.83 N 323 18 220, 130 1.85 1022 130 24.77 3.27, þ7.02
consider the simulation SGI: With relatively poor amplitudes results we obtain a better correlation better than for simulation SCh probably because we took into account some mode of initial wave triggering including time lags between the successive components of the segmented ruptured area (Table 2). The wave height overprediction of the simulation SGI is probably due to the fact that the authors calibrated their source with available tide gauges to the east with a 1 arc-minute-only grid-spacing for their computational grid. It is fair to say that their source would have been much improved if they used a finer gridding: Ioualalen et al. (2007a) showed for the Thailand case study a significant improvement of the simulated runup when from a 1arc minute grid spacing (mean over-prediction of 34%) to a 0.25arc minute (like here) (mean overprediction of 17% only) for the same GRI triggering source. However, since the FUNWAVE version we use does not allow multiple-grids we encounter difficulties to compute a finer gridspacing because of computing limitations. Since our validation conclusions are reasonable but less satisfactory than those of Ioualalen et al. (2007a) for the Thailand case study we need to perform sensitivity tests on the grid-spacing in order to estimate the convergence of our simulation and at least say if the simulations are in a coherent trend. For that purpose we have performed two simulations for each tsunamis source SGI and SCh with gridspacings 0.558 (e1040 m) and 1.112 (e2080 m) arc minute. Thanks to these grid tests, the basic statistics indicate that the SGI simulation is in a better convergence trend (Table 3): as far as the amplitude is concerned (std and rmse), the SGI simulations indicate a fair improvement when improving the grid spacing while no
M. Ioualalen et al. / Environmental Modelling & Software 25 (2010) 1874e1880
1879
Table 3 Comparison between the observed (OBS) and simulated runup or wave height maxima for coseismic sources SGI and SCh used for grid spacings 0.2780 (in bold), 0.5560 and 1.1120 . The basic statistics are the mean, the standard deviation, the root mean squared error (RMSE) (and their deviation in % from the observations), the cross-correlation qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P 2 coefficient and the norm L2. The usual RMSE representing the model skill is defined as i ðroi rsi Þ =n, where roi is the observed runup for each (i)-location of Table 1, rsi is for qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P 2 P 2 2 the associated simulated one and n is 40. The Norm L is defined as i ðroi rsi Þ = i ðroi Þ
Mean (m) STD (m) RMSE (m) Cross-correlation Norm L2
OBS
SGI 0.2780
SGI 0.5560
SGI 1.1120
Sch 0.2780
Sch 0.5580
Sch 1.1120
4.23 1.28
4.79 (D13%) 2.17 (D69%) 1.89 (45%) 0.55 0.43
4.37 (þ3%) 2.32 (þ81%) 1.98 (47%) 0.52 0.45
4.95 (þ17%) 2.32 (þ81%) 2.17 (51%) 0.47 0.49
3.99 (-6%) 1.51 (D18%) 1.58 (37%) 0.38 0.36
4.29 (þ1%) 1.66 (þ30%) 1.99 (47%) 0.10 0.45
4.07 (4%) 1.33 (þ4%) 1.48 (35%) 0.37 0.33
convergence trend appears for SCh ones. We obtain the same conclusion for the special distribution of the simulated runup: we improve the cross-correlation with finer grid spacing while no trend appears for SCh. We also tend to minimize the Norm L2 with SGI while there is no such feature for SCh. As a consequence, although the 0.278 arc minute grid, for SGI do not provide excellent results, we may assert that we place in a reasonable convergence trend and we may select this simulation for further process studies. This sensitivity test allowed us to select in a robust way the more coherent simulation between SGI and SCh while it is difficult to do so with only the absolute values of the 0.278 arc minute grid simulations. It is also fair to say that despite the relatively scarce grid spacing, the inherent uncertainties on the data used for the gridding (for example SRTM topography data may introduce errors on the slope of the beach and thus on the wave height), the possible gap between the positioning for observations and simulation locations (we compared heights at same locations or nearest point), we obtained a relatively coherent simulation using SGI triggering source. Based on that reasonable validation of the SGI simulation, we propose now a synoptic picture of the event (wave height distribution along the entire coast of Sri Lanka) and identify the physical processes that are responsible for the disparities. We propose here only the vulnerability picture regards to the physical event. Obviously, we do not discuss other vulnerability parameters suggested in the introduction (natural barriers, socio-economic features etc). Fig. 4 shows the predicted runup along entire coast of Sri Lanka. We use here the term “prediction” because the simulation has not been constrained by any observed runup or wave height values so that we place in a prognostic mode. In our case study, the SouthEastern coast of India exhibits a relatively pronounced coastal variability on tsunami impacting with clear risk areas and more sheltered ones for any tsunami event. This even can be considered as the worst case scenario for Sri Lanka because the threatening recurrent earthquakes happen in the considered subduction area and also because the Mw ¼ 9.1e9.3 event is extreme (based on an extreme size of the ruptured area). For that simple reason, the study of this event can be considered as a relevant for tsunami threat mapping. In the following we try to identify the features that are responsible for the runup distribution along the coast of Sri Lanka. It is obviously not the objective of this work to provide a GIS mapping of the tsunami risk. The wave amplification factors are relatively easy to identify. First the shoaling applies everywhere in the Eastern and Southern coasts where the wave reaches the coast at a pronounced incidence (Fig. 4) and propagates across a slope. This is not the case for the Western coast where the wave propagates along the slope, especially North of 7 N. Here, the flux of wave energy is conserved along a quasi-constant slope and thus does not yield a slope (or Green’s) effect.
Fig. 4 exhibits areas of wave refraction/focusing at the Eastern coast where spots of very high waves appear (noted FC in the Figure with red arrows). This is the case at around 9.2 N, 8.0 N and 7.3 N where we observe convex-type isobaths. At around 8.0 N, although Google Earth indicates the presence of deserted beaches, this non-occupied area is likely to be vulnerable. Because of its poor density in population, the 2004 event did not generate significant damages in that area. However the period of recurrence of such event (centenaries time scale), future development plan of the area should take into account the potential threat.
Fig. 4. Simulated maxima of wave height (in meters) recorded over the entire time integration for simulation SCh. Isobaths 30, 100 and 500 m are represented in while curves. FC stands for “focusing” and UC for underwater canyons. Blue arrows indicate wave divergence and red ones correspond to wave convergence (focusing). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
1880
M. Ioualalen et al. / Environmental Modelling & Software 25 (2010) 1874e1880
Fig. 4 also indicates areas where underwater canyons are present (noted UC in the figure with blue arrows). At their termination (the coast line), the wave is diverging (concave-type isobaths) and escapes from its sides (a defocusing). The effects of these canyons are to inhibit the wave amplification and eventually balance shoaling effects. This has been shown already for the Bangladesh coast and the same event where a multitude of underwater canyons have defocused the wave which encountered highest waves only offshore at the continental shelf (Ioualalen et al., 2007b). This explains why the coast of Bangladesh has been so little affected by the event. Acknowledgements The authors would like to thank the French Agence Nationale Pour la Recherche, ANR, for funding this work through the grant TSUMOD ANR-05-CATT-016-02. M. Ioualalen addresses his thanks to the Department Milieu et Environement of the Institut de Recherche pour le Développement, IRD, for having granted him a long stay at INOCAR, Ecuador, especially J. Boulegue and P. Soler and P. Charvis and J.-Y. Collot (Géosciences Azur) for their constant support. References Altinok, Y., Ersoy, S., 2000. Tsunamis observed on and near the Turkish coast. Nat. Hazards 21, 185e205. Chen, Q., Kirby, J.T., Dalrymple, R.A., Kennedy, A.B., Chawla, A., 2000. Boussinesq modeling of wave transformation, breaking, and run-up. II: 2D. J. Wtrwy, Port, Coast, Oc. Engrg., ASCE 126 (1), 48e56. Chlieh, M., Avouac, J.-P., Hjorleifsdottir, V., Song, T.-R.A., Ji, C., Sieh, K., Sladen, A., Hebert, H., Prawirodirdjo, L., Bock, Y., Galetska, J., 2007. Coseismic slip and afterslip of the great Mw 9.15 SumatraeAndaman earthquake of 2004. Bull. Seism. Soc. Amer. 97 (1A), S152eS173. Grilli, S.T., Ioualalen, M., Asavanant, J., Shi, F., Kirby, J.T., Watts, P., 2007. Source constraints and model simulation of the December 26, 2004 Indian Ocean tsunami. J. Wtrwy, Port, Coast, Oc. Engrg., ASCE 133 (6), 414e428.
Ilayaraja, K., Ioualalen, M., Chlieh, M., Krishnamurthy, M.V., Ramana Murthy, M.V., 2009. Numerical modeling of the 26th December 2004 India Ocean tsunami at Andaman and Nicobar Islands. J. Coast. Conserv.. doi:10.1007/s11853-008-00338. Ioualalen, M., Asavanant, J., Kaewbanjak, N., Grilli, S.T., Kirby, J.T., Watts, P., 2007a. Modeling of the 26th December 2004 Indian Ocean tsunami: case study of impact in Thailand. J. Geophys. Res. Oceans 112, C07024. doi:10.1029/ 2006JC003850. Ioualalen, M., Pelinovsky, E., Asavanant, J., Lipikorn, R., Deschamps, A., 2007b. On the weak impact of the 26 December Indian Ocean tsunami on the Bangladesh coast. Nat. Hazards Earth Syst. Sci. 7, 141e147. Kennedy, A.B., Chen, Q., Kirby, J.T., Dalrymple, R.A., 2000. Boussinesq modeling of wave transformation, breaking, and run-up. I: 1D. J. Wtrwy, Port, Coast, Oc. Engrg., ASCE 126 (1), 39e47. Liu, P.L.-F., Lynett, P., Fernando, H., Jaffe, B.E., Fritz, H., Higman, B., Morton, R., Goff, J., Synolakis, C., 2005. Observations by the international tsunami survey team in Sri Lanka. Science 308 (5728), 1595. doi:10.1126/science.1110730. Merrifield, M.A., Firing, Y.L., Aarup, T., Agricole, W., Brundrit, G., Chang-Seng, D., Farre, R., Kilonsky, B., Knight, W., Kong, L., Magori, C., Manurung, P., McCreery, C., Mitchell, W., Pillay, S., Schindele, F., Shillington, F., Testut, L., Wijeratne, E.-M., Caldwell, P., Jardin, J., Nakahara, S., Porter, F.-Y., Turetsky, N., 2005. Tide gauges observations of the Indian Ocean tsunami, December 26, 2004. Geophys. Res. Lett. 32, L09603. doi:10.1029/2005GL022610. Okada, Y., 1985. Surface deformation due to shear and tensile faults in a half-space. Bull. Seis. Soc. Am. 75 (4), 1135e1154. Pattiaratchi, C.B., Wijeratne, E.M.S., 2009. Tide gauge observations of 2004e2007 Indian Ocean tsunamis from Sri Lanka and Western Australia. Pure Appl. Geophys. 166, 233e258. Smith, W.H.F., Scharroo, R., Titov, V.V., Arcas, D., Arbic, B.K., 2005. Satellite altimeters measure tsunami: early model estimates confirmed. Oceanography 18 (2), 10. Socquet, A., Vigny, C., Simons, W., Chamot-Rooke, N., Rangin, C., Ambrosius, B., 2006. GPS determination of the relative motion between India and Sunda, and its accommodation in Myanmar. J. Geophys. Res. 111, B05406. doi:10.1029/ 2005JB003877. Stein, S., Okal, E., 2005. Speed and size of the Sumatra earthquake. Nature 434, 581e582. Wang, X., Liu, P.L.-F., 2006. An analysis of 2004 Sumatra earthquake fault plane mechanisms and Indian Ocean tsunami. J. Hydraulic Res. 44 (2), 147e154. Wei, G., Kirby, J.T., 1995. Time-dependent numerical code for extended Boussinesq equations. J. Wtrwy, Port, Coast, Oc. Engrg., ASCE 121 (5), 251e261. Wei, G., Kirby, J.T., Grilli, S.T., Subramanya, R., 1995. A fully nonlinear Boussinesq model for free surface waves. Part 1: highly nonlinear unsteady waves. J. Fluid Mech. 294, 71e92.