# Scenario-based tsunami hazard assessment at Kolumbo submarine volcano

Alessandro Tadini, Matteo Cerminara, Raphaël Paris, Augusto Neri, R.

Stephen J. Sparks, Georges Vougioukalakis, Anna Koutroulli, Benedetta Calusi

► To cite this version:

Alessandro Tadini, Matteo Cerminara, Raphaël Paris, Augusto Neri, R. Stephen J. Sparks, et al. Scenario-based tsunami hazard assessment at Kolumbo submarine volcano. Bulletin of Volcanology, 2025, 87 (52), 10.1007/s00445-025-01837-w. hal-05120126

HAL Id: hal-05120126

https://hal.science/hal-05120126v1

Submitted on 19 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Scenario-based tsunami hazard assessment at Kolumbo submarine volcano

Alessandro Tadini<fnref n="1" /> Ⓒ · Matteo Cerminara<fnref n="1" /> · Raphaël Paris<fnref n="2" /> · Augusto Neri<fnref n="1" /> · R. Stephen J. Sparks<fnref n="3" /> · Georges Vougioukalakis<fnref n="4" /> · Anna Koutroulli<fnref n="4" /> · Benedetta Calusi<fnref n="5" />

Received: 23 December 2024 / Accepted: 18 May 2025
© The Author(s) 2025

## Abstract

Volcanic-induced tsunamis have a potentially devastating impact, especially in densely populated and/or touristic coastal areas. Kolumbo submarine volcano (Greece) experienced in 1650 CE an explosive eruption with eyewitnesses' accounts of major tsunamis along the coasts of Santorini (Thera) and other islands. We present a scenario-based tsunami hazard assessment at this volcano based on existing simulations from literature and new simulations of tsunamis triggered by a less investigated but important mechanism, i.e., submarine landslides on the volcano flanks or within its crater. Simulation results show that the remobilization of a landslide volume of 150–300 Mm<fnref n="3" /> inside the crater can produce tsunami waves larger than 10 m high along the NE coast of Thera and of the order of 5 m along the E and SE coasts. The expected tsunami arrival time ranges from 2 to 3 min along the NE coast of Thera up to 8–10 min on its SE coast. Such scenarios produce inundation areas consistent with those reconstructed for the 1650 CE event and tsunami waves propagating inland at velocities from 2 to 12 m/s. Simulation results also suggest that, given the landslide parameters assumed, it is unlikely that such a large volume landslide would mobilize from the SW-facing Kolumbo crater slopes given the relatively gentle topo-bathymetry of this area. The study findings are part of a broader project which includes a structured expert judgment on different hazard-related topics, including tsunami. The outcomes of this latter, which indicate that the chances of having waves larger than 1 m high on the NE coast of Thera have median probabilities of 50–60%, are relevant to the results of the present study.

**Keywords** Kolumbo Central Volcano (Greece) · 1650 CE eruption · Tsunami hazard · Numerical simulation

## Introduction

Despite not being a frequent event, volcanogenic tsunamis represent one of the deadliest phenomena potentially associated with volcanic eruptions (Auker et al. 2013; Paris

et al. 2019; Schindelé et al. 2024), during which eruptive, volcano-tectonic, and gravitational phenomena can lead to the generation of tsunamis (Paris 2015). Tsunami hazard is potentially high for submarine volcanoes in shallow-water conditions (Day 2015; Syamsidik et al. 2020; Terry et al. 2022), where tsunamis can be triggered by explosions, the gravitational collapse of the subaerial eruptive column, pyroclastic density currents entering the sea, caldera collapse, or submarine flank instability (landslides). Moreover, violent explosive eruptions in shallow-water may also trigger meteo-tsunamis as observed during the 2022 Hunga Tonga Hunga Ha'apai eruption (Shen et al. 2024).

A shallow-water condition is met for Kolumbo, an active submarine volcano located ~ 7 km NE of Santorini (Thera in Greek language) island within the Aegean Sea (Fig. 1) which last erupted in 1650 CE (Vougioukalakis et al. 1996, 2025; Cantner et al. 2014). This explosive eruption has been characterized by several hazardous phenomena, the main ones being tephra fallout (Fuller et al. 2018), gas poisoning

Editorial responsibility: M. Edmonds

<figure><img src="image_2.png" /></figure>

Alessandro Tadini
alessandro.tadini@ingv.it

<footnote n="1">¹ Istituto Nazionale Di Geofisica E Vulcanologia, Sezione Di Pisa, Pisa, Italy</footnote>

<footnote n="2">² Laboratoire Magmas Et Volcans, CNRS, Université Clermont Auvergne, Aubière, IRD, OPGC, France</footnote>

<footnote n="3">³ School of Earth Sciences, University of Bristol, Bristol, UK</footnote>

<footnote n="4">⁴ Hellenic Survey of Geology and Mineral Exploration (HSGME), Athens, Greece</footnote>

<footnote n="5">⁵ Dipartimento Di Matematica E Informatica “Ulisse Dini, Università Degli Studi Di Firenze, Florence, Italy</footnote>
**Fig. 1** Location of the Santorini archipelago (Thera, Therasia, Nea Kameni, Palea Kameni, and Aspronisi) and the Kolumbo Central Volcano. Part of the Kolumbo Volcanic Field (which extends further to the NE) is enclosed within the red dashed ellipse. Coordinates expressed in the GGRS'87/Greek Grid- EPSG:2100 UTM projection. The inset in the upper left corner shows the location of the Santorini archipelago and the Kolumbo Central Volcano with respect to mainland Greece and the nearby islands (Ios, Amorgos, and Anafi)

<figure><img src="image_3.png" /></figure>

(Konstantinou 2020), and tsunami (Nomikou et al. 2014b; Ulvrova et al. 2016; Katsigera et al. 2024). Building on the heavy consequences that this eruption had on the population living in Thera island and on the proximity of Kolumbo Central Volcano to this highly touristic island, a comprehensive hazard assessment project has been carried out for Kolumbo Central Volcano under the auspices of the Greek civil protection authorities (Sparks et al. 2025). A core aspect of this project has been the setup of a performance-based structured

expert judgment (SEJ) exercise (Bevilacqua et al. 2025) to quantify the major uncertainties associated with specific aspects of the reconstruction of the 1650 CE eruption but also with potential hazards, including tsunami, associated with future explosive eruptions of Kolumbo Central Volcano. In brief, during performance-based SEJ elicitation, a group of experts is firstly calibrated through a set of appropriate questions and weighted according to their capability of correctly estimating uncertainties on their answers to

<figure><img src="image_4.png" /></figure>
those questions. Then, the answers of the experts to other questions of interest with respect to the problem under investigation are pooled using each expert's weight: The result is, for each question of interest, a group synthesized distribution (defined through three percentiles, usually the 5th, 50th, and 95th) that represents the uncertainty around each question. More details on this approach could be found in Cooke (1991), Aspinall (2006), and Bevilacqua et al. (2025).

The aim of this paper is to present a scenario-based tsunami hazard assessment for Kolumbo Central Volcano based on a review of existing modeling studies and new numerical simulations that were carried out to inform the expert elicitation target questions on tsunami hazard. The previous literature using numerical simulations includes a study by Ulvrova et al. (2016), investigating a variety of possible mechanisms of tsunamis in 1650 CE (including pyroclastic flows, underwater explosions, and caldera collapse) and a study by Karstens et al. (2023) on combining landslide and explosion mechanisms to reconstruct this event. These studies, together with the outcomes of our new simulations, can be used to investigate the tsunami hazard for future eruptions.

In particular, our new simulations aimed to investigate the dynamics of tsunamis generated by the, so far, less studied mechanism, that is, the formation of a submarine landslide on the volcano flanks or within the crater. The new simulations can describe the coupled dynamics of landslides and tsunamis thus illustrating the effect of the source process on the tsunami propagation. Simulation results also provide new insights on the potential of submarine landslides, in particular of those occurring inside the crater, to generate tsunamis that reach the coast of Thera, as well as on the interpretation of the genesis of the 1650 CE tsunamis. Modeling results are also compared to and discussed in relation to the tsunami-related outcomes of the expert elicitation procedure developed in the project (see Bevilacqua et al. 2025, target questions of block Q14).

This paper is organized as follows: We begin with a review of (i) the Kolumbo Central Volcano and the 1650 CE eruption with particular focus on tsunami observations (“The Kolumbo Central Volcano and its 1650 CE eruption” section) and (ii) the previous tsunami simulations (“Previous tsunami simulations” section). Then, we illustrate our simulation strategy (“Materials and methods” section), present the results of the new simulations of submarine landslides and associated tsunamis and relate these to the outcomes of the expert elicitation questions (“Results” section), and discuss the implications of the simulation results (both the new ones and those from literature) for the tsunami hazard assessment (“Discussion” section). Finally, we summarize the main findings and some future perspectives in the “Conclusions” section.

# Background

## The Kolumbo Central Volcano and its 1650 CE eruption

The Kolumbo Central Volcano is part of the Kolumbo Volcanic Field, situated ~10 km NE of Thera island (Fig. 1). A review of the Kolumbo Central Volcano and the related volcanic field is reported in Vougioukalakis et al. (2025). Here, we summarize features important in the context of tsunami hazard. The polygenetic and active Kolumbo Central Volcano is the southernmost edifice of the Kolumbo Volcanic Field (Fig. 2a) and covers a total area of ~29 km². The large central crater, with an approximate perimeter of 7.75 km, ranges in water depth from 18 m on the southwest side to the deepest point at 205 m in the N-NW sector. The inner crater wall has a mean slope of 29°, ranging from >15° to over 50° (Fig. 2b). The outer slopes of the pumice cone are in contrast lower than the inner ones (<10°) but show some scars (particularly to the E-SE sides) which could be interpreted as landslide scarps.

Several seismic profiles across the Kolumbo crater floor show outward-dipping reflectors beyond the crater rim within the pumice cone and inward-dipping reflectors (Karstens et al. 2023). These reflectors are caused by internal stratification within the pumice deposits of 1650 CE eruptions (Vougioukalakis et al. 2025). Reflectors representing crater infill deposits are also observed beneath the crater floor, and slope analysis in Vougioukalakis et al. (2025) also indicates that mass wasting deposits have partially infilled the crater. Failures of the crater walls are indicated too by the scalloped morphology of the crater rim of the pumice cone. Karstens et al. (2023) documented a basal detachment surface on the NW flank of the Kolumbo Central Volcano, based on detailed bathymetry and seismic data.

ROV observations (Cantner et al. 2014; Nomikou et al. 2014b) revealed two regions. The upper region consists of a thick (up to ~250 m thick) upper sequence of highly stratified pumiceous submarine pyroclastic deposits in the crater walls, which has been attributed to the 1650 CE eruption. The 1650 CE pumice deposits that make up the main cone edifice are weakly consolidated (Cantner et al. 2014) and prone to failure. The deposits of this latter are also observed to depths of between 275 and 290 m on the western, northern, and eastern flanks inside Kolumbo crater. Below that depth, the crater walls are covered by loose tephra scree derived from the instability of the pumice cone inner walls. Apart from these mostly unconsolidated deposits, lava, necks, and dykes occur on the western side beneath the 1650 CE deposit, mostly between the depths of 400 and 450 m. In terms of volume, the 1650 CE deposit is estimated as 2.5 km³ DRE (with an uncertainty interval between 0.9 and 4.9 km³; Vougioukalakis et al. 2025).
**Fig. 2** a Slope distribution of the topography and bathymetry of Santorini with indication of its main villages or settlements (yellow dots), ports, and the coastline (blue line). The submarine Kolumbo Central Volcano is shown on the NE side of Thera, about 7 km away from its coast, and it is also zoomed in (b). Slope angle values are derived from the topo-bathymetry described in the “Topo-bathymetry” section. The computational domain used in the numerical simulations is evidenced by the red box. Coordinates expressed in the GGRS’87/Greek Grid-EPSG:2100 UTM projection

<figure><img src="image_5.png" /></figure>

The Kolumbo Central Volcano experienced its most recent eruption from September to December 1650 CE with a major explosive eruption on 29 September. This powerful eruption formed a 500-m deep crater with a diameter of 1700 m (Fig. 2b; Nomikou et al. 2012).

