How well are aerosol–cloud interactions represented in climate models? – Part 1: Understanding the sulfate aerosol production from the 2014–15 Holuhraun eruption

For over 6 months, the 2014–2015 effusive eruption at Holuhraun, Iceland, injected considerable amounts of sulfur dioxide (SO₂) into the lower troposphere with a daily rate of up to one-third of the global emission rate, causing extensive air pollution across Europe. The large injection of SO₂, which oxidises to form sulfate aerosol (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$), provides a natural experiment offering an ideal opportunity to scrutinise state-of-the-art general circulation models' (GCMs) representation of aerosol–cloud interactions (ACIs). Here we present Part 1 of a two-part model inter-comparison using the Holuhraun eruption as a framework to analyse ACIs. We use SO₂ retrievals from the Infrared Atmospheric Sounding Interferometer (IASI) instrument and ground-based measurements of SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ mass concentrations across Europe, in conjunction with a trajectory analysis using the Hybrid Single Particle Lagrangian Integrated Trajectory (HYSPLIT) model, to assess the spatial and chemical evolution of the volcanic plume as simulated by five GCMs and a chemical transport model (CTM). IASI retrievals of plume altitude and SO₂ column load reveal that the volcanic perturbation is largely contained within the lower troposphere. Compared to the satellite observations, the models capture the spatial evolution and vertical variability of the plume reasonably well, although the models often overestimate the plume altitude. HYSPLIT trajectories are used to attribute to Holuhraun emissions 111 instances of elevated sulfurous surface mass concentrations recorded at European Monitoring and Evaluation Programme (EMEP) stations during September and October 2014. Comparisons with the simulated concentrations show that the modelled ratio of SO₂ to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ during these pollution episodes is often underestimated and overestimated for the young and mature plume, respectively. Models with finer vertical resolutions near the surface are found to better capture these elevated sulfurous ground-level concentrations. Using an exponential function to describe the decay of observed surface mass concentration ratios of SO₂ to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ with plume age, the in-plume oxidation rate constant is estimated as 0.032 ± 0.002 h⁻¹ (1.30 ± 0.08 d e-folding time), with a near-vent ratio of 25 ± 5 (µg m⁻³ of SO₂ $/$ µg m⁻³ of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$). The majority of the corresponding derived modelled oxidation rate constants are lower than the observed estimate. This suggests that the representation of the oxidation pathway/s in the simulated plumes is too slow. Overall, despite their coarse spatial resolutions, the six models show reasonable skill in capturing the spatial and chemical evolution of the Holuhraun plume. This capable representation of the underlying aerosol perturbation is essential to enable the investigation of the eruption's impact on ACIs in the second part of this study.

1 Introduction

The Bárðarbunga volcanic system in Iceland began experiencing noteworthy sequences of earthquakes on 16 August 2014 (Ágústsdóttir et al., 2016; Gudmundsson et al., 2014; Sigmundsson et al., 2015). This seismic activity created an effusive eruption at Holuhraun (64.85^∘ N, 16.83^∘ W) from 31 August 2014 to 27 February 2015. The resulting flow of lava is estimated to have covered 84 km² with an average discharge rate of 90 m³ s⁻¹, making it the largest effusive eruption in Iceland since the 1783–1784 Laki eruption (Pedersen et al., 2017). Ground-based observation estimates suggest the Holuhraun eruption released a total of 9.6–11.8 Mt of sulfur dioxide (SO₂) with little tephra (Gíslason et al., 2015; Pfeffer et al., 2018). These emissions represent up to 215 times the amount of Icelandic anthropogenic SO₂ emissions for 2019 (Keller et al., 2022) and approximately one-tenth of the current global annual anthropogenic SO₂ emissions (e.g. Aas et al., 2019). During these months of intense degassing, several studies using ground-based observations and remote sensing suggest that the volcanic plume reached altitudes of 1–4 km (Arason et al., 2015; Carboni et al., 2019a; Flower and Kahn, 2020; Pfeffer et al., 2018). This release of SO₂ adversely affected air quality over large parts of Iceland, often exceeding World Health Organization (WHO) surface concentration limits (Gíslason et al., 2015; Ilyinskaya et al., 2017; Schmidt et al., 2015). Such high rates of SO₂ emission into a relatively pristine, unpolluted environment provide an ideal opportunity to observe aerosol–cloud interactions (ACIs) and elucidate aerosol impacts on the climate system (e.g. Breen et al., 2021; Chen et al., 2022; Gettelman et al., 2015; Haghighatnasab et al., 2022; Malavelle et al., 2017; McCoy and Hartmann, 2015; Toll et al., 2017).

Remote sensing data estimates of SO₂ concentrations, a common tool to monitor the evolution of volcanic plumes, show that the September and October meteorological conditions transported the plume across Europe including the Fennoscandian Peninsula (Grahn et al., 2015; Ialongo et al., 2015), Belgium, northern France (Boichu et al., 2016), the UK, Ireland, the Netherlands (Schmidt et al., 2015), and Germany (Steensen et al., 2016). Previous studies have combined satellite data with ground-based observations and trajectory modelling to attribute local pollution events to the Holuhraun eruption and investigate the plume characteristics (e.g. Boichu et al., 2019; Schmidt et al., 2015; Twigg et al., 2016). Understanding such characteristics, particularly those that are hazardous to health (i.e. fine particulate matter), is important for air quality monitoring and attribution of any exceedances of statutory thresholds (Heaviside et al. 2021; Stewart et al., 2022). However, most studies assessing Holuhraun impacts on air quality focus only on pollution experienced at a few ground-based stations across small geographical areas, and the need to increase the quantity of air monitoring stations investigated has been noted previously (Schmidt et al., 2015; Twigg et al., 2016).

In addition to observational evidence, many studies have explored the impacts of the Holuhraun eruption using numerical models to simulate the transport of the volcanic SO₂ emissions. For instance, the Icelandic Meteorological Office employed the CALPUFF dispersion model for near-time probabilistic hazard mapping (Barsotti, 2020) and to support aviation safety decisions (Barsotti et al., 2020) following the eruption. Boichu et al. (2016) and Steensen et al. (2016) used the EMEP MSC-W and the CHIMERE chemical transport models (CTMs), respectively, to explore the far-range air pollution effects caused by the eruption, whereas Schmidt et al. (2015) used the NAME dispersion model to do so. Recently, Haghighatnasab et al. (2022) analysed the results from high-resolution simulations performed with the ICON model in its numerical weather prediction mode to explore the impact of the aerosol perturbation introduced by Holuhraun on cloud properties. This impact was also examined in earlier works using general circulation models (GCMs) of coarser resolutions (CAM5 – Gettelman et al., 2015; HadGEM3, CAM5, and a NorESM variant – Malavelle et al., 2017). Considering the opportunity the Holuhraun eruption presents to assess the modelling capability of GCMs in capturing aerosol–cloud interactions, the number of GCM studies of the Holuhraun eruption to date is surprisingly low.

