Demarcation of High Hazard Areas Based on Human Stability Considerations in Tsunami Overland Flow

On December 26, 2004, the coastal belts of Sri Lanka and several other countries bordering the Indian Ocean suffered massive loss of life and damage to property due to the tsunami unleashed by the third largest earthquake ever recorded. One way of mitigating potential loss of lives from a similar event in the future is through advance warning of an approaching tsunami and quick evacua­ tion of vulnerable coastal communities to safer areas. The detailed planning necessary for such evacua­ tion exercises requires information about the level of hazard in each locality, which is usually based on inundation maps showing the depth flooding. However, the ability of a person to survive a flood flow depends not only on the depth of inundation but also on the flow velocity. Accordingly, the present paper employs a numerical model to compute the temporal variation of both the flow depth and the flow velocity across the entire computational domain in three cities on the south coast of Sri Lanka as well as the results from an experimental study on the human stability in tsunami overland flow to develop a procedure to demarcate the areas that should be evacuated in case of a tsunami warning.


Introduction
On December 26, 2004, coastal belts of several Indian Ocean countries including Indonesia, Sri Lanka, India and Thailand suffered massive loss of life and damage to property due to the tsunami unleashed by the mega-thrust earth quake of Miv = 9.1 -9.3 in the Sumatra-Andaman subduction zone. In Sri Lanka, 13 of the 14 districts lying along the coastal belt were affected: the death toll was over 35,000 with 20,000 injured and about 100,000 housing units either completely or partially damaged leaving half a million people homeless and causing massive disruption to livelihoods [3]. The fisheries and tourism sectors were among the hardest hit with many fishing boats and beach front hotels damaged or destroyed. The enor mity of this unprecedented tragedy clearly showed the potential vulnerability of the coastal belt of Sri Lanka to the forces of nature, especially to the extreme geological and meteo rological happenings in the vast expanse of the Indian Ocean around the country.
The mega-tsunami of December 26, 2004 as well as the two subsequent major earthquakes on March 28, 2005 and on September 12, 2007 off Sumatra which led to tsunami warnings or alerts in Sri Lanka clearly underscored the need to have a proper system in place for tsunami early warning as well as for quick evacuation of vulnerable coastal communities to safer areas. One necessary requirement in planning-such evacuation is to identify which areas are not safe and which areas are safe. Often, such a classification of the tsunami hazard is based on inundation maps showing the depth of flood ing. However, the ability of a person to survive a flood flow depends not only on the depth of inundation but also on the flow velocity. Accordingly, the present paper employs a numerical model to compute the temporal variation of both the flow depth and flow veloc ity throughout the areas of interest, and then combines the model output with results of an experimental study of human stability in flood flows to delineate areas that are deemed unsafe. As a case study, we apply this method to three cities on the south coast of Sri Lanka, namely, Galle, Matara and Hambantota, devastated by the 2004 tsunami.

Tsunami Overland Flow
We briefly review here the main characteristics of the tsunami overland flow (Fig. 1).
The tsunami waves as they shoal into shallow water gain in height, steepen and either break before arriving at the shoreline or merely surge over the nearshore bathymetry without break ing, depending on the wave characteristics and the slope of the nearshore sea bottom. A brea king tsunam i wave front travels inland as a m m surging bore whilst a non-breaking wave will appear as a fast rising tide with a rapid increase of water level causing flooding onshore. This turbulent, fast-moving flow results in damage, collapse, or floating away of buildings and vegetation. People are drowned owing to the high water, the difficulty of withstanding the fluid forces, the difficulty in coping with the large turbulent eddies, or impact with debris. The amount of debris picked up by tsunami overland flow depends on the distance of incur sion, the nature of coastal construction and vegetation, as well as the characteristics of the wave. Inundation distance is the maximum horizon tal penetration of tsunami flooding in the direc tion normal to the shoreline whilst run-up is the difference between the elevation of maxi mum tsunami penetration and the sea level at the time of the tsunami. (Fig. 2). Theinundation distance and the run-up height of the overland flow depend on the wave character istics, the local bathymetry and topography, as well as the nature of any coastal structures and the type of land use. It must also be added that the water that ran up inland will then begin to run down (recede or drain away) to the sea under gravity causing further damage and carrying off loose objects and people in its path as it retreats. One factor that affects the dynamics of tsunami waves during both the run-up and run-down processes mentioned above is the bottom friction. Often, Manning's form of the bottom roughness coefficient is employed in character izing bottom friction in model simulations. The presence of buildings and other structures, if sufficiently strong, too retard the tsunami flood flow and help dissipate energy, however, inclu sion of such effects in two-dimensional depthaveraged model formulations remains a challenge.

