Natural Hazards and Earth System Sciences Examination of three practical run-up models for assessing tsunami impact on highly populated areas

This paper describes the examination of three practical tsunami run-up models that can be used to assess the tsunami impact on human beings in densely populated areas. The first of the examined models applies a uniform bottom roughness coefficient throughout the study area. The second uses a very detailed topographic data set that includes the building height information integrated on a Digital Elevation Model (DEM); and the third model utilizes different bottom roughness coefficients, depending on the type of land use and on the percentage of building occupancy on each grid cell. These models were compared with each other by taking the one with the most detailed topographic data (which is the second) as reference. The analysis was performed with the aim of identifying how specific features of high resolution topographic data can influence the tsunami run-up characteristics. Further, we promote a method to be used when very detailed topographic data is unavailable and discuss the related limitations. To this purpose we demonstrate that the effect of buildings on the tsunami flow can be well modeled by using an equivalent roughness coefficient if the topographic data has no information of building height. The results from the models have been utilized to quantify the tsunami impact by using the tsunami casualty algorithm. The models have been applied in Padang city, Indonesia, which is one of the areas with the highest potential of tsunami risk in the world.


Introduction
Estimating tsunami behavior on land is essential for mitigating tsunami damage and for designing optimum community evacuation plans along the coast (Imamura, 2009).The main goal of community evacuation plans is to complete the Correspondence to: A. Muhari (aam@tsunami2.civil.tohoku.ac.jp) evacuation before a tsunami comes (Koshimura et al., 2006).However, learning from the Japan experiences, even in cases when tsunami warning was well designed and quantitative information about tsunami arrival time and wave height was available, the activities necessary to mitigate damages were still difficult to be well implemented (Imamura, 2005).Even when tsunami warning and evacuation information were issued just after the earthquake, several problems before and during the evacuation were encountered by some communities that probably had inadequate time to complete the evacuation.One of the lessons learned was that in order to minimize the tsunami casualty toll, it is necessary to exclude from the potential evacuation route those roads where tsunami casualty is likely to happen.By using a detailed tsunami run-up model coupled with a tsunami casualty model, this problem can be solved.This is the main rational for the study presented in this paper.
The biggest obstacle that usually arises when performing detailed inundation calculations is the availability of coastal topographic and bathymetric data.Optimizing general open access data such as GEBCO (British Oceanographic Data Center, 2008) or other nautical charts in a particular area by performing several interpolations often does not reflect physical reality with the desired accuracy.On the other hand, making significant efforts to create a very high resolution topographic data set that also includes the building height information is costly and, furthermore, requires very fine computational grids that imply longer computational time.Such problems are mostly experienced by developing countries where neither high resolution spatial databases are yet available nor are very good performance computing facilities.To overcome these problems, the present study intends to provide a method to assess tsunami impact even if very detailed topographic data do not exist.
This paper is divided into two parts; the first part examines three practical run-up models differing in the way they handle topographic data to determine the one giving the most Published by Copernicus Publications on behalf of the European Geosciences Union.
A. Muhari et al.: Examination of three practical run-up models for assessing tsunami impact reliable results.The second part integrates the chosen runup model with an extended tsunami casualty model based on Koshimura et al. (2006).
Considering that the results of tsunami inundation analysis are used for further evacuation planning, annual updating is necessary to keep the analysis in line with the city development.In this context, a simple but sufficiently highaccuracy numerical model is needed.Also, the model should be able to at least describe the tsunami inundation characteristics on roads that will be used for evacuation.For this reason, in the first part we look at three practical run-up models originally based on the TUNAMI code (Imamura, 1996), which has already been used world-wide for tsunami calculation.The first model is called Constant Roughness Model (CRM).It applies a uniform bottom roughness coefficient throughout the model area.The second is named Topographic Model (TM) and includes the building height integrated with the topographic data.The third model is the Equivalent Roughness Model (ERM), which uses different bottom roughness coefficient that are function of land use and of the percentage of building occupancy on each grid cell.
When the tsunami inundates densely populated areas, Goto and Shuto (1983) stated that there are three important regions where major hydraulic phenomena occur during the flow passing through obstacles: The first is the entry region, including the first row of obstacles; the second is the rearmost region; and the third is the area between the first and the second.In the first region, the head loss is mainly caused by sudden contraction, while in the second region the sudden expansion is the most important occurrence.In the third region, the tsunami flows like through a narrow channel of roughened walls.However, based on their findings, the effect of the third region is relatively small compared with the first region.Therefore, in the first part of the paper we omit considering the third region in the examination of the three run-up models.Moreover, the present study uses two-dimensional models that do not allow the calculation of roughness on the vertical wall.
In the second part, the results of the method we selected out of the three models examined are utilized to calculate the hydrodynamic force acting on a human body.Once the human body cannot stand against the tsunami flow, it is assumed that a "tsunami casualty" occurs.The potential of tsunami casualty is then spatially quantified by the Tsunami Casualty Index (Koshimura et al., 2006) to establish the new term called the Road Risk Map.
To illustrate the model applications, Padang city in West Sumatra Province, Indonesia, was chosen as a case study.This city is exposed to tsunami threat due to the large amount of accumulated energy since the past events.Sieh et al. (2008) predicted that probably the locked area will be ruptured within the next several decades.To assess the potential impact of tsunami from this predicted earthquake, several studies related to tsunami inundation modeling were conducted.Borrero et al. (2006) produced a preliminary tsunami source scenario and inundation model with relatively coarse data (200 m).Goseberg and Schlurmann (2008) and Taubenböck et al. (2009), through the "Last-mile" project, provided very detailed surveyed topography (0.25 m resolution) and bathymetry (3 m pixel spacing) data sets that allowed them to image the buildings integrated along with the Digital Elevation Model (DEM).In principle, depending on the stability of the numerical method, such detailed data should give better results.But still, from the practical point of view, the use of such detailed data in numerical calculations is highly inefficient because it requires very long computational time.Gayer et al. (2010) used different data and methodology to optimize the bottom friction.In their work, the resistance factor of the buildings in residential areas was taken into account by assigning different values of roughness, depending on the land use conditions.They performed several numerical calculations using an "artificial city" to obtain an adequate Manning roughness coefficient in an urban area.The resulting various values of Manning coefficients were then set as time-independent variables along the simulation.
In view of future implementation and of the data and equipment required for further updating, it is difficult to suggest which model is more suitable for a given community because the previous efforts mentioned above were conducted under different input data accuracy.For this reason, the present paper first examines the performance of three run-up models at the same grid size, and then it exercises the method we propose on a larger grid size to demonstrate the capability and the accuracy of the model.The topographic data we use here represents the highest possible accuracy that probably can be obtained from other areas of Indonesia where the proposed method will be applied within the framework of the nation-wide tsunami inundation mapping project.

