Modeling SO₂ dispersion from future eruptions in the Auckland Volcanic Field, New Zealand

Abstract

Auckland city (pop. 1.7 M) is Aotearoa New Zealand’s largest city and an important economic hub. The city is built upon the active intraplate basaltic Auckland Volcanic Field (AVF). An AVF eruption would cause considerable impacts. An important component of volcanic risk management is assessing the likely volcanic hazards to help inform emergency planning and other preparedness activities. Previous volcanic hazard assessments for the AVF, particularly those for emergency planning scenarios, have modeled multiple volcanic hazards including lava flows, pyroclastic density currents, ballistic projectiles and tephra fall. Despite volcanic gas being an important and impactful hazard from intraplate basaltic field eruptions, there has been limited consideration of volcanic gas in AVF hazard assessment to date. This project is one of the first to quantitatively assess potential volcanic gas hazards for an explosive eruption scenario. For basaltic volcanism, sulfur dioxide (SO₂) gas is typically the most consequential volcanic gas emitted. The aim of this exploratory study was to model SO₂ dispersion from a high impact eruption during weather conditions conducive to high ground level pollutant concentrations. Since ground level SO₂ concentrations are influenced by complex wind patterns resulting from interactions of locally driven flow circulations and topographically influenced weather, we modeled SO₂ dispersion using the HYSPLIT model, a state-of-the art hybrid Eulerian and Lagrangian dispersion model widely used for volcanic gases, using high-resolution meteorological forcing fields given by the Weather Research and Forecasting (WRF) model.

Modeled air parcel trajectories and ground level SO₂ concentrations illustrate the effect of the converging sea breeze winds on SO₂ dispersion. Under worst-case dispersion conditions, extensive areas of up to hundreds of square kilometers to the north and northwest of the eruption location would exceed New Zealand short-term (24 h) air quality standards and guidelines for SO₂, indicating heightened health risks to downwind communities. Using this numerical modeling approach, this work presents a methodology for future applications to other AVF eruption scenarios, with a wider range of meteorological conditions that can help in exploring consequences for health services such as anticipated emergency department respiratory admissions.

Introduction

Tāmaki Makaurau/Auckland is the most populous city in Aotearoa New Zealand, with a population of ~ 1.7 million (Fig. 1). It is an important economic hub and contributes approximately 38% of New Zealand’s total Gross Domestic Product (Stats NZ: Auckland 2021; Reid et al. 2009). The city is built on the Auckland Volcanic Field (AVF), a distributed basaltic intraplate volcanic field with 53 recognized volcanic centers (Fig. 2A) (Hopkins et al. 2020). In the last ~ 190 ka, the volcanic field has erupted over 55 times, depositing roughly 2 km³ of volcanic deposits on the surface (Deligne et al. 2017). Visible features of previous AVF eruptions include maar craters with tephra (tuff) rings, scoria cones, and lava flows (Deligne et al. 2017; Hayes et al. 2018). These visible features suggest that for future eruptions there may be multiple hazards including pyroclastic surges, ballistic projectiles, tephra fall, lava flows, earthquakes, tsunami, and fires (Deligne et al. 2017; Hayes et al. 2018; Hopkins et al. 2020). While the annual probability of a new eruption is comparatively low (Bebbington et al. 2018), any new eruption would likely be catastrophic for Auckland and have substantial national impacts for New Zealand. This has led to the establishment of the Determining Volcanic Risk in Auckland (DEVORA) program in 2008 to improve the hazard assessment and emergency planning for Auckland (Hayes et al. 2018). An important component of volcanic risk management in Auckland has been the development of eight eruption scenarios (Fig. 2B) using a multidisciplinary approach to hazard planning, which incorporates different locations, eruption types and lengths, and emergency management responses (Deligne et al. 2017; Hayes et al. 2018, 2020). The scenarios are not intended to predict the next eruption, but to support disaster risk management (e.g. emergency planning) for the many eruption possibilities. As such, the scenario suite is a compromise between fully deterministic and fully probabilistic approaches (Hayes et al. 2020).

Fig. 1

Population distribution across the Tāmaki Makaurau/Auckland region, Aotearoa New Zealand

Fig. 2

A Volcanic centers of the Auckland Volcanic Field (AVF) B Scenario locations within Auckland developed for the DEVORA program (Hayes et al. 2018)