Methodology
In the following, we describe the methodology adopted in demarcating the high risk areas based on human stability considerations in tsunami overland flow, i.e., during both the run-up and run-down processes, in the three cities under consideration, namely, Galle Matara and Hambantota.
The stability of a person caught in a tsunami overland flow depends upon the local hydrody namic characteristics of the flow as well as the physical characteristics of the person, such as weight, height and strength. The instability occurs when the hydrodynamic load on a person due to the tsunami overland flow exceeds the resistance offered by the person through the friction between the land surface and the person concerned. The other forces acting on the human body under the above conditions are the self-weight and the buoy ancy.
The following criteria (Table 1) may be deduced from the experimental results reported in [9] in connection with a series of tests carried out to examine the stability of people caught in currents. These tests have been carried out over a range of flow depths (r|) and velocities (U) using male persons standing in a steady flow. It must be mentioned that the conditions given in Table 1 should be considered only as first approximations employed owing to the lack of more comprehensive data covering all variables that could contribute to the stability of a human body against a current.
In the present study, we employ a numerical model based on non-linear shallow water equa tions to compute the flow depths (q) and flow velocities (U) at every time step across the areas of interest.
The question that arises then is what initial conditions, i.e., the earthquake source param eters, that we should adopt to compute the flow depth and velocities in order to assess the level of tsunami hazard based on human stability considerations in tsunami overland flow. , Unfortunately, the short length of the historic record of tsunami in the Indian Ocean does not permit a comprehensive probabilistic tsunami hazard assessment for Sri Lanka. However, in such situations, the alternative procedure commonly employed in Japan and in the USA is to make use of the worst credible tsunami event from the available historic earthquake and tsunami records as the basis for hazard assess ments. Accordingly, the tsunami that occurred in December 2004 can be considered as the worst case scenario for Sri Lanka, so the results reported in the present study are based on that event.

Numerical Modelling
In the following, we outline the fault plane model and the hydrodynamic model employed to simulate the generation and the propagation of the 2004 tsunami. The nested grid set-up used in the hydrodynamic model is also described.

