Physical modeling of spikes during the volcanic tsunami generation

Volcanic tsunamis are water waves generated by several eruption processes, ground uplift, and slope instabilities at volcanoes. Certain amounts of mechanical or thermal energy from volcanic activities are transferred into water bodies such as oceans or lakes (Sigurdsson , 2015). The majority of volcanic tsunamis is generated by underwater eruptions, pyroclastic flows, and volcanic earthquakes, which sum up to around 60% of all volcanic tsunamis (Latter, 1981). The recorded historical volcanic tsunami event with most fatalities is the 1883 Krakatau eruption, which caused more than 36 000 victims (Simkin, 1983). The volcanic cone collapse during the 2018 Anak Krakatau volcanic eruption generated a tsunami in the Sunda Strait and caused localized runup heights exceeding 80 m on two neighboring islands (Borrero , 2020; Syamsidik , 2020). The tsunami waves generated by underwater explosion and atmospheric pressure shockwave in the 2022 Tonga eruption reached Tongatapu coastline with runup heights up to 15 m (Borrero , 2023). Submarine volcanic eruption generated tsunamis can be described as cylindrical water waves with a leading N-wave and a train of trailing dispersive waves (Le Méhauté and Wang, 1996; Paris, 2015). The source mechanisms and wave characteristics of volcanic tsunamis need to be investigated besides the existing source models of seismically generated tsunamis (Yeh , 2008; Liu , 2009; Behrens and Dias, 2015; and Wei , 2015).

Volcanic flank collapses and mass flow-induced tsunamis have been studied in the laboratory with solid and granular mass flows in flumes. Rigid bodies of various shapes, lengths, heights, and weights sliding down a slope were adopted to model the landslide tsunamis in laboratory (Law and Brebner, 1968; Watts, 1997; Watts, 2000; and Enet and Grilli, 2007). Extensive granular landslide generated tsunami experiments conducted in a flume yielded energy conversion rates from the slide to the first wave ranging from 2% to 30% (Fritz, 2002; Fritz , 2004). Highly supercritical slide velocity impacts generated an impact crater with the water displacement significantly exceeding the slide volume (Fritz , 2003a; Fritz , 2003b; Fritz , 2009; and Weiss , 2009). The relative wave amplitude is defined with the Froude number and relative slide thickness. A large-scale pneumatic landslide tsunami generator launched granular landslides into the 3D tsunami wave basin at Oregon State University (Mohammed and Fritz, 2012). Three-dimensional large-scale experiments on subaerial landslide generated tsunamis were conducted for planar slope, fjord, and conical island scenarios (McFall and Fritz, 2016; Kim , 2020). Experimental results including wave gauge recordings and PIV surface reconstruction were compared against field events including the 2007 Chehalis Lake landslide tsunami (McFall and Fritz, 2017; McFall , 2018). Laboratory experiments based on bottom uplift as tsunami generation mechanism were first introduced with a vertical piston motion in a two-dimensional flume (Hammack and Segur, 1974). Three-dimensional bottom uplift experiments have been limited to small scales with the presence of capillary effects in the wave generation and propagation processes (Jamin , 2015; Shen , 2021). Large-scale three-dimensional laboratory experiments can contribute to the understanding of the tsunami generation process, as well as the advancement and validation of numerical models.

Cylindrical water waves were studied theoretically and experimentally. The Laplace equation with the fully nonlinear free surface boundary conditions can be solved using the third-order approximation of the Dini expansions and iteration method to study the sloshing in a cylindrical containers (Mack, 1962). The Boussinesq equations in the cylindrical coordinates are solved numerically to simulate axial symmetrical waves propagating over an arbitrary bottom along the radial direction (Chwang and Wu, 1977). With the water depth decreasing along the direction of wave convergence and decreasing r, the cylindrical solitary wave amplitude grows much faster than rate of the $r − 1 / 2$ law. An improved kernel representation of one-dimensional long-wave dynamics on a uniform slope could be extended to study cylindrical long-wave propagation (Shimozono, 2020). The cylindrical standing wave can cause wave breaking and jet formation through the convergence of concentric waves at the center of a circular wave basin (McAllister , 2022). Various spike formations can be found in many fields of fluid dynamics, including the bubble bursting at the free surface and breaking wave (Longuet-Higgins, 1983; Duchemin , 2002; Singh and Das, 2021; and Nie , 2023), Worthington jets after water entry and cavity collapse (Gekle and Gordillo, 2010), air pulse jets (Ghabache , 2014), and microscale jets (Huang , 2022; Taraki and Said Ismail, 2022).

