# Cascading events during the 1650 tsunamigenic eruption of Kolumbo volcano

Jens Karstens, Gareth J Crutchley, Thor H Hansteen, Jonas Preine, Steven Carey, Judith Elger, Michel Kühn, Paraskevi Nomikou, Florian Schmid, Giacomo Dalla Valle, et al.

► To cite this version:

Jens Karstens, Gareth J Crutchley, Thor H Hansteen, Jonas Preine, Steven Carey, et al. Cascading events during the 1650 tsunamigenic eruption of Kolumbo volcano. Nature Communications, 2023, 14, 10.1038/s41467-023-42261-y. hal-04801322

HAL Id: hal-04801322

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

Submitted on 25 Nov 2024

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
<figure><img src="image_2.png" /></figure>

# Cascading events during the 1650 tsunamigenic eruption of Kolumbo volcano

Received: 27 October 2022

Accepted: 4 October 2023

Published online: 26 October 2023

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

Check for updates

Jens Karstens Ⓔ<fnref n="1" />, Gareth J. Crutchley Ⓔ<fnref n="1" />, Thor H. Hansteen Ⓔ<fnref n="1" />, Jonas Preine Ⓔ², Steven Carey³, Judith Elger Ⓔ<fnref n="1" />, Michel Kühn Ⓔ<fnref n="1" />, Paraskevi Nomikou⁴, Florian Schmid<fnref n="1" />,⁷, Giacomo Dalla Valle⁵, Karim Kelfoun⁶ & Christian Berndt Ⓔ<fnref n="1" />

Volcanic eruptions can trigger tsunamis, which may cause significant damage to coastal communities and infrastructure. Tsunami generation during volcanic eruptions is complex and often due to a combination of processes. The 1650 eruption of the Kolumbo submarine volcano triggered a tsunami causing major destruction on surrounding islands in the Aegean Sea. However, the source mechanisms behind the tsunami have been disputed due to difficulties in sampling and imaging submarine volcanoes. Here we show, based on three-dimensional seismic data, that ~1.2 km³ of Kolumbo's northwestern flank moved 500–1000 m downslope along a basal detachment surface. This movement is consistent with depressurization of the magma feeding system, causing a catastrophic explosion. Numerical tsunami simulations indicate that only the combination of flank movement followed by an explosive eruption can explain historical eyewitness accounts. This cascading sequence of natural hazards suggests that assessing submarine flank movements is critical for early warning of volcanogenic tsunamis.

Volcanogenic tsunamis have caused more than 55,000 fatalities since the late 18th century<fnref n="1" />. Tsunami generation is often due to a combination of processes, including submarine explosions, pyroclastic flows, caldera subsidence, and flank failures². The 2022 Hunga Tonga-Hunga Ha'apai (HTHH) eruption in Tonga³,⁴ and the 2018 Anak Krakatau sector collapse in Indonesia demonstrated the devastating tsunami hazard of volcanic eruptions in shallow marine environments and at the land-sea-boundary⁵. Both eruptions have become a focus of volcanological research in recent years. The 2018 tsunami of Anak Krakatau hit neighboring coasts without any warning and caused hundreds of casualties⁵. Satellite geodetic data from more than ten years before the event indicated deformation within the volcanic edifice⁶. However, it was only after the collapse that these observations were interpreted as revealing prolonged, deep-seated deformation in the buildup to the collapse⁷. Internal deformation of a volcanic flank preceding catastrophic failure has also been demonstrated for the 1888 Ritter Island sector collapse⁸,⁹,

and episodic gravity-driven movement along deep detachment surfaces is known from several large volcanoes¹⁰⁻¹². At Mt Etna, an average seaward motion of 3 to 5 mm per year is observed¹¹, while the southwestern flank of Kilauea shows both a transient cm-scale yearly seaward movement and m-scale slip events accompanied by major earthquakes¹³.

Kolumbo is a fully submerged volcano 7 km northeast of Santorini and part of the Hellenic Volcanic Arc. It is by far the largest volcano of the Kolumbo Volcanic Chain, consisting of 24 volcanic cones formed ~350 kyrs ago¹⁴,¹⁵ (Fig. 1). Kolumbo's edifice was created by at least five eruptive cycles, the earliest dating back more than 1 Myrs¹⁵,¹⁶. The most recent eruption, in 1650 CE, formed a cone consisting of up to ~260 m-thick stratified pumice deposits, which breached the sea surface before being destroyed by a violent explosive eruption that formed a 500-m-deep and 2500-m-wide crater¹⁷,¹⁸. The eruption deposits consist of highly vesicular rhyolitic pumice with cm-size mafic inclusions, likely reflecting an eruption triggered by mafic replenishment into a

<footnote n="1">¹GEOMAR Helmholtz-Zentrum für Ozeanforschung Kiel, Kiel, Germany. ²University of Hamburg, Institute of Geophysics, Hamburg, Germany. ³University of Rhode Island, Kingston, USA. ⁴National and Kapodistrian University of Athens, Athens, Greece. ⁵Italian National Research Council, Institute of Marine Science ISMAR, Bologna, Italy. ⁶Laboratoire Magmas et Volcans, Université Clermont Auvergne, OPGC, CNRS, IRD, F-63000 Clermont, Ferrand, France. ⁷Present address: K.U.M Umwelt und Meerestechnik Kiel GmbH, Kiel, Germany. ☒e-mail: jkarstens@geomar.de</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 1** | Kolumbo volcano in the Aegean Sea. a Map of the Hellenic Volcanic Arc (red) with active volcanoes (stars) and the location of this study area (blue box). b Topographic map of the Christiana-Santorini-Kolumbo volcanic zone and neighboring islands of Sikinos and Ios with historical tsunami eyewitness accounts. c Perspective view of high-resolution (2 m) bathymetry data<sup>44</sup> showing Kolumbo's</figcaption></figure>

