Tsunami inundation modelling based on detailed roughness maps of densely populated areas

An important part within the German-Indonesian Tsunami Early Warning System (GITEWS) project was the detailed numerical investigation of the impact of tsunamis in densely populated coastal areas of Indonesia. This work, carried out by the German Research Centre Geesthacht (GKSS), in co-operation with DHI-WASY, also provides the basis for the preparation of high resolution hazard and risk maps by the German Aerospace Center (DLR). In this paper a method is described of how to prepare very detailed roughness maps for scenario computations performed with the MIKE 21 Flow Model FM in three highly resolved ( ∼10 m) priority regions, namely Kuta (Bali), Padang (West-Sumatra), and Cilacap (southern coast of Java). Roughness values are assigned to 43 land use classes, e.g. different types of buildings, rural and urban subareas, by using equivalent coefficients found in literature or by performing numerical experiments. Comparisons of simulations using differentiated roughness maps with simulations using constant values (a widely used approach) are presented and it is demonstrated that roughness takes considerable influence on run-up and inundation. Out of all simulations, the results of the worst case scenarios for each of the three priority areas are discussed. Earthquakes with magnitudes of MW=8.5 or higher lead to considerable inundation in all study sites. A spatially distinguished consideration of roughness has been found to be necessary for detailed modelling onshore. Correspondence to: G. Gayer (gerhard.gayer@gkss.de)


Introduction
The Sumatra earthquake of 26 December 2004 was one of the largest ever detected rupture in the Earth's crust.Only a few minutes after the earthquake the first tsunami waves hit the coastline of Northern Sumatra.
Shortly after the tsunami disaster, Germany offered technical support for the installation and implementation of a tsunami early warning system in the Indian Ocean.In a joint cooperation between Germany and Indonesia a tsunami early warning system has been established in Indonesia.The implementation was mostly completed in 2009.Optimization will be necessary, operation and maintenance will continue.The German-Indonesian activities are fully integrated into the overall UN plans and strategies for the establishment of global and regional early warning systems (Rudloff et al., 2009).
The German-Indonesian Tsunami Early Warning System (GITEWS) consists of a variety of components.It includes a seismological network consisting of broadband seismometers which rapidly localizes the earthquake and determines its strength, as well as GPS stations monitoring the deformation of the ground.It also contains a network of GPS buoys additionally equipped with ocean bottom pressure sensors and a tide gauge network to detect sea level changes.The respective sensors are connected by satellite communication to the Early Warning and Mitigation Center operated by the Indonesian Meteorological Climatological and Geophysical Agency (BMKG) in Jakarta.
The main scientific and technological challenge for the setup of an Early-Warning System in Indonesia is the tectonic setting of the so-called Sunda-Arc-Structure, an active continental margin almost parallel and close to the Indian Ocean coastline of Indonesia resulting in tsunami arrival times of about 30-40 min after the occurrence of an earthquake.Therefore complete new technologies and scientific concepts have been developed to reduce early-warning times down to 5-10 min (Rudloff et al., 2009).
Because such short early-warning times do not allow numerical computations if required, a data base of thousands of pre-computed Indian Ocean Tsunami scenarios was established and substantially compiled from simulations performed by the Alfred-Wegener Institute of Polar and Marine Research (AWI) and the German Research Centre for Geosciences in Potsdam (GFZ).
In the GITEWS project, the task of GKSS and its partner institution DHI-WASY was to calculate the detailed tsunami run-up and inundation -for all of those scenarios with wave heights higher than 1 m at the coast -in three priority areas, namely Padang (Sumatra), Cilacap (Java), and Kuta (Bali) to provide an analogues basis (300 scenarios for Cilacap, 91 for Padang, and 137 for Kuta) for the preparation of high resolution hazard and risk maps by the German Aerospace Center (DLR).
A general model description is out-lined in Sect. 2. The model set-up, bathymetry and mesh characteristics, boundary and initial conditions, are described in Sect.3.Because of the importance of bottom roughness, a method to establish detailed roughness maps is presented in Sect. 4. In Sect. 5 the impact of variable roughness on run-up and inundation is discussed.Finally, in Sect.6 results of the worst case scenario for each of the three priority areas are shown.

