Assessment of tsunami threat to Sri Lanka from potential megathrust earthquakes in the Arakan subduction zone

This paper describes a numerical study carried out to assess the threat posed to Sri Lanka by potential tsunamigenic earthquakes in the Arakan subduction zone in the northeastern Bay of Bengal. The fault plane model adopted in the present simulations corresponds to a moment magnitude of M w = 8.8, which may be considered as a worst-case scenario. A hydrodynamic model based on linear shallow-water equations was employed to simulate tsunami propagation from the source to the shoreline around Sri Lanka. The model results extracted at an average water depth of about 5 m off the coastline were processed to obtain the spatial distribution of the maximum value of the ‘tsunami amplitude’ (i.e., the height of the crest of the tsunami waves above mean sea level) as well as the arrival time contours for the leading wave. The numerical simulations suggest that the coastal zone of the Northern Province is at most risk from a tsunami generated in the Arakan fault plane with maximum tsunami amplitudes of the order of 10 m, whilst the maximum amplitudes along the coastal belt of the Eastern and Southern Provinces could reach up to about 5.6 m and 2.6 m, respectively. On the other hand, the coastline in the Western Province in the shadow zone of direct tsunami impact will only receive waves of small amplitude up to 0.6 m. The numerical results also indicate that the tsunami waves will first hit the east coast about 160 minutes after the earthquake, followed by the coastal belts of the Northern and Southern Provinces. The information presented in this paper relating to likely tsunami amplitudes and arrival times around the coastline of Sri Lanka would help authorities responsible for evacuation to make better judgement as to the level of threat in different areas along the coastline, and to act accordingly, if a large earthquake as to occur in the Arakan subduction zone in future.


INTRODUCTION
The massive tsunami on 26 December 2004, generated by the mega-thrust earthquake of moment magnitude M w = 9.1-9.3 1 in the Andaman-Sunda trench, created an awareness of the potential for destructive tsunamigenic earthquakes to occur in the active subduction zones around the Indian Ocean.The December 2004 earthquake ruptured a 1200 km long segment of the Andaman-Sunda trench, starting near the northern part of the Sumatra Island and ending offshore of the Andaman Islands 1 (Figure 1).This massive earthquake also likely increased the stress on adjacent segments of the subduction zone, raising the seismic hazard at both ends of the rupture, i.e., further south in the Sunda trench as well as further north in the Arakan trench in the northeastern Bay of Bengal along the coast of Myanmar [2][3][4] .The M w = 8.5 Nias earthquake of 28 March 2005 that occurred immediately south of the epicentre of the December 2004 earthquake could be another example of such triggered seismicity 5 .
Recent GPS (Global Positioning System) velocity field data and past seismicity confirm that the Arakan trench, where the Indian plate underthrusts the South-East Asian plate, is still an active subduction zone [5][6][7][8] .Past seismic activity in the Arakan trench includes historical reports of large earthquakes in 1762 and in 1878 causing significant uplifts of 3-7 m and 6 m, respectively 9,10 .In a recent article, which received worldwide attention, Cummins 11 concluded that there is a high risk of occurrence of a giant tsunamigenic earthquake in the Arakan trench with potential for causing great loss of lives and destruction.His assessment was based on a detailed examination of the tectonic environment, stress and crustal strain observations as well as historical earthquake activity in the Arakan subduction zone.
If such an earthquake causing a tsunami were to occur, the location and the orientation of the Arakan fault plane such that much of the north and east coasts of Sri Lanka would likely be exposed to the path of direct wave energy from the source.The low-lying coastal terrain with high population density in many areas of the north and east could be highly vulnerable to such an event depending on the amplitude of the tsunami waves (i.e., the height of the crest of the tsunami waves above mean sea level) approaching the shore.
A coarse-grid model simulation carried out for the entire region 11 in general indicated significant offshore tsunami amplitudes in deep-water off the north and east coasts of Sri Lanka.However, the bathymetry of the shallow continental shelf around Sri Lanka could amplify these offshore wave amplitudes several times as a result of shoaling and energy focusing due to refraction.This paper presents a high-resolution numerical study of tsunami propagation, carried out in deep-water as well as over the shallow continental shelf, to assess the potential threat posed to the coastal belt of Sri Lanka by a tsunamigenic earthquake in the Arakan subduction zone.

