The 16 September 2023 Greenland Megatsunami: Analysis and Modeling of the Source and a Week‐Long, Monochromatic Seismic Signal

Large tsunamigenic landslides have occurred in Greenland in prehistoric and recent times (Korsgaard et al., 2023). The effects of global warming and changes in permafrost are likely to further reduce slope stability and increase the incidence of landslides and tsunamis (Gruber and Haeberli, 2007; Patton et al., 2019; Poli, 2024; Svennevig et al., 2023). The occurrence and propagation of tsunamis, especially in fjord locations, are considered to be one of the most devastating of all natural disasters (Blikra et al., 2006; Evans et al., 2007; Higman et al., 2018; Svennevig et al., 2023). A landslide‐triggered tsunami happens when substantial soil or rock formations abruptly give away and fall or glide into a body of water. Recent and devastating tsunami examples have been taking place and studied mainly in western Greenland, such as in 2017, when a rock avalanche of $50millionm3$ (Paris et al., 2019) impacted the Karrat Fjord; the landslide triggered a tsunami that flooded the village of Nuugaatsiaq, destroying 11 houses and killing 4 people (Paris et al., 2019; Svennevig et al., 2020). However, the East Coast of Greenland can also be the source of megatsunamis that exceed 100 m in height in the near field and with the effects recorded across Europe, as this study shows.

A large number of seismic signals from near‐surface processes have been observed in the Nordic and Polar regions. Very long‐period (VLP) seismic signals have been detected at regional and teleseismic distances (Nettles and Ekström, 2010), most of them originating along the Greenland coast at the location of large outlet glaciers. The high temporal correlation between glacial earthquakes and ice loss, with most of the recorded events occurring in late summer, suggests that most of these signals were related to iceberg collapse (Nettles and Ekström, 2010). Recent efforts to detect and locate the source of long‐period seismic signals on a global scale (Poli, 2024) have revealed tens of thousands of events. The related seismic sources have been identified primarily in the Polar regions, including the events associated with glaciers and landslides (Poli, 2024).

VLP, long‐lasting signals have also been observed in volcanic environments (Ohminato et al., 1998; Kumagai and Chouet, 2000). Long duration (∼20 min), monochromatic (∼16 s) VLP signals, produced by the resonance of a deep depleting magma reservoir at Mayotte, Comoros Island, were observed on a global scale (Cesca et al., 2020). Long oscillations were also observed following the Karrat Fjord landslide in 2017 (Paris et al., 2019): in that case, the seismic signals lasted ∼55 min, with the particle motion of horizontal components oriented perpendicular to the channel, which could be attributed to seiching of the channel.

On 16 September 2023, a megatsunami was triggered in a remote area of East Greenland. There were initial reports of a tsunami at Dickson Fjord on various social media channels (see Data and Resources), where the remote event was identified due to water waves hitting a military unit located on Ella Island.

In this study, we show satellite imagery in which we identify a missing cliff patch on the southern shore of Dickson Fjord (72.817° N, 26.965° W), ∼60 km farther west of Ella Island in the narrowing fjord system, suggesting that a rockslide caused the tsunami. We observe two types of seismic signals, first a short‐duration long‐period (LP) one, then a 7‐day‐long, monochromatic VLP signal. We focus on the seismological observations associated with the Greenland event of 16 September 2023, evaluating full waveforms of sparse networks and teleseismic arrays. We invert the source of the LP signal and the VLP to reconstruct the evolution of the rockslide and fjord resonance processes.

An analytical model of the sloshing physics is used to confirm the conceptual resonance model.

Seismic signals of the 16 September 2023 event were observed globally, sparking curiosity about triggering and resonance or sustainment mechanisms. The signals start on 12:35 UTC with a strong LP signal in the range of 0.02–0.06 Hz, which is best recorded at the closest stations, SCO, SUMG, and DAG, located at less than 600 km distance (Fig. 1a,b).

The LP signal is followed by a monochromatic oscillating VLP signal (Fig. 1b) lasting at least up to 1 week. This VLP signal is visible from regional to teleseismic distances, up to at least 5000 km. Its dominant frequency is ∼0.0109 Hz. This frequency is relatively stable over long time periods. A small decrease in the order of 0.1% of the dominant frequency is seen during the first day (Fig. S1, available in the supplemental material to this article).