Hydrodynamical model
The model used for the calculation of the tsunami wave propagation and subsequent run-up in coastal areas was MIKE 21 FM (flexible mesh module) developed by DHI (2008a,b).Based on the numerical solution of the nonlinear formulation of the two-dimensional shallow water equations (the depth-integrated incompressible Reynolds averaged Navier-Stokes equations) the model considers main factors important for the calculation of tsunamis such as bottom shear stress, Coriolis forces, momentum dispersion, flooding and falling dry.
In MIKE 21 FM the depth integrated non-linear shallow water (NLSW) equations are used including lateral stresses.These depend on viscous friction, turbulent friction and differential advection and are estimated using an eddy viscosity formulation based on depth averaged velocity gradients.The bottom shear stress is determined by a quadratic friction law, where the friction coefficient can be expressed with the Manning number or the Chézy number (see Sect. 4).
The spatial discretization of the NLSW equations is performed using a cell-centred finite volume method.The spatial domain is discretized by subdivision of the continuum into non-overlapping triangles.An approximate Riemann solver is used for the computation of convective fluxes, which allows handling of discontinuous solutions.For the time integration an explicit scheme is used.For the majority of all scenarios, a time step out of the range of 0.01-10 s was sufficient.In cases of numerical problems (especially, when high gradients of wave heights or flow velocities occurred) the time step range was reduced to 0.001-1 s.

Model set-up, boundary values and initial fields
A variety of data was used to establish the general model bathymetries/topographies. The bathymetries in shallow areas were derived from echo sounder measurements by the Agency for the Assessment and Application of Technology, Indonesia (BPPT), DHI-WASY, and University of Hanover.Furthermore, navigational charts (C-Map) and the General Bathymetric Chart of the Ocean (GEBCO) were used.The topographies were based on digital surface or terrain maps, street maps and building maps, provided by DLR, and differential GPS measurements performed by DHI-WASY.
Table 1 gives an overview of the priority areas with respect to size of the computational domain, numbers of mesh elements/nodes and mean mesh element sizes.
For the model area of Padang (Sumatra), the mean onshore terrain height is 12.41 m, but most parts of the city between beach and eastern mountain range have lower elevations in the order of 5 m.Just at the foot of the mountains there is a sharper increase of elevations to 20 m.
The city of Cilacap (Java) is sheltered in the south by a peninsula with steep coasts.In between, a shipping channel (∼30 m deep in the eastern waters and up to 19 m near the southern harbour) allows the tsunami waves to intrude to the west and after a sharp turn to the north (here, the channel has water depths less than 5 m).The general structure of the city ground (mean onshore terrain height = 6.55 m) consists of ridges (up to 8 m-9 m high) and valleys in between (2 m-3 m), parallel to the curved eastern coast.The first ridge along this coast has an average elevation of 4 m-5 m.
Kuta (Bali) is characterized by an airport area of low elevation (3 m-4 m).At the western coast and north of the airport area the terrain rises sharply up to 5 m-6 m, followed by a landscape with small hills (12 m high) and valleys some of them only 4 m-5 m a.m.s.l., allowing tsunami waves (higher than 6 m) to intrude far inland.Greater parts of the peninsula in the South-East might also get flooded from easterly directions but, together with a very shallow basin covered partly by mangroves this peninsula reduces the impact on Kuta.
In order to identify minimum requirements for mesh resolutions from a physical point of view, preliminary investigations have been done.Offshore mesh resolutions have been determined with the help of inverse solutions from tide gauge records for the 17 July 2006 Java Earthquake tsunami (Fujii and Satake, 2006), presented in Kongko et al. (2008).Onshore this topic is considered in Leschka et al. (2009a,b) based on the same event.The size of each of the computational domain reflects differences between the study areas.For example, the model for Kuta included areas of North Bali as well as part of the coast of Java in order to take into account possible wave reflections.These areas have been represented with a coarser mesh resolution than the area of interest itself.
In less sensitive parts (e.g.deep ocean or higher grounds never getting flooded), the mesh resolution of the model was approximately 1 km, which was sufficient for our purpose.The resolution was refined as approaching the area of interest.There, mesh element lengths down to 10 m enabled the mapping of finer structures such as sand banks, streets and buildings.
At open boundaries the models are driven by time series of sea surface elevations, which were extracted from the Indian Ocean tsunami scenario data base.The time series (delta = 1 min) were interpolated from mesh nodes of AWI's numerical mesh to straight lines with spatial resolutions of 100 m.During a model run, these input time series were interpolated to our mesh nodes and time steps, respectively.
Depending on the direction of the approaching tsunami wave, the time series of sea surface heights are taken as driving force at appropriate open boundaries.At all other boundaries, the condition is set to zero velocity gradients across these boundaries, enabling outflows but no inflows.
In addition to boundary values, each scenario needs an initial field of the sea surface heights (ssh) and a new bathymetry/topography, because due to the earthquake the bottom and the water columns above had moved vertically (thereby changing the distance to the mean sea level reference line).Both can be derived most easily from the initial bottom displacements (also taken from the tsunami scenario data base) (Babeyko, 2008).Interpolated to our mesh, the bottom displacements were added to the general bathymetry/topography, resulting in an adjusted new one, specific for a certain scenario.According to the bottom displacements, the water columns above were shifted simultaneously.With respect to mean sea level, the initial field of sea surface heights was formed.These displacements are not restricted to the area of the epicentre.The movement of the bottom plate at the source might also trigger the movement of adjacent plates.In some cases, this also has influence on distant coastal areas.For example, the impact of an earthquake of magnitude 9 with its epicentre ∼150 km to the south of Bali was still noticeable at the coast of Southern Bali.Bottom and sea surface were lowered by up to 0.9 m, instantaneously generating a wave trough.