While hazard footprints have been developed for lava flows, pyroclastic flows and surges, edifice building, ashfall, and ballistic projectiles (Hayes et al. 2020), no attempts have so far been made to quantitatively assess volcanic gas dispersion from future AVF eruptions. Gas emissions are in general an under investigated volcanic hazard (e.g., Edmonds et al. 2015). They are primarily comprised of water vapor (H₂O), carbon dioxide (CO₂) and sulfur dioxide (SO₂), with more minor quantities of hydrogen sulfide (H₂S), hydrochloric acid (HCl), and hydrofluoric acid (HF) (Hansell and Oppenheimer 2004). High emission fluxes of SO₂ are probable in a future AVF eruption due to the high sulfur content of basaltic magma (e.g., Oppenheimer et al. 2011) and the rapid magma ascent rates expected for future AVF eruptions (Brenna et al. 2018). The degassing episodes expected for the AVF are likely to be medial in scale: plumes are released into the troposphere and hence impacted by short-term meteorological variability associated with diurnal weather patterns and air flow interaction with local topography (Poulidis et al. 2018). The plumes are likely to be primarily SO₂ but may also contain minor HCl and HF (Stewart et al. 2019).

In the atmosphere, SO₂ reacts with oxygen and water vapor, in the presence of sunlight, to form sulfate aerosol particles of approximately 0.1–0.3 µm diameter (Longo et al. 2008), primarily composed of dilute sulfuric acid. The conversion rates of SO₂ to sulfate are highly variable. Businger et al. (2015) estimate SO₂ loss rates of 0.01–5% h⁻¹ in clear air and 3–50% s⁻¹ in cloudy air, and note that these rates vary diurnally, with faster removal in the daytime.

While SO₂ gas is colorless, the mixture of unreacted SO₂ gas and sulfate aerosol forms hazy air pollution known as vog (volcanic smog). Vog has been a well-documented phenomenon on the Island of Hawai’i, and beyond, since the long-term eruption of Kīlauea volcano commenced in 1983 (Pattantyus et al. 2018; Smithsonian 2023). There have been several epidemiological, clinical and descriptive studies of health impacts of vog on residents of both the Island of Hawai’i and further afield (summarized in Stewart et al. 2022). The SO₂ component of vog is greater near the sources (Nelson and Sewake 2008). SO₂ is a fast-acting respiratory irritant, as when it contacts the body’s moist mucous membranes it reacts to form the severe irritant sulfurous acid (Reno et al. 2015). Short-term exposure can irritate the eyes, nose, throat and lungs, leading to coughing, wheezing and bronchoconstriction (Stewart et al. 2022; National Research Council 2010). Asthmatics are particularly sensitive to SO₂ exposure (Reno et al. 2015). Other high-risk groups include people with respiratory or heart conditions, older adults and children (Longo et al. 2008; Longo 2013).

There is also growing interest in the international aviation sector about the impacts of volcanic SO₂ emissions on both the health of aircraft occupants and the service life of aircraft components (Schmidt et al. 2014; International Civil Aviation Organization (ICAO) 2021).

While many of the hazards associated with a future AVF eruption have been explored in depth, there is little understanding of gas emission rates and no attempts have been made to model volcanic gas dispersion for the DEVORA eruption scenarios. Given this gap in AVF hazard assessment, in this paper we present a new methodology for simulating volcanic gas dispersion and apply this model to one of the eight DEVORA scenarios, taking into account the complex interactions of locally driven flow circulations and topographically influenced weather patterns. Because of the exploratory nature of this study, we decided to model for a worst-case scenario by choosing an eruption scenario with high gas flux rates and meteorological conditions associated with minimal dispersion of air pollutants.

Methods

Choice of scenario

A single eruption scenario (Scenario C: Māngere Bridge; Fig. 2B) was selected from the DEVORA scenario suite for the following reasons: 1) this scenario is the most fully realized (Blake et al. 2017; Deligne et al. 2017); and 2) this scenario is at the upper end of eruptive volumes for the suite (0.1 km³) but has a relatively short duration (32 days), implying high gas fluxes (Hayes et al. 2018). Table 1 summarizes Scenario C eruption conditions used in this study.

Table 1 Chosen eruption parameters from DEVORA scenario suite, from Hayes et al. (2018)

Estimation of gas fluxes from petrologic method

There have been no historic accounts of eruptions in the AVF to provide direct inputs to gas dispersion models. However, as part of an ongoing study, Smid et al. (in prep) have estimated gas emissions for selected past eruptions in the AVF as well as gas emissions and daily fluxes for the DEVORA eruption scenario suite. In this approach, estimated magma volume and eruption duration from the DEVORA eruption scenarios, geochemical analyses of previously erupted volcanic material, and the petrologic method (e.g. Devine et al. 1984; Sigurdsson et al. 1985; Horn and Schmincke 2000; Sharma et al. 2004; Kereszturi et al. 2013), are used to estimate total SO₂ emissions of past eruptions. In the petrologic method, sulfur concentrations in melt inclusions and degassed groundmass glass are analyzed via electron probe microanalyzer. The sulfur concentration in the melt inclusion is assumed to represent the maximum concentration of sulfur found in the magma, while the sulfur concentration in the groundmass glass is taken as the post-eruption sulfur concentration still left in the volcanic deposits. After taking into account melt volume, density, and crystallinity, the difference between the two analyses is assumed to represent the minimum total degassed (emitted) amount of sulfur for that eruption. For the AVF, Smid et al. (in prep) assumed a melt density of 2,590 kg/m³ and a melt fraction of 0.8 (crystallinity of 0.2; Saito et al. 2005; Ferguson 2018). Dense rock equivalent erupted volumes for AVF volcanoes may be found in Kereszturi et al. (2013).