The VLP signal onset is emergent; its timing roughly corresponds to one of the LP events.

The VLP signal is slowly and smoothly attenuated but lasts for at least 1 week. The amplitude attenuation is best seen from the envelope of seismic beams at different arrays after the stacking of the velocity signals recorded at single stations (Fig. 1c–h).

Particle motions at regional and teleseismic distances confirm the excitation of both P–SV‐ and SH‐wave types but with a peculiar radiation pattern: SH‐type waves have the highest amplitude at stations located roughly toward east and west, whereas P–SV type waves are predominant to the north and south (Fig. S2).

We analyzed data available from the European Space Agency and from the 430+ Dove mini‐satellite swarm (PlanetLabs Inc.). Specifically, we consider Sentinel‐2 data acquired on 13 and 17 September 2023 to provide a before‐and‐after view of the area with a high spatial resolution of up to 10 m, as well as all Dove data acquired between 23 August 2023 and 29 September 2023 at 3 m resolution. We identify and extract the dimensions of the source area and consider a digital elevation model (ArcticDEM) to estimate morphologic and geometric constraints such as fall heights, slopes, and tsunami runups.

A Dove satellite image taken only ∼1 hr after the collapse (Fig. 2b) shows floating debris, ice fragments, and muddy discolored water in the Dickson Fjord. On the southern cliff of the Dickson Fjord, in a side valley oriented north–south (Fig. 3), we find a new rockslide scar enclosing an area of $1.6×105m2$ located on the upper‐eastern cliff (72.808° N, 26.944° W; Fig. 3a) and dark discoloration of the narrow glacier downstream from this point. Combined with the DEM this suggests that the detached rock mass initially collapsed westward from the 30° to 40° steep eastern cliff, with a fall height of ∼300–400 m, before hitting the western side, deflecting northward along the side valley and then becoming a mixed rock–ice avalanche that propagated north for ∼1.6 km where it reached the Dickson Fjord (Fig. 3a,c). The dark discolored coastal regions in post‐tsunami Sentinel‐2 and Dove data suggest peak tsunami runups exceeding 200 m near the entry point and an average of ∼60 m along a 10 km long profile (Fig. S3).

To study the LP seismic source, we perform centroid moment tensor (CMT) and single force (SF) inversions using the displacement records of 24 broadband seismic stations located between 300 and 2000 km from the rockslide region (Fig. 1). Source inversions are computed using a Bayesian bootstrap‐based probabilistic inversion scheme implemented in Grond (Heimann et al., 2018). Green’s functions were computed with QSEIS (Wang, 1999; Heimann et al., 2019). To obtain more accurate results, we combined crustal and oceanic velocity models (Kennett et al., 1995; Laske et al., 2001; Fig. S4). We invert Love and Rayleigh waves in time and frequency domains using the frequency ranges of 0.02–0.06 Hz when using a crustal model and 0.02–0.04 Hz with the oceanic velocity model. For the assessment of uncertainties, the inversions are run in 101 independent bootstrap reweightings: random weightings of the station‐component‐based misfits provide a measure of uncertainties along with mean solutions.

The CMT source is resolved with an $Mw$ 4.9 ± 0.14—a compensated linear vector dipole (CLVD) component of 5% ± 23% and a negative isotropic component of 27% ± 15% (Fig. S5). The inverted location of the CMT source is 11.9 ± 5.6 km south and 4.7 ± 7.1 km west from the rockslide origin. Because of the large non‐double‐couple components and a rockslide as a presumable source mechanism, we perform an SF inversion. Figure 3 shows a comparison between the SF results and remote sensing data. The force source is resolved 10 ± 3.6 km south and 4.4 ± 4.6 km east from the rockslide origin. The resolved horizontal direction is mainly westward oriented. The vertical component of the force (⁠$Fd$⁠) is oriented downward and $|F|=530$GN (GigaNewtons).

In a simple model, the dominant forces generated by a rockslide consist of two elements: (1) the mass acceleration associated with the detachment of mass uphill in the initial phase of the rockfall, and (2) the deceleration associated with the impact of the mass at the bottom of the slope (e.g., Allstadt, 2013). The sliding of the mass with a constant velocity in between is often assumed to generate negligible forces if the trajectory of the rockfall is straight.