deformed northwestern flank with terrace morphology, exposed part of the
detachment surface, and the 1650 crater with morphological expression of dykes.
Eye symbol in (b) shows the viewing direction. The cyan, magenta and purple lines
indicate the position of the profiles depicted in Fig. 3.

cooler silicic reservoir<sup>18,19</sup>. The 1650 eruption triggered a tsunami
causing major destruction on the neighboring islands of Santorini, Ios,
and Sikinos. Historical eyewitness accounts reported maximum run-up
heights of 20 m on the southern coast of Ios, 240 m inundation on
Sikinos<sup>20,21</sup>, and flooding of up to 2 km<sup>2</sup> of land on the eastern coast of
Santorini which caused major destruction in the towns of Kamari and
Perissa<sup>17,20,21</sup> (Fig. 1b). Flooding at Kamari and Perissa was preceded by
sea retreat, indicating a negative-amplitude leading wave<sup>20</sup>. Onshore
tsunami deposits have been identified up to 14 m above sea level on
the north coast of Santorini, 6 m above sea level in Kamari, and up to
360 m inland in Perissa as well as on Ios and Sikinos<sup>20,21</sup>.

The genesis of volcanogenic tsunamis is often complex and may
involve underwater explosions, earthquakes, caldera subsidence,
pyroclastic density currents, flank failures or a combination of these
processes². Even for prominent and comparably recent events, such as
the 1883 eruption of Krakatau, the tsunami source mechanisms are still
debated. Previous analyses of the 1650 tsunami tested the submarine
emplacement of pyroclastic flows, caldera collapse and underwater
explosions as potential singular source mechanisms²⁰.

However, the movement of Kolumbo's northwestern flank (Fig. 1c) has not been considered as a source mechanism in previous simulations. In this study, we use new high-resolution 3D reflection seismic data and numerical tsunami simulations to reconstruct the mechanisms that caused the tsunami waves recorded in 1650. We investigate Kolumbo's flank deformation and evaluate its tsunami- genic potential, as well as its role in triggering the violent explosions. Our analyses show that flank deformation initiated a cascade of events triggering the violent and tsunamigenic eruption of Kolumbo in 1650.

Results

Structural architecture of Kolumbo

High-resolution (2 m) bathymetry data reveal intricate morphological
details of the cone formed by the 1650 eruption of Kolumbo (Fig. 1c).

The northwestern flank of the volcano is characterized by a distinctive terraced morphology, in stark contrast to the smooth flanks to the south and east (Fig. 1c). 3D seismic reflection data reveal the subsurface architecture of this structural asymmetry (Figs. 2 and 3). While the southeastern flank comprises undeformed, sub-parallel strata, the northwestern flank displays pronounced internal deformation (Figs. 2 and 3a). This deformed sedimentary unit lies directly above the -1600 BCE Minoan eruption ignimbrites¹⁶. Aside from the 1650 eruption, there is no evidence of other recent explosive eruptions in sediment cores²²,²³. Therefore, the deformation affects solely volcanic sediments deposited by the 1650 eruption. Mapping of the top and base of the deformed sediments yields a volume of ~1.2 km³, assuming a seismic velocity of 1700 m/s (±50 m/s) based on refraction seismic data¹⁵. The bathymetric data show elongated ridges that strike slope parallel around the northwestern flank of the volcanic edifice (Fig. 1c). The analysis of the 3D seismic data reveals that they are the surface expression of anticlinal folds associated with downslope compression (Figs. 2 and 3). Their curved shapes in planform, their relative steepness, and interpreted internal thrusting indicate that the ridges are the result of compression and not contourite sediment depositional processes seen at some other submarine volcanoes²⁴. Deformation at the toe of the flank is complex and segmented (manifesting itself in the terrace-like morphology), while a distinct detachment surface can be observed toward the top of the remnant cone (hatched region, Figs. 1c and 4a). The exposed detachment surface has a slope angle of ~19°, while the underlying strata dip at ~11°, which is significantly steeper than strata within the southwestern flank that dip at ~6° (Fig. 3a). The steeper stratigraphic dips beneath the northwestern flank are likely due to deposition onto a pre-existing remnant cone (Figs. 2 and 3a), which is not present beneath the southeastern flank. Additionally, there is a morphological step between the two slopes, as the base of the 1650 cone lies 50 - 100 m deeper on the northwestern side of the edifice (Fig. 3). The 1650 eruptive products were deposited on top of
<figure><img src="image_5.png" /></figure>

**Fig. 2 | 3D seismic dataset of the Kolumbo volcano.** Chair-cut view of the 3D seismic data around the volcanic edifice of Kolumbo, showing the deposits of the 1650 eruption (orange semi-transparent overlay) and the crater formed by an explosive eruption. The northwestern flank is characterized by deformation as the

result of lateral displacement along a basal detachment surface (dashed pink line). Folded sediments above the detachment are consistent with compression. Dashed black lines are tentatively interpreted thrust faults associated with the seafloor ridges.

remnants of a cone from a previous eruption (K3 in ref. 15), which introduce structural heterogeneities within the northwestern flank of Kolumbo, in contrast to the relatively flat southeastern flank (Fig. 3). This complex pre-eruption topography, with a steep northwestern flank, likely controlled the instability of this flank segment as slope angle is a dominant factor for slope stability²⁵. Previous numerical simulations²⁰ suggested a violent explosion, pyroclastic flow emplacement, and caldera subsidence as potential tsunami source mechanisms during the 1650 eruption. However, the 3D seismic data show neither indications for sizable pyroclastic flow deposits (Figs. 2 and 3) nor for ring faults that would accommodate caldera collapse (Fig. 3b). Therefore, we deem pyroclastic flow emplacement and caldera collapse as unlikely tsunami source mechanisms. The existence of a 2500-m-wide and 500-m-deep crater, however, provides clear evidence for a highly explosive eruption phase. The seismic and bathymetric datasets suggest a link between this explosion and the flank deformation. At the top of the remnant cone, the detachment surface intersects prominent ridges that are the seafloor expression of feeder dykes (Fig. 4a, d). Their present exposure is probably due to their higher mechanical stability relative to the surrounding volcanoclastic deposits, which have been eroded since the 1650 eruption. While the explosion represents an obvious tsunami source mechanism, the previously unconsidered tsunamigenic potential of the flank deformation requires careful evaluation.

