The application of open-source software ANUGA installed on a BlueJ/Q computer for a tsunami modelling in costal areas of city of Trieste using high resolution laser-scanning and multibeam data

STERZAI, Paoloa, COREN, Francoa and CREATI, Nicolaa

a Istituto Nazionale di Oceanografia e di Geofisica Sperimentale, Italy, e-mail:

Abstract — Due to the characteristics of seismo-tectonic is an event of a tsunami in the northern Adriatic considered as unlikely. However many historical sources report that in 1511 a tsunami induced by an earthquake occurred devastating the city. The problems of the city in case of heavy rains about the coastal flooding is well known especially in the case of west - winds. This paper presents a dynamic hydraulic modeling of the city of Trieste in the event of a tsunami during a contemporary storm surge along the coastal areas using the open-source software ANUGA installed on a machine BlueJ/Q. The hydraulic modeling used for this purpose calculates the solutions from the fluid motion equation through the Navier-Stokes equations using some mathematical simplifications in the case of wave equations due to shallow water equations, as the fluid incompressibility, uniform vertical velocity distribution, uniform distribution of the hydrostatic pressure and a gentle slope of the seabed. A digital elevation model of the terrain has been used along the coast obtained by laser-scanning fused with a multibeam batimetry for the onshore data, with a final resolution of about 1 point per square meter. The data were corrected for the average height of the seas defined by the geoid ADBVE 2008. Due to the large amount of data using a high resolution digital terrain as well as batymetric model, the modelling computations have been done using the super-computer BlueJ/Q, which optimized the parallelization of the code. The hydraulic modeling high resolution required in urban areas which takes account of buildings, road structures, channels etc.. requires a computing capacity which is not possible to deal with PC seen the magnitude of the mesh used and the complexity considered to obtain reliable temporal evolution. Considering occurrence of an event of a tsunami in correspondence to an event of a flood wave 2 meters due to storm surges, then the worst possible case, by entering into a GIS all simulations it was possible to outline the areas of highest risk, indicating the areas that would require surgery in order to prevent problems of civil protection.

Keywords — risk management; hydraulic modeling, civil protection

1 Introduction

After the event of the tsunami in Southeast Asia of 26/12/2004 , the scientific community has begun to consider the event of these phenomena even in areas where previously were considered unlikely . One of these areas is considered the northern Adriatic (Lorito et al. 2008), in which due to the characteristics of the seismo- tectonic zone, it is considered an unlikely event. The main sismogenetic faults in this area are present on the Dinaric Alps and Julian Alps, located in areas far from being considered potential tsunami generators. An earthquake in these areas could cause tsunamis (Tinti et. al. 2004). From the year 1000 A.C. to now, five tsunamis have occurred in the Adriatic (among the 32 registered in Italy), cataloged by the National Seismological Service. The most important one took place in 1511 in the North Adriatic and caused a rise in the sea at Trieste (Tinti et. al. 2004); even though different authors seems to deny such an event (Camassi et . al. 2011). The catastrophic earthquake that broke out between Friuli and Slovenia in the afternoon of March 26, 1511 shocked all our neighboring regions to the Marche, (CPTI, 2004). This earthquake caused moreover one of the most tragic tsunami ever recorded in the Adriatic , which hit the towns of Venice and Trieste. The victims of the event overall , it is estimated, were approximately 10000 and the estimated sea level rise was of the order of one meter. In addition to the tsunami event one more frequent problem of the town of Trieste is the storm surge; this occurs in case of heavy rains especially in the case of south winds (scirocco), (Marone E., 2013), that can cause a water surge in the Gulf of Trieste of about 2 meters.

This work is based on modeling of a hypothetical tsunami similar to that occurred in 1511 , generated by an earthquake located in the position of 46.200 of latitude, 13.430 of longitude, see (CPTI, 2004) , that occurred contemporarily to an event of storm surge, that would be the worst case scenario. In order to optimize the modeling we have chosen to perform the calculations using a high-resolution terrain model of the city of Trieste obtained from laser scanning for the terrestrial part, while the seas part has been derived from multibeam.

2 The Case Study

