Initialising landslide-generated tsunamis for probabilistic tsunami hazard assessment in Cook Strait

Abstract

The Cook Strait Canyon is a submarine canyon which lies within 10 km of Wellington, the capital city of New Zealand. The canyon flanks are scarred with the evidence of past landslides that may have caused large local tsunamis. City planning and civil defence management require information on the magnitude and frequency of these tsunamis to adequately plan for them. Submarine-landslide-generated tsunamis are by nature local features. While they may be catastrophic in the near field, they are generally far smaller scales than co-seismic tsunamis and their energy does not travel very far. Including them within a comprehensive tsunami hazard assessment requires accounting for a large number of potential landslide sources. Unless we only use simple rules of thumb to approximate tsunami height, this requires considerable computing power. This article describes a technique for expanding two-dimensional vertical-slice tsunami generation by landslide modelling into a two-dimensional horizontal surface which can be used for tsunami propagation and inundation modelling. As such, it spans the gap between full three-dimensional modelling of the landslide and simple initialisation.

Keywords

Tsunami submarine landslides Cook Strait hazard assessment

Introduction

The 2004 Indian Ocean Tsunami was a wake-up call to tsunami scientists around the world. It reminded us that we should not assume that what has occurred in the recent past limits what we can expect in the future. Tsunami return periods are very long, there may be little record left of previous large events. Even in Japan, where tsunamis are relatively frequent, the more recent Tohoku-oki event showed how it is easy to underestimate the magnitude of possible events. This does not mean that we should give up trying. Rather, we need to make best use of all available evidence to actively move beyond the boundaries of what is currently known. Bringing together knowledge of what has occurred in other places with geologic knowledge of our own region, palaeo-tsunami knowledge and modelling, we can estimate hazards facing us that have not occurred within recorded history. This is especially important in places like New Zealand, where the historical tsunami record is less than 200 years long (Downes, 2014).

While co-seismic tsunamis (tsunamis caused by earthquakes) are the most common and well known type of tsunami, they are not the only source of tsunamis. Around 7% of tsunamis are associated with submarine landslides (Harbitz et al., 2014).

The Papua New Guinea tsunami of 1998 was a dramatic reminder of this when a tsunami after a moderate earthquake caused a wave over 15 m high in Sissano Lagoon killing over 2000 people (Tappin et al., 2001). Later research showed that the majority of this tsunami was not caused by the earthquake, but by a submarine landslide triggered by one of the larger aftershocks (Tappin et al., 2008).

New Zealand is a geologically complex, tectonically active region. It lies on the interface between the Pacific plate and the Indo-Australian Plate. To the north, the Pacific Plate is subducting beneath the Indo-Australian Plate in the Kermadec and then Hikurangi subduction zones. Just south of Cook Strait, this plate boundary becomes strike-slip along the alpine fault and then further south the Indo-Australian Plate is subducted beneath the Pacific Plate. Thus, the Wellington region is compressed and this has formed the North Island fault system. Numerous faults cross the region, many of them trending NNE-SSW parallel to the subduction zone. Many of these faults extend under Cook Strait and have the potential to generate tsunamis when they rupture. The last major earthquake in the Wellington region (the 1855 Wairarapa earthquake, M_w 8.2–8.4) caused a tsunami with local wave run-up over 10 m and affected the entire region (Downes, 2014). Palaeo-tsunami evidence also points to catastrophic past tsunamis in Wellington (Goff and Chague-Goff, 2009).

But it is not just co-seismic tsunamis that have the potential to affect the Wellington region. The Cook Strait Canyon is a deeply incised submarine canyon which cuts into the continental shelf to within 10 km of Wellington, the capital city of New Zealand (see Figure 1). The canyon walls reveal evidence for long-term landslide activity in the form of landslide scars. The geometry of the landslide scars, as well as evidence from seismic reflection data, indicates that the landslides are deep-seated failures occurring within ‘soft-rocks’ or over-consolidated tertiary marine sedimentary sequences as opposed to young surficial sediments (Micallef et al., 2012; Mountjoy et al., 2009).

Figure 1.

New Zealand with inset showing southern part of Wellington and Cook Strait Canyon. Landslide scars can be seen on the sides of the canyon.

The canyon’s location in Cook Strait means there is inhabited land in the path of both forward and backward propagating waves. To get a full understanding of the hazards affecting their region, city planners, civil defence managers and coastal communities need quantifiable information on how submarine-landslide-generated tsunamis compare with co-seismic tsunamis and other potential hazards, including information on both the magnitude and the frequency of this type of tsunami.

