Research ArticleGEOLOGY

Insights into the 9 December 2019 eruption of Whakaari/White Island from analysis of TROPOMI SO2 imagery

See allHide authors and affiliations

Science Advances  18 Jun 2021:
Vol. 7, no. 25, eabg1218
DOI: 10.1126/sciadv.abg1218

Abstract

Small, phreatic explosions from volcanic hydrothermal systems pose a substantial proximal hazard on volcanoes, which can be popular tourist sites, creating casualty risks in case of eruption. Volcano monitoring of gas emissions provides insights into when explosions are likely to happen and unravel processes driving eruptions. Here, we report SO2 flux and plume height data retrieved from TROPOMI satellite imagery before, during, and after the 9 December 2019 eruption of Whakaari/White Island volcano, New Zealand, which resulted in 22 fatalities and numerous injuries. We show that SO2 was detected without explosive activity on separate days before and after the explosion, and that fluxes increased from 10 to 45 kg/s ~40 min before the explosion itself. High temporal resolution gas monitoring from space can provide key insights into magmatic degassing processes globally, aiding understanding of eruption precursors and complementing ground-based monitoring.

INTRODUCTION

This paper presents SO2 flux and plume height time-series data measured before, during, and after the 9 December eruption of Whakaari/White Island volcano, New Zealand (110). Data were retrieved from SO2 imagery measured with Tropospheric Monitoring Instrument (TROPOMI) on board the polar-orbiting Sentinel-5P satellite (11) using an updated and improved back-trajectory analysis from that used by Pardini et al. (12, 13) and Queißer et al. (14).

Hydrothermal systems sitting atop volcanoes are supplied with heat and volatiles from underlying magmatic systems. Activity at the surface is usually expressed as relatively benign fumarolic gas emissions and elevated temperatures. However, there is a fine balance controlling the style of activity between the overpressure in the hydrothermal/magmatic system at depth driving gas to the surface and the permeability of the shallow crust above the hydrothermal system (15). In steady-state quiescent degassing, gas overpressure is released through permeable gas pathways to the surface. However, a phreatic explosion can occur if the overpressure in the hydrothermal system increases and exceeds the sum of the confining pressure and rock strength. This may happen due to magmatic fluid transport increasing gas or heat flux (1618) or due to a decrease in permeability between the hydrothermal system and the surface (19) arising, for example, from hydrothermal mineralization (20). Such explosions can be localized and small, but still produce a major risk if people are nearby (21). If magma intrudes into the hydrothermal system, larger, more sustained explosions can occur. There is a continuum of activity between phreatic and phreatomagmatic, and it can be challenging to distinguish between the two processes, leading to diverse interpretations (22). Studying the pattern of gas emission rate may yield insights into whether an explosion is phreatic or phreatomagmatic and help determine the future style of activity.

A recent example of a phreatic eruption is the 2014 explosion of Ontake, Japan, which killed 63 people (23). Ontake, a popular hiking and tourist destination, produced its first recorded phreatic eruption in 1979, followed by small explosions in 1991 and 2007 (24). Seismic data recorded the ascent of pressurized fluid to the surface in the minutes before the eruption (23). The explosion was driven by pressurized gas emissions from a shallow hydrothermal system (25). The heat and volatiles were sourced from magma emplaced 3 km below the summit some years before (26). While there were some volcano-tectonic swarms in the days before the eruption, these did not show a clear acceleration before the explosion (23), and therefore, a deterministic forecast of the explosion was impossible (27). Seismic data recording the ascent of pressurized fluids in the minutes before the eruption (23) suggest that gas emission rates may have increased before the explosion if the hydrothermal system were not completely sealed. The continuous fumarole degassing observed at both Ontake and Whakaari/White Island suggests a partially sealed/open system. Similar precursory gas-driven, very long–period seismic activity was observed before phreatic eruptions at Ijen, Indonesia and Whakaari/White Island, both in 2013 (28). Gas emission rate monitoring might therefore inform eruption forecasting in general in the minutes before eruption if there is any leak from the overpressured system. Detailed studies of the natural variability in emission rates are essential to understand if phreatic eruption precursors are distinguishable from the background.

