Research ArticleGEOPHYSICS

The 22 December 2018 tsunami from flank collapse of Anak Krakatau volcano during eruption

See allHide authors and affiliations

Science Advances  15 Jan 2020:
Vol. 6, no. 3, eaaz1377
DOI: 10.1126/sciadv.aaz1377

Abstract

On 22 December 2018, a devastating tsunami struck Sunda Strait, Indonesia without warning, leaving 437 dead and thousands injured along the western Java and southern Sumatra coastlines. Synthetic aperture radar and broadband seismic observations demonstrate that a small, <~0.2 km3 landslide on the southwestern flank of the actively erupting volcano Anak Krakatau generated the tsunami. The landslide did not produce strong short-period seismic waves; thus, precursory ground shaking did not provide a tsunami warning. The source of long-period ground motions during the landslide can be represented as a 12° upward-dipping single-force directed northeastward, with peak magnitude of ~6.1 × 1011 N and quasi-sinusoidal time duration of ~70 s. Rapid quantification of a landslide source process by long-period seismic wave inversions for moment-tensor and single-force parameterizations using regional seismic data available within ~8 min can provide a basis for future fast tsunami warnings, as is also the case for tsunami earthquakes.

INTRODUCTION

On 22 December 2018, a modest amplitude but devastating tsunami struck Indonesia’s Sunda Strait (Fig. 1) from around 14:27 to 15:00 UTC (21:27 to 22:00 local time), leaving at least 437 dead and tens of thousands injured. The tsunami swept onto the coastlines of southeastern Sumatra and western Java without warning from either earthquake ground shaking or any alert system. The source of the tsunami was along the southwest coast of Anak Krakatau, an active volcano on the rim of the caldera of the great 1883 Krakatau eruption (1). The 1883 Krakatau eruption and its subsequent major tsunamis caused more than 35,000 casualties at coastal cities around the Sunda Strait, but the mechanism generating the tsunami is still controversial, with proposed ideas of pyroclastic flow [e.g., (26)] and submarine explosion (79). Anak Krakatau first emerged from the ocean in 1929 and rapidly built a cone, rising ~300 m above sea level by 2018. The volcano had been in active Strombolian eruption since June 2018, including throughout the day of 22 December [e.g., (10, 11)].

Fig. 1 Affected area, tsunami amplitude, and arrival times at tidal gauge stations from the 22 December 2018 Sunda Strait tsunami.

Cyan zones indicate areas affected by the tsunami, adapted from the map of Tsunami Selat Sunda created by Badan Nasional Penanggulangan Bencana dated 14 January 2019. Zero-to-peak wave heights (H) and arrival times (hh:mm UTC on 22 December 2018) are from Joint Research Centre emergency reporting released by 24 December 2018. The origin time of the source is estimated to be around 13:57 UTC, as discussed in the text. The insert map locates the Sunda Strait and Anak Krakatau in Indonesia.