Tsunami source
The city of Padang is located in the center of the west coast of Sumatra.The population of the city is approximately 850 000 people.In the last decade, the subduction zone off the west coasts of Sumatra Island was identified as the most active tectonic region.As shown in Fig. 1, the northern part of this region ruptured in 2004 and 2005.The southern part was partly ruptured in 2007 (Konca et al., 2008), and the southernmost is considered ruptured in 2000 (e.g.Abercrombie et al., 2003).The above events left a gap centralized in the Mentawai Islands, which has stored the energy since 1797 and overlapped with 1833 event (Natawidjaja et al., 2006).This area is strongly coupled (Chlieh et al., 2008), and it is expected to cause a M w = 8.8 (Sieh et al., 2008) or a M w = 8.9 (Natawidjaja et al., 2009) earthquake in the near future.To start the numerical analysis, we use the tsunami source based on the work of Natawidjaja et al. (2009), implying the most reliable tsunami source scenario in Padang city (e.g.Muhari et al., 2010).The predicted earthquake with M w = 8.9 creates a maximum initial sea surface height of 3.7 m calculated using the seismic deformation model by Okada (1985).Previous results by Muhari et al. (2010) show that this source generates a maximum tsunami height of 6 m in the coastline, and approximately 20 min of arrival times is predicted along the coasts.In the present paper, the sea floor displacement is assumed to be equal with the sea surface deformation, and there is no effect of rise time nor is the rupture velocity included in the tsunami generation model.Some of the accumulated energy in the southernmost area was probably released after the M = 7.7 Mentawai earthquake in October 2010 (USGS, 2010).However, since the highest energy beneath the Siberut Island still remains, the scenario adopted in the present study is still considered reliable.

Tsunami propagation
We prepared a set of bathymetric and topographic data to construct a nested grid system consisting of five sub-regions.The grid size varies from 405 m for regional tsunami propagation to 5 m for local inundation in all run-up models.
The ratio between the larger to smaller region is set as 1:3 to satisfy the numerical stability (Fig. 2).We resampled the GEBCO 08 data set (British Oceanographic Data Center, 2008) from 30 arc-sec (about 926 m) down to 405 m for offshore tsunami propagation in region 1.Then the SRTM data (CGIAR-CSI), surveyed bathymetric data (MoMAF, 2006) with 200 m accuracy, and local bathymetric charts number 165 and 279 were compiled and interpolated to obtain the 135 m and 45 m grids for region 2 and region 3, respectively.Lastly, we used the Digital Terrain Model (DTM, Intermap, 2008a) to create 15 m and 5 m grids for region 4 and 5.All of the necessary extra-interpolations especially for bathymetric data in the fourth and fifth domains were conducted using a Kriging exponential method available in ArcGIS software (e.g.Naoum and Tsanis, 2004).
The DTM data in the present study has its original accuracy of 5 m.In a rectangular grid system, it produces the area of 25 m 2 for each grid cell.Concerning the applicability of this cell size to represent the building features in the study area, we analyze the size of around 7500 building in the central part of Padang city.From the GIS analyses, we found that only around 730 buildings have an area below 25 m 2 (with an average of 18 m 2 ), i.e. less than 5 % of the total buildings stock (with average of 174 m 2 ).From this finding, we assumed that the used grid size is able to represent the building shape in the study area, even though we acknowledge that several narrow roads with less than 5 m width are misrepresented and covered by buildings.Another concern of data accuracy is the position of the buildings that are not fixed in the rectangular grid system, with the consequence that the building shape is automatically adjusted according to the grid arrangement around them.The data processing for each model will be discussed furthermore in Sect.2.3.
Tsunami threat in Padang city is characterized by a nearfield tsunami case with unique bathymetric features (Fig. 2).Therefore, we use the non-linear shallow water equations to describe the tsunami propagation along continental shelf until it inundates inland.The non-linear shallow water equations are discretized by using the Finite Difference Method (FDM) with the leap-frog scheme (Imamura, 1996).The set of equations is shown below, In these equations, M and N denote the discharge flux in the x and y directions, respectively; η denotes the water elevation, n stands for the Manning coefficient, and h represents the water depth.