METHODS AND MATERIALS
In the following sections, the fault plane model and the hydrodynamic model employed to simulate the generation and the propagation of a tsunami that could be caused by a shallow, submarine earthquake in the Arakan subduction zone were outlined.

Fault Plane Model
The fault plane model of Cummins 11 corresponding to a mega-thrust earthquake of moment magnitude M w = 8.8 is adopted in the present study as the worst-case scenario for the Arakan coast.The fault line is oriented so that its upper edge lies along the deformation front (Figure 1) and a 10 m reverse slip results in the observed pattern of uplift and subsidence for the Arakan earthquake of 1762 mentioned earlier 9 .
The fault parameters proposed by Cummins 11 for the above earthquake are given in Table 1.Assuming that the sea surface follows the sea bed deformation instantaneously, we use Okada's dislocation model 12 to obtain the initial sea surface elevation for the above coseismic tsunami source (Figure 2).

Hydrodynamic Model
Grid set-up: A dynamically coupled system of two nested grids was employed to simulate the tsunami propagation from the Arakan subduction zone towards the shoreline of Sri Lanka.The bathymetry data for the largest grid employed in the simulations, i.e., Grid-1 shown in Figure 3, was obtained by interpolating GEBCO 13 (General Bathymetric Chart of the Oceans) data with a resolution of 1 arc-minute to a grid of 1.3560 arcminutes (~ 2500 m) spacing.Grid-2, which is embedded in Grid-1 for the simulation of tsunami propagation over the shallow continental shelf off Sri Lanka at a finer resolution of 0.2712 arc-minutes (~ 500 m), is also shown in Figure 3.The bathymetry for Grid-2 was at first interpolated from the GEBCO grid and was then updated with data from UK Admiralty charts.These navigation charts typically cover depths down to about 3000-4000 m at scales of 1:150,000 or 1:300,000.The nearshore  bathymetry at some localities was further updated with available data from higher resolution Admiralty charts at scales of 1:10,000 and 1:15,000.

Model formulation:
The mathematical model used in the present study was the Cornell Multi-grid Coupled Tsunami Model 14 (COMCOT, coded in FORTRAN 90), which solves the non-linear shallow water equations on a dynamically coupled system of nested grids using a modified leap-frog finite difference numerical scheme.This model had been validated by experimental data 15 and successfully used to investigate several historical tsunami events, including the 2004 Indian Ocean tsunami 14,16,17 .
In the nested grid system employed in COMCOT, the inner (finer grid) regions adopt a smaller grid size and time step and were nested inside an outer (larger grid) region.At the beginning of each time step along the interface of two different regions, the volume flux, which was the product of water depth and depth-averaged velocity was interpolated from the outer (larger grids) region into inner (finer grids) region.The water surface elevations and the volume fluxes in the inner (finer grids) region were calculated and the resulting free surface elevations were averaged to update those values in the larger grids, which overlaps the inner region.The volume fluxes in the outer (larger grids) region can also be updated.With this algorithm, features of tsunami dynamics over the shelf with a high spatial and temporal resolution could be captured and at the same time, a high computational efficiency could be maintained.
The amplitude of a trans-oceanic tsunami during its propagation was typically of the order of magnitude of 1 m, whilst the wavelength of the leading wave was of the order of magnitude of 100 km 17 .Thus, linear shallow water equations were adequate to solve tsunami propagation in Grids-1 and -2.
Furthermore, as the spatial extent was comparatively larger in Grids-1 and -2, the spherical coordinate system was used to solve the linear shallow water equations given in the following 14 : where, ζ was free surface elevation; (ψ, ϕ) denote the longitude and latitude of relevant location on the Earth; R was the Earth's radius; P and Q stood for the volume fluxes (P = hu and Q = hv , with u and v being the depth-averaged velocities in longitudinal and latitudinal directions, respectively); h was the still water depth; and f represented the Coriolis force coefficient.

