Research

The PMCC Algorithm Method

For my research, I conducted array processing using the Progressive Multi-Channel Correlation (PMCC) algorithm (Cansi, 1995) to detect hydroacoustic and infrasound signals at each station to determine the back-azimuth of arrival. PMCC is a time-domain cross-correlation method that calculates signal arrival time delays between the sensor pairs at a specified array and measures the consistency of the signal between the sensor triplets using the relation:

If a signal has a consistency r that is below a specified value between i,j,k sensors, the PMCC algorithm regards it as a coherent detection and temporarily records the wave characteristics. The algorithm repeats this calculation for a number of additional time windows and frequency bands until a full list of detections with corresponding arrival times, arrival azimuths, frequencies, velocities, amplitudes, and other characteristics is generated. PMCC progressively adds additional subsets of array stations to the calculation. If the addition of progressive subsets results in strong variation of the azimuth, velocity, or time for a detection, PMCC rejects the detection, but if all subset combinations are consistent the detection is logged as a ‘pixel’. The PMCC method is the signal detector in use by the IDC and has frequently been used for processing of infrasound data, but has not commonly been applied to hydroacoustic data. Therefore in preliminary work for this study we tested the method using published eruption case studies to calibrate the parameters specifically to hydroacoustic data (e.g., choosing appropriate wave velocity, frequency, number of bands, etc.), making this a good demonstration of the PMCC method for hydroacoustic data.

The International Monitoring System Network

The newly constructed International Monitoring System (IMS) infrasound and hydroacoustic networks are valuable sources of global infrasonic and hydroacoustic data. The Comprehensive Nuclear Test Ban Treaty (CTBT) of 1996 proposed a ban on nuclear detonations in all environments, with the mandate that a global system of monitoring stations would be installed as part of the compliance and verification regime, and would consist of seismic, infrasound, hydroacoustic, and radionuclide stations. The station locations were chosen for optimal global coverage to detect nuclear explosions as small as 1 kiloton anywhere on Earth (Christie and Campus, 2010), but are also useful for detecting volcanic explosions since these occur in a similar frequency band.

Each IMS infrasound station consists of 4 or more infrasonic sensors forming an array with 1-3 km aperture. The sensors are high-sensitivity microbarometers which record fluctuations in micro-pressure caused by infrasound waves propagating through the atmosphere at a sampling rate of 20 Hz. The IMS hydroacoustic network contains 11 hydroacoustic arrays composed of 1-2 sensor triplets, spaced throughout Earth’s oceans. The stations used in this study are moored hydrophones with a sampling rate of 250 Hz that are positioned at the depth of the SOFAR channel axis (approximately 0.5 to 1.5 km deep, depending on latitude) to intercept long-distance hydroacoustic signals.

Because the true extent of seafloor volcanism is unknown, there may be thousands to millions of undiscovered seafloor volcanoes which have not yet been identified. Volcanic infrasound is now regularly detected, and hydroacoustic signal is occasionally detected, at IMS stations and used to determine eruption characteristics such as locations, durations, and eruptive processes (e.g., Fee et al., 2010, Matoza et al., 2011a,b, Green et al., 2013, Metz et al., 2016, and Matoza et al., 2018).

Anak Krakatau (“child of Krakatau” in Indonesian) is a highly active volcano located on the rim of the caldera formed by the 1883 paroxysmal eruption of Krakatau in the Sunda Strait, Indonesia. The Krakatau volcanic complex is the result of arc volcanism caused by subduction of the Australian plate beneath the Sunda plate, directly below Krakatau, and crustal thinning caused by the extensional faulting and rifting of the Sunda Strait (Dahren et al., 2011). The 1883 eruption culminated in a sequence of large, magmatically-driven explosions that generated catastrophic pyroclastic flows and a tsunami in the Sunda Strait, which together are thought to have killed at least 36,000 people (Auker et al., 2013). The eruption destroyed most of the volcanic edifice, forming a partially-submarine caldera. Anak Krakatau began forming on the northeast rim of the caldera during an eruption in 1927, and first emerged above sea level in 1928 (Stehn et al., 1929). Since its rebirth in 1929, the volcano has erupted frequently, with at least 40 recorded episodes of Strombolian to Vulcanian style eruptive activity characterized by heightened seismic activity, explosions, ash plumes, incandescent ejecta, lava flows, and fire fountaining (GVP 2019). By October 2018, the cone had grown to 338 meters above sea level (GVP 2018a). Anak Krakatau was mainly active on its southwest side toward the center of the 1883 caldera, making the southwest volcanic flanks unstable. Numerical simulations published in 2012 indicated that a hypothetical flank collapse toward the southwest would generate a tsunami that would reach the islands surrounding Anak Krakatau within 1 minute, and would reach the coastal cities of Sumatra and Java within 35-45 minutes of the collapse (Giachetti et al., 2012).