Fault Plane Model
Seismic inversion models (e.g., [1]) suggest that the rupture propagated approximately north ward from the epicentre located at (95.51oE The fault plane model adopted in the present study is that proposed by Crilli et al. [4] who employed available seismic data as well as hydrodynamic data pertaining to the genera tion and propagation of the 2004 tsunami to iteratively develop a 5-segment co-seismic tsunami source for the December 26, 2004 event; see Fig. 3 and Table 2, respectively, for locations and source parameters of the 5 dislo cation segments.
The contours of initial surface elevation for the above co-seismic tsunami source computed using Okada's [8] dislocation model and reported in Grilli et al. [4] is shown in Fig. 4. Further, it is assumed here that the sea surface follows the sea bed deformation instanta neously. Note that, in Fig. 4, the continuous lines represent uplift and dashed lines repre sent subsidence, both at 1 m contour intervals in the range -5 m to +8 m. The background bathymetry is also plotted in grey at 500 m contour intervals.
In a subsequent study [5], the numerical simu lations based on this fault plane model showed good agreement with the measured run-up heights along the Andaman coast of Thailand. Table 2     A dynamically coupled system of nested grids was employed to simulate the tsunami propa gation from Andaman-Sumatra subduction zone towards Sri Lanka and the subsequent inundation of the onshore areas of the cities under study. The bathymetry data for the largest grid employed in the simulations, i.e., Region 1 shown in Fig. 5, was obtained by inter polating E T 0 P 0 2 data with a resolution of 2 arc minutes to a grid of 0.6765 arc minutes (-1250 m) spacing. Region 2, which is embedded in Region 1, is also shown in Fig. 5.

Hydrodynamic M odel Grid Set-Up
The bathymetry data for Region 2 with a grid resolution of 0.1353 arc minutes (-250 m) shown in Fig. 6 as well as for Region 3 (grid spacing 50 m) was at first interpolated from E T 0 P 0 2 data and was then updated with data from the available navigation charts. These navigation charts typically covered depths down to about 3000 m -4000 m at scales of 1:150,000 or 1:300,000. The nearshore bathym etry of some areas was updated with higher resolution navigation charts at scales of 1:10,000 and 1:15,000.
Three fourth level grids with a grid spacing of 10 m were used to cover the three cities (see Fig.  7 for location of this grid with respect to Region 3). The topographic data for this grid was obtained from high resolution LIDAR (Light Detection and Ranging) survey data (The precise digital earth model of the coastal areas of Sri Lanka, Project Director: Prof. Fabrizio Ferrucci, Italy) made available to the author by the Ministry of Disaster Management and Human Rights of the Government of Sri Lanka. These LIDAR data have been acquired at a horizontal resolution of 1 m and a vertical resolution of not more than 0.3 m and were originally projected on to the UTM WGS84 -Zone 44N coordinate system.

Hydrodynamic M odel Formulation
The mathematical model used in the present work is the Cornell Multi-grid Coupled Tsunami Model (COMCOT, coded in FOR TRAN 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 has been validated by experimental data [7] and has been successfully used to investi gate several historical tsunami events, includ ing the 2004 Indian Ocean tsunami [6,10].
volume flux is assigned to be zero. Once the water surface elevation at the wet grid is higher than the land elevation in its adjacent dry grids, the "shoreline" is moved one grid toward the dry grid and the volume flux is no longer zero and needs to be calculated by the governing equations.
The amplitude of the 2004 Indian Ocean tsunami during its propagation was in the order of magnitude of 1 m, whilst the wavelength of the leading wave was in the order of magnitude of 100 km. Thus, linear shallow water equations are adequate to solve tsunami propagation in Regions 1, 2 and 3.
On the other hand, in the inundation areas in Region 4, the water depth becomes very small and approaches zero at the tip of the surging bore. Thus, the non-linearity, i.e., the wave amplitude to depth ratio, could become signifi cant. However, owing to the shallowness of the water, the frequency dispersion effects, which are represented by the water depth to wave length ratio, can still be negligible. Considering these, we adopt non-linear shallow water equa tions inclusive of bottom frictional terms to describe the tsunami overland flow in fourth level grids.
Using a nested grid system, COMCOT is capable of simultaneously calculating the tsunami propagation in the ocean and the inun dation in the targeted coastal zones. In the nested grid system, the inner (finer grid) regions adopt a smaller grid size and time step and are 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 is product of water depth and depth-averaged velocity, is interpolated from the outer (larger grids) region into inner (finer grids) region. The water surface eleva tions and the volume fluxes in the inner (finer grids) region are calculated and the resulting free-surface elevations are 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, we can capture nearshore features of tsunami dynamics with a high spatial and temporal resolution, and at the same time, we can still maintain a high com putational efficiency. To simulate onshore flooding, a moving boundary scheme described in [2] was employed, in which the "shoreline" is defined as the interface between a wet grid and its adjacent dry grids. Along the "shoreline", the As the spatial extent is comparatively larger in Regions 1 and 2, we use the linear shallow water equations on the spherical coordinate system as given in the following: where, Z , is free surface elevation; (ip,cp) denote the longitude and latitude of the Earth; Ris the Earth's radius; P and Qstand for the volume fluxes (P=huandQ=hv, with u and ubeingthe depth-averaged velocities in longitudinal and latitudinal directions, respectively); h is the still water depth; and/represents the Coriolis force coefficient.
However, given the comparatively smaller spatial extents of Regions 3 and 4, we solve the governing equations in Cartesian coordinate system in these regions. The non-linear shallow-water equations in Cartesian coordi nate system can be expressed as: where, Z , is free surface elevation; (x,y) denote the horizontal coordinates in the Cartesian coordinate system; xx and yt are the bottom shear stress in x-axis (pointing to the east) and y-axis (pointing to the north); P and Q stand for the volume fluxes (P = hu and Q = hv, with wand v being the depth-averaged velocities in horizontal and vertical directions); H is the total water depth. Moreover, because of the smaller spatial extents of the second, third and fourth ' level grids, Coriolis force effect is not consid ered for these regions.

Results and Discussion
At the outset, the numerical model was calibrated by comparing the model simulated extent of overland flow for the city of Matara with the field measurements of the points of maximum penetration of inundation reported in [11], It was found that a Manning's bottom roughness value of n = 0.02 (which is the typical value of n for firm soil) gave good agreement of model predicted extent of inundation with the measured points of limit of inundation, as shown in Fig. 8. Similar comparisons carried out for the cities of Galle and Hambantota also showed good agreement. Now, we consider in Fig. 9 the high hazard areas delineated based on the methodology outlined in Section 3 for the cities of (a) Galle, (b) Matara, and (c) Hambantota. The topogra phy of the areas of interest as well as rivers, streams and other waterways present are also shown in these figures. We see in Fig. 9(a)for the city of Galle that the high hazard area extends to a distance of about 300 m -500 m inland from the shoreline in that part of the city fronting the bay to the east of longitude ~80.22oE. We also see that the main bus terminal and the surrounding area including the locality to the north to a distance of about 750 m in the valley between the two hilly areas also falls into high hazard category based on the classification method used in the present paper. There is another notable stretch of onshore lands belonging to high hazard category to the w est of longitude ~80.20oE.
Further, the wall around the historic Dutch Fort ( Fig. 9(a)) has prevented the tsunami from directly attacking the area inside the Fort. We, however, see in Fig. 9(a) that certain areas of the Fort on its eastern side are classified as high hazard. This is because the flooding outside had made its way into the Fort through horse shoe shaped road entrances, and also as a result of overtopping of a comparatively lower part of the perimeter wall on the eastern side with sufficient momentum to be classified high hazard. The model results also indicate that Mahamodara Lake (~80.200oE, 6.043oN) and its waterways have conveyed the flow of water further inland thus making certain areas on the periphery of the lake hazardous.
In the city of Matara ( Fig. 9(b)), the high hazard area extends, in general, up to about 500 m -750 m onshore from the coastline, except the areas to the east of longitude ~80.56oE due to the elevated ground with steep beach fronts there and around longitude 80.52oE. It is also clear that River Nilwala has helped spread the high hazard areas further inland.
The numerical results indicate that a surge of tsunami induced flow rushes through the comparatively low-elevated ground on either side of location A (see Fig. 9(c)) towards Saltern-1 thus making this area extremely hazardous. We also see that the sand dune stretching from around (81.133oE, 6.134oN) to (81.148oE, 6.140oN) has prevented the tsunami from directly attacking the settle ments behind thus making these areas safer. Longitude (in degrees)

ENGINEER
Note also that the coastal lands to the south of latitude ~6.127oN (except those adjacent to the shoreline) are largely safer from the tsunami threat owing to the hilly terrain with steep onshore slopes. However, a breach of sand dune at location Bhas facilitated the tsunami overland flow to rush through with consider able momentum.
Finally, it must be added that there are several limitations involved in modelling tsunami overland flow. One is that, since the initial condition for the modelling is determined by the dislocation of the ocean bottom along the fault line, the largest source of errors is the earthquake model. Another significant limita tion is that the resolution of the modelling is no greater or more accurate than the bathymetric and topographic data used. Although, high resolution topographic data were available for the present study, the resolution of the bathy metric data used is much less.Nevertheless, the results from the present study are useful for local authorities and emergency managers to identify areas that should necessarily be evacu ated in the event of a major tsunamigenic earth quake in the Andaman-Sumatra fault line as well as in public education and awareness activities. Because of the uncertainties inherent in this type of modelling, these results are, however, not intended for land-use regulation.

Conclusions
A numerical coding of the non-linear shallow water equations has been employed to compute the temporal variation of flow depths and flow velocities across areas of interest in three cities on the south coast of Sri Lanka. The flow parameters obtained in this way have been combined with the results of an experi mental study on the human stability in tsunami overland flow to demarcate high hazard areas that should be evacuated in the event of a tsunami warning due to an earthquake similar to that occurred in December 2004.