Probabilistic tsunami hazard assessments (PTHAs) tend to focus on tsunamis generated by earthquakes, for example, Power (2013) includes a PTHA of New Zealand but only including co-seismic sources. It is only recently that other sources, such as landslides, have started to be seriously considered (Grilli et al., 2009; Leonard et al., 2014). Often, this omission is a matter of practicality. Subduction faults are generally well known and mapped. There are a limited number of scenarios that cover the majority of the earthquake-based sources. By comparison, there are many more landslide sources that may affect at least part of the region of interest because landslide-generated tsunamis are by nature local features. While they may be catastrophic in the near field, they are generally far smaller scales than earthquake-generated tsunamis, and so their effects are limited to tens of kilometres from the source. The Sissano Lagoon tsunami that occurred in Papua New Guinea in 1998 caused wave heights over 10 m high over a stretch of coastline 10–20 km long, but outside of that region, the wave heights dropped off relatively rapidly with heights being around 2 m about 50 km away (Tappin et al., 2001). The Sissano Lagoon landslide tsunami is a good analogy for the type of event that could occur in Cook Strait, due to similar landslide size and water depth. Including landslide-generated tsunami within a comprehensive tsunami hazard assessment requires accounting for a large number of potential landslide sources. Unless we only use simple rules of thumb to approximate tsunami height, this requires considerable computing power.

As in many places in the world, the Indian Ocean Tsunami prompted a review of the tsunami hazard in New Zealand (Berryman, 2005). As part of this review, landslides in submarine canyons close to the New Zealand coast were identified as a potential hazard that required more research. A series of two projects were commissioned by New Zealand’s Natural Hazards Research Platform (NHRP) to investigate this hazard and to develop a methodology for calculating probabilistic wave heights at the coast to being it in line with Power (2013).

These NHRP research projects were multi-disciplinary, bringing together geological, physical and mathematical modelling knowledge to attempt to answer these questions. Scientific cruises carefully mapped the sea floor and identified landslide scars, core sampling to measure rock sample strength and thus to estimate failure coefficients, and seismic shaking maps were developed. All these were used to understand the probability of landslide failure and the likely magnitude frequency relationship. As well as understanding the likelihood of landslide occurrence, we need to know the size of tsunami a landslide can generate and what its impact is on the coastline. We can use mathematical modelling to fill this gap: initialising the tsunami based on the characteristics of the landslide and then modelling its propagation. This article looks at the process of modelling the landslide tsunami and in particular initialising the tsunami taking into account the characteristics of the landslide.

Modelling

The adaptive partial differential equation (PDE) solver Gerris was used for the modelling portion of the study (Popinet, 2003). Gerris is a very flexible numerical solver based on a quad-tree (oct-tree in three dimensions) mesh that can solve a range of equations including Navier–Stokes and Saint Venant. It has previously been used to model tsunamis (Popinet, 2011, 2012) and ocean dynamics (Popinet and Rickard, 2007) in addition to many other computational fluid dynamics (CFD) problems. The adaptive mesh allows for increased resolution in areas where it is needed, and this adaptivity can change to track a moving structure, for example, Gerris is able to track a tsunami front as it crosses the ocean, in turn refining and coarsening the cells as needed.

Unlike co-seismic tsunamis, there is not a golden standard method of initialising submarine-landslide-generated tsunami such as Okada (1985). Watts et al. (2003) developed TOPICS (Tsunami Open and Progressive Initial Conditions System), an initialisation method based on parameterising modelling, but this was developed for the continental slope. It is not always applicable to the situation in a canyon. Factors such as the canyon bottom and the steepness of far wall – not accounted for in TOPICS – can alter the size and shape of the tsunami wave generated. Furthermore, in the case of a landslide in a canyon, the landslide might be facing towards the land area of interest which is radically different from the situation where the failure happens on the continental shelf. For a landslide in a canyon, the landslide usually decelerates by hitting the bottom of the canyon (perhaps riding up the opposite side slightly), whereas in the case of a landslide on the continental shelf, the initialisation is done when the landslide reaches terminal velocity. Because of this, we decided to model the tsunami generation process explicitly.