The direct acquisition of a high density and accurate 3D point cloud has made lidar systems the preferred technology for the collection of topographic data in urban environment, (Wehr, 1999; Coren F. et al., 2004). In this study, Optech ALTM 3100 was used for the acquisition of the LiDAR data, with the flight performed on 30th April 2004. The study area was measured from an altitude of 1500 m, with a sampling density of about 4 points per square meter, and the radiometric resolution, scan frequency and scan width were 12bits, 70Hz and ±25°, respectively. The dataset has been processed using PosPac software for the trajectory computation. The final data were obtained using Optech DashMap software while TerraScan (Terrasolid Corporation) was used for data classification in order to produce a good map of the city buildings of the city, (Axelsson P., 1999). The lidar data have been integrated with batymetric coastal data collected with a multibeam system. A multibeam echosounder is a device typically used by hydrographic surveyors to determine the depth of water and the nature of the seabed. We used the Reason 6000 that is operated transmitting a broad acoustic fan shaped pulse from a specially designed transducer across the full swath acrosstrack with a narrow alongtrack then forming multiple receive beams (beamforming) that is much narrower in the acrosstrack (around 1 degree depending on the system). From this narrow beam a two way travel time of the acoustic pulse is then established utilizing a bottom detection algorithm. If the speed of sound in water is known for the full water column profile, the depth and position of the return signal can be determined. The usefull dataset covers the Trieste city costal-area for the first 500 meters from the shoreline.

A good data fusion between aerial laser-scan and multibeam data is usually hard to obtain due to the data gap caused because the laser-scan data cannot penetrate a water surface as well as the multibeam data cannot reach zones along the shoreline. To overcome this problem it was decided to fill the data gap using a kriging interpolator that allowed to obtain a smooth transition area between both data avoiding oscillations usually obtained using other interpolating techniques due to the triangulation.

With the aim to extend the bathymetric data of the Gulf of Trieste in order to obtain a total bathymetric cover until 2 km from the costal-area, about 50 seismic sections were digitized. The depth were obtained from the first reflection (estimated being the sea-bottom), multiplying the time for reflection with an average speed propagation of the seismic wave in the sea, by applying a correction due to the tide. The chosen grid, which well represents the spatial bathymetric resolution achieved by digitizing seismic sections, was 200 meters. Figure 1 shows the obtained bathymetric data fusion, while figure 2 represents the complete dataset.

Figure 1: The laser-scan elevation data fused with multibeam bathymetric map used in hydraulic modelling. With the pink color is identified the low-resolution bathymetry calculated from seismic lines digitalization. White color indicate no-data
Figure 2: Hypothetical location of the earthquake occured in 1512, the red point is the location of the historical earthquake, the square red-area is the tsunami investigated area where the computations has been done

Using the software Anuga is possible to import datasets of different resolutions, giving different weights (or priorities) to the individual areas. In this way was possible to reduce the computation time because the generated mesh is not uniformly spaced along the entire dataset but are dependent on the chosen resolution . The chosen resolutions was 1 meter in the coastal downtown area (where the water can reach) as well as in the first 500 meters from the coast, while in other areas a 200 meter mesh was considered as sufficient. This was done using the mesh generator of ANUGA defining the geometry of the problem in an interactive way. An almost common event in the Gulf of Trieste is the so-called high water, which can take place mainly in autumn and more precisely in November. The causes of this phenomenon are almost meteorological: a low atmospheric pressure and persistent south winds on the Adriatic basin can cause a rise in sea level in its northern part up to 70 cm. The same meteorological can generate another phenomenon in the Adriatic called the longitudinal sessa , which is a fluctuation in the level with a period of 21.5 h, and a width that can reach 50 cm in the Trieste area. Both phenomenon can also produce an overall sea level increase of 120 cm, E. Marone (2013). In the case that this events coincides with a high astronomical tide, a high water resulting can reach almost 2 meters. So far the maximum high water level recorded in Trieste was at 11/26/1969 with 193 cm up to the local mareographic zero called IGM. We adopted the value of 2 meters for the simulations. In order to simplify the problem it is necessary to make the following simplifications. The differential equations of a two-dimensional flow are considered according to the methodology proposed by Godunov: considering only those solutions to the first order of accuracy, however, apply only in the case of shallow-water, where the horizontal dimensions are considered as much greater than the column water. Finally the water surface is considered not disturbed in a body of water vertically homogeneous and isotropic.

