We developed and implemented our Machine Learning (ML) detection algorithm for pumice raft detection on the Google Earth Engine (GEE) platform (Gorelick et al., 2017). GEE is a web-based, publicly available platform that enables access to a vast catalog of satellite images and the resources to run global-scale analyses without the need to download or export large amounts of data. There are various satellite collections offered through GEE, such as low resolution (MODIS) and medium-high resolution imagery (Landsat, Sentinel-2). Although some other super-high-resolution image collections are available outside of GEE (e.g., Planet labs—3 m/pixel, Digital Globe—50 cm/pixel), they are typically not publicly available without commercial licenses. Thus, for this study, we have primarily focused on using GEE resources for the ML algorithm.

Specifically, we use GEE collections from the Sentinel-2 Multi-Spectral Instrument (MSI) as our baseline satellite product. Sentinel-2 (a pair of two satellites, each with MSI instrumentation) offers both high-resolution imagery (10–60 m/pixel), good coverage in regions of interest, and a relatively frequent repeat time (∼global 5-day revisits; Supplementary Text S6). Sentinel-2 data products are also freely available through the European Space Agency’s Copernicus Open Access Hub as well as other cloud environments. For our study, we chose to use Sentinel-2 as its high resolution imagery could be used to detect much smaller rafts than a lower resolution satellite (e.g., MODIS). In addition, Sentinel-2’s MSI collects data across 13 different spectral bands, with finer spectral coverage than other high resolution satellite image collections (e.g., Landsat 7 and 8) (See spectral response curve for Landsat 8 image of Puyehue-Cordón Caulle pumice in Figure 1B). An initial method using thresholds on only the visual bands to detect pumice rafts was insufficient, so the additional spectral bands are necessary in our ML algorithm (Supplementary Text S1). As illustrated by the variable importance in the Random Forest classifier (RF, Supplementary Figure S11), the multi-wavelength information is critical for accurate classification with a dominant role of the visible bands. In particular, the reflectance of pumice in the near-infrared (NIR) bands (700–900 nm) is much higher than the reflectance at those bands of most biological phenomena that are near the surface but still underwater (i.e., algal blooms, coral spawn, seaweed) (Biermann et al., 2020; Qi et al., 2020), so the NIR bands are critical for distinguishing between pumice and other visually similar-looking biological blooms. Our overall methodology is general and can be applied to other satellite collections in the future (Supplementary Text S10).

To identify spectral characteristics that can be used to classify Sentinel-2 image pixels as pumice rafts, we generated spectral response curves for pumice and other categories of interest in Figure 1B. Spectral response curves record the mean reflectance or brightness of an image pixel for a range of wavelengths. We used the Tonga pumice raft from 11 August 2019 to generate the spectral response curves (Figure 1B), as the particular eruption and the associated raft has been extensively analyzed by previous work (Brandl et al., 2020; Jutzeler et al., 2020). We also show the variance around the mean spectral response curve calculated for all of the pixels for each class (pumice, water, light clouds).

A key result from this analysis is that there is a significant difference between the spectral response curves of pumice, water, and light clouds (Figure 1B). Additionally, we find relatively minor (compared to differences with other classes) variation in the reflectance from pumice pixels within a single geo-temporal area, such as a specific day in Tonga (Figure 1B) or comparing across multiple days for the same raft (Supplementary Figure S16). Although there is some variation in pumice spectral response curves when comparing rafts from different chemical compositions, sources, and times (Figure 1B, comparison with Rabaul raft and Puyehue-Cordón Caulle raft), the general shape of the reflectance curve remains very similar. This characteristic shape of the spectral response curve for pumice pixels allows for an algorithm to identify pumice and differentiate from other classes (e.g., water, clouds) across a broad range of regions and time periods. Details for the Puyehue-Cordón Caulle raft are provided in Supplementary Text S13.

Our machine-learning algorithm uses a Random Forest classifier to read in an RGB Sentinel-2 image and return a classified image, where each pixel is colored according to the assigned class. The algorithm specifics are detailed in Supplementary Text S4. Since RF is a supervised learning algorithm, we need to train it on a set of manually demarcated and classified pixels. Our primary training data for pumice, ocean water, light cloud cover, and heavy cloud cover was sampled from the Tonga raft on 11 August 2019 (Figure 2A, only a small part of the raft pixels were used for training). We also included additional data from a Sentinel-2 scene of Rabaul, Papua New Guinea, on 20 April 2020 (Figure 2B, spectrally this is representative of the potential rafts from other days also). This image includes a large, distinct pumice raft as well as ocean water, light cloud cover, pumice mixture classes, and two different discolored water classes (additional information for the discolored water classes are included in Supplementary Text S9). Since the discolored water classes are not the primary focus of this study, our primary optimization for the RF algorithm was to ensure accurate detection of pumice rafts. We provide all the scripts, with step-by-step instructions for usage, used in our analysis in a publicly available repository (Supplementary Material—GEE Script Links).

We applied our classification algorithm to Sentinel-2 images from different geo-temporal regions to test model accuracy (Figure 2). In the Tonga area on 11 August 2019 (Figure 2A), the classifier displays pumice pixels in red, water in blue, light cloud cover in orange, and heavy cloud cover in white. The shape of the large raft is distinctly visible in the classified image. In the Rabaul region, on 20 April 2020 (Figure 2A), the classifier also includes additional classes: mixed/faint pumice—a mixture of water and pumice—shown in light blue, and two different classes of discolored water shown in turquoise and magenta. Overall, our algorithm is efficient at identifying pumice from other backgrounds though the accuracy of discolored water detection (mixing with corals) and shallow cloud is not great at present. Algorithm validation methods and results are included in Supplementary Text S5.