The potential for a flank collapse of rapidly growing Anak Krakatau was recognized by Giachetti et al. (12), who modeled the possibility of resulting tsunami generation and impact on the adjacent Sumatra and Java coastlines. Located on the margin of the 200- to 300-m below–sea level caldera floor of the 1883 eruption, the lack of buttressing support on the southwest side of Anak Krakatau was recognized as creating an unstable edifice. By modeling an abrupt 0.28-km3 flank collapse directed southwestward, it was anticipated that (i) immediately surrounding islands (Sertung, Panjang, and Rakata) would receive tsunami amplitudes of 15 to 30 m, and (ii) tide gauges at the coastal cities along Java and Sumatra would receive waves from 0.3 to 3.4 m high within 35 to 60 min after the onset of collapse. It was proposed that rapid detection of the collapse, together with an efficient alert system on the coast, could prevent such an event from being deadly. Unfortunately, the actual scenario played out close to the modeled case, with ~25-m tsunami run-up documented on Rakata (Smithsonian Institution Global Volcanism Program: https://volcano.si.edu/volcano.cfm?vn=262000), and tsunami zero-to-peak amplitudes at tide gauge stations on Java and Sumatra being 0.2 to ~1.0 m, with travel times from the source of ~30 to 60 min (Fig. 1). Post-tsunami survey by Tsunami and Disaster Mitigation Research Center (http://tdmrc.unsyiah.ac.id/the-latest-update-from-post-sunda-strait-tsunami-survey/) and another field survey found 10 to 15 m run-ups locally at Cipenyu beach in Pandeglang and up to ~5 m run-ups at coastal areas from western Banten to southern Lampung (13). However, there was no rapid detection of the collapse or any tsunami alert. The impact was particularly severe in western Java, where 317 fatalities occurred as the waves arrived beginning around 21:27 in the evening (14:27 UTC). The observed tsunami can be accounted for by a flank collapse, as shown by Grilli et al. (14).

During the daylight and early evening hours of 22 December before the tsunami, photographers such as O. Andersen (https://www.oysteinlundandersen.com) captured glowing eruptions at the summit of Anak Krakatau, readily visible from the touristic Java and Sumatra coasts. Before the collapse, the island misted over according to various photographs, presumably due to magma reaching the ocean, but there was no visible sign of unusual behavior of the volcano. The infrared window of Himawari-8 images show evidence that two pulses began at ~13:40 (UTC; 20:40 local time) and at ~15:40 (UTC; 22:40 local time) on 22 December and the late infrared pulse lasted until the end of 23 December (https://cimss.ssec.wisc.edu/goes/blog/archives/31176). Low levels of ground shaking throughout the day detected by nearby seismic instruments appear to include minor motions from the Strombolian eruptions. When the tsunami struck, it was not immediately apparent how it had been generated. No associated seismic event was reported by the U.S. Geological Survey National Earthquake Information Center (USGS-NEIC). In the daylight of 23 December, it became clear that there was ongoing, greatly increased phreatomagmatic eruption with dense ash, strongly suggesting that some process associated with Anak Krakatau had produced the tsunami. We quantify the source of the tsunami using geodetic and seismological observations and consider strategies for seismic monitoring of such sources to mitigate future impacts.

RESULTS

Indonesia has a high-quality broadband seismic network with stations distributed along Sumatra and Java and to the north in Kalimantan (Fig. 2). A subset of the network is openly available, but the full network has restricted access. Seismic monitoring using a subset of the Indonesian stations is routinely performed by the GEOFON Program of GFZ (GeoForschungsZentrum) in Potsdam, and a moment magnitude (Mw) 5.1 event on 22 December 2018 was detected with a manually revised location of 105.44°E, 6.15°S at 13:55:48.7 UTC (GFZ: https://geofon.gfz-potsdam.de/eqinfo/event.php?id=gfz2018yzre). A body and surface wave inversion for the mechanism indicates Mw 5.1 strike-slip faulting with steeply dipping nodal planes (strike, ϕ1 = 124°; dip, δ1 = 81°; rake, λ1 = 7°; ϕ2 = 33°; δ2 = 83°; λ2 = 171°; centroid location: 105.67°E, 6.16°S). The geometry and size are not indicative of strong tsunami excitation. Other than a long-period source detection [seismic magnitude (Ms) 5.4 at 13:56:08 UTC at 105.75°E, 6.25°S] by the LDEO (Lamont Doherty Earth Observatory) Surface Wave Event Locations procedure (https://www.ldeo.columbia.edu/~ekstrom/Research/SWD/current/RADB_SWD_grd.html), to the best of our knowledge, no other seismic detection of the event was reported. This indicates that only local seismograms or long-period teleseismic seismograms had signals with sufficient signal-to-noise ratio to detect the source process. Given the relatively large magnitude estimates, the absence of routine teleseismic detection suggests that teleseismic body waves were unusually small (Mw = 5.1 earthquake typically produce clear short-period signals at teleseismic distances), indicating a source distinct from typical earthquake faulting.

Fig. 2 Regional broadband three-component ground motions from the Indonesian (IA) and GSN (AU.XMI, MS.BTDF) networks.

The map shows the source location (red star) at Anak Krakatau and the regional broadband seismic network, with stations to the northwest labeled in blue and stations to the southeast labeled in purple. The lower figures show vertical (left), north-south (center), and east-west (right) components of ground displacement in the frequency band of 0.01 to 5.0 Hz. For each station, distance (Δ) and azimuth (ϕ) from the source (star) are listed in parenthesis.

Important information about the tsunami source included tsunami back-tracing from four tide gauge recordings [Kota Agung and Pelabuhan Panjang in Sumatra, and Marina Jambu and Banten (Ciwandan) on Java] (Fig. 1), which triangulates a source in the area of Anak Krakatau. The tsunami travel times inferred from the GFZ origin time are from 31 min (Marina Jambu) to 57 min (Pelabuhan Panjang), with the shallow water depth of the Sunda Strait being responsible for the slow tsunami velocity. The inferred travel times are smaller than the predicted travel times by Giachetti et al. (12) by about 5 min, but are very consistent with recent modeling by Grilli et al. (14). Allowing for inaccuracy of shallow bathymetry models used in the back-tracing, this can be readily reconciled. The combined information from the seismic event detection, the tsunami back-tracing, and the visible eruption strongly indicate that a volcanic deformation process was responsible for generating the tsunami.

The likely explanation for the tsunami was provided by before and after synthetic aperture radar (SAR) images from the Advanced Land Observing Satellite 2 (ALOS-2) and Sentinel-1A/1B systems. SAR intensity images of ALOS-2 data display dramatic geomorphic changes in the southwestern part of Anak Krakatau between 20 August 2018 and 24 December 2018 (2 days after the tsunami) (https://www.gsi.go.jp/cais/topic181225-index-e.html). The Geospatial Information Authority of Japan estimated that an area of 2 km2 on the southwestern side of the island had been affected by flank collapse by 5 p.m. on 24 December (UTC).

Improved time resolution is provided by the Sentinel data, with Sentinel-1A images at 22:33:45 UTC on 10 December and 22:33:44 UTC on 22 December, ~8 hours after the tsunami first reaching the Java coast. The difference in radar backscatter intensity between the images on 10 December and 22 December indicates that the western and southern parts of the volcano collapsed in that time frame (Fig. 3). Using the digital elevation map (DEM) (DEMNAS, spatial resolution of 0.27 arc sec; freely available from http://tides.big.go.id/DEMNAS/#Metode) collected before the event, we estimate upper bounds on volume changes of 0.007 and 0.04 km3 for the western and southern subaerial failures, respectively. Estimation of the possible slide volume below sea level is very difficult, but likely smaller than 0.28 km3 in the scenario of Giachetti et al. (12). The 28 December image shows major subsequent modification of the volcano, which was not tsunamigenic. This image corresponds quite well with the 24 December ALOS-2 image, indicating that not much further modification of the edifice occurred after 24 December. Similar images of the Sentinel data are presented by Williams et al. (15), who estimate a 0.004-km3 subaerial failure volume. This is much smaller than the collapse volume inferred by Grilli et al. (14), who assert that the landslide removed 50% of the subaerial extent. Gouhier and Paris (10) estimate that the flank collapse removed 0.094 km3 of subaerial lavas, based on their interpretation of the eruptive sequence.

Fig. 3 Sentinel-1A images showing the changes of Anak Krakatau from 10 to 28 December 2018.

(A) SAR image acquired on 10 December 2018. (B) SAR image acquired on 22 December 2018, ~8 hours after the landslide generated tsunami. Possible flank failure occurred in the southern (outlined by dashed red line) and western (difference between red outlines and white outline to the west of Anak Krakatau). (C) SAR image acquired on 28 December 2018, including further modification of the edifice by large eruptions on 23 and 24 December. (D) Interpretation of the evolving geomorphology of Anak Krakatau superimposed on the SAR image on 22 December. Colored lines in (D) indicate the changes in island surface area with the color coding used in (A) to (C). The looking directions of these images are similar, around 44°.

The SAR images provide a strong case for moderate collapse of the western and southern portions of Anak Krakatau between the time of photographic images on the evening of 22 December and about 8 hours later. Flank collapse that was not accompanied by strong short-period seismic wave generation is thus the putative source of the tsunami, as inferred by many observers by the day after the tsunami. With quite small subaerial volume loss having occurred and large uncertainty in the volume of slide material below sea level, we seek to further constrain the overall process by examination of seismic recordings.

Inspection of Global Seismic Network (GSN), Federation of Digital Seismic Network (FDSN), and F-net [National Research Institute for Earth Science and Disaster Resilience (NIED), Japan)] stations guided by the GFZ event location confirmed that long-period surface waves were produced by the event at 13:55:48.7 UTC. Examples of the F-net recordings are shown in fig. S1. The body wave phases in these recordings are very weak, barely visible in broadband traces, clarifying why too few short-period detections were available for the USGS-NEIC to form an event. However, long-period surface wave arrivals are very clear in fig. S1. The broadband signals from the Indonesia network (IA) show similarly low levels of short-period body wave energy, but clear move-out of longer-period phases across Indonesia (Fig. 2). These ground displacement waves swept across the network within 8 min, with many, but not all, stations providing good three-component recordings. Note that the primary ground shaking is concentrated within a 100-s interval, with particularly simple waveforms on the horizontal components, which indicate strong transverse component motion aligned with the northeast-southwest direction.

We incorporate these signals into a regional W-phase moment tensor inversion using the basic algorithm of Kanamori and Rivera (16), adapting to bandwidth (40 to 200 s) of the regional data (Fig. 4A). This solution, representative of what would be obtained in conventional operational application of the regional W-phase inversion procedure, has a best double couple with very shallowly dipping (1.4° to the southwest) normal faulting (or nearly vertical dip-slip motion), with a seismic moment of 8.0 × 1017 Nm (Mw = 5.9) for a source depth of 3.5 km. The centroid location is 6.11°S, 105.42°E, with a centroid time shift of 25 s. The overall solution yields a near vertical nodal plane, but the strike and dip of the shallow-dipping plane are not tightly resolved. Waveform fits for this solution (fig. S2) are fairly good for some stations, but poor at others, and there is substantial long-period noise in many waveforms. If the screening criterion used in the W-phase inversion to remove noisy traces [see (16)] is too restrictive, then only vertical components are retained in the inversion and the sense of vertical dip-slip becomes indeterminate due to the limited azimuthal coverage provided by the data. Narrowing the bandwidth to 30 to 83 s to eliminate long-period noise and performing a moment tensor inversion still indicate a shallow-dipping normal-faulting solution (fig. S3), with a lower estimated seismic moment of 1.3 × 1017 Nm (Mw = 5.3), and better waveform stability and matches are found, although some signals are still not well matched and some have low signal-to-noise ratios. There is a strong trade-off between dip and moment, as expected for a shallow dip-slip source. The passband used for the latter solution is shorter period than for conventional W-phase inversions of larger events (periods of 200 to 600 s are typical), making the validity of a one-dimensional (1D) Earth model tenuous and accounting for the propagation effects of the shorter-period energy in the passband necessary. However, the overall dip-slip nature of the solution is quite stable and is different from the strike-slip solution obtained by GFZ.

Fig. 4 Regional waveform inversion solutions for the double-couple and single-force sources.

(A) Best double couple from the W-phase moment-tensor inversion of regional waveforms in the period band of 40 to 200 s. Corresponding waveform comparisons are shown in fig. S2. W-phase inversion for the period band of 30 to 83.3 s is shown in fig. S3. (B) Single-force source representation, showing azimuth and force strength with dip of 0.1° downward from the horizontal indicated by the negative sign, from inversion of regional waveforms in the period band of 30 to 83.3 s. Corresponding waveform comparisons are shown in fig. S4. For each solution, the number of traces used is indicated along with the source geometry and centroid location information. (C) Comparisons of centroid locations of the W-phase moment tensor solution from (A), single-force solution from (B), and the moment tensor solution from GEOFON (GFZ), along with the epicenter locations from GFZ using regional IA network and LDEO using global long-period surface waves. The magenta star shows an approximate location of Anak Krakatau, and dashed magenta arcs indicate distances of 10, 20, and 30 km from Anak Krakatau.

A southward directed slump on Anak Krakatau could conceivably involve normal faulting, but the depletion of short-period energy in the seismic records suggests a minor role for brittle faulting. As discussed by Dahlen (17), if the sliding material in the hanging wall disaggregates so that it does not internally convey elastic waves like a solid, then the double-couple force representation for conventional faulting may be better replaced with a single force, representing the reaction to the slide mass accelerating and decelerating on the Earth’s surface. This concept has been successfully applied to many landslide events, notably the 1980 Mt. St. Helens landslide [e.g., (18, 19)] and many subsequent events [e.g., (2024)]. Okal (25) discussed the difference in ground motion and tsunami excitation between single forces and double couples, noting in particular that the physical requirement that the integral of the single force-time history must vanish at the end leads to fall-off of the low-frequency source spectrum.

The improved waveform fit of the single-force source (fig. S4 versus S3) is an indication that a single force is a better kinematic representation of a landslide than a double couple. We extend the W-phase inversion algorithm to include a single-force source representation (see details in Materials and Methods) and apply it to 30- to 83-s period ground motion recordings of the regional IA network. We assume the single force-time history to be a one cycle sine wave. A nearly horizontal force (dip, 0.1°) with an azimuth of 43.2° and a peak force amplitude of 3.93 × 1011 N (Fig. 4B), applied on the Earth surface, accounts for many waveforms very well (fig. S4), better than for the double couple source (fig. S3). The half duration of 27 s corresponds to one-half of the period of the sine waves.

The double-couple and single-force solutions are consistent in terms of representing nearly horizontal displacement of hanging wall material toward the southwest. Assuming that the landslide effectively disaggregated the sliding mass, the moment tensor double-couple solution is not quite physically correct, and the single-force solution may be a more sensible representation of the equivalent force system for elastic wave generation process. Both solutions can be obtained very quickly as soon as the seismic waves have swept across enough stations to provide a stable solution. With automation, the entire IA network could be processed to give corresponding solutions within about 10 min after the event, as long as a source location is rapidly obtained, as was the case here. For other events, if the body waves are even more depleted in high-frequency energy so that an event location cannot be determined, then it may be possible to adapt the automatic analysis to detect and locate the hypocenter using long-period waves alone, before automatic W-phase inversion for moment tensor and single-force representations. Given the location and mechanism information, the likelihood of tsunami excitation can be evaluated, and an informed warning can be issued.

The regional inversion for a single-force source representation of the long-period ground motions (Fig. 4B) provides a first-order approximation of the force-time history of a one-cycle sine function, with a total duration of 54 s. The specified form of the force-time history imposes a specific peaked shape of the source spectrum that may affect the peak force and force orientation in the inversion. To better resolve the actual shape of the force-time history and to ensure that the single-force representation is correct to first order, we directly extract the force-time history from the IA recordings by deconvolving the observed waveforms with the Green’s functions for a given single-force geometry. We use a selected set of signals with very good signal quality for the passband of 8 to 125 s. The dip of the single force estimated by the regional waveform inversion has significant uncertainty (fig. S5), so we search for the favored dip angle by minimizing the amplitude ratio between the stacked deconvolved force-time histories from vertical and tangential records. We find a preferred dip angle of −12° (upward from horizontal) for a force azimuth of 42° (Fig. 5B and fig. S6) and a peak force of 6.1 × 1011 N acting on the Earth’s surface. We estimate at least 10% uncertainty in the peak force due to trade-off with dip and due to uncertainty in the Green’s functions. This is our preferred representation of the source, involving the reaction force from the slide mass moving down the southwestward flank with an average 12° dip.

Fig. 5 Stacked force-time history of the landslide and constraint on the dip angle.

(A) Red and blue curves are linear averages of force-time histories obtained by deconvolving tangential (T comp.) and vertical (Z comp.) broadband data shown in (C) and (D), respectively. The force geometry is a single force with dip −12° (dipping upward) and azimuth 42°, used to calculate Green’s functions for each station for the deconvolutions. Red dashed lines indicate the approximate time window of the force-time history. (B) Variation of peak force amplitudes for the tangential (red) and vertical (blue) components inferred from the first peak of the stacked deconvolutions as a function of the force dip angle. The force azimuth is held fixed at 42°. A dip of −12° gives the optimal agreement between the force estimates and is the preferred case. (C) and (D) show the observed (black) and reconstituted (magenta) waveform comparisons of broadband displacements in the period band of 8 of 125 s for tangential and vertical components, respectively. The right panels show the deconvolved force-time history for each trace. Inferred peak-force values (in 1011 N) are listed for each trace.

Using the preferred force geometry of azimuth 42° and dip −12°, the deconvolved force-time histories consistently show a leading ~30-s duration positive triangular pulse followed by a somewhat broader, 40- to 50-s duration negative pulse at all stations (Fig. 5). Shorter period oscillations in the force-time histories of some stations likely represent limitations of the 1D Green’s functions used in the deconvolutions. Stacking the estimates from the tangential and vertical components gives similar force-time function shapes and peak amplitudes. The overall duration of the force is ~80 s, somewhat longer than the 54-s duration indicated by the single-force waveform inversion. A gradually increasing amplitude ratio of vertical to horizontal average force-time histories evident in the stacks from ~ 15 to 70 s in Fig. 5 and fig. S6 suggests that the dip angle of the reaction force (sliding surface) may decrease following initial strong acceleration into the negative amplitude deceleration phase. Computation of Green’s function for a specific single-force strength and geometry is very rapid using a precomputed library (see details in Materials and Methods), so determination of detailed force-time histories can also be performed rapidly, once the regional waveform source inversion is completed.

DISCUSSION AND CONCLUSION

The reaction force, F, acting on the ground due to the downhill movement of the landslide mass estimated from the seismic wave analysis is a valuable quantity that can help to infer the sliding volume and basal friction [e.g., (26)]. For F = 6.1 × 1011 N, g = 9.8 m/s, ρ = 2000 kg/m3, and γ = 12°, we get an estimated lower bound of sliding volume V = 0.15 km3 by ignoring any basal friction. This estimate involves large uncertainties in F, ρ, and γ, so the uncertainty in estimating sliding volume is large. With most of the slide being below water, the effect of sea-water buoyancy affects ρ, so we use a relatively low effective density, but this is uncertain. If one can independently constrain V, then perhaps by imaging the submarine slump directly, it will be possible to estimate the basal dynamic coefficient of friction, μ, by viewing the force amplitude as the difference between the driving force (Mg sinγ) and the frictional shear force μMg cosγ [e.g., (26); fig. S7]. Current constraints on V are too uncertain to make a robust determination, but assuming V = 0.2 km3, an estimate of μ = 0.05 is obtained; if V = 0.3 km3, then μ = 0.1. Given the small subaerial failure volume and the edifice shape of Anak Krakatau (12), the underwater sliding volume is probably no more than 0.2 km3, so the effective basal friction is likely smaller than 0.05. This is much smaller than values of 0.3 to 0.4 for most landslides in Japan, with much smaller volumes of 0.002 to 0.008 km3 (24). Saturated conditions and deposit of the landslide mass onto sediments in the preexisting caldera may account for low average frictional resistance. The geometry of the caldera basin generated during the 1883 Krakatau eruption precludes long run-out of the slide, which constrains the overall duration to be shorter than most inland landslides (fig. S8B) as well.

Numerous other seismic modeling procedures have been developed to characterize seismic radiation for landslides or to invert for point-force representations of landslide sources [e.g., (20, 21, 2729))], so the use of the W-phase approach is primarily motivated by the widespread use of this code for regional earthquake faulting determinations and tsunami warning. Comparing the findings for the Anak Krakatau flank collapse with the compilation of landslide measurements for other events (21) indicates that the force amplitude is near the center of those observations, but the runout duration is near the lower end of observations (fig. S8). On the basis of the Ms – log F pattern from Ekström and Stark (21), a surface wave magnitude of about 5.3 is predicted, consistent with the 5.4 value from LDEO surface wave detection. Among the history of large landslides that have produced tsunami discussed by Giachetti et al. (12), the Anak Krakatau event involves a small slide volume, even allowing for uncertainty in the extent of below sea level deformation. Smaller slides of this size occur more frequently than huge slides, so improving the warning capabilities for future events is very desirable.

The Anak Krakatau tsunami disaster presents difficult challenges for mitigation of the tsunami hazard for such events. Figure 6 provides a nonlinear timeline of critical points in the sequence. Given the 2012 recognition of a potential landslide hazard, the most straightforward approach to a reliable tsunami warning system seems to be to deploy a continuous ocean bottom pressure sensor near the potential source region, with telemetry to land to detect any large tsunami wave generated by landslide and/or faulting. However, this is a fairly difficult type of system to operate and maintain, and the option of putting it on an undersea cable to provide power and real-time communications is definitely expensive. Deployment of telemetered seismometers on the targeted volcano can provide helpful information, and a seismic station was operating on Anak Krakatau before the collapse, but given the ongoing eruptions for several months, local recording can be difficult to interpret and may not resolve long-period energy well. Continuous analysis of a regional seismic network is perhaps the most cost-effective approach, but as in the case of the 2018 collapse, it is important to rapidly apply robust long-period ground motion analysis when an event is detected. As shown here, rapid analysis of the long-period recordings from regional seismic stations can resolve a reliable source representation of the event. If a shallow dip-slip moment-tensor source with depleted short-period body waves is observed, then the single-force inversion can be performed to evaluate the possibility that a landslide is involved rather than a tsunami earthquake. In either case, the source geometry, strength, and duration can provide a basis for deciding whether to issue a tsunami warning. It is then critical to have a means of conveying the warning to remote areas and having a population informed of how to react when a warning is broadcast.

Fig. 6 Time sequence of geological activities on Mt. Anak Krakatau.

A nonlinear timeline is shown, extending from the first emergence of Anak Krakatau above sea level in 1929 through its construction of a ~300-m-high edifice by 22 December 2018 and its subsequent landslide and eruption resulting in large mass loss by 28 December 2018. Key times and intervals of processes related to the buildup and follow-on to the landslide along with key observation times from seismic network and Sentinel satellite are noted. The interval in which infrared satellite imagery of eruptive perturbation of the atmosphere occurred is indicated by the horizontal red line.

The 2018 Anak Krakatau volcano flank collapse appears to have involved a relatively small slide volume (<0.2 km3), yet it produced a deadly tsunami affecting local coastlines due to the lack of warning. The collapse did not produce strong short-period ground motions, but is very well recorded by broadband seismic stations throughout Indonesia. These data, all available within 8 min of the event, provide sufficient information to rapidly determine the source process by locating the event and inverting the long-period (>40 s) ground motions for equivalent body force systems that can represent shear faulting or a single-force geometry and force strength in the case of shallow-dipping landslide. Finding a 50- to 80-s duration force-time history for a surface horizontal single-force solution with Ms ~5.4 to Mw 5.9 would be immediately suggestive of an anomalous source process that can be tsunamigenic depending on location. Thus, use of modern seismic analysis procedures can contribute to tsunami warning even for submarine slumping events depleted in short-period energy. Similarly, rapid determination of moment tensors for shallow tsunami earthquakes is more viable and inexpensive compared to real-time offshore seafloor monitoring systems. In Japan, automated centroid moment tensor solutions using W-phase recorded on regional broadband network are obtained in 15 min after the event and used to update the tsunami warning (30). Ideally, both approaches would be used given the loss of life caused by this modest volume collapse.

MATERIALS AND METHODS

Analysis of the sentinel data

Good time resolution was provided by the Sentinel data, with Sentinel-1A images at 22:33:45 UTC on 10 December, 22:33:44 UTC on 22 December (~8 hours after the tsunami first reaching the Java coast), and 22:33:06 UTC on 28 December, all on descending track 47. Choosing the image on 10 December as a master image and co-registering the subsequent two, the images were multi-looked with a factor of 1 in azimuth and 4 in range, resulting in a spatial resolution of ~5 m. We also obtained one Sentinel-2A optical image collected on 16 November 2018. The multi-looked SAR images (Fig. 3) were geocoded by selecting ground control points on the Sentinel-2A image. The difference in radar backscatter intensity between the images on 10 December and 22 December indicates that the western and southern parts of the volcano had collapsed. The areas were estimated to be 0.4 km2 in the western part (difference between the yellow and magenta lines) and 0.6 km2 in the southern part (area outlined by the red dashed curve). Lacking a precise immediate post-event DEM, it is difficult to precisely determine the corresponding subaerial volume changes, but using the DEM collected from before the event, we can estimate upper bounds on the volume changes. This gives estimated upper limit values of 0.007 km3 for the western failure and 0.04 km3 for the southern failure.

Source inversion of regional broadband seismic data

We incorporated regional broadband ground motions from the IA network into a regional moment tensor inversion using the basic W-phase algorithm of Kanamori and Rivera (16). The bandwidth used for the regional data is 40 to 200 s (Fig. 4B and fig. S2), which was routinely used for the regional W-phase inversion. Green’s functions were calculated by a frequency–wave number method (31) using a regional 1D average velocity model for the Japan subduction zone. The solution has a best double couple, with very shallowly dipping (1.4° to the southwest) normal faulting (or nearly vertical dip-slip motion), with a seismic moment of 8.0 × 1017 Nm (Mw = 5.9) for a source depth of 3.5 km. The centroid location is 6.11°S, 105.42°E, with a centroid time shift of 25 s. The overall solution is stable for a range of data screening parameters, but the strike and dip of the shallow-dipping plane are not tightly resolved. If the screening criterion is too restrictive, then only vertical components are retained in the inversion, and the sense of vertical dip-slip becomes indeterminate due to the limited azimuthal coverage provided by the data.

To avoid drift of the ground motion records, a narrower passband of 30 to 83 s was also used for inverting the regional data. The regional moment tensor inversion then has a best double couple with shallowly dipping (14.1° to the west) normal faulting (or nearly vertical dip-slip motion), with a seismic moment of 1.3 × 1017 Nm (Mw = 5.3) and a source depth of 3.5 km. The centroid location is 6.106°S, 105.42°E, with a centroid time shift of 20 s. The solution and corresponding waveform fits are shown in fig. S3.

Modification of the W-phase inversion method for a single-force source

The regional waveform inversion algorithm was extended to include a single-force source representation, and we briefly describe how we modify the W-phase inversion method (16) for a single force source. We replace equation 8 of Kanamori and Rivera (16) byUi,k(t)=j=13Fjgi,kj(t)S(t)(1)

Here, U1,k(t), U2,k(t), and U3,k(t) are, respectively, north-south (NS), east-west (EW), and the vertical component of displacement as a function of time t at station k.g1,k1(t), g2,k1(t), and g3,k1(t) are, respectively, NS, EW, and the vertical component of displacement at station k due to a delta function force in NS direction at the source. gi,k2(t) and gi,k3(t) are similarly defined for a delta function force in EW and vertical directions, respectively. We call gi,kj(t) the Green’s functions. F1, F2, and F3 are the magnitude of the force in NS, EW, and vertical directions, respectively. The symbol ⊗ means convolution. S(t) gives the shape of the source time function given byS(t)=sin(πt/td) for 0t2td and S(t)=0 elsewhere(2)where td is the centroid time. This form of S(t) is appropriate for a landslide type source with acceleration and deceleration stages.

With this replacement, if we compute the Green’s functions for a single force, then we can use the standard W-phase code for single-force source to determine Fj (j = 1, 2, and 3), the centroid time, and location by the least-squares method and grid search (16, 32).

In practice, however, we use an integrated form of Eq. 1.Ui,k(t)dt=j=13Fjgi,kj(t)[S(t)dt](3)

The reason for using this form is that in the standard W-phase inversion code, we use a positive triangular moment rate function for the source, and using a positive source function S(t)dt=tdπ[cos(πtdt)1] requires a minimum modification of our code. Theoretically, using Eqs. 1 or 2 should make no difference, except that Eq. 3 emphasizes longer period component than Eq. 1.

With this modification in the method, we applied it to the IA ground motion recordings. The displacement records are first filtered in the period band of 30 to 83 s and integrated before inversion. The source function so obtained should then be differentiated to recover the force-time history. Green’s functions are obtained from the online precomputed responses for the IASPEI91 model in the IRIS Synthetics Engine (Syngine: http://ds.iris.edu/ds/products/syngine/). The inversion of 51 channels from 33 stations gives the single-force solution shown in Fig. 4A. A nearly horizontal force (dip 0.1°) with azimuth of 43.2° and a peak force amplitude of 3.93 × 1011 N, applied on the Earth surface, accounts for the waveforms well (fig. S4). The centroid location is 6.14°S, 105.42°E after a grid search on candidate locations [e.g., (32)]. The half duration of 27 s corresponds to the duration of each half-sine arch. To evaluate uncertainty of this inversion, we performed a bootstrap analysis. By randomly selecting subsets of data from a master set of 98 waveform traces with high signal-to-noise ratios, 1,000,000 inversions were performed, yielding solution distributions for a force amplitude of (4.4 ± 0.21) × 1011 N, a dip angle of 2.3° ± 3.02°, and a strike angle of 41.8° ± 2.55°, with a 95% confidence level.

Single-station deconvolutions for estimating force-time history

To better resolve the force-time history and to ensure that single-force representation is correct to first order, we directly extract the force-time history from the IA recordings by deconvolving the observed waveforms using the Green’s functions computed for a specified force geometry. Green’s functions for unit forces acting in the vertical, NS, and EW directions were precomputed using a frequency–wave number method (31) with a time sampling of 0.25 s. For a unit single force with a specific geometry, we can calculate corresponding Green’s functions using vector summation of the precomputed Green’s function library. A regional 1D velocity structure, adapted from Crust 1.0 (33), was used in computing the Green’s functions, so some variation among the deconvolved force-time histories is expected because of incorrect propagation effects, but the results are quite stable for a large number of stations. Tangential and vertical component displacement records filtered in the period band of 8 to 125 s were integrated for each station and deconvolved by the Green’s functions with a positivity constraint, and then the resulting integrated force-time history was differentiated and convolved with the Green’s function to give the comparisons shown in Fig. 5 (C and D). Note the difference in crustal model and single-force dip from those obtained from the single-force waveform inversion described above.

Estimation of slide volume

Considering the rock mass, M, to have moved down a slope with dip angle of γ, the reaction force, F, acting on the ground during the motion isF=Mg sinγμ Mg cosγ=ρVg(sinγμcosγ)(4)

Thus, we can estimate the slide volume, V, byV=F/[ρg (sinγμcosγ)](5)

Using the force amplitude of 6.1 × 1011 N from seismic data analysis and assuming an effective density ρ = 2000 kg/m3, allowing for reduction due to the buoyancy of sea water, we compute the variations of sliding volume with basal friction and dip angle of the force (sliding surface) shown in fig. S7. The lower bound of sliding volume is ~0.15 km3 for our preferred dip angle of 12° without basal friction. If the basal friction is 0.05, then the sliding volume is ~0.2 km3.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/3/eaaz1377/DC1

Fig. S1. F-net ground displacements in Japan for time windows of 370 to 2000 s after 13:55:48.7 on 22 December 2018 (UTC).

Fig. S2. Regional W-phase waveform fits for the moment-tensor source in Fig. 4A.

Fig. S3. Moment-tensor inversion using regional waveforms in the passband 30 to 83 s.

Fig. S4. Waveform fits for the single-force source model in Fig. 4B.

Fig. S5. Bootstrap results for 1,000,000 single-force regional waveform inversions.

Fig. S6. Stacked force-time history for variable dip angles.

Fig. S7. Estimated slide volume as functions of average basal friction and average detachment dip angle.

Fig. S8. Comparison of key estimates of the 2018 Anak Krakatau landslide with other landslides.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: We thank three anonymous reviewers for their comments on the manuscript, which helped to improve the presentation. Funding: This study was supported by Fundamental Research Funds for the Central Universities (to L.Y.), National Natural Science Foundation of China (no. 41874056 to L.Y.; no. 41874020 to Y.Z.), and the NSF (grant EAR1802364 to T.L.). Author contributions: L.Y. and H.K. conceived the project; L.Y., H.K., and T.L. designed the single-station deconvolution procedure; H.K., L.R., and L.Y. conducted moment-tensor and single-force inversions; Y.Z. performed SAR image analysis of Sentinel data; D.S. collected IA data; and all coauthors wrote the manuscript collaboratively. Competing Interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. The Incorporated Research Institutions for Seismology (IRIS) data management service (DMS) (http://www.iris.edu/hq/) was used to access the long-period (LH) seismic data from Global Seismic Network and Federation of Digital Seismic Network stations. Japan F-net broadband data are available from the National Research Institute for Earth Science and Disaster Resilience (NIED) (http://www.fnet.bosai.go.jp). Indonesia seismic data (IA-Net) are available from the Agency for Meteorology Climatology and Geophysics (BMKG) upon request. The Sentinel-1 images were obtained from the European Space Agency (https://sentinel.esa.int). The pre-event DEM is the Indonesian National DEM (DEMNAS), which is freely available from http://tides.big.go.id/DEMNAS/#Metode. Additional data related to this paper may be requested from the authors.
View Abstract

Stay Connected to Science Advances

Navigate This Article