Roughness coefficient maps
Because of high mesh resolution and physical effects (e.g.wave attenuation) the bottom roughness (in terms of Manning coefficients) was of main interest.For example shallow reefs, dykes, shoreline stabilisations, forests, or buildings might change the propagation of waves considerably due to being solid obstacles or due to roughness induced energy dissipation and therefore should be included in the final model mesh either as features with appropriate elevation heights or as adequate roughness elements mapped to the ground.
Generally, in case of 3-D features (e.g.buildings, trees), several options are available how to take them into account: 1.As obstacles, forcing the water to flow around.In this case, mesh elements covered by these obstacles would simply be declared non-calculation elements, disregarding any effects caused by the roughness at the sides of the obstacle.
2. According to their height, features are included as land elevations.This would make sense in cases where features are overtopped.
3. Substituted by ground elements with a certain roughness (e.g. a mangrove forest, disregarding the heights of the trees).
4. As a combination of 2 and 3 (e.g. a harbour wall with a roughness element placed on top of it).
After careful inspection it has been decided for option 3, because this allowed taking into account material structures.collapse.This differentiation would not be possible with option 1, when buildings are excluded from the modeling process.Also it has been distinguished between different vegetation forms: for example, a dense mangrove forest exerts a much higher roughness than palm trees in a park.Option 4 was chosen only for single solid objects like dykes and channels.In this case it has been ensured that the shapes of these objects were preserved by specifying the grid elements and their elevations manually.
For each of the priority areas, roughness maps (Manning values at each node of the computational mesh) were established in three steps: 1. Localising areas with a specific land use.

