Articles | Volume 22, issue 2
Research article
 | Highlight paper
25 Feb 2022
Research article | Highlight paper |  | 25 Feb 2022

Multilayer modelling of waves generated by explosive subaqueous volcanism

Matthew W. Hayward, Colin N. Whittaker, Emily M. Lane, William L. Power, Stéphane Popinet, and James D. L. White

Theoretical source models of underwater explosions are often applied in studying tsunami hazards associated with subaqueous volcanism; however, their use in numerical codes based on the shallow water equations can neglect the significant dispersion of the generated wavefield. A non-hydrostatic multilayer method is validated against a laboratory-scale experiment of wave generation from instantaneous disturbances and at field-scale subaqueous explosions at Mono Lake, California, utilising the relevant theoretical models. The numerical method accurately reproduces the range of observed wave characteristics for positive disturbances and suggests a relationship of extended initial troughs for negative disturbances at low-dispersivity and high-non-linearity parameters. Satisfactory amplitudes and phase velocities within the initial wave group are found using underwater explosion models at Mono Lake. The scheme is then applied to modelling tsunamis generated by volcanic explosions at Lake Taupō, New Zealand, for a magnitude representing an ejecta volume of 0.1 km3. Waves reach all shores within 15 min with maximum incident crest amplitudes around 0.2 m at shores near the source. This work shows that the multilayer scheme used is computationally efficient and able to capture a wide range of wave characteristics, including dispersive effects, which is necessary when investigating subaqueous explosions. This research therefore provides the foundation for future studies involving a rigorous probabilistic hazard assessment to quantify the risks and relative significance of this tsunami source mechanism.

1 Introduction

Subaqueous eruptions are poorly understood volcanic phenomena that can generate hazardous tsunamis. When fully submerged, eruptions can transfer energy into generating impulsive waves by the displacement of water, either by regular or periodic jetting, flank instability, or explosion (Duffy1992; Egorov2007; Paris et al.2014). These processes can expand the hazard footprint of an eruption far beyond the classical primary hazards, posing danger along the shores of caldera lakes and coastal areas near explosive volcanoes. While volcanoes are estimated to be responsible for just 5 % of all noted tsunamis since 1600 CE, they can be particularly dangerous in that they account for 20 %–25 % of all recorded fatalities resulting from volcanic activity (Mastin and Witter2000). However, tsunamis are often neglected from hazard maps of volcanoes. Recent events such as at Anak Krakatau in late 2018 have underlined the need to consider them in any disaster risk response (Grilli et al.2019; Williams et al.2019; Ye et al.2020). This is especially pertinent for communities that are not exposed to or are less familiar with seismogenic tsunamis and may not consider themselves at risk.

While experimental data and detailed field studies of these phenomena are rare, some reliable observations exist. For example, at Karymskoye Lake, 1996, a Surtseyan-style eruption was partially witnessed from the air including six explosions followed by tsunamis and base surges. Later ground investigation revealed run-up along the lake ranging from 19–1.8 m at distances 0.5–3 km from the vent, debris flows down the Karymskaya River, and boulder transportation up to 60 m inland (Belousov et al.2000; Torsvik et al.2010; Ulvrová et al.2014; Falvard et al.2018). The Ritter Island volcano generated one of the largest known tsunamigenic flank collapses in 1888, leaving only a small remnant above the water surface, and has since experienced occasional submarine eruptive activity and small local tsunamis in 1972, 1974, and 2007 (Johnson1987; Dondin et al.2012). In the Caribbean near Grenada, Kick'em Jenny volcano was discovered mid-eruption in 1939 and has been regularly active and progressively shoaling, with eruption columns breaching the surface in 1939 and 1965 generating minor tsunami waves (Smith and Shepherd1993, 1996). Many other candidate eruptions are historically documented with small amplitude waves such as Kavachi, Solomon Islands, or lack detailed proximal observations, which leaves uncertainty as to the source mechanism responsible, for example, the 1952 Myojin-Sho submarine eruption which destroyed a naval research vessel or the 1883 eruption of Krakatau (Dietz and Sheehy1954; Nomanbhoy and Satake1995).

Explosive volcanic eruptions are characterised by a directional gas-driven jet from the source, release of water vapour, and, in subaqueous settings, potentially violent vaporisation of sea or lake water on interaction with hot magma. This violent vaporisation is a characteristic of phreatomagmatic eruptions, and it leads to rapid expansion of the resultant water vapour at depth leading to disturbance of the water surface and propagation of waves. To investigate the potential hazard range, we need to understand the relationship between source parameters of the eruption and the nature of waves they generate.

Underwater explosions are well documented (Cole1948; Le Méhauté1971; Mirchina and Pelinovsky1988; Kedrinskiy2006; Egorov2007) primarily owing to military reports and research in blast mitigation and structural response in, for example, ship hulls and other coastal or off-shore structures (Klaseboer et al.2005; Aman et al.2012; Liu et al.2018). As a result, significant research efforts have usually been focused on non-linear fluid–structure interactions such as pressure loading from shock waves rather than any wave generation relationships. Still, some tests were conducted on this matter during the nuclear-testing age and led to the development of theoretical models describing explosion–surface interaction and dynamics of the resultant wave field (Le Méhauté and Wang1996). Physical experimentation since the end of nuclear testing has been rare due to cost, practicalities, environmental concerns, and the challenges of scale experienced by previous tests. In their place, numerical investigations are now the predominant area of research and offer the most likely route to advance the understanding of these processes.

The current theoretical models summarised by Le Méhauté and Wang (1996) have been used in recent years to simulate the wavefield generated from events that produce analogous water surface cavitation such as subaqueous volcanic explosions (Torsvik et al.2010; Ulvrová et al.2014; Paris et al.2019; Paris and Ulvrová2019) and asteroids impacting in ocean basins (Ward and Asphaug2000). However, numerical solutions often either utilise the empirically derived relations without validating their use in a numerical scheme against a suitable explosive physical experiment or test a generation mechanism in the local spatial range only at the cost of neglecting investigation of the generated wave field. Often, models such as those based on non-linear shallow water equations are applied to these problems without considering how dispersive the resultant waves may be (Paris and Ulvrová2019). Non-hydrostatic (NH) multilayer models such as NHWAVE (Ma et al.2012) and SWASH (Zijlema et al.2011) have been developed over the last 15 years to allow a wide range of spatial scales of flows such as surface waves in complex environments to be simulated in the same framework, and they have been applied to dispersive and non-linear tsunamis (Ma et al.2013; Grilli et al.2015; Glasbergen2018; Schambach et al.2019; Ruffini et al.2019; Grilli et al.2019, 2021; Schambach et al.2021).

This work uses a recently developed non-hydrostatic multilayer solver for free-surface flows to model the physical problem. Firstly, the method is validated against a laboratory-scale experiment of released columns of water to ascertain the numerical solution's robustness in resolving a simplified comparable wave generation mechanism. Secondly, data from one of the last military explosive test series focused on surface wave observations is compared with results produced by implementing the theoretical model's initial conditions in the numerical method. These tests are to establish the fitness of the underlying models, which are then applied to a hypothetical explosive subaqueous eruption at Lake Taupō, New Zealand, to provide an example of how these models can be used.

2 Methodology

2.1 Underwater eruption model

An explosive subaqueous volcanic eruption is a dynamic and complex event involving abrupt fragmentation, volume change, and numerous high-energy interactions between pressurised magma, volatiles, and water. Their wave generation capability depends on numerous physical parameters, including eruptive energy, depth, duration, and vent geometry (Egorov2007; Paris2015). The scarce availability of field observations combined with practical limitations both in field and laboratory necessitates simplifications to be made for an explosive eruption model such as considering it as a point-source explosion, as proposed and utilised by Torsvik et al. (2010), Ulvrová et al. (2014), and Paris and Ulvrová (2019).

The models developed for subaqueous explosions and their waves are derived from experimental data and visual observations from chemical and nuclear explosive testing during the 20th century. As documented at the time, water disturbances are born from the generation and rapid expansion of a gas bubble that interacts with the free-surface by collapsing into a crater-like cavity, accompanied by central jets of water and the initial dissipative cylindrical bore, which radially expands outward. The resultant cavity rapidly fills under gravity to produce the second, larger jet which produces a further cylindrical bore, after which the disturbance oscillates until rest, precipitating waves of decreasing amplitude. This free-surface interaction is strongly linked with the depth of explosion relative to its energy; small-yield or deep detonations lead the explosive bubble to transfer a large portion of its energy to the surrounding water through rapid oscillations, which significantly reduces wave-making efficiency. (Le Méhauté1971; Le Méhauté and Wang1996)

Bubble dynamics is a very active area of research in computational fluid dynamics (CFD), although, in the explosive realm, the focus is usually on pressure waves and solid interactions (Wang et al.2018). These studies are usually short in temporal range and are very computationally expensive as modelling the full problem requires accounting for compressibility and multi-phase flow; thus, researchers have generated specialised codes as a solution (Hallquist1994; Li et al.2018). Only in recent years have studies directly simulated expanding explosive bubbles interacting with deformable beds and a free-surface (Petrov and Schmidt2015; Daramizadeh and Ansari2015; Xu et al.2020). However, there is minimal focus on subsequent surface waves, let alone relations tested for their generation mechanisms or far-field propagation. Recently, physical experiments that utilise underwater gas jets by Shen et al. (2021a, b, c) quantify wave generation regimes and establish relationships between jet intensity, duration, depth, and maximum wave height, showing that for a given jet intensity, there is a critical depth and minimum duration limit that produces maximum wave heights.