A detailed review of the available evidence and literature together with an analysis of the chronicles of the 1650 CE eruption of Kolumbo Central Volcano is reported in Vougioukalakis et al. (2025). In the following, we will briefly recall the main facts, with specific reference to the tsunamis.

Most of the observations of the 1650 CE eruption and its related phenomena were compiled by Fouqué (1879) from earlier reports and descriptions of the phenomena. Between January 1649 CE and March 1650 CE, precursory earthquakes (with an intensity of VII on Mercalli scale) were felt on Thera Island. In September (14–26), subterranean roaring and green seawater attested to the beginning of submarine volcanic activity, building up a submarine pumice cone. Then, the intensity of the earthquakes increased, with several events reaching an intensity of VIII (Mercalli). Earthquakes persisted throughout the eruption, but their influence on the generation of the largest observed tsunami during the paroxysmal stage is unclear, as also evidenced by numerical simulations (Ulvrova et al. 2016). The islet emerged from seawater on September 26, although its nature is not yet identified (Vougioukalakis et al. 2025). During the following days, phreatomagmatic explosions frequently broke the sea surface with intermittent jets of gas and ash. The sea was progressively covered

with pumice. The following shallow-water phreatomagmatic phase (September 27–29) did not generate major tsunamis, but on September 27, a wave pushed a boat into the sea and then brought it back to the shore.

The paroxysmal phase started in the early morning of September 29 and lasted up to a maximum of 2 days. The major explosions were heard as far away as 400 km in the Dardanelles. Fine ash was deposited up to western Turkey (150–200 km). Gas affected Thera’s inhabitants (20–50 deaths) and animals. Earthquakes accompanying the paroxysmal phase were felt in Crete (120–140 km away), and one of the earthquakes was strong and long enough to cause damage to buildings in Thera. The largest tsunami occurred at some time on September 29 during this paroxysmal phase. Available testimonies did not report the timing and number of tsunamis; thus, it is difficult to associate the tsunami with a particular source mechanism. Waves were observed on the coast of Thera, where about 2 km² of land was eroded (revealing Hellenistic and Byzantine ruins at Kamari and Perissa), many trees were uprooted, and five churches were destroyed (Dominey-Howes et al. 2000). The spatial distribution of the tsunami deposits allows an estimation of the minimum wave runup in the order of 15–20 m a.s.l. on the eastern coast of Thera (> 630 m inundation near Monolithos settlement) and > 3.5 m a.s.l. on the southern coast (> 360 m inundation near Perissa) (Ulvrova et al. 2016). The coast inside the Thera caldera was not affected by the tsunami, consistent with numerical simulations (Ulvrova et al. 2016). The timing of the tsunami in the framework of the eruption chronology

<figure><img src="image_6.png" /></figure>
remains unclear, even if historical sources place it around nighttime.

The tsunami also impacted nearby islands, especially Ios and Sikinos, located more than 20 km northwest of the volcano. On Ios, a local wave runup of 14–20 m a.s.l. associated with pumice deposition was observed on a rocky shore, but the precise location remains unknown. On Sikinos, the tsunami penetrated 240 m inland (Fouqué 1879: “350 pas”, erroneously translated as feet in many publications). Sea agitation and damage to ships are mentioned in Kea Island (150 km northwest of Kolumbo) and Crete (towns of Dia and Chania, 120 and 165 km to the south, respectively).

Between October 1 and December 20, the activity was characterized by periodic explosions (waning phase), gas release (e.g., 20 deaths on November 4), and small tsunamis whose origin is unknown (possibly following a large explosion on November 4 and renewed activity on December 6). On October 2, nine men approaching the eruption site by boat were found swollen and burnt.

## Previous tsunami simulations

All numerical simulations published so far aimed to understand the possible source mechanism of the main tsunami observed at some time on September 29 during the paroxysmal phase of the 1650 CE eruption (Ulvrova et al. 2016; Karstens et al. 2023). The full list of parameters of the simulations presented in these two papers is provided in Table 1. Ulvrova et al. (2016) tested three different source mechanisms of tsunami:

a. Submarine explosions with energies ranging from $3 \times 10^{14}$ to $5.4 \times 10^{16}$ J (plausible range of energies for shallow-water volcanic explosions). The depth of the explosion was set at 150 m, which corresponds to the approximate present-day rim of the crater.

b. Caldera collapses of different geometries (full collapse of the central crater, collapse of the upper half only, or deepening/collapse of the lower half) and different durations (from 1 min to 1 h). Recent examples of caldera collapse (e.g., Pinatubo 1991; Hunga Tonga 2022) suggest that this source mechanism takes at least 30 min (Scott et al. 1996; Gupta et al. 2022).

c. Pyroclastic density currents resulting from the gravitational collapse of the subaerial eruptive column, with different flow densities (from 1100 to 1500 kg/m$^3$), different flow velocities (from 5 to 30 m/s), different volumes (from $5 \times 10^6$ to $100 \times 10^6$ m$^3$), and different volume fluxes (from $10^4$ to $10^7$ m$^3$/s).

Karstens et al. (2023) proposed a hybrid scenario, with a submarine landslide of the north-western flank of the

volcano (thus accounting for a possible initial retreat of the sea on the eastern coast of Thera), followed by a powerful submarine explosion. They tested different flow densities (from 1250 to 1750 kg/m$^3$) and yield strength (from 5 to 10 kPa) for the landslide, but only one volume (1.2 km$^3$). Explosion energies range from $3 \times 10^{14}$ to $2.2 \times 10^{16}$ J, following the methodology proposed by Ulvrova et al. (2016). They proposed that the slow movement of the landslide on the north-western flank of the volcano (of about 500–1000 m occurring in about 4 min) caused a major depressurization of the magmatic system and the associated major explosion.

Ulvrova et al. (2016) also traced the limits of sedimentary deposits from the 1650 CE tsunami, giving an idea of the extent of the inundation. The waves thus reached minimum altitudes ranging between 3.5 m a.s.l. (Perissa, southern coast) and 20 m a.s.l. (Monolithos, eastern coast), corresponding to a minimum inundation of 360 and 630 m, respectively.

## Materials and methods

The new simulations of submarine landslides and associated tsunamis have been carried out using the Multilayer-HySEA model (Fernández-Nieto et al. 2018; Macías et al. 2020a, b; Esposti Ongaro et al. 2021). We investigated in detail three scenarios: the case of a “small-scale” landslide occurring along the outer SW flank of Kolumbo Central Volcano described by simulation 1 with a volume of about 1.8 Mm$^3$ and two “large-scale” slope failures inside its central crater described by simulations 2 and 3 with volumes of about 147 and 300 Mm$^3$, respectively.

## Numerical model

The Multilayer-HySEA code is a multilayer, non-hydrostatic model in which the three-dimensional model equations are depth-averaged across a number of vertical layers. The governing equations correspond to a semi-discretization of the vertical variables of the Euler equations. The total pressure is decomposed into a sum of hydrostatic and non-hydrostatic pressures. In this process, the horizontal and vertical velocities are assumed to have a constant vertical profile in each layer. The model ensures an exact energy balance, and as the number of layers increases, its linear dispersion relation progressively converges toward that of Airy's theory, the most accurate wave dispersion relationship derived directly from the linearization of the Navier-Stokes equations (Fernández-Nieto et al. 2018). Moreover, the Multilayer-HySEA model can simulate the two-way interaction of the tsunamis with a landslide: In this case, the motion of the bottom surface is represented by an additional layer described by the shallow-water equations of granular material. In this way, the
<table><caption>**Table 1** Parameters used in simulations performed by Ulvrova et al. (2016) and Karstens et al. (2023) to reproduce the 29 September 1650 CE tsunami event</caption><thead><tr><th rowspan="2">Source type</th><th rowspan="2">Location</th><th colspan="3">Input parameters</th><th rowspan="2">Reference</th></tr><tr><th>Explosion depth (m)</th><th>Explosion energy (J)</th><th>Initial surface displacement (m)</th></tr></thead><tbody><tr><td>Submarine explosion</td><td>Central crater</td><td>150</td><td>3 × 10<sup>14</sup></td><td>50</td><td>Ulvrova et al. 2016</td></tr><tr><td>Submarine explosion</td><td>Central crater</td><td>150</td><td>5.6 × 10<sup>14</sup></td><td>100</td><td>Ulvrova et al. 2016</td></tr><tr><td>Submarine explosion</td><td>Central crater</td><td>150</td><td>3.2 × 10<sup>15</sup></td><td>150</td><td>Ulvrova et al. 2016</td></tr><tr><td>Submarine explosion</td><td>Central crater</td><td>150</td><td>9.8 × 10<sup>15</sup></td><td>200</td><td>Ulvrova et al. 2016</td></tr><tr><td>Submarine explosion</td><td>Central crater</td><td>150</td><td>2.6 × 10<sup>16</sup></td><td>250</td><td>Ulvrova et al. 2016</td></tr><tr><td>Submarine explosion</td><td>Central crater</td><td>150</td><td>5.4 × 10<sup>16</sup></td><td>300</td><td>Ulvrova et al. 2016</td></tr><tr><td>Submarine explosion</td><td>Central crater</td><td>150</td><td>3 × 10<sup>14</sup></td><td>100</td><td>Karstens et al. 2023</td></tr><tr><td>Submarine explosion</td><td>Central crater</td><td>150</td><td>1.5 × 10<sup>15</sup></td><td>125</td><td>Karstens et al. 2023</td></tr><tr><td>Submarine explosion</td><td>Central crater</td><td>150</td><td>3.2 × 10<sup>15</sup></td><td>150</td><td>Karstens et al. 2023</td></tr><tr><td>Submarine explosion</td><td>Central crater</td><td>150</td><td>2.2 × 10<sup>16</sup></td><td>240</td><td>Karstens et al. 2023</td></tr><tr><td></td><td></td><td>Volume (m<sup>3</sup>)</td><td>Flow density (kg/m<sup>3</sup>)</td><td>Flow velocity (m/s)</td><td>Volume flux (m<sup>3</sup>/s)</td></tr><tr><td>Pyroclastic flow</td><td>Radial—all flanks</td><td>5 × 10<sup>6</sup></td><td>1500</td><td>20</td><td>1 × 10<sup>7</sup></td></tr><tr><td>Pyroclastic flow</td><td>Radial—all flanks</td><td>5 × 10<sup>6</sup></td><td>1500</td><td>20</td><td>1 × 10<sup>7</sup></td></tr><tr><td>Pyroclastic flow</td><td>Radial—all flanks</td><td>5 × 10<sup>6</sup></td><td>1500</td><td>20</td><td>1 × 10<sup>7</sup></td></tr><tr><td>Pyroclastic flow</td><td>Radial—all flanks</td><td>10 × 10<sup>6</sup></td><td>1500</td><td>20</td><td>1 × 10<sup>7</sup></td></tr><tr><td>Pyroclastic flow</td><td>Radial—all flanks</td><td>10 × 10<sup>6</sup></td><td>1500</td><td>5</td><td>1 × 10<sup>7</sup></td></tr><tr><td>Pyroclastic flow</td><td>Radial—all flanks</td><td>10 × 10<sup>6</sup></td><td>1500</td><td>20</td><td>1 × 10<sup>7</sup></td></tr><tr><td>Pyroclastic flow</td><td>Radial—all flanks</td><td>10 × 10<sup>6</sup></td><td>1500</td><td>30</td><td>1 × 10<sup>7</sup></td></tr><tr><td>Pyroclastic flow</td><td>Radial—all flanks</td><td>50 × 10<sup>6</sup></td><td>1500</td><td>20</td><td>1 × 10<sup>7</sup></td></tr><tr><td>Pyroclastic flow</td><td>Radial—all flanks</td><td>50 × 10<sup>6</sup></td><td>1500</td><td>20</td><td>1 × 10<sup>7</sup></td></tr><tr><td>Pyroclastic flow</td><td>Radial—all flanks</td><td>50 × 10<sup>6</sup></td><td>1500</td><td>30</td><td>1 × 10<sup>7</sup></td></tr><tr><td>Pyroclastic flow</td><td>Radial—all flanks</td><td>50 × 10<sup>6</sup></td><td>1500</td><td>20</td><td>1 × 10<sup>7</sup></td></tr><tr><td>Pyroclastic flow</td><td>Radial—all flanks</td><td>50 × 10<sup>6</sup></td><td>1500</td><td>30</td><td>1 × 10<sup>7</sup></td></tr><tr><td>Pyroclastic flow</td><td>Radial—all flanks</td><td>100 × 10<sup>6</sup></td><td>1500</td><td>20</td><td>1 × 10<sup>7</sup></td></tr><tr><td>Pyroclastic flow</td><td>Radial—all flanks</td><td>100 × 10<sup>6</sup></td><td>1500</td><td>20</td><td>1 × 10<sup>7</sup></td></tr><tr><td>Pyroclastic flow</td><td>Radial—all flanks</td><td>100 × 10<sup>6</sup></td><td>1500</td><td>20</td><td>1 × 10<sup>7</sup></td></tr></tbody></table>