Constant Roughness Model (CRM)
The CRM is the conventional type of tsunami inundation model.It considers only ground elevation influences the inundation characteristic.The basic assumption used in this model is that buildings are not able to withstand the tsunami, and the difference of the bottom friction coefficient without building is very small to be able to influence the tsunami flow.
In this model, we used the DTM for topographic data.The bottom friction coefficient is set as constant along the computational grid with the value of 0.025.This value is adequate for natural channels in good condition (Linsley and Franzini, 1979).The topographic data used in this model is shown in Fig. 3.No change in Eqs.(1) to ( 4) is needed for this run-up model.

Topographic Model (TM)
In the Topographic Model (TM), the building shape and height is integrated in the Digital Elevation Model.It means that the buildings are assumed to withstand against the earthquake and remain standing during the tsunami.In this framework, the use of the shallow water theory to simulate the runup with such detailed grid was shown and verified by Goto We first extracted the building shape from the highresolution satellite images (Fig. 4a and b).Then the GIS file of the building shape is used to extract only the building height information from DSM data (Fig. 4c).Next, the building height data from (c) is overlaid to DTM (Fig. 4d).At this stage the "new" DTM has the building height information as is shown in Fig. 4e.It is remarked that, when a grid cell consists of both building and land, the feature in the center of the cell is automatically taken as the cell attribute (land or building).It implies that the area of some buildings is increased or decreased by less than 25 m 2 for each cell unit as is mentions in the previous sub-chapter.
In the numerical model, the existence of these "physical buildings" blocks the tsunami and changes its direction if the flow depth is less than the building height (Fig. 5a).As a consequence, the resulting inundation area has dry cells in the area covered by buildings that are not submerged by the tsunami.If the tsunami flow depth is higher than the building height, it overtops and submerges the buildings in the following time step.The calculation when the flow overtops the structures in the model is made by means the Homma formula as is described in Imamura (1996).
One should note that even the finest grid of 5 m does not allow us to simulate micro-phenomena such as the free surface flow when tsunami passes through the buildings, including the water splash effect.The present grid size is also not sufficient to determine the turbulence or the vortex shedding surrounding the structures since most of the width of the narrow roads is represented by only one or two cell(s).However, simulating such micro-phenomena is not the aims of this paper.Our interest in the examination is substantially focused on assessing the change of flow depth and velocity pattern due to the effect of large obstacles in the city scale.

Equivalent Roughness Model (ERM)
The work of ERM was formerly initiated by Goto and Shuto (1983), and it was continued by Aburaya and Imamura (2002).The basic assumption is similar to TM that the entire buildings are able to withstand against earthquake and tsunami.The main idea is to represent building resistance without exactly having the building features integrated in the Digital Elevation Model.The resistance is then represented by an equivalent roughness.The conception of ERM is shown in Fig. 5b.
The formulation of the methodology differs depending on building size relative to the grid size.It is explained in the following subsection.

The case of the building width (k) smaller than the grid size (dx)
Consider a control volume V formed by rectangular grid dx = dy and a tsunami flow depth D as shown in Fig. 6.The force acting on the flow (G) in the control volumeV is equal to bottom shear force (R 1 )+ resistance of obstacle (R 2 ) As previously mentioned, the ERM does not have the building height information on its topographic data, so the resistance of buildings and other man-made structures (R 2 ) should be taken into account through an equivalent roughness.
Assuming that A 1 is the bottom area occupied by water and that A 2 is the projection area of the obstacle underwater, then the above parameters R 1 and R 2 can be obtained by the following equations, Where dx and dy are the grid size in the x and y directions, respectively, and θ is the percentage of the bottom area occupied by the building in a grid on the numerical domain.  .The projection areas of the obstacles under water in a grid in the condition of D > H , and D < H consequently are given as follows: However, Eq. ( 8a) requires the building height information that it is not owned by the ERM.Therefore, it was simplified by assuming that the height of the projection area is equal to the flow depth.Next, we calculate the force exerted on water in the control volume V in Eq. ( 5), starting from the expression based on the friction slope concept (I ), which is G = pgV I (Goto and Shuto, 1983).Introducing the equivalent roughness n e in the new form of Eq. ( 5), now we have the relationship as follows: Considering that n 0 is the bottom roughness without obstacles, the friction force in the area occupied by water can be written as: The resistance due to the obstacle is assumed to depend on the square of velocity, and thus the drag equation can be represented as follows: where C D is the drag coefficient.FEMA 55 (FEMA, 2003) suggests drag coefficient values from 1.25 to 2, depending on the ratio between the width of building and the flow depth.In case of highly populated areas where the building size varies from small to large, it will be very difficult to assign a different value for each building.Therefore, in this case we choose a fixed average value of 1.5 applied in the residential area.G, R 1 and R 2 are forces and have same unit of kg m s −2 .Substituting Eqs. ( 10) and ( 11) into Eq.( 9), the equivalent roughness n e is obtained as follows: Substitution of the Eq. ( 12) into the Eqs.( 2) and (3) yields the following equations: The methodology was tested by using the experimental data by Aburaya and Imamura (2002).The hydraulic experiment was carried out using a 12 m wave tank with two obstacles assumed as buildings (Fig. 7).The comparison between experimental and model data yields the conclusion that finer grid in ERM yields better agreement with the experimental data.
It is worth recalling that ERM in the version presented in this subsection was applied by Suppasri et al. (2011) who used a 52 m grid size in the region with the highest resolution, and by Koshimura et al. (2009) who used a grid size of 23 m to model the 2004 tsunami impact in Khaolak and Patong Beach in Thailand, and in Banda Aceh, Indonesia.They obtained a good agreement with the field data survey.However, it is to be stressed that it is impossible to model the tsunami flow along the roads in residential areas with a grid size as big as above, because it is bigger than the width of the road and the average building size.This is the reason why we propose the use of the ERM by using an even smaller grid size with several enhancements as described below.