Here we present results from an inter-model comparison two-part study of the volcanic plume and its interactions with clouds within the vicinity of the Holuhraun eruption (44 to 80^∘ N, 60^∘ W to 30^∘ E) during September and October 2014. Participation in the study was organised through the AeroCom initiative (Schulz et al., 2006). We begin with a comparison of the volcanic SO₂ plume evolution between remote sensing observations and simulations of the eruption from five GCMs and a CTM. The analysis specifically focuses on the plume spatial distribution, plume altitude, and the total SO₂ mass burden. We further investigate the numerical models' capability to simulate the Holuhraun eruption by comparing modelled SO₂ and sulfate (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$) concentrations with in situ surface observations from 25 long-term monitoring stations from the European Monitoring and Evaluation Programme (EMEP) network. By using remote sensing and ground-based observations, in conjunction with trajectory modelling, we attribute sulfurous pollution events to Holuhraun emissions and assess the skill of the numerical models in capturing these episodes. Finally, this refined catalogue of volcanically influenced pollution events is used to estimate the rate at which SO₂ oxidises to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ within both observed and modelled volcanic plumes. We conclude with a discussion of whether the models simulate the observed Holuhraun aerosol perturbation with sufficient fidelity, a prerequisite for using a model to understand the climatic impacts caused by the eruption via ACIs (see Part 2 of this study, currently in preparation).

2 Methodology

We provide a brief description of the remote sensing and in situ observations that are used to assess the numerical models, the numerical models themselves, and the HYSPLIT trajectory modelling framework used to evaluate the pollutant transport of the observed local sulfurous events.

2.1 Satellite observations

Retrievals of volcanic SO₂ from satellite instrumentation typically use either measurements in the ultraviolet (e.g. Ozone Mapping and Profiler Suite – Nadir Mapper (OMPS-NM); Carn et al., 2015; Li et al., 2017; Wells et al., 2023; TROPOspheric Monitoring Instrument (TROPOMI); de Leeuw et al., 2021; Theys et al., 2017; Global Ozone Monitoring Experiment-2 (GOME-2); Twigg et al., 2016) or infrared (e.g. Infrared Atmospheric Sounding Interferometer (IASI); Clarisse et al., 2008, 2010; Haywood et al., 2010; de Leeuw et al., 2021) region of the electromagnetic spectrum. Here we use IASI measurements as they have proved valuable in monitoring the evolution of volcanic plumes in both the stratosphere (e.g. Haywood et al., 2010; de Leeuw et al., 2021) and the troposphere (e.g. Athanassiadou et al., 2016; Malavelle et al., 2017). Specifically, we use data from IASI retrievals on the MetOp-A and MetOp-B satellites produced by the University of Oxford as part of the NERC Centre for the Observation and Modelling of Earthquakes, Volcanoes and Tectonics (COMET) (Carboni et al., 2019a, b).

SO₂ column load and plume height are derived by applying the IASI retrieval algorithm of Carboni et al. (2012, 2016) to level 1C data from EUMETSAT and the CEDA Archive. The IASI SO₂ retrieval is performed only on pixels where the underlying SO₂ detection scheme returns a positive result. The detection scheme is a linear retrieval where a positive result is defined as when the free parameter, the SO₂ column load, exceeds a defined threshold. This threshold is set substantially greater than the standard deviation of the assumed Gaussian distribution describing the background atmospheric concentration of SO₂. Consequently, a positive result is exceedingly likely to be significantly different to the background and not a consequence of instrumental noise or climatological variations (see details in Walker et al. 2011, 2012). The threshold defined for the Holuhraun eruption in Carboni et al. (2019a) is 0.49 effective DU.

An iterative optimal estimation retrieval using forward modelled spectra is applied to pixels with a positive detection result. This retrieval uses all channels within 1000–1200 and 1300–1410 cm⁻¹ (the 7.3 and 8.7 µm SO₂ bands, respectively) and assumes a Gaussian vertical SO₂ profile to return the SO₂ column load (DU) and height (mb) which the retrieval algorithm subsequently converts to kilometres using European Centre for Medium-Range Weather Forecasts (ECMWF) meteorological profiles. The algorithm provides a comprehensive pixel-by-pixel error estimate on the retrieved parameters that is derived from an error covariance matrix computed using the differences between the measured IASI spectra and the simulated spectra (driven by ECMWF data). This means uncertainty due to imperfect knowledge of non-SO₂ atmospheric conditions (e.g. cloudiness, vertical distribution of constituents) and imperfect radiative transfer simulations are addressed (see details in Carboni et al., 2012). The thermal contrast between the plume and surface heavily influences the retrieval error such that retrievals of SO₂ plumes centred at lower altitudes have higher uncertainties. Note that the IASI retrieval algorithm can miss parts of the SO₂ plume, such as when overlaying clouds are present or under conditions of negative thermal contrast, and so the IASI SO₂ column load and mass burdens presented here should be considered an approximate minimum.

This study maps data from individual IASI overpasses to a regular 1.0^∘ × 1.0^∘ latitude–longitude grid using a nearest-neighbour approach with Gaussian weighting. The decision to weight closer neighbouring pixels allows retention of plume characteristics which can change abruptly over small spatial scales. The individual gridded overpass data are grouped into bidaily intervals (AM: 03:30–15:30 UTC, PM: 15:30–03:30 UTC) with overlapping cells averaged. Linear interpolation is used to estimate missing values in the gridded output that result from orbital gaps and/or pixels failing quality control. Each bidaily regridded IASI SO₂ column load and altitude map is visually inspected to ensure no obvious artefacts exist within the Holuhraun vicinity.

Table 1 Details of the 25 EMEP stations explored in this study. These stations are shown geographically in Fig. 1.

Table 2 Details of the numerical models used in this study.

2.2 Surface observations