<figure><img src="image_7.png" /></figure>
<table><caption>**Table 1 (continued)**</caption><thead><tr><th rowspan="2">Source type</th><th rowspan="2">Location</th><th colspan="3">Input parameters</th><th rowspan="2">Reference</th></tr><tr><th>Explosion depth (m)</th><th>Explosion energy (J)</th><th>Initial surface displacement (m)</th></tr><tr><th></th><th></th><th colspan="3">Collapse model</th><th></th></tr><tr><th></th><th></th><th colspan="3">Collapse duration (min)</th><th></th></tr></thead><tbody><tr><td>Caldera collapse</td><td>Central crater</td><td>full collapse</td><td>1</td><td></td><td>Ulvrova et al. 2016</td></tr><tr><td>Caldera collapse</td><td>Central crater</td><td>full collapse</td><td>5</td><td></td><td>Ulvrova et al. 2016</td></tr><tr><td>Caldera collapse</td><td>Central crater</td><td>full collapse</td><td>10</td><td></td><td>Ulvrova et al. 2016</td></tr><tr><td>Caldera collapse</td><td>Central crater</td><td>full collapse</td><td>20</td><td></td><td>Ulvrova et al. 2016</td></tr><tr><td>Caldera collapse</td><td>Central crater</td><td>full collapse</td><td>60</td><td></td><td>Ulvrova et al. 2016</td></tr><tr><td>Caldera collapse</td><td>Central crater</td><td>upper half</td><td>1</td><td></td><td>Ulvrova et al. 2016</td></tr><tr><td>Caldera collapse</td><td>Central crater</td><td>upper half</td><td>5</td><td></td><td>Ulvrova et al. 2016</td></tr><tr><td>Caldera collapse</td><td>Central crater</td><td>upper half</td><td>10</td><td></td><td>Ulvrova et al. 2016</td></tr><tr><td>Caldera collapse</td><td>Central crater</td><td>upper half</td><td>20</td><td></td><td>Ulvrova et al. 2016</td></tr><tr><td>Caldera collapse</td><td>Central crater</td><td>lower half</td><td>1</td><td></td><td>Ulvrova et al. 2016</td></tr><tr><td>Caldera collapse</td><td>Central crater</td><td>lower half</td><td>2</td><td></td><td>Ulvrova et al. 2016</td></tr><tr><td>Caldera collapse</td><td>Central crater</td><td>lower half</td><td>10</td><td></td><td>Ulvrova et al. 2016</td></tr><tr><td>Caldera collapse</td><td>Central crater</td><td>lower half</td><td>20</td><td></td><td>Ulvrova et al. 2016</td></tr><tr><td></td><td></td><td>Landslide volume (m³)</td><td>Flow density (kg/m³)</td><td>Yield strength (kPa)</td><td></td></tr><tr><td>Landslide</td><td>NW flank</td><td>1200 × 10<sup>6</sup></td><td>1250</td><td>5</td><td>Karstens et al. 2023</td></tr><tr><td>Landslide</td><td>NW flank</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>5</td><td>Karstens et al. 2023</td></tr><tr><td>Landslide</td><td>NW flank</td><td>1200 × 10<sup>6</sup></td><td>1750</td><td>5</td><td>Karstens et al. 2023</td></tr><tr><td>Landslide</td><td>NW flank</td><td>1200 × 10<sup>6</sup></td><td>2000</td><td>5</td><td>Karstens et al. 2023</td></tr><tr><td>Landslide</td><td>NW flank</td><td>1200 × 10<sup>6</sup></td><td>1250</td><td>7.5</td><td>Karstens et al. 2023</td></tr><tr><td>Landslide</td><td>NW flank</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>7.5</td><td>Karstens et al. 2023</td></tr><tr><td>Landslide</td><td>NW flank</td><td>1200 × 10<sup>6</sup></td><td>1750</td><td>7.5</td><td>Karstens et al. 2023</td></tr><tr><td>Landslide</td><td>NW flank</td><td>1200 × 10<sup>6</sup></td><td>2000</td><td>7.5</td><td>Karstens et al. 2023</td></tr><tr><td>Landslide</td><td>NW flank</td><td>1200 × 10<sup>6</sup></td><td>1250</td><td>10</td><td>Karstens et al. 2023</td></tr><tr><td>Landslide</td><td>NW flank</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>10</td><td>Karstens et al. 2023</td></tr><tr><td>Landslide</td><td>NW flank</td><td>1200 × 10<sup>6</sup></td><td>1750</td><td>10</td><td>Karstens et al. 2023</td></tr></tbody></table>
<table><caption>**Table 1 (continued)**</caption><thead><tr><th rowspan="2">Source type</th><th rowspan="2">Location</th><th colspan="3">Input parameters</th><th rowspan="2">Reference</th></tr><tr><th>Explosion depth (m)</th><th>Explosion energy (J)</th><th>Initial surface displacement (m)</th></tr></thead><tbody><tr><td>Landslide</td><td>NW flank</td><td>1200 × 10<sup>6</sup></td><td>2000</td><td>10</td><td>Karstens et al. 2023</td></tr><tr><td>Landslide</td><td>NW flank</td><td>1200 × 10<sup>6</sup></td><td>1250</td><td>20</td><td>Karstens et al. 2023</td></tr><tr><td>Landslide</td><td>NW flank</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>20</td><td>Karstens et al. 2023</td></tr><tr><td>Landslide</td><td>NW flank</td><td>1200 × 10<sup>6</sup></td><td>1750</td><td>20</td><td>Karstens et al. 2023</td></tr><tr><td>Landslide</td><td>NW flank</td><td>1200 × 10<sup>6</sup></td><td>2000</td><td>20</td><td>Karstens et al. 2023</td></tr><tr><td>Landslide</td><td>NW flank</td><td>1200 × 10<sup>6</sup></td><td>1250</td><td>50</td><td>Karstens et al. 2023</td></tr><tr><td>Landslide</td><td>NW flank</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>50</td><td>Karstens et al. 2023</td></tr><tr><td>Landslide</td><td>NW flank</td><td>1200 × 10<sup>6</sup></td><td>1750</td><td>50</td><td>Karstens et al. 2023</td></tr><tr><td>Landslide</td><td>NW flank</td><td>1200 × 10<sup>6</sup></td><td>2000</td><td>50</td><td>Karstens et al. 2023</td></tr><tr><td></td><td></td><td>Landslide volume (m<sup>3</sup>)</td><td>Flow density (kg/m<sup>3</sup>)</td><td>Yield strength (kPa)</td><td>Explosion energy (J)</td><td>Time delay (min)</td></tr><tr><td>Landslide + explosion</td><td>NW flank + crater</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>7.5</td><td>3 × 10<sup>14</sup></td><td>1</td></tr><tr><td>Landslide + explosion</td><td>NW flank + crater</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>7.5</td><td>1.5 × 10<sup>15</sup></td><td>1</td></tr><tr><td>Landslide + explosion</td><td>NW flank + crater</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>7.5</td><td>3.2 × 10<sup>15</sup></td><td>1</td></tr><tr><td>Landslide + explosion</td><td>NW flank + crater</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>7.5</td><td>2.2 × 10<sup>16</sup></td><td>1</td></tr><tr><td>Landslide + explosion</td><td>NW flank + crater</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>7.5</td><td>3 × 10<sup>14</sup></td><td>2</td></tr><tr><td>Landslide + explosion</td><td>NW flank + crater</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>7.5</td><td>1.5 × 10<sup>15</sup></td><td>2</td></tr><tr><td>Landslide + explosion</td><td>NW flank + crater</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>7.5</td><td>3.2 × 10<sup>15</sup></td><td>2</td></tr><tr><td>Landslide + explosion</td><td>NW flank + crater</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>7.5</td><td>2.2 × 10<sup>16</sup></td><td>2</td></tr><tr><td>Landslide + explosion</td><td>NW flank + crater</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>7.5</td><td>3 × 10<sup>14</sup></td><td>3</td></tr><tr><td>Landslide + explosion</td><td>NW flank + crater</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>7.5</td><td>1.5 × 10<sup>15</sup></td><td>3</td></tr><tr><td>Landslide + explosion</td><td>NW flank + crater</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>7.5</td><td>3.2 × 10<sup>15</sup></td><td>3</td></tr><tr><td>Landslide + explosion</td><td>NW flank + crater</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>7.5</td><td>2.2 × 10<sup>16</sup></td><td>3</td></tr><tr><td>Landslide + explosion</td><td>NW flank + crater</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>7.5</td><td>3 × 10<sup>14</sup></td><td>4</td></tr><tr><td>Landslide + explosion</td><td>NW flank + crater</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>7.5</td><td>1.5 × 10<sup>15</sup></td><td>4</td></tr><tr><td>Landslide + explosion</td><td>NW flank + crater</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>7.5</td><td>3.2 × 10<sup>15</sup></td><td>4</td></tr><tr><td>Landslide + explosion</td><td>NW flank + crater</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>7.5</td><td>2.2 × 10<sup>16</sup></td><td>4</td></tr></tbody></table>

<figure><img src="image_8.png" /></figure>
<table><caption>**Table 1 (continued)**</caption><thead><tr><th rowspan="2">Source type</th><th rowspan="2">Location</th><th colspan="3">Input parameters</th><th rowspan="2">Reference</th></tr><tr><th>Explosion depth (m)</th><th>Explosion energy (J)</th><th>Initial surface displacement (m)</th></tr></thead><tbody><tr><td>Landslide + explosion</td><td>NW flank + crater</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>7.5</td><td>3 × 10<sup>14</sup></td><td>5</td><td>Karstens et al. 2023</td></tr><tr><td>Landslide + explosion</td><td>NW flank + crater</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>7.5</td><td>1.5 × 10<sup>15</sup></td><td>5</td><td>Karstens et al. 2023</td></tr><tr><td>Landslide + explosion</td><td>NW flank + crater</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>7.5</td><td>3.2 × 10<sup>15</sup></td><td>5</td><td>Karstens et al. 2023</td></tr><tr><td>Landslide + explosion</td><td>NW flank + crater</td><td>1200 × 10<sup>6</sup></td><td>1500</td><td>7.5</td><td>2.2 × 10<sup>16</sup></td><td>5</td><td>Karstens et al. 2023</td></tr></tbody></table>

Multilayer-HySEA code incorporates the possibility of simulating the generation of tsunami produced by subaerial or submarine granular landslides. The motion of the landslide is described by a granular landslide model (Fernández-Nieto et al. 2008), in which the Pouliquen and Forterre (2002) friction law is implemented. This law characterizes the dependence of static and dynamic friction coefficients on landslide velocity and thickness through three parameters (three friction angles) $\delta_1$, $\delta_2$, and $\delta_3$ which represent characteristic angles of the material in start and stop conditions ($\delta_1$ and $\delta_3$) and the dependence of the rheology from flow thickness and velocity ($\delta_2$) (Macías et al. 2020b; Esposti Ongaro et al. 2021, 2025). The granular landslide model is in turn weakly coupled with the non-hydrostatic multilayer model. To close the system of equations of the tsunami propagation model, in Macías et al. (2020b), appropriate boundary conditions are imposed on both the water free surface and the granular layers. For water layers, a sponge edge is applied at the vertical boundaries of the computational domain to prevent artificial wave reflections, and a no-slip boundary condition is enforced at the bottom surface, whether it corresponds to the non-erodable bathymetry or the granular layer. No mass transfer is allowed across the water free surface or the granular-water interface. Regarding the granular layer, the friction between the granular slide and the underlying non-erodable surface is parametrized following Pouliquen and Forterre (2002). Momentum exchange at the granular-water interface is governed by a friction term. For more details on the model description and validation tests, we refer the reader to Macías et al. (2020a, b). The Multilayer-HySEA numerical code is designed to run on graphic processing unit (GPU)-accelerated high-performance computing (HPC) architectures (Escalante et al. 2018, 2019). The model in this configuration has been already applied to tsunamis generated at volcanic islands such as Stromboli (Italy) as described in Esposti Ongaro et al. (2021, 2025).

