Operational tsunami modelling with TsunAWI – recent developments and applications

In this article, the tsunami model TsunAWI (Alfred Wegener Institute) and its application for hindcasts, inundation studies, and the operation of the tsunami scenario repository for the Indonesian tsunami early warning system are presented. TsunAWI was developed in the framework of the German-Indonesian Tsunami Early Warning System (GITEWS) and simulates all stages of a tsunami from the origin and the propagation in the ocean to the arrival at the coast and the inundation on land. It solves the non-linear shallow water equations on an unstructured finite element grid that allows to change the resolution seamlessly between a coarse grid in the deep ocean and a fine representation of coastal structures. During the GITEWS project and the following maintenance phase, TsunAWI and a framework of pre- and postprocessing routines was developed step by step to provide fast computation of enhanced model physics and to deliver high quality results.


Introduction
The devastating tsunami in December 2004 triggered numerous activities aiming at the installation of a tsunami early warning system in the Indian Ocean.Since the time available for early warning is very short, especially along the Indonesian coast facing the Sunda Trench, first estimates of the tsunami impact after an earthquake are based on precomputed tsunami simulations.These can rely only on a model thoroughly tuned and verified following the guidelines based on Synolakis et al. (2007Synolakis et al. ( , 2008)).Furthermore, for warning systems in the Indian Ocean, the implementation plan (IOC, 2008) demands to benchmark the operational models according to these guidelines.
In the framework of the German-Indonesian Tsunami Early Warning System (GITEWS), the tsunami modelling group at AWI (Alfred Wegener Institute) has developed the operational model code TsunAWI and provided the warning centre at the Indonesian Agency for Meteorology, Climatology, and Geophysics (BMKG), Jakarta, Indonesia, with a repository of 3470 prototypic tsunami scenarios.Since the start of GITEWS in 2005, the warning system was constructed and improved step by step, integrating the daily experiences at the warning centre, and so was the TsunAWI code developed.In this article, we want to highlight some aspects of this process of TsunAWI developing from a research model to an operational code with a framework of preand postprocessing routines.The development of TsunAWI started as a spin-off of the ocean model FESOM (Finite Element Sea-Ice Ocean Model) (Wang et al., 2012b, and references therein), from which it inherited the main structure, the finite element discretisation, and some core routines.The first step was to drastically reduce the model physics, e.g. to eliminate the vertical dimension, the temperature and the salinity.On the other hand, inundation, that can be neglected in large scale ocean modelling, had to be implemented.
The code TsunAWI simulates all stages of a tsunami from the source to propagation and run-up by solving the nonlinear shallow water equations discretised with finite elements on a two-dimensional unstructured triangular grid.This discretisation is very flexible with respect to resolution and allows for an excellent representation of complicated coastlines and bathymetry.Other tsunami models also take the advantages of unstructured grids, though usually pared with finite volume discretisations (e.g.ANUGA, Roberts et al., 2007).As the price for the unstructured mesh is a less efficient implementation, other models rely on regular N. Rakowsky et al.: Tsunami modelling with TsunAWI grids, which may be nested to resolve areas of interest, the finite difference model TUNAMI-N3 (Tohoku University's Numerical Analysis Model for Investigation of Near-field Tsunami Number 3) (Imamura et al., 2006) being a well known example.
In Sect.2, the basics of TsunAWI are introduced, namely the underlying shallow water equations, the discretisation with finite elements, and some implementation issues.Section 3 presents results on the 2004 Indian Ocean tsunami.The tsunami scenario repository built for GITEWS is described in Sect. 4. Section 5 concludes with a look at ongoing and planned research.

The governing equations
The non-linear shallow water equations (SWE) The tsunami simulation code TsunAWI is based on the rotating non-linear SWE given by the vertically averaged equations of motion and continuity.We consider the boundary value problem in Cartesian coordinates (x, y) ∈ in the plane domain with boundary ∂ at the time t ∈ [0, T ] ∂v ∂t for the horizontal velocity vector v = (u, v) and the total water depth H = h + ζ > 0 as the sum of the unperturbed water depth h and the surface elevation ζ .Furthermore, ∇ = ∂ ∂x , ∂ ∂y denotes the gradient operator, f the Coriolis parameter, k the vertical unit vector, r the bottom friction coefficient, and K h the eddy viscosity coefficient.

Boundary conditions
On the solid part of the boundary ∂ 1 , the condition Oliger and Sundström (1978) derived (v, ζ ) | ∂ 2 = 1 with the operator of the boundary condition and a vector function 1 determined by the boundary regime, differing for inflow and outflow.In practice, the full information on the open boundary is unavailable, and one option is to impose the boundary condition on the elevation ζ | ∂ 2 = (x, y, t).The accuracy of this approach was analysed by Androsov et al. (1995).For most tsunami simulations, the radiation boundary condition is implemented in TsunAWI providing free linear wave passage through the open boundary, given the Coriolis acceleration plays only a small role (for a review of several open boundary conditions see Røed and Cooper, 1986).Furthermore, as TsunAWI operates on unstructured grids, the open boundary can be replaced by a coarse global grid with a finely resolved focus on the region of interest.

Initial conditions
Initial conditions can be taken from various sources, e.g.Okada parameters (Okada, 1985) from geological surveys for hindcasts, or prototypic ruptures for the Sunda Arc provided by the German Research Centre for Geosciences (GFZ) for GITEWS (RuptGen, Babeyko et al., 2008, 2010).
Especially for hindcasts, TsunAWI allows to elevate the ocean bottom and surface repeatedly within several time steps, following the time frame of the source model.Usually, only an initial sea surface height ζ 0 is given, while the velocity is assumed to be zero, but it is also possible to supply an initial velocity field additionally, e.g.derived from a land slide model as provided by Androsov and Voltzinger (2005).

Bottom friction
Manning bottom friction with a constant parameter is chosen for most calculations.Alternatively, a value varying with the surface structure can be provided.

Finite element method
Finite element methods are widely used in studies of wave generation and propagation in different fields of fluid dynamics.They are often employed to simulate propagation of long waves such as ocean tides and tsunamis in the ocean in the framework of shallow-water equations, e.g.Kienle et al. (1987); Greenberg et al. (1987); Baptista et al. (1993).The main reason to choose finite elements is the computational mesh that can be adapted to cover basins with irregular bottom topography and coastlines.The spatial discretisation of TsunAWI is based on the finite element approach by Hanert et al. (2005) with modifications like added viscous and bottom friction terms, corrected momentum advection terms, radiation boundary condition, and nodal lumping of mass matrix in the continuity equation.The basic principles of discretisation follow Hanert et al. (2005) with linear conforming elements P 1 for sea surface height ζ and water depth h, and linear non-conforming elements P NC 1 for the velocity v. Thus, values for sea surface are computed at element nodes, while the velocity components are regarded at the edges.

Leap-frog time stepping scheme
Simulation of tsunami wave propagation benefits from an explicit time discretisation.Indeed, numerical accuracy requires relatively small time steps, which reduces the main advantage of implicit schemes.Furthermore, modelling the inundation processes requires very high spatial resolution in coastal regions and consequently a large number of nodes, drastically increasing computational demands for an implicit time discretisation.
The leap-frog discretisation was chosen as a simple and easily implemented method.We rewrite Eqs. ( 1) and (2) in time discrete form with the time step length t and the time index n.The leapfrog three-time-level scheme provides second-order accuracy and is neutral within the stability range.This scheme, however, has a numerical mode, which is removed by the standard Robert-Asselin filter.Notice that friction and viscosity contributions deviate from the usual leap-frog method.

Momentum advection scheme
Because of the discontinuous character of velocity representation, special care is required with respect to the implementation of the momentum advection.Earlier experiments with P NC 1 − P 1 code revealed problems with spatial noise and instability of the momentum advection when the discretisation is used as described by Hanert et al. (2005).A modified implementation without upwinding terms was found to work well when paired with rather high viscous dissipation to remove small-scale noise in the velocity field.In addition, the implementation of the momentum advection for P NC 1 velocities involves cycling over edges in the numerical code, in addition to cycling over elements to assemble the elemental contributions.This leads us to a simpler approach, which provides a smoothing of the velocity field while removing edge contributions.
According to this approach, prior to calculating the advection term in the momentum equation, we project the velocity from P NC 1 to P 1 space in order to smooth it.This projection becomes numerically efficient by nodal quadrature (lumped mass matrix).The projected velocity is then used to estimate the advection term.Finally, the result is multiplied with a P NC 1 basis function and integrated over the domain (Maßmann et al., 2008).This results in a very stable behaviour.

Viscosity
The coefficient K h for the horizontal viscosity is determined by a Smagorinsky parameterization (Smagorinsky, 1963).Optionally, the standard linear viscosity K h = c lin x y may be chosen for comparison or to save computation time.The advantages of Smagorinsky viscosity are demonstrated in the simulation of the Okushiri channel experiment test case (Androsov et al., 2008).

Wetting and drying
For modelling wetting and drying, we use a moving boundary technique based on the "dry node concept" developed by Lynett et al. (2002).The idea is to exclude dry nodes from the solution, but to extrapolate the elevation to the dry nodes from their wet neighbours.As this scheme is neutrally stable, it requires horizontal viscosity in places of large gradients of the solution.The Smagorinsky viscosity (7) fulfills this task while keeping the dissipation on a level that does not affect the quality of the solution.
The implementation in TsunAWI showed that in order to retain stability it has to be decided which values to extrapolate and how.Finally, linear extrapolation was chosen for the gradient ∇ζ and the sea surface height ζ at dry elements and nodes.Furthermore, all other values at dry nodes, elements, and edges are excluded from further computations.
In the first versions of TsunAWI, the extrapolation scheme was not always stable at complicated coastlines when the mesh did not resolve the geography smoothly enough.This could lead to local numerical artefacts like spurious spikes in the amplitude or initial waves of high frequency in the far field of the source.Due to the fractal nature of many coastlines, this problem cannot be solved by mesh generation alone and was handed over to the postprocessing routines.With growing experience, the problem could be solved for TsunAWI, which now runs stably if the time step and mesh resolution are chosen appropriately.

Mesh generation
The quality of the triangulation of the model domain is crucial for the model results.The meshes used in the following studies were generated with a mesh generator based on the freely available software Triangle by Shewchuk (1996).Starting from a model domain defined within a topography/bathymetry data set (in our case GEBCO 30"); Triangle generates a mesh based on a refinement rule depending on the water depth and prescribed by the corresponding wave phase velocity and the Courant-Friedrichs-Lewy criterion (CFL).The triangulation will be refined until the edges in (8) As the Triangle output is not smooth enough for numerical experiments, several smoothing iterations are applied.Each iteration consists of edge swapping, torsion smoothing, and linear smoothing.Torsion smoothing tries to adjust the angles around each node, linear smoothing acts on the distance between nodes.These strategies are described in Frey and Field (1991).

Code optimisation and parallelisation
The operational version of TsunAWI is parallelised with OpenMP (Open Multiprocessing) by defining parallel regions in the numerically most demanding parts of the code and splitting up the corresponding loops, thus sharing the load to the cores involved.The implementation in this approach needs smaller changes than a fully parallel MPI (Message Passing Interface)-implementation, and proves to be numerically efficient, if data locality is ensured.
In fact, already in a serial code, data locality is crucial for the performance.In particular the representation of variables in an unstructured mesh, if not constructed carefully, can result in many cache misses during calculation and considerably slow down the code.For TsunAWI simulations we therefore order the variables along a space filling curve (SFC, for triangular meshes with regular structure see Bader et al., 2008), which guarantees good data locality on all levels of the memory hierarchy by construction.The ordering is visualised in Fig. 1, the effect on the computation time is shown in Table 1.Results for one commonly used ordering technique, reverse Cuthill-McKee (RCM, George and Liu, 1981), are given for comparison.

Pre-and postprocessing
Preprocessing for TsunAWI runs is dominated by producing a smooth computational grid that resolves the domain in the desired accuracy.In the first model runs, the time step length is tuned to be as long as possible while keeping the scheme stable.For a small number of runs in one setting, e.g. for a hindcast, the preliminary steps pass over into the main simulation.However, for repository calculations with several thousand scenarios, scripts have to be provided to automate the runs on large Linux clusters.While in the early days of GITEWS, TsunAWI computation was expensive compared to data storage, the ratio is now vice versa.We therefore now reduce the amount of data beforehand, e.g. by deciding at which locations a timeline is really needed in addition to global values like maximum amplitude, arrival time, and maximum flux.
TsunAWI output is written in netCDF (Common Data Format), which is standard in oceanography, but cannot be displayed by GIS applications.A set of postprocessing scripts  and routines extracts the data products and converts them to shapefile format needed by the warning system.Special care has to be taken for isochrones, because for smaller tsunamis, a fixed threshold can result in a rather disturbing picture as illustrated in Fig. 2. Though it is correct from the physical point of view, a map like this cannot be distributed to the media.In TsunAWI postprocessing, Fig. 2. Arrival time with a threshold of 2 cm for a minor tsunami scenario, M w = 8.0.In the regions marked in black, the fixed threshold is never reached, at some, it is reached by the second or third wave crest.
we therefore chose a relative threshold of 10 % of the local maximum amplitude, see Fig. 3.Still, islands and other complicated geographical settings can produce discontinuous arrival times and there is no other way but checking the isochrones of each scenario manually.

Application: hindcast of the Indian Ocean tsunami in December 2004
TsunAWI is still under development and the performance and model results must be constantly evaluated to ensure a consistent and stable code.Validation needs to be performed on several levels, like consistency and convergence, benchmarks in idealised settings with well defined results as proposed by Synolakis et al. (2007Synolakis et al. ( , 2008)) For the tsunami generated by the great Sumatra-Andaman earthquake on 26 December 2004, a large amount of data are available.In this section, model comparisons with observations from satellite altimetry, tide gauge records, and field surveys are discussed.The results are closely related to the studies in Harig et al. (2008), where additional information can be found.

Topography and bathymetry data
Topography and bathymetry data are based on several sources.The General Bathymetric Chart of the Oceans (GEBCO) data are globally available with 30 arc seconds resolution and are used as initial topographic and bathymetric information.It is replaced by more precise data wherever available.Bathymetric data are locally improved by ship measurements and nautical charts.Topographic information is improved by the Shuttle Radar Topography Mission (SRTM) data, which are freely available at a resolution of 90 m, however the vertical accuracy is usually not sufficient for model studies with inundation.For Indonesia, the SRTM data were additionally processed by the German Aerospace Agency (DLR) and provided at a resolution of 30m.In the Aceh region, additional bathymetric and topographic data were provided by BPPT Jogyakarta.

Source model
In tsunami modelling, the exact knowledge of the source, i.e. the initial conditions of the model, is crucial for comparisons with data.Usually, the source parameters are optimised with respect to a choice of all available measurements.The source model used in the following studies is based on the results of Tanioka et al. (2006).Their objective was to optimise the subfaults of the earthquake with respect to certain tide gauge records.The resulting subfaults are shown in Fig. 4. To demonstrate the impact of source parameters on the matching with data, the orientation of the southern faults has been modified as indicated in that figure.The strike angle of subfaults A/C is changed from 340/340 • as proposed by Tanioka et al. (2006) to 290/320 • .This change is motivated by Hirata et al. (2006), who optimised source parameters with respect to satellite altimetry.The geometry differs from Tanioka's source.The subfaults E1 and E3 vaguely correspond to subfaults A and C in Fig. 4 and have strike angles of 290 and 330 of the 2004 tsunami optimised with respect to selected data.
In the present study we are not aiming at source optimisation, rather we would like to highlight the sensitivity of model results with respect to small changes of source parameters.

Model set-up
The simulation code TsunAWI was employed on a computational mesh with 5 million nodes and 10 million triangles.
The whole Indian Ocean is covered with a resolution range between 15 km in the deep ocean and 40 m in the Aceh region in northern Sumatra, see Fig. 5.The finest resolution enforces a time step of 0.5 s, and the model is integrated for 10 h.This rather short integration time was just enough to compare with most measurements, and the focus was to perform a larger number of runs for different combinations of finite fault solutions and model parameters.In the positions  along the satellite tracks, the model state is written to a file every second to allow for a comparison between model and satellite altimetry data.Additional information on model setup and initialization can be found in (Harig et al., 2008).

Satellite altimetry
The Indian Ocean tsunami was the first tsunami to be observed by satellites.Jason 1 (J1) and Topex/Poseidon (T/P) were above the bay of Bengal about 2 h after the earthquake.By comparison to the previous average heights, the signature of the wave was determined.These data are available from NOAA (Smith et al., 2005).Table 2 summarises the cycles that were taken into account.Figure 6 contains the groundtracks of J1 and T/P, whereas Fig. 7 displays the extracted tsunami signal.The model results are interpolated in time and space and extracted from a Hovmöller diagram as shown in Fig. 8.The J1 signal in Fig. 7 shows a double peak in the leading wave crest, which is due to the partial waves generated by the southern subfaults, whereas in the T/P measurements these partial waves overlap.This behaviour is reproduced by the model, however, the matching strongly depends on the fault parameters.Changing the strike angles in fault A from 340 to 290 and in fault C from 340 to 320 improves the matching and lowers the root-mean-square (RMS) errors as indicated in Table 2.However, it has to be noted that Allen and Greenslade (2012) have recently pointed out that the RMS error is not necessarily the best parameter to compare model results, because it may be to sensitive to outliers.

Inundation
After the event in December 2004, several field surveys examined the run-up and maximum flow depth in the affected regions.In the area of Banda Aceh, which was most heavily hit, eight locations with well documented field measurements were chosen to be compared with model run-up.As flow  depth and run-up depend heavily on the prescribed roughness parameterization, several experiments with identical initial conditions and different values for the Manning parameter were conducted.In all cases, a constant manning number was applied in the whole domain.
Figure 9 shows the locations of measurements in the Aceh region, the boundary of inundation as it was obtained from satellite images (provided by DLR), and a TsunAWI simulation result for the modified source with strike angles 290/320 in subfaults A/C.Here, the best fitting Manning roughness parameter n = 0.03 was used.This value was determined by comparing the corresponding RMS errors for experiments with varying parameter n, as visualised in Fig. 10.The effect  of the different strike angles in the southern subfaults on the inundation simulation RMS error is shown in Fig. 11.Also in this respect, the modified orientation improves the results considerably.In fact, the unsatisfactory inundation result obtained with the original source was the reason to test the strike angles proposed by Hirata et al. (2006) for the southern subfaults, because they have the largest impact on the inundation in Banda Aceh.

Tide gauge records
The Indian Ocean tsunami was registered by tide gauges worldwide.Since arrival time and estimated amplitude are among the most important warning products, good matching between model results and tide gauge measurements of arrival times and the height of the leading wave crest is important.Figure 12 summarises the time series in some locations throughout the Indian Ocean.The comparisons between tide gauge records and model results generally match well with respect to arrival time.The height of the leading crest is sometimes underestimated, however, in most of the locations the agreement is very good.The orientation of the southern subfaults does not influence the results at far distances such as Salalah and Lamu, and the matching in Colombo is slightly better with the uniform values (340/340).This is consistent with the derivation of these strike values, as this station was used for optimisation in Tanioka et al. (2006).

Summary of the Indian Ocean study
Summing up, TsunAWI is able to reproduce observational results on several scales of wave propagation.In our approach with unstructured meshes, both the large scale propagation throughout the Indian Ocean and the inundation in a selected area with high resolution can be combined within a single model run.However, the results strongly depend on the quality of the source model.

Application in the German-Indonesian Tsunami Early Warning System (GITEWS)
The German-Indonesian Tsunami Early Warning System was founded after the devastating Indian Ocean tsunami, in 2004, as a joint project of German research institutes as well as Indonesian institutes and authorities.In 2005, Indonesia and Germany agreed in a joint declaration to develop a warning system coordinated by the UNESCO Intergovernmental Oceanographic Commission.The system was inaugurated in 2008, evaluated in 2010 by an international commission, and handed over to the Indonesian government in 2011.In the ongoing project for training, education and consulting for the tsunami early warning system (PROTECTS), the Indonesian staff are trained in the use and maintenance of all components.

Tsunami scenario repository for GITEWS
The Indonesian warning system meets the challenge of near field warning with extremely short warning time.As tsunamigenic earthquakes originate close to the shore, the time between a seismic event and the issue of a warning is limited to just a few minutes.The early warning system therefore relies on a repository of tsunami scenarios precomputed with the tsunami model TsunAWI.To speed up the processing in a warning situation, an index database is extracted offline.It contains data products for all scenarios like virtual sensor data and shapefiles with estimated maximum amplitude and arrival times at coastal forecast zones, and isochrones.

Warning system architecture
Data of various sensors, such as the seismic stations analysed by SeisComP (Seismological Communication Processor), global positioning system (GPS) sensors as well as data from tide gauges and buoys, are collected via the Tsunami Service Bus (TSB, Fleischer et al., 2010), and distributed to the Decision Support System (DSS, Steinmetz et al., 2010), which triggers the Simulation System (SIM) with the available sensor data.On this basis, the SIM delivers a set of best matching scenarios.The DSS performs a worst case aggregation over these scenarios and visualizes expected maximum amplitudes and arrival times for the Indonesian coasts.The system is installed at the warning centre in Jakarta, Indonesia, and assists the chief officers on duty to assess the potential tsunami risk and to disseminate warnings to governmental institutions, local disaster management, action forces, and media.

Tsunami scenario repository
In the beginning of the GITEWS project, emphasis was put on getting the warning system operational as fast as possible and enhancing the running system step by step.For the tsunami scenario repository (TSR), this approach was followed and as soon as good bathymetry data and a rupture model (RuptGen by Babeyko et al., 2010) had become available, first scenarios were computed with a first version of TsunAWI based on the linear SWE and with linear viscosity.As computation and quality control took about half a day per scenario, the repository first covered only a smaller number of epicentres and magnitudes M w = 8.0, 8.5, 9.0.In the daily routine at the warning center, it soon turned out that scenarios for lower magnitudes M w ∈ [7.2, 8.0] were most urgently needed, because smaller earthquakes occur far more frequently and their tsunami potential is much harder to assess.Therefore, the TSR was extended to contain lower magnitudes M w .
In early 2011, all scenarios in the TSR were replaced in one step by new ones that were calculated with the enhanced TsunAWI code including non-linear advection, Smagorinsky viscosity, improved inundation scheme, and with initial conditions derived from an extended version of RuptGen (Rupture Generator) offering also ruptures closer to the coast.The TSR now contains 3470 scenarios for prototypic ruptures with magnitudes evenly spaced in the range of M w = 7.2, 7.4, . . ., 9.0 on 528 different epicentres, see Fig. 13.
These scenarios are used in the warning situation as well as for risk analysis and hazard maps, as long as no simulations with higher resolution are available.For earthquakes with epicentres outside the area covered by the TSR, e.g. in the Banda Sea, a tsunami simulation is performed within 1 min computation time by the model EasyWave (Babeyko, 2012).EasyWave discretises the linear SWE in deep water on a regular coarse grid and estimates the amplitude at the coast by Green's law.Currently, we are preparing TsunAWI scenarios for epicentres in the eastern Sunda Arc.
For the TSR, TsunAWI operates on a grid with a resolution of 14 km in deep sea, 150 m at the coast and down to 50 m in regions of special interest, e.g.densely populated areas and at tide gauges.The grid contains 2.3 million nodes, model time step is 1s, and model time 3 h.One scenario with time steps written every minute occupies a disk space of 2 GB, which is reduced in a post processing step to 1.1 GB by deleting all grid nodes on land that are never inundated in any scenario.The total disk space required by the TSR adds up to 4 TB.The time series are kept to easily extract mareograms at new tide gauge or buoy locations.The dense net of scenarios allows to evaluate the effect of small changes in epicentre and magnitude as shown in Fig. 14 for artificial mareograms for the harbour of Padang, Sumatra.The logarithmic energy scale of the magnitude is clearly illustrated by the upper mareogram.The lower graph shows more complex features due to varying epicentres on a line perpendicular to the trench.With the magnitude fixed to M w = 8.0, the tsunamis originating closest to the coast have the lowest impact.On the one hand, this is due to the initial conditions generated by RuptGen, which assumes a source depth of 0 km at the trench increasing up to 100 km under the main islands.A second factor determining the strength of the tsunami is the height of the water column at the origin, which is small in coastal regions.
Epicentres to the southwest off the coast yield increasing wave heights in Padang, until the Mentawai islands are reached.Farther west, the corresponding tsunami is no longer trapped in reflections between Sumatra and the islands, relieving the situation for Padang.

SIM selection algorithm
The selection algorithm implemented in the Simulation Module (SIM) combines the available sensor data to acquire a set of best matching scenarios to an earthquake event.By basing the selection on different sensor types, uncertainties resulting from inaccurate measurements and errors in the tsunami model are reduced.A more detailed overview of the SIM is given in Behrens et al. (2010).
Originally, the so-called matching values generated for each sensor type for a scenario were accumulated in an overall weighted sum.However, the experiences with real sensor data in the GITEWS project showed that each sensor group has to be regarded separately with its characteristics in mind.The weighted sum has therefore been replaced by a stepped approach starting with seismic data delivering the most robust values, followed by GPS data.
The algorithm takes the tectonic structure of the Sunda Arc into account.For the preselection of scenarios based on seismic data, an elliptic area around the measured epicentre is constructed, see Fig. 15.The dimension of the ellipse depends on the measured magnitude M w with the long ellipse axis given by r L = max 10 0.5[M w +0.3]−1.8km, 180 km . (9) The lower value of 180 km ensures that at least one scenario is covered by the ellipse.The long axis is constructed parallel to the trench.
The second important sensor class in GITEWS are GPS dislocation vectors.The GPS measurements are fast and reliable, thus allowing for a better estimate of the tsunami hazard in the first few minutes after the earthquake.For each sensor and each scenario that remains after preselection by the seismic data, the lengths of the dislocation vector from the measurement and from the corresponding scenario are compared.Scenarios for which a defined minimum number of similar GPS values has been reached are taken into account for the final set of best fitting scenarios.This allows to narrow down the preselected area and to correct overestimated magnitudes.With more GPS sensors available, it will become possible to correct underestimated magnitudes and enlarge the seismic preselection.
Once several scenarios have been chosen for one epicentre, only the one with the largest magnitude within the uncertainty range [M w − 0.5, M w + 0.3] is kept in the list, because the scenarios with lower magnitudes will have no relevance for the worst case aggregation in DSS processing.Hence, disregarding them reduces the amount of data to be processed without changing the result.

Scenarios for Indonesia as Regional Tsunami
Service Provider (RTSP) In 2011, the warning centre at BMKG, Jakarta, became a regional tsunami service provider (RTSP) for the Indian Ocean together with the warning centres in India and Australia.Therefore, a second set of scenarios was required, one which covers the whole Indian Ocean with a model time of 24 h and provides a predefined set of data products like arrival times and maximum amplitudes in the forecast zones along the coast of the Indian Ocean.
Though the resolution was reduced to a minimum of 150 m along the coast and 25 km in the deep ocean, the grid consists of three times as many nodes as the regional grid.Together with longer model time, a full scenario would end up with a file size of 40 GB.Though the number of scenarios, 1900 in total, is smaller, because only those posing a threat outside the area already covered by the GITEWS TSR are required, the amount of disk space became impracticable.Therefore,  8.8,8.6,8.4,8.2,8.0,7.8,7.6,7.4only timelines at nodes with a water depth of at least 50 m were kept to provide smooth isochrones, and after post processing, the timelines were completely deleted.Global fields of maximum amplitude and arrival times are permanently stored, together with the all preliminaries to recalculate the scenarios, if necessary.

Inundation studies for risk assessment
The focus of this section is to highlight the importance of the bottom roughness parameter.Simulations of tsunamis with identical initial conditions were carried out in a mesh with a resolution down to 5 m.The study area is Mataram, Lombok, Indonesia, and as a tsunami source, the Okada parameters for an assumed earthquake with M w = 8.5, epicentre at 8.15 • S, 116.15 • E, in the back-arc region north of Lombok was provided to us by the Australia-Indonesia Facility for Disaster Reduction (AIFDR).
The topography data are based on the Intermap data set with 5 m horizontal resolution, provided to us by the DLR.Two versions, a digital surface map (DSM) and a digital terrain map (DTM), were used.The DSM is based on a firstreflection data set and contains elevations of vegetation and buildings.In the DTM, all these features have been removed.
The scenarios generated for the tsunami database are based on a coarser mesh with a resolution down to 50 m in priority regions.This resolution is a compromise between the significance of the simulation and the computational efforts.In particular, evacuation maps are based on data sets with much higher resolution, containing streets and other infrastructure.The corresponding model set-up should contain a comparable mesh, and in GITEWS, high resolution simulations were performed for large cities like Bengkulu, Cilacap, or Denpasar with MIKE 21 FM (Flexible Mesh) by the German Research Center Geesthacht (GKSS) (Gayer et al., 2010).In these simulations, the Manning roughness parameter varies spatially according to detailed roughness maps.Details on this approach and the model can be found in Gayer et al. (2010) and citations therein.
The simulated inundation in the TSR is used for the other coastal regions.Initially, the wetting and drying scheme was only included in TsunAWI to prevent artificial wave reflections at ocean basin walls.However, it soon turned out that the implemented scheme was also well suited for realistic inundation simulation.
Figure 17 displays the flow depth based on the DSM topography, Fig. 16 for the DTM with different values of Manning's n.Additionally, in case of a tsunami, buildings and infrastructure will partly collapse and result in a very  formed with the MPI-parallel research branch of TsunAWI, while the operational branch is kept as lean as possible.For the latter, with the code optimisation so far achieved and the advances in hardware near real time computing comes into reach -without reducing model physics.
First real time models, e.g., EasyWave (Babeyko, 2012) and RIFT (Wang et al., 2012a), allow to react in an early warning situation on the seismic information in addition to the epicenter and magnitude.Fast momentum tensor derivation and GPS measurements can narrow down the potential tsunami source in just a few minutes after the earthquake.To save computing time, EasyWave and RIFT employ reduced model physics based on linear SWE in deep water with the extrapolation of the maximum amplitude to the shore by Green's law.For the tsunami inundation hazard in the far 690 field, Tang et al. (2009) used the example of Hawaii to con-duct a sensitivity study with respect to source parameters, bathymetry data, and grid resolution.
The computational efficiency of TsunAWI is very high regarding model physics and resolution, but the code still has 695 to be adapted to the early warning situation.One strategy will be to coarsen grid resolution by evaluating the local CFL numbers reached in model runs in the given geographical setting.Another option is to employ general purpose graphics processing units (GPGPU) that make the massive computational power of a modern graphics accelerator available for general computing.But this option requires the maintainance of two code branches, because optimisation strategies for cache efficiency on typical CPUs and for the limited registers on GPGPUs differ extremly.

6 Conclusions
During seven years of development, TsunAWI matured from a reseach model to an operational code with a framework of pre-and postprocessing routines.In particular, the implementation of the non-linear advection and the wetting and 710 drying required the careful comparison of different schemes to achieve an algorithm that is stable without being too diffusive.The model has since proven its robustness, correctness, and efficiency in several hindcasts of tsunami events world wide and in more than 5000 scenario calculations for 715 the Indonesian Seas and the Indian Ocean.However, the experience of the user in model setup and evaluation remains an important requirement for good results in tsunami modelling.
Sidorenko, Widodo Pranowo, and Chaeroni.We thank the colleagues at DLR, Matthias Mück and Ralph Kiefl, for providing us with Intermap topography for Materam and with inundation areas estimated from satellite images, respectively.Widjo Kongko  heterogeneous fluid flow, which may be only poorly approximated by a bottom friction parameterization.The result based on DSM contains small scale features, however, since the destruction of buildings and vegetation is not taken into account, the extent of the inundation may be underestimated.Gayer et al. (2010) suggests the use of DTM data with detailed roughness maps.According to Gayer et al. (2010), the right panel corresponds to the situation of an urban area with partly collapsing buildings whereas the left panel corresponds to the situation of the whole area covered by coarse sand.
Given the vast differences displayed in Fig. 16, the results are meant to raise the awareness to the fact that the simulated inundation strongly depends on the friction parameters.The quality of topography data plays an even bigger role and the best suited data set and model parameterization must be chosen before hazard maps are produced.Furthermore, the derivation of the shallow water equations (SWE) (1) and ( 2) is based on the assumption of a small ratio of water height to wave length h/λ 1.If the wave is determined by very fine scale topographic features, the validity of the SWE becomes questionable.As long as the focus of such model studies lies on the inundated area alone, the limitations of the SWE might still be acceptable.In addition to the work presented here, it would be desirable to assess the limitations of TsunAWI compared to models with e.g. three dimensional simulation, non-hydrostatic flows, wave breaking, or sediment transport.

Outlook
TsunAWI is under constant development and our current research focusses on two modular additions, one for the nonhydrostatic pressure correction, the second for tidal forcing.Both aim to improve the tsunami hazard assessment.
The non-hydrostatic correction follows the approach of Walters (2005) and can be activated on demand, when coastline geometry or steep bathymetric gradients necessitate it for realistic simulation results.
Tides are usually neglected in tsunami modelling, because a typical tsunami wave has a much shorter wave length.However, simulation with TsunAWI (Androsov, 2010;Androsov et al., 2011) showed that the inclusion of tides can have a significant effect on the amplitude and phase of tsunami waves in regions with strong tidal currents or shallow bathymetry.For these experiments, a small simulation area was chosen and the tidal forcing was imposed on the boundary conditions.To simplify the set-up for simulations in arbitrary regions, TsunAWI is now extended to a global tidal model, and the region of interest is refined in the unstructured grid.
To make these extensions feasible, emphasis is put on an efficient implementation on modern computer architectures.The compute-intensive extensions named above are performed with the MPI-parallel research branch of TsunAWI, while the operational branch is kept as lean as possible.For the latter, with the code optimisation so far achieved and the advances in hardware, near-real-time computing comes into reach -without reducing model physics.
First real-time models, e.g.EasyWave (Babeyko, 2012) and RIFT (Real-time Inundation Forecasting of Tsunamis) (Wang et al., 2012a), allow to react in an early warning situation on the seismic information in addition to the epicentre and magnitude.Fast momentum tensor derivation and GPS measurements can narrow down the potential tsunami source in just a few minutes after the earthquake.To save computing time, EasyWave and RIFT employ reduced model physics based on linear SWE in deep water with the extrapolation of the maximum amplitude to the shore by Green's law.For the tsunami inundation hazard in the far field, Tang et al. (2009) used the example of Hawaii to conduct a sensitivity study with respect to source parameters, bathymetry data, and grid resolution.
The computational efficiency of TsunAWI is very high regarding model physics and resolution, but the code still has to be adapted to the early warning situation.One strategy will be to coarsen grid resolution by evaluating the local CFL numbers reached in model runs in the given geographical setting.Another option is to employ general purpose graphics processing units (GPGPU) that make the massive computational power of a modern graphics accelerator available for general computing.But this option requires the maintainance of two code branches, because optimisation strategies for cache efficiency on typical CPUs and for the limited registers on GPGPUs differ extremly.

Conclusions
During seven years of development, TsunAWI matured from a research model to an operational code with a framework of pre-and postprocessing routines.In particular, the implementation of the non-linear advection and the wetting and drying required the careful comparison of different schemes to achieve an algorithm that is stable without being too diffusive.The model has since proven its robustness, correctness, and efficiency in several hindcasts of tsunami events worldwide and in more than 5000 scenario calculations for the Indonesian seas and the Indian Ocean.However, the experience of the user in model set-up and evaluation remains an important requirement for good results in tsunami modelling.

N.
Rakowsky et al.: Tsunami modelling with TsunAWI all triangles fulfill the criterion x ≤ min c CFL gh, c g h ∇h .

Fig. 1 .
Fig.1.Three different orderings of variables in the regional Indonesian grid, from dark blue for the first node to dark red for the last one (2.3 million).On the left for the whole domain, on the right a zoom at Bali region reduced to 16 colours emulating the distribution on 16 OpenMP threads.

Fig. 3 .
Fig. 3. Arrival time for the same scenario as in Fig. 2 but with a relative threshold of 10 % of the local maximum amplitude.

Fig. 5 .
Fig. 5. Tide gauge locations and mesh density in the model domain.The green areas are land nodes which are initially dry.The resolution ranges from 15 km in the deep ocean to 40 m in Aceh, Sumatra.

Fig. 8 .
Fig.8.Hovmöller diagram for the modelled sea surface elevation of the Jason 1 track (Fig.6).The model results corresponding to the satellite observations (Fig.7) are interpolated in space and time.

Fig. 9 .
Fig. 9. Positions of field measurements and flowdepth obtained with Manning number n = 0.03 and original Tanioka source.The red line indicates the inundated area as it was determined from satellite images by DLR.

Fig. 10 .
Fig. 10.Flow depth comparison in locations marked in Fig. 9 for different roughness parameters n and the modified source (strike angle 290/320).

Fig. 11 .
Fig. 11.Flow depth comparison in locations marked in Fig. 9 for strike angles of plates A/C set to 340/340 and 290/320.

Fig. 12 .
Fig. 12. Tide gauge (TG) records and corresponding model results in locations displayed in Fig. 5.

Fig. 13 .
Fig. 13.Tsunami scenario repository for GITEWS.Scenarios for up to ten different magnitudes M w = 7.2, 7.4, . . ., 9.0 are calculated at each of the 528 positions marked with red dots, resulting in 3470 scenarios in total.

Fig. 14 .Fig. 15 .
Fig. 14.Influence of the magnitude M w and of the source location on the strength of a tsunami: mareograms at Padang Harbor for TSRscenarios originating at the same epicentre with different magnitude (above) and with equal magnitude, but different origin (below).

Fig. 16 .
Fig. 16.Flowdepth in the model study for Mataram/Lombok (background © OpenStreetMap contributors) based on Intermap DTM with three different Manning parameters n=0.02, 0.04, 0.06, 0.08 from the left to the right panel.

Fig. 16 .
Fig. 16.Flowdepth in the model study for Mataram/Lombok (background © OpenStreetMap contributors) based on Intermap DTM with four different Manning parameters n = 0.02, 0.04, 0.06, 0.08 from the left to the right panel.

Table 2 .
Satellite missions Jason 1 (J1) and Topex/Poseidon (T/P).The last columns quantify the RMS error obtained in the TsunAWI scenarios with different strike angles in subfaults A and C (Fig.4) during the approx.10 min of the satellite overpasses.