Assuming that magmatic sulfur concentrations and degassing patterns in future eruptions will be similar to those in past eruptions, we used the petrologic method to estimate total sulfur emissions for Scenario C, using the assumed bulk erupted volume (0.1 km³, Table 1). The total estimated SO₂ emission for Scenario C is 0.41 Mt (Smid et al. in prep). Daily fluxes were estimated by assuming an ideal log-normal distribution for magma effusion, wherein magma effusion (and consequently, the SO₂ degassing rate) is high at the start of eruptions and then wanes following a simple logarithmic decay curve (Wadge 1981). The total SO₂ emissions for Scenario C were distributed over its 32-day time period to simulate a Wadge-type curve (Fig. 3). This analysis yields a maximum daily flux of ~ 50,000 t (0.05 Mt) on Day 2 of the eruption (Smid et al. in prep). In comparison, a maximum daily flux of 0.03 ± 0.01 Mt was estimated for the 1973 Eldfell eruption, Iceland, with a similar eruptive volume to Scenario C (Wallace 2001; Stewart et al. 2019). The four-day window with the highest SO₂ fluxes was chosen as the modeling time period, in line with our worst-case scenario approach (Fig. 3).

Fig. 3

Estimated SO₂ flux over the active eruption period for a DEVORA scenario C, with the four days of highest flux highlighted in green

Meteorological scenario

The scenario for meteorological modeling (June 27th-30th, 2019) was selected after analyzing data from air quality stations at the AVF area (Auckland Council Environmental Data Portal 2021) as well as synoptic regimes associated with a high probability of significant air pollution in Auckland (Jiang et al. 2005, 2014; Jiang 2008; Griffiths et al. 2019). High pollution events (hourly aggregates) from PM2.5 pollution data from 2010–2020 at the Penrose and Takapuna stations (Fig. 6A) were analyzed. The high pollution threshold in our case study was the same used in the two-site analysis: the 95% percentile, or PM2.5 > 15.7 μg/m³. The data was then grouped by date, and frequency of each date occurring was counted. In order to identify persistent, high pollution weather regimes, dates with twelve or fewer hourly exceedances were removed. Finally, date clusters of three days or longer were identified. In total, there were 9 clusters identified. These clusters were then labelled based on their Kidson Synoptic Class, from the Kidson Type time series (Renwick 2021). The synoptic class identification matched previous literature analysis showing anticyclones correlated to high pollution events (Jiang et al. 2005, 2014). The regime “H”, representing a high-pressure regime centered over the North Island, accounted for 46% of the days, followed by regimes “HSE” and “HW” at 17%. The high-pressure regimes were in the blocking and zonal categories, both of which are associated with low wind speeds in the Auckland region and reduced precipitation (Kidson 1999). In order to select one high pollution anticyclone, a series of constraints were applied to the 9 scenarios. First, given that the predominate synoptic type was H, the event had to be majority H regimes, which three episodes fit. Second, the event needed to be validated in terms of persistence of pollution based on the pollution persistence event identifier tool in openair (Carslaw and Ropkins 2012). This was run with a threshold of 15.7 μg/m³ with a minimum persistence of eight hours. All three episodes appeared on the persistence event identifier, but only one event appeared for more than one day of persistence. This event was a primarily Kidson type H anticyclone occurring on June 26th-30th, 2019. It was a calm few days, with windspeeds under 3 m/s. There is a clear pattern of diurnal sea breezes, with wind peaking around midday at both Penrose and Takapuna sites. While there is a diurnal cycle, the sea breezes in the winter months are much weaker than the summer due to a reduced thermal gradient. Therefore, the late June 2019 event chose here represents a high-risk meteorology for reduced dispersion and increased pollution episodes due to the synoptic meteorology and season.