## Reassessment of the 1650 tsunami

Historical eyewitness accounts of tsunamis can be used to test the feasibility of tsunami simulations. For Kolumbo, the observations of

higher tsunami run-up at los than at Santorini and sea retreat preceding inundation in eastern Santorini¹⁷,²⁰,²¹ suggest a northwesterly directed wave propagation. Previous numerical simulations of the explosive eruption assumed a radial-symmetric initial waveform with a steep positive amplitude leading wave, followed by a trough that mimics the crater²⁰. By assuming a peak wave amplitude of 250 m, these simulations could reproduce tsunami run-up heights at Santorini. However, the simulations could not reproduce the initial sea retreat at Kamari and Perissa, and they underestimated the 20 m run-up height at los by a factor of four. Therefore, it appears unlikely that an underwater explosion was the only source of the 1650 tsunami.

To evaluate the tsunami potential of flank deformation, we extrapolated the present-day southeastern flank slopes toward the sea surface to reconstruct a cone that just breaches the sea surface (Fig. 5) to be consistent with eyewitness accounts¹⁷,¹⁸. This pre-explosion topography (Fig. 5a) and the seismically-mapped detachment surface (Fig. 5b) were used as input surfaces for numerical tsunami simulations using VolcFlow²⁶. Kolumbo's cone consists of poorly consolidated pumice with a median dry density of 0.725 g/cm³ and median vesicularity 67.5 (ref. 18), which, depending on the inter-pore porosity, results in a bulk deposit density of 1250–1500 kg/m³. This density is significantly lower than values used in previous volcanic landslide tsunami simulations, e.g., 1900 kg/m³ at Anak Krakatau⁵, 2000 kg/m³ at Ritter Island²⁷ and 2600 kg/m³ at Tenerife²⁸. In VolcFlow simulations, the slide material density, in combination with its yield strength, controls the dynamics of flank deformation. We tested a large parameter space (Supplementary Information) with yield strengths between 2000 and 75,000 Pa and densities between 1250 and
<figure><img src="image_6.png" /><figcaption>**Fig. 3** | Reflection seismic profiles of the Kolumbo edifice and flanks. a–c Semi-transparent orange regions represent material deposited by the 1650 eruption. Folded and disrupted seismic reflections within the northwestern flank of the volcano are the result of internal deformation, while parallel reflections within the</figcaption></figure>

southeastern flank of the volcano indicate the absence of deformation. A tentative interpretation of a thrust beneath compressional folding is shown by the broken black line in (a), soling into the basal detachment surface (broken pink line). The locations of the profiles are shown in Fig. 1c.

2000 kg/m³. We compared the resulting flank deformation and tsunami heights with the geophysical observations and the eyewitness reports.

Our analyses reveal that a range of yield strength-density pairs can reproduce the observed tsunami run-up heights at Ios, Sikinos, and northern Santorini (Fig. 6a, Supplementary Figs. S1–S25). Given the deposit density range defined by the pumice sample analysis (1250–1500 kg/m³), it is possible to reproduce the tsunami observations with simulations using a yield strength between 5000 and 20,000 Pa, which is the same value range used in previous VolcFlow simulations of Anak Krakatau²⁹ and Ritter Island²⁷. At Perissa and Kamari, the important observation of sea retreat is reproduced, but run-up heights are underestimated (blue dashed line in Fig. 6d–h).

In addition, we simulated an explosion-induced tsunami, which requires an estimate of the peak wave amplitude, previously assumed to be in the range of 150–240 m, to be able to reconstruct tsunami observations at Santorini²⁰. However, considering the crater rim lies in an average water depth of only 135 m, the initial wave height by the explosion is unlikely to have exceeded this. Our simulations of an explosion-induced tsunami with a peak wave height of 150 m can neither reproduce the wave height observed at southern Ios and Sikinos nor the reported initial sea retreat at Kamari and Perissa (green dotted line in Fig. 6d–h).

Since the seismic data show that both source mechanisms could have contributed to tsunami genesis and neither the slope failure nor the explosion-induced tsunami simulations alone could reproduce eyewitness accounts, we performed additional model runs to simulate the effects of a combination of both source mechanisms, superimposing the wave field of an explosion onto that caused by flank deformation. The lag time between flank deformation and the explosion is unknown. However, we interpret that the explosion was caused

by depressurization of a critical, erupting system or exposure of the feeder system, which would require lateral movement of the overlying sediment by a distance on the order of 500–1000 m. Given a simulated velocity of 3 to 4 m/s for the flank movement, this equates to a lag time of -2 to 5 min from the initiation of the slope instability. Our results show that assuming a lag time of 4 min from the beginning of flank deformation to the explosion, and an explosion-derived initial wave height of 150 m, we are able to simulate a tsunami pattern that is consistent with all known historical eyewitness accounts (red solid line in Fig. 6). In addition, various combinations of lag time and explosion-derived initial wave height show good agreement with the run-up heights at Ios as well as the initial sea retreat and subsequent run-up heights on the eastern Santorini coastline (Supplementary Figs. S1–S25). The observed inundation of 240 m on Sikinos cannot be directly transferred to a run-up height due to the simplified topography used in the simulations and the lack of a precise location of the eyewitness account. However, the reported inundation indicates that the southeastern coast of Sikinos was affected by a tsunami with significant run-up height. The good agreement between our simulations and the eyewitness accounts gives confidence in the combined simulation approach and highlights the complexity of tsunami generation related to the cascading chain of events during volcanic eruptions (Fig. 7). While a lateral explosion might be able to lead to a regional tsunami pattern consistent with historical eyewitness accounts, we do not have any evidence for such an explosion. We argue that the simplest explanation for tsunami genesis, supported by our data, is flank movement followed by explosive eruptions.