We used Gerris for both the landslide and the tsunami modelling components. For landslide modelling, we used a volume of fluid (VoF) formulation of the Navier–Stokes equations (Hirt and Nichols, 1981; Scardovelli and Zaleski, 1999). This enables us to solve multiphase flows (Agbaglah et al., 2011) such as modelling the air–sea interface explicitly. For the tsunami modelling, we used Gerris to solve the Saint Venant equations (also known as the non-linear shallow water equations). These equations assume the tsunami waves are long enough so that essentially the depth is shallow everywhere and that there is no dispersion. Although not strictly true in this case, studies have shown that dispersion generally affects the timing of the waves more than run-up heights (Pederson, 2008). Basilisk, the successor to Gerris which is currently being developed, is able to solve the Serre–Green–Nagdhi equations (also known as the non-linear Boussinesq equations; Popinet, 2015), which include weakly dispersive effects. It could be used in future studies to investigate the importance of these effects.

Initial investigations of landslides in canyons showed that factors such as the position and steepness of the far wall of the canyon could play a strong role in the generation of a tsunami. The Cook Strait Canyon varies widely in steepness, width and depth. Landslides can face any direction, and in some places there is land both in front of and behind the landslide. There is no small set of typical scenarios that cover the entire hazard. We need to consider landslides occurring throughout the canyon. To achieve this, the walls of the entire Cook Strait Canyon were divided into 176 locations, approximately 1 km apart, shown in Figure 2. Furthermore, landslides in the Cook Strait Canyon can range in size up to 1–2 km³. In order to understand how the volume of the landslide affects the tsunami generated, we modelled three landslide volumes (0.1, 0.3 and 1 km³) at each location: 528 scenarios in total. Thus, we needed a computationally efficient initialisation technique to allow this number of simulations.

Figure 2.

Close-up of the Cook Strait Canyon with cross-sections where landslide scenarios were modelled.

While tsunami modelling is standard within Gerris, the solver had not previously been used to model the process of landslides generating tsunamis. To ensure that Gerris was able to correctly solve these problems, we first modelled test cases from Enet et al. (2003), Enet and Grilli (2007) and Liu et al. (2005). Figure 3 shows time steps from modelling of scenario taken from Liu et al. (2005). The adaptive resolution of the mesh can be seen projected onto the slope. In this, the landslide is modelled as a solid block with a prescribed motion and the VoF formulation solves for the air and the sea as fluids with different densities; we also use a three-phase VoF formulation where the landslide was also modelled as a semi-rigid, dense fluid. In the latter case, the landslide is able to evolve under gravity and does not need a prescribed motion. Although we are able to solve these three-dimensional (3D) landslides generating tsunamis, they are computationally expensive and we needed a more efficient method in order to model all the scenarios.

Figure 3.

Modelling landslide as solid block sliding down slope causing tsunami. The grey top layer is the water surface. The mesh on slope shows the adaptive refinement, with higher levels of refinement around the sliding block and the shoreline.

Our initial investigations were modelled as two-dimensional (2D) vertical slices. The bathymetry was represented as a simplified version with a flat shelf on either side of the canyon and a near-side slope which the landslide would slide down and then a far-side slope. This 2D modelling is equivalent to assuming that the canyon and the landslide extend indefinitely without change in the cross-ways direction. While these give many insights into the characteristics of landslide into canyons and are far more computationally efficient, we cannot use this directly to initialise our tsunami propagation model as they only have a one-dimensional water surface. To understand the relationship between the 2D and 3D modelling, we investigated a specific example of a submarine landslide generating a tsunami shown in Figure 4. The left-hand-side of this figure shows the scenario it is based on and the right-hand-side shows the 2D and 3D formulations of this scenario. We considered three different landslide widths, 1000, 1700 and 2500 m, to explore how the finite width of the landslide affects the results.

Figure 4.

Simplified landslide modelling. The left-hand-side shows the landslide taken as the model scenario. The top right-hand-side is the 2D vertical-slice representation of this landslide, showing the sea floor with the canyon and the wedge-shaped landslide. Bottom right-hand-side shows the 3D representation.

Figure 5 shows the results from the 2D vertical-slice modelling (which we shall refer to as η_2D) along with the surfaces produced from scenarios assuming a finite. The top left-hand-side shows the wave height in cross-section through y = 0. The overall height of the wave decreases as the width of the tsunami decreases. As the width increases to infinitely wide, the solution tends to the 2D vertical-slice solution. The other three figures show the wave surfaces for the three cases with the 2D vertical-slice solution drawn in the y = 0 plane for comparison. The 2D vertical slice captures the basic outline of the 3D modelling. However, if we are to use the 2D modelling to initialise the tsunami modelling, we need to convert the results from the 2D vertical-slice model into a 2D horizontal surface that will be the initial condition for the tsunami modelling.