Following development of the physical theory of underwater explosions, mathematical models were developed by applying inverse methods to experimental time series and simplifying the result to a two-parameter model corresponding to initial conditions on the free-surface elevation (η0) representing the maximum surface displacement from the explosive disturbance (Le Méhauté and Wang1996). ηc corresponds to the maximum depth of the disturbance below equilibrium and R to its radius extent. These parameters physically represent the size of the initial cavity and are functions of explosive yield E, water depth h, burst depth z, and bed characteristics for which calibration is made with empirical data. The initial disturbance can therefore be described analytically by one of a number of candidates for the general profile as Eqs. (1)–(2) (Le Méhauté and Wang1996). Note that Eq. (1) is discontinuous at its edge, while Eq. (2) returns back to zero. A schematic diagram illustrates the problem and initial profiles in Fig. 1a–c.


Figure 1Illustrations of the subaqueous explosion problem. (a) An explosion of yield E at depth z in water of depth h. (b) The initial surface elevation profile is given by Eq. (1). (c) The initial surface elevation profile is given by Eq. (2). (d) Schematic of a volcanic scenario in which such an explosion would occur at maximum depth where z=h, and crater diameter CD can be measured or calculated using estimated ejecta volume V with Eqs. (9)–(10).


The relations between parameters ηc, R, and explosive characteristics described here are derived empirically after many series of small- and larger-scale experimental observations and are well described and reviewed by Le Méhauté and Wang (1996). These relations depend on classifications of water depth h and charge depth below water surface z relative to explosive energy released E. In terms of a depth parameter,

(3) D = c h E 3 ,

where c=406.2 is an imperial unit conversion constant, three categories are specified when considering wave generation:

(4) Classification = Deep , D > 14 Intermediate , 1 < D 14 Shallow , D 1 .

For deep and intermediate cases, the cavity parameters are defined as follows:


where constants a and b vary as follows:


For shallow cases, it is implicitly assumed that the explosion develops a cavity that extends through the entire water column and exposes the bed, and therefore its radius is larger than the water depth (R>h). In this instance, cavity radius is defined as follows:

(9) R = 0.03608 E 1 4 .

The data calibrating these models include charges ranging from small (<500 lb of TNT or <9.5×108 J) to a handful larger (<9500 lb of TNT or <1.8×1010 J) and further include a 23 kT nuclear test (Le Méhauté and Wang1996).

Volcanic context

For a volcanic case, illustrated in Fig. 1d, such explosions would occur on or in an edifice, meaning that the charge depth is equivalent to the water depth at that point (z=d); therefore events that are capable of hazardous wave generation fit into the intermediate or shallow depth classes.

An approach first demonstrated by Torsvik et al. (2010) and expanded on by Ulvrová et al. (2014) and Paris and Ulvrová (2019) identifies the use of the parabolic initial deformation (Eq. 1) and its compatibility with the shallow class radius definition (Eq. 9) when paired with the intermediate condition of Eqs. (5) and (7). This is valid when the released explosion energy is estimated from volcanic crater diameter CD using the following empirical relationship by Sato and Taniguchi (1997):

(10) E = 4.45 × 10 6 C D 3.05 .

This method has recently been used for probabilistic hazard analysis of volcanogenic tsunamis at the Campi Flegrei caldera, Italy (Paris et al.2019), and at Taal Lake, Philippines (Pakoksung et al.2021). Some land studies suggest that the size of a volcanic crater cannot be assumed to directly reflect the size of the largest explosion causing it (Valentine and White2012), so a further relation from Sato and Taniguchi (1997) relates the ejecta volume V with the released explosion energy:

(11) E = 4.055 × 10 6 V 1.1 .

2.2 Numerical method

To compute simulations of the models described earlier, we use the open-source CFD framework Basilisk (Popinet2013). The software is widely used in studies involving multi-phase problems from jet dynamics to viscoelastic and surface tension investigations and includes several free-surface solvers with application to tsunamis, wavefield transformation, and other hydrological areas (e.g. Lane et al.2017; López-Herrera et al.2019; Berny et al.2020). The benefits of the Basilisk framework include OpenMP/MPI parallelism capability and an adaptive-tree-based Cartesian grid with refinement algorithms that accept user-defined criteria. These allow the appropriate resolution of areas of irregular shape or distant from one another within a domain with high-scale separation while increasing the overall efficiency and speed of computation. (Popinet2011; van Hooft et al.2018)

2.2.1 Multilayer Scheme

The majority of the following work utilised the non-hydrostatic multilayer free-surface solver developed and described by Popinet (2020). A brief outline is given here. The scheme of n layers is a horizontally gridded and vertically discrete approximation of the incompressible Euler equations with a free-surface and gravity. It is described by the following system:


where, in the x-z reference frame, k is the layer index, hk the layer thickness, g the gravitational acceleration, uk and wk the horizontal and vertical velocity components, ϕk the non-hydrostatic pressure, η the free-surface height (sum of layer thicknesses and bathymetry height zb), and

(16) z k + 1 / 2 z b + l = 0 k h l

the height of layer interfaces.

Between them, the set expresses the evolution of layer thickness (Eq. 12), conservation of momentum (Eqs. 13, 14), and conservation of volume and mass (Eq. 15). The framework allows the model to be built modularly, starting from the hydrostatic case where ϕ=0 and vertical momentum conservation (Eq. 14) is removed; these effectively become the generalised multilayer Saint-Venant (SV) or stacked shallow water equations. More components are added, for example, vertical remapping, adaptivity, non-hydrostatic and Keller box vertical projections, and a wave breaking method. The latter is implemented by introducing dissipation. This is handled by limiting the maximum vertical velocity by setting

(17) w k n + 1 sgn w k n + 1 min w k n + 1 , b g H ,

where g|H| is the characteristic horizontal velocity scale of the wave, b is a specified breaking parameter smaller than 1, and sgn and min are sign and minimum functions, respectively. The breaking algorithm used throughout the multilayer scheme is described in greater depth in the model's defining paper by Popinet (2020). Lastly, terrain is handled by looking up a pre-processed k-dimensional-tree-indexed database of heights, and the model is able to discriminate and resolve areas of wetting and drying.

The scheme has been tested against numerous benchmark cases, including standing waves, sinusoidal wave propagation over a bar, the Tohoku tsunami of 2011 and its dispersive features, viscous hydraulic jumps, and breaking Stokes waves (Popinet2020). Source code of these examples can be found at (last access: : 17 February 2022). Able to provide accurate solutions for surface gravity waves across a wide range of scales, the multilayer scheme is ideal for resolving the evolution of an explosive source initial condition, resulting wave propagation, and shore run-up (Popinet2020).

2.2.2 Other schemes

Accompanying models used for validation and comparative purposes include a Navier–Stokes volume-of-fluid (N-S/VOF) method in the same framework which solves the two-phase Navier–Stokes equations for interfacial flows, including variable density and surface tension (Popinet2009, 2018), a solver for the shallow water or Saint-Venant (SV) equations, and finally another for the Serre–Green–Naghdi equations (SGN), a Boussinesq higher-order approximation for non-linear and weakly dispersive flows (Popinet2015). Their inclusion is to support and inform the evaluation of the multilayer scheme against well-known and commonly used methods.

In this work's context, the main discriminations between the schemes used are both the hydrostatic assumptions involved and their resolution of vertical gradients such as the velocity profile: in the hydrostatic SV solver, a constant velocity profile is assumed; the non-hydrostatic SGN equations represent a Boussinesq-type analytical approximation of the vertical structure; the multilayer method resolves the vertical to n layers with the capability to include non-hydrostatic terms; and the N-S/VOF scheme fully resolves the vertical as an additional dimension. The addition of having to fully solve a 2D slice through the vertical versus 1D means the N-S/VOF has a larger domain to compute for the same region when considering water waves.

Only recently have non-hydrostatic multilayer methods advanced towards a point where their stability and efficiency has approached or challenged Boussinesq-type schemes. These are typically used for applications such as landslides which can generate highly dispersive (high kh) initial waves and significant vertical accelerations. A preferred current method is often to initialise the wave generating source, usually a short duration compared to the total tsunami simulation length, using either a direct numerical method (e.g. Abadie et al.2012) or a NH multilayer scheme (e.g. Grilli et al.2012, 2019; Schambach et al.2021; López-Venegas et al.2015) and then to interpolate or otherwise insert the results into an initialisation of a Boussinesq-type scheme to further propagate the waves across a larger domain.

3 Laboratory-scale validation