Identifying typical roughness elements and land cover
characteristics for each land use, and 3. assigning a Manning coefficient for each land use.
Areas of land use classes were localised with the help of maps provided by Bakosurtanal (National Coordinating Agency for Surveys and Mapping of Indonesia), Lapan (Indonesian National Institute of Aeronautics and Space) and DLR.Since the land use classification of these maps was quite general in some cases (e.g.airport area), satellite images (Quickbird), street maps, and building maps were used for refined subdivisions (e.g.airport asphalt, airport meadow, airport settlement, airport terminal).Figure 1 presents the airport area of Kuta, Bali.There, yellow lines mark areas with specified land use classes.
Typical roughness elements are for example buildings, shrub, lakes, mangroves, forest.Typical land cover characteristics are meadow, paddy fields, garden, sand (beach), bare soil or asphalt.Representative Manning values can be found in literature, e.g.Acrement et al. (1989), Latif and Hadi (2007) or Weichel et al. (2007).With the help of satellite images, roughness coefficients were assigned to delimited areas depending on averaged fractions of roughness elements contained in the land use areas as classified in Table 2.For example, an area with land use class "airport meadow" was uniformly covered to 100% with only one type of roughness elements, namely meadow, whereas for non-uniform areas with land use type "hotel" (Fig. 2) the average percentages of main roughness elements were estimated as 40% coverage with buildings, 40% coverage with meadows, using a Manning value M=40 m 1/3 /s (Kouwen, 1992), and 20% coverage with trees, represented by M=14.3 m 1/3 /s (Latief and Hadi, 2007).
Manning values for such non-uniform land use classes then have been calculated as the coverage weighted average of values specific for each of the roughness element fractions.
Not all types of roughness elements (and their Manning values) occurring in the land use classes could be found and had to be determined otherwise.In case of blocks of houses with streets in between, an appropriate value was determined by performing a series of numerical scenarios for an artificial city as shown in Fig. 3. Inundation results of a run with stable houses, represented by non-calculation mesh elements were compared with inundation results where houses were represented by roughness elements.According to typical slopes and buildings found in the study areas the set-up of the numerical grid (covering an area of 4 km×0.8km) was as follows: a spatial resolution of 10 m, a bathymetry slope from a water depth of -40 m to 0 m (over a distance of 2 km), an adjacent beach elevation of 1 m, an elevation slope in the city from 2 m to 7 m (over a distance of 1 km), and from 7 m to 12 m behind the city (over a distance of 1 km).
The city itself started at the beach side with a row of bigger houses (representing hotels), followed by blocks of smaller houses (width 10 m) and single bigger houses (width 20 m).All the larger streets were parallel or perpendicular to the beach, some of the smaller streets could cross in an arbitrary manner.No distinction was made concerning the quality of the houses and streets.Generally, the Manning value was set to M=32 m 1/3 /s (coarse sand) on the sea side and in mesh cells not covered by houses.
Figure 4 shows four final results of (from top to bottom) case a, a scenario without any city, included only to show the difference (and the importance to take roughness into account) to the other cases, of case b, the reference run with houses as non-calculation cells (solid objects), and of cases c and d, derived after a series of simulations, varying the Manning value until the same inundation line was reached as in case b.Case c was included to estimate a Manning value M=7.7 m 1/3 /s for areas for which only general information like "covered with stable houses" would be available, but no detailed information about the spatial distribution and size of houses.Case d resulted in a very high roughness of M=2.5 m 1/3 /s for the elements representing houses.
Furthermore, during a tsunami event, water might cause a building to collapse, depending on water level, current speed, floating objects on one hand and building stability on the other hand.Presently, the model does not take these criteria into account explicitly.For simplified application purposes it has been estimated that 50% of normal houses would partly collapse during a big tsunami event.Since the Manning coefficient for collapsed buildings was unknown, this type was represented by an equivalent Manning-Strickler value of k S =20 m 1/3 /s for big rocks in a mountain torrent with supercritical flow (Schneider, 1992).Using the relation 1/M = n = k 1 / 6 S /25.4 (DHI, 2008a) led to a Manning value of M=15.2 m 1/3 /s.Finally the coefficient for buildings was determined as M=11.1 m 1/3 /s.
In case of land use class "hotel areas" (buildings are taken into account to 40%) the weighted average of all roughness values resulted in the final Manning value of M=25 m 1/3 /s.This kind of evaluation resulted in roughness maps of Manning values for each of the priority areas.As an example, Fig. 5 shows the final roughness map of Kuta.Blue colours depict high roughness (e.g.mangroves, shrub), yellow colours represent mean roughness (e.g.sand, rice fields) and red colours represent low roughness (e.g.lakes, streets, and airport).Mangrove swamps (blue colour) along the eastern coast help to protect the land from high flooding.Further information on roughness map generation with the focus on the area of Cilacap is given in Leschka et al. (2009c).

