Modelling with Volna-OP2: Towards tsunami threat reduction

Accurate and efficient tsunami modelling is essential for providing tsunami forecasts and hazard assessments. Volna-OP2 is a finite volume solver of the nonlinear shallow water equations and its capabilities of producing both faster than real time ensembles and high resolution inundation studies are presented here. The code is massively parallelised and can utilise various high performance computing architectures. When an earthquake is detected there is always some uncertainty on the source parameters. Generating a faster than real time ensemble for maximum wave heights which captures this uncertainty would be of great benefit to tsunami warning centres. The 2003 Boumerdes earthquake (Algeria) acts as a test case for showing Volna-OP2's ability at rapidly forecasting regional maximum wave heights. Drawing on various earthquake sources proposed in the literature and scaling the magnitudes to mimic uncertainty on the source, 20 separate earthquake realisations are simulated for 4 hours real time in 97s on two Nvidia V100 GPUs. Further a reduced ensemble of the Lisbon 1755 tsunami with an emphasis on the effects to the Irish coastline is presented. Where again various earthquake sources have been drawn from the literature and simulated on a regional scale. Finally, a pilot study which builds upon the reduced ensemble results investigates the inundation of a Lisbon tsunami on key sections of the Irish coastline. The results of this pilot study highlight that the inundation is constrained to low-lying areas with maximum run-up heights of $\approx 3.4m$ being found.


Introduction
Tsunamis are infrequent events that have the capacity of being extremely destructive. Owing to their rarity the tsunami community relies heavily on mathematical and computational modelling to provide accurate forecasts and risk assessments. Tsunamis can travel at extremely high speeds in deep water, thus arrival times at the coastline are often on the O(mins). Tsunami Warning Centres (TWC) are therefore under severe time constraints when deciding the level of warning and areas at risk from an event. To aid in this decision process, TWCs currently deploy a variety of tools which produce regional levels of warnings.
However, it is well documented that the local bathymetry can have a substantial effect on the local levels of maximum wave heights, inundation and induced currents. To gain an understanding of how the local bathymetry affects an incoming tsunami, prior high resolution risk assessments are often carried out. These studies require extensive computational resources and are time consuming. Epistemic uncertainty involved in the modelling and source parameters has lead to the development of probabilistic tsunami hazard assessments (PTHA)[GP06, GBB + 17]. One of the key limiting factors for both TWCs and risk assessment studies is the inherent tsunami modelling software. Extremely efficient tsunami codes, which utlitize modern computational techniques, have been developed to reduce the run times involved.
The software package Volna-OP2 is a tsunami code which has been optimised to run on various high performance architectures, CPUs, Xeon-Phis and GPUs. To date, Volna-OP2 has been used by various research groups around the world. In particular its efficiency has been leveraged in conjunction with statistical studies to carry out computationally expensive tasks like sensitivity studies and building statistical surrogates [SGGD17, GSD + 18, GHG20]. The code solves the depth-averaged nonlinear shallow water equations, utilising modern reconstruction and shock capturing techniques. The original code was introduced in 2011 [DPD11] and has since undergone extensive parallelisation [RGG + 18]. Both the original and newly updated version have been verified against the traditional tsunami mitigation benchmark problems and an extensive error analysis of the newly parallelised version has been carried out [GKS + 20].
Volna-OP2 is presented as a candidate code which could be implemented into a TWC workflow. It is shown here to be capable of rapidly simulating tsunamis on a regional scale which can then be incorporated into 'on the fly' ensemble warnings. Further, to demonstrate the code's capability at performing high resolution inundation studies, a pilot mapping study of the Irish coastline has been carried out. This paper is thus split into three parts; in section(2) an outline of the numerical scheme and parallelisations used in Volna-OP2 are described. In section(3) a faster than real time tsunami ensemble is presented for the Mediterranean Sea and a reduced ensemble for the North East Atlantic Ocean. In section (4) a high resolution inundation mapping pilot study of the Irish coastline is shown. The main tsunami threat for the Irish coastline lies with a Lisbon type earthquake or a mid-Atlantic submarine landslide. This pilot study focuses on the Lisbon 1755 tsunami as a reasonable worst case scenario and builds upon the reduced ensemble results.