Three different teleseismic arrays in southeast, southwest, and northwest directions (Fig. 1c–h) are used to decipher the amplitude decay over 1 week. Velocity records are filtered (band‐pass 0.008–0.015 Hz), down‐sampled to 1 s, the instrument response is removed, and northeast (NE) components are rotated to radial‐transversal components (see Figs. S8–S11). Cross‐correlation (cc)‐based time shifts are used to compute a common apparent phase velocity for each array. We only consider stations with cc >0.9 to remove bad‐quality traces. We stack the waveforms and compute the envelope for each array (Fig. 4).

Figure 4 shows the envelope fits of the stacked records for each array using equations (3), (5), and (7). Because of the low signal‐to‐noise ratios (SNRs) at later times, we fit the curves in a time window of ∼35 hr, starting 10 min after the event origin time and extrapolate the fit. The energy loss of the oscillation is about 1% per cycle in the first hour and about 0.1% after 2.5 days.

The joint analysis of seismological and remote sensing data resolved three phases of the chain of events (see Fig. 2b). Phase 1 corresponds to a rockslide, which produced a short‐duration LP seismic signal visible at regional distances. The approximate location of the rockslide signal from the CMT and SF inversion generally agrees with the very accurate location retrieved from the satellite records. A view comparison before and after the landslide (Fig. 3a) illustrates a difference in the upper‐eastern cliff of a side valley. The rockslide fell a westward‐dipping slope impacting the narrow valley. We estimate a slope of 30°–40° and a fall height of 300–400 m for the initial rockfall originating on the upper‐eastern cliff. This scenario is independently confirmed by modeling LP seismological data using a single force, a model that has been used to reproduce seismic signals from the landslides (e.g., Allstadt, 2013; Chao et al., 2018). The west and downward orientation of the resolved SF is in good agreement with the slope orientation and the landslide path; therefore, we assume that the LP signal is caused by the impact of the collapsing material on the western wall of the side valley.

During phase 2, the slide continued northward along the narrow canyon down to the Dickson Fjord. The narrow canyon hosted a glacier, darkly discolored from the point of collapse, indicating a mixed rock–ice mass movement and possible volumetric gain before reaching the ocean.

Phase 3 corresponds to the generation of a tsunami, triggered by the impact of the rock–ice avalanche into the Dickson Fjord. The tsunami was reported by its impact on Ella Island. Further evidence for a large tsunami is given by remote sensing imaging. The consequent water redistribution built up a seiching process in the Dickson Fjord. The seismological evidence for this stage is given by the long‐lasting VLP signal. Indeed, the source of the VLP signal, which will last more than 1 week, begins at the time of or shortly after the landslide event. This timing suggests a direct connection between the landslide occurrence and the following tsunami. A second argument supporting the role of the tsunami in the generation of the VLP signal is the location and the orientation of the single force obtained by waveform inversion, which is dominantly horizontal and perpendicular to the Dickson Fjord. A plausible source mechanism should explain the peculiar features of the VLP signal, such as its dominant frequency, long duration, and very slow, and quadratic attenuation.

We hypothesize here a possible source model and support it with analytical computations. The initial excitation is the flow of the ice–rock debris into the narrow Dickson Fjord and the resulting tsunami trigger, which might have exceeded 200 m in height in its initial stage.

The landslide plunged into a segment of the Dickson Fjord, which is almost straight on a length of about 20 km. The western end of this segment is a turn to the northwest where the fjord is terminated by a glacier after farther 5 km. The eastern end of the segment is a sharp southward corner. As a first‐order model, we can treat this segment as a box of $L×W=20×2km$ (Fig. 2b). Considering this geometry we estimate that for the fundamental mode and a wave height of 1 m, an oscillating horizontal force with an amplitude of 64 GN is produced. In addition, alternating gravitational loading results in a vertical force pair with an amplitude of 255 GN (see details in the supplemental material). From the SF inversion, we obtain a horizontal force amplitude of 64 GN about 5 hr after the onset of the signal. The initial amplitudes are at about 160 GN, which would translate to an average fundamental‐mode wave height over the fjord segment of about 2.6 m. Assuming a standard global velocity model, the expected maximum vertical ground velocity at 1500 km distance is ∼41 and 3.5 nm/s for the horizontal and vertical oscillator models, respectively, in agreement with observed ground velocities (Fig. S12).