## Topo-bathymetry

For the topo-bathymetry of the study area, we considered a 23 × 24 km domain including the Kolumbo crater and the whole Thera Island (Fig. 1). The employed topo-bathymetry is the result of a merging of two different datasets: (1) the bathymetry of the surroundings of Kolumbo Central Volcano and Thera Island and (2) the topography of Thera Island. The bathymetry was collected during several multi-beam echosounder surveys by R/V Aegeo during 2001, 2006, and 2017 cruises (Alexandri et al. 2003; Sigurdsson et al. 2006; Sakellariou et al. 2010; Nomikou et al. 2014a). The pixel size of the original file is 0.03125 arc minutes (X Y coordinate system WGS 84/EPSG:4326), which corresponds to a pixel size of ~48 × 51 m in the GGRS87' coordinate system. The original file was provided by the Hellenic Center for Marine Research, Hellenic National Oceanographic Data Center. The missing bathymetry of the continental shelf and slope around Santorini was complemented by using the bathymetric maps of Nomikou et al. (2013, 2014a). Bathymetric contours between depths of −5 and −180 m were digitized and edited at 10-m intervals to recreate the shallow bathymetry of the Santorini-Christian complex and the Santorini caldera. The resulting files were extracted using the GGRS87' coordinate system, with a cell size of ~48 × 51 m.

The topography derives from a dataset of 5-m pixel size orthophoto maps, created from color countrywide air photos acquired during the “Large Scale Orthophotos” (LSO) Project (2007–2009). This dataset features a vertical accuracy/ RMSEz of +2.0 m or <3.92 m at a 95% confidence level and uses the GGRS’87/Greek Grid coordinate reference system (EPSG:2100) with TM projection.

The merging of the bathymetry and the topography was performed in three steps: (i) by individually resampling the two datasets to a common cell size of 10 m using the nearest
neighbor method, (ii) by mosaicking the two datasets to a new raster, and (iii) by performing a smoothing of the junction zones between the two original datasets to avoid sharp changes in the topo-bathymetry that could cause instabilities in the numerical code. These steps were performed within the ArcGIS10© platform, resulting in a final topo-bathymetry with a 10-m cell size.

We note that the employed topo-bathymetry is obviously different from that in 1650 CE prior to the tsunami-generating event. The simulations described in the following sections are not specifically aimed at reproducing the 1650 CE tsunami, although they can provide useful insights for the interpretation of the tsunami source mechanism of this event.

## Simulation source conditions

To deal with a realistic source condition and to get an accurate dispersion modeling of the waves, the Multilayer-HySEA model was used by adopting a 3-vertical layer configuration: The lower granular layer, representing the tsunamis-generating landslide, is two-way coupled with water, represented by the remaining two layers.

A crucial aspect of the simulations conducted involves delineating the volume involved in triggering the tsunamis

through the generation of the landslide. This delineation was achieved by intersecting the current bathymetric data of the volcano with an ellipsoidal-shaped volume centered at a specified spatial coordinate and defined by its three semi-axes. This process entails carving out a region from the bathymetric surface representing the potentially unstable granular layer that could lead to landslide formation. Following the delineation of the ellipsoid-shaped region, the simulation further determines the actual mobilized volume primarily on the basis of the friction angle used to model the rheological properties of the granular material constituting the landslide. Importantly, both the amount of mobilized material and its dynamics are not prescribed a priori but are evaluated based on the shallow water equations governing the granular layer and its interaction with the water layers. The proposed methodology allows the investigation of dynamics induced by instabilities within the actual bathymetry, facilitating a deeper understanding of potential tsunami generation mechanisms. While the present set of simulations does not explicitly explore the effect of the internal stratigraphic structure of the Kolumbo cone and flanks or the addition and flow of juvenile material resulting from a volcanic eruption onto the bathymetry, we believe that useful insights on the potential impact of such a generation mechanism can be inferred from the outcomes of these new simulations.

<figure><img src="image_9.png" /><figcaption>**Fig. 3** **a** Plain and (b) vertical views of the ellipsoid (dotted lines) used to carve out the bathymetry. The black arrows indicate the line of sight from the opposite panel. The black solid line in panel (b) represents a schematic topo-bathymetry</figcaption></figure>

<table><caption>**Table 2** Parameters defining the ellipsoid used to carve out the topo-bathymetry for the three simulations performed in the study</caption><thead><tr><th rowspan="3">Simulation</th><th colspan="4">Ellipsoid parameters</th></tr><tr><th>Point L</th><th>Point A</th><th>a</th><th>b</th></tr><tr><th>(m UTM, m UTM, m b.s.l)</th><th>(m UTM, m UTM, m b.s.l)</th><th>(m)</th><th>(m)</th></tr></thead><tbody><tr><td>1</td><td>630,548; 4,040,116; 287.0</td><td>631,830; 4,042,023; 6.9</td><td>1384</td><td>1000</td></tr><tr><td>2</td><td>633,008; 4,042,928; 503.2</td><td>631,830; 4,042,023; 6.9</td><td>1723</td><td>1000</td></tr><tr><td>3</td><td>633,008; 4,042,928; 503.2</td><td>630,585; 4,042,928; 236.3</td><td>3397</td><td>1700</td></tr></tbody></table>
The source ellipsoid is horizontally oriented with semi-axis *a*, *b*, and *c* (Fig. 3). Its position, angle of rotation, and size are defined by (see Table 2): (i) the ($x_c, y_c$) horizontal coordinates of the center C; (ii) the ($x_a, y_a$) coordinates of the point A at the end of horizontal semi-axis *a*; and (iii) the length of horizontal semi-axis *b*. In addition, the vertical coordinate of points C and A ($z_a$) is defined by the elevation of point A on topo-bathymetry plus 20 m of vertical offset to avoid a vertical cut in correspondence of point A. Finally, the vertical semi-axis *c* is prescribed by imposing that the lowest point L (which has the same horizontal coordinates of C) is on the topo-bathymetry: $L(x_c, y_c, z_l)$. Points A and L should be selected so that the former is shallower than the latter. A sketch summarizing the above-mentioned geometrical elements is provided in Fig. 3.

The resulting ellipsoid uniquely intersects the topo-bathymetry, and the intersection space represents the involved volume that could potentially slide, but it does not necessarily coincide with the real landslide volume mobilized. Such a mobilized volume of the landslide is dependent on the chosen rheology, depending in turn on the landslide density and the assumed friction angle (see the next section for more details on the volumes involved).

By varying the above-mentioned parameters (Table 2), we carried out three simulations: simulations 1, 2, and 3. As already mentioned and as we will see in detail in the next subsection, simulation 1 corresponds to a “small-scale” (about 1.8 $Mm^3$ in volume) landslide occurring along the outer SW flank of Kolumbo, whereas simulations 2 and 3 correspond to “large-scale” (of 147 and 300 $Mm^3$ in volume, respectively) landslides generated by slope failures of the internal walls of the Kolumbo crater.

We remark that these three simulations were chosen in order to investigate two potentially dangerous generation mechanisms of tsunamis at Kolumbo, such as (1) a small-scale submarine landslide on the SW external slopes of the volcano directly facing the island of Thera (simulation 1) and (2) a large-scale submarine landslide produced by the internal failure of the crater slopes (simulations 2 and 3). We assumed that a small-scale landslide on the internal slopes would not produce a tsunami large enough to impact on Thera. We anticipate that these three simulations were aimed at having a first quantification of the scale of the tsunamis generated by a potential failure of the external slope of Kolumbo Central Volcano toward Thera and of an internal failure of the crater walls and therefore cannot be considered representative of the most likely tsunamigenic submarine landslides potentially occurring at Kolumbo.

In addition to these values, for all the three simulations, a landslide density of 2000 $kg/m^3$ and friction angles $\delta_1 = 10^\circ$, $\delta_2 = 15^\circ$, and $\delta_3 = 10^\circ$ were considered. The chosen friction angles are assumed equal to those used in similar studies carried out at Stromboli volcano (Italy) (Esposti

Ongaro et al. 2021). Esposti Ongaro et al. (2025) provided a preliminary sensitivity analysis on wave amplitude by varying alternatively $\delta_2$ (from 15 to 32°) and $\delta_1 = \delta_3$ (from 7 to 20°). Results indicate that, by increasing both angles, there is a decrease (up to 25%) in wave amplitude. Following the conservative approach of Esposti Ongaro et al. (2025), the above-mentioned values of $\delta_1 = 10^\circ$, $\delta_2 = 15^\circ$, and $\delta_3 = 10^\circ$ were therefore considered. The density of the porous material underwater should be considered larger than that measured in air because the pores of the lapilli of the landslide and interstitial space between lapilli are partially or totally filled with water. The present study does not investigate the effects of density and friction angle variations on tsunami wave dynamics. However, initial findings on simulations at Stromboli volcano suggest that altering landslide density from 2500 to 1700 $kg/m^3$ leads to a variability in tsunami wave heights that is considerably less significant when compared to the substantial uncertainty associated with parameters such as instability volume and position (Cerminara et al. 2024; Trolese et al. 2024).

Simulations were carried out at INGV-Pisa on the “hpc-gpu” server (equipped with 20 Intel-Xeon CPU cores and 3 NVIDIA P100 GPUs). The computational time to simulate 1000 s of dynamics (corresponding to 16.7 min) is 13 h (for simulation 1), 20 h (for simulation 2), and 29 h (for simulation 3).

## Results

Here, we present the main outcomes of the three simulations. Some refer to the landslide generating the tsunami, whereas others refer to the propagation of the tsunami itself.

In Table 3, we summarize, for each simulation, the volume available for mobilization as defined by intersecting the bathymetry with the source ellipsoid, the actual mobilized volume, and the time of movement of the submarine landslide (i.e., the time necessary for 99% of the mobilized volume to stop). The dynamics of the mobilized volumes are influenced by the shape and

<table><caption>**Table 3** Maximum landslide volume (as defined by intersecting the bathymetry with the source ellipsoid), mobilized volume, and landslide movement time (i.e., the time necessary for 99% of the mobilized volume to stop) for the three simulations</caption><thead><tr><td>Simulation</td><td>Maximum landslide volume (Mm³)</td><td>Actual mobilized volume (Mm³) (percentage with respect to the maximum volume)</td><td>Movement time (s)</td></tr></thead><tbody><tr><td>1</td><td>75</td><td>1.82 (2.5%)</td><td>48</td></tr><tr><td>2</td><td>245</td><td>147 (61%)</td><td>62</td></tr><tr><td>3</td><td>2,200</td><td>300 (14%)</td><td>48</td></tr></tbody></table>
geometry of the source volume, the slope of the bathymetry in the source area, and the rheological properties of the landslide material. In particular, only the portion of the volume with an initial slope steeper than $δ_1$ can be mobilized. Vice versa, the movement times (i.e., how long the landslides keep moving before coming

to rest) are quite similar across the three simulations, despite the very different volume scales and geometries. Such a scale-independent property is mostly related to the adopted rheology of the landslide and represents an interesting topic for future research. Figure 4 shows the evolution in time of the ratio between the integral of the

<figure><img src="image_10.png" /><figcaption>**Fig. 4** Graph showing the evolution in time of the ratio between the integral of mobilized volume at a given time and the total volume remobilized for the three simulations carried out. The asterisks indicate, for each simulation, the time necessary for 99% of the mobilized volume to stop</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 5** Simulation 1: a plain view of the ellipsoid (in blue) with the bathymetry variation (at the end of the landslide propagation) in color scale. The two independent landslides mentioned in the main text are labeled as landslide 1 and landslide 2 and are shown in the zoomed area; b vertical section (red line in panel a) of the bathymetry of the volcano along the main horizontal axis of the ellipsoid carving out the bathymetry (note that the vertical axis of the upper panel of (b)</figcaption></figure>

has a vertical exaggeration; the true scale is provided in the lower panel). The light blue area is a section of the maximum volume of the landslide, whereas the continuous and dashed lines correspond to the profile of the bathymetry at the beginning and end of the simulation, respectively. The orange dot in (b) is the lowest point on the bathymetry intersected by the ellipsoid. Coordinates expressed in the GGRS'87/Greek Grid-EPSG:2100 UTM projection