The four-day period selected demonstrated persistent poor air quality and a synoptic regime associated with high ground level pollution. In Auckland, most high pollution days occur with low windspeeds, but are not correlated with a specific wind direction (Jiang 2008) (Supplementary S4). Seasonally, the climatology of Auckland shows decreased windspeeds and turbulent mixing during the autumn and winter months, partially due to the reduction in diurnal cycles (Khan 2010). As noted previously, these conditions were selected in order to model a worst-case scenario whereby a period of high SO₂ fluxes coincides with meteorological conditions promoting limited dispersion, and therefore high ground-level concentrations, of SO₂.

The surrounding water, complex topography and land use of the Auckland region has a significant impact on the mesoscale wind patterns. The mountainous ranges of the region, including the Waitākere ranges to the west, and Hunua ranges to the southeast, shelter inland regions from the synoptic winds and sea breezes (Fig. 4 and Supplementary S1 and S2). High variation in terrain heights, cooler land surface moisture and lower temperatures act to reduce sea breeze-induced wind speed intensities and delay their onset (Khan 2010). Coastal areas of the city also tend to be windier compared to inland areas as weather fronts pass through (Chappell 2013). The volcanic hills within the city also act as wind shelters at low elevation and produce higher wind speeds at higher elevation (Supplementary S3). As a result, it is very difficult to accurately simulate multiscale processes (from the microscale to mesoscale) but it has been previously shown (Khan 2010; Soltanzadeh et al. 2017; Pretorius et al. 2019) that sea-breeze and other terrain-induced wind systems can be modeled accurately using a numerical weather model setup with a grid scale of 1 km or less (in our case study, 300 m). This ensured realistic meteorological forcing parameters were carried through to the SO₂ dispersion model.

Fig. 4

Auckland region sea breeze convergence zone and topography. Sea breezes (blue) occur primarily from the southwest and northeast, meeting in the convergence zone (enclosed by dotted lines). Early in the day, harbor/bay breezes (red) also occur. Adapted from Fig. 8, Chappell 2013

Meteorological modeling

The meteorological conditions of the chosen study period were reproduced with the mesoscale numerical Weather Research and Forecasting (WRF) model (Skamarock et al. 2021), configured to 300 m of spatial resolution in the vicinity of Scenario C eruption location (Fig. 5). For this study, WRF was configured with 4 nested domains, from an outer domain with 8.1 km spatial resolution (domain 1 (d01) in Fig. 5) to an inner domain with 300 m spatial resolution (d04 in Fig. 5). This progressive downscaling technique is typically used in numerical modeling and results from a compromise between the necessary detail of relevant weather features (localized wind patterns, complex topography), which are more relevant near the eruption location, and the heavy computational needs of a high-resolution simulation. The spatial resolution of 300 m is necessary to resolve the meteorologic impact of the complex topography in the AVF area, such as the diurnal development of an atmospheric boundary layer and the thermal wind systems associated with land/sea contrast and topographic slopes (Fig. 6).

Fig. 5

WRF model configuration with progressive resolution increases, from 8.1 km (d01) to 2.7 km (d02), 900 m (d03) covering AVF region, and finally 300 m (d04) in the nearby vicinity of the eruption location

Fig. 6

A Terrain and station locations and B Land Use Index for the study area. Note the highly urbanized landscape. Land Use Index: Urban = 1, Ocean = 16, 12 = Deciduous

In our case study application, WRF was configured with the options listed in Table 2. WRF simulations were forced with the publicly available global reanalysis from ERA5 data set (Hersbach et al. 2020), with 0.25 degrees resolution (approximately 25 km) and daily sea surface temperature from GHRSST Analysis with 0.1 degrees resolution (Toshio et al. 2017).

Table 2 WRF model configuration options used in this study (details in Skamarock et al. 2021)

Atmospheric dispersion modeling

Atmospheric dispersion modeling is a useful tool for understanding how emissions from natural hazard events can affect air quality by providing a prognosis of the spatial and temporal distribution of airborne gases and particulate matter concentrations. Typically, a mesoscale meteorological model, such as WRF, is used to model the relevant phenomena (happening between 10 and 1000 km scales) near the site of the event. The meteorological fields can then be coupled with a Lagrangian particle dispersion model, which simulates the movement of “tracer particles” in the supplied meteorological conditions (Hegarty et al. 2013).

There are a number of Lagrangian dispersion models. In volcanology, one of the most widely used for dispersion of volcanic ash and gases is the Hybrid Single Particle Lagrangian Integrated Trajectory (HYSPLIT) model (Davidson et al. 2009; Stein et al. 2015a, b; Hurst and Davis 2017; Trancoso et al. 2022). The Vog Measurement and Prediction (VMAP) project in Hawai’i, USA, uses HYSPLIT to produce comprehensive vog forecasts (Summary of the Vog Measurement and Prediction (VMAP) Project (2021)). In South Africa, HYSPLIT was utilized to model long-range transport of SO₂ from a 2015 eruption occurring in South America (Sangeetha et al. 2018). In New Zealand, HYSPLIT is operationally run by MetService, for airborne volcanic ash, as part of one of the nine international Volcanic Ash Advisory Centers, and since 2014 MetService has also run HYSPLIT to produce pre-emptive ashfall forecasts in a joint effort with GNS Science (Trancoso et al. 2022). Both models are driven by calibrated high resolution WRF forecasts developed at MetService (Hurst and Davis 2017).