The case of the building's width (k) larger than grid size (dx)
Application of ERM using a finer grid (5 m) implies that some grid cells are fully occupied by a building.To illustrate the application of ERM in this case, consider a control volume that consists of several grid cells with building A and B hit by a tsunami with velocity u as shown by Fig. 8.
In the numerical calculation, the tsunami flow passes across buildings.The building resistance (R 2 ) is given by Eq. ( 11).The value of R 2 may differ from cell to cell, depending on the flow depth on each cell in the respective time step.For example, the values of R 2 for the cells in the second row of a building may be lower than the ones in the first row because the flow depth already decreased when it passed the first row.Also, the values of R 2 may vary spatially, depending on the percentage of building occupancy in each cell.For the whole area of a building, the total resistance undergone by the flow in a specific time step can be written as R = N f × R 2 , where N f is the number of grid cells opposing the wave front.
High percentage of occupancy in Eq. ( 12) causes the equivalent roughness value to significantly increase, and reduces the flow speed when a tsunami passes through the building areas.In case of long wave run-up, deceleration in the wave front caused by bottom friction will cause the amplification of the flow depth.If it occurs over a relatively long time, it will change the direction of the flow to the surrounding alleys.However, there will be no dry cells in ERM.The flow will keep inundating the occupied grid,  and the amplified flow depth will be distributed over the surrounding grid cells.Therefore, all the cells in the inundated area in ERM will be submerged.Consequently, no change is needed on the mass continuity.For numerical stability, it is set that the maximum allowed value for the percentage of building occupancy is 99, so Eqs.( 6), ( 7) and (12) can still be utilized.
The ERM requires more input data than the previous models.In addition to topographic data from DTM (Fig. 9e, bottom figure), ERM requires the bottom roughness coefficient obtained from the land use information, and also the percentage of building occupancy data that is processed by using GIS.The overview of the data preparation is explained as follows with reference to Fig. 9.
Land use data derived from the high-resolution satellite imagery (Fig. 9a) by German Aerospace Center (DLR, 2008) was reclassified according to the values given in Table 1.The roughness data is then gridded with the same cell size of topographic data (Fig. 9d), and saved as separated file input (Fig. 9e, center figure).
As to the percentage of building occupancy data, the extracted building shape (Fig. 9b) is overlaid to the grid file, and the occupancy ratio is calculated by means of GIS tools (Fig. 9c).The resulting output is then saved as an independent input file (Fig. 9e, top figure).
To explore the performance and accuracy of the ERM in a larger grid size, using the same data sources, we used a grid size of 15 m as the smallest region in the nested grid system.This size is a little bit larger than the average building size in study area and two consequences arise.The first is that the domain model only consists of four regions.The second consequence is the increasing number of small streets that are not described as a road because they are smaller than the size of the grid, and that can be interpreted as a building area as result of the existence of houses along them.
One should note that the lack of "physical building" and the selected grid size in the present model implied that the approach may not be perfectly able to reproduce the drag due to building.However, as it was mentioned before, the purpose of our examination is to promote a practical method to use when very high accuracy topographic data including the building height is not available.Indeed ERM is devised for cases when the building height information does not exist.
Ideally, the examination should be performed referring to the field observation data from an actual tsunami as a benchmark.Unfortunately, there is no field measurement data from the past tsunamis available in the study area.Therefore, we use TM as the benchmark since it was run with the most complete topographic data including building height information.
The models were run for two hours with 0.05 s time step.The model execution on a quad-core workstation with 20 GB memory requires the maximum computational time of four days if three of them are running at the same time.Notice that for the ERM with 15 m grid size, the computing time required is only less than one day.