Since the early 1970s, the EMEP network has monitored air pollution and surface deposition across Europe at ground-level stations outside of notable conurbations where significant sources of local pollution are minimised, thus creating a comprehensive database useful for assessing long-range transportation of a plethora of pollutants (Tørseth et al., 2012). The use of EMEP stations to evaluate model output has proven fruitful previously (e.g. Hardacre et al., 2021; Mulcahy et al., 2020). This study only considers EMEP stations that provide both SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ surface mass concentration measurements at the same temporal sampling frequency during September and October 2014. This criterion results in 25 stations located across 12 countries being selected for this study (see Table 1). The observations include hourly and daily measurements made using online ion chromatography and filter-pack measurements, respectively, with the former to a precision of 0.001 µg m⁻³ and the latter to either 0.01 or 0.001 µg m⁻³ depending on the station. The hourly and daily sampling midpoints are centred on 30 min past the hour and on the hour, respectively. Further details on the instruments and sampling techniques are provided in the EMEP Standard Operating Protocol (NILU, 2014). This study screens out invalid and missing measurements in accordance with the EMEP data quality flags (NILU, 2020). For each station monthly surface mass concentration climatologies for SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ are calculated from the temporal coverage listed in Table 1. For a given station this coverage may differ for the two chemical species. Subsequently, the combined total sulfur content climatologies are only calculated across periods where the temporal coverages align (e.g. 1988–2017 for Aspvreten and 2006–2020 for Irafoss). Here we define a significant sulfurous pollution event as when the surface mass concentration of the total sulfur content observed exceeds the 90th percentile of the corresponding monthly climatological value. Note that the number of EMEP stations carrying out SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ measurements has fallen since the late 2000s due to the reduced need to monitor the declining sulfur emissions from anthropogenic sources (Boichu et al., 2019; Schmidt et al., 2015).

2.3 Numerical model simulations

Included in this study are Holuhraun eruption simulations by five GCMs: UKESM1.0, HadGEM3-GA7.0, MIROC6.1-SPRINTARS, ECHAM6.3-HAM2.3, and ECHAM6.3-HAM2.3-P3. Simulations are performed using the atmosphere-only component at a global scale (AMIP-style). To help clearly discriminate between signal and noise, the modelled horizontal winds and potential temperature are constrained (“nudged”) to ERA-Interim reanalysis data (Dee et al., 2011) on a 6-hourly timescale and use monthly observational datasets to prescribe sea surface temperature and sea ice boundary conditions (e.g. HadISST, Rayner et al., 2003). All other modelled variables evolve physically and dynamically as their setup dictates and are subject to the parameterisations in play. Also included in our inter-model comparison is OsloCTM3, a global CTM. Unlike GCMs, CTMs do not simulate atmospheric dynamics explicitly; instead OsloCTM3 uses pre-calculated 3-hourly meteorological fields from ECMWF forecasts produced daily with a 12-hourly spin-up starting from ERA-Interim reanalysis. All numerical model simulations assume the eruption starts on 31 August 2014 and that the Holuhraun SO₂ emissions are distributed equally in the vertical for grid cells between 0.8 and 3 km in the column containing the eruption vent following the magnitude and altitude profile of emissions described in Malavelle et al. (2017). All models include additional background SO₂ emissions from anthropogenic and natural sources. The simulations are continued from multiyear control simulations. All model output is regridded to a common regular 1.0^∘ × 1.0^∘ latitude–longitude grid using linear interpolation. Details specific to individual numerical models and key references can be found in Table 2.

2.4 Backward trajectories

Lagrangian modelling has been used previously to study the long-range transport of Holuhraun pollutants (e.g. Boichu et al., 2019; Schmidt et al., 2015). Here the origin, age, and travel distance of air parcels associated with sulfurous pollution events detected in the EMEP network observations are estimated using backward single-particle trajectories generated by the Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model developed by the National Ocean and Atmospheric Administration (NOAA) Air Resources Laboratory (Stein et al., 2015). We use 6-hourly ERA-Interim reanalysis data interpolated to an hourly resolution and regridded onto a 1.0^∘ × 1.0^∘ latitude–longitude grid as the meteorological input to the HYSPLIT model, a choice made to keep the driving meteorological data consistent between the trajectory analysis, and the nudging of the GCMs and CTM. Beginning on 1 September 2014 00:00 UTC, every hour at each EMEP station a new 27-member ensemble of 10 d backward trajectories is initiated at the coordinates and starting heights listed in Table 1 until 31 October 2014 (a total of 1464 ensembles for each station), with locations along each trajectory saved hourly. The ERA-Interim reanalysis data for each ensemble are offset by a fixed grid factor, a maximum of 1.0^∘ of latitude/longitude in the horizontal and 0.01σ units in the vertical, and so all the possible meteorological offsets result in the 27 members within each ensemble. As our ensembles are initiated at the beginning of every hour, pollution events observed at a daily resolution with a sampling midpoint centred on the hour have 25 ensembles available for analysis (12 h either side of the midpoint and the midpoint itself), whilst those events observed hourly with a sampling midpoint centred 30 min past the hour will have two ensembles (the bounding hours of the midpoint). This equates to a total of 675 and 54 individual trajectories, respectively, to evaluate the pollutant transport of each event.

A limiting factor with using a backward trajectory analysis is that the trajectories are not expected to arrive exactly at the eruption vent. Subsequently, a domain must be defined that sets the bounds as to whether a trajectory is deemed close enough to be attributed to the volcano. Defining this domain can be done visually through satellite imagery (e.g. Pardini et al., 2017) or by a statistical analysis (e.g. Hughes at al., 2012). Here we adopt the latter, defining multiple 3-D bounding cylinders centred on 64.85^∘ N, 16.83^∘ W with a height above mean sea level of 4.5 km and various bounding radii (see dotted regions in Fig. 1). The cylinder height is based on the maximum plume altitudes within the literature (Arason et al., 2015; Carboni et al., 2019a; Flower and Kahn, 2020; Pfeffer et al., 2018), whilst bounding radii are dependent on the distance from the eruption that the trajectories are initiated at. Trajectories released from stations distanced 1200–1500, 1500–1800, 1800–2100, and 2100–2400 km from Holuhraun are subject to radii of 320, 380, 440, and 500 km, respectively. These values are based on the positional error of a trajectory being approximately 10 %–30 % of the total distance travelled (Stohl, 1998). A special case is made for the Irafoss station due to its close proximity to Holuhraun (∼ 200 km). Over this distance, due to the finite hourly resolution of the trajectories and magnitude of local wind speeds, a trajectory is likely to travel further in a single time step than the estimate obtained from the “10 %–30 % distance travelled” method. Hence, even an ideal trajectory passing directly over the eruption may not be outputted within the bounding radius. Consequently, to ensure a near-100 % likelihood these ideal trajectories are captured, we define the Irafoss radius as the 99th percentile of the September and October ERA-Interim reanalysis horizontal wind speeds of the grid cells containing the horizontal eruption location and with midpoints below the cylinder height (see Fig. S1 in the Supplement).