In our study, HYSPLIT was used to model dispersion of SO₂ over the first four days of the eruption period (June 27th-30th, 2019), using the meteorological parameters from the WRF modeling described in the previous section, and utilizing the parameters of the VMAP standard (Businger et al. 2015). The HYSPLIT model was set up with a plume of constant vertical dispersion from 10 to 2500 m above ground level. The plume was released for the entirety of the 96 h of simulation to model the highest emission period of the eruption, with a constant flux for simplicity. The initial emission is 100% SO₂, however there is a 1% per hour conversion to H₂SO₄ as a sulfate aerosol (Businger et al. 2015). While the conversion rate does vary diurnally, the complex inputs required to the model to simulate any changes in the conversion were not considered here for simplicity. Both SO₂ and H₂SO₄ were removed through dry deposition, again for simplicity. The boundary layer turbulent mixing parameter was given by the forcing WRF model’s Turbulent Kinetic Energy field. 20,000 particles were released per cycle (hour) for a high-resolution output. The model output was integrated in vertical layers from surface up to 100 m above ground level, with concentrations for both pollutants averaged hourly to best visualize ground level concentrations and remain consistent with standard observations.

Finally, HYSPLIT was run in deterministic mode (single output), as well as in ensemble mode varying the turbulence (randomly) and meteorology parameters (one grid spacing in each direction) to obtain uncertainty estimates.

Results and discussion

Meteorology modeling

The WRF model was validated using hourly surface wind speed, temperature, and wind directions in seven observation stations around Auckland (Fig. 6A). The diurnal cycle of the boundary layer stability (or the tendency for vertical motion represented by the thermal gradient profile; Stull et al. 2017) was well represented as shown in Fig. 7, where we compared day and night vertical thermal gradient profiles at Penrose stations. As the day progresses with radiative warming of the surface, the boundary layer vertical structure transitions from stable (nighttime conditions; Fig. 7) to near neutral and unstable in the late afternoon (Fig. 7, daytime conditions) which is conducive for enhanced vertical mixing. This also results in increased surface wind speed associated with the sea breeze inflow. In contrast, the colder early morning conditions show a positive lapse rate in potential temperature for the first 200 m above ground level, which is representative of a temperature inversion and negative buoyancy, suppressing vertical mixing of pollutants and reducing surface wind speeds.

Fig. 7

Vertical potential temperature profile at night (3AM) and day (3PM) from the WRF model at the Penrose station location. Height is above ground level

Atmospheric dispersion modeling

Sulfur dioxide (SO₂)

Modeled 24-h average ground-level concentrations of SO₂ over Auckland for the first four days of the eruption scenario are shown in Fig. 8. Under these dispersion conditions, SO₂ is dispersed primarily to the northwest quadrant. The effect of topography on dispersion is evident, with the plume trapped in regions of low elevation bordered by the Waitākere Ranges to the west and higher ground to the north of Auckland. Very high 24-h average SO₂ concentrations of > 1000 µg/m³ extend over hundreds of square km, persisting throughout the duration of the scenario. The area affected includes Auckland’s city center (CBD, see Fig. 1), and at times, Auckland’s central and western suburbs and the North Shore (Fig. 8).

Fig. 8

24-h averages of SO₂ concentration over Auckland in Domain 4, from 0–100 m above ground level at a 500-m resolution. The grey line encloses areas exceeding 100 µg/m³

One-hour average SO₂ concentrations, shown as an animation over the 96 h of the scenario, are included as Additional File 1. As for the 24-h averages, SO₂ is dispersed primarily to the northwest quadrant, with occasional dispersion to the southwest. Smaller areas are sometimes subjected to hourly concentrations exceeding 10,000 µg/m³, even at distances > 10 km from the vent. These concentrations are far in excess of typical ambient SO₂ concentrations in New Zealand. Even at the seven sites routinely monitored due to their proximity to known industrial sources of SO₂, 24-h average concentrations between 2017 and 2020 ranged from 0.4–10.8 µg/m³ (Stats NZ 2022).