Tsunami impact assessment
In this section, we integrate a module of an extended tsunami casualty model into the proposed run-up model.Based on casualty data from 2004 Indian Ocean tsunami in Banda Aceh, Koshimura et al. (2009) shows that the ratio of death is correlated with the inundation depth of tsunami flow.This is in line with the observations of Nishikiori et al. (2006a), whose research in Ampara District, Sri Lanka showed that more than 75 % of tsunami victims died due to drowning.And this is also confirmed by the autopsy result of 13 135 bodies of tsunami victims of Japan 2011, where 12 143 of them (92.4 %) died due to drowning (NPA, 2011).
On the other hand, Nishikiori et al. (2006b) and Doocy et al. (2009) showed that the tsunami casualties in Sri Lanka and Banda Aceh from the 2004 tsunami vary depending on gender and age.Based on the above facts, if we can quantify the possibility of tsunami impact on various age and gender based on hydraulic properties of tsunami inundation, safer evacuation route can be determined.
In the past, numerous studies were conducted to obtain empirical relationship of human-flow interaction based on hydraulic experiment (e.g.Karvonen et al., 2000;Jonkman et al., 2008), and on numerical assessment (e.g.Endoh and Takahashi, 1995;Lind et al., 2004;Koshimura et al., 2006, andYeh, 2010).Yeh (2010) initiated the use of anthropometric data in a tsunami casualty model.However, a better representative of the human body factor is probably needed in order to utilize the anthropometric data in more detail.This is the main point of the present model.Oversimplifications such as representing the human body as a simple single cylinder or as two elliptic cylinders will not perfectly describe the different sizes of each part of the human body that yield to different projection area against the hydrodynamic force.
This study enhances the human casualty model developed by Koshimura et al. (2006) through the completion of the human fall mechanism based on Endoh and Takahashi's paper (1995).It is supplemented by a preliminary verification using experimental data by Endoh and Takahashi (1995).The actual human body was represented by a set of cylinders with the composition as shown in Fig. 10.
Two types of fall are considered: the first is if the hydrodynamic force exceeds the friction force on the human feet soles that yields to the slipping fall type; the second is if the moment in the bottom back of heel due to the hydrodynamic force is bigger than the resistance moment due to human body weight, which leads to the tumbling fall type (Endoh and Takahashi, 1995).We described the hydrodynamic force incorporating the tsunami flow through the Morrison equation (Morison et al., 1950).The governing equations for the slipping fall type and the tumbling fall type are expressed below: 16) In these equations, µ is the friction coefficient in the evacuee's feet soles, W 0 is the body weight, W is the buoyancy, h G is the vertical distance from the floor where the resultant force acts as the function of flow depth, and I G is the distance between the center gravity of the body and the bottom back of the heel.In the Morrison equation, ρ s denotes the sea water density, u is the velocity, S is the projection area according to flow depth, ∂u ∂t is the local acceleration in the respective point, V is the volume of submerged human body.C D and C M are the drag and inertia coefficients.
In the equation above, C D is determined as 1.1(1-Lf/D), where Lf is the distance between feet, and D is the tsunami flow depth.I G is determined as the 1/4 of the feet (Endoh and Takahashi. 1995).Lf and I G are obtained from the anthropometric data.C M is assumed as 0.5 (Koshimura et al., 2006).The friction coefficient in the evacuee feet soles is assumed as 0.7 (Yeh, 2010).
Equations ( 15) and ( 16) carry the implication that in a specific human body size, slipping fall type will be mostly determined by the friction coefficient on the human feet soles.The slipping fall type usually occurs in the areas where the flow depth is not so high, but the velocity is large.On the other hand, the distance from gravity center to the bottom back of the heel becomes a critical component on the tumbling fall type, a condition that usually occurs when the flow depth is high but low velocity.Figure 11 shows the influence of these aspects to determine the critical condition of human instability.The graph was created numerically using the human body criteria as it gives in Table 2 for adult male.
The size of each cylinder obtained from the anthropometric data varies for each country.For Padang city, we use anthropometric data for for adult male and female Indonesians  from Chuan et al. (2010).The size of each cylinder according to the gender classification is given in Table 2.The calculation of hydrodynamic forces acting on a projection area depends on the tsunami flow depth.Therefore, the computation is divided into three sections (refer to Fig. 10).First, when the flow depth is less than H 1 (arm height); second, when the flow depth is between H 1 and H 2 (leg height); and the third when the flow depth is less than H 3 (shoulder height).The condition where the flow depth is higher than H 3 is not considered because the human body will tend to float rather than fall.
The quantification of potential tsunami casualty in spatial distribution uses the so-called Tsunami Casualty Index (TCI) introduced by Koshimura et al. (2006) as follows: T C is the time interval of run-up condition that satisfies Eqs. ( 15) or ( 16), and T I is the total time of tsunami runup during the simulation (the duration of the interval during which the cell is submerged by water).TCI is an index to show the tsunami casualty potential spatially during the calculation.So, it can be extracted at the end of the simulation or until a specific time of tsunami occurrence, i.e. until the expected time for the completion of tsunami evacuation operations.
Before applying this model to the study area, we conducted verification using experimental data from Endoh and Takahashi (1995).The experimental data was obtained from the stability test of somebody who is 164 cm tall and 64 kg weight in an open channel as shown in Fig. 12.In this experiment, the person standing at point B faces the direction of the flow.We arrange a set of cylinders of the same size in the model and run it on similar condition of flow depth and velocity as the experiment.The verification result is shown in Fig. 12.