Major 2018 Eruption and Recent Activity

After 15 months of quietude, a new eruptive phase began on June 18, 2018, indicated by an increase in seismic activity (GVP 2018a). Ash plumes from the volcano were first visible on June 21, 2018, followed by the first incandescence ob- served on July 1, 2018. Lava flows reached the sea on September 22, 2018, and frequent explosions, ash plumes, and incandescent material were observed through December of 2018. On December 22, 2018, there was a significant increase in eruptive activity, with 423 explosions observed in one six-hour period from 12:00 to 18:00 local time by the Volcanological Survey of Indonesia. At 21:03 local time, the southwest side and summit of the volcano collapsed into the ocean as a result of the eruption, resulting in a tsunami first detected at 21:27 local time. The tsunami hit the coastlines of Bantan and Lampung, traveling up to 330 m inland with a maximum surveyed runup of 13.5 m (Muhari et al., 2019). Estimates by Ye et al. (2020) indicate that the flank collapse slide volume was small (<∼0.2 km3), and did not produce short-period seismic waves strong enough to trigger a tsunami warning from the existing Indonesian tsunami warning system which was designed to detect tsunamis from earthquake sources (Strunz et al., 2011). Indonesian authorities reported 437 deaths, and tens of thousands of injuries as a result of the tsunami (GVP, 2019). As a result of the December 22 partial collapse, the vent was submerged under tens of meters of seawater. Surtseyan activity and base surges were observed from December 23, 2018 through January 9, 2019 when a rim of tephra formed around the vent. The rim created a barrier between the vent and the ocean, forming a crater lake above the vent. The flank collapse reduced the height of the volcano from 338 m to 110 m. Since the collapse and tsunami, eruption activity has continued at a smaller scale. Intermittent ash plumes, explosions, and seismic events recorded on Anak Krakatau’s seismic network were reported throughout 2019 to April 2020, and in mid-April 2020 Strombolian activity accompanied a lava flow from the crater which filled the crater lake, covering the vent (GVP, 2020). The eruption is continuing at the time of writing (January, 2021).

Anak Krakatau 2018 Infrasound Data and Chronology

For my research, I analyzed data from the three closest infrasound stations to Anak Krakatau (stations IS06, IS07, and IS52, at distances of 1,156 km, 3,475 km, and 3,638 km from the volcano, respectively). Data was analyzed from the start of the new eruptive phase in June 2018, through January 2019 in this study (Rose, 2020). The climactic eruption phase beginning on December 22, 2018 is clearly visible in the waveforms, spectrograms, and array processing detections at the closest station, IS06.

To investigate waveform similarity, I used a network STA/LTA method to detect events within the continuous IS06 infrasound signal. Network STA/LTA calculates the ratio between the short-term average and the long-term average amplitude of the signal at each sensor in the network in moving time windows to automatically pick events when amplitude ratios that exceed a specified threshold value. All event triggers which did not fall within the time window of a PMCC bulletin detection ±10 seconds at IS06 were discarded. This method was used to identify all events which could be associated with eruption activity at Anak Krakatau, without regard for event waveform shape similarity. This network STA/LTA plus PMCC cross-check method was applied to 10 days of picked infrasound events from December 19 - 29, 2018.

This method produced an event detection list of 10,324 events, which were then used for waveform cross correlation analysis. A period of weak, intermittent activity was present from December 19 to December 22, during which the satellite imagery shows the volcanic edifice and cone was still intact. Then on December 22 the events became strong and consistent with increased amplitudes, corresponding to a period of intense Strombolian activity reported at the volcano. Event amplitudes decreased sharply with the flank collapse, with the first post-collapse satellite imagery 8 hours later showing that at least one large landslide had carried a large part of the flank and the cone of the volcano into the ocean. After the flank collapse, the eruption entered a period with a continuous, pulsing ash plume, explosions, base surges, and Surtseyan phreatomagmatic activity (GVP, 2019).

I then calculated event similarity using cross-correlation coefficient (pxy) for each of the 10,324 waveforms included in the IS06 event detection list, using the cross correlation function:

Using this function as a measurement of similarity, we employ an algorithm which cross-correlates each waveform with every other waveform in the event list, normalizing the waveforms to quantify waveform shape but disregard the amplitude of the events. This allows quantitative comparison of the source variability of the events, but not of their magnitude. I log the maximum correlation coefficient value between each pair of signals and the lag l at the maximum correlation, and use the cross-correlation values to build a maximum correlation coefficient matrix with m rows and m columns, where m is the number of events in the list.