Numerical scheme and parallelisation
In the deep ocean tsunamis exhibit long wavelengths when compared to the depth. Thus, the simplest approach in modelling tsunami waves is to invoke the shallow water limit and to model the tsunami dynamics using the linear shallow water equations. These linear equations can be simply dealt with using computational techniques and have been shown to work well for modelling tsunamis across a deep ocean basin. However, once the ratio between the depth and wave height becomes comparable, nonlinear dynamics play an important role and the nonlinear shallow water equations are needed to represent the tsunami dynamics. Continuing on this trend, if one is to include dispersion or non-hydrostatic pressures one gets the Boussinesq equations. There are plenty of numerical solvers that exist which tackle these higher order Boussinesq equations. Undoubtedly, physical dispersion can play an important role but in most cases the nonlinear shallow water equations are sufficient.
Volna-OP2 is a finite-volume solver which solves the depth-average nonlinear shallow water equations, where H = (b + η) is the total water depth, b(x, y, t) is the time-dependent bathymetry, η(x, y, t) is the free surface deformation, u(u, v) is the depth averaged fluid velocity in the x and y horizontal directions, I is the identity matrix and g is the acceleration due to gravity. Volna-OP2 utilises unstructured triangular meshes with dependent variables assigned to the triangle's barycenters. In order to achieve 2 nd order accuracy in space, which is essential in modern methods, a MUSCL reconstruction scheme is implemented. For the temporal discretisation an option of using either a 2 nd or 3 rd order strong stability preserving Runge-Kutta scheme is present. Extensive details on the numerical scheme used in Volna-OP2 can be found in previous works [DPD11, GKS + 20].
When Volna was first developed in 2011 [DPD11] it was naturally a serial code. However, the code underwent a series of parallelisations culminating in the present Volna-OP2 version. The main parallelisation of the serial code occurred when it was built on the OP2 domain-specific language (DSL) for unstructured mesh computations [MGR + 12]. This DSL enables unstructured mesh calculations to be expressed at a high level, with a suite of automated tools to translate the scientific code into a wide range of targeted high performance implementations. OP2 makes use of various parallelisation approaches: MPI, OpenMP and CUDA. It allows for the Volna numerical algorithms to be written once, which in turn is then automatically parallelised to use multiple CPUs and GPUs. This parallelisation of Volna resulted in a huge improvement in the codes scalability and performance [RGG + 18, GKS + 20]. It has enabled the code to be used for various extremely computationally expensive tasks like building statistical emulators [BG16, LG17, GSD + 18, BG16, GHG20], carrying out sensitivity analysis and stochastic inversions [GVR + 17]. Extensive details on the parallelisation approaches used in OP2 can be found here [GMS + 11, MGR + 12] and details on the integration of Volna into OP2 are given in [RGG + 18].

Faster than Real Time Ensemble
As stated TWCs are tasked with deciding the level of threat posed by a tsunamigenic event and the areas at risk. Warning centres utilise various different tools to aid in this decision process. Leveraging modern high performance computing resources and improvements in the accuracy of earthquake detection has allowed one to simulate tsunami's across an ocean basin with run times which are faster than real time (FTRT) [LLM + 19]. The results of these faster than real time simulations are supplemental to the traditional decision matrix and precomputed database techniques that are currently utilised [GHLH13]. However, when an earthquake is detected, there is initially uncertainty with respect to various physical parameters of the event. Thus, to capture this uncertainty, the main goal is to produce faster than real time ensemble techniques. With arrival times on the order of minutes it is imperative that these ensemble results are completed in a timely fashion. Workflows are currently being developed on the implementation of these highly efficient codes into an automatically triggered response to seismic activity [LLM + 19]. In the following section Volna-OP2's ability to provide extremely rapid regional forecasts and a brief example of a potential ensemble prediction system is shown with the Mediterranean as a test case.