The software we used is ANUGA that records the evolution of the water depth and the horizontal time in the interior of each cell in function of time by solving the equations of shallow water starting from a finite element approach, for the mathematical background see Nielsen, O. (2007) . The software is written in the object-oriented language Python, allowing to interface easily the Anuga library functions. In order to increase the code efficiency as well as linking directly with the structures of Python numpy some more complex components are written in the C language. The core of ANUGA is the hydrodynamic module, called shallow-water, which is based on the finite volume method using triangular cells meshes to solve the shallow water equations. The finite volume method is a very robust and flexible numerical method technique; however it is certainly not the fastest. If the geometry is simple, a finite difference method can be more efficient.

The frictional resistence is applied using the Manning formula, although the software was not yet tested in relation to the roughness of the sea-bottom. The Manning’s coefficient is set at 0.1 at each point of the mesh, (M. H. Dao and P. Tkalich, 2007). In order to make the problem more simple, Anuga uses a computational approach applying the so-called Dirichlet boundary which is often used in hydro- dynamics, the no-slip condition for viscous fluids states that at a solid boundary have zero velocity relatively to the boundary.

The calculations have be done considering the tide at 2 meters, considering the mean-sea level of the sea referred to the ADBVE 2006 geoid, (Sterzai et al., 2008). Regarding the tsunami wave; the modeling has been performed a modeling considering an option called “fixed- wave" that consider a wave applied to a chosen boundary. Therefore we do not modelled the tsunami source parameters in terms of an hypothetical sesimic event but we used planar wave approximation; this assumption is justified by the fact that the source of the tsunami is located onshore and far away from the basin (almost 100 km; see figure 2). This also imply that we cannot model the arrival time of the tsunami wave cosidering that we didn’t effectively modelled the seismi source and therefore the earthqwuacke that generate the tsunami. In addition the modeling does not take into account other external factors that could bring the system out of equilibrium as rains, water production and external stress due to wind and atmospheric pressure gradients.

2.1 Data analysis

As can be seen from the Figure 5, the Trieste downtown is affected the most compared to the other harbour area in the bay. The reason is the water depth is almost shallow there and that part is an enclosed area which traps and amplifies the tsunami energy. By using these simulation results and topographical data, inundation maps for Trieste area have been prepared in a GIS environment.

From the data of temporal evolution of the tsunami was possible to calculate the maximum flood wave of the tsunami, the absolute moment and the speed of impact.

The momentum can be a good variable to point out the tsunami damage risk in a area, H. G. Loomis (2002) . It contains in the largest positive wave crest is a measure of its damaging potential and might be used as a measurement of tsunami magnitude. The momentum flux represents the most suitable damage indicator and can be computed from the equation:

MT(t) =
zFT(t)dz, (1)

where F(t) is the impact force computed from the integral over the whole depth z and n and h are the local amplitude and undisturbed water depth at the pile. It can be noted that the absolute momentum is a variable not depending from the wave period, see the figure 3.

Figure 3: Absolute momentum of Trieste tsunami modelling.

Since it is considered to an earthquake place in the north-west, from where the wave front could move to the city of Trieste, it is noted that the areas of higher values of the absolute moment can be found on the coast near the city downtown, area where cruise ships dock, see figure 3. Particularly vulnerable seems to be also the areas behind the dams and the port area. For the protection of the dams is to be specially protected area to the east.

2.2 Risk Map for Trieste city

The simulation duration is taken as 30 min of real-time computation with the time step of 0.1 sec. The software computes the propagation of the tsunami wave at every 0.1 sec and gives the outputs in every 60 sec. In order to be able to compare the maximum positive and maximum negative wave amplitudes near Trieste, artificial computation gauges are placed to specific locations in the study domain. In the simulation all tsunami parameters in the near shore area of Trieste town are computed. Initial wave is selected as a site specific tsunami, and its characteristics represent a tsunami characteristics based on the tsunami generated in the offshore area in the direction from the earthquake direction to Trieste. After a 30 minutes long simulation, the software gives the propagation of the tsunami, run-up in addition to sea states at specified time steps and maximum amplitude during the simulation duration. According to the simulation results, the maximum elevation of the tsunami to Trieste from simulation source is with the wave height of 2.0 m, figure 4.

Figure 4: Maximum innundation map of Trieste downtown. The colorscale is in meters. Here is possible to note the inundation map of the city. The water drowns most of the old town and the coastal area that by the way is a very populated area.