Figure 5.

Comparison of 3D and 2D vertical-slice modelling for scenario shown in Figure 4 for three different widths: 1000, 1700 and 2500 m.

The characteristic timescale of the landslide, t₀, is assumed to be the time it takes the landslide to reach the bottom of the canyon (measured here as when the y component of the centre of mass reaches a minimum). The characteristic wavelength, λ₀, is then calculated from t₀ by

λ_{0} = t_{0} \sqrt{g H_{s h e l f}}

(1)

where H_shelf is the depth of the shelf and g is acceleration due to gravity. We take the water surface, η_2D, at either twice the characteristic time length or just before the wave reaches the edges of the modelling domain (in which case the wave would be getting close to the shore), whichever occurs first. This gives the starting cross-section to use as our initialisation for the tsunami propagation modelling.

The 2D vertical-slice modelling is equivalent to assuming that the landslide and the canyon extend indefinitely in the along-canyon direction. When the landslide has finite width, the height of the wave is decreased as can be seen in Figure 5. Watts et al. (2005) give a formula relating the 2D maximum wave height with the 3D maximum wave height based on the width of the landslide and its characteristic wavelength

W_{2 D - 3 D} = \frac{w i d t h}{(w i d t h + λ_{0})}

(2)

where width is the width of the landslide and λ₀ is the characteristic wavelength as calculated in equation (1). Note that as the width goes to infinity, this limits to 1. Thus, $η_{m a x, 3 D} = W_{2 D - 3 D} \times η_{m a x, 2 D}$ . We use equation (2) to adjust our η_2D values to take the width into account.

As the landslide slides into the canyon, it pushes a bow wave of water in front of it. It also creates a hole above its initial starting point where the water is unable to flow in fast enough to account for the submarine movement. The disturbance then radiates outwards from the starting point of the landslide. To create the initial water disturbance from our water heights in the along-transect direction, we first split the water heights into two pieces η₊ and η₋ depending on whether it is that part in front of or behind the centre of the starting point of the landslide. The part in front of the landslide is assumed to have been pushed by the landslide and so its shape is dependent on the width of the landslide in three dimensions. The part behind the landslide is assumed to be predominantly due to the radiating out of the water disturbance. Thus, we set

R_{+} (r, θ) = W_{2 D - 3 D} \times η_{+} (r) \times e^{- (\frac{θ^{2}}{{(\tan^{- 1} (w i d t h / l e n g t h))}^{2}})}

(3)

and

R_{-} (r, θ) = W_{2 D - 3 D} \times η_{-} (r)

(4)

We combine equations (3) and (4) to get our initial water disturbance, namely

η_{3 D} = R_{+} (r, θ) + R_{-} (r, θ)

(5)

Once we have our map view representation of the water surface disturbance, we then rotate the axes of the surface so the x-axis aligns with the along-transect axis and shifts the disturbance in space so that it is correctly oriented and occurs in the right part of the domain.

Results

In order to check out this initialisation, we looked at cases where we used the initialisation from the full 3D VoF modelling and compared this with the 2D initialisation. When we did the full 3D VoF modelling, we were also able to output the vertically integrated velocity at the initial time-point as well as the wave surface, so we considered three cases: initialisation using both the wave surface and the vertically integrated velocity from the 3D modelling, initialisation using only the wave surface from the 3D modelling and initialisation using the extended wave surface from the 2D modelling. We present the results comparing the first case (seen as the ‘truth’ here) with the third case for different widths of landslide: 1000, 1700 and 2500 m. Figure 6 shows the water surfaces for the 2D approximation and the 3D modelling. Note that these have been rotated so that they align with the landslide direction in the case shown in Figure 4.

Figure 6.

Initialisation for tsunami from submarine landslide. Top row: 2D approximation; bottom row: 3D for width (left to right) 1000, 1700 and 2500 m. Grey scale ranges from −30 m (white) through to 30 m (black). Initialisation has been rotated so that it faces in the direction of the landslide (see Figure 4).