To determine the suitability of the numerical method for use in modelling a subaqueous disturbance, validation against a suitable case study is required. Prins (1958) conducted a flume experiment investigating surface waves produced from an instantaneously raised or depressed column of water. This case is replicated using the multilayer scheme and, additionally, the N-S/VOF, Saint-Venant, and SGN schemes for full numerical comparison. Unless otherwise specified, all length units in this section only are in imperial standard for consistency and ease of comparison with the Prins (1958) dataset (to approximately convert to metric units, from feet to metres, divide by 3.281).

Figure 2(a) Schematic diagram of the experimental set-up of Prins (1958), with gauge locations and the initial disturbance. (b) ηQ relation for first arrival crest or trough amplitude at x=5 ft. (c) ηQ relation for first arrival crest or trough amplitude at x=25 ft. Regression lines through origin are plotted for each L-model combination. For approximate metric conversion of feet to metres, divide by 3.281.


The experimental set-up involved a column of water at rest of extent L, height from equilibrium Q and in a 1 ft wide flume of water depth h; these parameters were varied by steps: L=2, 1, and 1/3 ft; Q=±0.1, ±0.2, ±0.3 ft; h=2.3, 0.5, 0.35, 0.2 ft. Figure 2a illustrates a schematic diagram of the set-up. The height-adjusted column was held by a difference in air pressure and released by means of a sliding door opened by hand in under 0.03 s. The resulting waves are measured at five gauges positioned every 10 ft starting at 5 ft. A sloped “wave absorber” dampened reflections from the flume end.

The numerical models to replicate this are relatively simple. For the multilayer, SGN, and SV schemes, a 1D domain of length 100 ft is initialised with the origin translated by L, the positive part of the domain with initial water depth η0=h, and the negative with η0=h+Q. The displaced column is released instantaneously. While a reflective boundary condition is applied to both ends, extending the positive domain far beyond the experimental case to prevent reflections has negligible computational cost. In the multilayer simulations, 20 regularly spaced layers were used, and the breaking parameter, following guidance and methods used in Popinet (2020), is set to b=0.38. The N-S/VOF scheme's geometry is built similarly but in 2D vertical (2DV) such that the origin is translated by (-L,-h). Material properties (i.e. density, dynamic viscosity, surface tension coefficient) of the two phases are also specified, as is gravitational acceleration. For all presented cases, unless otherwise specified, the maximum level in the grid is 14, which, for a domain of length 100 ft, gives a maximum horizontal mesh resolution of 6.1×10-3 ft.

3.1 Model comparison

Figures 2b–c and 3 present direct comparisons between the numerical models and the published experimental data. The N-S/VOF and multilayer approaches each conform well when varying L and Q at h=2.3 ft and generally confirm the experimental findings for initial wave amplitude throughout the flume range. Some experimental data points at high Q are missing; however, the model output highlights a slight divergence in first arrival amplitudes positive (higher) and negative (lower) disturbances at larger Q.

Figure 3Model comparison of time series of waves produced by a disturbance of Q=0.3 ft, and L=13 ft when h=0.5 ft. For approximate metric conversion of feet to metres, divide by 3.281.


Figure 3 plots the full time series for a mid-range parameter run and includes solutions from the non-linear shallow water (SV) and Boussinesq-type (SGN) schemes. While the SV solution exhibits a good approximation of the phase velocity of the initial wave (4.39 ft s−1), it is, unsurprisingly, unable to resolve any of the trailing wave field. The SGN scheme resolves this element at far-field well; however, it often overestimates initial wave amplitude and the overall dispersivity. The multilayer scheme is shown to be very accurate compared to the N-S/VOF result and the experimental trace and reinforces the good fit found in Fig. 2. This suggests that the scheme faithfully replicates the process of a collapsing water column or infilling of a uniform depth and the associated wave generation, along with reflections from behind the initial disturbance.

Figure 4Comparison of numerical solutions for the disturbance collapse from initialisation to t=0.5 s. Distances are normalised with water depth h and are vertically exaggerated.


The generation process is shown in Fig. 4, which illustrates how the four models resolve an initial disturbance and reveals significant variations in handling the vertically critical jump in the first half-second. The SV source quickly develops into the classic steep-fronted crest often seen in dam-break problems in which the amplitude is proportional to the initial disturbance height. The higher-order SGN scheme is able to resolve frequency dispersion and therefore handles the source development better; however, there are clear exaggerations such as at 0.05 s near x=-L where, intriguingly, the water height temporarily increases above Q. Qualitative comparison between the development of both the N-S/VOF and multilayer solutions shows a high consistency across all stages involving the collapse of initial disturbance and breaking of the generated wave. Notable variation includes the varying maximum amplitude of the resultant wave's intermediate stage and any bubbles/droplets produced, demonstrating the restriction of a single value for function on the 1D multilayer scheme not present on a 2D multi-phase VOF solver, meaning that there are many neglected phenomena in the multilayer scheme such as bubbles and plunging breaks that, while this may be negligible in this lab-scale experiment, could be more significant at larger scale. Also, note graphical artefacts on free-surface height at the leading gradient in the N-S/VOF solution which occur at coarser regions of the adaptive grid. Despite these minor variations, it is clear that the multilayer model has greater validity in application to this case than either of the single layer models and is remarkably consistent with the physical experiment, as well as the directly simulated approach.

Finally, Table 1 presents performance metrics for the numerical schemes across three maximum refinement levels for an example run case. All were performed with OpenMP parallelism on eight CPU cores. The multilayer scheme fits between the SV and SGN methods in terms of wall time processing and offers a vast improvement in computational efficiency compared with the N-S/VOF solver considering result similarity. It is also faster than the dispersive SGN method primarily due to the computation time required to solve the higher-order approximation.

Table 1Table of numerical model performance for run h=2.3 ft, L=0.2 ft, and Q=0.3 ft (h=0.701 m, L=0.061 m, and Q=0.0914 m). Runtime in seconds and speed in points × time steps × layers per second per core.

Download Print Version | Download XLSX

3.2 Wavefield classification

All simulated cases are plotted in Fig. 5 by the size of initial disturbance relative to water depth using a non-linearity parameter |Q|h against a dispersivity parameter kh, where k=2π2L. Six additional simulations were run outside of the original experimental parameter space to expand the model dataset reach. As done by Prins (1958), the +Q runs are categorised into groups with similar wavefield characteristics starting with strong oscillatory properties (blue) where kh>10|Q|h, tending through increasingly solitary wave properties (purple and green) once kh 10|Q|h until a succession of diminishing amplitude solitary waves result where kh<|Q|h (orange). Beyond this region, the initially generated bore survives far enough down the numerical flume before it would likely separate (black). For the Q domain, all resultant wave fields were similar except for the length of initial trough relative to the following periodic wave group. This ratio initially remains approximately unity in the same region as the +Q oscillatory character group and grows larger towards higher |Q|h and lower kh.

Figure 5Plots of numerical runs by dispersivity parameter against non-linearity parameter. All data from time series are at x=35 ft. (a) Classification of +Q disturbances by their generated wave field. The legend underneath is illustrated by plots of example runs. (b) Q disturbances coloured by the ratio of initial trough length divided by the following wavelength.


Results of the +Q disturbance (Fig. 5a) corroborate the experimental wave field descriptions of Prins (1958), including the transition of an oscillatory field through to solitary initial waves. Bore formations in the first stages were also observed during the experiment, and those which last a considerable length of the flume match model results. Accounting for tolerance in qualitative descriptions, groups match closely, and the additional results beyond the experimental scope further confirm these definitions in the studied range. Such an analysis of the Q part was not attempted in the original research; however, effort in this area can be made with the numerical results (Fig. 5b). The length ratio of initial trough to the following oscillatory waves increases with higher Q/h and lower kh. This matches the trend towards solitary characteristics with +Q. Intriguingly, this pattern holds regardless of the length ratio that defines the initial disturbance (i.e. Q/L).

These results bear similarity to those found by Hammack and Segur (1974, 1978a, b) in experiments involving a piston producing vertical bottom motions described by Hammack (1973) and modelling using the Kortewig–de Vries equation, also with an initialised rectangular wave source. Notable similitudes include the generation of potentially multiple solitons of decreasing amplitude for shallow water positive initialisations (as seen in orange region of Fig. 5a), followed by a train of dispersive oscillatory waves, and that no solitons are generated from negative vertical motions, instead producing a wave train of the type illustrated in Fig. 5b of an initial “triangular” wave of greater speed than the trailing modulated oscillatory waves.

In suitable replication of the experimental findings, the present numerical scheme is seen to be fit for generating accurate waves from initialised disturbances and modelling their near-field propagation across a significant regime range where non-linear and dispersive effects may be prevalent. The method also demonstrates suitability for further investigations either beyond or in complement, such as to widen parameter spaces with relatively low computational expense.

4 Field-scale validation

The next stage is to assess the use of the underwater explosion models of Sect. 2.1 within the numerical scheme. To do this, we utilise datasets from the Mono Lake test series in 1965, conducted by the Waterways Experiment Station, and documented by Walter (1966), Wallace and Baird (1968), Whalin et al. (1970), and Pinkston et al. (1970). This was one of the largest chemical explosive test series designed to investigate subsequent water wave generation and shore effects. A series of 10 approximately 9250 lb ( 4196 kg) spherical TNT charges were detonated off the south shore of Mono Lake, California. The test area is illustrated in Fig. 6a, and this also shows the location of wave gauges arranged in four radials directed away from ground zero (GZ), along with contoured terrain. GZ was located at a site of approximate water depth h=39 m.