## Chain of events during the 1650 eruption

The combination of high-resolution marine geophysical datasets and historical eyewitness accounts offers a unique opportunity to
<figure><img src="image_7.png" /><figcaption>**Fig. 4** | Kolumbo's deformed edifice. a 3D view of Kolumbo's northeastern crater wall with exposed dykes and the exposed detachment surface (slide plane). b, c ROV images of pumice layering representative of the 1650 cone's lithology. d ROV image of an exposed dyke within the crater.</figcaption></figure>

reconstruct the sequence of volcanic processes that triggered the 1650 explosion and tsunami. The 1650 cone consists of highly vesicular pumice, which was deposited as fallout from the eruption column, where many of the large pumice clasts floated at the sea surface before becoming water-saturated and sinking¹⁸ (Fig. 7a, b). The 1650 cone had a volume of approximately 4 to 6 km³ (ref. 16), which was deposited in a very short timeframe of only 2 weeks, based on the onset of visually detectable eruption products in eyewitness accounts¹⁸,²⁰,²¹. This rapid deposition has likely precluded significant consolidation of the pumice deposits covering the pre-1650 topography. The pre-1650 topography was characterized by the remnants of a cone on the northwestern side of Kolumbo from the previous eruptive cycle (K3 in Fig. 3), as well as a topographic step cutting through the volcanic edifice. Consequently, deposits from the 1650 eruption were deposited on steep slopes on the northwestern side of Kolumbo, while the southeastern side was relatively flat (Figs. 2 and 7b). ROV-dives targeting neighboring cones from the Kolumbo chain have revealed pronounced biogenic cementation of volcanic material¹⁴,³⁰, thereby forming a stark material contrast with the overlying 1650 eruption products, which is indicated by a positive polarity reflection (Figs. 2 and 3).

Slope stability is governed by sub-seafloor effective stress (total stress minus pore pressure) and material strength properties. The 1650 pumice deposits were sensitive to slope failure triggered by dynamic loading because (1) rapid deposition precluded the development of cohesion and (2) the exceptionally low deposit density (-1250–1500 kg/m³) resulted in low effective stress within the cone. We suggest that the pronounced material contrast between the cone and the older, underlying volcanic material, as well as the steeper north-western slope (19° compared to ~6° elsewhere), explains why the failure occurred on the northwestern flank. Eruptive activity was accompanied by earthquakes noted in the eyewitness reports¹⁸, which

<figure><img src="image_8.png" /><figcaption>**Fig. 5** | Cone and slide plane reconstructions used as simulation input. a 3D view of the reconstructed 1650 eruption cone of Kolumbo. b 3D view showing the slide plane used for the numerical tsunami simulations.</figcaption></figure>

will have caused ground acceleration and a likely associated dynamic pore pressure increase. This scenario is consistent with the development of a listric failure surface that cuts through the northwestern flank of the cone and soles out onto the pre-1650 material, defining the detachment surface of the slope instability (Figs. 2, 3 and 6c).

The water displacement caused by the movement of the volcanic flank not only triggered a tsunami, as demonstrated by our numerical simulations, but also directly affected the dynamics of the ongoing eruption. Considering that the active vent had emerged above sea level at this stage¹⁸, failure of the northwestern flank probably caused the subaerial vent area to slide into the sea along with the flank. Based on the geometry of the detachment surface, the slope failure removed up to 200 m of material, thereby unroofing the underlying magmatic system. Assuming a bulk density of 1300–1500 kg/m³, this unroofing may have resulted in a pressure reduction of up to 2.6 to 3 MPa, affecting the underlying feeder system and magma reservoir at a depth as shallow as 2 km, as indicated by seismic full-waveform inversion results³¹. The failure will have exposed deeper levels in the vent directly to seawater, while crack formation and possible re-activation of existing failure planes will have allowed phreatomagmatic interactions within the edifice (Fig. 7d). Thus, both decompression through unroofing and interaction with seawater must have affected the upper parts of the magma feeding system within the timeframe of 4 min, as suggested by our model results. Decompression of the magma by unroofing will have led to an expansion of the magmatic fluid phase within the already critical system, directly followed by enhanced
<figure><img src="image_9.png" /><figcaption>b Simulation of submarine explosion (150 m peak wave height)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>c Simulation of slope failure and submarine explosion after 240 s</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 6** | Tsunami simulations of the 1650 flank deformation. a Map showing the maximum tsunami heights of the simulation assuming a slide density of 1500 kg/m³ and a yield strength of 7500 Pa. b Map showing the maximum tsunami heights of the simulation assuming an explosion with a peak tsunami height of 150 m. c Map showing the maximum tsunami heights of the simulation assuming a slide density</figcaption></figure>

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

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

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

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

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

of 1500 kg/m³ and a yield strength of 7500 Pa and an additional explosion at 4 min after the initiation of the slide with a peak tsunami height of 150 m. d-h Virtual tide gauges showing the tsunami waveforms up to 30 min after the beginning of the flank deformation.

internal (closed-system) degassing and associated further pressure buildup.

The flank failure took place during an ongoing eruption, meaning that there was a hydraulic connection in the form of a highly vesicular magma¹⁸ between the summit region and the deeper magma feeding system. A decompression of 3 MPa in the conduit at, e.g., 10 MPa ambient pressure (or about 600 m depth within the volcano base, assuming a magma density of 1500 kg/m³) would lead to a sudden volume increase of about 50% for bubbles existing in the magma, and thus produce efficient fragmentation in the upper vent region³². The pressure released upward in the open conduit would be sufficient to increase the internal bubble pressure above the typical yield strength of about 5 MPa for a felsic magma³³, leading to magma fragmentation at shallow conduit levels. The magma strain rate would suddenly exceed a critical limit³⁴, and efficient magma fragmentation would