While the main components of volcanic gases are usually H2O and CO2, gas flux monitoring has focused on measurements of SO2 flux since the 1970s (29). This is due to the relative ease of data acquisition, as scattered sunlight is used to measure the strong ultraviolet (UV) absorption of SO2, which is readily detectable at wavelengths of between 305 and 320 nm, and because SO2 has very low background atmospheric concentrations, unlike CO2 and H2O. Ground-based monitoring is now done routinely at approximately 50 volcanoes using networks of scanning spectrometers (30), but many more volcanoes produce an SO2 plume than those that are monitored with networks (31). Polar-orbiting satellites offer global coverage, but their sensitivity to SO2 has previously been insufficient to routinely measure quiescent degassing (32). Carn et al. (31) provided a global inventory of volcanic emissions from annually averaged SO2 fluxes from data acquired by the Ozone Monitoring Instrument (OMI), gaining sensitivity at the cost of temporal frequency. The launch of TROPOMI aboard Sentinel-5P provided a step change in sensitivity and resolution compared with OMI and is currently the most sensitive daily global satellite observation platform for SO2 emission monitoring (33). Intraday SO2 flux time series quantification was demonstrated for Calbuco, Chile (12), Etna, Italy (14), and Fuego, Guatemala (13), and this paper builds on and extends that work to demonstrate the capabilities of TROPOMI for an explosion with modest gas emission.

Whakaari/White Island volcano (37.52°S, 177.18°E) is an active basaltic-andesitic composite stratovolcano well known for its high-temperature fumarole degassing (34, 35). It is located 48 km to the northeast (NE) of New Zealand’s North Island, in the Bay of Plenty (Fig. 1), at the northern end of the Taupo Volcanic Zone (TVZ) (36). The TVZ is the manifestation of volcanic activity arising from convergence between the Pacific and Australian plates, where the Pacific plate is subducting beneath the Australian plate at a rate of ~50 mm/year (34). Whakaari/White Island has an area of 3.3 km2 and a maximum height of 321 m above sea level. The crater floor is made up of variably consolidated volcanic ash deposits, with permeabilities ranging from ~10−19 to ~10−11 m2, which is typical of a material with 40 to 50% porosity (37). Multiple hot springs and fumaroles in the main crater release high-temperature (up to ~700°C) gases, producing a sustained plume. Regular airborne plume measurements of SO2, H2S, and CO2 flux have been conducted since 2003 using Correlation Spectrometer (COSPEC) and Licor sensors (20, 38). Miller et al. (39) reported approximately monthly airborne SO2 fluxes from March 2018 to November 2019, with emissions ranging from 2.5 to 22 kg/s, mean flux of 7.4 kg/s, and SD of 4.6 kg/s.

Fig. 1 Location and photograph of Whakaari/White Island.

(A) Geographical location (denoted by red triangle) and (B) photograph of Whakaari/White Island viewed from the southwest. Photo credit: Stephen James O’Meara.

Multiple phreatic, phreatomagmatic, and magmatic eruptions occurred from 1970 to 2000 (36, 40). After a little more than a decade of reduced activity from 2000, a period of higher activity resulted in phreatic and phreatomagmatic eruptions in 2012, 2013, and 2016 (27, 4143) and a small dome extrusion in 2012. The last eruptive activity before the December 2019 eruption consisted of a series of six small, discrete, phreatic explosions over a 35-min interval on 27 April 2016 (41). Because these explosions occurred at night, there were no visitors on the island and therefore no casualties. However, eruption deposits covered the crater floor, so if this had occurred during the day when visitors were present, multiple casualties may have resulted. An increase in steam-driven, geyser-like activity within the lake was noted in late September 2019 (1), similar to that seen in the 2012–2013 unrest and eruption episode. Increases in the SO2 flux and the seismic tremor over several weeks from mid-2019 led GeoNet, the geological hazard monitoring agency in New Zealand, to raise the Volcanic Alert Level from 1 to 2 (on a scale of 0 to 5) and the Aviation Color Code to Yellow on 18 November 2019 (2). This heightened tremor and flux continued up to the eruption on 9 December 2019.

At 14:11 NZDT (New Zealand Daylight Time), 01:11 UTC (all times reported in UTC hereafter), on 9 December 2019, a short-lived eruption began, lasting only 1 to 2 min. In response, GeoNet increased the Volcanic Alert Level to 4 and the Aviation Color Code to Orange (3). The eruption produced a gas and ash plume that was observed visually rising an estimated 4 km above the vent. At the time of the explosion, 47 people were present on the island and the eruption tragically resulted in 22 fatalities and numerous severe injuries. In the hours after the eruption, the seismic and eruptive activity level decreased, with GeoNet lowering the Volcanic Alert Level to 3 (4). Search-and-Rescue teams noted strong sulfur smells and irritation to their eyes and any exposed skin (44). Seismic tremor increased again on 11 December, accompanied by steam venting and mud jetting at the active vent sites (5). The Volcanic Alert Level was decreased to 2 on 12 December (6), in response to the lack of continuing eruptive activity. Small-scale gas jetting was observed on 13 December (7), and vigorous gas emissions were observed on 15 December (8). During post-eruption overflights on 12 and 19 December 2019, SO2 gas fluxes of ~20 kg/s and ~15 kg/s were measured, respectively (9). A multiple-lobed lava dome was observed in January 2020, but its time of emplacement is uncertain. Throughout 2020, vent temperatures continued to decline and the Volcanic Alert Level was lowered to 1 on 16 June 2020 (10). A brief period of weak, nonexplosive ash emission saw the Volcanic Alert Level raised again to 2 on 13 November 2020 (45) before being lowered back to 1 on 7 December 2020 (46).