Figure 6(a) Setting of test series on the south of Mono Lake with bathymetry and near-shore terrain elevation. The yellow star marks ground zero (GZ), and crosses mark gauge locations. Gauge radials are coloured by red (one), purple (two), blue (three), and white (four). (b) Numerical results of maximum crest amplitude for Shot 3 are overlaid by observed maximum crest amplitudes at gauge locations.


The numerical model was built with the multilayer scheme, including a digital elevation model of the Mono Lake bathymetry by Raumann et al. (2002), in which the lake water level was set at elevation above mean sea level 1945.7 m. This left approximately 2 m vertically of the bathymetric model dry to act as the shore surrounding the lake. Numerical gauges were defined at distances from GZ described in the experiment. The maximum horizontal resolution across the domain was 1.46 m, and simulations used five layers.

Table 2Numerical values used in Mono Lake models.

Download Print Version | Download XLSX

Data from Shots 3 and 9, individual explosive tests within the experimental series, are chosen to test the numerical method against. These experiments resulted in the highest-quality data and were detonated at different depths: z3=0.427 m and z9=6.401 m. Following a similar method used by Paris et al. (2019), the initialised disturbance takes the form of Eq. (2), and, for these charges, D=6.1; therefore the depth classification (Eq. 4) is intermediate. The resultant values using Eqs. (4)–(7) are given in Table 2. Considering doubts regarding the charge magnitudes raised in the report of Walter (1966), additional lower-yield simulations were added.

4.1 Generated wavefield

Figure 6b illustrates the generated field of maximum crest amplitudes from the numerical simulation paired with data transcribed from the experimental gauge records for Shot 3. This and all other included simulations are presented quantitatively in Fig. 7 to show maximum crest amplitudes for all gauges. The maximum amplitude of the initial envelope decreases with radial distance from GZ as expected, matching experimental observations. Shoaling is most pronounced at the closest incident shoreline, in the region 500 m east of Radial 1, and becomes far less significant with distance as seen on the western shore of the region. Maximum crest amplitudes at gauge locations follow similar patterns in all runs; however, the higher-energy-yield simulations produce greater amplitudes (additional 0.11–0.05 m) throughout and experience more significant shoaling in all shallow areas. The lower-yield simulations are a closer match to experimental observations for both shots, especially in shallow zones; however, the experimental data have noticeably greater variation at shore.

Figure 7Time series comparison at closest gauge to GZ for both considered shots.


Figure 8Maximum crest amplitudes at gauge locations for full yield model, fitted model, and experimental tests.


The experimental gauge time series at locations closest to GZ are plotted alongside the lower-yield model traces in Fig. 8 and are useful to compare the phase arrival times and the initial development of the wave train. For both simulations, the first arriving phases match the experimental record very well, with the exception of minor noise immediately following detonation in Shot 9 which is likely to be shock or debris related. The time of maximum envelope amplitude also conforms well, in which the differences are 9 and 7 s for Shots 3 and 9 respectively. The latter part of the initial wave group maintains higher amplitudes in the experimental trace for both records, whereas the envelope decay is sooner in the numerical model. Shot 3 also seems to exhibit a positive amplitude shift in the early part of the experimental envelope. These could be a genuine underestimation of wave amplitudes towards the end of the initial group which could be due to variations in dissipation by breaking waves of the initially generated energy.

In terms of parameters pertaining to dispersion and non-linearity, kh and ka, waves in the initial group near the source at the nearest gauges on each radial were in the ranges 1.34<kh<3.56 and 3×10-3<ka<1.181, and those across the gauges beside the shore were in the ranges 0.105<kh<0.235 and 3.6×10-3<ka<1.2×10-2. As would be expected, moderately non-linear waves are generated, and kh decreases as the waves approach the shore and become shallow. Wave steepness ka immediately decreases when propagating away from source across near-flat bathymetry but then increases on average as water depth h decreases towards shore.

4.2 Model implications

While the models used fit the experimental results and trends well overall, it is significant that the lower-yield data are a much better fit. In reports following the test series (e.g. Whalin et al.1970) it is noted that, except for two charges, all shots delivered below-average expected maximum crest amplitudes as predicted by earlier experimentation, particularly the deep water shots. Further similar deep tests at Mono Lake in the following year, which were conducted for other investigations, delivered much greater amplitude waves in line with expectations, leading to the suggestion by Wallace and Baird (1968) that, beyond scaling effects or measurement issues, the charges used in the series may have been faulty and delivered a lower yield. This is supported by the numerical data as, when energy is reduced, the maximum crest amplitudes can be accurately predicted in addition to the resultant early wave group and individual phases. Moreover, further energy dissipation not accounted for in the physical model may be responsible for the greater fit of the reduced yield simulation, such as losses from a higher amount of dissipation from breaking of the initial waves from the explosion or some of the energy transferred to the nearby bed as elastic deformation. While the initialisation model is calibrated to charge depths relative to water depth, bed characteristics were not strongly considered.

Many data from the experimental series were unreported or discarded following the series due to various problems, including excess noise generated from the explosion itself or wind-driven waves, and are thus missing from comparisons in this work. Instrumental issues meant gauges near the shore were often unable to be acceptably calibrated for amplitude. Many experienced positive noise due to wave breaking and bores, responsible for some spikes in measurements; however, these were kept for reliable arrival times and periods (Wallace and Baird1968).

Despite the experimental challenges, the multilayer scheme is shown to be good at resolving the initialised disturbance into a wave field, which is very consistent with that generated by an intermediate depth submarine explosion. With the capability to accurately propagate such a source in the relatively deep (kh≈3) near-field through to the shore, it demonstrates suitability to model such events at these spatial scales and timescales. However, this test did not contain any shallow depth underwater explosion tests, for which there are no case examples at this scale.

5 Taupō scenario

We now apply the multilayer model presented earlier to a hypothetical sublacustrine eruption to demonstrate a case-use example for a potential hazard study. Lake Taupō is New Zealand's largest freshwater lake with an area of approximately 616 km2. It lies in the southern section of the Taupō Volcanic Zone and conceals most of Taupō Volcano and the caldera of Earth's youngest super-eruption (Oruanui) at ca. 25.5 ka, ejecting >1100 km3 of pyroclastic material (Davy and Caldwell1998; Wilson2001; Vandergoes et al.2013; Allan2013). Taupō Volcano and its post-super-eruption magmatic system have been extensively studied because of how recently the event occurred, along with the volcano's relatively high activity. As a result, the eruptive history of Taupō is well constrained through radiocarbon dating methods, stratigraphy, and composition analysis (Wilson1993; Wilson et al.2009; Barker et al.2015).

Figure 9Lake Taupō setting with lake bathymetry, surrounding terrain, marked fault zones, and built-up areas. Also highlighted are the caldera zones of the Hatepe/Taupō eruption and the Oruanui super-eruption (Leonard et al.2010). The dashed red area encloses area of activity within 5 kyr from Wilson (1993) and Barker et al. (2015). The star indicates simulated eruption location at 38.813 S, 176.016 E. Below is a cross-sectional profile of high-run initial condition with terrain elevation in a north–south direction in line with source location.

The Oruanui event itself was the evacuation of a rapidly formed (<3000 years) and well-mixed high-silica rhyolite-melt-dominant body (Wilson et al.2006; Allan2013; Wilson and Charlier2009; Allan et al.2013) which left the remaining mush source heavily modified by intrusion of hot mafic magmas (Barker et al.2014, 2015). Following a notably brief 5 kyr quiescence, three further, much smaller (0.004 to 0.05 km3 dense rock equivalent, DRE) dacitic events occurred from  20.5–17 ka which preceded another 5 kyr gap before an additional 25 rhyolitic eruptions from  12–1.8 ka, which are all remarkable in their high variation in eruption volumes (0.015 to 35 km3 DRE) and concentration of vent locations, shown in Fig. 9 (Gelman et al.2013; Wilson1993; Sutton et al.2000; Barker et al.2015). Compositions of these eruptions reveal a thermal reset of the magmatic system, with more crustal material incorporated into the younger magmas in the series (Barker et al.2014, 2015). The largest of these, the Taupō Plinian eruption in  232 CE, was one of the most powerful events globally in the past 5000 years (Wilson and Walker1985; Houghton et al.2010) and resulted in the further collapse of the caldera beneath the lake and, afterwards, the formation of the Horomatangi Reefs (Davy and Caldwell1998).

After the last events, the lake underwent refilling and eventually breakout down a main outlet, the Waikato River (Manville et al.1999; Manville2002). In modern times, the lake has impeded monitoring and investigation of the volcanic system state, but this has improved with developments in seismology, geodesy, and other geophysical techniques, including identifying numerous seismic swarms and ground deformation in 1983, 1997, 2008, and 2019 (Barker et al.2020; Illsley-Kemp et al.2021). This, accompanied with a good understanding of the varying eruptive history, informs the current state of Taupō. However, the lack of relationship between repose and a potential eruption size makes it challenging to infer future activity. Therefore, any attempt at scenario-based modelling involving Taupō Volcano should consider plausible cases across the full span of statistical likelihood. Here we model a single eruption from Taupō to demonstrate the multilayer model but highlight that future work will be required to provide an assessment of tsunami hazards from this volcano.