General results
Here we discuss the spatial features of run-up heights and flow velocities obtained by three run-up models.Maximum run-up height distribution of three models is given by Fig. 13.In terms of the maximum inundation line, significant discrepancy appears when we compare the results of TM and ERM with CRM.This indicates that the influence of flow resistance caused by the existence of "physical building" in TM, and by the equivalent roughness coefficient in ERM, have a significant impact for the determination of the tsunami inundation area.These results show the limitations of CRM: since the bottom friction is set constant along the simulation, hence the characteristics of the tsunami flow are only influenced by the ground elevation.The CRM produces an inundation area of about 21.46 km 2 , with an average flow depth of 2.77 m.On the other hand, the TM estimates that 11.36 km 2 will get inundated with an average 2.55 m deep flow, while ERM obtains 12.71 km 2 and 2.29 m average flow depth.
We observe that the TM has a maximum inundation line that is greater than the ERM, but the flooded area of the TM is smaller than the ERM.This is due to the tsunami in the TM inundating only non-building area plus the area of buildings with a height less than the tsunami flow depth.
For practical reasons, the ERM can be used on a largercell grid since the value of bottom roughness is influenced by the percentage of building occupancy in a grid.However, one should note that a lower resolution grid may not be able to represent in detail the shape of topographic features in the study area.It can be seen from the result of the ERM with 15 m cell size that there is a general underestimation of flow depth compared with the one from the finer grid, especially in the area near to the bend of the river flow.Nevertheless, differences in flow depth in these areas are not more than 0.5 m.
Meanwhile, in another distant part of the river, it appears that the difference in the size of the cells is not very influential in determining the maximum inundation line.In this region, the maximum inundation line of the ERM with 15 m cell size is very identical to that obtained by using the ERM with 5 m cell size.So, here we return to common understanding that finer grid size will yield better representation of topographic features.

Run-up height and velocity distribution
Here, we will discuss the abilities and the limitations of three run-up models in more detail.Along the coastline of Padang, the government built sea walls and groins to protect the shore from erosion, the groins were constructed both parallel and perpendicular to the coastline (Fig. 14, upper photos).The existence of these structures affects the tsunami propagation from offshore to onshore.In this area, after passing through the sea walls, the maximum flow depth drops drastically.In order to show the impact of these structures in influencing the flow depth distribution in the study area, we extract the maximum flow depth from the control points in Fig. 13 (TM) that are marked with a sequence of consecutive numbers from one to six from top to bottom.All three models show similar results as shown by Fig. 14 (bottom figure).These facts imply that our desire to compare the differences of the inundation behavior of the three models when the flow enters the first row (waterfront) of the residential-area cells cannot be fulfilled.Differences will appear later that are more onshore, especially in correspondence of intersections and crossroads that separate the previous and the next blocks in the residential area.
The maximum run-up heights and velocity profiles from the cross-sections one and two in Fig. 13 (ERM figure) are extracted and shown in Fig. 15.These cross-sections were chosen because they represent two different characteristics.Cross-section one has a higher topography landward.It has two intersections at a distance of 380 m and 780 m from the coastline (see bottom panels -left hand side of Fig. 15).In cross-section two, the direction is not perpendicular to the coastline.The direction forms an angle of about 40 degrees to the coastline from the north.It also has two intersections at a distance of 150 and 350 m from the coastline in the direction of the road.The ground is relatively flat with a height difference of less than 1 m.In addition, the crosssection is characterized by a reverse slope to within 100 m from the coastline (see the bottom panel -right hand side of Fig. 15).

Run-up height distribution
In cross-section one, TM shows that the run-up height is relatively constant in the region from the coastline until the first intersection.In this area, the run-up height of TM is higher than the ones from the ERM and of the CRM.Next in the first intersection, the sudden expansion phenomenon is clearly visible and influential on the sudden decrease in run-up height profile of the TM.As said by Goto and Shuto (1983), in the narrow road, the water flows as if threw a narrow channel: it means that buildings that are aligned along the narrow road in the TM create the channel effect that perfectly blocks and amplifies the tsunamis until it meets the first intersection.At the intersection, the amplified tsunami is spread out through the roads in different directions, causing a sudden drop of the run-up height.
The same characteristic of the TM is visible between the first and second intersection.Along the narrow road, the run-up height is relatively constant until it then drops significantly at the second intersection.In this region, the run-up height of TM is lower compared to the CRM.This is because the amplified run-up height in the first region of the TM is already scattered in other direction after it passed the first intersection, whereas the CRM does not experience the sudden expansion phenomena because it does not take into account the effect of buildings on its calculation.
Different characteristics are shown by the two ERM, with 15 m and 5 m cell size.As mentioned earlier, this cross section has a higher topography landward.This leads to a decrease in flow velocity at the wave front.Supposedly, in the long wave phenomenon, the decrease of flow velocity at the wave front results in increased run-up height because the velocity of the second wave, the third and so forth still remains.However, there is no vertical wall in the ERM that would make a perfect channel effect.The resistance of the buildings is represented by the high value of bottom friction, depending on the percentage of building occupancy in each cell.But it does not prevent the flow to remain flooding the building area.The amplified run-up height at the wave front will then be distributed over the surrounding cells, whether they have the same or lower bottom friction values.Therefore, the amplification of tsunami run-up when running down the narrow roads makes the elevation of the TM to be greater than the ERM which allows it to flow further.This is the cause that makes the maximum value in the ERM to be not as big as in the TM, and becomes one of the limitations that should be noted when using the equivalent bottom roughness to represent the resistance when tsunami flows in residential area.It also causes the maximum distance of tsunami inundation produced by the ERM shorter than in the TM in some areas.
Another phenomenon of TM that is not observed in the ERM and in the CRM is the difference between the run-up height in front of and behind the building (indicated by significant difference of run-up height at several points along the cross section in Fig. 15, left hand side-top panel).But one should keep in mind that the extracted data is the maximum value during the simulation.Perhaps the ERM has the similar phenomena as the TM when the wave front passes through a building in the first attack.However, the accumulation of the run-up height when all the cells are already inundated during the simulation make these values not to be stored as the maximum value during the computational time.
In cross-section number two, better results were shown by both of the ERM.Although the run-up that results is not as big as the one produced by the TM, flow characteristics that indicate the sudden expansion phenomenon is successfully demonstrated by the ERM with a pattern similar to the TM.This indicates that the application of ERM on representing the building resistance would be better in a gentle topography.