We then clustered the IS06 infrasound events into groups of similar signals, or “multiplets” based on their cross-correlation values, using a clustering algorithm based on the method of Matoza and Chouet (2010). The threshold for cross-correlation was set at 0.65, to balance identifying as many waveforms belonging to a multiplet as possible, but excluding those with poor signal- to-noise ratios. Using the cross-correlation threshold of 0.65, the algorithm generated a total of 574 multiplets which contained 8,722 (84.5%) of the 10,324 IS06 infrasound events.

To ensure that similar waveforms were classified together despite poor signal-to-noise ratios, the algorithm further grouped the multiplets into families. For this procedure, the mean waveform (or “stack”) for each multiplet is calculated, and each multiplet stack is cross-correlated with each other stack in the same manner used to cross-correlate individual events. Multiplets with stacks that correlated with other multiplet stacks above a cross-correlation threshold of 0.75 were merged into “families”. The waveforms in the family were ordered chronologically and aligned to the most correlated event. The families classification procedure generated 487 families, with the four largest families containing 5394 (61.8%) of the 8722 multiplet events.

Anak Krakatau 2018 Hydroacoustic Data and Chronology

Of the three IMS hydroacoustic stations located in the Indian Ocean, we used only H01W and H08S in this study, located at distances of 3,307 km and 3,720 km from Anak Krakatau, respectively. Both stations consist of a triplet of hydrophones in a triangular configuration with ∼2 km horizontal separation between each station, moored at or near the SOFAR channel axis. We excluded hydrophone station H04 due to preliminary analysis showing no identifiable hydroacoustic signal from the Anak Krakatau direction recorded at this station.

We analyzed hydroacoustic and infrasound data for the 8 months from the beginning of the new eruptive phase in June 2018 through January 2019 using waveforms, spectrograms, and PMCC array processing to detect coherent signal from the direction of Anak Krakatau. PMCC results confirmed that these stations did not record the main eruption or flank collapse. However, detection results show a 12 day long swarm of nearly continuous hydroacoustic signals arriving at H08S from the back-azimuth of Anak Krakatau (89.4 degrees back-azimuth) beginning 24 days before the flank collapse and tsunami, from November 29 through December 11, 2018. At station H01W, the majority of PMCC bulletin detections were attributed to southern Indian Ocean iceberg cracking and breakup, but hydroacoustic detections were present in the direction of Anak Krakatau from November 30 - December 4, 2018, consistent with the detections at the closer hydroacoustic station, H08S. A shorter signal duration may have been captured at station H01W due to lower SNR or scattering by bathymetric obstructions such as seamounts, ridges, and the Cuvier Plateau in the path from volcanic source to receiver.

We also investigated hydroacoustic event similarity to determine whether repetitive signals could be identified and interpreted for information on submarine eruption source processes at Anak Krakatau. Twelve days of hydroacoustic data from November 29 to December 11, 2018 hydroacoustic swarm were beamformed and filtered with a 1-4 Hz Butterworth bandpass filter. Hydroacoustic events were picked manually from the filtered and beamformed signal, and picks were only made where there were concurrent PMCC detections from Anak Krakatau.

Using the same procedure as with the infrasound data, we cross-correlated hydroacoustic events to determine event waveform similarity. We then grouped the hydroacoustic signals into multiplets based on their cross-correlation values, setting the threshold for cross-correlation to 0.6. Multiplets were further grouped into families to show more broadly which multiplets were similar to each other, and to reduce the effect of the lower signal-to-noise ratio of the hydroacoustic events. Multiplets with stacks correlating with other multiplet stacks above a threshold of 0.6 were assigned into families. The events in each family were again cross-correlated to determine which event was most similar to the most other events, and the waveforms were ordered temporally and aligned to the most correlated event in the family. The families classification procedure yielded 153 signal families.

Due to the challenges of data collection and visual observation in a submarine environment, many questions remain regarding source mechanisms and conditions of submarine explosive volcanism. Therefore, the discussion presented here is speculative and represents one set of possible interpretations of eruption events at the volcano. Several factors must be considered when investigating the physical source mechanisms which produced the recorded hydroacoustic signal, which include the frequency, duration, continuity, recurrence interval, waveform similarity, and visual observations at the ocean surface. There is no evidence of harmonic tremor in the H08S hydroacoustic swarm signal. Since most of the hydroacoustic swarm events are < 30 seconds long, we classify them as impulsive signals. These impulsive events may be phreatomagmatic explosions, bursting magmatic gas bubbles at the opening of a conduit on the flank, or volcanogenic earthquakes coupled into the water column. Based on the variability in shape and duration of the hydroacoustic waveform families, it is possible that any or all of these processes were occurring within the same timeframe. Since tectonic earthquakes are shown not to be a plausible cause of these impulsive events and a volcanic source is likely, this analysis shows that remote hydroacoustic signal analysis may be useful in identifying volcanic events that were not flagged by other geophysical methods.