A few wave generation mechanisms have been applied to generate water waves in wave basins to model the tsunami due to landslides and bottom deformation. However, no large-scale physical model with controllable cylindrical wave generation has been implemented to investigate submarine volcanic eruption generated tsunamis. A new wave generation method is carried out to model the submarine volcanic eruption-induced tsunami and gain insight into the relevant wave generation mechanisms. First, the novel volcanic tsunami generator (VTG) and experiment setup information are presented in Sec. II. The DIC method and water surface reconstruction are discussed in Sec. III. The wave generation processes and water surface spike are recorded by linear potentiometer and high-speed cameras and discussed in Sec. IV. Finally, the relevant results are summarized in Sec. V.

The experiments were specifically designed to study tsunami waves generated by underwater volcanic eruptions. The wave maker is a novel volcanic tsunami generator with pneumatic actuators, which enables controlled vertical acceleration. The laboratory setup, including the VTG, the instrumentation bridge, and the wave gauges in the wave basin, is shown in Fig. 1. The volcanic tsunami generator experiments were performed in the tsunami wave basin at the O. H. Hinsdale Wave Research Laboratory at Oregon State University.

The layout of the tsunami wave basin, shown in Fig. 1(a), features the deployment of the volcanic tsunami generator (VTG) and resistance wave gauges including vertical wire gauges, cantilever wave gauges and runup gauges. The dimensions of the tsunami wave basin are 48.8 × 26.5 m² with the filled water depth ranging from 0.73 to 1.5 m for the VTG experiments. The featured nearfield wave gauge No.4 at $r=1.73 m$ is installed in the wave basin to measure the tsunami wave profile generated by the VTG. The wave gauges were strategically positioned to measure the radial propagation, dispersion, and decay of the generated cylindrical waves.

The volcanic tsunami generator is designed as a single stage telescopic column driven by eight FESTO pneumatic cylinders supplied by air tanks. The initial and end positions of the VTG vertical stroke are shown in Fig. 2. This design includes eight pneumatic cylinders with 0.080 m diameter and 0.305 m stroke length, which allows to cover a wide spectrum of motion patterns. The generalized Froude similarity is applied to the volcanic eruption and tsunami generation based on combinations of water depths and vertical eruption velocity controlled by the driving tank pressure. The spike Froude number $F s p= v p g h s p$ is defined as the dimensionless eruption velocity, which relates the vertical eruption velocity $v p$ of the VTG to the shallow water depth wave celerity $g h s p$⁠. The Froude scaling allows for comparison with real-world events from slow mud-volcanoes to explosive eruptions. The nearfield tsunami waves are measured as the variation of the water surface $η$ from the still water level with resistance wave gauges.

The axial symmetry of the wave propagation was considered during the deployment and the angular independence checked during data analysis. The wave gauge measurements can characterize the radial wave propagation. This physical model can be considered as an ideal problem of the generation and propagation of a cylindrical wave in water of uniform depth h. It represents a simplified volcanic tsunami model isolating the mechanical eruption process as well as a unique benchmark for validating mathematical models of the cylindrical wave propagation and runup.

The motion of the VTG and the deformation of the free surface are recorded by multiple cameras. In total, eight cameras were deployed above and under water including two PIV cameras. Two Delta Vision HD underwater cameras were installed in the tsunami wave basin on the side opposing the deployment of the wave gauges. They were oriented laterally toward the VTG to record the vertical column motion and its interaction with the surrounding water body. Two Panasonic PTZ cameras mounted on the bridge were oriented toward the VTG to observe the water surface during the tsunami generation. Videos from both the PTZ and underwater cameras were acquired and recorded in the recorder Sound Device PIX 260i.