Flow velocity distribution
This section explores the characteristic and the comparison of the flow velocity of the TM and the ERM.We purposely omit the CRM in this section because in the previous analysis it was seen that CRM does not exhibit specific behaviors when the tsunami flows in residential areas.
According to Goto and Shuto (1983), the tsunami flow depth and its velocity will change when the tsunami flows at the entrance depending on the coefficient of contraction that is a function of the ratio of contraction and the Froude number.Based on the equation that they recommend on their research, it can be predicted that the tsunami flow velocity at the entrance cells will be larger than the previous cells.
The results of the TM in both cross-sections show a consistent pattern, that is an enlarged tsunami flow velocity when starting to penetrate the first row of the residential area (Fig. 15).However, the results of the TM are not perfectly replicated by the ERMs.The increasing tsunami flow velocity of the ERM at the observation points as shown in crosssection one (Fig. 15, left) is not as significant as given by the TM.Especially when compared to the ERM with 15 m cell size, the resulting velocity profile is not able to represent the local effects since cells cannot fully describe the geometry of the cross section.In the larger cells grid it occurs that some part of the roads are not considered to be a road, but instead a building with a high value of building occupancy.with 15 m cell size are smaller, the pattern generally can represent an increase in flow velocity, especially at the second intersection.
Expanding the area of analysis into the city scale (Fig. 16), referring to the TM, satisfactory results are obtained by the ERMs compared with the CRM.The higher velocity distribution along the roads compared with the ones in the residential area is clearly visible in the ERM.Both ERM models which have 5 and 15 m cell size showed the same results.It indicates that basically the use of the equivalent roughness in the ERM could make a significant difference between the tsunami velocities in the settlement blocks (the grid that is occupied by building) and the tsunami velocity on the roads in these settlements.This is one of the expected results from the application of the ERM in a rigorous grid.
It is realized that modeling the tsunami run-up with fine grid yields a further challenge to consider the building damage due to hydrodynamic force during the tsunami inundations in the residential areas.Based on the building damage analysis of the 2004 tsunami in Banda Aceh, Koshimura et al. (2009) found that the buildings are destroyed with more than 50 % probability if the flow depth is deeper than 3 m or flow velocity is faster than 2 m s −1 .The maximum flow depth and maximum flow velocities produced by the CRM are respectively 7.88 m and 12.8 m s −1 .The TM produced 7.83 m and 16.5 m s −1 , while the ERM with 5 m cell is 7.82 m and 12.2 m s −1 , and the ERM with 15 m cell yielded 7.78 m and 10.1 m s −1 .In this context, for evacuation planning purposes especially to avoid the underestimation on the determination of the evacuation zone for future events, the results of the CRM should be used as a limit of the maximum extent of tsunami inundation area.Experience from the 2004 tsunami in Banda Aceh, and the 2011 tsunami in Japan showed that most of the buildings may not be able to withstand the tsunami.Although some of them can survive the surge force in the wave front, but very rarely the building can withstand the drag accumulation in the zone where the flow velocity is more than 2 m s −1 .However, on the other hand, the results given by the TM and ERM will be very useful to see more detailed flow characteristic in residential areas, and the assessment of the potential accumulation of hydrodynamic forces associated with the building and human vulnerability in tsunami affected area.

Tsunami impact to humans
In accordance with the purpose of the study mentioned at the beginning of this paper, we would like to propose a practical method to estimate tsunami impacts on human beings, particularly in regards to the preparation of evacuation plans.So, the tsunami casualty model was applied to the TM and to the ERM.We also analyze the performance of the ERM on representing the TCI on the roads that may be used to evacuate from the tsunami inundation zone.The tsunami casualty algorithm that was applied in the ERM with 15 cell size used an average adult male body size.Next, we utilize the spatial distribution of the TCI to generate the road risk map that shows the potential occurrence of tsunami casualty along the computational time, or during the first two hours of tsunami attack in reality.
The first impression that can be seen from the results shown by Fig. 17 is that not all regions in the tsunamiaffected areas are threaten in the same way.Some places, especially close to the inundation limit, may still be flooded by a tsunami up to 0.5 m.However, very low tsunami velocity makes the hydrodynamic force that accompanies tsunami flow unable to affect the balance of adult humans.Note that we extract the TCI distribution of the ERM only in an area that is inundated by the tsunami on TM for color clarity reasons.
The average TCI value obtained from the TM for adult males and adult females were respectively 0.59 and 0.61.In the same order, the ERM resulted 0.48 and 0.52.The ERM with 15 cell size produced the average TCI for adult male of 0.47.Even though the ERM produces a lower result compare to the TM, the differences are not more than 0.12 (∼10 %) for each of classification.So it can be said that the ERM can provide the predictive result for the TCI with the margin error not more than 10 %.