RESULTS AND DISCUSSION
The numerical model results were processed to obtain the spatial distribution of the maximum values of the amplitude of the tsunami as well as the arrival time contours.The following results refer to a potential tsunamigenic earthquake of moment magnitude M w = 8.8 in the Arakan oblique subduction zone off the west-coast of Myanmar in the northeastern Bay of Bengal.
First, Figure 4 shows three 'snapshots' of the simulated tsunami wave front approaching Sri Lanka at:  coast of Bangladesh in the neighbourhood of the fault plane.The Andaman and Nicobar Islands to the south of the fault line too appeared to receive a significant amount of wave energy.However, the geographical location of the fault plane ensured that the entire west and northwest coasts as well as part of the south coast of Sri Lanka were in the shadow zone of direct impact of tsunami waves.
Prior information relating to the time it took for the leading wave of a tsunami to arrive at a given coastline was useful in early warning and emergency planning.Accordingly, the computed arrival time contours for the coastal belt of Sri Lanka corresponding to a tsunami generated by an earthquake in the Arakan trench were shown in Figure 5.Note that the arrival times given are in minutes after the occurrence of the earthquake and were based on the first 1 cm rise of the mean water level.As the tsunami propagation speed depends only on the water depth, these calculated arrival times ought to be applicable to a tsunamigenic earthquake of any magnitude in that part of the Arakan trench.slowed down the speed of the tsunami waves considerably.Furthermore, the south and the west coasts in the shadow zone of direct impact would receive tsunami waves, albeit of lesser magnitude as seen later, in about 180-210 minutes, and 210-240 minutes, respectively, after the earthquake.
In the following sections, the distribution of the maximum amplitude of the tsunami would be considered in detail.Figure 6 shows a shaded plot of the spatial variation of the maximum values of the tsunami amplitudes extracted from the computational domain 1 (grid spacing~2500 m).It could be seen that the generally shallower source bathymetry limits the amount of tsunami energy generated even though the moment magnitude of the earthquake was as high as M w = 8.8.The resulting tsunami was thus mainly confined to the Bay of Bengal with only marginal penetration to the Indian Ocean Basin.The distribution of computed tsunami amplitudes in the near-field also indicated that the amount of tsunami energy derived from the relatively deeper southern part of the fault plane was comparatively more than that from the shallower northern part.
Moreover, the location and orientation of the fault was such that, in the far field, the zone of maximum energy release was directed towards the southern part of the east coast of India between 10 o N~15 o N, and to a lesser extent, towards the north coast of Sri Lanka, particularly, the Jaffna peninsula (Figure 6).It also appears that, although oriented lateral to the zone of maximum energy release, some stretches of the east coast of India north of 15 o N too could be hit by significantly large wave heights owing to the effects of shoaling and energy focusing caused by the bathymetry.However, it was not intended to examine the tsunami amplitudes along the east coast of India in detail in the present analysis.This was because the model results for the shoreline off the east coast of India were available only from Grid-1 with a coarser spatial resolution; moreover, the bathymetry that was employed in the present study for the eastern seaboard was based only on those of 1 arc-minute resolution and was not updated with data from higher resolution navigation charts.
In order to have a closer look at the maximum water levels expected along the coastal belt, the computed maximum tsunami amplitudes at water points nearest the shoreline, viz., at an average water depth of about 5 m, were extracted from Grid-2 (grid spacing~500 m).Accordingly, Figure 7 shows the variation of the maximum tsunami amplitude with longitude ( o E) or latitude ( o N) along the various coastal sectors of the island.Also, the mean as well as the range of the computed tsunami amplitudes for each segment of the coast at provincial level are given in Table 2 whilst the same at the district level in Table 3.It is seen from Table 2 that the tsunami amplitudes are the largest along the northern coastline with a mean value of 4.3 m and a maximum of 10 m.The location and the alignment of the source as well as the effects of the shelf bathymetry were the primary reasons for the higher values of tsunami amplitudes here.Several prominent peaks in the variation of tsunami amplitude could be seen in this sector, as for instance, those at ~80.37 o E (A) and at ~80.8 o E (B) (Figure 7a).In order to further examine the factors contributing to the higher wave amplitudes at these locations, the wave fronts spaced at 2 minute time intervals showing the effects of refraction caused by the bathymetry are depicted in Figure 8.The computed tsunami amplitudes are also shown in this figure for convenience.Clearly, at both locations, A (~80.37 o E) and B (~80.8 o E), the peaks in tsunami amplitudes are due to the focusing of wave energy by the bathymetry as indicated by the arrows.The computed tsunami amplitudes along the eastern sector in Figure 7(b) also shows considerable spatial variation, ranging from 0.8 m to 5.6 m with a mean value of 2.4 m.Furthermore, the mean value of the tsunami amplitudes is seen to decay from north to south along the east coast.Table 3 indicates that, in the eastern sector, the coastal belts of Trincomalee and Batticaloa are likely to experience larger wave amplitudes of up to about 5.6 m than the Ampara coast where the highest amplitude is about 3.5 m.There are also several prominent troughs and peaks in the computed tsunami amplitudes shown in    3 indicates that the tsunami amplitudes along the shorelines of Galle and Matara districts are likely to be significantly less than those along the coast of the Hambantota district.The tsunami amplitudes on the western and northwestern coastlines are small, ranging from 0.1 m to 0.6 m with a mean value of 0.3 m.These peak values of the computed amplitudes along the north and east coasts due to a potential tsunami generated in the Arakan trench are about 50-60% of those due to the 2004 tsunami that originated in the Northern Sumatra-Andaman trench; however, the corresponding values for the west and the south coasts are only about 10% and 20%, respectively.
The computed tsunami amplitudes shown in Figure 7 indicate that it may not be necessary to evacuate the entire coastline of Sri Lanka if a large earthquake capable of generating a major tsunami were to occur in the Arakan trench.A warning requiring immediate evacuation could be issued to all coastal communities in the north and east, and probably in the southeast (Hambantota district) as well; however, given the comparatively smaller waves anticipated, a 'tsunami watch' would be sufficient, in general, for coastal areas of the west and northwest, and perhaps, also for the districts of Galle and Matara in the south.
It may be mentioned that tsunami warnings for evacuation had to be issued to the entire coastal belt of Sri Lanka owing to a lack of prior threat assessments, when two large earthquakes occurred on March 28, 2005 and on September 12, 2007 off the west coast of Sumatra.Had such threat assessments as the present study been available to the authorities, evacuation orders could have been limited only to areas at risk, thus ensuring optimum utilization of resources and minimizing the economic losses as well as potential injuries during evacuation.
It is emphasised that the numerical results such as those presented in Figure 7 only show which parts of the coastline of Sri Lanka are likely to be exposed to greater risk in terms of oncoming tsunami waves near the shoreline.However, factors such as the onshore topography, the population density and the construction standards could further enhance or reduce the vulnerability of a given community.For example, the elevation of onshore lands significantly influences the inundation pattern as well as the flow depth and the run-up heights whilst ordinary brick and mortar housing are more susceptible to tsunami loading than reinforcedconcrete buildings.
Furthermore, there are several limitations and approximations inherent in numerical modelling of tsunami.Since the initial condition for the modelling is determined by the displacement of the ocean bottom along the fault line, the largest source of errors is the earthquake model itself.Another significant limitation is that the resolution of the modeling is no greater or more accurate than the bathymetric data used.It must also be added that shallow water models such as that employed in the present study assume a uniform velocity profile across the flow depth and neglect vertical accelerations.
Finally, the next tsunamigenic earthquake in the Arakan subduction zone need not necessarily be as large as the M w = 8.8 used in the present simulations.If the magnitude of the next earthquake is smaller compared to that used in the present simulations, the tsunami amplitudes reaching the shoreline of Sri Lanka will be smaller than those computed; nevertheless the spatial variation of the tsunami amplitudes around the coastline is expected to be qualitatively similar.