commence at shallow levels. The increased fragmentation rates lead to a rapid downward migration of the fragmentation level, thereby initiating a chain reaction of accelerated decompression, producing sustained eruption rates, which further decompresses the deeper parts of the magma feeding system leading to a runaway process of magma ascent and subsequent depressurization. Decompression experiments³⁵ suggest a lag time of ~180 s between initial rapid decompression and accelerated degassing caused by bubble coalescence, which would be in agreement with the assumed slope failure dynamics. Magma-seawater interactions leading to phreatomagmatic explosions³⁶ must have played a central role at this stage, adding to the pressure buildup. The combined result of bubble expansion and phreatomagmatic explosions is rapid degassing in the magma feeding system, triggering a chain of highly explosive events and sustaining a runaway explosive eruption. Potential evidence of enhanced
**Fig. 7** | **Schematic of the Kolumbo volcanogenic tsunami cascade.** a Eruption initiates, producing large amounts of pumice that float on the sea surface before settling to the seafloor. b Ongoing eruption and pumice settling results in the establishment of the 1650 cone, which breaches the sea surface. c Northwestern flank deforms and moves downslope, causing a directed tsunami and resulting in depressurization of the underlying volcanic system. d Depressurization leads to degassing and expansion, which triggers a violent explosion. e High-resolution bathymetry and 3D seismic data reveal the present-day architecture of Kolumbo.

multiple caldera-forming eruptions at Tenerife³⁸, and initial results suggest that this process may have triggered the 2022 HTHH eruption⁴. The up to 2500-m-wide crater of Kolumbo is the result of one or multiple violent explosions, which agrees well with historical reports that the explosion was heard at least 400 km away at the Dardanelles¹⁸. The analysis of marine tephra in sediment cores²² indicates that the 1650 eruption produced ~4.4 km³ (dense-rock equivalent (DRE) volume of ~2.5 km³) of tephra. Considering the large amounts of pumice forming the 1650 cone (~4-6 km³, ref. 16) and the reported pumice rafts floating throughout the southern Aegean Sea¹⁸, the total eruption volume was larger than 10 km³ (5 km³ DRE). In comparison, the 1980 Mount St. Helens eruption produced only 1.1 km³ of tephra, while first mass flow rate-based estimates suggest an eruption volume of ~1.9 km³ for the 2022 HTHH eruption³⁹.

## Toward improved flank stability monitoring and volcanogenic tsunami early-warning systems

The 3D seismic analysis of Kolumbo highlights the complexity of tsunami genesis during volcanic eruptions, with flank deformation being the nucleus of a chain of events in a critical system. Internal deformation within volcanic flanks can progress over long periods, as seen by 3D seismic data of Ritter Island's failed flank⁸ and pre-collapse InSAR data for Anak Krakatau.⁷ The most prominent historical sector collapse triggered the 1980 eruption of Mount St. Helens, which was preceded by intrusion-controlled flank deformation. Although not all deformation of a volcanic flank will lead to a catastrophic collapse, deformation was a key precursor to these historical volcanic collapse events. Current monitoring strategies targeting the detection of flank deformation are based on satellite data and repeated photogrammetric imaging,⁷ methods that cannot be used to survey submerged volcanic edifices. Submarine geodetic networks have proven capable of detecting periodic slip events at Mount Etna¹¹, indicating that submarine deformation monitoring is possible. However, neither remote sensing-based subaerial approaches nor submarine geodetic monitoring techniques allow real-time observations of flank deformation. This may not be a problem for detecting prolonged deformation as a precursory warning signal but would be unsatisfactory in the case of rapidly developing systems like Kolumbo or for detecting and characterizing the transition from slow to fast (tsunamigenic) flank movement. It is not yet known whether unroofing by a slope failure played a key role in triggering the 2022 HTHH eruption, nor whether the edifice underwent longer-term (weeks to months) deformation prior to the eruption.

There are several active (sub) marine volcanoes that have been affected by repeated slope failures in the past but which are currently not monitored. Kick'em Jenny (Lesser Antilles) is highly active and has been affected by repeated small-scale cone collapses in recent decades⁴⁰. While these collapses have not resulted in significant tsunami events in recent history, the shallow water depth of the active vent (~190 m), and its formation on the scar of a major sector collapse, highlight the potential hazard. Slope stability analyses and numerical simulations indicate a significant tsunami hazard associated with a potential sector collapse⁴¹, which would be amplified further by a phase of explosive eruptive activity, as indicated by the 1650 Kolumbo and 2022 HTHH eruptions. Another candidate for future submarine

phreatomagmatic fragmentation of Kolumbo magma may lie in the abundance of very fine-grained tephra found in close proximity to the volcano in box cores and chemically correlated to the products of the 1650 event³⁷.

Depressurization by unroofing is known to have triggered violent explosive eruptions such as the 1980 Mount St. Helens eruption and

<figure><img src="image_17.png" /></figure>
flank monitoring is the Myojinsho submarine volcano in Japan⁴². Its 1952–1953 explosive eruptions bore many similarities to the evolution of Kolumbo by forming a ~4 km diameter symmetrical volcaniclastic cone, which breached the sea surface before being destroyed by violent explosions⁴². While there is no evidence for past slope deformation in the available low-resolution bathymetry data, the steeply sloped cone (21°) formed on the northeastern rim of a larger pre-existing caldera. Renewed activity could result in a flank collapse of the volcaniclastic cone and depressurization of the underlying magmatic system, similar to the Kolumbo eruption, with the potential for the generation of much larger tsunamis than those caused by the 1952–1953 cone building activity. Finally, Kolumbo itself poses a significant hazard to the densely populated northern and western coasts of Santorini. A low P-wave velocity zone at 2 to 4 km depth indicates a shallow melt reservoir beneath Kolumbo, and repeated earthquake swarms may indicate fracturing associated with ascending melts³¹,⁴³. Renewed volcanic activity or a strong earthquake may reactivate the deep-seated deformation of the northwestern flank, with the potential to trigger tsunami waves that would arrive at Santorini within 5 minutes (Fig. 6). Bottom waters in Kolumbo's crater are highly acidic⁴⁴, and even a minor slope failure could potentially trigger a limnic eruption releasing large amounts of CO₂ and other toxic trace gases like SO₂. Given that most fatalities from the 1650 eruption of Kolumbo were associated with toxic gases¹⁸, this scenario deserves particular attention. However, currently available submarine monitoring systems require long lead-in times (often several years) and costly infrastructure, which would have precluded monitoring efforts even if subaerial precursors had been identified. Local populations, decision-makers and scientists are currently unprepared for the threats posed by submarine eruptions and slope failures, as has been demonstrated by the recent 2018 sector collapse of Anak Krakatau and the 2022 HHTH eruption. Therefore, new shore-line crossing monitoring strategies (as tests as part of the SANTORY project⁴⁵) are required that are capable of being deployed as part of rapid response initiatives during volcanic unrest and which enable real-time observation of slope movement.