Modeled ground-level SO₂ concentrations in this study are broadly comparable with those recorded for recent basaltic eruptions worldwide. For instance, in the 2014–2015 Bárðarbunga-Holuhraun fissure eruption (Iceland), described as ‘one of the most intense, large-scale volcanogenic air pollution events in centuries’ (Ilyinskaya et al. 2017), 15-min average SO₂ concentrations of 43,000 µg/m³ were recorded in the airborne plume close to the vent, and 3100 µg/m³ in the grounding plume 4–20 km from the vent (Ilyinskaya et al. 2017). Similarly, for the 2018 lower East Rift Zone eruption of Kīlauea volcano (Hawaiˈi, USA), SO₂ concentrations spiking at over 100 ppm (260,000 µg/m³) were recorded near the erupting fissures, with concentrations of approximately 5 ppm (13,000 µg/m³) recorded in residential areas several km away (Stewart et al. 2019).

New Zealand and international guidelines for SO₂ are shown in Table 3, for different exposure periods. The NZ National Environmental Standards for Air Quality (NES-AQ) are legally enforceable regulations that set lower and upper limits for SO₂ concentrations for 1-h exposures, designed to be protective against short-term exposures to elevated concentrations. The NZ National Ambient Air Quality Guidelines (NAAQG) are based on longer exposure timeframes and are designed to protect communities from longer-term health impacts but are not legally enforceable. In Table 3, we also include the US Environmental Protection Agency’s Acute Exposure Guideline Levels (USEPA-AEGLs), which describe human health effects from occasional exposures to airborne pollutants, for SO₂.

Table 3 USEPA Acute Exposure Guideline Levels (AEGL1-3) and New Zealand (Land Air Water Aotearoa 2023) and international (WHO 2021; US EPA 2023) guideline levels for SO₂ concentrations (µg/m³)

Modeled SO₂ concentrations (Fig. 8) suggest that life-threatening concentrations (79,000 µg/m³ for exposure periods between 10 min and 1 h, and 25,000 µg/m³ for an exposure period of 8 h) may occur close to the vent, generally within 1 km. The Auckland Volcanic Field Contingency Plan (Auckland Emergency Management 2015) stipulates that an evacuation zone will be set up around the inferred location of a new eruption consisting of a Primary Evacuation Zone (PEZ) extending 3 km radially from the vent, and a Secondary Evacuation Zone extending a further 2 km radially from the PEZ (Deligne et al. 2017). These emergency management measures should mitigate the life safety hazard from SO₂, but regular monitoring of ground-level SO₂ concentrations will be necessary to inform the management approach.

Under the modeled scenario, large areas of Auckland will be subjected to sustained concentrations of SO₂ above the 24-h NAAQG recommendation (120 µg/m³) throughout the four-day duration of the eruption, indicating an increased risk of adverse health effects, particularly for sensitive groups such as asthma sufferers. The level of risk is heightened relative to the more protective 24-h guideline of 40 µg/m³ set by the World Health Organization.

Within the above area of heightened risk, hundreds of square km of residential areas will be exposed to 24-h average concentrations > 1000 µg/m³ (Fig. 8), and extensive residential areas will be exposed to shorter-term one-hour average concentrations exceeding both AEGL-1 (the threshold for the onset of discernible effects) and AEGL-2 (the threshold for the onset of more serious or irreversible effects) (Additional File 1). Strategies for managing exposure to these concentrations of SO₂ include the following: advising the public to remain indoors, keeping doors and windows tightly closed and turning off heat pumps/air conditioning units pulling in outdoor air, minimizing sources of indoor air pollution such as smoking and burning candles or incense, and reminding individuals with respiratory or heart conditions to keep their prevention and relief medications at hand and use as prescribed (DEVORA 2022).

In addition to the well-characterized irritant effects of SO₂, Orellano et al. (2021) concluded, in a recent review and meta-analysis, that a short-term rise in SO₂ concentrations (from hours to days) increases the risk of all-cause and respiratory mortality.

Risk ratios for both 1 h and 24-h exposure, and all-cause and respiratory mortality, are shown in Table 4. These refer to the ratio by which deaths increased for every 10 µg/m³ increase in SO₂ concentration. For example, the daily risk ratio for SO₂ of 1.0059 means that for every 10 µg/m³ increase in the 24-h average concentration of SO₂, deaths in the general population from all causes increased by 0.59%, with 95% of the data in the studies analyzed lying between 0.46% and 0.71%.

Table 4 Risk ratios for short term exposure to SO₂ (from Orellano et al. 2021)

While risk ratios may be small, the health burden can still be significant when applied across a large population. For example, a simple first-order approach can be taken for the areas with concentrations ≥ 1000 µg/m³, for which a risk ratio of ≥ 1.59 (1 + 100*0.0059) will apply. For a population of 200,000, the base daily mortality rate of 1.9 would increase by ≥ 59%, or ≥ 1.1 additional deaths per day, attributable to SO₂ exposure (Mortality web tool 2020).