5.1 Model set-up

Figure 9 illustrates the lake bathymetry and surrounding terrain accompanied by geological features and settlement locations. The lake is fed by multiple rivers and notably from Tongariro hydroelectric power station via the Tokaanu Tailrace Canal. The sole outflow is the Waikato River, controlled by gates at the largest settlement on the lake, Taupō, which leads to numerous further hydroelectric dams downstream. Surrounded by abundant geothermal resources, strong trout fishing, and agricultural industries, the area also boasts plentiful tourism opportunities, hosting over 1 million tourists each year.

In building a model representing an example tsunamigenic explosive eruption in Lake Taupō, building a terrain dataset by combining a bathymetric model of the lake was required (Rowe et al.2002) with an elevation model generated using lidar datasets from the Waikato Regional Council which cover the entire foreshore. The limiting resolution of the resultant digital terrain model is that of the bathymetric model (10 m); therefore the simulations were performed at a grid refinement level of 11 resulting in a horizontal resolution of 16 m.

The eruption site chosen is within the region where most Holocene vents are located and in an active geothermal field (Bibby et al.1995; de Ronde et al.2002). Considerations made for the size of eruption involve selecting a size which is plausible with respect to the eruptive history. A single, gigantic blast corresponding to an instantaneous release of material equivalent to a larger eruption is not realistic as eruptions of such size typically endure for hours or days (Pyle2015). Therefore, for this example, a simulation is run with an ejecta volume corresponding to a plausible mass eruption rate (MER) of 1.2×107 kg s−1, which is typical of a volcanic explosivity index (VEI) 4 intensity eruption (Barker et al.2019). Over the simulation time (1000 s), the estimated ejecta volume V is 0.004 km3, and, using Eqs. (9)–(10), the equivalent explosive energy can be obtained. At this depth (h=65 m), all events of V>0.0012 km3 yield D<1, falling into the shallow depth class. Following the rest of the method in Sect. 2.1.2, model parameters for both models are detailed in Table 3. A further simulation of the scenario was performed with the Saint-Venant scheme with identical terrain and model geometry for comparison.

Table 3Numerical values used in Lake Taupō models.

Download Print Version | Download XLSX

5.2 Generated wavefield and shoreline impacts

Figure 10 describes the travel times and maximum crest amplitudes. As demonstrated earlier, a succession of waves radially propagate outwards from the centre of the explosion. The initial phase velocity typically starts from  40 m s−1 until they heavily interact with the nearby Horomatangi Reefs. Due to the limited area of the lake, the entire shoreline experiences wave phenomena from these sources within 15 min. Long wavelength oscillations (up to seiche) linger throughout the lake beyond simulation length (17 min).

Figure 10Multilayer simulation of potential volcanic explosion in Lake Taupō. (a) First arrival travel times in minutes. (b) Field of maximum crest amplitudes. (c) Snapshot of wave amplitudes at t=4.5 min. (d) Equivalent of (c) but using the Saint-Venant scheme.

The initial wave at all locations is a crest, reflecting the positive amplitude lip of the initial condition. Earliest arrivals at shore occur about the closest point east of the source at 4 min, after which arrival time generally scales with radial distance except for areas with extended shallow zones. The highest wave heights incident to the shore are, unsurprisingly, located nearest the event on neighbouring eastern and northern shores where crest amplitudes reach over 0.11 m. The lowest are found in the further area of the south-west beside Gauge 8 (Tokaanu) and in sheltered parts away from direct paths such as by Gauges 10 and 11 (Kinlock). Taupō township is relatively sheltered compared to the surrounding shoreline due to shielding from the lake morphology.

Figure 11Time series for numerical gauges with gauge locations highlighted in red on inset map including ground zero signified by the star. Time series have varying y-axis scales and truncated x axes.

Figure 11 presents numerical gauge time series for both runs. Throughout the domain, the high ejecta volume run returns significantly higher maximum crest amplitudes than the low run. Wave periods vary from  40 s early in the group to  15 s towards its end, 10 min later. The first arrival is generally the longest period wave but rarely the greatest amplitude at the gauge locations. Gauges 4 and 5, positioned in shallow zones near the eruption, initially record initial waves of amplitude 0.017 and 0.04 m before waves of 0.02 and 0.11 m respectively. Wave generated by this source are initially highly dispersive within 5 km, with waves through the group in the ranges 2.17<kh<5.57 and 4.3×10-4<ka<8.5×10-3, and approaching nearby shores (e.g. most numerical gauges shown in Fig. 9) they become more non-linear, exhibiting values in the ranges 0.03<kh<1.35 and 1.8×10-4<ka<5.7×10-3.

The resultant waveforms produced from this source are similar to those in Sect. 3 and indicate that the multilayer scheme is suitable for modelling their propagation. The combinations and transformations between highly dispersive, non-linear, and shallow waves across the whole domain require treatment that would likely lack validity if only modelled using the shallow water equations. A snapshot of this comparison is illustrated in Fig. 10c–d and clearly shows the heavy dispersion modelled in the multilayer scheme which is missing in the Saint-Venant method. The benefit of non-hydrostatic pressure terms and numerical wave breaking therefore adds robustness to the near-shore solution for this type of source. Notably, the explosion model could prescribe an initial condition that intersects the bathymetry in a higher ejecta volume case, which would result in additional mass added via the volume of lip surrounding the cavity. It would be expected that a higher-energy explosion in similarly shallow water depth would transmit less energy into the water and towards wave-making. To rectify mass imbalance, the lip height could be lowered to better match the excavated volume.

5.3 Model implications

While re-emphasising this is a model case-use example, these preliminary results suggest some potential implications for modelling wave hazard from explosive subaqueous eruptions. If an eruption of sufficient magnitude at Lake Taupō produces an initial explosion, there could be a threat posed to nearby shores. As with most lacustrine tsunami hazards, there is minimal time from source to shore impact; no existing warning system would ever be able to respond to an eruption with sufficient speed.

The underlying caldera has frequently experienced minor unrest (Potter et al.2015), and current thought suggests eruption probabilities in the near-term are not negligible; for instance, an event of magnitude above that considered here (0.1 km3) is estimated to have 5 % probability within 100 years (Bebbington2020). Taupō Volcano can produce far greater magnitude events than considered here, including the aforementioned 232 CE Taupō eruption (35 km3 DRE) and the ca. 25.5 ka Oruanui super-eruption (530 km3 DRE). Events of such magnitude, even when considering they may consist of multiple smaller episodes, undoubtedly carry wide-ranging hazards well beyond the lake's proximity, and even those towards the lower end of the scale could produce numerous volcanic dangers, including ashfall and pyroclastic density currents. Therefore, an additional complexity of modelling the suite of hazards posed by subaqueous volcanism will be to determine the relative weight each source component possesses with events of varying location and magnitude.

Only a brief effort is made at present on modelling this hazard as an example application of the multilayer numerical scheme and its benefits on resolving the resultant wavefield and significant outputs. Further work should incorporate this or similar numerical methods into conditional probabilistic hazard models for assessing the relative significance of this tsunami source mechanism, using wide parameter space comprising all likely sublacustrine eruption locations and magnitudes. Such a model would be able to take advantage of this scheme's broad wave regime validity and computational efficiency and potentially be able to investigate inundation in detail at various points onshore with, for instance, building data, key infrastructure impacts (e.g. control gates at outlet to the Waikato River), or similar additional model layers.

6 Conclusions

The non-hydrostatic multilayer scheme used in this paper has been shown to accurately replicate the collapse of various initial disturbances into a resultant wavefield that exhibits varying degrees of non-linear properties and frequency dispersion. By capturing some depth aspects of the model without fully resolving the vertical, it is superior in accuracy to shallow water equation-based schemes while being far more computationally efficient than direct numerical methods. The method was used to verify experimental results of positive amplitude disturbance wave generation and probe their validity at a wider parameter range while also further investigating negative amplitude disturbances, revealing the extension of the leading trough relative to trailing oscillations for larger size disturbances and smaller water depths.

Initialisations of wave generation via underwater explosion were tested for use with the multilayer scheme by simulating detonations of explosives as done in a US Army test series at Mono Lake. Including consideration of uncertainty with experimental data, the combination of empirically derived underwater explosion model and the numerical scheme was able to capture the significant elements of the generated waves as measured experimentally and help validate the use of the underlying empirical relations.

The multilayer scheme was then used to simulate an example volcanic explosion under Lake Taupō based on estimated eruptive energy by ejecta volume and MER. Implications for the developed model are suggested in the context of a small-magnitude phreatomagmatic eruption, with data revealing areas of varying exposure to waves – of above 0.1 m at locations near the source – and waves reaching throughout the lake within 15 min. A probabilistic investigation is likely required to assess the full range of possible scenarios at this location, including eruption geometry and size, while potentially considering further complexity, such as any syneruptive variation in initial conditions. In addition, any hazard model could usefully incorporate the numerical outputs in investigating cascading hazards such as one involving wave impact on the inlet dam at the Waikato River. Such a modelling effort would help resolve the significance of this hazard source compared to alternative tsunamigenic sources and volcanic hazards across varying magnitudes of eruption. Further work could include numerical investigation of non-detonation eruptions involving sustained jetting.