Mediterranean Tsunami
By way of example the major earthquake of Boumerdes 2003 (6.8 Mw) is chosen. The earthquake has been extensively studied in the literature [BDA + 04, MMC + 04, YLM + 04, SCC05]. The exact source parameters from various papers vary and will be used to highlight the ability of Volna-OP2 to provide faster than real time ensemble results. As stated in the immediate aftermath of an earthquake there is always uncertainty on the source parameters and magnitude. To capture this uncertainty, the earthquake sources presented in the literature have been scaled using Wells and Coppersmith [WC94] relations to provide initial earthquakes with magnitudes varying between 6.7-7.1 Mw. In total, there are 20 different initial earthquake sources. The initial water surface elevation is obtained by using the Okada Model [Oka85].
Maximum wave heights at each grid point are outputted after 4 hours of simulated real time. Absolute maximum and mean wave heights which incorporate the complete set are then calculated and outputted. The resolution of the grid is 2arc minutes and is thus too coarse for direct comparison to tidal gauge measurements. This grid resolution results in ≈ 385,000 cells. The run time for the complete set of 20 separate 4 hour simulations is 97s, using 2 Nvidia V100 GPUs. It should be noted that this is the run time for just the simulations themselves and doesn't include the time taken to post-process the data. Further work on these workflows is required and has been previously highlighted [LLM + 19]. From the results shown it can be seen that the larger magnitude earthquakes produce larger maximum wave heights. The regional areas which are worst affected are consistent across all the realisations, ie. mainly the local Algerian and the Balearic coastlines. The mean maximum wave heights can be deemed as the most probable distribution of maximum wave heights from this event while the absolute maximum wave heights can be seen as a worst case scenario. In the absolute maximum wave height case, non-zero wave heights are seen along the coastline of southern France and Sardinia. These 'on the fly' ensemble results could augment the warnings deduced from the traditional decision matrix or precomputed databases. Further, the absolute and mean maximum wave heights can be normalised to a traditional non-dimensional warning level scale, which in turn can be given to the relevant authorities [GHLH13].

North East Atlantic Tsunami
The last major tsunami to effect the Atlantic basin is the Lisbon 1755 event. This tsunami was generated by a ≈ 8.5 Mw earthquake which occurred on the morning of the 1 st November 1755 roughly 100km off the coast of Lisbon. The earthquake and subsequent tsunami caused a large number of fatalities and severe damage along the coasts of Portugal, Spain and Morocco. The effects of the tsunami were felt across the Atlantic basin, with historical reports of the tsunami waves arriving in the Caribbean, United Kingdom and Ireland [ODD13]. There is still no general consensus on the exact earthquake parameters for the event. With the effect of this event on the Irish coastline in mind a 'reduced' faster than real time ensemble will be shown. This is considered a 'reduced' ensemble as only six different earthquake sources were selected from the literature and only a section of the north east Atlantic was chosen. The initial sea surface was again generated using the Okada model [Oka85]. The resolution of the grid is 2 arc minutes and contains ≈ 610,000 cells. Each individual simulation is run for 8 hours and then the maximum wave heights at each grid point are outputted. The complete run time for the six different tsunami simulations is 93s on two Nvidia V100 GPUs. The plots of the maximum wave heights all highlight the impact on the local coastlines of Portugal, Spain and Morocco. They however also show a substantial variation in the directionality of the wave energy outside the gulf of Cadiz, with some sources showing the wave energy propagating towards the Irish Coastline. The absolute and mean maximum wave heights provide a rough estimate at a worst case and most likely scenario respectively. They are considered a rough estimate as the ensemble contains so few members. To investigate the effect of this event on the Irish coastline the sources introduced in the above figure will be investigated further.