Figure 7 shows the maximum wave heights for tsunamis generated by landslide with widths 1700 and 2500 m initialised using the full 3D initialisation including velocity and by the 2D approximation. Grey scale ranges from 0 to 20 m in all cases. It can clearly be seen that the 2D approximation lacks some of the details of the full 3D initialisation, but the overall pattern is a good representation. Specifically, we are interested in the wave heights at the coast, so in Figure 8, we are only looking at points around the Wellington coastline. The inner colour represents the 2D approximation and the outer colour the 3D initialisation. Where only one colour is apparent shows a good match. The same information is shown without the spatial locations but with maximum wave height on the y-axis for comparison in the bottom right-hand-side. Again, the general pattern is well represented, although there are some discrepancies in the exact details. Table 1 shows the percentage of wave heights at coast for the 2D approximation that lie within 0.5 and 1 m of the 3D simulation. Time series for three points are shown in Figure 9. Again, there is generally good agreement, although the timing is not correct for Lake Ferry. But even in this case, the wave height is generally in agreement.

Figure 7.

Comparing maximum wave height from scenarios initialised using full 3D and 2D vertical-slice VoF modelling. Two slide widths (1700 and 2500 m) are considered.

Figure 8.

Maximum wave height at coast plotted on positions around Wellington coastline for widths 1000 m (top left), 1700 m (top right) and 2500 m (bottom left); outer colour is 3D value and inner colour is 2D value; grey scale ranges to 10, 15 and 20 m, respectively. The same information is plotted bottom right without the spatial information and with maximum wave height on y-axis.

Table 1.

Percentage of errors of maximum wave height at coast.

Width (m)	\|Error\| < 0.5 (%)	\|Error\| < 1 (%)
1000	93	98
1700	86	96
2500	83	96

Figure 9.

Time series for 3D and 2D initialisation at three points around Cook Strait.

Discussion and conclusion

While tsunami simulations initialised with the 2D approximation rather than the full 3D representation do not capture all the details of the propagation, they do capture the first-order patterns and effects. In addition, they are considerably faster to calculate, making it possible to model for all 528 scenarios needed to capture the variability of tsunamis generated by landslide in the Cook Strait Canyon. Given that the results of the tsunami propagation are being aggregated to provide probabilistic wave heights at the coast and that the variability between landslides from different regions in the canyon is far greater than the differences between the 2D and 3D initialisation, we judge this to be an acceptable method for initialising the tsunami scenarios.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the New Zealand Natural Hazards Platform (contract CS-GNS-25).

Author biographies

Dr Emily M Lane works at the National Institute of Water and Atmospheric Research (NIWA, Taihoro Nukurangi) in New Zealand studying coastal hazards - primarily inundation caused by tsunamis and storm surge. She has a specific interest in tsunamis generated by submarine landslides. She is a member of the Tsunami Expert Panel who provide expert advice to Civil Defence during tsunami alerts.

Dr Joshu J Mountjoy is a marine geologist at the National Institute of Water and Atmospheric Research (NIWA, Taihoro Nukurangi) in New Zealand. He has worked extensively on submarine canyon systems in New Zealand and overseas, researching quaternary activity of sedimentary systems and the long term influence of climate and tectonics on canyon development.

William L Power is a tsunami scientist working at GNS Science in New Zealand with a specific interest in probabilistic tsunami hazard. He compiled the 2013 Review of Tsunami Hazard in New Zealand and is a member of the Tsunami Expert Panel.

Stephane Popinet is a directeur de recherche at Centre National de la Recherche Scientifique (CNRS) based at Institut Jean le Rond d’Alembert of Université Pierre et Marie Curie (UPMC), Paris, France. He is an internationally renowned researcher in quad-tree adaptive grids and has created two state-of-the-art multiphase computational fluid dynamics models using these techniques: Gerris and its successor, Basilisk.

References

Agbaglah

Delaux

Fuster

. (2011) Parallel simulation of multiphase flows using octree adaptivity and the volume-of-fluid method. Comptes Rendus Mécanique 339(2–3): 194–207.

Berryman

(2005) Review of tsunami hazard and risk in New Zealand. Institute of Geological & Nuclear Sciences client report 2005/104, p. 139 Available at: http://www.civildefence.govt.nz/assets/Uploads/publications/GNS-CR2005-104-review-of-tsunami-hazard.pdf

Downes

(2014) New Zealand tsunami database: Historical and modern records. Available at: http://data.gns.cri.nz/tsunami/

Enet

Grilli

(2007) Experimental study of tsunami generation by three-dimensional rigid underwater landslides. Journal of Waterway Port Coastal and Ocean Engineering-ASCE 133(6): 442–454.

Enet

Grilli

Watts