<figure><img src="image_12.png" /></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 6** Simulation 2: a plain view of the ellipsoid (in blue) with the bathymetry variation (at the end of the simulation) in color scale; b vertical section (red line in panels a, c, and d) of the bathymetry of the volcano along the main horizontal axis of the ellipsoid carving out the bathymetry (note that the vertical axis of the upper panel of (b) has a vertical exaggeration; the true scale is provided in the lower panel). The light blue area is a section of the maximum volume of the landslide, whereas the continuous and dashed lines are the profile of the bathymetry at the beginning and end of the simulation, respectively; plain view of the topo-bathymetry (c) before and (d) after the landslide occurrence. The orange dot in (b) is the lowest point on the bathymetry intersected by the ellipsoid that corresponds to Kolumbo Central Volcano crater bottom. Coordinates expressed in the GGRS'87/Greek Grid-EPSG:2100 UTM projection</figcaption></figure>

mobilized volume at a certain time and the total mobilized volume at the end of the simulation.

In the following subsections, we present, for each simulation, the main results describing the source geometry and landslide dynamics (Figs. 5, 6, and 7 for simulations 1, 2, and 3, respectively) and the tsunami propagation and impact (Figs. 8, 9, and 10, again for simulations 1, 2 and 3, respectively).

In particular, Figs. 5, 6, and 7 describe different views of the initial geometry of the landslide volume as well as the initial and final bathymetry. Vice versa, Figs. 8, 9, and 10 illustrate the tsunami propagation over the whole computational domain by showing (a) the 1-cm wave arrival times (assumed representative of the arrival time of the tsunami), (b) the maximum free surface heights, and (c) the inundation areas (only for simulations 2 and 3) of the Santorini archipelago in terms of water depth values (meters above ground level). Finally, in Fig. 11, we illustrate a comparison between the maximum free surface height and the 1-cm arrival times of the three simulations along the NE to SE coast of Thera, whereas in Fig. 12 and Table 4, we show the same comparison for the maximum free surface height at several control points and related profiles perpendicular to the coastlines.
<figure><img src="image_14.png" /><figcaption>**Fig. 7** Simulation 3: **a** a plain view of the ellipsoid (in blue) with the bathymetry variation (at the end of the simulation) in color scale; **b** vertical section (red line in **a**, **c**, and **d**) of the bathymetry of the volcano along the main horizontal axis of the ellipsoid carving out the bathymetry (note that the horizontal and vertical axes of the upper panel have vertical exaggeration; the true scale is provided in the lower panel). The light blue area is a section of the maximum vol-</figcaption></figure>

ume of the landslide, whereas the continuous and dashed lines are the profile of the bathymetry at the beginning and end of the simulation respectively; plain view of the topo-bathymetry (**c**) before and (**d**) after the occurrence of the landslide. The orange dot in (**b**) is the lowest point on the bathymetry intersected by the ellipsoid that corresponds to Kolumbo crater bottom. Coordinates expressed in the GGRS'87/Greek Grid-EPSG:2100 UTM projection

## Simulation 1

The source geometry and simulation results (Fig. 5a) show a submarine landslide scenario from the outer SW flank of Kolumbo Central Volcano, which is directly facing the NE coast of Thera Island (see Fig. 1). A more direct description of the dynamics of the landslide and associated tsunamis wave is provided by a video animation attached as Supporting Information 1. Despite the large source area and volume encompassed within the ellipsoid, a relatively small volume is mobilized due to the low slope angle of the bathymetry

and consequently the propagation of the landslide is limited (see also Fig. 1). The results show that with this source configuration, two small independent landslides with a combined volume of 1.82 Mm³ are produced: one moving from NE to SW along the Kolumbo external crater walls (landslide 1, Fig. 5a) and another one moving from NW to SE along a relief on the topo-bathymetry (landslide 2, Fig. 5a). This dynamic is clearly visible in the video of Supporting Information 1.

The tsunami propagation and impact results (Fig. 8) reflect the small size of the two submarine landslides,
<figure><img src="image_15.png" /><figcaption>**Fig. 8** Simulation 1: **a** 1-cm wave first arrival times (isolines values are in s). The black asterisk is the point along the NE to SE coast of Thera with the shortest value of arrival time (the green triangle in Fig. 11a); **b** maximum free surface (m a.s.l.), the area with the maximum free surface height > 1.5 m is located on the outer rim of the Kolumbo crater and is evidenced by the black polygon. The</figcaption></figure>

green asterisk is the point along the NE to SE coast of Thera with the highest value of maximum free surface height (the green diamond in Fig. 11b). Coordinates expressed in the GGRS'87/Greek Grid-EPSG:2100 UTM projection. The computational domain is evidenced by the red box

as the maximum free surface elevation is ~1.5 m near the outer rim of the Kolumbo crater (Fig. 8b) and along the NE coast of Thera (Fig. 11b). As a consequence, the resulting inundation for Thera is negligible and, for this reason, has not been illustrated. Arrival times for the 1-cm waves range from ~180 s along the NE coast of Thera to ~540 s near Perissa, along the SE coast of Thera (Figs. 8a and 11a).

## Simulation 2

Simulation 2 presents a potential slope failure inside the central crater of Kolumbo, as the bathymetry variation is all contained within this crater (Fig. 6 and video included as Supporting Information 2).

Such a scenario involves a much larger mobilized volume (about 147 Mm³) and therefore a larger impact on the Santorini archipelago compared to simulation 1. In particular, while the 1-cm wave arrival times do not differ significantly from those computed for simulation 1 (Fig. 9a), the maximum free surface height can reach values up to ~23 m near the Kolumbo crater (Fig. 9b) and up to 17 m near the settlement of Cape Kolumbo's on Thera (see also Fig. 11b). The resulting inundation could affect several areas on Thera Island, such as the NE coast but also large parts of the E coast, especially around the settlement of Monolithos (Fig. 9c).

## Simulation 3

Simulation 3 (Fig. 7) presents a similar scenario to simulation 2, although at a larger scale (mobilized volume now of about 300 Mm³, i.e., about double that of simulation 2). As shown in Fig. 7a, the area enclosing the maximum landslide volume is, in this case, large enough to also include some parts of the outer Kolumbo crater walls, although the collapse is all internal at the crater, as for simulation 2. The video animation of this simulation is provided as Supporting Information 3.

The impact of such a scenario on the Santorini archipelago is the largest among the three simulated. Arrival times of the 1-cm wave do not differ significantly with respect to the other simulations, especially for the NE coast of Thera (Fig. 10a), but maximum free surface heights can now reach > 70 m near the Kolumbo crater (Fig. 10b) and up to about 36 m along the NE coast of Thera (near the settlement of Cape Kolumbo, see also Fig. 11b), resulting in large inundation areas from NE Thera (wide areas on both sides of settlement of Cape Kolumbo) to E Thera (from Monolithos settlement to Kamari) and even to SE Thera (Perissa and southern coast, see Fig. 10c).

## Hazard variables along the E coast of Thera

For a more quantitative analysis of the hazard associated with the simulation outcomes, we also present two additional figures illustrating some key variables along the E coast of Thera:
<figure><img src="image_16.png" /><figcaption>**Fig. 9** Simulation 2: a 1-cm wave first arrival times (isolines values are in s). The black asterisk is the point along the NE to SE coast of Thera with the smallest value of arrival time (the violet triangle in Fig. 11a); b maximum free surface (m a.g.l.), the areas with the maximum free surface height > 20 m are shown by the black polygons. The green asterisk is the point along the NE to SE coast of Thera with the</figcaption></figure>

highest value of maximum free surface height (reported as violet diamond in Fig. 11b); c inundation area with water depth values. Coordinates expressed in the GGRS'87/Greek Grid-EPSG:2100 UTM projection. The computational domain is evidenced by the red box

1. The distributions, along the NE to SE coast of Thera (roughly from the settlement of Cape Kolumbo to Perissa), of the arrival times of the 1-cm wave and

of the maximum free surface height of the tsunami (Fig. 11)

2. A comparison between the maximum surface height of the three tsunami simulations at eight control points

<figure><img src="image_17.png" /></figure>
<figure><img src="image_18.png" /><figcaption>**Fig. 10** Simulation 3: a 1-cm wave first arrival times (isolines values are in s). The black asterisk is the point along the NE to SE coast of Thera with the smallest value of arrival time (the light blue triangle in Fig. 11a); b maximum free surface (m a.s.l.), the area with the maximum free surface height > 70 m is evidenced by the black polygon. The green asterisk is the point along the NE to SE coast of Thera</figcaption></figure>

with the highest value of maximum free surface height (the light blue diamond in Fig. 11b); c inundation area with water depth values. Coordinates expressed in the GGRS'87/Greek Grid-EPSG:2100 UTM projection. The computational domain is evidenced by the red box

and along their associated profiles perpendicular to
the NE-SE coastline of Thera and the NW coastline of
Therasia (Fig. 12)

Regarding Fig. 11, we have chosen to investigate the dis-
tribution of arrival time and maximum surface height along
this profile as it represents the closest area to Kolumbo
<table><caption>**Table 4** Maximum free surface heights and arrival times at the 8 control points of Fig. 12 for the three new simulations carried out in this study</caption><thead><tr><th>Control point number along profile</th><th>Control point depth (m)</th><th>Simulation</th><th>Max free surface height (m)</th><th>1 cm wave arrival time (s)</th></tr></thead><tbody><tr><td rowspan="3">1</td><td rowspan="3">-13.4</td><td>1</td><td>0.07</td><td>424</td></tr><tr><td>2</td><td>1.26</td><td>378</td></tr><tr><td>3</td><td>4.08</td><td>372</td></tr><tr><td rowspan="3">2</td><td rowspan="3">-11.3</td><td>1</td><td>0.28</td><td>172</td></tr><tr><td>2</td><td>3.94</td><td>168</td></tr><tr><td>3</td><td>8.48</td><td>150</td></tr><tr><td rowspan="3">3</td><td rowspan="3">-11.8</td><td>1</td><td>0.51</td><td>140</td></tr><tr><td>2</td><td>5.89</td><td>148</td></tr><tr><td>3</td><td>12.98</td><td>150</td></tr><tr><td rowspan="3">4</td><td rowspan="3">-10.1</td><td>1</td><td>0.27</td><td>196</td></tr><tr><td>2</td><td>2.46</td><td>200</td></tr><tr><td>3</td><td>6.98</td><td>172</td></tr><tr><td rowspan="3">5</td><td rowspan="3">-13.8</td><td>1</td><td>0.28</td><td>268</td></tr><tr><td>2</td><td>2.37</td><td>264</td></tr><tr><td>3</td><td>4.68</td><td>228</td></tr><tr><td rowspan="3">6</td><td rowspan="3">-13.9</td><td>1</td><td>0.11</td><td>394</td></tr><tr><td>2</td><td>1.27</td><td>372</td></tr><tr><td>3</td><td>4.08</td><td>336</td></tr><tr><td rowspan="3">7</td><td rowspan="3">-11</td><td>1</td><td>0.07</td><td>574</td></tr><tr><td>2</td><td>0.97</td><td>520</td></tr><tr><td>3</td><td>4.81</td><td>480</td></tr><tr><td rowspan="3">8</td><td rowspan="3">-17</td><td>1</td><td>0.04</td><td>726</td></tr><tr><td>2</td><td>0.41</td><td>634</td></tr><tr><td>3</td><td>1.95</td><td>594</td></tr></tbody></table>

Central Volcano and the portion of Thera Island (especially from Monolithos toward S) with plain shores and relatively large settlements and critical infrastructures. For the two variables (i.e., arrival time and max free surface height) and the whole E coastline, we also show the minimum of the 1-cm arrival times (triangles in Fig. 11a) and the maximum of the free surface (diamonds in Fig. 11b). These outcomes have relevant applications from a civil protection point of view to identify areas with highest expected waves and shortest time delays between the triggering of the phenomenon at source and the associated effects along the coast.