Figure 1 Map of the 25 EMEP stations explored in this study. Stations marked with filled triangles experienced at least one sulfurous pollution event between 1 September and 31 October 2014 attributed to Holuhraun emissions, whereas stations marked with unfilled triangles did not. A red star indicates the location of the Holuhraun eruption (64.85^∘ N, 16.83^∘ W) with the surrounding dashed lines outlining the horizontal boundaries of the Holuhraun bounding areas defined in this study. From the inner circle outwards, the radii are 72, 320, 380, 440, and 500 km. Colouring links a station to the bounding area it is subject to.

This study attributes an observed sulfurous pollution event to volcanic emissions from Holuhraun when at least 25 % of the released trajectories pass through the relevant 3-D bounding cylinder, equating to a minimum of 167 and 14 trajectories for daily and hourly sampled events, respectively. Whilst this threshold could be considered low, allowing for other sources to contribute to the sulfurous pollution detected, the sheer volume of Holuhraun emissions within the region versus other sources during September and October 2014, and the rural location of the surface stations, gives confidence that this threshold is sufficient. For a given pollution event, we average the transport time and travel distance of the individual trajectories attributed to Holuhraun at their point of closest approach to the eruption to estimate the age and distance travelled by the plume. The error in the plume age is estimated as the larger value of either the standard error of the trajectories sampled or the trajectory temporal resolution (1 h).

Like all frameworks based on single-particles trajectories, our trajectory analysis is subject to the inherent uncertainty associated with individual trajectories (Stohl, 1998), with the uncertainty in the input meteorology often regarded as the dominant contribution (Engström and Magnusson, 2009; Gebhart et al., 2005; Harris et al., 2005). Here we minimise this main cause of uncertainty by perturbing the meteorology for each ensemble member. However, other uncertainties, such as the choice of meteorological dataset and/or trajectory model, and the exclusion of turbulence, have not been accounted for in our trajectory framework. Although these uncertainties are relevant, the focus of this study is to inter-compare numerical models consistently rather than through a rigorous dispersion exercise, and so they will not be considered further. Despite the simplicity of our single-particle trajectory framework, similar methods have been applied successfully to surface monitoring stations (e.g. Nieminen et al., 2015; Räty et al., 2023; Väänänen et al., 2013).

3 SO₂ plume spatial distribution

The SO₂ column load from IASI retrievals and model simulations on 23 September 2014 AM are displayed in Fig. 2, and an animation spanning September and October 2014 is given in Fig. S2. We mask IASI SO₂ column load retrievals below 0.5 DU to identify the observed horizontal extent of the plume. This threshold sufficiently exceeds the region's typical background SO₂ column load of approximately 0.1 DU, a background derived using the 2007–2009 September mean SO₂ mass burden of a similar geographical region given in Schmidt et al. (2015). Applying such a threshold ensures enough SO₂ from other sources is screened out whilst not removing data associated with the main volcanic plume. To enable a comparison with the IASI retrievals, the model output is regridded onto the same 1.0^∘ × 1.0^∘ latitude–longitude grid and is sampled within the observed plume extent at 09:30 UTC and 21:30 UTC for the AM and PM grouped retrievals, respectively. Areas within the observed plume extent that the models fail to capture (i.e. values below 0.5 DU) are shown in a distinct blue.

Figure 2 SO₂ column load from (a) Infrared Atmospheric Sounding Interferometer (IASI) retrievals and (b–g) model simulations of the Holuhraun eruption on 23 September 2014 AM. IASI data below 0.5 DU (Dobson units) are masked to identify the observed plume extent. Model output is sampled at 09:30 UTC. Non-hatched areas show model output within the observed plume extent, with the distinct blue highlighting model values below the 0.5 DU threshold. Hatched areas show model output outside the observed plume extent that is above the 0.5 DU threshold. Both IASI and model data have been regridded onto the same 1.0^∘ × 1.0^∘ latitude–longitude grid. Dashed red areas represent the Holuhraun region defined in this study. An animation of September and October 2014 is given in Fig. S2.

Within the Holuhraun region (dashed red area) in Fig. 2, UKESM1.0, HadGEM3-GA7.0 and OsloCTM3 perform well in capturing the observed plume extent with minimal blue areas present, whilst the ECHAM variants and MIROC6.1-SPRINTARS roughly capture 50 % and 30 %, respectively. These performances largely hold true across September and October as evident in the animation. Although, due to the binary nature of the observed plume extent masking, no magnitude on how far the modelled values lie below 0.5 DU is given, and so this metric should not be considered a sole indicator of model performance. In addition, we sample the models outside the observed plume extent when the modelled SO₂ column load exceeds 0.5 DU (hatched areas). In Fig. 2 all models except MIROC6.1-SPRINTARS simulate the plume outside the observed area over parts of western Europe. These areas potentially arise due to IASI retrieval limitations causing parts of the plume to be missed (e.g. cloud cover, high latitude, swath width) or due to high levels of background SO₂ emissions in the models (e.g. volcanic activity from Mt. Etna, anthropogenic activity). In September, for UKESM1.0, HadGEM3-GA7.0, and OsloCTM3 the hatched areas tend to dominate the blue areas, suggesting that their modelled plume areas are greater than the IASI retrievals, whereas the opposite is true for the remaining models. In October, all models largely show a greater modelled plume area than observed, yet this is partly due to the low IASI coverage across this period.

Overall, Fig. 2 and the animation make it apparent that the Holuhraun eruption is observed and modelled as the main source of SO₂ in our region of interest. Contributions from high background SO₂ sources are minimal relative to the total regional SO₂ and occur either outside or just within the outer bounds of the region and so are unlikely to substantially influence this study. Both visualisations show that the models capture the general features of the observed plume, particularly the dispersion over the Fennoscandian Peninsula and the UK during September, suggesting that nudging the models to ERA-Interim reanalyses gives credence to the models' ability to accurately simulate the plume dispersion despite their coarse resolution. It is worth noting that the animation shows possible regridding artefacts in the IASI retrievals on 11, 12, and 17 September PM. These artefacts occur outside our region of interest and so will not be considered further here.