Waveform records of the VLP signal show that the seiching starts with an incoherent phase from which a fundamental oscillation mode develops within the first minutes. The tsunami wave becomes coherent, propagating back and forth across the narrow width of the Dickson Fjord, alternately hitting the southern and northern shores. The build‐up of the fundamental‐mode oscillation can be qualitatively reproduced by finite‐difference modeling of shallow water waves in a rectangular resonator box (aquarium) (Figs. S13–S14 and videos). Beating patterns due to interference of the multiple reflected waves are suggested by the model, especially in the early phase of the oscillation. Such beating patterns can also be seen in the VLP observation (Fig. S9). Waves are efficiently trapped in the short dimension of the resonator, even when one of the ends in the long dimension is modeled as completely open. Longer‐period eigenfrequencies along the long dimension of the fjord decay more quickly because of the dissipation of wave energy toward the ocean. In addition, these modes are slower by a factor of 10, given the aspect ratio of the fjord segment, and would be difficult to observe with the seismometers.

The north–south‐directed motion of the standing wave involves the periodic application of vertical, gravitational, and horizontal forces to account for periodic water mass shifts from north to south and vice versa. The amplitude of the resulting signal may be damped with time by different processes: (1) continuous radiation of water wave energy toward the ocean, (2) friction loss at the fjord seafloor and around obstacles, (3) internal friction in the water body, especially at geometrical irregularities, in which oscillation patterns due to different eigenmodes meet, and (4) excitation of seismic waves. Radiation of tsunami waves toward the ocean (1) may contribute to the exponential part of the decay. Seismic waves (4) must be a minor contributor to the attenuation of the seiche: less than 100 W of seismic energy is radiated in the case at hand (rough estimate: $P=1/2ρv2$ × $A=14W$ with $ρ=3000kg/m3$⁠, $A=2πRD$ with $v=100nm/s$⁠, $R=1500km$⁠, and effective surface‐wave depth D of 100 km). Remarkably, this power is sufficient to be globally observable but it does not significantly contribute to the damping of the seiche. Boundary friction associated with a laminar flow (2) would be confined to a thin boundary of a few centimeters at the seafloor and would be proportional to flow velocity (e.g., Batchelor, 1967). Numerical modeling of a laminar boundary layer with periodically changing flow direction suggests an energy loss on the order of $10−5–10−4$ per cycle for an assumed water depth of 160 m (supplemental material). Internal friction in the water wave (3) is a minor contributor to the energy loss of the seiching: a simple analytical model of the fundamental‐mode oscillation suggests an energy loss of $10−13$ per cycle (supplemental material). Although the exact determination of fluid dynamics acting within the seiching fjord is beyond the scope of our study, we here compare curve fits to the decay of the VLP signal as a seismically derived proxy.

In Figure 4, we showed three possible decay fits to explain the observed amplitude decrease. Equation (5) fits the signal but overestimates the total duration. Equation (3) reproduces the damping only at the beginning and significantly underestimates the duration. A combination of both terms (equation 7) reproduces the decay during the total time period best. The terms included in equations (3) and (5) are physically interpreted as a loss of energy and friction, respectively, meaning that friction plays a major role in producing long‐lasting oscillation amplitude decay.

Water level changes may lead to small changes in the oscillation frequency because wave speed c is related to water depth H by $c=gH$⁠. Such water level changes may be introduced by the outflow of the water, which is initially displaced by the injected landslide volume. During the first hours, the observed peak frequency changes in the order of ∼0.1% could correspond to 20–40 cm of water level change, equivalent to a displaced water volume of $8−10×106m3$ or 50–100 m of rock height on the landslide area. Alternatively or additionally, small frequency changes during the early phase of the oscillation can also be attributed to geometrical irregularities of the fjord. Changes in frequency are observed in the numerical modeling of the oscillation build‐up. When the water depth or the width varies along the major axis of the channel, the laterally shifting interference patterns that occur when the fundamental‐mode oscillation stabilizes lead to small changes in the frequency of the average oscillation (see description and video links in the supplemental material). A minor variation of maximum frequencies between 0.010874 and 0.010879 Hz with a rhythm of about 6 hr is observed after the initial frequency decrease, which seems to be correlated to the tides (Fig. S1c). Interestingly, the observed modulation of the dominant eigenfrequency is opposed to expectations: larger water depths at high tides would suggest higher velocities of gravity waves in the water column and, therefore, higher eigenfrequencies, as opposed to observing lower frequencies at high tides. Other factors influencing the water height and frequency can be diurnal changes in wind directions or meltwater inflow from the surrounding glaciers. In addition, in a complex fjord system, changes in the cross‐sectional area between connected fjord basins due to tidal changes in water depth could explain shifts in eigenfrequencies (e.g., Dahm, 1992). Because of the decreasing SNRs, observations of changes in the oscillation frequency are limited to three tidal cycles.