The travel time between the triggering landslide and the wave arrival at the E coast is between about 3 min, at Koloumbus settlement, and about 10 min, at Perissa. The coastline of Thera Island inside the Santorini caldera is mostly characterized by high cliffs, and it is significantly more protected from large tsunami waves generated by the Kolumbo Central Volcano. However, two important infrastructures are located along this coast; these are the touristic ports of Thera

and Athinios (see Fig. 2). In addition, another touristic port is located in front of the town of Therasia on Therasia Island (see Fig. 2). For these three localities, our largest simulation (simulation 3) indicates.

* 1-cm wave arrival times of ~7 min, ~8 min, and ~6 min at Thera port, Athinios port, and Therasia port, respectively

* Maximum free surface heights of ~0.75 m, ~1.2 m, and ~0.8 m at Thera port, Athinios port, and Therasia port, respectively

The islands inside the Santorini caldera (Nea Kameni, Palea Kameni, and Aspronisi) are not populated and therefore not considered in this analysis, although the density of people can be relatively high in Nea Kameni during tourist visits. In comparison, based on the results of Ulvrova et al. (2016), the tsunami travel times (i.e., the arrival time of the first wave generated at the Kolumbo Central Volcano) for the surrounding islands (see Fig. 1) are of the order of ~6 min for the S coasts of the island of Ios, 8 min for the island of Anafi, and ~10 min for the island of Amorgos.

Figure 12 analyzes the maximum surface height computed in the three simulations in selected locations along the outside Thera coast. The control points selected are similar to those used by Ulvrova et al. (2016). In particular, we drew a profile line passing through each point and perpendicular to the coastline and show the topo-bathymetry (from -30 to 10 m a.s.l.) and the maximum free surface height along it for the three simulations (Fig. 12). For such control points, we also show in Table 4 the corresponding depth from our topo-bathymetry, the 1-cm wave arrival time, and the maximum free height.

To have a better appreciation of the potential hazard associated with tsunami waves generated by Kolumbo Central Volcano, some semi-quantitative comparisons between available data of the 1650 CE event and some selected results of our simulations were performed. For instance, Fig. 13 shows a comparison between the inundation computed by simulation 3 and the field evidence of the 1650 CE tsunami. These latter data have been documented by Ulvrova et al. (2016) on several outcrops along the NE and SE Thera coasts where tsunami deposits have been recognized and linked to the 1650 CE event. The figure clearly shows that simulation results are semi-quantitatively consistent with the historical inundation of the 1650 CE event as reconstructed from field observations. Similarly, Fig. 14 provides an estimation of the potential inundation flow velocities of simulation 3. To get the flow velocity (u), we used the empirical correlation given by Matsutomi and Okamoto (2010), namely

$$u = Fr\sqrt{gh_f} \quad (1)$$
<figure><img src="image_19.png" /><figcaption>**Fig. 11** **a** Arrival times (s) of the 1-cm wave and (b) maximum free surface height (m a.s.l.) along the whole E coast of Thera (highlighted in blue in the small insets). Main towns/settlements along the coast are also indicated with their total inhabitants (2011 census) and some main infrastructures or attractions. The two plots also show the</figcaption></figure>

values of the minimum (for arrival times) and maximum (for free surface height) for the three simulations along the coast (triangles and diamonds, respectively). These latter are also reported, as respectively black and green asterisks, in Figs. 8a, b; 9a, b; and 10a, b

where $h_f$ is the water depth measured at the front of an impacted element (e.g., building, hill), $Fr$ is the Froude number (assumed equal to 0.66), and $g$ is the gravitational acceleration. Values obtained with this relation, based on the largest potential fluid force exerted by the flow on the impacted element, range from 0.2 to 12 m/s highlighting that such flows can produce major devastation of the coastal areas.

## Outcomes of the expert elicitation on tsunami hazard

To have a better appreciation of the values associated with the new tsunami simulation described above, in this section, we briefly report some of the outcomes of the expert elicitation procedure presented in Bevilacqua et al. (2025) with specific reference to questions about tsunamis hazard at Kolumbo (block 14; questions 14a–c, 14d1–14d3). The

tsunami simulations presented here were carried out at the same time as the expert elicitation sessions and therefore the outcome of one informed and guided the other, thus favoring a more robust and consistent description of the phenomenon (Sparks et al. 2025).

The elicitation outcomes, summarized in Table 5, clearly imply that tsunamis are a significant threat for Thera and the other islands around Kolumbo. Moreover, experts’ answers indicate that there is a medium to high probability of having a major tsunami in the circumstances of an explosive eruption occurring in the next 30 years. Tsunamis are more likely to occur from the main crater of Kolumbo, particularly during the paroxysmal phase of the eruption. A complementary result is that tsunamis are less likely to be generated from the other edifices of the volcanic field, although the uncertainty of these estimates is large. In terms of tsunami wave height, median probabilities (for the [5th-50th-95th] percentile values, see Table 5) of having a wave higher than
<figure><img src="image_20.png" /><figcaption>**Fig. 12** Map of the location of the control points (pale blue circles) and associated free surface elevation profiles (purple segments) perpendicular to the Thera coastline, showing the maximum free surface heights as computed for simulations 1, 2, and 3. Vertical red dashed</figcaption></figure>

lines on profile graphs indicate the location of the control points along each profile. Coordinates expressed in the GGRS'87/Greek Grid-EPSG:2100 UTM projection
1 m along the NE coast of Thera are, for the CM (classical model)-EW (equal weight), 52–60%, whereas they decrease to about 19–27% for waves higher than 5 m and to 7–12% for waves higher than 10 m. All these probabilities refer to potential events occurring in a future period of 30 years.

The three new simulations presented above were able to generate maximum wave heights of about 1.5, 1.7, and 36 m (for simulations 1, 2, and 3, respectively, see Fig. 11b) along the NE coast of Thera and values of the order of 1, 5, and 10 m along approximately the central part of the NE coast of Thera (i.e., the area between Cape Kolumbo and Monolithos). As a consequence, the three simulations can be considered approximately representative of the three different scales associated with the elicitation questions (questions 14d1–14d3) and, in turn, to the probabilities of occurrence in case of an eruption in the next 30 years (Table 5).

# Discussion

## New simulated scenarios and associated hazard

The new simulated scenarios allow us to discuss the dynamics of tsunamis triggered by submarine landslides at Kolumbo Central Volcano and their potential in terms of hazard. In this area, most of the modeling efforts so far focused on tsunamis triggered by PDCs, caldera collapse, and underwater explosions (see the “Previous tsunamis simulations” section), while submarine landslides have only been recently considered by Karstens et al. (2023). Our choice to investigate the effect of this generation mechanism is largely dependent on the current geological and bathymetric conditions of the Kolumbo volcanic system, which lacks a subaerial part.

The submarine pumice cone formed in the 1650 CE eruption could indeed become unstable in the case of a future explosive eruption or a regional earthquake such as the 1956 Amorgos earthquake (see e.g. Papazachos and Kkalas 2025). A potential internal failure of the crater rim, as well as a landslide on the external rim of the crater, could be generated during a future explosive eruption (Katsigera et al. 2024) or large-magnitude regional earthquake. In particular, with respect to the simulated volumes of simulations 2 and 3 (i.e., 147 and 300 Mm³), we remark that the Kolumbo inner crater walls are predominantly composed of screees of unconsolidated pumiceous material related to the 1650 CE eruption at about their repose angle.

For the assumed rheological properties of the landslide described above, results of simulation 1 suggest that the SW-facing slope of Kolumbo Central Volcano has low potential to generate submarine landslides resulting in > 1 m tsunami waves on the NE coast of Thera Island. This is interpreted as due to the low slope angles along the SW side of the

Kolumbo Central Volcano, which are generally well below (about 10°) the friction angles (10°–15°–10°) considered in our simulations (see Fig. 1). However, the definition of realistic values of the friction angles of such landslides is still the subject of debate (see, e.g., Poulain et al. 2023). For the Kolumbo case, the definition of appropriate friction angles is furthermore complicated due to the complex stratigraphy of its outer slopes, where fine beds of pumice clasts could act as weak failure surfaces (especially if shaken by a regional earthquake) resulting in an increase in pore pressure and reduced friction angles (see, e.g., d’Acremont et al. 2022; Tran et al. 2024). Therefore, further investigation should be undertaken to evaluate the effect of lower friction angles (< 10°) on the remobilized volume of the landslide and, thus, on its tsunamigenic potential.

Vice versa, simulations 2 and 3, representative of the medium and largest scenarios here considered, describe the propagation of a tsunami triggered by a landslide occurring on the inner slopes of the Kolumbo crater with volumes of about 150 and 300 Mm³, respectively, and have the potential to represent a major hazard for the E coast of Thera (see Figs. 9 and 10). We note that the first wave arrival times at the control points (see Table 4) are ~2–3 min on the NE coast of Thera (e.g., profiles 2–3) and ~8–10 min near Perissa on the SE coast of Thera (profile 7). Moreover, if we consider the whole E coast of Thera, we observe that both the shortest arrival times (Fig. 11a) and the highest waves (Fig. 11b) of our simulated scenarios are expected near the settlement of Cape Kolumbo, but significantly short times (~3–4 min) and high waves (~1–35 m) are also computed along the whole E coast. Such timings and potential wave heights should be properly taken into consideration when preparing emergency/evacuation plans and designing alert systems linked to the detection of the triggering phenomenon.

Finally, although the aim of our simulations was not to reproduce the tsunamis that occurred during the 1650 CE eruption, we note that the internal crater failure scenarios investigated (e.g., simulation 2 and especially simulation 3, see Fig. 10) produce an inundation of the Monolithos/Perissa areas which is semi-quantitatively consistent with the historical accounts (see the “Introduction” section) and the field evidence of the 1650 CE tsunami (see Fig. 13). In addition, we observe that simulation 3 produces an inundation along the E coast of Thera Island comparable with that produced by the simulations of Karstens et al. (2023) but with about 20% of their mobilized volume, a simpler triggering mechanism, and a shorter emplacement time of the mass movement (4 min from Karstens et al. (2023) compared to the 48–62 s of the present simulations, see Table 3 and Fig. 4). With these observations, we are not claiming that the 1650 CE tsunami was certainly generated by an internal collapse of the Kolumbo crater, but this triggering
<figure><img src="image_21.png" /></figure>
►Fig. 13 Comparison between the inundation area (color legend refers to water depth in m a.g.l.) of simulation 3 and locations with evidence of 1650 CE tsunami deposits (white squares, from Ulvrova et al. 2016). (a–c) Enlargements of the main panel. Coordinates expressed in the GGRS’87/Greek Grid–EPSG:2100 UTM projection

mechanism is certainly plausible and compatible with this
event. Further ad hoc simulations should be done to better
reconstruct the 1650 CE tsunami event as well as to forecast
the most likely scenarios in case of tsunamis triggered by
submarine landslides.

Simulated scenarios from literature and hazard
assessment

The considerable number of simulations done in previous
work (Ulvrova et al. 2016; Karstens et al. 2023) and in the
current study can provide an evaluation of tsunami hazard
at Kolumbo. To this aim, we developed a simplified matrix
diagram shown in Fig. 15 that links, from one side, the tsu-
nami-triggering mechanism such as submarine explosion,
pyroclastic density current, caldera collapse, and submarine
landslide and, on the other side, the probabilities of having
major tsunamis with different wave heights (i.e., 1, 5 and 10
m) along the NE coast of Thera in case of an eruption in the
next 30 years as estimated by expert elicitation (Bevilacqua
et al. 2025).

For each of the source mechanisms of tsunami mentioned
in Fig. 15, we identified scenarios for which waves on NE
Thera exceed the values used in the expert elicitation (i.e.,
1, 5, and 10 m). These scenarios are based on the sensitivity
of the simulation outcomes to key source parameters that
have the highest influence on wave amplitude: explosion
energy, pyroclastic density current volume flux, caldera col-
lapse duration, and submarine landslide volume (Paris 2015;
Schindelé et al. 2024). In the matrix diagram of Fig. 15,
all the simulations reported are related to a source location
inside the current Kolumbo crater. Moreover, even consider-
ing the new simulations done assuming such a source, both
the simulations from previous works and those from the pre-
sent one do not reproduce all the scenarios corresponding to
the elicited wave height along NE Thera (i.e., crater internal
failure generating > 1 m wave height and landslide along the
crater external slopes generating > 10 m wave height). These
scenarios should be considered in future developments.