Table 3 Monthly IASI-retrieved and IASI-simulated SO₂ plume heights and SO₂ mass burdens across the Holuhraun region (44 to 80^∘ N, 60^∘ W to 30^∘ E) for September (S) and October (O) 2014. Plume heights for MIROC6.1-SPRINTARS and ECHAM6.3-HAM2.3-P3 are derived from monthly resolution sampling, rather than bidaily as is done for other model estimates. The OsloCTM3 simulation does not contain the required diagnostics for a plume height estimate.

4 SO₂ plume height and mass burden

The average observed and modelled SO₂ plume heights across the Holuhraun region for September and October 2014 are shown in Figs. 3a and S3, respectively, with monthly values provided in Table 3. The IASI retrievals show that the observed SO₂ plume height, specifically the central height of a Gaussian SO₂ vertical profile, exists primarily (∼ 75 %) between 0.8 and 2.5 km above mean sea level, very rarely exceeding 3 km, showing that the volcanic perturbation to the region is contained well within the lower troposphere. The modelled SO₂ plume vertical extents are determined from model output that has been subject to the same sampling used in Sect. 3 and an additional masking in the vertical of grid cells with SO₂ mass concentrations below 4.5 µg m⁻³, a threshold based on the clean air SO₂ concentration of 1 ppbv (roughly 3 µg m⁻³ at 2–3 km) given in Theys et al. (2013). From the remaining grid cells with sufficient concentrations in the observed plume extent, the modelled central SO₂ plume height, represented by a solid line, is calculated as the regional mean of the heights of the grid cells containing the maximum SO₂ concentration in each column. The associated error is based on the typical vertical resolution of the model between 2 and 4 km above mean sea level. The top and bottom of the modelled SO₂ plume, illustrated by the outer dotted lines enclosing the envelope, are the regional averages of the maximum and minimum heights of the sufficiently polluted grid cells. Figure 3c shows the underlying number of data points contributing to the regional means. Regional means calculated from a number of data points below the minimum threshold of 35 (dashed red line) are deemed inadequate for a fair comparison and do not contribute to the rolling means shown in Fig. 3b. Note that from the output provided to this experiment, ECHAM6.3-HAM2.3-P3 and MIROC6.1-SPRINTARS SO₂ plume heights can only be compared with IASI retrievals at a monthly resolution, whereas no comparison is possible for OsloCTM3.

Figure 3 (a) Bidaily temporal evolution and (b) 5 d rolling mean of IASI-retrieved and IASI-modelled SO₂ plume heights across the Holuhraun region (44 to 80^∘ N, 60^∘ W to 30^∘ E) for September 2014. Black crosses and error bars represent the regional mean IASI SO₂ plume height and associated error. Models are sampled only within the observed plume extent. The coloured lines and error bars represent the regional mean central height of the modelled SO₂ plume and associated error, whilst an envelope indicates the regional mean height of the top and bottom of the modelled SO₂ plume (see text). The 5 d rolling means tolerate a maximum of two missing data points. (c) Number of grid cells used to calculate the regional means. The dashed red line indicates the minimum number of grid cells deemed sufficient to enable a fair comparison between observations and the model. Both IASI and model data have been regridded onto the same 1.0^∘ × 1.0^∘ latitude–longitude grid.

Generally, UKESM1.0, HadGEM3-GA7.0, and ECHAM6.3-HAM2.3 overestimate the central plume height across September, as clearly evidenced in the 5 d rolling mean model-to-observed ratio. This overestimation is greatest during the third week, particularly on 16 and 18 September. This feature may be a consequence of the vertical winds in the models not being constrained and/or additional variability in the momentum flux during the eruption, which is not accounted for in the prescribed emission profile used in the models. The observed variability in the plume height is largely well represented in the three models with all peaks, aside from 11 September PM, captured within error, whilst the observed central height is very rarely found outside the modelled vertical extents. The bottoms of the modelled vertical extents are close to the surface, suggesting that ground-based stations within the region are likely to experience moments of sulfurous pollution due to the eruption. The performance of the three models during October is similar (see Fig. S3), although the limited IASI retrievals make comparison harder. On a monthly scale, Table 3 indicates that all models where comparison is possible agree with the mean September and October observed heights within error, providing confidence that the models adequately capture the plume height within the Holuhraun region at this temporal resolution.

The IASI-retrieved and IASI-modelled SO₂ mass burdens across the Holuhraun region for September and October 2014 are shown in Figs. 4a and S4, respectively, with monthly values provided in Table 3. Both the observed and modelled mass burdens are derived by summing the product of the SO₂ column load and surface area of the individual grid cells within the observed plume extent for each bidaily interval. The same method is applied to the SO₂ column load IASI retrieval error to estimate an SO₂ mass burden observational error. For September and October 2014, we estimate an average bidaily SO₂ mass burden of 52 ± 15 kt of SO₂ and 34 ± 9 kt of SO₂ from IASI retrievals, respectively, which is in excellent agreement with Malavelle et al. (2017), who report corresponding mass burdens of 52 kt of SO₂ and 30 kt of SO₂ using an independent IASI dataset for the same geographical region. September mass burdens derived in Schmidt et al. (2015) using the Ozone Monitoring Instrument (OMI) are considerably higher, averaging 99 ± 49 kt of SO₂, across a slightly smaller area (45 to 75^∘ N, 60^∘ W to 30^∘ E). There exists substantial bidaily variation in the observed mass burden evident by the peaks of 180 ± 60, 130 ± 30, and 130 ± 30 kt of SO₂ on 5, 20, and 23 September, respectively, and the low values below 15 kt of SO₂ (e.g. 13–16 September). This variation is likely caused by a combination of the plume passing in and out of the defined region, changing IASI retrieval coverage (see Fig. 3c), and fluctuations in the volcanic SO₂ emission flux (Thordarson and Hartley, 2015).

Figure 4 (a) Bidaily temporal evolution and (b) 5 d rolling mean of IASI-retrieved and IASI-modelled SO₂ mass burdens across the Holuhraun region (44 to 80^∘ N, 60^∘ W to 30^∘ E) for September 2014. Black crosses and error bars represent the mass burdens and associated error derived from the IASI retrievals. Models are sampled only within the observed plume extent, and their derived mass burdens are given by the coloured lines. The 5 d rolling means tolerate a maximum of two missing data points. Both IASI and model data have been regridded onto the same 1.0^∘ × 1.0^∘ latitude–longitude grid. The same grid cells used in Fig. 3 are used to derive the mass burdens shown here.