Code and data availability

The numerical model framework code is free software and available at (last access: 17 February 2022). The latest version date used in this work is 12 November 2021. Code used in Sects. 3 and 4 are available at (Hayward2022), and code used in Sect. 5 is available on demand. Simulation data for all sections are available on demand.

Video supplement

Animation of the free-surface elevation of the Lake Taupō multilayer run is available from the TIB AV-Portal at (Hayward2021) and for high quality on request.

Author contributions

MWH undertook manuscript preparation, model designing, validation, and running, background research, conceptualisation, and formal analysis. CNW, EML, and WLP designed the research concept, were responsible for funding acquisition, provided supervision, and contributed to revisions and editing of the manuscript. SP and JDLW contributed to reviewing and revisions of the manuscript.

Competing interests

The contact author has declared that neither they nor their co-authors have any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


We are grateful for the comments and suggestions of three reviewers, Stephan Grilli, Simon Barker, and an anonymous reviewer, which have markedly improved this work. We wish to thank Sanjay Wadhwa for use of the bathymetric model of Lake Taupō courtesy of MBIE/NIWA contract C01X1005 (Cumulative Effects of Stressors on Aquatic Ecosystems) and Waikato Regional Council for lidar datasets of the Lake Taupō foreshore and surrounding topography. We acknowledge the use of New Zealand eScience Infrastructure (NeSI) high-performance computing facilities as part of this research.

Financial support

This research has been supported by the Marsden Fund Council, Royal Society Te Apārangi (grant no. 17-NIW-017), awarded to NIWA.

Review statement

This paper was edited by Animesh Gain and reviewed by Simon J. Barker, Stephan Grilli, and one anonymous referee.


Abadie, S. M., Harris, J. C., Grilli, S. T., and Fabre, R.: Numerical modeling of tsunami waves generated by the flank collapse of the Cumbre Vieja Volcano (La Palma, Canary Islands): Tsunami source and near field effects, J. Geophys. Res.-Oceans, 117, C05030,, 2012. a

Allan, A. S. R.: The Oruanui eruption: Insights into the generation and dynamics of the world's youngest supereruption, PhD thesis, Victoria University of Wellington, (last access: 21 February 2022), 2013. a, b

Allan, A. S. R., Morgan, D. J., Wilson, C. J. N., and Millet, M.-A.: From mush to eruption in centuries: assembly of the super-sized Oruanui magma body, Contrib. Mineral. Petr., 166, 143–164,, 2013. a

Aman, Z., Wen-shan, Y., and Xiong-liang, Y.: Numerical simulation of underwater contact explosion, Appl. Ocean Res., 34, 10–20,, 2012. a

Barker, S. J., Wilson, C. J. N., Smith, E. G. C., Charlier, B. L. A., Wooden, J. L., Hiess, J., and Ireland, T. R.: Post-supereruption magmatic reconstruction of Taupo volcano (New Zealand), as reflected in zircon ages and trace elements, J. Petrol., 55, 1511–1533,, 2014. a, b

Barker, S. J., Wilson, C. J. N., Allan, A. S. R., and Schipper, C. I.: Fine-scale temporal recovery, reconstruction and evolution of a post-supereruption magmatic system, Contrib. Mineral. Petr., 170, 1–40,, 2015. a, b, c, d, e

Barker, S. J., Van Eaton, A. R., Mastin, L. G., Wilson, C. J. N., Thompson, M. A., Wilson, T. M., Davis, C., and Renwick, J. A.: Modeling ash dispersal from future eruptions of Taupo supervolcano, Geochem. Geophy. Geosy., 20, 3375–3401,, 2019. a

Barker, S. J., Wilson, C. J. N., Illsley-Kemp, F., Leonard, G. S., Mestel, E. R. H., Mauriohooho, K., and Charlier, B. L. A.: Taupō: an overview of New Zealand's youngest supervolcano, New Zeal. J. Geol. Geop., 64, 1–27,, 2020. a

Bebbington, M. S.: Temporal-volume probabilistic hazard model for a supervolcano: Taupo, New Zealand, Earth Planet. Sc. Lett., 536, 116141,, 2020. a

Belousov, A., Voight, B., Belousova, M., and Muravyev, Y.: Tsunamis generated by subaquatic volcanic explosions: unique data from 1996 eruption in Karymskoye Lake, Kamchatka, Russia, Pure Appl. Geophys., 157, 1135–1143,, 2000. a

Berny, A., Deike, L., Séon, T., and Popinet, S.: Role of all jet drops in mass transfer from bursting bubbles, Physical Review Fluids, 5, 033605,, 2020. a

Bibby, H. M., Caldwell, T. G., Davey, F. J., and Webb, T. H.: Geophysical evidence on the structure of the Taupo Volcanic Zone and its hydrothermal circulation, J. Volcanol. Geoth. Res., 68, 29–58,, 1995. a

Cole, R. H.: Underwater Explosions, Princeton University Press, ISBN-13: 978-0486613840, 1948. a

Daramizadeh, A. and Ansari, M.: Numerical simulation of underwater explosion near air–water free surface using a five-equation reduced model, Ocean Eng., 110, 25–35,, 2015. a

Davy, B. W. and Caldwell, T. G.: Gravity, magnetic and seismic surveys of the caldera complex, Lake Taupo, North Island, New Zealand, J. Volcanol. Geoth. Res., 81, 69–89,, 1998. a, b

de Ronde, C. E. J., Stoffers, P., Garbe-Schönberg, D., Christenson, B. W., Jones, B., Manconi, R., Browne, P. R. L., Hissmann, K., Botz, R., Davy, B. W., Schmitt, M., and Battershill, C. N.: Discovery of active hydrothermal venting in Lake Taupo, New Zealand, J. Volcanol. Geoth. Res., 115, 257–275,, 2002. a

Dietz, R. S. and Sheehy, M. J.: Transpacific detection of Myojin volcanic explosions by underwater sound, Geol. Soc. Am. Bull., 65, 941–956,[941:tdomve];2, 1954. a

Dondin, F., Lebrun, J. F., Kelfoun, K., Fournier, N., and Randrianasolo, A.: Sector collapse at Kick 'em Jenny submarine volcano (Lesser Antilles): Numerical simulation and landslide behaviour, B. Volcanol., 74, 595–607,, 2012. a

Duffy, D. G.: On the generation of oceanic surface waves by underwater volcanic explosions, J. Volcanol. Geoth. Res., 50, 323–344,, 1992. a

Egorov, Y.: Tsunami wave generation by the eruption of underwater volcano, Nat. Hazards Earth Syst. Sci., 7, 65–69,, 2007. a, b, c

Falvard, S., Paris, R., Belousova, M., Belousov, A., Giachetti, T., and Cuven, S.: Scenario of the 1996 volcanic tsunamis in Karymskoye Lake, Kamchatka, inferred from X-ray tomography of heavy minerals in tsunami deposits, Mar. Geol., 396, 160–170,, 2018. a

Gelman, S. E., Deering, C. D., Gutierrez, F. J., and Bachmann, O.: Evolution of the Taupo Volcanic Center, New Zealand: petrological and thermal constraints from the Omega dacite, Contrib. Mineral. Petr., 166, 1355–1374,, 2013. a

Glasbergen, T.: Characterization of incoming tsunamis for the design of coastal structures: A numerical study using the SWASH model, master's thesis, published by Delft University of Technology, (last accessed: 17 February 2022), 2018. a

Grilli, S. T., Harris, J. C., Shi, F., Kirby, J. T., Bakhsh, T. S. T., Estibals, E., and Tehranirad, B.: Numerical modeling of coastal tsunami impact dissipation and impact, Coastal Engineering Proceedings, 1,, 2012. a

Grilli, S. T., O’Reilly, C., Harris, J. C., Bakhsh, T. T., Tehranirad, B., Banihashemi, S., Kirby, J. T., Baxter, C. D. P., Eggeling, T., Ma, G., and Shi, F.: Modeling of SMF tsunami hazard along the upper US East Coast: detailed impact around Ocean City, MD, Nat. Hazards, 76, 705–746,, 2015. a

Grilli, S. T., Tappin, D. R., Carey, S., Watt, S. F., Ward, S. N., Grilli, A. R., Engwell, S. L., Zhang, C., Kirby, J. T., Schambach, L., and Muin, M.: Modelling of the tsunami from the December 22, 2018 lateral collapse of Anak Krakatau volcano in the Sunda Straits, Indonesia, Sci. Rep., 9, 1–13,, 2019. a, b, c