A pair of LaVision Imager Pro X cameras was deployed in a stereo optical setup on top of the bridge to enable 3D surface reconstruction of the water surface above the eruptive column to characterize the initial condition of the wave source. The Edgertronic high-speed camera with 1280 × 1024 frame size and 495 fps is installed on the sidewall of the wave basin to record the free surface spike motion. The cameras were calibrated before experiment trials to enable stereo optical surface reconstruction and obtain quantitative data from the recordings.

The calibration process ensures accurate projection from the real-world dimensions to the image pixel dimensions. The calibration process is carried out with two plywood boards as shown in Fig. 3. The above water calibration plate (3.66 × 2.44 m²) and the underwater one (2.44 × 1.22 m²) were painted black with white calibration dots. The dots on the larger calibration board were 38 mm in diameter and 300 mm apart. The dots on the calibration board for the underwater cameras were 38 mm in diameter and 200 mm apart. Figure 4 shows the above water calibration plate analyzed by dot recognition and stereo-calibration algorithms in DaVis 8.0 software (LaVision Inc.).

The digital image acquisition, stereo PIV analysis, and surface reconstruction are performed with the FlowMaster PIV and StrainMaster DIC packages in the DaVis 8.0 software. Synchronized water surface image pairs are acquired at a frame rate of 13 fps in the experiments. The deformation and velocity of the water surface are determined with an iterative multi-pass cross-correlation-based DIC algorithm applied to the tracer particle pattern on the water surface.

The DIC algorithm is applied only to the generation and propagation of the leading wave due to the initial distribution of surface tracer particles and discontinuous in subsequent waves, which ensures the necessary texture for image reconstruction (McFall , 2018). The surface reconstruction algorithms include several approaches like cross-correlation analysis, fractional window offset, and iterative multi-grid processing with window distortion (Westerweel, 1997; Fincham and Delerce, 2000; and Scarano and Riethmuller, 2000). The surrounding light reflections and other instruments in the camera view can interfere the DIC algorithm applied for the surface reconstruction. The measuring area of 9  $m 2$ around the VTG is processed with the DIC algorithm. The pixel size is 0.002 m in the surface reconstruction. The cross-correlation analysis is conducted on the dyed water surface by customizing the StrainMaster DIC package. A standard cross-correlation interrogation is performed with a relatively large interrogation window size of 128 to 128 pixels and a mean initial window shift. In the next round iteration of the velocity and deformation, the calculated velocity field is applied as the initial displacement for the multi-pass cross-correlation algorithm. The velocity prediction determines the window shift for the higher resolution level with a refined interrogation window size. The iteration is repeated until the final window size is reached. The DIC surface reconstruction can only be applied on image sequences with a sufficient density of tracer particles in the water surface region to identify the peak shift in the correlation domain.

The uncertainty in experimental measurements includes systematic and random errors. The errors in the surface reconstruction measurements can be summarized as $ε tot= ε ν+ ε optics$⁠, where $ε ν$ is the random error and $ε optics$ is the optical imaging error. The optical image error generally comes from the image recording, rectification, and calibration process. Additional optical image errors in the DIC method can be due to the out-of-plane motion of the water surface. The maximum elevation change at the center of VTG is approximately $0.15 m$⁠. This corresponds to an estimated $2.5%$ variation in the observation distance between the water surface and the camera position. The uncertainty in the surface reconstruction has been investigated based on several benchmark cases with DIC algorithm (Garcia , 2002; Hwang , 2013), and the correlation algorithm is accurate to 1/100 pixel and the grid extraction algorithm accuracy is 1/30 pixel. Considering the joint effects of multiple error sources, the vertical accuracy in the surface reconstruction with this approach is estimated at around 5%.