From a hazard perspective, values of tsunami wave
heights at the coast are essential for estimating flow veloci-
ties (c.f. (1)) during the inundation stage. This is critical to
determine how fast humans within the potential footprint can
respond (see Sparks et al. 2025) and to estimate what forces
will be exerted on structures (buildings, harbor piers, etc.).
During the inland inundation phase, tsunami flows typically
have Froude numbers between 0.66 and 2 (Matsutomi and

Okamoto 2010). As an example, a tsunami with initial wave
heights of 1 m, 5 m, and 10 m at the coast will propagate
inland (using (1)) at velocities higher than 2 m/s, 4.5 m/s,
and 6.5 m/s, respectively. For a 1650-like event, inundation
flow velocity in Perissa and Monolithos will likely exceed,
respectively, 4 m/s and 7 m/s. Such estimates are more dif-
ficult to obtain for the Cape Kolumbo area because of the
steep topography in that part of the island. For comparison,
a tsunami with flow velocities > 5 m/s can severely damage
stone and concrete walls with breakage of most columns
(Suppasri et al. 2011).

Conclusions

Based on the findings of the expert elicitation carried out
(Bevilacqua et al. 2025; Sparks et al. 2025), a 30-year future
eruption in the Kolumbo Volcanic Field has the potential to
generate a major tsunami, with high chances of having a >
1 m wave on the NE coast of Thera. The presence of a large
unstable tephra cone and a high probability that the next
eruption at Kolumbo will be dominantly explosive mean that
the volcano meets several conditions conducive to tsunami
formation.

In this paper, we reviewed the results of the tsunami sim-
ulations carried out in previous works (Ulvrova et al. 2016;
Karstens et al. 2023) and new simulations produced during
this project. This analysis allows us to outline a description
of the tsunami scenario–based hazard at Kolumbo. Moreo-
ver, the new simulations carried out in this study allowed
us to describe the dynamic evolution of two new triggering
mechanisms at the Kolumbo crater: a landslide involving
just the outer, SW-facing crater slopes and an internal fail-
ure of the crater rim (simulated by assuming two different
volumes).

The simulation outputs allowed us to obtain the first
important insights regarding hazard associated with tsuna-
mis triggered by submarine landslides. A scenario involv-
ing a significant (i.e., volume of the order of 150–300
Mm³) internal failure of the Kolumbo crater (not necessar-
ily related to an eruption) is capable of producing tsunami
waves up to > 10 m high along the NE coast of Thera. Such a
scenario is made more likely by the steep slopes of the inner
crater wall in agreement with Katsigera et al. (2024). Such
volumes are also associated with wave heights > 5 m along
the E coast of Thera and ~ 5 m along the SE coast of Thera.
Tsunami risk associated with these scenarios is relatively
low along the NE coast of Thera due to the low density of
settlements, whereas it is potentially much larger along the
E and SE coasts, given the presence of main infrastructures
(i.e., Thera airport near Monolithos on the E coast) and
highly touristic areas (near Perissa and Kamari on the E and
SE coasts). The expected arrival times of the first waves are
<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>
►Fig. 14 Comparison between the calculated flow speed (from Eq. 1) of simulation 3 and locations with evidence of 1650 CE tsunami deposits (white squares, from Ulvrova et al. 2016). (a–c) Enlargements of the main panel. Coordinates expressed in the GGRS’87/ Greek Grid–EPSG:2100 UTM projection

also remarkably short—ranging from 2 to 3 min in the most exposed coasts up to 8–10 min on the south-eastern coast of Thera. Such tsunami waves would propagate inland at velocities typically ranging from 2 to 12 m/s. We refer to Sparks et al. (2025) for a first quantitative analysis of the individual and societal risk associated with Kolumbo tsunamis.

We have also shown how simulation results suggest that, given the specific set of landslide rheological parameters assumed, it is unlikely a landslide with a large volume would mobilize just from the SW-facing Kolumbo crater slopes (i.e., without involving the inner crater), given the relatively gentle topo-bathymetry of this area. Nevertheless, more work is needed to investigate the effect of landslide rheological properties on this outcome. In comparison, a scenario involving an internal failure of the Kolumbo crater could

produce an inundation on Thera Island semi-quantitatively consistent with that observed during the paroxysm of the 1650 CE eruption, although we acknowledge that during such a complex eruption, multiple mechanisms could have generated the tsunamis.

Finally, while we underline that the present study can contribute to the improvement of tsunami hazard assessment at Kolumbo by identifying a new effective mechanism of tsunami generation, we believe that more research is needed to fill important knowledge gaps. In particular, additional simulations and analyses are required to investigate tsunami events triggered by landslides from (i) different areas along the outer slopes of Kolumbo crater (e.g., to the NW) including also part of the inner crater, (ii) smaller volumes due to partial internal crater failure, (iii) different areas from the Kolumbo Volcanic Field, and (iv) different physical and rheological properties of the collapsing material (in particular frictional parameters and density). In addition, a larger computational domain including the nearby islands (Ios, Anafi, and Amorgos) is required to better characterize the volcanic tsunami hazard assessment in this region.

<table><caption>**Table 5** Summary of key elicitation outcomes for tsunami-related questions. CM (classical model) and EW (equal weight) decision-maker solutions have been rounded to the first digit with respect to the values provided in Bevilacqua et al. (2025)</caption><thead><tr><th>Question ID</th><th>Content of question</th><th>Numerical outcomes—CM and EW solutions (5th-50th-95th)</th><th>Key elicitation outcomes</th></tr></thead><tbody><tr><td>14a</td><td>Next eruption within the next 30 years (main crater)—major tsunami probability</td><td>CM: [18.7–57.3-92.5]<br>EW: [15.0–64.2-93.4]</td><td>●Medium to high likelihood (~ 60%) of a major tsunami associated with a future eruption, but large uncertainty<br>●Tsunami is a significant threat in the case of a future eruption of Kolumbo main crater</td></tr><tr><td>14b</td><td>Next eruption within the next 30 years (outside main crater)—major tsunami during paroxysmal eruption</td><td>CM: [1.3–9.8-43.4]<br>EW: [1.5–12.4-37.4]</td><td>●Chances of a major tsunami for an eruption outside main crater has medians between 9.8 and 12.4%<br>●Tsunamigenic processes such as explosions and gravitational processes (landslides, pyroclastic flows) are viewed as unlikely in the case of an eruption outside the main crater (i.e., in the Kolumbo volcanic field)</td></tr><tr><td>14c</td><td>Next eruption within the next 30 years (main crater)—tsunami during paroxysmal stage</td><td>CM: [14.3–60.1-89.1]<br>EW: [14.6–70.0-92.0]</td><td>●High likelihood (60–70%) of a tsunami in the case of an explosive paroxysm in the central crater<br>●Low level of uncertainty on this scenario according to the answers of the experts</td></tr><tr><td>14d1</td><td>Next eruption within the next 30 years (main crater)—tsunami wave height > 1 m on the NE coast of Thera</td><td>CM: [15.1–51.7-92.7]<br>EW: [12.5–59.8-94.3]</td><td>●Chances of > 1 m wave on the NE coast of Thera have medians between 51.7 and 59.8%, even if there is a large variability in the answers given by the experts</td></tr><tr><td>14d2</td><td>Next eruption within the next 30 years (main crater)—tsunami wave height > 5 m on the NE coast of Thera</td><td>CM: [4.8–18.8-63.8]<br>EW: [5.2–27.2-74.6]</td><td>●Chances of > 5 m wave on the NE coast of Thera have medians between 18.8 and 27.2%</td></tr><tr><td>14d3</td><td>Next eruption within the next 30 years (main crater)—tsunami wave height > 10 m on the NE coast of Thera</td><td>CM: [1.5–7.4-29.7]<br>EW: [1.3–12.5-56.9]</td><td>●Chances of > 10 m wave on the NE coast of Thera have medians between 7.3 and 12.5%</td></tr></tbody></table>
<figure><img src="image_24.png" /><figcaption>**Fig. 15** Matrix diagram for volcanic tsunami hazard at Kolumbo volcano. Values in square brackets correspond to the [5th-50th-95th] percentiles of the CM decision-maker solutions (rounded to the first decimal digit) of the corresponding questions of the expert elicitation (see Table 5 and Bevilacqua et al. 2025). Colored areas indicate different scenarios produced by different tsunami-triggering mechanisms and report the value of the most important parameters required to</figcaption></figure>

generate waves of a given height on the NE coast of Thera, based on numerical simulations. We remark that each colored area refers just to one specific triggering mechanism which is here considered independent from all the others. “n.a.” indicates that a simulation describing such a specific scenario has not yet been conducted. <sup>a</sup>Ulvrova et al. (2016), <sup>b</sup>Karstens et al. (2023)

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s00445-025-01837-w.

