Multilayer modelling of waves generated by explosive submarine volcanism

. Theoretical source models of underwater explosions are often applied in studying tsunami hazards associated with submarine volcanism; however, their use in numerical codes based on the shallow water equations can neglect the signiﬁcant dispersion of the generated waveﬁeld. A non-hydrostatic multilayer method is validated against a laboratory-scale experiment of wave generation from instantaneous disturbances and at ﬁeld-scale submarine explosions at Mono Lake, California, utilising the relevant theoretical models. The numerical method accurately reproduces the range of observed wave characteristics for 5 positive disturbances and suggests a previously unreported relationship of extended initial troughs for negative disturbances at low dispersivity and high nonlinearity 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¯o, New Zealand, for a magnitude range representing ejecta volumes between 0.04 – 0.4 km 3 . Waves reach all shores within 15 minutes with maximum incident crest amplitudes around 4 m at shores near the 10 source. This work shows that the multilayer scheme used is computationally efﬁcient and able to capture a wide range of wave characteristics, including dispersive effects, which is necessary when investigating submarine explosions. This research therefore provides the foundation for future studies involving a rigorous probabilistic hazard assessment to quantify the risks and relative signiﬁcance of this tsunami source mechanism. with data revealing areas of varying exposure to hazardous waves, of above 4 m near source to 0.3 m in the most sheltered areas, and waves reaching throughout the lake within 15 minutes. A probabilistic investigation is 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. This will help resolve the signiﬁcance of this hazard source compared 410 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.

. However, tsunamis are often neglected from hazard maps of volcanos -recent events such as 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. 25 While experimental data and detailed field studies of these phenomena are rare, some reliable observations exist, for example, at Karymskoye Lake, 1996, where 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 30 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 1972in , 1974in and 2007in (Johnson, 1987Dondin 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 Shepherd, 1993, 1996). Many other candidate eruptions are historically documented with 35 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 Sheehy, 1954;Nomanbhoy and Satake, 1995).
Explosive volcanic eruptions are characterised by a directional gas-driven escape from the source, exsolution of water vapour and, in submarine settings, potentially violent vaporisation of sea or lake water on interaction with hot magma. This 40 phreatomagmatic eruption can lead 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 is needed.
Underwater explosions are well documented (Cole, 1948;Mirchina and Pelinovsky, 1988;Le Méhauté, 1971;Kedrinskiy, 2006;Egorov, 2007) primarily owing to military reports and research in blast mitigation and structural response in, for example, 45 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 nucleartesting age and led to the development of theoretical models describing explosion-surface interaction and dynamics of the resultant wave field (Le Méhauté and Wang, 1996). Physical experimentation since the end of nuclear testing has been rare 50 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 submarine volcanic explosions 55 (Torsvik et al., 2010;Ulvrová et al., 2014; and asteroids impacting in ocean 2 https://doi.org/10.5194/nhess-2021-109 Preprint. Discussion started: 7 July 2021 c Author(s) 2021. CC BY 4.0 License.
basins (Ward and Asphaug, 2000). 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 60 waves may be .
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 65 theoretical model's initial conditions in the numerical method. These tests are to establish fitness of the underlying models, which are then applied to hypothetical explosive submarine eruptions at Lake Taupō, New Zealand.

Underwater eruption model
An explosive subaqueous volcanic eruption is a dynamic and complex event involving abrupt fragmentation, volume change 70 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 (Egorov, 2007;Paris, 2015). 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 .

75
The models developed for submarine 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 80 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 and significantly reduce wave-making efficiency.
(Le Méhauté, 1971;Le Méhauté and Wang, 1996) Bubble dynamics is a very active area of research in computational fluid dynamics (CFD), though, in the explosive realm, 85 the focus is usually on pressure waves and solid interactions . These studies are usually short in temporal range and are very computationally expensive as modelling the full problem requires accounting for compressibility and multiphase flow; thus, this has spawned specialist codes for their solution (Hallquist, 1994;Li et al., 2018). Only in recent years have studies appeared that directly simulate expanding explosive bubbles interacting with deformable beds and a free-surface 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). (Petrov and Schmidt, 2015;Daramizadeh and Ansari, 2015;Xu et al., 2020). However, there is minimal focus on subsequent 90 surface waves, let alone relations tested for their generation mechanisms or far-field propagation.
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 Wang, 1996). η c corresponds to the maximum depth of the disturbance below equilibrium and R to its radius 95 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 4 https://doi.org/10.5194/nhess-2021-109 Preprint. Discussion started: 7 July 2021 c Author(s) 2021. CC BY 4.0 License. described analytically by one of a number of candidates for the general profile as Eq. (1-2). 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 Figure 1a-c.