Therefore, before defining the risk, first hazard and vulnerability of the case should be specified clearly. Mitigation plan for a disaster management includes reconstruction and preparedness stages for the disaster event. Since we know the extent of the hazard in terms of source of the waves caused by earthquake, affected zones can be determined by using GIS. The inundation map formed in GIS environment is used for this purpose. The intersections of the building and road layers and inundation. Map layer have been used to detect the buildings affected by tsunami waves. Besides this, zonal statistics tool has been used to extract the elevation information from DEM. Finally, by considering the inundation possibility and elevation parameter, a risk map has been prepared for Trieste after the tsunami, see the figure 5.

Figure 5: GIS of the risk map for Trieste city with identification of the building interested by inundation.

3 Conclusion

In conclusion, a tsunami hydraulic modelling has been computed the inundation map for a hypothetical event caused by a earthquake that has occurred in 1511 in the area of the town of Trieste superimposed of a surge storm sea level rise of 2 meters. Using a high- resolution terrain model obtained by laser-scanning for the terrestrial and marine multibeam for the coastal data and using a machine with high computing power as is the Blue J / Q was possible to model the possible consequences of a tsunami even of moderate intensity. In addition, through the measurement of the absolute momentum that is associated with the intensity of a tsunami was possible to point out the most dangerous zones: these can be identified. The future developments will be focused trying to extend the study in other areas where high-resolution digital terrain models are available.


Axelsson, P. (1999): Processing of laser scanner data—algorithms and applications, ISPRS Journal of Photogrammetry & Remote Sensing 54, 138–147

Bellotti, G.and Cecioni C. (2010): Propagation of Tsunamis over Large Areas using COMSOL, Excerpt from the Proceedings of the COMSOL Conference 2010 Paris

Camassi, R. ; Caracciolo, C.H., Castelli, V., Slejko d. (2011) •The 1511 Eastern Alps earthquakes: a critical update and comparison of existing macroseismic datasets. J Seismol (2011) 15:191–213 DOI 10.1007/s10950-010-9220-9

Catalogo Parametrico dei Terremoti Italiani, versione 2004 (CPTI04), INGV, Bologna.

Coren, F. and Sterzai, P., (2006): Radiometric correction in laser scanning. International Journal of Remote Sensing, 27 (15-16), 3097-3104.

Dao, M.H., Tkalich P. (2007): Tsunami propagation modelling - a sensitivity study, Nat. Hazards Earth Syst. Sci., 7, 741-754

Evaluating Tsunami Impact Metrics, Tsunami Pilot Study Working Group Seaside, Oregon Tsunami Pilot Study

Marone, E., et al. (2013): Harmonic tidal analysis methods on time and frequency domains: similarities and differences for the Gulf of Trieste, Italy, and Paranagu Bay, Brazil, Vol. 54, n.2, June 2013 pp. 183-204

Loomis, H.G. (2002): Science of Tsunami Hazards, Vol. 24, No. 5, page 317 (2006)

Loomis, H.G. (2002): The Momentum of a Tsunami, Science of Tsunami Hazards, Vol.20, No.1

Lorito, S. et al. (2008): Earthquake-generated tsunamis in the Mediterranean Sea: Scenarios of potential threats to Southern Italy, Journal of Geophysical Research, Vol. 113, B01301, doi:10.1029/2007JB004943

Salap, S. et al.: Tsunami Risk Analysis and Disaster Management by Using GIS: A Case Study in Southwest Turkey, Goecek Bay Area

Nielsen, O. (2007): ANUGA v1.0 User Manual. Geoscience Australia and the Australian National University, pp. 71-72

Sterzai , P. et al.(2008): An Improved Geoid in North Eastern Italy, Observing our Changing Earth , International Association of Geodesy Symposia , Vol. 133, 2008, pp 427-430

Tinti S., Maramai A. and Graziani L. (2004). The New Catalogue of Italian Tsunamis. Natural Hazards, Vol.33, n.3, 437-463.

Wehr, A. and Lohr, U. (1999): Airborne laser scanning—an introduction and overview, ISPRS Journal of Photogrammetry & Remote Sensing 54, 68–82


Sterzai, P., Coren F. and Creati, N. (2015): The application of open-source software ANUGA installed on a BlueJ/Q computer for a tsunami modelling in costal areas of city of Trieste using high resolution laser-scanning and multibeam data. In: Planet@Risk, 3(1): 121-125, Davos: Global Risk Forum GRF Davos.