## Methods

### 3D seismic survey

The high-resolution P-Cable 3D seismic dataset was collected during research cruise POS538 aboard R/V Poseidon in October 2019 (ref. 46). The receiver array consisted of 16 streamers (each with 8 channels, spaced 1.5625 m apart) attached to a cross cable of ~200 m length. The source was a 45in³–45in³ generator-injector airgun (harmonic mode), fired every 4 s, with a nominal ship speed of 3.5 knots through water. Processing included geometry definition, trace editing (removing dead and very noisy traces), static time corrections, anomalous amplitude attenuation, frequency-wavenumber (FK) filtering and trace interpolation. We then applied adaptive subtraction multiple suppression on common receiver gathers using a 1D pattern-matching algorithm. The resulting data were sorted to common midpoints on a grid with cell sizes of 3.125 m by 3.125 m, corrected for normal moveout using a water velocity of 1517 m/s, and then stacked. Gaps were filled via post-stack trace interpolation in both the inline and crossline directions. Signal to noise ratio was further improved with a post-stack 3D F-K coherency filter. We migrated the data with a constant velocity (1517 m/s) 3D Stolt migration, followed by residual 2D finite difference migrations applied in the inline and then crossline direction. We estimated velocities for the finite difference migrations via semblance analysis of long-offset 2D seismic data that crossed our 3D survey. The final data volume has a maximum spatial resolution of 3.125 m, and a maximum vertical resolution (half wavelength) near the seafloor of ~7 m, based on a dominant frequency of ~100 Hz.

## Numerical tsunami simulations

We used the numerical simulations code VolcFlow for all presented simulations (details in refs. 26–29,47). VolcFlow simulates the slide as a fluid and couples the interaction between the slide and the water based on a depth-averaged approximation on a topography-linked coordinate system. VolcFlow applies a plastic behavior with turbulence as rheology, which is primarily controlled by yield strength and density. The code uses general depth-averaged equations of mass and momentum conservation¹⁴,⁴⁸. The slide was defined by reconstructing a pre-slide topography (Fig. 5), which was achieved by using the geometry of the remaining cone and constraints about the pre-slide geometry (i.e., that the cone breached the sea surface) and mapping the detachment surface in the 3D seismic datasets as slide plane, which had to be converted into the time domain using a seismic velocity of 1.7 km/s (ref. 15). The simulations had a lateral resolution of 100 m and simulated the wave field for 30 min. We performed 20 slide simulations using various densities (1250, 1500, 1750 and 2000 kg/m³) and yield strengths (5000, 7500, 10,000, 20,000 and 50,000 Pa) to test the parameter space (see. Supplementary Figs. S1–S25). To simulate the effects of an explosion, we followed the approach by Ulvrova et al.²⁰ and assumed a radial-symmetric initial explosion-induced water undulation for four different peak wave heights (100, 125, 150 and 240 m). Finally, we performed 20 simulations combining a slide and explosion using the same explosion-induced water undulations and five different time gaps between slide and explosion (1, 2, 3, 4, 5 min) combined with a slide with a representative density of 1500 kg/m³ and yield strength of 7500 Pa (as described in ref. 26). For this, the initial wave field from the explosion-derived tsunami simulations was added to the wave field of the landslide tsunami simulations at one specific time step (1, 2, 3, 4, 5 min) and the code used this superimposed wave fields for the rest of the simulations. To compare simulated wave heights with eyewitness accounts, we created maps showing the maximum wave height as well as virtual tide gauges at locations where direct tsunami observations are available. The virtual tide gauges were placed in shallow water (5–10 m and not at the coastline) to allow analyzing a potential sea retreat at each location.

## Autonomous underwater vehicle (AUV) mapping

During cruise POS510 (ref. 44), detailed bathymetric data were acquired through 7 planned dives by AUV Abyss, which is equipped with the dual frequency (200/400 kHz) RESON Seabat 7125 multibeam. During the acquisition, the nominal frequency of 200 kHz was applied, and all the logged data were further processed using the open-source software package MB-System. First, the navigation was adjusted using the interactive graphical interface program MBnavadjust and then merged with the bathymetric data. The following step included identifying and cleaning erroneous beams; therefore, dedicated filters were applied and manual cleaning through 2D and 3D editors. Finally, the processed bathymetric data were gridded, yielding a spatial resolution of 2 m, using the function mbgrid and, more specifically, the Gaussian weighted mean algorithm.

## Data availability