The importance of roughness: comparisons
In Fig. 6 two model runs of the same scenario but different roughness are compared.The picture to the left shows the maximum sea surface heights (sea side) and total water     depths (in Kuta) of a simulation using the detailed roughness map shown in Fig. 5. On the right hand side a constant roughness with a Manning coefficient of M=32 m 1/3 /s for the whole model area is assumed.This value represents land covered with coarse sand and is used in many applications when no detailed roughness information is available.In both simulations the wave heights at the western shore are almost identical.With the exception of the airport area, the colour patterns suggest that differences in water depths in flooded areas seem to be low.Locally, differences of 0.5 m or less occur (see red circles in Fig. 6).In general, in the first simulation the inundated area is smaller, suggesting on average a higher roughness and slowed down water flows as shown in Fig. 7.In the airport area, the roughness used in the first scenario is much lower than in the second one, enabling very high flow velocities of about 10 m/s.
In the previous chapter it has been shown that the Manning value for areas for which only a general information like "covered with stable houses" is available (numerical experiment case c) is M=7.7 m 1/3 /s, a roughness value comparable to that of Mangrove forests.Figures 8 and 9 show the inundation results and comparisons for the same scenario as discussed above, indifferently applying this high roughness value everywhere on land.Expectedly, the differences (compared to the model run using a differentiated roughness map) are huge.
Of course the inundation not only depends on roughness but also on local topographical conditions.Therefore this kind of comparison is repeated for a scenario affecting the city of Cilacap on Java.Again, differences are most pronounced for the case of using a Manning value of M=7.7 m 1/3 /s everywhere.Here, only the comparison of velocities (Fig. 10) is shown.As can be seen, not only that there are big differences in inundation but also that flow velocity distributions differ significantly.Whereas in case of the simulation with constant roughness coefficient (right panel) the velocities are rather low (∼2 m/s at the eastern coast, decreasing gradually further in land), the velocities of the simulation using a detailed roughness map show a high variability (left panel, red stripes ∼10 m/s, yellow stripes ∼7 m/s, dark green areas ∼4 m/s), in agreement with topographical features and high variability of roughness.

Worst case scenarios
Using highly resolved bathymetries/topographies (Sect.3) and differentiated roughness maps as described in Sect.4, the run-up and inundation was modeled in the priority areas Padang, Cilacap, and Kuta.The results of 91 scenarios for Padang, 300 for Cilacap, and 137 for Kuta provided input data for high resolution hazard and risk maps (prepared by DLR).All scenarios were selected from the Indian Ocean tsunami scenario data base and fulfilled the condition of wave heights higher than 1 m at the coast of the respective area.
In this chapter inundation results of worst case scenarios are shown to give an impression of the hazard in each of the priority areas.Probabilities of occurrence of such severe events and relations to goods like human beings and buildings are not discussed here.They are taken into account (estimated from historical records) and described in the process of establishing hazard and risk maps.Here, the worst case scenarios are presented to contribute high resolution computations when the accuracy of various tsunami model results will be discussed and validated.Figures 11,12,and 13 show the maxima of sea surface heights (maximum water elevation during a 3-hourly simulation) of the worst of all scenarios of earthquake magnitude 9.0, overlaid with inundation lines of the worst  of all scenarios of magnitude 8.5 (black lines) and of magnitude 8.0 (white lines) for each of the priority areas Padang, Cilacap, and Kuta, respectively.The scenarios are named according to their epicentres.The first digits give the location along the Sunda Trench (starting with 1 at the south-eastern end), the following digits give the location with respect to the As expected, the inundation increases with earthquake magnitude.Tsunamis generated by earthquakes of magnitude 8.0 are not high enough to cause remarkable inundation.Only lower grounds are flooded as in the case of the eastern coast of Kuta, an area mainly covered by mangroves (Fig. 13).In the case of Padang (Fig. 11), very local inundations are visible at river deltas and along rivers and channels (please note that river banks are not overtopped).
With increasing incoming tsunami heights as in magnitude 8.5 scenarios, the situation gets much worse as can be seen from the flooded areas indicated by black lines.In case of Padang, this line roughly corresponds to the 5 m elevation contour line.In Cilacap (Fig. 12), the water intrudes ∼2 km inland.
During worst case events of magnitude 9.0, almost entire cities are flooded.In Padang, the inundation reaches the 8 m contour line at the foot of the mountain range.In Cilacap and in Kuta, only a few spots (city areas on higher grounds, elevations can be estimated from adjacent ssh values) are spared from being flooded.