Conclusions
We have conducted the examination of three practical tsunami run-up models on highly populated areas.Significant difference between the TM and the ERM to CRM confirms the important aspect of the building existence when tsunami flows in residential areas.Our examination results show that if the building height information is not integrated to the digital elevation model, one can use the equivalent bottom roughness depending on the percentage of building occupancy.Here we must state that the selection of the grid size for the ERM should be based on to what purpose the analysis is carried out.If the analysis is intended to look at the general characteristic of tsunami inundation in the area with a topography which is relatively less varied, then one can use a cell size that is larger than the building size.On the contrary, if the ERM is used to look in more detail into the characteristic of tsunami inundation in populated area with more varied topography, we recommend the modeler to use a cell size smaller than the average building size in the study area.This is because in a numerical scheme with rectangular cell type, cells smaller than the average building size provide a better description for specific features in settlements such as building, roads, river embankment, and etc.
Let's now consider the significant differences between the models.If the purpose is to make estimates of future events, then due to the high uncertainty of the potential bulding damage by the tsunami attack, the CRM results should be interpreted as the maximum extent of that can be reached by the tsunami inundation.As for the TM and the ERM, because we ran the model for 2 h simulation time, the results should be interpreted as the maximum condition that can be reached within the first two hours of tsunami attack if no buildings were destroyed by the tsunami.They can be utilized to identify the distribution of inundation parameters in residential area associated with the human and infrastructure vulnerability assessment.However, the maximum extent of tsunami inundation produced by them should not be used for further mitigation purposes as likely smaller than expected due to the resistance force of the building works along the simulation.
By integrating the tsunami casualty algorithm into the tsunami run-up model, we have identified the spatial level of tsunami risk on human through the values of Tsunami Casualty Index (TCI), and have established the road risk map in the study area.Further research on evacuation model is needed for this analysis to become more useful in determining the percentage of the community that can complete the evacuation before tsunami comes, in estimating the number of victims, and in planning the appropriate countermeasures.

Fig. 1 .
Fig. 1.Tsunami source scenario (left hand side, Natawidjaja et al., 2009) and the calculated sea floor displacement (right hand side).Inset on the right hand side is the position of numerical tide gauges, and the computed tsunami wave forms at the bottom figure.White dotted lines are the historical faults, indicated by year in italics(Natawidjaja et al., 2006); transparent rupture zones are the result of recent seismic activity, indicated by year in bold(Sieh et al., 2008).

Fig. 2 .
Fig. 2. Computational area from region 1 to 5 in a nested grid system.

Fig. 3 .
Fig. 3. Digital Terrain Model for topographic data in CRM.

Fig. 4 .
Fig. 4. The process of input data preparation for Topographic Model.

Fig. 5 .
Fig. 5. Schematic descriptions of TM and ERM in residential area, (A) TM changes the direction of the flow due to the existence of the building.(B) ERM will significantly reduce the flow velocity (dotted line) when the tsunami passes through the building as the effect of the equivalent roughness.

Fig. 6 .
Fig. 6.Descriptive explanation of the force exerted by the tsunami flow.

Fig. 7 .
Fig. 7. Open channel experiment scheme (upper figure) conducted by Aburaya and Imamura (2002), and the resulting comparison with the ERM model (bottom figure) in front of the first obstacle.

Fig. 8 .
Fig. 8. Description of the ERM application if the building size is larger than the grid size.

Fig. 9 .
Fig. 9.The process of input data preparation for Equivalent Roughness Model.Color legend in panel D refers to the entire roughness data in the middle figure on panel E.

Fig. 10 .
Fig. 10.Schematic description of human body and the type of fall condition: (A) Slipping fall type, and (B) tumbling fall type.

Fig. 11 .
Fig. 11.Influence of the friction coefficient on the human feet soles (µ) and of the distance of the projection of the center of mass to the bottom back of heel (I G ) to determine the critical condition of human instability.

Fig. 12 .
Fig. 12. Verification of the proposed model using data from a hydraulic experiment, upper figure shows the schematic wave tank used in the experiment by Endoh and Takahashi (1995).Bottom figure is the verification result with present model.Solid line is the critical condition for the slipping fall type obtained from the model.

Fig. 13 .
Fig. 13.Comparison of the maximum flow depth.From left to the right hand side: CRM, TM, ERM (5 m), and ERM (15 m).Cross-sections and control points on TM are used for elevation check in Figs. 14 and 15.Dashed line in CRM and in the ERM (15 m) indicate the maximum inundation line of the ERM, while the solid line in CRM is the maximum inundation line from the TM.

Fig. 14 .
Fig. 14.Small seawalls and groins along the coast of Padang city (top figures), and the reduction of tsunami flow depth because of their existence (bottom figure).

Fig. 15 .
Fig. 15.Run up height and velocity comparison between TM and ERM from the cross-section 1 (left hand side), and the cross section 2 (right hand side).The circles denote the road junctions and corresponding values of inundation height and flow velocity.
The ERMs give results similar to each other at the second intersection.Even though the values produced by the ERM A. Muhari et al.: Examination of three practical run-up models for assessing tsunami impact

Table 2 .
Dimensions of the human body model derived from anthropometric data (in m).