1. Introduction
After the devastating 2004 Indian Ocean tsunami, there was a concerted effort by the scientific community to mitigate the damages posed by these natural hazards on coastal communities. After the event Bernard et al. [
1] wrote that the greatest contribution that science can make to society is to serve it purposefully. Prior to the 2004 tsunami, which touched virtually all the global coastlines [
2], there existed a earthquake-centric tsunami warning network in the Pacific Ocean [
3]. However, the global threat posed by these natural hazards was brought into sharp focus by the 2004 event. Warning centres were set up across all the major ocean basins, with numerous international collaborations being formed. The scientific community’s efforts were focused on the continued development of tsunami early warning centres, carrying out prior hazard assessments, deploying tsunami wave gauges and increasing public awareness [
4]. Despite the global average of one tsunami occurring per year [
1], on a regional level tsunamis are inherently rare events. The tsunami community thus relies on high-fidelity computational/mathematical models to carry out simulations of numerous potential tsunami sources in order to quantify the associated hazards.
The last major tsunami to effect the Atlantic basin is the Lisbon 1755 event. This tsunami was generated by a
Mw earthquake which occurred on the morning of the 1 November 1755 roughly 100 km 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 [
5]. 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 [
6]. There is still no general consensus on the exact earthquake parameters for the event. There exists a substantial catalogue of work which focuses on the earthquake source for this event [
7,
8,
9] and numerous studies on the effects of the resultant tsunami waves in the near and far field areas [
5,
10,
11].
Focusing on the Irish coastline there exists a handful of historical accounts on the effects of this event. One notable account is given by the harbour master in Kinsale (South Coast of Ireland) [
6]. There is however a severe lack of studies quantifying the level of hazard associated from the event on a high resolution level. The Defra study [
12] utilised a minimum mesh element size of 1 km along the southern Irish coastline and did not capture inundation as a cut off bathymetry depth was used. However, one of the relevant findings from this study was that tsunami amplification occurred along the South East coastline of Ireland. Further, some studies from the French tsunami warning centre (CENALT) found that the west coast was vulnerable to a Lisbon-type event, with the tsunami waves refracted by the Porcupine bank into Galway Bay and amplified. Therefore the efforts presented here in quantifying the hazard for the Irish coastline on a high resolution level are the first of its kind. The high resolution simulation results will aid in the disaster management planning for the identified vulnerable areas.
However, in the immediate aftermath of a tsunami occurrence, it is a tsunami early warning centre’s (TWC) responsibility to deduce the level of hazard posed by the event. TWCs’ tasks include detecting the tsunami source, deducing the level of threat posed, deciding on the areas most at risk and then notifying the relevant authorities. As tsunamis can travel at extremely high speeds in deep water arrival times at the coastline are often on the
(mins). TWCs are therefore under severe time constraints when carrying out their duties. To aid in this decision process, they currently deploy a variety of tools which produce regional levels of warnings. These tools are primarily based on extensive pre-computed databases and pre-calculated scaling laws [
13]. However, as the computational efficiencies of tsunami codes increase, faster than real time techniques are now being developed. These promise to be extremely beneficial as multiple tsunami simulations constrained by the tsunami source can be simulated in an ensemble. This ensemble would capture in a faster than real time setting the uncertainty on the resultant tsunami waves. A reduced ensemble for the North-East Atlantic is presented in
Section 3.
The local bathymetry can have a substantial effect on the local levels of maximum wave heights, inundation and induced currents, something which is missed by the regional warnings. So in order to gain an understanding of how the local bathymetry affects an incoming tsunami, prior high resolution risk/hazard assessments are often carried out, such as the results presented in
Section 4. These studies can require extensive computational resources and are time consuming. Epistemic uncertainty involved in the modelling and source parameters has led to the development of probabilistic tsunami hazard assessments (PTHA) [
14,
15]. One of the key limiting factors for both real time tsunami forecasts and risk assessment studies are the runtimes required to run the modelling software. Extremely efficient tsunami codes, which utilise 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 [
16,
17,
18]. The code solves the depth-averaged nonlinear shallow water equations, utilising modern reconstruction and shock capturing techniques. The original code was introduced in 2011 [
19] and has since undergone extensive parallelisation [
20]. Both the original and newly updated version have been verified against the traditional tsunami mitigation benchmark problems [
20] and an extensive error analysis of the newly parallelised version has been carried out [
21]. Volna-OP2 is presented here 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.
This paper is split into three parts; in
Section 2 an outline of the numerical scheme and parallelisations used in Volna-OP2 are described. In
Section 3, two faster than real time reduced ensembles are presented for the North East Atlantic Ocean, with the difference in ensembles being the level of resolution focused along the Irish coastline. Building upon these faster than real time simulations, in
Section 4 a high resolution inundation mapping study for vulnerable sections of the Irish coastline is shown.
2. 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
is the total water depth,
is the time-dependent bathymetry,
is the free surface deformation,
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 2nd 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 2nd or 3rd 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 [
19,
21].
When Volna was first developed in 2011 [
19] 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 [
22]. 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 are then automatically parallelised to use multiple CPUs and GPUs. This parallelisation of Volna resulted in a huge improvement in the codes scalability and performance [
20,
21]. It has enabled the code to be used for various extremely computationally expensive tasks like building statistical emulators [
17,
18,
23,
24], carrying out sensitivity analysis and stochastic inversions [
25]. Extensive details on the parallelisation approaches used in OP2 can be found in References [
22,
26] and details on the integration of Volna into OP2 are given in Reference [
20].
4. High Resolution Inundation Mapping
In order to provide accurate information on the inundation extent of such an event on the Irish coastline, higher resolution meshes must be generated and used. At present it is unfeasible to carry out these high resolution simulations with Volna-OP2 with run times which would be beneficial to a warning centre. The explicit temporal scheme is the main constraint on this for Volna-OP2, as the time step is governed by the minimum mesh element size.
In order to generate the high resolution meshes one must first obtain bathymetry/topography data (≈ 10 m) for areas of interest (Galway Bay and Dunmore East). The topography for the grids was obtained from raw LiDAR data. In order to ensure that the high resolution data was continuous with the regional data, all the LiDAR/near shore data was spatially (latitude/longitude) referenced to WGS84 (a global coordinate system). The raw LiDAR and near shore bathymetry data were vertically referenced to OD Malin (vertical datumn used in Ireland), thus an assimilation step was required to remove the discontinuity with the regionally defined data. Further, due to the difficulty in obtaining bathymetry data in shallow regions, gaps between the near shore and the onshore LiDAR data were present. Some smoothing filters were applied to remove these gaps and to produce continuous bathymetry/topography grids (
Figure 6 and
Figure 7). The resultant grids were used in the generation of the high resolution mesh (
Figure 8).
It should be noted that the relevant tidal levels, which could have a substantial effect on the extent of inundation, are not captured in this study. The authors acknowledge that to provide a fully comprehensive hazard assessment, these tidal effects along with other factors discussed in
Section 5 need to be investigated.
The nonuniform unstructured meshes (
Figure 3 and
Figure 8) were generated using a customised mesh sizing function and the GMSH software [
30]. The mesh sizing function splits 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
. 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 Reference [
18].
In the interest of minimising the runtime a ’hotstart’ technique was employed. From the simulations in
Section 3.2, the dependent variables (
) were outputted after 4 h of simulated time. These variables were interpolated onto the fine mesh (
Figure 8) and therefore acted as the initial conditions. The simulations on the fine mesh were then carried out for a further 6 h, with the runtime for all four fine mesh simulations taking a total of 1 h 59 min, using 2 Nvidia V100 GPUs.
4.1. Inundation Maps
The following inundation maps highlight land areas which were inundated during the high resolution simulations (10 m resolution). The maximum run-up height was calculated to be ≈ 3.4 m 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 this parameter.
4.1.1. Galway Bay
Figure 9 highlights the areas of inundation for the Galway Bay region from each of the four sources simulated. The grey points in
Figure 9 indicate cells which were inundated during the various simulations. It can be seen that the model predicts inundation at Kinvarra, Clarinbridge and Oranmore (marked on map) (
Figure 10). 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. Source 1 causes the least amount of inundation while there is a consistency in the areas and extent of inundation for Sources 2–4.
4.1.2. Dunmore East
Figure 11 highlights the areas of inundation for the Dunmore East region from each of the four sources simulated. The grey points in
Figure 11 indicate cells which were inundated during the various simulations. The areas of inundation include Tramore, Saltmills and near the Great Island Power Station (marked on map) (
Figure 12). It should be noted that the flow in the River Suir 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.2 m. The maximum run-up height was calculated to be ≈ 3.4 m. Similar to the Galway bay results presented above, Source 1 causes the least amount of inundation while there is a consistency in the areas and extent of inundation for Sources 2–4.
5. Discussion and Concluding Remarks
In line with the goals and directives set out by the tsunami community after the devastating tsunami events of the 2004 Indian Ocean and 2011 Tohuku tsunamis, this paper presents new attempts to quantify the tsunami hazard associated for sections of the Irish Coastline. As the 1755 Lisbon tsunami was the last major tsunami to affect the Irish coastline, it serves here as a worst case scenario. Multiple earthquake scenarios presented in the literature are simulated, with the inundation levels at vulnerable sections of the Irish coastline investigated for the four most impactful sources. Faster than real time tsunami simulations are achieved with the massively parallelised Volna-OP2 tsunami code. These faster than real time simulations will greatly improve the warning capabilities of tsunami early warning centres, with the Volna-OP2 code presented as a candidate code which could be implemented into a tsunami warning centre’s workflow.
In the immediate aftermath of a tsunamigenic source being detected there is always some uncertainty on the physical parameters of the source. Thus, multiple tsunami realisations are needed to capture the spread of potential tsunami wave distributions. Any potential sources 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 tsunami warning centres’ pre-computed database. It is shown here with a uniform mesh resolution of ≈ 3.3 km that 6 separate tsunami simulations can be simulated using 2 Nvidia V100 GPUs in 93 s. However in order to refine these regional forecasts a non-uniform mesh with a minimum element size of 100 m is used. The four most impactful sources identified in the initial ensemble are simulated again for 10 h with a runtime of ≈ 21 min (using 2 Nvidia V100 GPUs). Although this runtime is outside the scope of being beneficial for a tsunami warning centre, more computing resources would contribute to reducing these times. The refined regional forecasts would allow a warning centre to issue more accurate warnings. These regional results corroborate previous studies on what sections of Irish coastline (the west and south east) are particularly vulnerable to a repeated Lisbon-type event. The authors acknowledge that the ’on the fly’ ensembles presented here are 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.
There exists a severe lack of studies in the literature on the hazard assessment associated with a Lisbon-type tsunami on the Irish coastline. Therefore, the high resolution tsunami hazard assessment presented here is considered a pilot study on the extent of inundation for the most vulnerable sections of the Irish coastline from a Lisbon-type event. High resolution simulations of the areas of interest were achieved by utilising a non-uniform unstructured mesh and simulating the tsunamis 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 assessment found that there was a consistency in the areas and levels of inundation across the four sources simulated, with only source 1 producing a notably less amount of inundation. The inundation extent was constrained to low-lying areas in the Dunmore East and Galway Bay regions, with a maximum run-up of m for both regions. The depth of water above the underlying topography (inundation depth) reached 2 m 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, 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.