Most Impactful Earthquake Source
From an Irish perspective the main tsunami threats lie in a Lisbon type event or a tsunami originating from a Mid-Atlantic submarine landslide. Extensive studies have been carried out on the collapse of the Rockall bank slide complex (located off the north-west coast of Ireland) and the subsequently generated tsunami [SGGD17]. However, this inundation mapping study focuses on earthquake generated tsunamis and from an Irish perspective the Lisbon 1755 event can be considered a worst case scenario.
Previous work on the effects of the Lisbon tsunami have highlighted the vulnerable sections of the Irish coastline to be the west (Galway Bay) and south east (Dunmore East) coastlines (areas highlighted in Fig 3). The south-east area was chosen based on a previous study [HR 06] which highlighted that this area was prone to tsunami amplification. The west coast was selected based on tsunami simulations carried out by the French tsunami warning Centre (CENALT). The CENALT simulations showed that the tsunami wave is directed into Galway Bay and amplified because of refraction from the Porcupine Bank. The reduced ensemble results (Fig  2) corroborate these findings, with the main tsunami energy impacting the west and south east coastline of Ireland. To build on the reduced ensemble, the sources which exhibited higher impact (maximum wave height) on the Irish coastline were chosen (Table 4.1). These sources were then re-run on a non-uniform mesh which is refined around the areas of interest (Fig 4). The resultant maximum wave height plots (Fig 6) and outputs at various wave gauges (Fig 5) allowed the most impactful source to be deduced. This source was then used for the high resolution inundation mapping (section 4.3) and can represent a reasonable worst case scenario.   • Source 2: The second source involves the rupture of the Horseshoe Fault (HF) and has been purposed by Baptista et al. [BMOA11] for a tsunami inundation study of the Lisbon downtown area. The source was scaled in order to obtain rupture parameters comparable to those outlined by Solares and Arroyo [ML04].
• Source 3: The third source has been described by BRGM (Bureau de recherches géologiques et minières) in a pilot study for the TANDEM project (http://www-tandem.cea.fr) and also involves the HF.
It can be seen from the maximum wave heights (Fig 6) that the overall wave energy directionality of the tsunami is consistent across sources (2 − 4). Source 1 fails to produce signs of refraction of the wave by the Porcupine Bank at the wave height scale used (Fig 6). It is interesting to note the presence of wave amplification as it propagates onto the continental shelf. This phenomenon has been observed elsewhere and has been reproduced in experimental set ups [KY18]. Based on the selected wave gauge plots (Fig 5) and the maximum wave heights (Fig 6), the candidate source was chosen to be source 3, which has been proposed by BRGM. It can be seen from the wave gauge plots (Fig 5) that this source produces the maximum wave heights along the Irish coastline, with the maximum wave height of ≈1m being produced at Dunmore East. This source will now be used for the high resolution inundation mapping simulation.
It should be noted that from these results the south coastline was also highlighted for its vulnerability to tsunami impact. Unfortunately, high resolution topography data was unavailable to the authors at the time that this pilot inundation study was carried out.

High Resolution Meshes
High resolution bathymetry/topography grids (≈ 10m) for the Galway Bay and Dunmore East areas were used in the generation of the unstructured meshes for Volna-OP2. The topography for the grids was obtained from raw LiDAR data. In order to ensure that the bathymetry and LiDAR data was continuous, all the LiDAR data was referenced to WGS84 (a global coordinate system) and then locally vertically referenced to OD Malin. Two separate high resolution (25km-10m) meshes for the Galway and Dunmore East coastlines were then generated.
The nonuniform unstructured meshes (Fig 4 and Fig7) were generated using a customised mesh sizing function and the GMSH software [GR09]. The mesh sizing function spilts the domain into three separate regions: offshore, onshore and port region (area of interest). In the offshore region the mesh size is calculated based on the bathymetry value b(x, y). Onshore cell sizes are dependent on the distance to the coastline while in the area of interest a fixed cell size is used. Full details on the mesh sizing function are given in [GHG20].

Inundation Maps
The following inundation maps highlight land areas which were inundated during the high resolution simulations (≈ 10m resolution). The grey points indicate cells which were inundated during the simulation. The maximum run-up height was calculated to be ≈ 3.4m for both the Dunmore East and Galway Bay areas. It should be noted that the simulations with these highly non-uniform meshes were extremely dependent on the Courant-Friedrich-Levant (CFL) parameter used. Too high of a CFL condition and the numerical simulation became unstable, while too low of a value and the tsunami signal became diffused and dispersed. Thus, multiple preliminary runs were needed to optimise for this parameter.

Galway Bay
Figure (9) highlights the areas of inundation for the Galway Bay region. It can be seen that the model predicts inundation at Kinvarra, Clarinbridge and Oranmore (marked on map). Further it predicts little inundation at Galway itself, the lack of high resolution bathymetry/topography in the Galway city area (Mutton Island) could be responsible for this.

Dunmore East
Figure (10) highlights the areas of inundation for the Dunmore East region. The areas of inundation include Tramore, Saltmills and near the Great Island Power Station. It should be noted that the flow in the River Barrow was not modelled in these simulations. This river flow could lead to the formation of a bore and thus allow the tsunami wave to propagate further up the river. The maximum water height at the Waterford coastline is predicted to be ≈ 1.2m. The maximum run-up height was calculated to be ≈ 3.4m. Figure 10: Inundation maps for Dunmore East: the grey areas on the map are land values which were inundated during the simulation.

Concluding remarks
The massively parallelised Volna-OP2 code is shown here to be a candidate code which could be implemented into a Tsunami Warning Centre's workflow. It is shown to produce extremely rapid regional forecasts for the Mediterranean Sea. The 2003 Boumerdes earthquake was chosen to act as an example for this process. Most tsunami warning centres now have access to substantial pre-computed databases of potential earthquake sources (1,000s of unit sources) and resultant tsunami wave heights [GHLH13]. Once a tsunami is detected, any relevant sources in the database can become an input into an 'on the fly' ensemble with Volna-OP2 being the tsunami code kernel. The resultant absolute maximum and mean wave height plots can augment the results from the pre-computed database. The results presented here for the 'on the fly' ensembles are acknowledged to be preliminary and major effort in developing efficient workflows around the code is needed. However, it is clear that Volna-OP2 has promising potential when it comes to this respect.
The codes capability at producing high resolution inundation studies is explored in section (4). This pilot inundation study highlights recent efforts to increase the tsunami hazard assessment for the Irish coastline. As stated, the main tsunami hazard for Ireland lies with a repeated Lisbon 1755 type event or a mid-Atlantic submarine landslide. The inundation study presented here assumes the Lisbon 1755 event as a worst case scenario and examines the potential inundation of sections of the Irish coastline from such a scenario. The impacts of various earthquake sources proposed in the literature are explored and the one which produced the highest maximum wave heights for the Irish coastline and at various gauge points was chosen for the inundation study. Utilising non-uniform unstructured meshes a high resolution simulation of the areas of interest was achieved with Volna-OP2. It was found that these high resolution inundation simulations were extremely dependent on the CFL condition, and thus required multiple preliminary runs to optimise for.
The tsunami hazard assessment of the Irish coastline is considered a pilot study on the extent of inundation for the most vulnerable sections of the Irish coastline from a Lisbon type event. The results shown corroborate previous studies on what sections of coastline (the west and south east) are particularly vulnerable. The assessment found that inundation was constrained to low-lying areas in the Dunmore East and Galway Bay regions, with a maximum run-up of ≈ 3.4m for both regions. The depth of water above the underlying topography (inundation depth) reached 2m in parts. This level of inundation could pose a risk to various infrastructures (Great Island Power Station) and members of the public located near the shoreline. To date the above pilot study provides the best estimate of assessing the hazard associated with a Lisbon type tsunami event on the Irish coastline. However, it is acknowledged that to provide a full hazard assessment of the Irish coastline, some further sites (south coast) and factors should be included, such as tidal forcing, 'street level' inundation and the effect of wave dispersion. Thus, the inundation plots produced should only be used as a guidance on the potential tsunami hazard of an area.