The radiation pattern of the VLP horizontal components signal illustrates a clear difference between the stations located east–west and north–south from the landslide origin, where transversal components are dominant in the east–west direction, whereas vertical components are dominant in the north–south direction (Fig. S2). This pattern is in agreement with the radiation produced by an SF source. The observed particle motions support that the VLP signal triggered by the landslide could be produced by the seiching effect of water waves in the channel after the impact of material in the north direction.

On 16 September 2023, a mass‐wasting event occurred at 12:35 UTC in East Greenland, where a rockmass on an area of $1.6×105m2$ and a height of 50–100 m detached and impacted initially to the west, falling around 300–400 m in a slope of 30°–40°. The initial west direction is well resolved with an SF model, which retrieves $|F|=530$GN. After the first impact, the rockslide material is redirected to the north and interacts with the glacier until it enters the Dickson Fjord, where it triggers a megatsunami with a peak run‐up height locally exceeding 200 m. After the LP signal produced by the landslide, a VLP monochromatic oscillation was generated and lasted at least one week. This long signal and its amplitude decay are well explained with a damped oscillator model that contains a linear and quadratic contribution, in which the major role is played by the quadratic term, but a combination of both factors is needed to fit the amplitude decay. A fundamental‐mode standing wave with an average initial amplitude of about 2.6 m in the narrow direction of a 2 × 20 km fjord segment explains the dominant frequency of the VLP oscillation and gives a plausible mechanism for the strength, radiation pattern, and long duration of the seismic VLP signal.

The supplemental material includes additional figures, videos, and text that provide more detailed insights into the methodology and the data sets. Software packages Pyrocko (Heimann et al., 2017), NumPy (https://numpy.org/; Harris et al., 2020), SciPy (https://docs.scipy.org/doc/), Generic Mapping Tools (GMT) (Wessel et al., 2013), and Matplotlib (Hunter, 2007) were used for data processing and plotting. Social media reports were available at https://twitter.com/OJoelsen/status/1704839193623400519 and https://twitter.com/K_Svennevig/status/1705161265625157733. All station sensor orientations were checked with AutoStatsQ (Petersen et al., 2019). Shuttle Radar Topography Mission (SRTM) (Farr et al., 2007) and 1 Arc‐Minute Global Relief Model (ETOPO1) (Amante and Eakins, 2009) topographic data were used. Sentinel‐2 and Dove satellite data were analyzed with the open‐source software QGIS version 3.36 (https://www.qgis.org/). Sentinel‐2 data are free to access from the European Space Agency’s (ESA) Copernicus program (https://dataspace.copernicus.eu/) and the images can be available at https://sentinelshare.page.link/sYLK. Dove data are free to use for academics and research and were provided to T. R. W. at the University Potsdam, Image © 2023 Planet Labs PBC. The digital elevation model (ArcticDEM) is freely available at https://livingatlas2.arcgis.com/arcticdemexplorer/. The overview image is taken from Google Earth Pro (Landsat/Copernicus). Seismic waveform data used in this study are publicly available at GEOFON (https://geofon.gfz-potsdam.de/waveform/webservices/fdsnws.php), Orfeus (https://www.orfeus-eu.org/data/eida/webservices/), and Incorporated Research Institutions for Seismology (IRIS) (https://service.iris.edu/fdsnws/) data centers. Citations of used seismic networks are provided in Table 1. All websites were last accessed in April 2024.

The authors acknowledge that there are no conflicts of interest recorded.

A. C. P. received funds from the National Agency for Research and Development (ANID) /Scholarship Program/DOCTORADO BECAS CHILE/2019‐72200544.