With respect to the models, HadGEM3-GA7.0, ECHAM6.3-HAM2.3, and OsloCTM3 simulate average bidaily SO₂ mass burdens that lie close to those of the IASI retrievals for September and October 2014. UKESM1.0 overestimates the observed mass burdens, particularly during the early stages of September, which is potentially due to overpredicting total column SO₂, a bias that has been noted previously (Hardacre et al., 2021). As the IASI instrumentation is not able to sample the full intricacies of the plume, the observed SO₂ mass burden presented here is to be considered a lower estimate, and so UKESM1.0 exceeding this total may not necessarily be an indicator of poor performance. ECHAM6.3-HAM2.3-P3 and MIROC6.1-SPRINTARS largely underestimate the IASI-derived mass burdens, yet as the models are only sampled within the observed plume extent, these two models, as well as the remaining models, may simulate considerable mass in regions outside this extent (hatched areas in Fig. 2 and animation in Fig. S2). All models capture the observed variability, simulating larger mass burdens during early September when the eruption is most powerful and prescribed emission rates the highest, before decreasing during October. Correcting the IASI retrievals for parts of the SO₂ plume potentially missing has proved valuable (e.g. Carboni et al., 2019a) and could improve the comparison of the modelled heights and mass burdens presented here, yet as the general variability of both characteristics is well captured and no significant defects exist, using a cloud-adjusted correction is deemed unnecessary here.

5 Sulfurous surface mass concentrations

Observed time series of surface mass concentrations of total sulfur content, SO₂, and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ from 1 September to 31 October 2014 at EMEP stations Anholt (Denmark) and Kårvatn (Norway) are shown in Fig. 5a–c and d–f, respectively. Time series of the remaining EMEP stations are provided in Figs. S5–S27. Anholt and Kårvatn feature numerous peaks in sulfurous concentrations that exceed the climatological monthly statistics suggesting that these concentrations are significant and, given the rural locations of the sites, are likely caused by far afield sources of pollutants. Using the definition given in Sect. 2.2, we find that Anholt experienced 10 pollution events during September and October whilst Kårvatn experienced 14. We see merit in defining a pollution event using the total sulfur content concentration, rather than the commonly used SO₂ concentration (e.g. Boichu et al., 2019), as additional events are identified due to their high ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations which otherwise would have been missed (e.g. 20 and 23 September at Anholt, 12 September at Kårvatn). The number of sulfurous pollution events observed across the 25 EMEP stations during September and October 2014 is shown in Fig. 5g. Birkenes II (Norway), Irafoss (Iceland), Pallas Matorova (Finland), and Risoe (Denmark) all experienced roughly an event every 2 d, whilst only Leba (Poland), Preila (Lithuania), and Rucava (Latvia) did not experience any pollution episodes. In total, 283 pollution events are observed at 22 EMEP stations, indicating that widespread sulfurous pollution occurred across Europe in the months following the eruption.

Figure 5 Observed surface mass concentration time series of total sulfur content, SO₂, and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ at (a–c) Anholt (Denmark) and (d–f) Kårvatn (Norway) between 1 September and 31 October 2014. Climatological monthly mean and 90th percentile values are shown by the solid and dashed red lines, respectively. (g) Sulfurous pollution events (see Sect. 2.2) identified across the 25 EMEP stations highlighting those attributed to Holuhraun emissions (solid green filled) and those not (hashed orange filled).

The likelihood of Holuhraun being a main source of pollution for the 283 events can be established qualitatively using the IASI-retrieved and IASI-modelled SO₂ column load animations or more robustly using the trajectory framework outlined in Sect. 2.4. Using the latter approach, the main source of pollution for 111 (39.2 %) of the events can be attributed to Holuhraun emissions (see Fig. 5g for a station-by-station breakdown). Of the 22 EMEP stations experiencing a sulfurous pollution episode between September and October 2014, 19 stations endured at least one event influenced by the eruption. Note that other sources of pollutants may contribute to the mass concentrations observed at these 111 events, yet these contributions are likely minor given the rural setting of EMEP stations and that Holuhraun is the dominant sulfurous source in the region covering this period. None of the combined 17 events observed at Bredkälen (Sweden) and Tustervatn (Norway) are attributed to Holuhraun emissions which, given that the plume has been shown to pass this area (e.g. Grahn et al., 2015; Ialongo et al., 2015), suggests an inconsistency in the trajectory analysis. This inconsistency could be resolved by revising the heights the trajectories are released at, by incorporating additional meteorological datasets and/or trajectory models, or by using a more comprehensive trajectory framework (e.g. dispersion modelling). As assessing our trajectory framework is not the focus here, these inconsistencies have not been explored further. Nevertheless, the trajectory analysis shows that Holuhraun brought about significantly elevated sulfurous surface mass concentrations across Europe in September and October 2014, testament to the sheer volume of SO₂ emitted into the region by the eruption.