The AUV bathymetry dataset is available at the PANGAEA data repository [https://doi.org/10.1594/PANGAEA.958275]. The entire 3D seismic dataset will become available in September 2025 at the PANGAEA data repository [https://doi.org/10.1594/PANGAEA.960973]. For previous access on reasonable request, please contact the authors.

## Code availability

The used numerical tsunami simulation code VolcFlow can be downloaded at https://lmv.uca.fr/volcflow/.
## References

1. Auker, M. R., Sparks, R. S. J., Siebert, L., Crosweller, H. S. & Ewert, J. A statistical analysis of the global historical volcanic fatalities record. *J. Appl. Volcanol.* **2**, 1–24 (2013).

2. Paris, R. et al. Volcanic tsunami: a review of source mechanisms, past events and hazards in Southeast Asia (Indonesia, Philippines, Papua New Guinea). *Nat. Hazards* **70**, 447–470 (2014).

3. Seabrook, S. et al. Pyroclastic density currents explain far-reaching and diverse seafloor impacts of the 2022 Hunga Tonga Hunga Ha'apai eruption. Preprint at ResearchSquare https://doi.org/10.21203/rs.3.rs-2395332/v1 (2023).

4. Cronin, S. et al. *The 15 January 2022 Hunga Eruption, Tonga—First Petrographic and Geochemical Results (No. EGU22-13584)* (Coper-nicus Meetings, 2022).

5. Grilli, S. T. et al. Modelling of the tsunami from the December 22, 2018 lateral collapse of Anak Krakatau volcano in the Sunda Straits, Indonesia. *Sci. Rep.* **9**, 1–13 (2019).

6. Chaussard, E. & Amelung, F. Precursory inflation of shallow magma reservoirs at west Sunda volcanoes detected by InSAR. *Geophys. Res. Lett.* **39**, 21311 (2012).

7. Walter, T. R. et al. Complex hazard cascade culminating in the Anak Krakatau sector collapse. *Nat. Commun.* **10**, 1–11 (2019).

8. Karstens, J. et al. From gradual spreading to catastrophic collapse—reconstruction of the 1888 Ritter Island volcanic sector collapse from high-resolution 3D seismic data. *Earth Planet. Sci. Lett.* **517**, 1–13 (2019).

9. Watt, S. F. et al. From catastrophic collapse to multi-phase deposition: flow transformation, seafloor interaction and triggered eruption following a volcanic-island landslide. *Earth Planet. Sci. Lett.* **517**, 135–147 (2019).

10. Puglisi, G. & Bonforte, A. Dynamics of Mount Etna Volcano inferred from static and kinematic GPS measurements. *J. Geophys. Res. Solid Earth* **109**, B11404 (2004).

11. Urlaub, M. et al. Gravitational collapse of Mount Etna's southeastern flank. *Sci. Adv.* **4**, eaat9700 (2018).

12. Segall, P., Desmarais, E. K., Shelly, D., Miklius, A. & Cervelli, P. Earthquakes triggered by silent slip events on Kīlauea volcano, Hawaii. *Nature* **442**, 71–74 (2006).

13. Patrick, M. R. et al. The cascading origin of the 2018 Kīlauea eruption and implications for future forecasting. *Nat. Commun.* **11**, 5646 (2020).

14. Nomikou, P. et al. Submarine volcanoes of the Kolumbo volcanic zone NE of Santorini Caldera, Greece. *Glob. Planet. Change* **90**, 135–151 (2012).

15. Preine, J. et al. Spatio-temporal evolution of the Christiana-Santorini-Kolumbo volcanic field, Aegean Sea. *Geology* **50**, 96–100 (2022).

16. Hübscher, C., Ruhnau, M. & Nomikou, P. Volcano-tectonic evolution of the polygenetic Kolumbo submarine volcano/Santorini (Aegean Sea). *J. Volcanol. Geotherm. Res.* **291**, 101–111 (2015).

17. Fouqué, F. *Santorin et Ses Eruptions* (G. Masson, 1879).

18. Cantner, K., Carey, S. & Nomikou, P. Integrated volcanologic and petrologic analysis of the 1650 AD eruption of Kolumbo submarine volcano, Greece. *J. Volcanol. Geotherm. Res.* **269**, 28–43 (2014).

19. Klaver, M. et al. A distinct source and differentiation history for Kolumbo submarine volcano, Santorini volcanic field, Aegean arc. *Geochem. Geophys. Geosyst.* **17**, 3254–3273 (2016).

20. Ulvrova, M. et al. Source of the tsunami generated by the 1650 AD eruption of Kolumbo submarine volcano (Aegean Sea, Greece). *J. Volcanol. Geotherm. Res.* **321**, 125–139 (2016).

21. Dominey-Howes, D. T. M., Papadopoulos, G. A. & Dawson, A. G. Geological and historical investigation of the 1650 Mt. Columbo (Thera Island) eruption and tsunami, Aegean Sea, Greece. *Nat. Hazards* **21**, 83–96 (2000).

22. Kutterolf, S. et al. The medial offshore record of explosive volcanism along the central to eastern Aegean Volcanic Arc: 2. Tephra ages and volumes, eruption magnitudes and marine sedimentation rate variations. *Geochem. Geophys. Geosyst.* **22**, e2021GC010011 (2021).

23. Pope, E. L. et al. Origin of spectacular fields of submarine sediment waves around volcanic islands. *Earth Planet. Sci. Lett.* **493**, 12–24 (2018).

24. Miramontes, E., Garziglia, S., Sultan, N., Jouet, G. & Cattaneo, A. Morphological control of slope instability in contourites: a geo-technical approach. *Landslides* **15**, 1085–1095 (2018).

25. Kelfoun, K., Giachetti, T. & Labazuy, P. Landslide-generated tsunamis at Réunion Island. *J. Geophys. Res. Earth Surf.* **115**, F04012 (2010).

26. Karstens, J., Kelfoun, K., Watt, S. F. & Berndt Combining 3D seismics, eyewitness accounts and numerical simulations to reconstruct the 1888 Ritter Island sector collapse and tsunami. *Int. J. Earth Sci.* **109**, 2659–2677 (2020).

27. Giachetti, T., Paris, R., Kelfoun, K. & Pérez-Torrado, F. J. Numerical modelling of the tsunami triggered by the Güimar debris avalanche, Tenerife (Canary Islands): comparison with field-based data. *Mar. Geol.* **284**, 189–202 (2011).

28. Giachetti, T., Paris, R., Kelfoun, K. & Ontowirjo, B. Tsunami hazard related to a flank collapse of Anak Krakatau Volcano, Sunda Strait, Indonesia. *Geol. Soc. Lond. Spec. Publ.* **361**, 79–90 (2012).

29. Tucker, M. E. et al. Carbonate crusts around volcanic islands: composition, origin and their significance in slope stability. *Mar. Geol.* **429**, 106320 (2020).

30. Chrapkiewicz, K. et al. Magma Chamber Detected Beneath an Arc Volcano With Full-Waveform Inversion of Active-Source Seismic Data. *Geochem. Geophys. Geosystems* **23**, e2022GC010475 (2022).

31. Cashman, K. & Blundy, J. Degassing and crystallization of ascending andesite and dacite. *Philos. Trans. R. Soc. Lond. A Math. Phys. Eng. Sci.* **358**, 1487–1513 (2000).

32. Papale, P. Strain-induced magma fragmentation in explosive eruptions. *Nature* **397**, 425–428 (1999).

33. Gonnermann, H. M. Magma fragmentation. *Annu. Rev. Earth Planet. Sci.* **43**, 431–458 (2015).

34. Burgisser, A., Bergantz, G. W. & Breidenthal, R. E. Addressing complexity in laboratory experiments: the scaling of dilute multi-phase flows in magmatic systems. *J. Volcanol. Geotherm. Res.* **141**, 245–265 (2005).

35. Cas, R. A. & Simmons, J. M. Why deep-water eruptions are so different from subaerial eruptions. *Front. Earth Sci.* **6**, 198 (2018).

36. Fuller, S., Carey, S. & Nomikou, P. Distribution of fine-grained tephra from the 1650 CE submarine eruption of Kolumbo volcano, Greece. *J. Volcanol. Geotherm. Res.* **352**, 10–25 (2018).

37. Hunt, J. E., Cassidy, M. & Talling, P. J. Multi-stage volcanic island flank collapses with coeval explosive caldera-forming eruptions. *Sci. Rep.* **8**, 1–11 (2018).

38. Yuen, D. A. et al. Under the surface: pressure-induced planetary-scale waves, volcanic lightning, and gaseous clouds caused by the submarine eruption of Hunga Tonga-Hunga Ha'apai volcano. *Earthq. Res. Adv.* **2**, 100134 (2022).

39. Allen, R. W. et al. 30 years in the life of an active submarine volcano: a time-lapse bathymetry study of the Kick-'em-Jenny Volcano, Lesser Antilles. *Geochem. Geophys. Geosyst.* **19**, 715–731 (2018).

40. Dondin, F., Heap, M. J., Robertson, R. E. A., Dorville, J. F. & Carey, S. Flank instability assessment at Kick-'em-Jenny submarine volcano (Grenada, Lesser Antilles): a multidisciplinary approach using experiments and modeling. *Bull. Volcanol.* **79**, 1–15 (2017).

41. Fiske, R. S., Cashman, K. V., Shibata, A. & Watanabe, K. Tephra dispersal from Myojinsho, Japan, during its shallow submarine eruption of 1952–1953. *Bull. Volcanol.* **59**, 262–275 (1998).

42. Schmid, F. et al. Heralds of future volcanism: swarms of micro-seismicity beneath the submarine Kolumbo volcano indicate opening of near-vertical fractures exploited by ascending melts. *Geochem. Geophys. Geosyst.* **23**, e2022GC010420 (2022).
43. Carey, S. et al. CO₂ degassing from hydrothermal vents at Kolumbo submarine volcano, Greece, and the accumulation of acidic crater water. *Geology* **41**, 1035–1038 (2013).

44. Hannington, M. D. RV POSEIDON Fahrtbericht/Cruise Report POS510-ANYDROS: rifting and hydrothermal activity in the Cyclades Back-Arc Basin, Catania (Italy)—Heraklion (Greece) 06.03.–29.03. 2017 (GEOMAR Helmholtz Centre for Ocean Research, 2018).

45. Nomikou, P. et al. SANTORY: SANTORini's Seafloor Volcanic ObservatorY. *Front. Mar. Sci.* **9**, 421 (2022).

46. Karstens, J. et al. R/V Poseidon Cruise Report 538—*THESEUS Tsunami hazard of explosive submarine eruptions, 7th October–28th October, Cartagena (Spain)*—Heraklion (Greece) (GEOMAR Helmholtz Centre for Ocean Research, 2019).

47. Kelfoun, K. & Druitt, T. H. Numerical modeling of the emplacement of Socompa rock avalanche, Chile. *J. Geophys. Res. Solid Earth*, 110 (2005).

48. Kutterolf, S. et al. The medial offshore record of explosive volcanism along the central to eastern Aegean Volcanic Arc: 1. Tephra ages and volumes, eruption magnitudes and marine sedimentation rate variations. *Geochem. Geophys. Geosyst.* **22**, e2021GC010010 (2021).

## Acknowledgements

We thank the masters and crews of RV Poseidon for their support during research cruise POS538. We thank Mark Hannington and Sven Petersen for providing access to the AUV bathymetry dataset. Thanks also to Dirk Klaeschen for giving valuable advice at several key stages of seismic processing. We also thank Schlumberger and IHS for granting educational licenses to their software.

## Author contributions

J.K. and C.B. were responsible for organizing and designing the reflection seismic experiments. J.K., G.J.C., J.P., F.S., P.N., M.K., G.D.V., and J.E. acquired, processed and interpreted the reflection seismic data. S.C. provided the ROV footage from within the crater. J.K. and K.K. performed the numerical tsunami simulations. J.K., G.J.C., T.H., and S.C. integrated the geophysical results into the volcanological framework of Santorini. J.K., G.J.C. and T.H. drafted the manuscript, while J.P., S.C., P.N., M.K., G.D.V., F.S., K.K., J.E., and C.B. discussed the dataset and provided comments and corrections to the manuscript.

## Funding

Open Access funding enabled and organized by Projekt DEAL.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-023-42261-y.

**Correspondence** and requests for materials should be addressed to Jens Karstens.

**Peer review information** *Nature Communications* thanks Thomas Walter, Alberto Armigliato, and the other anonymous reviewer(s) for their contribution to the peer review of this work. A peer review file is available.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**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 license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license 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 license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023