The water surface reconstruction is performed for select daylighting scenarios. Given the erupted VTG column height of 1.03 m, experimental runs with still water depth h > 1.03 m are fully submerged and the shallower ones with h < 1.03 m are surface piercing. Among the submerged runs, the cases between 1.03 and 1.27 m are considered as daylighting scenarios as the thickness of the water layer submerging the volcano is comparable to the amplitude of the first wave trough above the VTG top surface. The domes are only sustained in the submerged scenarios (h > 1.03 m) due to the existence of VTG column in the surface piercing scenarios. With increasing water depth, the corresponding dome height decreases. Therefore, the daylighting scenario is the most critical case in terms of wave generation. The basin water is dyed white and 0.25-in.-diameter polypropylene plastic spheres are used as tracer particles as shown in Fig. 5. Due to the outward spreading of tracer particles during the VTG motion, the surface reconstruction is limited to initial image sequences with a suitable density of tracer particles in the water surface region. The identification of the peak shift in the correlation domain is only valid with adequate tracer particle patterns as shown in Figs. 5(a) and 5(b). The surface reconstruction fails after the tracer particles have been pushed outward by the VTG motion as shown in Fig. 5(c).

Six sample frames of the surface reconstruction at water depth h = 1.17 m (⁠ $d h=1.03, d=1.20 m$⁠) are shown in Fig. 6. The surface reconstruction results are radially interpolated given the axial symmetry of the wave generation and propagation, which has been confirmed with the wave gauges and PTZ camera recordings. The 120° sector of the axisymmetric pattern is applied to the interpolation. The center of the circular dome profile is first determined for each frame, and the interpolation is calculated as the average of all actual reconstructed data points at a given radial coordinate. The averaging process during interpolation also increases the accuracy of the axisymmetric dome elevation. The smooth dome shape is observed during the wave generation, which is followed by the rebounding and the leading trough at the source. The dome profile reaches peak height in 0.5 s and the peak height at the center is about 0.15 m, which is 50% percent of the stoke length. The correlation between the stroke length and the dome height is matched with the rapid motion regime introduced in Jamin (2015).

The parameters defining the wave generation and propagation are illustrated in Fig. 7 with the dimensions of the VTG and water depth definitions, where the VTG column diameter d is 1.2 m and the vertical stroke length s is 0.305 m. The generated waves are propagating outward in axisymmetric pattern. The leading wave generated by VTG eruption can be classified as a leading elevation N-wave based on video observations and wave gauge recordings. Circular wave breaking around the VTG is observed at the erupted platform and attributed to radial currents and mantel surface reflections.

The vertical displacement is recorded with a FESTO linear encoder type potentiometer (MLO-POT) installed in the VTG. The resulting velocity time series is found as the time derivative of the recorded displacement. The stroke velocities are controlled by the driving tank pressure settings and water depths in experimental trials. The peak velocity $v p$ is applied to define the Froude number $F= v p/ g h$ and characterize the initial condition of wave generation.

In both completely submerged and surface piercing scenarios, a major peak value followed by lower amplitude fluctuation appears in each velocity time series as shown in Fig. 8. The peak value is positively correlated with the driving tank pressure at a given water depth. The fluctuation observed in the latter half of the motion is generated by transient air pressures on the ventilating and exhaust sides of the moving pistons as well as the shock absorbers installed at the end of the stroke to prevent equipment damage due to impact at the cylinder heads. The displacement and velocity time series are compared and verified with underwater camera recordings.

Under certain water depths from h = 1.10 to 1.23 m, a concentric vertical spike can be observed during bore formation after the rising phase of motion. The spike is formed by concentric superposition of the first inward propagating wave crest. The surface spike shows up as a vertical spike with stable axial symmetry during the life cycle from formation to scattering. At the later stage of the spike projectile motion, the initially continuous jet of the spike breaks up into fragments. The spike eventually scatters into droplets beginning from the top. This process is mainly controlled by the velocity of the VTG and counteracting gravity as well as capillary effects during the fragmentation. At laboratory scale, the top of the rough spike at h = 1.17 m water depth can exceed 3 m above the water surface. At larger water depth from h = 1.20 to 1.23 m, the spike pattern becomes continuous and smooth compared to the cases under relatively shallower water depth.

The above and under water cameras capture the characteristics of the surface spike. There are several different patterns of the surface spike as shown in Fig. 9, which are controlled by the spike Froude number $F s p= v p g h s p$ and the relative spike VTG diameter $D s p= d h s p$⁠, where $v p$ is the peak velocity of the VTG, $d=1.2 m$ is the VTG column diameter, and $h s p=h− h c−s=h−1.03 m$ is the submerged depth between the still water surface and the top surface of the erupted VTG. At a certain water depth, there may be a transition between various patterns under different Froude numbers.