**Acknowledgements** We gratefully acknowledge the EDANYA group from the University of Malaga (Jorge M. Sánchez, Manuel J. Castro, Cipriano Escalante, Sergio Ortega and José M. González-Vida) for developing and sharing the Multilayer-HYSea code (https://edanya.uma.es/hysea/index.php/models/multilayer-hysea) and for their support in its use. The Large Scale Orthophotos (LSO) Project was implemented by the Hellenic Cadastre and co-funded by the European Union within the framework of the Operational Program “Information Society” of the 3rd Community Support Framework and NSRF 2007-2013. We thank two anonymous reviewers and the editor Christopher Gregg for useful comments that improved the quality of this paper.

**Funding** This research has been supported by the project “Hazard and risk assessment for Kolumbo Volcano, Greece,” Hellenic Survey of Geology & Mineral Exploration (HSGME). We thank Hellenic Survey of Geology & Mineral Exploration (HSGME) for hosting the first project workshop on Santorini, and Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Pisa, for hosting the second project workshop.

**Data availability** The online version contains supplementary material available as Supporting Information 1, Supporting Information 2, and Supporting Information 3.

## Declarations

**Competing interests** The authors declare no competing interests.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated

otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

* Alexandri M, Papanikolaou D, Nomikou P (2003) Santorini volcanic field - new insights based on swath bathymetry. Abstracts Iugg 30. Sapporo, Japan

* Aspinall WP (2006) Structured elicitation of expert judgement for probabilistic hazard and risk assessment in volcanic eruptions. In: Statistics in volcanology special publications of IAVCEI No. 1. The Geological Society for IAVCEI, pp 15–30

* Auker MR, Sparks RSJ, Siebert L, Crosweller HS, Ewert J (2013) A statistical analysis of the global historical volcanic fatalities record. J Appl Volcanol 2:1–24. https://doi.org/10.1186/2191-5040-2-2

* Bevilacqua A, de’ Michieli Vitturi M, Tadini A, Neri A, Aspinall WP, Sparks RSJ, (2025) The Kolumbo Volcanic Field, Greece: expert elicitation findings supporting volcanic hazard and risk assessment. Bull Volcanol 87:41. https://doi.org/10.1007/s00445-025-01822-3

* Cantner K, Carey S, Nomikou P (2014) Integrated volcanologic and petrologic analysis of the 1650 AD eruption of Kolumbo submarine volcano, Greece. J Volcanol Geotherm Res 269:28–43. https://doi.org/10.1016/j.jvolgeores.2013.10.004

* Cerminara M, Esposti Ongaro T, De’ Michieli Vitturi M, Tadini A, Trolese M, Fornaciai A, Nannipieri L, Rodriguez Galvez JF (2024) Simulated scenarios of volcanic mass movements and associated tsunamis at Stromboli (Aeolian archipelago, Tyrhenian Sea, Italy). version 1. Istituto Nazionale di Geofisica e Vulcanologia (INGV). https://doi.org/10.13127/stromboli/sciara_del_fuocohardt

* Cooke RM (1991) Experts in uncertainty: opinion and subjective probability in science. Oxford University Press, New York, p 321
d'Acremont E, Lafuerza S, Rabaute A, Lafosse M, Castelot MJ, Gorini C, Alonso B, Ercilla G, Vazquez JT, Vandorpe T, Juan C, Migeon S, Ceramicola S, Lopez-Gonzalez N, Rodriguez M, El Moumni B, Benmarha O, Ammar A (2022) Distribution and origin of submarine landslides in the active margin of the southern Alboran Sea (Western Mediterranean Sea). Mar Geol 445:106739. https://doi.org/10.1016/j.margeo.2022.106739

Day SJ (2015) Volcanic tsunamis. In: Sigurdsson H (ed) The encyclopedia of volcanoes, 2nd edn. Academic Press, pp 993–1009. https://doi.org/10.1016/B978-0-12-385938-9.00058-4

Dominey-Howes DTM, Papadopoulos GA, Dawson AG (2000) Geological and historical investigation of the 1650 Mt. Columbo (Thera Island) eruption and tsunami, Aegean Sea, Greece. Nat Haz 21:83–96. https://doi.org/10.1023/A:1008178100633

Escalante C, de Luna TM, Castro MJ (2018) Non-hydrostatic pressure shallow flows: GPU implementation using finite volume and finite difference scheme. Appl Math Comput 338:631–659. https://doi.org/10.1016/j.amc.2018.06.035

Escalante C, Fernández-Nieto ED, de Luna TM, Castro MJ (2019) An efficient two-layer non-hydrostatic approach for dispersive water waves. J Sci Comput 79:273–320. https://doi.org/10.1007/s10915-018-0849-9

Esposti Ongaro T, de' Michieli Vitturi M, Cerminara M, Fornaciai A, Nannipieri L, Favalli M, Calusi B, Macías J, Castro MJ, Ortega S, González-Vida JM, Escalante C (2021) Modeling tsunamis generated by submarine landslides at stromboli volcano (Aeolian Islands, Italy): a numerical benchmark study. Front Earth Sci 9:628652. https://doi.org/10.3389/feart.2021.628652

Esposti Ongaro T, Cerminara M, de' Michieli Vitturi M, Tadini A, Trolese M, Fornaciai A, Nannipieri L, Calusi B, Macías J, Castro MJ, Escalante C, Ortega S, González-Vida JM, Rodriguez Galvez JF, (2025) Modeling and numerical simulation of volcanic mass movements and associated tsunamis at stromboli (Aeolian archipelago, Tyrrhenian Sea, Italy). Front Earth Sci 13:1548961. https://doi.org/10.3389/feart.2025.1548961

Fernández-Nieto ED, Bouchut F, Bresch D, Castro Díaz MJ, Mangeney A (2008) A new Savage-Hutter type model for submarine avalanches and generated tsunami. J Comput Phys 227:7720–7754. https://doi.org/10.1016/j.jcp.2008.04.039

Fernández-Nieto ED, Garres-Díaz J, Mangeney A, Narbona-Reina G (2018) 2D granular flows with the μ (I) rheology and side walls friction: a well-balanced multilayer discretization. J Comput Phys 356:192–219. https://doi.org/10.1016/j.jcp.2017.11.038

Fouqué F (1879) Santorin et ses éruptions. Masson, Paris, p 588

Fuller S, Carey S, Nomikou P (2018) Distribution of fine-grained tephra from the 1650 CE submarine eruption of Kolumbo Volcano, Greece. J Volcanol Geotherm Res 352:10–25. https://doi.org/10.1016/j.jvolgeores.2018.01.004

Gupta AK, Bennartz R, Fauria KE, Mittal T (2022) Eruption chronology of the December 2021 to January 2022 Hunga Tonga-Hunga Ha'apai eruption sequence. Commun Earth Environ 3:314. https://doi.org/10.1038/s43247-022-00606-3

Karstens J, Crutchley GJ, Hansteen TH, Preine J, Carey S, Elger J, Kühn M, Nomikou P, Schmid F, Dalla Valle G, Kelfoun K, Berndt C (2023) Cascading events during the 1650 tsunamic eruption of Kolumbo Volcano. Nat Commun 14(1):6606. https://doi.org/10.1038/s41467-023-42261-y

Katsigera A, Nomikou P, Pavlopoulos K (2024) A preliminary hazard assessment of Kolumbo Volcano (Santorini, Greece). GeoHazards 5(3):816–832. https://doi.org/10.3390/geohazards5030041

Konstantinou KI (2020) Magma chamber evolution during the 1650 AD Kolumbo eruption provides clues about past and future

volcanic activity. Sci Rep 10(1):15423. https://doi.org/10.1038/s41598-020-71991-y

Macías J, Escalante C, Castro MJ (2020) Multilayer-HySEA model validation for landslide generated tsunamis – part 1: rigid slides. Nat Hazards Earth Syst Sci 21:775–789. https://doi.org/10.5194/nhess-21-775-2021

Macías J, Escalante C, Castro MJ (2020) Multilayer-HySEA model validation for landslide generated tsunamis – part 2: granular slides. Nat Hazards Earth Syst Sci 21:791–805. https://doi.org/10.5194/nhess-21-791-2021

Matsutomi H, Okamoto K (2010) Inundation flow velocity of tsunamis on land. Isl Arc 19:443–457. https://doi.org/10.1111/j.1440-1738.2010.00725.x

Nomikou P, Carey S, Papanikolaou D, Bell KC, Sakellariou D, Alexandri M, Bejelou K (2012) Submarine volcanoes of the Kolumbo volcanic zone NE of Santorini Caldera, Greece. Global Planet Change 90:135–151. https://doi.org/10.1016/j.gloplacha.2012.01.001

Nomikou P, Carey S, Croff Bell K, Papanikolau D, Bejelou K, Alexandri M, Cantner K, Martin J (2013) Morphological analysis and related volcanic features of the Kolumbo submarine volcanic chain (NE of Santorini Island, Aegean Volcanic Arc). Z Geomorphol 57(3):029–104. https://doi.org/10.1127/0372-8854/2013/S-00142

Nomikou P, Parks MM, Papanikolau D, Pyle DM, Mather TA, Carey S, Watts AB, Paulatto M, Kalninš M, Livanos I, Bejelou K, Simou E, Perros I (2014) The emergence and growth of a submarine volcano: the Kameni islands, Santorini (Greece). GeoResJ 1:8–18. https://doi.org/10.1016/j.grj.2014.02.002

Nomikou P, Carey S, Bell K, Papanikolaou D, Bejelou K, Cantner K, Sakellariou D, Perros I (2014) Tsunami hazard risk of a future volcanic eruption of Kolumbo submarine volcano, NE of Santorini Caldera, Greece. Nat Haz 72:1375–1390. https://doi.org/10.1007/s11069-012-0405-0

Papazachos C, Kkallas C (2025) Deterministic seismic hazard assessment for the Kolumbo area. Bull Volcanol (in review)

Paris R (2015) Source mechanisms of volcanic tsunamis. Philos Trans R Soc Lond A 373:20140380. https://doi.org/10.1098/rsta.2014.0380

Paris R, Ulvrova M, Selva J, Brizuela B, Costa A, Grezio A, Lorito S, Tonini R (2019) Probabilistic hazard analysis for tsunamis generated by subaqueous volcanic explosions in the Campi Flegrei caldera, Italy. J Volcanol Geotherm Res 379:106–116. https://doi.org/10.1016/j.jvolgeores.2019.05.010

Poulain P, Le Friant A, Mangeney A, Viroulet S, Fernandez-Nieto E, Castro Diaz M, Peruzzetto M, Grandjean G, Bouchut F, Pedreros R, Komorowski JC (2023) Performance and limits of a shallow-water model for landslide-generated tsunamis: from laboratory experiments to simulations of flank collapses at Montagne Pelée (Martinique). Geophys J Int 233(2):796–825. https://doi.org/10.1093/gji/ggac482

Pouliquen O, Forterre Y (2002) Friction law for dense granular flows: application to the motion of a mass down a rough inclined plane. J Fluid Mech 453:1–19. https://doi.org/10.1017/S0022112001006796

Sakellariou D, Sigurdsson H, Alexandri M, Carey S, Rousakis G, Nomikou P, Georgiou P, Ballas D (2010) Active tectonics in the Hellenic volcanic arc: the Kolumbo submarine volcanic zone. Bull Geol Soc of Greece XLIII/2:1056–1063

Schindelé F, Kong L, Lane EM, Paris R, Ripepe M, Titov V, Bailey R (2024) A review of tsunamis generated by volcanoes (TGV) source mechanism, modelling, monitoring and warning systems.
Pure Appl Geophys 181:1745–1792. https://doi.org/10.1007/s00024-024-03515-y

Scott WE, Hoblitt RP, Torres RC, Self S, Martinez ML, Nillos TJ (1996) Pyroclastic flows of the June 15, 1991, climatic eruption of Mount Pinatubo. In: Newhall CG, Punongbayan S (eds) Fire and mud: eruptions of Mount Pinatubo, Philippines. Philipp. Inst Volcanol and Seismol, Quezon City, pp 545–570

Shen Y, Whittaker CN, Lane EM, White JD, Power W (2024) Physical experiments of waves generated by submerged steam eruptions with applications to volcanic tsunamis. Sci Adv 10(45):eadk6208. https://doi.org/10.1126/sciadv.adk6208

Sigurdsson H, Carey S, Alexandri M, Vougioukalakis G, Croff KL, Roman C, Sakellariou D, Anagnstou C, Rousakis G, Ioakim C, Gogou A, Ballas D, Misaridis T, Nomikou P (2006) Marine investigations of Greece’s Santorini volcanic field. Eos 87(34):337–342. https://doi.org/10.1029/2006EO340001

Sparks RSJ, Vougioukalakis G, Aspinall WP, Neri A, Antonakos A, Baxter PJ, Bevilacqua A, Cerminara M, de’ Michieli Vitturi M, Francalanci L, Koutroulli A, Maeno F, Mastroianni F, Papazachos K, Pardini F, Paris R, Tadini A, Vaselli O (2025) Future eruptions of the Kolumbo volcanic field: prognosis with hazard and risk assessment. Bull Volcanol (in press)

Suppasri A, Koshimura S, Imamura F (2011) Developing tsunami fragility curves based on the satellite remote sensing and the numerical modeling of the 2004 Indian Ocean tsunami in Thailand. Nat Hazards Earth Syst Sci 11:173–189. https://doi.org/10.5194/nhess-11-173-2011

Syamsidik B, Luthfi M, Suppasri A, Comfort LK (2020) The 22 December 2018 Mount Anak Krakatau volcanogenic tsunami on Sunda Strait coasts, Indonesia: tsunami and damage characteristics. Nat Hazards Earth Syst Sci 20:549–565. https://doi.org/10.5194/nhess-20-549-2020

Terry JP, Goff J, Winspear N, Bongolan VP, Fisher S (2022) Tonga volcanic eruption and tsunami, January 2022: globally the most significant opportunity to observe an explosive and tsunamigenic submarine eruption since AD 1883 Krakatau. Geosci Lett 9(1):24. https://doi.org/10.1186/s40562-022-00232-z

Tran QA, Sørlie E, Grimstad G, Eiksund G, Takahashi H, Sassa S (2024) Influence of sediment permeability in seismic-induced submarine landslide mechanism: CFD-MPM validation with centrifuge tests and analysis. Comput Geotech 174:106588. https://doi.org/10.1016/j.compgeo.2024.106588

Trolese M, Cerminara M, Esposti Ongaro T, de’ Michieli Vitturi M, Tadini A (2024) Modeling tsunami generation and propagation: insights from sensitivity analysis of landslide parameters at Stromboli, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-19260, https://doi.org/10.5194/egusp here-egu24-19260

Ulvrova M, Paris R, Nomikou P, Kelfoun K, Leibrandt S, Tappin DR, McCoy FW (2016) Source of the tsunami generated by the 1650 AD eruption of Kolumbo submarine volcano (Aegean Sea, Greece). J Volcanol Geotherm Res 321:125–139. https://doi.org/10.1016/j.jvolgeores.2016.04.034

Vougioukalakis G, Francalanci L, Mitropoulos D, Perissoratis K (1996) The 1649–1650 eruption of the Kolumbo submarine volcanic center, Santorini. In: Second Workshop on European Laboratory Volcanoes, Thira, Santorini, Greece

Vougioukalakis G, Koutroulli A, Sparks RSJ, Francalanci L, Mastroianni F, Laurenzi MA, Schaen G, Papazachos C, Aspinall WP, Kanellopoulos C, Baxter PJ, Neri A, Tadini A, Bevilacqua A (2025) The Kolumbo volcanic field, Greece: a review of evidence to reconstruct the 1650 CE eruption and to support hazard and risk assessment. Bull Volcanol (in press)