To assess the utility of our algorithm for new submarine eruption detection, we applied the classifier over a single region for an extended period of time. We focused on Rabaul, a partially submarine volcano located on the Gazelle Peninsula’s tip at the northeast end of New Britain in Papua New Guinea (Figure 3A). The Rabaul caldera (∼8 × 14 km size) was formed as a consequence of multiple large explosive eruptions in the past few hundred thousand years, with the present day shape due to an eruption ∼1,400 years ago (GVP and wunderman, 1994). The caldera is mostly shallow submarine ( < 200 m water depth) and is connected to the sea on the east through a wide opening (Blanche Bay). The main raft-forming eruptions for this volcano occurred in 1878, 1937, and 1994, and no raft formation has been recorded since 1994 (GVP and Wunderman, 1994; GVP and Wunderman, 2006). No activity has been recorded at either of the main vents (Vulcan and Tavurvur) since 2014 (Bernard and Bouvet de Maisonneuve, 2020). More detailed eruption history is provided in Supplementary Text S8.

In the Rabaul area, we applied our algorithm from November 2015 (start of the Sentinel-2 coverage for the Rabaul region) to August 2020—a total of 239 distinct days with images. More details on our algorithm application method are included in Supplementary Text S7.

Of these 239 days, we found that 74 days were too cloudy for the classifier to detect any pumice meaningfully. Cloudy days were filtered out by manually examining classified images and removing images in which every pixel was labeled as heavy or light cloud cover. In the future, this step can be automated by explicitly filtering the images based on the classified heavy cloud fraction. We detected potential rafts in 28 (red lines, Figure 3B) of the remaining 165 days (gray lines, Figure 3B), leading to a detection rate of 16.97%. As illustrated in Figure 3B, most of our raft detections were after January 2018 (Figure 3B). This is likely a consequence of increased revisit frequency (∼5-day) after the second Sentinel-2 satellite launch. Before 2018, when only one Sentinel-2 satellite was in operation, there are significantly fewer images available. It is noteworthy that none of the rafts detected in our analysis had been previously reported in the scientific literature (to the best of our knowledge) or the Smithsonian Global Volcanism Catalog (Global Volcanism Program, 2013). The sizes of our detected rafts varied greatly, with raft areas as small as 20 km² to as great as 10,000 km². At present, we do not have any direct ground truthing of our detections. However, considering the distinctive spectral features of pumice compared to other biological sources, especially in the Near Infrared (Biermann et al., 2020; Qi et al., 2020), we interpret the classified features as rafts of pumiceous material - pumice rafts (Supplementary Text S16 for discussion regarding difference between “true” pumice and more general pumiceous material).

Considering our interpretation that the detected rafts in Rabaul are secondary rafts, an important question to consider is what potential physical mechanisms are responsible for the mobilization of pumiceous material. One possibility is that resuspension is a consequence of local climatological conditions, e.g., high rainfall events, high wind conditions that dislodge pumice along coastlines and riverbanks back into the water. Local weather can lead to landslides and dislodgement of small pumice rafts [e.g., local pumice raft from rockslide in the Askja caldera, Iceland on 21 July 2014 (Icelandic Meteorological, 2014)]. Using ERA5 Daily Aggregate Reanalysis Product (Hersbach et al., 2020) (directly accessible through GEE), we generated time series of various atmospheric properties—daily mean air temperature, wind magnitude, wind direction, and precipitation. These time series were all sampled from the same location, directly on top of one of Rabaul’s vents, and the time series spanned the entire Sentinel-2 coverage period in the area. We did not observe any significant correlation between the daily mean air temperature and the detection of pumice rafts in the area (Supplementary Figure S3). We also explored potential correlations with weather parameters up to 10 days before raft detection to allow for some unknown advection time (Supplementary Data). Overall, we did not find significantly different results across these windows. The main statistically robust relationships in our dataset are between raft detection and wind and precipitation.

Although our ML algorithm is reasonably successful for pumice raft detection, it is not fully automated. The classification process requires manual checks to filter out incorrect classifications of pumice and cloud cover. In particular, the light cloud cover with a flat spectral response curve can at times be misclassified as pumice (and vice versa). Also, the satellite’s viewing geometry may create a “Sun glint” in certain images, where all of the pixels in the RGB rendering are affected and off-colored. The classifier subsequently has difficulty correctly identifying the correct class of each pixel. There are some ways these issues can be addressed. Better atmospheric corrected products, specifically for oceanic regions, would help improve detection. For instance, in some cases, using the atmosphere corrected Surface Reflectance (Level-2A) product can allow us to detect pumice rafts on images discolored due to atmospheric effects (Supplementary Figure S10). Alternatively, more stringent data filtering for satellite viewing angle and cloudiness bounds can help reduce potential false positives. Finally, incorporating additional satellite imagery data e.g., geostationary imagery with high temporal frequency (more detailed raft tracking as well as better cloud detection based on motion) and radar imagery (sensitive to surface roughness and reduced sensitivity to atmospheric effects) could help improve detection accuracy. Additional potential options for algorithmic improvement are described in Supplementary Text S11.