Grilli, S. T., Zhang, C., Kirby, J. T., Grilli, A. R., Tappin, D. R., Watt, S. F. L., Hunt, J. E., Novellino, A., Engwell, S., Nurshal, M. E. M., Abdurrachman, M., Cassidy, M., Madden-Nadeau, A. L., and Day, S.: Modeling of the Dec. 22nd 2018 Anak Krakatau volcano lateral collapse and tsunami based on recent field surveys: Comparison with observed tsunami impact, Mar. Geol., 440, 106566,, 2021. a

Hallquist, J. O.: LS-DYNA3D theoretical manual, Livermore software technology corporation Livermore, Calif., 1994. a

Hammack, J. L.: A note on tsunamis: their generation and propagation in an ocean of uniform depth, J. Fluid Mech., 60, 769–799,, 1973. a

Hammack, J. L. and Segur, H.: The Korteweg-de Vries equation and water waves. Part 2. Comparison with experiments, J. Fluid Mech., 65, 289–314,, 1974. a

Hammack, J. L. and Segur, H.: The Korteweg-de Vries equation and water waves. Part 3. Oscillatory waves, J. Fluid Mech., 84, 337–358,, 1978a. a

Hammack, J. L. and Segur, H.: Modelling criteria for long water waves, J. Fluid Mech., 84, 359–373,, 1978b. a

Hayward, M.: Simulation of Lake Taupō explosive tsunami, TIB AV-Portal [video],, 2021. a

Hayward, M. W.: Multilayer modelling of waves from disturbances and explosions: Cases at laboratory-scale and field-scale at Mono Lake, California, Zenodo [code],, 2022. a

Houghton, B. F., Carey, R. J., Cashman, K. V., Wilson, C. J. N., Hobden, B. J., and Hammer, J. E.: Diverse patterns of ascent, degassing, and eruption of rhyolite magma during the 1.8 ka Taupo eruption, New Zealand: evidence from clast vesicularity, J. Volcanol. Geoth. Res., 195, 31–47,, 2010. a

Illsley-Kemp, F., Barker, S. J., Wilson, C. J. N., Chamberlain, C. J., Hreinsdóttir, S., Ellis, S., Hamling, I. J., Savage, M. K., Mestel, E. R. H., and Wadsworth, F. B.: Volcanic unrest at Taupō volcano in 2019: Causes, mechanisms and implications, Geochem. Geophy. Geosy., 22, e2021GC009803,, 2021. a

Johnson, R. W.: Large-scale volcanic cone collapse: the 1888 slope failure of Ritter volcano, and other examples from Papua New Guinea, B. Volcanol., 49, 669–679,, 1987. a

Kedrinskiy, V. K.: Hydrodynamics of Explosion: Experiments and Models, Springer, ISBN: 978-3-642-06130-1, 2006. a

Klaseboer, E., Hung, K. C., Wang, C., Wang, C. W., Khoo, B. C., Boyce, P., Debono, S., and Charlier, H.: Experimental and numerical investigation of the dynamics of an underwater explosion bubble near a resilient/rigid structure, J. Fluid Mech., 537, 387–413,, 2005. a

Lane, E. M., Borrero, J., Whittaker, C. N., Bind, J., Chagué-Goff, C., Goff, J., Goring, D., Hoyle, J., Mueller, C., Power, W. L., Reid, C. M., Williams, J. H., and Williams, S. P.: Effects of inundation by the 14th november, 2016 Kaikōura tsunami on banks Peninsula, Canterbury, New Zealand, Pure Appl. Geophys., 174, 1855–1874,, 2017. a

Le Méhauté, B.: Theory of explosion-generated water waves, Adv. Hydrosci., 7, 1–79, 1971. a, b

Le Méhauté, B. and Wang, S.: Water Waves Generated By Underwater Explosion, World Scientific, ISBN: 978-981-02-2083-9, 1996. a, b, c, d, e, f, g

Leonard, G. S., Begg, J. G., and Wilson, C. J. N.: Geology of the Rotorua Area: Institute of Geological & Nuclear Sciences 1 : 250 000 Geological Map 5, 1 Sheet + 102 pp., GNS Science, Lower Hutt, New Zealand, ISBN: 978-0-478-19778-5, 2010. a

Li, T., Wang, S., Li, S., and Zhang, A.-M.: Numerical investigation of an underwater explosion bubble based on FVM and VOF, Appl. Ocean Res., 74, 49–58,, 2018. a

Liu, Y. L., Zhang, A. M., Tian, Z. L., and Wang, S. P.: Numerical investigation on global responses of surface ship subjected to underwater explosion in waves, Ocean Eng., 161, 277–290,, 2018. a

López-Herrera, J., Popinet, S., and Castrejón-Pita, A.: An adaptive solver for viscoelastic incompressible two-phase problems applied to the study of the splashing of weakly viscoelastic droplets, J. Non-Newton. Fluid, 264, 144–158,, 2019. a

López-Venegas, A. M., Horrillo, J., Pampell-Manis, A., Huérfano, V., and Mercado, A.: Advanced tsunami numerical simulations and energy considerations by use of 3D–2D coupled models: The October 11, 1918, Mona passage tsunami, Pure Appl. Geophys., 172, 1679–1698,, 2015. a

Ma, G., Shi, F., and Kirby, J. T.: Shock-capturing non-hydrostatic model for fully dispersive surface wave processes, Ocean Model., 43, 22–35,, 2012. a

Ma, G., Kirby, J. T., and Shi, F.: Numerical simulation of tsunami waves generated by deformable submarine landslides, Ocean Model., 69, 146–165,, 2013. a

Manville, V.: Sedimentary and geomorphic responses to ignimbrite emplacement: readjustment of the Waikato River after the AD 181 Taupo eruption, New Zealand, J. Geol., 110, 519–541,, 2002. a

Manville, V., White, J. D. L., Houghton, B. F., and Wilson, C. J. N.: Paleohydrology and sedimentology of a post–1.8 ka breakout flood from intracaldera Lake Taupo, North Island, New Zealand, Geol. Soc. Am. Bull., 111, 1435–1447,<1435:pasoap>;2, 1999. a

Mastin, L. G. and Witter, J. B.: The hazards of eruptions through lakes and seawater, J. Volcanol. Geoth. Res., 97, 195–214,, 2000. a

Mirchina, N. and Pelinovsky, E.: Estimation of underwater eruption energy based on tsunami wave data, Nat. Hazards, 1, 277–283,, 1988. a

Nomanbhoy, N. and Satake, K.: Generation mechanism of tsunamis from the 1883 Krakatau eruption, Geophys. Res. Lett., 22, 509–512,, 1995. a

Pakoksung, K., Suppasri, A., and Imamura, F.: Probabilistic Tsunami Hazard Analysis of Inundated Buildings Following a Subaqueous Volcanic Explosion Based on the 1716 Tsunami Scenario in Taal Lake, Philippines, Geosciences, 11, 92,, 2021. a

Paris, R.: Source mechanisms of volcanic tsunamis, Philos. T. Roy. Soc. A, 373, 20140380,, 2015. a

Paris, R. and Ulvrová, M.: Tsunamis generated by subaqueous volcanic explosions in Taal Caldera Lake, Philippines, B. Volcanol., 81, 14,, 2019. a, b, c, d

Paris, R., Switzer, A. D., Belousova, M., Belousov, A., Ontowirjo, B., Whelley, P. L., and Ulvrova, M.: Volcanic tsunami: A review of source mechanisms, past events and hazards in Southeast Asia (Indonesia, Philippines, Papua New Guinea), Nat. Hazards, 70, 447–470,, 2014. a

Paris, R., Ulvrová, M., Selva, J., Brizuela, B., Costa, A., Grezio, A., Lorito, S., and Tonini, R.: Probabilistic hazard analysis for tsunamis generated by subaqueous volcanic explosions in the Campi Flegrei caldera, Italy, J. Volcanol. Geoth. Res., 379, 106–116,, 2019. a, b, c

Petrov, N. V. and Schmidt, A. A.: Multiphase phenomena in underwater explosion, Exp. Therm. Fluid Sci., 60, 367–373,, 2015. a

Pinkston, J., Skinner, F. W., and Strange, J. N.: Results of Surface Wave Experiments: Mono Lake Explosion Test Series, 1965, Waterways Experiment Station, 1970. a

Popinet, S.: An accurate adaptive solver for surface-tension-driven interfacial flows, J. Comput. Phys., 228, 5838–5866,, 2009. a

Popinet, S.: Quadtree-adaptive tsunami modelling, Ocean Dynam., 61, 1261–1285,, 2011. a

Popinet, S.: Basilisk, Darcs [code], (last access: 9 April 2021), 2013. a

Popinet, S.: A quadtree-adaptive multigrid solver for the Serre–Green–Naghdi equations, J. Comput. Phys., 302, 336–358,, 2015. a

Popinet, S.: Numerical models of surface tension, Annu. Rev. Fluid Mech., 50, 49–75,, 2018. a

Popinet, S.: A vertically-Lagrangian, non-hydrostatic, multilayer model for multiscale free-surface flows, J. Comput. Phys., 148, 109609,, 2020. a, b, c, d, e

Potter, S. H., Scott, B. J., Jolly, G. E., Johnston, D. M., and Neall, V. E.: A catalogue of caldera unrest at Taupo Volcanic Centre, New Zealand, using the volcanic unrest index (VUI), B. Volcanol., 77, 1–28,, 2015. a