CONCLUSION
A numerical model based on shallow-water equations were employed to evaluate the potential threat to Sri Lanka from a tsunamigenic earthquake of moment magnitude M w = 8.8 11 in the Arakan subduction zone off the west-coast of Myanmar.The tsunami propagation across the shallow and narrow continental shelf around Sri Lanka was modelled using a high-resolution innergrid to enable simulation of essential nearshore processes such as shoaling and refraction.
The numerical simulations indicated that, in the far-field, the east coast of India and the north and east coasts of Sri Lanka lie in the direct path of the tsunami wave front.Further, the geographical location of the fault plane ensured that the entire west and northwest coasts as well as part of the south coast of Sri Lanka were in the shadow zone of direct impact of tsunami waves.The model results also indicated that the tsunami waves would first hit the eastern shores about 160 minutes after the earthquake, followed by the northern and southern shorelines.
The maximum tsunami amplitudes at water points nearest the shoreline around the coastline of Sri Lanka were extracted in order to assess the spatial variation.The computed peak tsunami amplitudes along the northern shoreline ranged from 1.9 m to 10.0 m with a mean value of 4.3 m, whilst in the eastern sector, from 0.8 m to 5.6 m with a mean of 2.4 m.The tsunami amplitudes along the western and northwestern coastlines were small: the mean value was 0.3 m and the maximum was 0.6 m.
The information presented in this paper relating to likely tsunami amplitudes and arrival times around the country would help authorities responsible for evacuation to make better judgements as to the level of threat in different areas along the coastline, and to act accordingly, if a large earthquake were to occur in the Arakan subduction zone.