The time trajectory of the surface spike during the wave generation process is synchronized with the VTG velocity profile and wave gauge recordings as shown in Fig. 10. The wave time series measured by the wave gauge is a leading elevation N-wave with axial symmetry. The wave generator has completed its stroke at t = 0.9 s and the dome in the free surface at source is formed. The free surface above the VTG column forms a dome shape and reaches its maximum height at the center when the wave generator stopped. After the wave generator motion stopped, the hump begins to collapse into a trough at t = 1.75 s. Once the trough is formed, the inward propagation of the leading wave collides at the center of the wave source at t = 2.3 s, and the spike flow is thus formed. The spike shoots up aggressively and the main jet reaches its top height at t = 2.55 s. The spike is formed by superposition of the first inner propagation wave crest during the bore formation.

There are several types of surface spikes depending on the wave generator motion and water depth, which can be classified according to Fig. 11. The spikes appear in a limited range of water depth in the experiment. The pattern of the surface spike is very sensitive on the relative water depth and spike Froude number as shown in Fig. 11. Considering the spike generation process by VTG motion and the scaling law of Worthington jets proposed by Ghabache (2014), an empirical VTG motion parameter $γ= F s p 0.25 D s p$ is defined to represent both initial spike tip velocity and spike height. The parameter $γ$ is based on the Froude number similarity and the condition that the spike height and spike tip velocity is controlled by inertia and gravity in the laboratory scale. At a certain water depth, there may be a transition between two patterns under different VTG motion parameter $γ$ defined with spike Froude number $F s p$ and relative spike VTG diameter $D s p$⁠.

The first pattern is a splash spike at $D s p=17.1, γ∈[15.8,17.0]$ where the driving tank pressure and the VTG motion parameter are low. At $D s p=17.1, γ∈ 17.7 , 21.0$ and $D s p=10.9, γ∈ 11.8 , 12.5$ with higher tank pressure or larger water depth, the spike pattern fragments and transforms to a splash with relatively low height. The splash pattern is very random in its shape and height, while the vertical spike pattern has a visually matching repeatability.

At $D s p=10.9, γ∈[9.5,11.4]$ and $D s p=8.6, γ∈[8.7,9.6]$⁠, the pattern represents a strong and rough spike. This pattern includes the largest spike height among all scenarios, which is around 3 m above the water surface. The height is defined by the continuous water jet, while individual water droplets are sprayed much higher. Scattered droplets were observed but excluded from spike height measurements.

At $D s p=8.6, γ∈[7.6,8.3]$ and $D s p=7.1, γ∈[6.7,7.5]$⁠, the spike can still form a smooth and slower pattern compared to the previous rough scenario. The height of the spike is reduced compared to the previous one.

For scenarios with water depth outside of the range $5.5<γ<21.0$⁠, the water spike eventually reduced to a hump, which is similar to a theoretical model. The concentric collision of a cylindrical solitary wave was studied by Chwang and Wu (1977) with the Boussinesq wave model. This model can predict a superposition of a cylindrical solitary wave. However, the Boussinesq wave model cannot reproduce local violent deformation such as the spike formation. The spikes observed in the experiments belong to the violent free surface flows and, as a singularity, beyond the theoretical predictions.

Figure 12 shows the classification of three different patterns by the relative column diameter and the Froude number. The three patterns are clustered with different ranges of dimensionless VTG parameter $γ= F s p 0.25 D s p$⁠. With increasing $γ$⁠, the spike pattern transforms through the three distinct regimes: smooth spike, rough spike, and splash spike. At extreme large $γ$ where either the water depth is relatively small or the column diameter is relatively large, the inward propagation process is similar to an inward collapsing cylindercal dambreak wave resulting in a random splash as shown in Figs. 10(c) and 10(d). In the experimental trials with $γ<5$⁠, the characteristic spike disappears and a residual hump forms at the center.