Zheng et al. (2021) evaluated evidence on the effects of short-term exposure to SO₂ and other air pollutants and reported an association between 24-h average SO₂ concentrations and asthma emergency department (ED) visits and hospitalizations, with a risk ratio of 1.010 and a 95% confidence interval of 1.001 – 1.020. For a population exposed to 24-h average SO₂ concentrations ≥ 1000 µg/m³, the baseline rate of asthma ED visits and hospitalizations would increase by ≥ 200%.

We emphasize here that these simple calculations are indicative only and have substantial uncertainties. In particular, while the relationships between risk ratios and SO₂ concentrations were linear for the range of studies reviewed by Orellano et al. (2021) and Zheng et al. (2021), the concentrations modeled here are higher than for those studies, introducing uncertainty into the practice of scaling up the risk ratios linearly. A further source of uncertainty is the potential for interactions between SO₂ and other urban air pollutants such as ozone (O₃) and nitrogen dioxide (NO₂).

Sulfate aerosol (H₂SO₄)

In the atmosphere, SO₂ reacts with oxygen and water vapor, in the presence of sunlight, to form sulfate aerosol particles of ~ 0.1–0.3 µm diameter. Modeled daily average ground-level concentrations of sulfate aerosol (H₂SO₄) formed by oxidation of SO₂ (see “Methods” section) are shown in Fig. 9. One-hour average sulfate aerosol concentrations, shown as an animation over the 96 h of the scenario, are included as Supplementary Video 1 and 2 (Additional File 2).

Fig. 9

Daily average concentrations of sulfate aerosol over Auckland in Domain 4, from 0–100 m above ground level at a 500-m resolution. The grey line shows areas exceeding 25 µg/m³

As for SO₂, sulfate aerosol is dispersed primarily to the northwest quadrant (Fig. 9). However, higher concentrations are slightly offset from the eruption location, consistent with the slow conversion of SO₂ gas to aerosol (Businger et al. 2015). For the first day of the scenario, an extensive area of hundreds of square km will experience sulfate aerosol concentrations exceeding 25 µg/m³, although this area will decrease over the following days. Similarly, the area experiencing > 100 µg/m³ decreases substantially over the following days.

The most appropriate metric to assess the significance of the modeled sulfate aerosol concentrations is PM_2.5 (particulate matter smaller than 2.5 µm in diameter). New Zealand and WHO guidelines for PM_2.5, for time averaged periods of 24 h and one year, are shown in Table 5. There is currently no legally enforceable NES-AQ for PM_2.5 in New Zealand, although the Ministry for the Environment has proposed establishing NES-AQ for PM_2.5 based on the current NAAQG (Ministry for the Environment 2020).

Table 5 New Zealand and International (WHO) guideline levels for PM_2.5 concentration (µg/m³)

Modeled daily average concentrations of sulfate aerosol (Fig. 9) show substantial areas exceeding NZ and WHO 24-h guidelines, indicating heightened risk. Here we do not further investigate the health implications of the modeled sulfate aerosol concentrations for the following reasons. Firstly, there are existing sources of anthropogenic PM_2.5 in Auckland, which will be additive with sulfate aerosol from a new AVF eruption. Conditions that favor limited dispersal of SO₂ and sulfate aerosol will also favor the buildup of anthropogenic PM_2.5. Secondly, while strong causal links between even short-term exposure to fine PM and adverse health effects are now well-established (Orellano et al. 2020), it is currently unknown whether this includes PM of volcanic origin (Stewart et al. 2022).

Model uncertainty and ensemble modeling

Ensemble modeling is a method of assessing uncertainty and improving accuracy, constructed by combining multiple runs that vary one or more input and/or dependent parameters. HYSPLIT has three built-in ensemble modeling capabilities to measure uncertainty: meteorology, turbulence and physics (Stein et al. 2015a, b). In this study, meteorology and turbulence ensembles were run based on what ensembles are utilized in methodology in current literature (Stein et al. 2007; Holland et al. 2020). Turbulence ensemble is utilized in order to determine the uncertainty of the model’s solution to the variability of random motions of turbulence by varying the random number generator used to simulate random motion. The meteorology ensemble tests the uncertainty with the meteorological input by offsetting the meteorological grid input by one or more grid points and thus accounting for the uncertainty in the spatial predictions (Stein et al. 2015a, b).

Meteorology ensemble

In the meteorology ensemble, for the probability of exceeding the daily NAAQG for SO₂, the highest likelihood is to the west of the eruption location, with a mostly symmetrical plume (Fig. 11). The highest percentage area is 50%, suggesting a degree of variability within the model members. This is supported by the coefficient of variation mapping for SO₂, which shows a relatively high variation for a large amount of the study area. As expected from a near-surface plume in complex topography, a variation in eruption location, however slight, has a significant impact on plume dispersion patterns. For example, given the eruption location on a narrow body of water, a minor shift in the eruption location on the grid may result in the eruption location being on land and therefore affected by different meteorology patterns.