(2003) Laboratory experiments for tsunamis generated by underwater landslide: Comparison with numerical modeling. In: Proceedings of the thirteenth (2003) international offshore and polar engineering conference, 25–30 May, 2003 pp. 372–379.

Goff

Chague-Goff

(2009) Brief Communication: Cetaceans and tsunamis – Whatever remains, however improbable, must be the truth? Natural Hazards and Earth System Sciences 9(3): 855–857.

Grilli

Taylor

O-DS

Baxter

CDP

. (2009) A probabilistic approach for determining submarine landslide tsunami hazard along the upper east coast of the United States. Marine Geology 264(1–2): 74–97.

Harbitz

Lovholt

Bungum

(2014) Submarine landslide tsunamis: How extreme and how likely? Natural Hazards 72(3): 1341–1374.

Hirt

Nichols

(1981) Volume of fluid (VOF) method for the dynamics of free boundaries. Journal of Computational Physics 39(1): 201–225. Available at: https://dx-doi-org.web.bisu.edu.cn/10.1016/0021-9991(81)90145-5

10.

Leonard

Rogers

Mazzotti

(2014) Tsunami hazard assessment of Canada. Natural Hazards 70(1): 237–274.

11.

Liu

PL-F

T-R

Raichlen

. (2005) Runup and rundown generated by three-dimensional sliding masses. Journal of Fluid Mechanics 536: 107–144.

12.

Micallef

Mountjoy

Canals

. (2012) Deep-Seated Bedrock Landslides and Submarine Canyon Evolution in an Active Tectonic Margin: Cook Strait, New Zealand. In: Yamada

Kawamura

Ikehara

. (eds) Submarine Mass Movements and Their Consequences, vol. 31. Dordrecht: Springer, pp. 201–212.

13.

Mountjoy

Barnes

Pettinga

(2009) Morphostructure and evolution of submarine canyons across an active margin: Cook Strait sector of the Hikurangi Margin, New Zealand. Marine Geology 260(1–4): 45–68.

14.

Okada

(1985) Surface deformation due to shear and tensile faults in a half-space. Bulletin of the Seismological Society of America 75(4): 1135–1154.

15.

Pederson

(2008) Modeling runup with depth integrated equation models. In: Liu

PL-F

Yeh

Synolakis

(eds) Advances in Coastal and Ocean Engineering: Advanced Numerical Models for Simulating Tsunami Waves and Runup. Singapore: World Scientific, pp. 3–41.

16.

Popinet

(2003) Gerris: A tree-based adaptive solver for the incompressible Euler equations in complex geometries. Journal of Computational Physics 190(2): 572–600.

17.

Popinet

(2011) Quadtree-adaptive tsunami modelling. Ocean Dynamics 61(9): 1261–1285. Available at: https://dx-doi-org.web.bisu.edu.cn/10.1007/s10236-011-0438-z

18.

Popinet

(2012) Adaptive modelling of long-distance wave propagation and fine-scale flooding during the Tohoku tsunami. Natural Hazards and Earth System Sciences 12(4): 1213–1227.

19.

Popinet

(2015) Available at: http://Basilisk.basilisk.fr/

20.

Popinet

Rickard

(2007) A tree-based solver for adaptive ocean modelling. Ocean Modelling 16(3–4): 224–249.

21.

Power

(2013) Review of tsunami hazard in New Zealand (2013 update). GNS Science Consultancy Report 2013/131. Lower Hutt: GNS Science.

22.

Scardovelli

Zaleski

(1999) Direct numerical simulation of free-surface and interfacial flow. Annual Review of Fluid Mechanics 31: 567–603.

23.

Tappin

Watts

Grilli

(2008) The Papua New Guinea tsunami of 17 July 1998: Anatomy of a catastrophic event. Natural Hazards and Earth System Sciences 8(2): 243–266.

24.

Tappin

Watts

McMurtry

. (2001) The Sissano, Papua New Guinea tsunami of July 1998-offshore evidence on the source mechanism. Marine Geology 175(1–4): 1–23.

25.

Watts

Grilli

Kirby

. (2003) Landslide tsunami case studies using a Boussinesq model and a fully nonlinear tsunami generation model. Natural Hazards and Earth System Sciences 3(5): 391–402.

26.

Watts

Grilli

Tappin

. (2005) Tsunami generation by submarine mass failure. II: Predictive equations and case studies. Journal of Waterway Port Coastal and Ocean Engineering-ASCE 131(6): 298–310.