Depth classification
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 105 released E. In terms of a depth parameter D = ch 3 √ E , where an imperial unit conversion constant c = 406.2, three categories are specified when considering wave generation: For deep and intermediate cases, the cavity parameters are defined as: 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: 120 The data calibrating these models include charges ranging from small (< 500 lb or < 9.5 × 10 8 J) to a handful larger (< 9500 lb or < 1.8 × 10 10 J) and further include a 23-KT nuclear test (Le Méhauté and Wang, 1996).

Volcanic context
For a volcanic case, illustrated in Fig. 1d, such explosions would occur on or near 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 . This is valid where the released explosion energy is estimated from volcanic crater diameter C D using the following empirical relationship by Sato and Taniguchi (1997): This method has recently been used for probabilistic hazard analysis of volcanogenic tsunamis at the Campi Flegrei caldera, Italy  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 White, 2012), so a further relation from Sato and Taniguchi (1997) relates the ejecta volume V with the released explosion energy:

Numerical method
To compute simulations of the models described earlier, we use the open source CFD framework, Basilisk (Popinet, 2013).
The software is widely used in studies involving multiphase problems from jet dynamics to viscoelastic and surface tension investigations and includes several free-surface solvers with application to tsunamis, wavefield transformation and other hy-

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 system: where, in the x-z reference frame, k is the layer index, h k layer thickness, g gravitational acceleration, u k , w k the horizontal and vertical velocity components, φ k the non-hydrostatic pressure, η the free-surface height (sum of layer thicknesses and bathymetry height z b ), and 160 the height of layer interfaces.
Between them, the set expresses the evolution of layer thickness (Eq. 11), conservation of momentum (Eq. 12, 13), and conservation of volume/mass (Eq. 14). The framework allows the model to be built modularly, starting from the hydrostatic case where φ = 0 and vertical momentum conservation (Eq. 13) is removed; these effectively become the generalised multilayer Saint-Venant (SV) or stacked shallow water equations. More components are added, for example, vertical remapping, adaptiv-165 ity, 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 where g |H| ∞ is the characteristic horizontal velocity scale of the wave and b is a specified breaking parameter smaller than one. sgn and min are sign and minimum functions respectively. Lastly, terrain is handled by looking up a pre-processed 170 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 propagation and shore run-up. (Popinet, 2020)