Prins, J.: Characteristics of waves generated by a local disturbance, Eos, Transactions American Geophysical Union, 39, 865–874,, 1958. a, b, c, d, e

Pyle, D. M.: Chapter 13 – Sizes of Volcanic Eruptions, in: The Encyclopedia of Volcanoes (Second Edition), edited by: Sigurdsson, H., Academic Press, Amsterdam, 2nd edn., 257–264,, 2015. a

Raumann, C. G., Stine, S., Evans, A., and Wilson, J.: Digital bathymetric model of mono lake, california, Miscellaneous field studies map MF-2393, 2002. a

Rowe, D. K., Shankar, U., James, M., and Waugh, B.: Use of GIS to predict effects of water level on the spawning area for smelt, Retropinna retropinna, in Lake Taupo, New Zealand, Fisheries Manag. Ecol., 9, 205–216,, 2002. a

Ruffini, G., Heller, V., and Briganti, R.: Numerical modelling of landslide-tsunami propagation in a wide range of idealised water body geometries, Coast. Eng., 153, 103518,, 2019. a

Sato, H. and Taniguchi, H.: Relationship between crater size and ejecta volume of recent magmatic and phreato-magmatic eruptions: Implications for energy partitioning, Geophys. Res. Lett., 24, 205–208,, 1997. a, b

Schambach, L., Grilli, S. T., Kirby, J. T., and Shi, F.: Landslide tsunami hazard along the upper US East Coast: effects of slide deformation, bottom friction, and frequency dispersion, Pure Appl. Geophys., 176, 3059–3098,, 2019. a

Schambach, L., Grilli, S. T., and Tappin, D. R.: New high-resolution modeling of the 2018 Palu tsunami, based on supershear earthquake mechanisms and mapped coastal landslides, supports a dual source, Front. Earth Sci., 8, 598839,, 2021. a, b

Shen, Y., Whittaker, C. N., Lane, E. M., White, J. D. L., Power, W., and Melville, B. W.: Waves generated by discrete and sustained gas eruptions with implications for submarine volcanic tsunamis, Geophys. Res. Lett., 48, e2021GL094539,, 2021a. a

Shen, Y., Whittaker, C. N., Lane, E. M., White, J. D. L., Power, W., and Nomikou, P.: Laboratory experiments on tsunamigenic discrete subaqueous volcanic eruptions. Part 1: Free surface disturbances, J. Geophys. Res.-Oceans, 126, e2020JC016588,, 2021b. a

Shen, Y., Whittaker, C. N., Lane, E. M., White, J. D. L., Power, W., and Nomikou, P.: Laboratory experiments on tsunamigenic discrete subaqueous volcanic eruptions. Part 2: Properties of generated waves, J. Geophys. Res.-Oceans, 126, e2020JC016587,, 2021c. a

Smith, M. S. and Shepherd, J. B.: Preliminary investigations of the tsunami hazard of Kick'em Jenny submarine volcano, Nat. Hazards, 7, 257–277,, 1993. a

Smith, M. S. and Shepherd, J. B.: Tsunami waves generated by volcanic landslides: an assessment of the hazard associated with Kick’em Jenny, Geol. Soc. Spec. Publ., 110, 115–123,, 1996. a

Sutton, A. N., Blake, S., Wilson, C. J. N., and Charlier, B. L.: Late Quaternary evolution of a hyperactive rhyolite magmatic system: Taupo volcanic centre, New Zealand, J. Geol. Soc., 157, 537–552,, 2000. a

Torsvik, T., Paris, R., Didenkulova, I., Pelinovsky, E., Belousov, A., and Belousova, M.: Numerical simulation of a tsunami event during the 1996 volcanic eruption in Karymskoye lake, Kamchatka, Russia, Nat. Hazards Earth Syst. Sci., 10, 2359–2369,, 2010. a, b, c, d

Ulvrová, M., Paris, R., Kelfoun, K., and Nomikou, P.: Numerical simulations of tsunamis generated by underwater volcanic explosions at Karymskoye lake (Kamchatka, Russia) and Kolumbo volcano (Aegean Sea, Greece), Nat. Hazards Earth Syst. Sci., 14, 401–412,, 2014. a, b, c, d

Valentine, G. A. and White, J. D.: Revised conceptual model for maar-diatremes: Subsurface processes, energetics, and eruptive products, Geology, 40, 1111–1114,, 2012. a

Vandergoes, M. J., Hogg, A. G., Lowe, D. J., Newnham, R. M., Denton, G. H., Southon, J., Barrell, D. J. A., Wilson, C. J. N., McGlone, M. S., Allan, A. S. R., Almond, P. C., Petchey, F., Dabell, K., Dieffenbacher-Krall, A. C., and Blaauw, M.: A revised age for the Kawakawa/Oruanui tephra, a key marker for the Last Glacial Maximum in New Zealand, Quaternary Sci. Rev., 74, 195–201,, 2013. a

van Hooft, J. A., Popinet, S., van Heerwaarden, C. C., van der Linden, S. J. A., de Roode, S. R., and van de Wiel, B. J. H.: Towards adaptive grids for atmospheric boundary-layer simulations, Bound.-Lay. Meteorol., 167, 421–443,, 2018. a

Wallace, N. and Baird, C.: Explosion-Generated Waves 1965 Mono Lake Test Series, Tech. rep., Oceanographic Services Inc., Santa Barbara, CA, US Defence Technical Information Center, Accession number: AD0682962, 1968. a, b, c

Walter, D.: Explosion-generated Wave Tests, Mono lake, California: Ground and Aerial Photography, Tech. rep., URS Systems Corp, Burlingame, CA, US Defence Technical Information Center, Accession number: AD0627354, 1966. a, b

Wang, S.-P., Zhang, A.-M., Liu, Y.-L., Zhang, S., and Cui, P.: Bubble dynamics and its applications, J. Hydrodyn., 30, 975–991,, 2018. a

Ward, S. N. and Asphaug, E.: Asteroid impact tsunami: a probabilistic hazard assessment, Icarus, 145, 64–78,, 2000. a

Whalin, R. W., Pace, C. E., and Lane, W. F.: Mono Lake Explosion Test Series, 1965: Analysis of Surface Wave and Wave Runup Data, Waterways Experiment Station, US Army Technical Report N-70-12, 1970. a, b

Williams, R., Rowley, P., and Garthwaite, M. C.: Reconstructing the Anak Krakatau flank collapse that caused the December 2018 Indonesian tsunami, Geology, 47, 973–976,, 2019.  a

Wilson, C. J. N.: Stratigraphy, chronology, styles and dynamics of late Quaternary eruptions from Taupo volcano, New Zealand, Philos. T. Roy. Soc. A, 343, 205–306,, 1993. a, b, c

Wilson, C. J. N.: The 26.5 ka Oruanui eruption, New Zealand: an introduction and overview, J. Volcanol. Geoth. Res., 112, 133–174,, 2001. a

Wilson, C. J. N. and Charlier, B. L. A.: Rapid rates of magma generation at contemporaneous magma systems, Taupo Volcano, New Zealand: insights from U–Th model-age spectra in zircons, J. Petrol., 50, 875–907,, 2009. a

Wilson, C. J. N. and Walker, G. P. L.: The Taupo eruption, New Zealand I. General aspects, Philos. T. Roy. Soc. A, 314, 199–228,, 1985. a

Wilson, C. J. N., Blake, S., Charlier, B. L. A., and Sutton, A. N.: The 26.5 ka Oruanui eruption, Taupo volcano, New Zealand: development, characteristics and evacuation of a large rhyolitic magma body, J. Petrol., 47, 35–69,, 2006. a

Wilson, C. J. N., Gravley, D. M., Leonard, G. S., and Rowland, J. V.: Volcanism in the central Taupo Volcanic Zone, New Zealand: tempo, styles and controls, Studies in volcanology: the legacy of George Walker, Special Publications of IAVCEI, 2, 225–247,, 2009. a

Xu, L.-Y., Wang, S.-P., Liu, Y.-L., and Zhang, A.-M.: Numerical simulation on the whole process of an underwater explosion between a deformable seabed and a free surface, Ocean Eng., 219, 108311,, 2020. a

Ye, L., Kanamori, H., Rivera, L., Lay, T., Zhou, Y., Sianipar, D., and Satake, K.: The 22 December 2018 tsunami from flank collapse of Anak Krakatau volcano during eruption, Sci. Adv., 6, eaaz1377,, 2020. a

Zijlema, M., Stelling, G., and Smit, P.: SWASH: An operational public domain code for simulating wave fields and rapidly varied flows in coastal waters, Coast. Eng., 58, 992–1012,, 2011. a

Short summary
Volcanic eruptions can produce tsunamis through multiple mechanisms. We present validation cases for a numerical method used in simulating waves caused by submarine explosions: a laboratory flume experiment and waves generated by explosions at field scale. We then demonstrate the use of the scheme for simulating analogous volcanic eruptions, illustrating the resulting wavefield. We show that this scheme models such dispersive sources more proficiently than standard tsunami models.
Final-revised paper