Our catalogue of 111 sulfurous pollution events attributed to Holuhraun emissions is used to assess model performance in capturing the surface level behaviour of the plume. Figure 6 displays the modelled versus observed surface mass concentrations of total sulfur (Fig. 6a–f), SO₂ (Fig. 6g–l), and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ (Fig. 6m–r) and the SO₂-to-${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ ratio (Fig. 6s–x) of the volcanically influenced episodes. Colouring is used to highlight the plume age at the time of sampling (see Sect. 2.4). In terms of reproducing the observed total sulfur content, MIROC6.1-SPRINTARS, HadGEM3-GA7.0, and UKESM1.0 perform well, with 76.6 %, 75.7 %, and 71.2 % of values within the 3-to-1 range, respectively, whilst the two ECHAM variants capture just below 60 %. These lower values are largely due to underestimating the higher sulfur content observed in the young plume (0–96 h), which itself is due to a considerable underestimation of the underlying SO₂ during these early stages of the plume. The SO₂ surface concentrations are also underestimated in MIROC6.1-SPRINTARS, yet for a slightly more mature plume (72–120 h), whilst the opposite is evident in UKESM1.0 and HadGEM3-GA7.0, predominantly for plume ages exceeding 96 h. The overestimation of European surface SO₂ concentrations in UKESM1.0 has been noted previously (Hardacre et al., 2021), albeit over longer timescales. OsloCTM3 simulates the SO₂ well, showing no apparent overprediction or underprediction. All models depict the observed decrease in SO₂ concentration with increasing plume age.

Figure 6 Modelled versus observed surface mass concentrations of total sulfur, SO₂, and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and the SO₂-to-${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ ratio for sulfurous pollution events attributed to Holuhraun emissions across the EMEP network for September and October 2014. Solid and dashed black lines represent parity and the 3-to-1 region, respectively, with counts of points within the latter given. Colouring illustrates the plume age at the time of sampling. Observational errors are too small to discern.

With respect to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, all models improve on their SO₂ performance, with each model having at least 61 % of simulated concentrations within the 3-to-1 range. The models show no obvious overestimation or underestimation, aside from UKESM1.0 and MIROC6.1-SPTRINARS where concentrations for plume ages above and below 96 h are generally underpredicted and overpredicted, respectively. This underprediction of surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ across Europe by UKESM1.0 has been stated previously (Hardacre et al., 2021; Mulcahy et al., 2020). In addition, the models struggle most in capturing the observed SO₂-to-${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ ratios, with only OsloCTM3 and UKESM1.0 having more than 50 % of simulated values within a 3-to-1 ratio of the observed. Broadly, all models, except for the ECHAM variants, simulate lower ratios than observed for the young plume (0–96 h) and higher ratios for the mature plume (144–216 h). The two ECHAM variants largely underestimate the observed ratio across all plume ages. There exists a notable underprediction in the ratio by MIROC6.1-SPTRINARS for plume ages roughly between 30 and 60 h. Both the observed and modelled ratios decrease with increasing plume age, suggesting that SO₂ oxidation to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ is occurring within the observed and simulated plumes.

Possible causes of differences between observed and simulated surface level behaviour of far-afield Holuhraun pollutants, such as vertical resolution, source emission profile, and sub-grid turbulence parameterisations, have been explored in depth previously (e.g. Boichu et al., 2016; Schmidt et al., 2015) and so will not be explored further here. As these challenges are not specific to Holuhraun and feature extensively in most numerical dispersion problems, the discrepancies that arise from them should not act as evidence against the use of these models for the ACI investigation in Part 2 of this study. In fact, given the relatively fine spatial and temporal resolution that these coarse models are being assessed against here, they perform commendably in capturing the surface level behaviour of the plume.

6 In-plume SO₂ oxidation to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$

We have demonstrated that Holuhraun emissions affected the troposphere over long distances, triggering SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ pollution events across Europe, and that the ratio of SO₂ to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ during these episodes decreases as the plume ages. This suggests that SO₂ oxidation to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ is occurring as the plume matures. There are two main pathways for this conversion in the troposphere: gas-phase reactions, largely with the hydroxyl radical (OH-), and aqueous-phase reactions with dissolved ozone (O₃) and hydrogen peroxide (H₂O₂) (e.g. Calvert et al., 1978; Stevenson et al., 2003). The ratio of SO₂ to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ is therefore useful in assessing whether oxidation processes are being accurately represented in the models; a ratio greater than that observed suggests that the overall oxidation processes are too slow, whilst a ratio less than that observed suggests that the overall oxidation processes are too fast. This assessment is carried out here using the 111 pollution events attributed to Holuhraun emissions. By focusing on the ratio of the two pollutants as opposed to the total sulfur content and assuming that volcanic SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ coexist, the variation in the absolute Holuhraun daily sulfurous emission flux can be ignored.

Table 4 Summary of the observed and modelled in-plume SO₂ oxidation to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ using the sulfurous pollution events attributed to Holuhraun emissions for September and October 2014.

Figure 7a shows the observed surface mass concentration ratio of SO₂ to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ on a logarithmic scale of the 111 sulfurous pollution events attributed to Holuhraun emissions versus the age of the plume at the time of sampling, with colouring highlighting the plume's travel distance from the eruption (see Sect. 2.4). A variety of plume ages and SO₂-to-${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ ratios make it apparent that the plume is sampled during different stages of maturity. The linear characteristics of Fig. 7a, along with the curve depicted in the equivalent linear scale figure in Fig. S28), imply an exponential decay of the SO₂-to-${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ ratio with plume age, a relationship commonly associated with the depletion of volcanic SO₂ and one that assumes first order kinetics with respect to the SO₂ concentrations (e.g. Ilyinskaya et al., 2021; McGonigle et al., 2004; Oppenheimer et al., 1998; Pattantyus et al., 2018). Consequently, here we fit the observed data to an exponential decay function using non-linear total least-squares regression as this considers both the ratio and plume age uncertainties. The derived SO₂-to-${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ oxidation rate constant is 0.032 ± 0.002 h⁻¹, corresponding to an e-folding time of 1.30 ± 0.08 d. Using IASI-retrieved SO₂ column load, Carboni et al. (2019a) estimate Holuhraun SO₂ depletion as having a mean 6-month e-folding time of 2.4 ± 0.6 d, whilst Schmidt et al. (2015) derive a mean September SO₂ e-folding time of 2.0 ± 0.8 d using NAME simulations of the eruption. Whilst not directly comparable, as these studies have not estimated the oxidation rate explicitly and focus on different time periods, both estimates are of similar magnitude to the SO₂ oxidation e-folding time found here. Assuming our exponential decay relationship holds close to the eruption vent, this study estimates a near-vent SO₂-to-${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ ratio of 25 ± 5. This result agrees with Ilyinskaya et al. (2017), who report SO₂-to-${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ ratios of 2 to 250 and 4 to 94 at 100 and 250 km from the vent, respectively. Boichu et al. (2019) estimate a slightly lower near-vent SO₂-to-${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ ratio of 19.7 using a linear model created from five observations to describe the evolution of the SO₂-to-${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ ratio.

Figure 7 (a) Observed and (b–g) modelled SO₂-to-${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ surface mass concentration ratios of sulfurous pollution events attributed to Holuhraun emissions with respect to the plume age at the time of sampling for September and October 2014. Red lines represent the exponential decay fits. Dashed black lines are the observed exponential fit overlaid onto the modelled ratios. Colouring illustrates the plume's travel distance at the time of sampling. Observational errors are too small to discern.

Figure 7b–g depict the modelled surface mass concentration ratio of SO₂ to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ versus the age of the plume that these ratios are sampled in. All models display an exponential relationship between the SO₂-to-${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ ratio and plume age, albeit only for the young plume (< 70 h) in the case of MIROC6.1-SPRINTARS. Each model is fitted to an exponential decay function, which is given by the solid red line with the observed fit overlaid by the dashed black line for comparison. Due to how the plume age error increases with plume age, and how the magnitude of the ratio error is negligible in comparison, ratios sampled in the mature plume have a larger total error and so are weighted less in the fitting than those sampled in the younger plume. The exponential decay parameter estimates for the models are displayed in Table 4. The modelled near-vent ratios are all smaller than those derived from observations, yet they still agree with those found in Ilyinskaya et al. (2017). Except for MIROC6.1-SPRINTARS, all model-derived oxidation rates are slower than those derived from observations, ranging from being roughly just under twice as slow in ECHAM6.3-HAM2.3-P3 and OsloCTM3 to 4.5 times as slow in UKESM1.0, HadGEM3-GA7.0, and ECHAM6.3-HAM2.3. The seemingly poor fit of MIROC6.1-SPRINTARS is due to the underestimation of the ratios between approximately 30 and 60 h (Fig. 6w), which, for reasons stated previously, have a larger influence on the fit than ratios sampled in the more mature plume. If ratios across this time range were better captured by MIROC6.1-SPRINTARS, a fit that results in a slower derived oxidation rate more in keeping with the other models and observations could be expected. Interestingly, there is no apparent correlation between a model's vertical resolution and a model's ability to capture the in-plume SO₂-to-${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ oxidation.

7 Summary and conclusions

By releasing 9.6–11.8 Mt of SO₂ into the lower troposphere across nearly 6 months, the 2014–2015 Holuhraun eruption offers an opportunity to challenge the capability of GCMs in capturing the characteristics of tropospheric sulfate aerosol intricacies resulting from effusive eruptions and assess the potential impact of subsequent aerosol–cloud interactions. A model inter-comparison effort has been initiated to leverage this opportunity, and the results from Part 1 of the two-part analysis are presented here. Remote sensing data of SO₂ and surface level SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ mass concentration measurements are used, in conjunction with trajectory modelling, to evaluate the performance of five GCMs and a CTM in simulating the spatial and chemical evolution of the SO₂ plume across the North Atlantic and Europe.

A comparison against IASI SO₂ retrievals shows that the models capture the evolution of the volcanic plume within the surrounding region well during September and October 2014. Holuhraun emissions are the dominant source of SO₂ in the models, and the spatial transport of the associated SO₂ plume is well replicated. The SO₂ plume height is slightly overestimated by the models, whereas there is no general overestimation or underestimation in simulating the SO₂ mass burdens; it is model dependent. The temporal variability of both these plume characteristics is well captured. Discrepancies with the IASI retrievals could be due to several factors including the limitations of the IASI retrievals (e.g. Carboni et al., 2019a) and discrepancies between the idealised volcanic emission profile used by the models and the real emissions (e.g. Steensen et al., 2016). A comparison against retrievals of volcanic SO₂ from other satellite instrumentations may yield different conclusions, yet the descriptions of the plume spatial distribution made with other remote sensing products are similar (e.g. OMI: Ialongo et al., 2015; Schmidt et al., 2015, Steensen et al., 2016; OMPS-NN: Ialongo et al., 2015; GOME-2: Twigg et al., 2016). Even though the model spatial representations of the eruption are not perfect, our intent here is rather to identify the variations in the models' transport of Holuhraun SO₂ as this will help discern the impact on cloud properties and assess the ACIs in Part 2.

By combining the surface mass concentration measurements of SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ made during September and October 2014 across the EMEP network with single-particle trajectories calculated using the HYSPLIT model, the simulated surface level behaviour of the plume was assessed. Of the 283 sulfurous pollution events identified, 111 are attributed to Holuhraun emissions. Generally, the models reproduce the measured elevated surface level concentrations during these volcanically influenced events, yet they struggle in simulating the correct magnitude, notably the ratio of SO₂ to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, which is often underestimated and overestimated for the young and mature plume, respectively. However, this should not be discouraging as capturing volcanic sulfurous pollutant surface mass concentrations far afield at a specific location and time is challenging, even for CTMs of finer scales. We note that the models with finer vertical resolutions, UKESM1.0, HadGEM3-GA7.0, and OsloCMT3, describe ground-level concentrations of Holuhraun pollution episodes best, a feature that has been noted previously by Boichu et al. (2016). Given the relatively coarse scale of the simulations discussed here, the surface level performance of the models is admirable.

Both the observed and modelled ratios of SO₂-to-${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ surface mass concentrations sampled within the plume are shown to decrease with increasing time and distance from the eruption vent, suggesting SO₂ oxidation to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ is occurring. To explore this further, the ratios as a function of plume age have been analysed, revealing an exponential decay. By fitting this decay to an exponential function, observed and modelled rate constants for the volcanic SO₂ oxidation are estimated. Aside from MIROC6.1-SPRINTARS, the in-plume SO₂-to-${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ oxidation is shown to be slower in the models than observed. This implies that the volcanic SO₂ introduced into the simulations may not be chemically converted fast enough relative to what is derived from surface measurements (if deposition effects are ignored). The considerable underestimation of the ratios sampled at plume ages between 30 and 60 h is the suggested reasoning for MIROC6.1-SPRINTARS exhibiting opposing behaviour. No correlation between a model's vertical resolution and a model's derived SO₂ oxidation rate is found.

The oxidation rate constants explored here are generalised values representing both the gas-phase and aqueous-phase pathways. This study attempted to help elucidate the complexity of volcanic SO₂ oxidation by fitting the SO₂-to-${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ ratios to a bi-exponential function, a sum of two individual exponential decay components, to distinguish between the two pathways by estimating individual gaseous and aqueous oxidation rate constants. Despite the success of previous studies in estimating and applying multiple rate constants to describe the depletion of volcanic SO₂ (Ilyinskaya et al., 2021; Pattantyus et al., 2018), this study found no significant improvement in the fitting function versus a standard exponential decay (one exponential decay component). Categorising the SO₂-to-${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ ratios in terms of the conditions the air parcels are subject to during transport, such as time spent in cloud, relative humidity, cloud pH, oxidant concentrations, time of day, and deposition rates, could reveal the dominant oxidation pathway/s affecting a particular group. The subsequent fitting versus plume age would then provide rate constant estimates of the mechanism/s in play. However, such a method would likely require a sophisticated Lagrangian framework, and it is beyond the scope of this work to explore the intricate chemical kinetics of volcanic SO₂ oxidation.

Overall, the six models considered here provide reasonable simulations of the spatial and chemical evolution of the Holuhraun plume and are considered competent enough to be used to explore the impacts of the eruption on ACIs in the region (see Part 2 of this study). It is important to acknowledge, and is possibly relevant to the wider ACI community, that this analysis has also highlighted that the models do not perfectly capture the secondary ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ aerosol production during a large degassing event. This pitfall may contribute to the well-documented disagreements between model ACIs estimates as the underlying aerosol perturbations could differ. We hope that our application of in situ sulfurous surface measurements to assess numerical models helps bolster the case to retain and extend air monitoring networks of volcanic pollutants for use in future studies.