The spike formation and wave generation from the cylindrical VTG satisfies the Froude number similarity. Considering the complex processes often superimposed during volcanic tsunami generations, this work focuses solely on the mechanical aspect of the volcanic eruption. The experimental data from this isolated process can serve as benchmark data for numerical tsunami model validation from laboratory to real scales. In the gravity- and inertia-dominated flow regime, this large-scale physical model will deepen understanding of the nonlinear free surface responses generated by the volcanic eruption. The initial tip velocity $v t$ is defined as the spike tip velocity evaluated at the spike generation moment and near the still water level. The initial tip velocity and the spike height are determined by the image tracking of the leading fluid droplet with diameter larger than 10 mm. The initial tip velocity is the time derivative of the spike top height. The spike tip velocity is related to the vertical motion of the wave maker with the spike Froude number defined as $F s p= v p g h s p$⁠, where $h s p$ is the submerged depth between the still water surface and the VTG erupted top surface, and $v p$ is the peak velocity of the VTG. The dimensionless spike height $ζ h s p$ and spike tip velocity $v t g h s p$ depend on the spike Froude number and relative VTG diameter $D s p$⁠.

As shown in Fig. 15, the rough spike has the peak value of the dimensionless maximum spike height and spike tip velocity among all tested water depths. The dependence of the maximum spike height on the submerged depth and VTG parameters is similar with the Worthington jets of the relaxation flow in the rough spike scenario. However, the transition from rough spike to smooth spike shows the VTG motion Froude number $F s p$ is also critical. The maximum spike height as well as spike tip velocity in the smooth spike scenario is positively correlated with the VTG peak velocity $v p$ at certain submerged depth $h s p$⁠.

Large-scale physical model experiments were performed to study tsunami waves generated by underwater volcanic eruptions. The volcanic tsunami generator (VTG) is designed as a single stage telescopic column driven by eight pneumatic cylinders supplied by internal and external air tanks, which enables controlled vertical acceleration. The VTG was installed on the floor of the tsunami wave basin and completely submerged. The vertical displacement of the wave generator is recorded by the linear encoder potentiometer installed inside the VTG and connecting the base with the uplifting cap. The optical and gauge measurements of the free surface elevations are deployed to record the near-source volcanic tsunami wave. In the tsunami wave basin, eight cameras were deployed above and under water including two PIV cameras in a stereo-optic configuration to capture the nearfield surface deformation.

The wave characteristics are described by the VTG Froude number, relative wave maker stroke, submerged depth, relative wave maker diameter, and the radial propagation distance. The wave generator performances are evaluated with the nearfield wave recordings. A concentric vertical spike is observed in a specific range of water depth, which is generated by superposition of the circular bore from the wave generator motion. The spike may be clustered with different ranges of a dimensionless similarity parameter. With variation of the spike parameter, the pattern transforms through three distinct categories: smooth spike, rough spike, and splash spike. The maximum dimensionless spike height appears in the rough spike regime, which can reach 30 times the submerged water depth on top of the erupted VTG. The measured spike tip height and velocity are compared with the existing theories of Worthington jets, and the spike under various categories has different correlation with the VTG motion parameters. The characteristics of the rough spikes are matched with the typical Worthington jets, while the smooth spikes are more dependent on the wave generator velocity instead of the submerged depth. The splash spikes seem to lose the axial symmetry after generation process and become more random and complicated than other scenarios.

The deformation of free surface and the dynamics of spike formation are suited for further investigations with sophisticated numerical and theoretical methods. The experimental data reported in this paper can serve as a benchmark to validate and improve three-dimensional volcanic tsunami models.

This work was supported by the National Science Foundation (NSF), under ENH Award Nos. CMMI-1563217 and CMMI-1519679. Support and assistance from the staff of the O.H. Hinsdale wave research laboratory, Oregon State University, the Natural Hazards Engineering Research Infrastructure (NHERI), and the NSF REU students are acknowledged. Supporting data for this paper can be found on the DesignSafe-CI web-based platform at www.designsafe-ci.org.

The authors have no conflicts to disclose.

Yibin Liu: Data curation (lead); Formal analysis (lead); Investigation (equal); Methodology (equal); Writing – original draft (lead); Writing – review & editing (equal). Hermann Fritz: Conceptualization (lead); Funding acquisition (lead); Investigation (equal); Methodology (equal); Supervision (lead); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.