Figure 1 :
Figure 1: Active subduction zones in the Indian Ocean: Arakan, Andaman and Sunda Trenches

Figure 2 :
Figure 2: The initial sea-surface displacement calculated from the fault model

Figure 3 :
Figure 3: Grid-1 of the computational domain;Location of Grid-2 is also shown.

Figure 6 :
Figure 6: Distribution of computed maximum tsunami amplitude across Bay of Bengal

Figure 5 :
Figure 5: Computed tsunami arrival times for the coastal belt of Sri Lanka.The times given are in minutes after the occurrence of earthquake

Figure 7 (
b).The troughs at around 7.3 o N, 7.5 o N, 7.8 o N and 8.6 o N are all as a result of energy defocusing along the valleys of the submarine canyons at these locations.On the other hand, the ridges on either side of these canyons help focus tsunami energy, resulting in increased tsunami amplitudes seen at around 7.4 o N, 7.60~7.75o N and 7.9 o N in Figure 7(b).

Figure 7 :
Figure 7: Computed maximum tsunami amplitudes at water points nearest the shoreline: (a) Northern coastline, (b) Eastern coastline, (c) Southern coastline, and (d) Western coastline and part of North-Werstern coastline.

Figure 8 :
Figure 8: Pattern of tsunami wave fronts spaced at 2 min.time intervals off the north coast.Inset: computed tsunami amplitudes.

Table 1 :
Fault parameters adopted in the present study.

Table 2 :
Statistics of computed tsunami amplitudes along the coast at provincial level.

Table 3 :
Statistics of computed tsunami amplitudes along the coast at district level.