RESULTS

Insights into the 9 December eruption dynamics arise from analysis of a single TROPOMI Level 2 SO2 image [details of the SO2 retrieval are reported by Theys et al. (11)], collected at 02:09 UTC (15:09 NZDT), fortuitously only 58 min after the eruption onset at 01:11 UTC (Fig. 2A). After removing high-noise SO2 pixels (those with SO2 less than three times reported SO2 error), we further refined the data using a nearest-neighbor approach, which required that each pixel must have at least two other adjacent low-error SO2 pixels, removing spurious individual pixels not associated with a plume (Fig. 2B). This may result in an underestimation of SO2 emissions but greatly decreases the chance of noise being erroneously included in our analysis. We calculate backward trajectories using Hybrid Single Particle Lagrangian Integrated Trajectory (HYSPLIT) (47) and three-dimensional (3D) wind fields from the National Oceanic and Atmospheric Administration (NOAA) National Centers for Environmental Prediction (NCEP) Global Forecast System (GFS) to determine the age and height of each pixel in the image (see Materials and Methods) and correct the height-sensitive vertical column density (VCD) of SO2 in each pixel for the appropriate retrieved height (Fig. 2C). VCDs are reported in Dobson units (DUs; 1 DU = 2.69 × 1016 molecules cm−2). The sensitivity and spatial resolution of TROPOMI is remarkable, with a detection limit of ~2 DU in this image at 1- to 2-km altitude (sensitivity is strongly dependent on altitude and increases with increasing plume altitude), and pixel size of 3.6 km × 5.6 km (for details on TROPOMI and operations, see https://sentinel.esa.int/web/sentinel/user-guides/sentinel-5p-tropomi). The data quality is helped by the fact that the image was collected in the middle of a cloud-free day during the Southern Hemisphere summer, when solar insolation is maximum, and that the plume was close to the swath center. Significantly lower resolution and sensitivity would be achieved at the edge of the swath (swath width is 2700 km and orbit height is 824 km, so pixel size at the swath edge is ~13 km) and in midwinter at high latitudes.

Fig. 2 Results from the analysis of the TROPOMI image collected on 9 December 2019 at 02:09 UTC, 58 min after the Whakaari/White Island eruption.

(A) Unfiltered L2 TROPOMI SO2 1-km VCD data. (B) SO2 VCD, with the VCD more than three times the reported random error, and with two neighboring pixels also above this threshold. (C) SO2 VCD, corrected to the retrieved plume altitude for that pixel; back trajectories for each successful pixel, colored by (D) pixel age and (E) pixel altitude; pixel maps colored by (F) closest approach distance for the successful trajectory, (G) pixel age, (H) pixel altitude, and (I) instantaneous pixel SO2 flux. The eruption time is denoted by a white line on the color bar of (D) and (G).

Retrieved plume back trajectories for each pixel, selected by closest approach to the volcanic source as a function of height (Fig. 2F; see Materials and Methods), are shown in Fig. 2 (D and E), colored by pixel age and pixel height, respectively. Maps of pixel age and height are shown in Fig. 2 (G and H). Our results show that the oldest part of the plume propagated at an altitude of 0.5 to 1.0 km toward the NE, and the maximum observed pixel age was 2.5 hours, containing gas emitted 1.5 hours before the eruption. The youngest part of the plume is attached to the volcano and is also at an altitude of ~1 km and was emitted approximately 1 hour after eruption. We therefore capture in a single image the gas emissions produced by Whakaari/White Island volcano in the period from 1.5 hours before the eruption to 1 hour after eruption. The highest plume emission reached 5 km above sea level and propagated to the southeast, reaching the coastline at the moment of image capture. The greater distance traveled by this SE lobe compared to the older plume region indicates that the wind velocity at 5 km was greater than that at 1 to 2 km. Together with the diverse wind directions, this reflects the wind shear in the atmosphere at the time (see Materials and Methods). Information on the wind conditions is critical in this analysis and was provided by NOAA NCEP GFS 0.25-degree meteorological data (ftp://arlftp.arlhq.noaa.gov/pub/archives/gfs0p25/).

In Fig. 2I, we report the map of pixel SO2 flux, calculated by dividing the mass of SO2 in each pixel (which is the product of VCD and pixel area) with the time taken for wind to advect across the pixel. The time associated with the pixel flux is the time of emission. This is a new data product that has not previously been calculated but, as shown below, is effective in revealing the dynamics of the eruption. We see in Fig. 2I that the low-altitude plume propagating NE has SO2 fluxes of 2 to 4 kg/s, while the higher altitude plume propagating SE has a peak pixel SO2 flux of 8 to 10 kg/s.

Temporal evolution of plume height, SO2 mass, and SO2 flux

We define the atmospheric injection time for each pixel as the time at which its back trajectory reached minimum distance to the volcano. We know the time from which to calculate the back trajectory, as this is the acquisition time for each pixel and is reported in the TROPOMI SO2 product. As we see in Fig. 3A, the distance of closest approach to the volcano varies between 0 and 7 km, i.e., approximately the dimension of one pixel. During an explosive eruption, a degree of plume spreading takes place due to a range of initial ballistic trajectory angles and due to spreading at neutral buoyancy, which is strongly dependent on mass injection rate (48). For these reasons, the trajectories do not all converge at zero distance from the volcano. In other words, the eruption column is not a point source, but is instead a somewhat distributed one. However, keeping distances of closest approach within the size of a pixel is an indication of good performance. In a strongly explosive eruption forming an umbrella cloud, this distance can become >50 km. We must also consider that the time of closest approach is the time at which the gas began to advect with the wind field and therefore does not include the transport time from vent to altitude where advection begins. Lastly, the HYSPLIT model has a fixed 10-min time step, introducing a 10-min temporal resolution to the back trajectories.

Fig. 3 Results from analysis of the TROPOMI image collected on 9 December 2019 at 02:09 UTC, 58 min after the eruption of Whakaari/White Island (indicated by the red line, reflecting the observed explosion duration).

All data are shown versus the pixel injection time. (A) Bottom, pixel injection altitude, colored by SO2 mass; top, mean pixel altitude for pixels within the 10-min bin, colored by total SO2 mass of those pixels. Black diamonds show peak pixel altitude within the bin. (B) Bottom, pixel SO2 mass, colored by pixel injection altitude; top, total SO2 mass of pixels within the 10-min bin, colored by the mean pixel injection altitude. Binned flux (calculated using the 10-min bin time) is scaled on the right axis. (C) Bottom, instantaneous SO2 pixel flux, colored by the pixel’s injection altitude; top, average transit time through pixels within the 10-min bin, colored by the mean pixel injection altitude.

Using the time of closest approach to the volcano for each pixel, we present the plume pixel height as a function of injection time in Fig. 3A. We have calculated the error on time and height, and these are shown with error bars. The time error is determined by calculating back trajectories for the vertices of the pixel as well as the center to determine the minimum and maximum time by which the gas could advect from the volcano to the pixel. The height error is determined by the minimum and maximum heights from which the edge of the pixel can be reached from the location of closest approach (see Materials and Methods for full error analysis). The plume injection height was stable at <1 km until 00:00, when we see three pixels just under 2 km, and then from 00:40 the injection height reached up to 2.5 km. The time and duration of the eruption is shown with a gray band at 01:11 to 01:15, and we see a single pixel at 3-km altitude in the time between 01:10 and 01:20, and three pixels close to 5-km altitude in the following two time bins. Injections at lower altitudes continued throughout the time up to image acquisition time, in two clusters at 1 and 2 km, respectively. In the upper part of Fig. 3A, we report these heights in histogram form, showing the mean and SD of the measured heights, with diamonds representing peak heights. We observe an onset of plume emission at 1 km at 23:40, followed by a step increase to 2 km at 00:40, with a peak in altitude in the 10-min bin after the eruption occurred. This lag between plume height peak time and explosion time reflects the transport time from vent to altitude where wind advection begins. This lag time is 5 to 10 min, and the maximum height is 6 km, reflecting an average ascent speed of 10 to 20 m/s. The color of each bin in the histogram and the color of pixels reflect the mass of SO2 in metric tons and show the increase in SO2 emission around the explosion time.

Figure 3B shows the time series of pixel SO2 mass, with each point colored by injection altitude. Errors on SO2 mass reflect the random error reported by the TROPOMI retrieval and the height error, as the mass calculation is height dependent. We observe a wide range of masses that increase in variability during the sequence. A source of this variability is the geometry of each pixel and distance from plume center, and the fact that different injection heights have different wind speeds. The upper section on Fig. 3C shows a histogram of integrated mass and colored by average altitude, and this shows a clear increasing trend up to the explosion and decreasing after. This histogram also shows with the right-hand scale the integrated SO2 flux in kg/s for each bin, calculated by dividing the integrated mass of each bin by the duration of each bin (10 min).

Figure 3C shows the time series pixel flux, which follows a clear increasing trend before the explosion followed by a decline. This plot appears to be more physically representative of the eruption dynamics than that of individual SO2 pixel masses shown in Fig. 3, reflecting the key additional information contained in the wind speed. This is reflected in the time taken to traverse the pixel, termed the pixel duration, and this is shown in histogram format at the top of Fig. 3C. The trend here is affected by higher wind speed at higher altitude, resulting in shorter pixel duration during the explosion. The color of each histogram bar reflects the average height of pixels in the time bin. The typical time error for each pixel is ±10 min, meaning that there is some mixing of SO2 between time bins producing a smoothing of the gas emission signal, but the precursory increases in height and flux start ~40 min before the eruption, indicating that the spread in flux is not wholly attributable to averaging but instead reflects an observable precursory increase.

Overall, we observed a precursory increase in SO2 flux and plume height from 00:40, ~30 min before the explosion. The average SO2 flux for the whole image is calculated as the integrated mass of all pixels (165 [+11/−18] metric tons) divided by the age of the oldest pixel (161 [+0/−40] min), and is 17.1 [+1.1/−1.8] kg/s. The peak SO2 flux observed is 49 [+2/−6] kg/s.

In Fig. 4A, we report the SO2 flux for each time bin and the seismic amplitude recorded from a seismic station on Whakaari/White Island, bandpass-filtered between 2 and 5 Hz, and averaged into 10-min bins (raw data for the entire time series are shown in fig. S4). The tremor signal also demonstrates an increased amplitude compared with background before the explosion, in agreement with our pre-explosion elevated flux observations, and the post-eruption trend of decreasing flux is also reflected in the tremor signal. The actual peak in SO2 flux is the time bin before the explosion, suggesting perhaps that the high aerosol content of the explosion produced a degree of underestimation in SO2 columns. The overall agreement between seismic amplitude and SO2 flux is consistent with similar relationships observed on other volcanoes [e.g., (49)]. The temporal coincidence of the tremor and gas signals suggests that the wind field accuracy results in timing errors <10 min. Figure 4B shows pixel height plotted against pixel flux, showing three clusters of emissions at heights of 1-, 2-, and 5-km altitude, consistent with a background degassing, elevated pre- and post-eruptive degassing, and explosive degassing, respectively.

Fig. 4 Tremor and SO2 mass time series and scatter plot of pixel flux versus injection altitude.

(A) SO2 mass as a function of time, with Real-time Seismic Amplitude Measurement (RSAM) tremor ground velocity in nanometers per second from a seismic station on Whakaari/White Island. (B) Pixel injection altitude against SO2 flux. Three clusters are defined by altitude, for the 1-km altitude initial phase, 2-km elevated pre- and post-explosion phases, and the 5-km explosion phase.

To evaluate whether SO2 degassing observed on 9 December 2019 was exceptional or typical, we analyzed TROPOMI data for 1 month before and after the explosion. The time series of average fluxes and plume heights for each image are shown in Fig. 5, and the corrected VCD maps are shown in Fig. 6. We detected an SO2 plume on 3 days of the 30 days before 9 December on 26, 29, and 30 November, and on 7 days of the 31 days after explosion. The average fluxes seen after eruption (3.62 ± 1.58 kg/s) are higher than those seen before eruption (2.02 ± 2.14 kg/s). The pre-eruption plume heights were higher (2.59 ± 1.07 km) than the post-eruption period (1.29 ± 1.58 km). The plume from 30 November (fig. S2) was the largest noneruptive event in the series, with an SO2 mass of 97 [+5/−8] metric tons. The 10-min flux for that day did not exceed 4 kg/s (fig. S3B)—the flux exceeded this value 30 min before the eruption on 9 December, although only regularly in the 20 min before. The injection altitude, however, was the lowest of the pre-eruptive events, not exceeding 2 km (fig. S3A).

Fig. 5 Time series of average to daily SO2 flux and the average injection altitude for 30 days before and after the 9 December eruption (eruption denoted by red line) from this work and COSPEC measurements.

The highest flux of the period is recorded on the day of the eruption (9 December 2019). The flux is calculated using the altitude-corrected SO2 masses for TROPOMI pixels that pass the quality control divided by the age of the oldest observed pixel. The error bars reflect the mass uncertainty but do not include the systematic error arising from gas pixels below the detection limit; thus, they are a lower limit.

Fig. 6 Corrected SO2 VCD plots for all detected plumes within the study period (9 November 2019 to 9 January 2020).

The observation of multiple days with detectable SO2, high fluxes, and plume heights without explosions suggests that sporadic intense degassing occurs on Whakaari/White Island without leading to eruption but may also be influenced by meteorological conditions that occasionally favor a higher advection of the plume. We also see spikes in tremor unassociated with detectable SO2 emissions (fig. S4). The maximum age of each observed plume is less than 12 hours, so chemical processing of SO2 to H2SO4 is unlikely to be a significant factor in our observations (50).

DISCUSSION

Measuring gas emissions from small volcanic explosive events is rarely achieved, because the focus on the volcanoes where gas flux is regularly monitored is usually quiescent degassing, meaning the high gas/ash fluxes that typify explosions make proximal measurements extremely challenging and strongly affected by light dilution (51). In our case, the TROPOMI image was collected 58 min after eruption, allowing solid ash sedimentation and aerosol evaporation to occur, reducing the impact of aerosols, but perhaps not completely excluding the impact. This work demonstrates that we can now analyze the dynamics of small eruptions at high temporal resolution and gas mass accuracy through the combination of PlumeTraj and TROPOMI data, opening up a new frontier in volcanological research. Globally, the frequency of explosive eruptions increases exponentially with decreasing magnitude (52), so this new capacity to measure gas fluxes and plume heights from small/modest eruptions greatly increases the number of eruptions for which quantitative SO2 flux and plume height time series can be made.

Our results allow us to reflect on the processes that triggered the 9 December 2019 eruption of Whakaari/White Island, and this volcano’s gas flux. Christenson et al. (20) report that Whakaari/White Island exhibits pulsatory degassing, and we detected only a few days where emissions were visible with TROPOMI. The sensitivity of TROPOMI is strongly dependent on plume height, and processes that increase the plume height will increase the probability of detection. Such processes include higher gas emission rates, elevated gas temperatures, low wind speeds, and a strong adiabatic lapse rate. The observation of SO2 emissions with a relatively high SO2 flux and integrated mass in the absence of eruptive activity on 30 November is consistent with previous observations of pulsatory gas fluxes on Whakaari/White Island, and that the transition from intense degassing to eruption is challenging to forecast accurately from gas flux alone. One implication is that during an intense degassing event, the system may be close to producing a phreatic explosion, and a trivial change in permeability/overpressure may trigger the eruption (a “leaky pipe”). This would explain the continued increase in the SO2 flux before the 9 December event and why the 30 November event did not trigger an explosion—the gas continued to vent, but on 9 December, the system was simply not able to expel the gas quickly enough to reduce the overpressure, leading to the explosion. A further distinction is that the plume height was significantly lower on 30 November, consistent with a lower gas flux and potentially lower temperature of emission. The maximum flux and plume height observed on 30 November was 4 kg−1 and 2 km, respectively (fig. S3), suggesting that these might be threshold values beyond which an eruption becomes more likely.

The pattern of degassing observed in our dataset, where the gas flux decreases rapidly after the explosion, is consistent with a shallow subsurface gas overpressure. This is also consistent with increased pressure and fluid migration revealed by geophysical observations before explosions at White Island (18) and increased gas flux before phreatic activity at Poas in Costa Rica (16, 17). There are several processes that may have produced the subsurface overpressure on 9 December 2019: a decrease in the permeability of the degassing pathway (53, 54), possibly due to external processes such as meteorological or tidal effects (38); an increase in fluid supply from deeper in the magmatic plumbing system; and a shallow magma intrusion that cooled on emplacement and degassed, producing an additional gas source. The extrusion of a small lava dome days to weeks after the eruption suggests that a mostly degassed magma batch was present at shallow levels in the conduit. The intrusion of magma in the weeks to months prior likely drove the initial unrest including elevated tremor, gas flux, and geysering activity (39). Processes that can explain the observed increase in gas flux ~40 min before the explosion include gradual failure of a seal holding pressurized gas or that a gas pulse ripped through the system leaking out through whatever channels available but ultimately led to failure of the critically stressed seal. We can, however, exclude a primarily magma ascent–driven process; in this case, an increasing gas flux following eruption initiation would be expected as the magmatic system depressurized and degassed, and this was not observed.

In terms of eruption forecasting, our observations of an increase in SO2 flux, plume height, and tremor ~40 min before the explosion of 9 December suggest that a significant increase in heat flux may have occurred, providing a potential precursory signal that indirectly reflects on gas flux. Coupling of gas flows with wall-rock can produce volcanic tremor, exemplified by observed correlations between degassing rate and tremor intensity (49), and this appears to be the case for this explosion. The degassing event on 30 November produced no such correlation though (fig. S5). TROPOMI observations of Whakaari/White Island can provide insights into quiescent degassing, but plume detection is also influenced strongly by the meteorological conditions, which may promote higher plume altitudes, playing a key role in allowing a plume to be detected. Once detected, we can reconstruct a period of emission rates whose duration increases with plume height and the magnitude and duration of the gas flux. Our results are therefore complementary to ground-based gas flux networks, because for such weak plumes routine daily monitoring is not readily achievable from space, but the satellite observations are useful when stronger degassing or explosions occur and higher plumes form, which may produce challenges for quantification from ground-based systems. In these cases, TROPOMI and PlumeTraj provide insights into degassing mechanisms, and data can be produced in near real time, 3 to 4 hours after overpass. A plume containing just 165 [+11/−18] metric tons of SO2 is one of the smallest volcanic gas masses to be quantified from space, highlighting the sensitivity of TROPOMI.

Our work demonstrates that high–temporal resolution observations of SO2 flux are possible even from very modestly degassing volcanoes. That magmatic degassing from more volcanoes than ever before could be observed regularly from space heralds a new frontier in observational volcanology. When combined with geophysical observations of deformation, tremor, and thermal emission, these new data will provide a new window into the processes controlling volcanic activity and potentially aiding in the detection of eruption precursors.

MATERIALS AND METHODS

ESA’s (European Space Agency) Sentinel constellation is an integrated series of satellite platforms and instruments designed to monitor the Earth system as a whole. Launched in October 2017, Sentinel-5P orbits with a local equatorial overpass time of 13:30 LT (ascending). It carries TROPOMI, a hyperspectral instrument composed of four spectrometers, including one UV spectrometer, in push-broom imaging configuration. It has a spatial resolution of 3.6 km × 5.6 km (nadir, along-track × across-track; 7.2 km × 5.6 km before August 2019) and a swath width of 2600 km, providing daily global coverage (55).

This work uses the Level 2 Offline SO2 dataset [S5P_OFFL_L2__SO2__, version 1.01.07 from the ESA-Copernicus Sentinel-5P Pre-Operation Data-Hub (https://s5phub.copernicus.eu/)], providing SO2 total VCDs in DU. Averaging kernels (AKs) are used to convert the slant column concentrations measured by the instrument into the VCDs used in the analysis. The AKs are produced for a modeled atmosphere, with SO2 located in one of three altitude bands, each 1 km thick and centered on 1, 7, or 15 km (33).

To calculate an emission flux from a single satellite image of a volcanic plume, the age (time since emission) and mass of SO2 within each pixel must be determined. Using the 1-km L2 TROPOMI dataset, plume pixels are selected by passing two criteria. First, they must have an SO2 VCD higher than three times the random error; second, there must be at least two other pixels in the surrounding eight pixels that also pass this threshold, reflecting the connected plume morphology. High ash or aerosol loading in a volcanic plume can lead to high error values, meaning pixels in the core of the plume are sometimes excluded. The retrieved SO2 VCD is dependent on the height of the plume, and this is generally a major source of uncertainty in SO2 quantification from space. We developed an improved version, focusing on quantitative error budgets, of the PlumeTraj method, adapted from the method proposed by Pardini et al. (56) and Queißer et al. (14), to calculate the measurement and plume altitudes as well as the emission time.

For each selected pixel, backward trajectories are calculated using the HYSPLIT dispersal model (47). The NOAA NCEP GFS 0.25-degree meteorological dataset (www.ready.noaa.gov/archives.php) is used as the 3D wind field for the trajectory calculations. The center of each flagged pixel is used as the initiation location, and the back trajectories are performed over a user-selected range of altitudes. Natural wind shear in the atmosphere means that air at different altitudes moves at different speeds and in different directions (Fig. 7, A, C, and D). The point of closest approach to the volcano for each trajectory is calculated and used to find the altitude that generates the minimum approach distance (Fig. 7B). This method of trajectory selection is the major difference from the original method (56) and was found to be more robust. The altitude uncertainty is calculated using the range in modeled initiation altitudes that would produce back trajectories that pass within half the pixel’s width of the minimum approach distance to the volcano. For example, for a measurement made at nadir, with a pixel width of 3.6 km, the altitude uncertainty would correspond to back trajectories with an approach distance of ±1.8 km from the minimum approach distance. To be classified as a successful trajectory, the minimum approach distance must be less than the threshold distance of 25 km (approximately the size of one 0.25° wind data grid box) to account for discrepancies within the wind field.

Fig. 7 Trajectory analysis for one pixel from TROPOMI orbit 11161 (9 December 2019).

(A) Trajectories for all altitudes, from 500 m to 7.5 km. (B) Trajectory initiation altitude versus distance of closest approach. The minimum distance is the red line, giving an altitude of 1.66 km. The black dashed line represents minimum distance + 0.5 × pixel width, giving the uncertainty on the altitude [−0.21, +0.19 km]. (C and D) Vertical profile of (C) wind speed and (D) wind direction at the moment of explosion on 9 December 2019. Note that vertical wind shear in both speed and height helps to identify the correct plume height.

For each pixel, trajectories are run again using the altitudes corresponding to the minimum approach distance and the two half-pixel width distances. These trajectories give the time of flight from the point of measurement to the point of closest approach; subtracting this from the measurement time gives the pixel’s emission time. The SO2 concentration for each pixel is linearly interpolated to the minimum approach altitude from the three L2 dataset altitudes (1, 7, and 15 km) and then multiplied by the pixel area to give the SO2 mass. Binned SO2 flux is calculated by combining the SO2 mass and the emission time for each pixel, grouping all arrivals within a specified time interval. Pixel SO2 flux is calculated by dividing the SO2 mass of each pixel with the time taken for air to travel through the pixel using the pixel shape, wind speed, and wind direction at the retrieved height of that pixel. This time interval is limited by the time step of the trajectory model (how often the trajectory location is calculated). For HYSPLIT, this is 10 min, meaning it is possible to produce a 10-min flux or any higher multiple thereof. The calculations are performed on all three trajectories (the optimal altitude and the half-pixel width error altitudes), giving the error bounds on each of the values derived from the optimal altitude.

The method is dependent on the vertical stratification of the wind field (Fig. 7, C and D). When the wind shear is weak, trajectories from numerous altitudes can return to close to the volcano, meaning there is a greater uncertainty in the altitude and, therefore, the SO2 mass. On days with strong wind shear, a small change in the altitude can result in a significant difference in the speed and direction of the trajectory, leading to a greater uncertainty in the emission time and emission rate.

TROPOMI data were downloaded for the period 1 month either side of the eruption (9 November 2019 to 9 January 2020). To identify orbits containing a plume from Whakaari/White Island, the raw 1-km SO2 VCD data for each orbit were plotted and inspected for a distinct plume. Plume structure and absolute concentration were manually examined to ensure that orbits with only noise pixels were excluded (see fig. S1 for 1-km raw VCD plots for all orbits within the study period). Once all orbits containing a plume were identified, PlumeTraj, as described above, was initiated on altitudes between 0.1 and 7.5 km (every 100 m between 0.1 and 1.0 km, then every 250 m from 1.0 to 7.5 km). The SO2 column density for each pixel was then calculated using linear interpolation between the three altitude bands (1, 7, and 15 km) to the altitude returned by the trajectory analysis (Fig. 6). The individual orbit fluxes were then combined to produce a time series (Fig. 5).

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/7/25/eabg1218/DC1

https://creativecommons.org/licenses/by-nc/4.0/

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 S. Sherburn, the New Zealand GeoNet project, and its sponsors EQC, GNS Science, and INZ for providing the seismic data used in this study. We gratefully acknowledge constructive reviews from E. Liu and an anonymous reviewer. Early reviews of this manuscript were provided by G. Jolly and S. Sherburn. Funding: We acknowledge the NERC-funded V-PLUS project NE/S004106/1 and COMET, the NERC Centre for the Observation and Modelling of Earthquakes, Volcanoes and Tectonics, a partnership between UK Universities and the British Geological Survey (https://comet.nerc.ac.uk/). The research leading to these results has received funding from the European Union Horizon 2020 project EUROVOLC (agreement no. 731070). C.M. and B.C. are supported by the MBIE-funded GNS Science Volcanic Hazards Programme. Author contributions: M.B. and C.H. conceived this study. C.H. analyzed data and produced figures. M.B. and C.H. led manuscript production, with all authors contributing to interpretation, writing, and figure preparation. 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. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article