Other Schemes
Accompanying models used for validation and comparative purposes include a Volume-of-Fluid method (N-S/VOF) in the same framework which solves the two-phase Navier-Stokes equations for interfacial flows, including variable density and surface tension (Popinet, 2009(Popinet, , 2018, a solver for the shallow water or Saint Venant equations (SV), and finally another 180 for the Serre-Green-Naghdi equations (SGN), a Boussinesq higher order approximation for non-linear and weakly dispersive flows (Popinet, 2015). Their inclusion is to support and inform 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 185 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 capability to include non-hydrostatic terms; 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.
3 Laboratory-scale validation 190 To determine suitability of the numerical method for use in modelling a submarine 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-    (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 220 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.
The generation process is shown in Figure 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 where the amplitude is proportional to the initial disturbance  solver. Also, note graphical artefacts on free-surface height at the leading gradient in the N-S 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.

235
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 11 https://doi.org/10.5194/nhess-2021-109 Preprint. Discussion started: 7 July 2021 c Author(s) 2021. CC BY 4.0 License. with the N-S/VOF solver considering result similarity. It also demonstrates superior speed than the dispersive SGN method, primarily due to the computation time required to solve the higher order approximation.

Wavefield classification
All simulated cases are plotted in Figure 5 by the size of initial disturbance relative to water depth using a nonlinearity parameter |Q| h against a dispersivity parameter k*h 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 k*h > 10 |Q| h , tending 245 through increasingly solitary wave properties (purple and green) once k*h ≤ 10 |Q| h until a succession of diminishing amplitude solitary waves result where k*h < |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 k*h. 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, those of 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 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 255 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 k*h. 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).
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-260 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.

Field-scale validation
The next stage is to assess 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 265 Walter (1966); Wallace and Baird (1968); Whalin et al. (1970); 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 ten approximately 9,250 lb (∼4,196 kg) spherical TNT charges were detonated off the south shore of Mono Lake, California. The test area is illustrated in Figure 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. are given in Table 2. Considering doubts regarding the charges raised in the report of Walter (1966), additional lower yield simulations were added.     Figure 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 285 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 has noticeably greater variation at shore.
The experimental gauge time series at locations closest to GZ are plotted alongside the lower yield model traces in Figure 8 and are useful to compare the phase arrival times and the initial development of the wave train. For both simulations, the first 290 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, where the differences are 9 and 7 seconds 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.

Model implications
While the models used fit the experimental results and trends well overall, it is significant that the lower yield data is 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. A further similar deep test at Mono Lake in the following year delivered much greater amplitude waves in line with expectations, 300 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.
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.

305
Instrumental issues meant gauges near the shore were often unable to be acceptably calibrated for amplitude. Many experienced +ve 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 Baird, 1968) Despite the experimental challenges, the multilayer scheme is shown to be excellent at resolving the initialised disturbance into a wave field very consistent with that generated by an intermediate depth submarine explosion. With the capability to 310 accurately propagate such a source in the relatively deep near-field through to the shore, it demonstrates suitability to model such events at these scales. However, this test did not contain any shallow depth underwater explosion tests, for which there are no case examples at this scale.

Taupō Scenario
Lake Taupō is New Zealand's largest freshwater lake of area approximately 616 km 2 . It lies in the southern section of the Taupō   315 Volcanic Zone and conceals most vent sites and features of Taupō volcano; the lake itself is the result of caldera collapse after the c. 25.5 ka Oruanui supereruption, with modifications in following events. The lake has experienced numerous volcanic episodes post-Oruanui, with some 21 events occurring from 7.05-1.8 ka, including the ∼ 232 CE Taupō eruption which is globally one of the most powerful eruptions in the Holocene. (Barker et al., 2020) Figure 9 illustrates the lake bathymetry and surrounding terrain accompanied by geological features and settlement locations.

320
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 one million tourists each year, and is of great cultural significance accompanied by the surrounding land.

325
In building a model representing an example tsunamigenic explosive eruption in Lake Taupō, it was required to build a terrain dataset by combining a bathymetric model of the lake (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 (Barker et al., 2020) and, for this example, two simulations are run corresponding to the range of most eruptions occurring during this period, excepting three larger episodes. The estimated ejecta volume V is between 0.04 -0.4 km 3 (Wilson, 1993) and, using Eqs. (9-10), the equivalent resulting models were run for 1000 seconds of simulated time. A further simulation of the scenario was performed with the Saint Venant scheme with identical terrain and model geometry for comparison. The initial wave at all locations is a crest, reflecting the positive amplitude lip of the initial condition. Earliest arrivals at 345 shore occur about the closest point east of the source at 3 min, after which arrival time generally scales with radial distance excepting 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 4 m. The lowest are found in the further area of the south-west besides 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 350 lake morphology. 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 ∼65 s early in 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. Gauge 4 and 5, positioned in shallow zones near the eruption, initially record bores of amplitude 1.6 m and 3.5 m 355 respectively in the high run. While phase velocities are very similar between runs, the generated group velocity is slower for the low run. This is best seen in Gauges 2, 3 and 6 which are positioned with relatively unobstructed direct paths from the source.  numerical wave breaking therefore adds robustness to the near-shore solution for this type of source. Notably, the explosion model prescribes an initial condition that intersects the bathymetry in the present high ejecta volume case, resulting in additional 365 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; however, no alteration to the explosion model is made in present work.

Hazard implications 370
These preliminary results suggest many potential implications for wave hazard from explosive subaqueous eruptions. If an eruption of sufficient magnitude at Lake Taupō produces an initial explosion there is clearly a threat posed to nearby shores.
As with most lacustrine tsunami hazards, there is minimal time from source to shore impact; no possible warning system would ever be able to respond to an eruption with sufficient speed. Instead, resilience should come from preparation for events including disaster management and exclusion zone planning by local authorities, with ongoing monitoring of the volcanic 375 system such as seismicity, ground deformation, changes to the geothermal system, and geophysical imaging.
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 at least in the upper range considered here (0.2 km 3 ) is estimated to have 5% probability within 100 years (Bebbington, 2020). Taupō volcano can produce far greater magnitude events than considered here, including the aforementioned CE 232 Taupō eruption (24 km 3 ) and the c. 25.5 380 ka Oruanui supereruption (> 1100 km 3 ). Events of such magnitude 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 submarine volcanism will be to determine the relative weight each source component possesses with events of varying location and magnitude.

385
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, potentially able to 390 investigate inundation in detail at various points onshore with, for instance, building data or similar additional model layers.

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 395 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 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 400 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 use of the underlying empirical relations.
These were then used to simulate a volcanic explosion under Lake Taupō based on estimated eruptive energy by ejecta 405 volume. Implications for tsunami hazard under lakes are suggested for small to medium magnitude phreatomagmatic eruptions, and surrounding topography. We acknowledge the use of New Zealand eScience Infrastructure (NeSI) high performance computing facilities as part of this research.