Conclusions
The great number of simulations of different tsunami scenarios enabled us to get a detailed insight of possible hazards during flood events in the priority regions Padang, Cilacap, and Kuta which were investigated by GKSS and DHI-WASY within the GITEWS project.Prerequisite for a successful simulation of tsunami run-up and inundation was the careful preparation of bathymetries/topographies of high resolution and detailed roughness maps to take into account the influence of structures smaller than features considered in the topography by using spatially varying Manning coefficients.
In this paper a method for detailed roughness map preparation has been described.Implications of using such differentiated roughness information were shown by comparing inundations with results of simulations using constant Manning coefficients like M=32 m 1/3 /s (a value commonly used in many applications) or M=7.7 m 1/3 /s (a value representing a settlement, derived from numerical experiments).Other Manning coefficients have been chosen based on literature.
It has to be noted, that in the pilot regions no data is available for validation purposes.A validation with the help of data sites, e.g.Banda Aceh, could not been fully performed due to the lack of highly detailed input data dated to the period before the 2004 Indian Ocean tsunami.In additional studies, an estimation of uncertainty for flux simulations due to different roughness representations for the area of Cilacap has shown that maximum water levels onshore vary considerably, but the difference regarding velocities is much higher (Leschka et al., 2009c).As one would expect from mean onshore terrain heights and their standard deviations in Table 1, a comparison of uncertainties for settlement areas in all study areas confirms, that flow depths uncertainties differ among the areas, depending on the slope of the terrain (Leschka et al., 2010).When interpreting the results it has to be mentioned, that other energy loss such as inertia and important effects such as debris flow to date are not fully understood and cannot be considered adequately for application purposes.It should also be noted that discretization technique and solution algorithm of the non-linear shallow water model has an influence on the results, which lies in the same order like the differences in roughness consideration onshore.Furthermore, temporal changes of the roughness during a tsunami event are not considered adequately with respect to real forces acting on structures.However, using detailed roughness maps for tsunami inundation modelling has been shown to be applicable to relatively large areas such as cities.Especially in terms of computed velocities onshore, the method increases the degree of detail within the results remarkably, which is required for hazard and risk mapping.Extensive work is still required in order to achieve a physically fully validated method to determine tsunami inundation in large areas.

Fig. 1 .
Fig. 1.Quickbird satellite image of airport area of Kuta, Bali.Yellow lines indicate areas with specific land use classes.Coordinate system WGS 1984 UTM 50 S.

Fig. 2 .
Fig. 2. Quickbird satellite image of typical hotel areas of Kuta, Bali.Blue lines indicate the hotel area, red lines mark streets.Coordinate system WGS 1984 UTM 50 S.

Fig. 3 .
Fig. 3. 3-D view of the artificial "roughness city", consisting of stable houses and streets.

Fig. 4 .
Fig. 4. Inundations (max.water depths): (a) No city.(b) Houses are non-calculation cells.(c) Unknown distribution of houses, Manning value of M=7.7 m 1/3 /s on the land side.(d) Houses are substituted by roughness elements with a Manning value of M=2.5 m 1/3 /s.

1684G.
Gayer et al.: Tsunami inundation based on detailed roughness maps

Fig.
Fig. Roughness map (Manning value M at mesh nodes) of Kuta, Bali.

Figure 7 :
Figure 7: Maximum velocities for the same scenario as in figure 6.

Fig. 7 .
Fig. 7. Maximum velocities for the same scenario as in Fig. 6.

Fig. 9 .
Fig. 9. Maximum velocities for the same scenario as in Fig. 8.

Table 1 .
Mesh characteristics of study areas.