The high degree of variation in the coefficient of variation in the meteorology ensemble (Fig. 10) also suggests that errors in the weather driving fields, either in time or in space can have a significant impact on plume dispersion. Therefore, the scenario modeled is limited to the specific location of the eruption and the quality of the WRF model simulations, reinforcing the need of having probabilistic model results from full meteorological ensemble suites (Trancoso et al. 2022). For this scenario, the meteorology ensemble modeling further supports the impact of the complex terrain in Auckland and the resulting mesoscale meteorology on pollutant dispersion. This can be seen in the distinct impact of the converging sea breezes on the plume direction (Supplementary Video S1 and S2).

Fig. 10

Meteorology ensemble for the eruption period (June 27th-30th, 2019). A An average probability of exceeding the daily average NAAQG of 120 µg/m³ of SO₂, and B) coefficient of variation

Turbulence ensemble

The turbulence ensemble over the full eruption period (Fig. 11) shows a larger variability in where the plume could exceed 120 ug/m3 in the meteorological vs the turbulence ensemble, which suggests that meteorology is having a larger impact on the model than the turbulence number (Fig. 10). This suggests that meteorology accuracy has greater impact on dispersion patterns than varying the model’s turbulence parameters. The plume extends in the same W-NW direction for both types of ensembles, but the higher probability (> 75%) in the turbulence ensemble shows a reduced variation between model members. This is supported further by the coefficient of variation mapping (Fig. 11) which shows much lower coefficients of variation (< 25) for the turbulence ensemble compared to the meteorology ensemble (< 125). The highest (< 500) coefficient of variation is at the edges of the overall plume, suggesting that the main direction of the plume is consistent throughout the model members with the greatest variation on the spread of the plume to the east and south. The high coefficient of variation values suggest that the standard deviation is large compared to a small mean.

Fig. 11

Turbulence ensemble for the eruption period (June 27th-30th, 2019). An average probability of exceeding the daily average NAAQG of 120 µg/m³ SO₂, and B coefficient of variation

Concluding remarks

In this contribution, we present a new methodology for simulating volcanic gas dispersion for a new eruption at the Auckland Volcanic Field and apply this model to a credible ‘worst-case’ eruption and meteorologic scenario. Under the ‘worst-case’ dispersion and eruption scenario, modeled ground level concentrations of SO₂ greatly exceed typical ambient SO₂ concentrations in the Auckland region and are comparable to concentrations measured for recent basaltic eruptions in Iceland and Hawai’i. Life-threatening concentrations may occur near the vent; this hazard will need to be managed by establishing specific evacuation zones around the vent, ideally informed by real-time monitoring of SO₂ concentrations. Large areas of Auckland will be subjected to concentrations exceeding effects thresholds. Public health strategies will include advising the public to remain indoors with doors and windows closed and turning off heat pumps/air conditioning units pulling in outdoor air. Simple risk ratio calculations indicate that some areas of the city may experience a doubling of baseline rates for asthma emergency department visits and hospitalizations. Acknowledging that this ‘worst-case’ scenario is comparatively unlikely, it does highlight that gas hazards from a future AVF eruption could create substantial health and associated emergency management challenges.

The major uncertainty in predicting ground-level SO₂ dispersion is the location, size and style of a future eruption in the AVF. Future work should address this limitation by modeling various eruption types in varying locations to better understand the impact of size, style and location of eruption on SO₂ dispersion patterns. The other limitation of the study is that the eruption is only modeled for a specific meteorological period, and therefore does not represent all possible meteorological conditions. In future studies, modeling each DEVORA scenario for each synoptic weather type will provide a more comprehensive assessment of the SO₂ risk to the greater Auckland area in the event of an AVF eruption. Finally, future work should also address how volcanic gas dispersion mapping could be utilized to inform a population analysis of health impacts across the Auckland region.

This work presents a promising approach for rapidly supporting public health and emergency management operations in the event of an imminent eruption. The development of the coupled weather and SO₂ dispersion model simulations enables dispersion models to be run for any AVF eruption scenario in a variety of synoptic weather conditions, including all DEVORA eruption scenarios and (potentially) in a future AVF eruption. Given the relatively fast computation time (provided that good-quality meteorological forcing is available) of less than an hour for an ensemble model, gas dispersion modeling is an effective tool that could be used before, during and after an eruption event for health agencies to evaluate risk of volcanic gas exposure to communities and inform appropriate emergency management strategies.