Evaluation of atmospheric sulfur dioxide simulated with the EMAC (version 2.55) Chemistry–Climate Model using satellite and ground-based observations

Sulfur dioxide (SO₂) is a key atmospheric pollutant, primarily emitted through human activities such as fossil fuel combustion. In atmospheric models, accurate representation of SO₂ emission sources, transport, and removal processes are essential for evaluating air quality and radiative forcing.

In this study, we present, for the first time, a comprehensive examination of atmospheric SO₂ simulated by the ECHAM/MESSy Atmospheric Chemistry (EMAC) model, here operated under the Chemistry–Climate Model Initiative (CCMI-2022) protocol. First, the tropospheric sulfur budget simulated by EMAC is verified to be closed. This closure means that all sulfur sources and sinks are balanced and no artificial gain or loss occurs over time due to numerical or conceptual errors. This budget closure is a prerequisite for any further analysis. Second, the results of EMAC simulations are compared with observations from three ground-based networks (the Clean Air Status and Trends Network (CASTnet), the European Monitoring and Evaluation Program (EMEP), and the Acid Deposition Monitoring Network in East Asia (EANET)), mainly over polluted regions, and with vertical column densities retrieved from a TROPOspheric Monitoring Instrument (TROPOMI) on board the Copernicus Sentinel-5 Precursor mission (Sentinel-5P) satellite. The EMAC simulated SO₂ concentrations near the Earth's surface for the year 2019 are, depending on the region, between 1.4 and 1.8 times larger than observed. This discrepancy aligns well with the differences between simulated and retrieved satellite-based measurements of SO₂ vertical column densities over the same regions. It indicates that the prescribed SO₂ emissions used for the EMAC simulations might be overestimated. Over a longer time period (2000–2019), the EMAC simulation reproduces the measured declining trends of SO₂ concentrations and deposited sulfur fluxes in the USA and Europe, but fails to simulate the observed trends in East Asia. This is most likely attributable to the prescribed SO₂ emission inventories. Furthermore, sensitivity simulations are performed to assess the emitted amount of SO₂ following the Raikoke and Ulawun volcanic eruptions in 2019. The results show a very good agreement of the simulated temporal evolution of the amount of atmospheric SO₂ after the eruptions with that retrieved from satellite-based observations.

1 Introduction

Air pollution remains a significant global challenge, affecting both, human health and the Earth's climate (Wood et al., 2024; Arias et al., 2021). Among various pollutants, SO₂ plays a key role due to its strong influence on atmospheric chemistry, air quality, and climate processes (Seinfeld and Pandis, 2016; Myhre et al., 2013). Anthropogenic SO₂ emissions primarily originate from the combustion of sulfur-containing fossil fuels, oil refining, and metal smelting (Smith et al., 2011; Klimont et al., 2013), while natural sources include volcanic eruptions, the oxidation of dimethyl sulfide (DMS) emitted from the ocean, and minor biogenic contributions from land ecosystems (Lana et al., 2011; Fioletov et al., 2016; Quinn et al., 2011).

SO₂ is the dominant precursor of sulfate aerosols, which influence the Earth's radiation balance by scattering incoming solar radiation and acting as cloud condensation nuclei (CCN) (Charlson et al., 1992; Lohmann and Feichter, 2005). These processes contribute to short-term climate cooling, partially offsetting warming caused by greenhouse gases (Arias et al., 2021; Albrecht, 1989). At the same time, SO₂ contributes to adverse environmental effects such as acid deposition and impacts on the stratospheric ozone layer (Seinfeld and Pandis, 2016; Solomon, 1999).

While global SO₂ emissions have declined in many industrialized regions due to regulatory efforts (Klimont et al., 2013; Liu et al., 2018), emissions remain high in rapidly developing countries like India and China (Dahiya et al., 2020), despite recent reductions (Zhang et al., 2019). Furthermore, episodic volcanic eruptions introduce large amounts of SO₂ into the atmosphere, affecting its distribution on regional and global scales (Carn et al., 2017). Despite improvements in satellite monitoring and emission inventories, uncertainties remain regarding the atmospheric lifetime, transport, and transformation of SO₂.

Accurately simulating the complex processes governing SO₂ behavior in the atmosphere is essential for understanding its role for air quality, climate forcing, and environmental impacts such as acid rain. Chemistry–climate models (CCMs) like EMAC provide a comprehensive framework to represent emissions, chemical transformations, transport, and deposition of sulfur compounds within the coupled atmosphere system. EMAC, in particular, incorporates detailed tropospheric and stratospheric chemistry schemes, making it well suited to investigate the sulfur cycle from emissions to atmospheric sinks (Jöckel et al., 2010).

The current study investigates the distribution and budget of tropospheric SO₂ using the EMAC model and observational datasets. The study evaluates the model's ability to reproduce SO₂ spatial and temporal distributions by comparing model simulation results with observations retrieved from a satellite instrument and with ground-based measurements. Furthermore, this paper investigates the tropospheric sulfur chemistry within the employed EMAC model. This is done by examining the SO₂ emissions, the sulfur related chemical processes, and the sink processes (including wet and dry deposition, and sedimentation) to verify the model's ability to conserve sulfur mass. This conservation is a prerequisite for the comparative analysis and inter-comparison with results from other models and with observational data, thereby showing the numerically correct representation of chemical processes simulated within the model.

This paper includes the following sections: The used EMAC model setup is illustrated in Sect. 2. A detailed study of the tropospheric sulfur budget in the EMAC model is presented in Sect. 3. Section 4 shows the evaluation of the simulated global distribution of SO₂, as well as the variations of SO₂ following eruptive volcanic events, using data retrieved from the TROPOMI instrument. Section 5 assesses the comparison of simulated SO₂ concentrations and sulfur deposition fluxes at the Earth's surface with ground-based measurements. Last but not least, the Conclusions and an Outlook of this study are presented in Sect. 6.

2 Model description

2.1 The EMAC model

In the present study, a detailed investigation and evaluation of SO₂ simulated by the global EMAC model (Jöckel et al., 2016) integrated within the Modular Earth Submodel System (MESSy) framework (Jöckel et al., 2010), is undertaken. The EMAC model is a comprehensive global CCM that represents physical and chemical processes in the troposphere and middle atmosphere, along with their interactions with the land surface, ocean systems, and human-induced changes such as emissions and land-use (Jöckel et al., 2010, 2016). The second version of the Modular Earth Submodel System and the 5th generation European Center Hamburg general circulation model (ECHAM5, Röckner et al., 2006) make up the EMAC model used in this study. The physics-related submodels within the MESSy framework have been adapted from the physics routines of ECHAM5 (Jöckel et al., 2016). Only the spectral dynamical core, the flux form semi-Lagrangian (FFSL) large scale advection scheme (Lin and Rood, 1996), the time integration loop, and the Newtonian relaxation methods retain their original structure from the ECHAM5 base model.

The results analysed here stem from the RD1SD-base-01 EMAC simulation (Jöckel et al., 2024b) that has been performed under the CCMI-2022 protocol (CCMI, 2023). Here, the emissions of SO₂ and other sulfur species are calculated by the submodels OFFline EMISsions (OFFEMIS, formerly called OFFLEM) for prescribed emission fluxes (Kerkweg et al., 2006b), and AIRSEA (calculating the air-sea exchange of chemical species Pozzer et al., 2006). Details about the emission setup used here are described in a separate subsection (Sect. 2.2).

Chemical reactions in the gas phase are computed by the submodel Module Efficiently Calculating the Chemistry of the Atmosphere (MECCA, Sander et al., 2019), while the Scavenging Submodel for Regional and Global Atmospheric Chemistry Modeling (SCAV) simulates the aqueous phase kinetics and scavenging processes in the atmosphere (Tost et al., 2006). The dry deposition of gases and aerosols is calculated by the dry deposition submodel (DDEP, formerly called DRYDEP) (Kerkweg et al., 2006a), and aerosol sedimentation is calculated by the aerosol sedimentation submodel (SEDI, Kerkweg et al., 2006a). The sampling along sun-synchronous satellite orbits submodel (SORBIT) to sample model results on-line along orbits of sun-synchronously orbiting satellites, as described by Jöckel et al. (2010), has been applied to facilitate a direct comparison between simulated Vertical Column Density (VCD) of trace gases such as SO₂, with observations from satellite instruments.

The RD1SD-base-01 EMAC simulation results analysed in this study cover the years 1979 to 2019. The simulation was performed at a resolution of T42L90MA with output of results every 5 h of the simulated period. The spectral resolution (triangular truncation) T42 corresponds to a quadratic Gaussian horizontal grid of roughly 2.8° × 2.8° in both, longitude and latitude coordinates, and L90 denotes 90 vertical layers (with a median lowest level height of 60 m) between the surface and the uppermost model layer centered around 0.01 hPa (Jöckel et al., 2010). For the RD1SD-base-01 simulation, the gas phase chemistry is calculated throughout the entire atmosphere using the Mainz Isoprene Mechanism (MIM1) based on Pöschl et al. (2000). This mechanism accounts for hydrocarbons up to 4 carbon atoms, along with isoprene (5 carbon atoms). SO₂ in the gas phase is oxidized by the hydroxyl-radical OH or direclty photolysed. The major sink, however, is by transition into the aqueous phase, mainly cloud water, and further oxidation in the liquid phase. Since gaseous SO₂ is not released on evaporation of cloud and rain droplets, the sulfur contents is in these cases transferred into a so-called residual (res) pseudo-aerosol tracer ${{\mathrm{SO}}_{\mathrm{4}}}_{\mathrm{res},\mathrm{cs}}$ with characteristics of a coarse mode soluble (cs) aerosol. This tracer is treated as aerosol tracer, e.g. by the sedimentaion submodel (SEDI), by the dry deposition submodel (DDEP), and by the wet scavenging submodel (SCAV). The details of the gas phase and liquid phase chemistry are documented as supplemantary material. Other aerosol effects were prescribed in both, the troposphere and the stratosphere, to consider their impact through heterogeneous chemistry and radiative forcing (Jöckel et al., 2016). To allow for a direct comparison of the simulation results, in particular chemical tracers, between the simulated and observational data, the RD1SD-base-01 simulation was operated in “specified dynamics” (SD) mode, for which the prognostic variables temperature, divergence, vorticity and the logarithm of surface pressure were “nudged” by Newtonian relaxation towards the fifth generation of European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis data (ERA5; Hersbach et al., 2020). The model dynamics of the SD simulations are then aligned with the observed dynamics, aiming a good reproduction of real meteorological situations.

2.2 Description of the used sulfur emissions

Sulfur emissions of both, anthropogenic and natural sources, need to be taken into account. The following prescribed emission inventories were used:

• Throughout this study the Coupled Model Intercomparison Project Phase 6 (CMIP6) inventory is selected as the standard inventory for global EMAC simulations, since it was recommended by the experimental protocol for participation in the CCMI-2022 model intercomparison initiative (Eyring et al., 2016). The CMIP6 inventory has a horizontal resolution of 0.5° × 0.5° and it primarily combines bottom-up inventories to provide emission data for climate models. Bottom-up inventories involve estimating emissions based on detailed data about specific sources and activities, such as energy consumption and industrial processes. The CMIP6 inventory contains historical emissions from 1850 to 2014, provided by the Atmospheric Chemistry and Climate-Model Intercomparison Project (ACCMIP) developed by Lamarque et al. (2010). The historical data are then combined with the shared socioeconomic pathways (SSPs) for projected future emissions from the IPCC Sixth Assessment Report (AR6, Calvin et al., 2023). The SSPs used within the CMIP6 inventory provide a range of future scenarios based on varying levels of greenhouse gas emissions and societal changes, such as SSP1-1.9, SSP1-2.6, SSP2-4.5, SSP3-7, SSP4-6, and SSP5-8.5 (Riahi et al., 2017). These SSPs present different emission scenarios, in order to explore different future climate outcomes based on varying levels of greenhouse gas emissions and societal changes. For the RD1SD-base-01 simulation the SSP2-4.5 scenario has been used to prescribe trace gases emissions, including SO₂ emissions after 2014. The SSP2-4.5 is a middle-of-the-road scenario with moderate emissions, leading to a radiative forcing of 4.5 W m² by 2100 (Riahi et al., 2017).

• The terrestrial Dimethyl Sulfide (DMS_terrestrial) emissions are based on the global inventory developed by Spiro et al. (1992). This inventory was mainly developed to examine gaseous sulfur emissions. Over the years, this inventory has been evaluated by other studies, such as Chin et al. (2000); Vallina and Simó (2007), and Lana et al. (2011). DMS_terrestrial emissions originate from both, vegetation and soils, and are available as a monthly resolved annual climatology at a resolution of 1° × 1° (Bates et al., 1987).

• Volcanic sulfur emissions from both, continuously degassing and explosive volcanoes are represented by an inventory of the Aerosol Inter Comparison (AeroCom) project as a zonal mean climatology (Dentener et al., 2006). Volcanic sulfur is emitted as 97.5 % SO₂ and 2.5 % SO₄. The data are based on the bottom-up Global Emissions Inventory Activity (GEIA) for the years 1750 and 2000 (Andres and Kasgnoc, 1998). Continuously degassing sulfur in the AeroCom inventory is equally distributed over the grid points with GEIA volcano locations and amounts to a multi-annual total emission of 12.6 Tg (S) yr⁻¹ (Teragrams of Sulfur per year) over all the years (Dentener et al., 2006). The height of these emissions is defined in the upper third of the volcano altitudes, simulating the degassing processes that occur predominantly at the volcano flanks. Explosive volcanic emissions are quantified at approximately 2 Tg (S) yr⁻¹ over all the years. This estimation is based on the Aerosol Index (AI) provided by the Total Ozone Mapping Spectrometer (TOMS) satellite sensors (Dentener et al., 2006). The emissions data are distributed evenly across grid boxes that include volcanoes, which were active in the last century (Halmer et al., 2002). It is important to note that these emissions are treated as being continuously released rather than episodic, due the fact that only about one-third of such emissions occur during violent explosive events (Dentener et al., 2006). Furthermore, these emissions are typically defined to occur between 500 and 1500 m above the peaks of the volcanoes, to accurately represent their dispersal in the atmosphere. The injection height, time resolution and the sulfur flux of the different volcano types are listed in Table 1.

Table 1 Parameters of the AeroCom explosive and continuous volcanic emissions. V_T (volcano top) corresponds to the altitude of the top of the volcano.

In addition to the prescribed emissions, sulfur from oceanic Dimethyl Sulfide (DMS_airsea) and from Carbonyl Sulfide OCS are calculated using the submodels AIRSEA (for gas exchange between air and sea) (Pozzer et al., 2006) and TNUDGE (Kerkweg et al., 2006b) for Newtonian relaxation towards prescribed mixing ratios, respectively.

Since the concentration simulated by the model is affected by the prescribed emissions, it is important to understand the differences between the used CMIP6 emission inventory and other emission inventories. In this study, the Emissions Database For Global Atmospheric Research (EDGAR) emission inventory (Solazzo et al., 2021) is used for this comparison (See Sect. 5). Same as for the CMIP6 inventory, EDGAR is also considered a bottom-up inventory, but with a finer horizontal resolution on grid-maps at 0.1° × 0.1°. EDGAR is developed using a bottom-up approach combining internationally available statistics on activity data with emission factors derived from scientific literature and guidelines (e.g., IPCC) (Crippa et al., 2019). The data is available as yearly and monthly mean and is emitted into 7 vertical tropospheric levels (0, 20, 92, 184, 324, 522 and 781 m), as described by Bieser et al. (2011). Version 5.0 of EDGAR (EDGAR5) contains solely historical data about anthropogenic emissions from different sectors such as fossil fuels, agricultural waste burning, ships and roads, starting from 1970 till present. Other emissions from large scale biomass burning, forest fires and sources from land-use forestry are excluded (Crippa et al., 2019). Here, the EDGAR 5.0 and the EDGAR8 inventories are used for inter-comparison to provide an estimate of the uncertainty of the magnitude and variability of SO₂ emissions over time from the CMIP6 emission inventory.

3 Tropospheric sulfur budget simulated with EMAC

This section provides a comprehensive evaluation of the global tropospheric sulfur budget in the RD1SD-base-01 EMAC simulation by examining prescribed sulfur emissions and the removal of sulfur-containing species via deposition over the years 2010 to 2019. These years were selected based on the availability of corresponding observational datasets used for later evaluation. The goal is to verify the internal consistency of the model's sulfur budget: the sulfur emitted into the atmosphere must either remain in the atmosphere (as part of the sulfur burden) or be removed through deposition processes. This ensures that the model conserves mass and accurately represents the sulfur cycle. The principle can be formulated as follows for each year:

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }B\left(t\right)=E\left(t\right)-D\left(t\right)\end{array}$$

Here, $\mathrm{\Delta }B\left(t\right)=B\left(t_{\text{end-of-the-year}}\right)-B\left(t_{\text{start-of-the-year}}\right)$ is the annual change in the atmospheric sulfur burden (in units of mass), E(t) is the total sulfur emission over the year, and D(t) is the total sulfur deposition over the same year. All quantities are integrated over the year. In this context, the burden B(t) represents the total mass of sulfur in the atmosphere (summed over all sulfur-containing species in the model domain) at a given time t. The difference ΔB(t) reflects the net accumulation (or loss) of sulfur in the atmosphere over the year.

Prescribed sulfur emissions, as applied in EMAC, arise from both, anthropogenic and natural sources. In the present study and for the year 2010, fossil fuel consumption, DMS from the ocean (denoted as DMS_airsea), DMS from terrestrial sources (denoted as DMS_Terrestrial), volcanic activity, and maritime shipping collectively contribute to nearly 95 % of the sulfur emissions released into the EMAC model atmosphere. Other sources, such as OCS, agricultural waste burning, and road emissions, constitute the remaining 5 % of the emitted sulfur.

Table 2 Detailed list of the emitted and deposited sulfur species for the year 2010 in the EMAC model. The first column represents the sulfur emission sectors and the third column (Tracers) shows the sulfur species deposited within the EMAC model. Suffixes _cs and _l denote species in coarse mode aerosol and liquid phase, respectively.

The released sulfur from these sectors becomes oxidized and is removed from the atmosphere through dry deposition, sedimentation (of sulfuric particles), and wet deposition/scavenging, which rinses sulfur through convective and large scale precipitation (cv + ls). The emitted, deposited and remaining sulfur species are exemplarily examined for the year 2010 in Table 2.

Table 3 Sulfur budget closure in Tg (S) yr⁻¹ between 2010 and 2019 in the EMAC model.

Table 4 A comparison of emission/deposition fluxes in Tg (S) yr⁻¹ between EMAC, STOCHEM-Ed model and IPCC AR3 for the year 1990.

Following the same analysis done for the year 2010, an evaluation of the sulfur budget for the years between 2010 and 2019 is shown in Table 3. Here, the tropospheric sulfur budget is nearly perfectly closed with a value near 0 (sulfur deficit) for the other years as well. In other words, this shows that the tropospheric sulfur budget is effectively balanced, accounting for the contributions of various sulfur species and their interactions over the specified time period.

Following the examination of the sulfur budget closure in the EMAC model, a comparison with other atmospheric chemistry models and relevant literature data helps assessing the model performance and identification of potential areas for improvement (see next Sect. 3.1).

3.1 How does the tropospheric sulfur budget in the EMAC model compare to that of other atmospheric chemistry models?

In their studies, Stevenson et al. (2003) and Penner et al. (2001) address the tropospheric sulfur budget within their atmospheric chemistry models for the year 1990. Stevenson et al. (2003) used the STOCHEM-Ed model, which is a global three-dimensional Lagrangian Chemistry-Transport Model (CTM), while the IPCC Third Assessment Report (IPCC AR3, Penner et al., 2001) incorporates the average results from 11 distinct models providing a comprehensive overview of the global tropospheric sulfur budget. Since the RD1SD-base-01 EMAC simulation covers the years 1979 to 2019, we directly compare the year 1990, which is used in both references. Note that in Stevenson et al. (2003), sulfur anthropogenic emissions were presented as a single sector, unlike the detailed breakdown presented in our study, which categorizes the prescribed emissions into sectors such as fossil fuels, aircraft and shipping emissions, among others. Additionally, the models in both studies applied an interactive aerosol model, unlike the EMAC model for the RD1SD-base-01 simulation. Therefore, some sulfur species present in Stevenson et al. (2003) have not been considered in the EMAC model, such as the Methane Sulphonic Acid (MSA) and vice versa (such as OCS).

For a clear comparison of the processes used by Stevenson et al. (2003) with those in EMAC, Table 4 provides the values for both, the prescribed emission and the deposition rates, for the year 1990. Despite discrepancies arising from the applied emission inventories between EMAC and the literature, as well as differences in used chemical reactions and physical processes, this comparative analysis reveals a good agreement between our model results and existing literature.

Consequently, the magnitudes of sulfur emissions across different sectors are consistent between EMAC and the compared models. The total emissions from EMAC align perfectly with those reported in the IPCC AR3, with a slight difference of about 1 Tg (S) yr⁻¹ in EMAC. Regarding the STOCHEM-Ed model, the EMAC simulation utilises larger sulfur emissions by approximately 12 %, related to the applied emission inventories. This overall good agreement in magnitude is also apparent in the deposition rates. In EMAC, the deposited sulfur through wet and dry processes is 93.79 Tg (S) yr⁻¹ compared to 95.3 and 96.5 Tg (S) yr⁻¹ for the STOCHEM-Ed model and the IPCC AR3, respectively. The remaining deposited sulfur from sedimentation is present only in EMAC and not in the cited literature.

This evaluation highlights the validity of the EMAC model in comparably capturing sulfur emission and deposition rates, despite minor differences in the chemical species and mechanisms used by the compared models. This discrepancy is attributed to the sulfur emissions inventories utilized as input for the models, which play an important role in shaping both, sulfur mixing ratios and deposition processes. Further refinement and validation of these emissions inventories may help to improve the models performance in representing the tropospheric sulfur cycle.

4 Comparison of SO₂ simulated by the EMAC Model with TROPOMI data

The primary objective of this section is to conduct a comprehensive global-scale comparison between the EMAC simulations and retrievals based on data from the TROPOMI instrument on board the Sentinel-5P satellite, which is the first Copernicus mission specifically designed for atmospheric monitoring, as mentioned by ESA (2017).

TROPOMI is performing atmospheric measurements, particularly for the quantification of VCDs of various gases and aerosols (Veefkind et al., 2012). These include ozone, formaldehyde, nitrogen dioxide, carbon monoxide, methane, aerosols, and sulfur dioxide, which is of specific importance for this study. Additionally, TROPOMI, notable for being the first satellite instrument to measure SO₂ columns with the highest spatial resolution (3.5 km by 7 km) among other satellites, serves as a pivotal dataset for this study. The analysis is focused on the year 2019, chosen for its status as, at the time of our study, being the first complete year of data available in both datasets (model and observations).

Performing a direct comparison between the EMAC model results and the TROPOMI/Sentinel-5P datasets presents inherent challenges due to their different methodologies and resolutions. Therefore, Appendix A2 describes the methodologies adopted in this study to enable a fair and detailed comparison of SO₂ between both datasets. To compare with the TROPOMI-retrieved VCDs, EMAC's VCDs are calculated by applying the standard/polluted averaging kernel (for the so called “standard case”) labeled “AK_polluted” from the Covariance-Based Retrieval Algorithm (COBRA; Theys et al., 2021). The COBRA product represents an advanced SO₂ retrieval technique from TROPOMI. Its improved sensitivity to low SO₂ levels enables the detection of previously undetected emission sources, including weakly emitting volcanoes and power plants (Theys et al., 2021).

This section is divided into three subsections. A comparison of the EMAC results on the global scale, using the standard (“polluted”) AK, against TROPOMI retrievals for the standard case is described in Sect. 4.1. The effects of eruptive volcanoes on the simulated SO₂ total column are discussed in Sect. 4.2. Last but not least, a comparative evaluation of SO₂ from both, anthropogenic and outgassing volcanic emissions, is detailed in Sect. 4.3.

Figure 1 Geographical representation of SO₂ VCDs from EMAC, calculated with the standard AK (left panels), against TROPOMI retrievals (right panels) in DU for the first 6 months of 2019. The grey zones represent areas with no valid measurements.

Figure 2 Continuation of Fig. 1 for the last 6 months of 2019.

4.1 SO₂ Vertical Column Densities in 2019: EMAC vs. TROPOMI

A comprehensive analysis of the global distribution and discrepancies of SO₂ VCDs as derived from the RD1SD-base-01 EMAC simulation and from TROPOMI/Sentinel-5P satellite observations for the year 2019 is performed. Figures 1 and 2 show the monthly spatial distribution of SO₂ VCDs, with EMAC results on the left and TROPOMI data on the right.

Overall, both datasets display low SO₂ VCDs (< 0.5 DU) over regions like Western Europe, Africa, and Australia. However, notable discrepancies appear in areas such as Southern Italy, Northeast China, India, Central America, and Papua New Guinea, where EMAC simulates persistently higher values (> 1 DU) compared to TROPOMI (< 0.5 DU) throughout most of the year, except in June and July. In August, elevated TROPOMI signals over Papua New Guinea are caused by the Ulawun eruption (Kloss et al., 2021).

In June and July, TROPOMI detects higher SO₂ over northern latitudes, particularly following the June 2019 Raikoke volcanic eruption in Russia (De Leeuw et al., 2021). These enhanced signals are absent in the EMAC simulation due to the omission of such episodic volcanic emissions in the CMIP6 inventory.

Figure 3 SO₂ VCD differences between the EMAC model results and TROPOMI retrievals in DU as a 10 months (without June and July) mean.

To disregard the missing volcanic SO₂ emissions from Raikoke and Ulawun in the EMAC simulation, a 10 months average excluding June and July 2019 is used for the comparison with TROPOMI observations. The resulting differences in SO₂ VCDs are shown in Fig. 3.

The comparison between EMAC simulations and TROPOMI observations reveals significant discrepancies in Northeastern China, India, Southern Italy, and at localized sources in the Southern Hemisphere, where EMAC consistently overestimates SO₂ VCDs by 1.5 to 2 DU.

Conversely, a better agreement (0.1–0.5 DU differences) is found in Fig. 3 over Western Europe, the USA, Africa, Australia, and most of South America, where SO₂ emissions are relatively low. The most significant differences arise in regions with large anthropogenic emissions or active volcanoes, such as Etna (Italy), Mt. Fuji (Japan), Ulawun (Papua New Guinea), Nevado Ojos del Salado (Chile), and Kunlun (Tibetan Plateau). Elevated SO₂ VCDs in Beijing and India highlight the impact of industrial emissions.

To ensure accurate comparisons, different SO₂ sources require appropriate TROPOMI product types and Averaging Kernels (AKs). For instance, AK_15km is used to assess the impact of explosive volcanic eruptions (e.g., Raikoke, Ulawun) on SO₂ emissions and deposition, as detailed in Sect. 4.2. For anthropogenic and outgassing volcanic emissions, a relative comparison of SO₂ VCD magnitudes is performed, without directly assessing the absolute emitted and deposited SO₂ mass from individual point sources.

The differences between the AKs and the corresponding VCDs are explained in detail in Appendix A2.

4.2 Effects of volcanic eruptions on the simulated atmospheric SO₂ in the EMAC Model

TROPOMI detected significant SO₂ signals in June and July 2019, which were not visible in the EMAC RD1SD-base-01 simulation due to the lack of eruptive events in the applied emission inventories (Fig. 3). The Raikoke (48.29° N, 153.25° E) and Ulawun (5° S, 151° E) eruptions injected volcanic SO₂ into the stratosphere, increasing stratospheric Aerosol Optical Depth (sAOD) across both hemispheres (Kloss et al., 2021). Raikoke, the largest SO₂ injection into the Upper Troposphere and Lower Stratosphere (UTLS) since Nabro (2011), released about 1.5 ± 0.2 Tg SO₂ (Muser et al., 2020; De Leeuw et al., 2021), while TROPOMI estimated 0.14 Tg (SO₂) from Ulawun in June and 0.2 Tg (SO₂) in early August (Kloss et al., 2021).

Table 5 Input parameters of the three sensitivity simulations used in this study. The table lists the prescribed volcano emissions for each sensitivity simulation, with their eruption time and the injected SO₂ mass into the stratosphere and into all model layers.

To investigate the effects of these two eruptions, we performed additional sensitivity simulations (June–December 2019, initialized with results end of May 2019 of the RD1SD-base-01 simulation), in which the volcanic SO₂ emissions were taken into account with the submodel TREXP (Jöckel et al., 2010). Two different SO₂ vertical emission profiles (named StratProfile and VolRes1.5, respectively, adopted from De Leeuw et al., 2021) were utilized for the Raikoke eruption. In EMAC, the StratProfile has been applied in the RD1SD-raik-02 (raik-02) sensitivity simulation, whereas the VolRes1.5 injection profile was used in the RD1SD-raik-03 (raik-03) simulation. In both simulations the SO₂ release is approximately 1.5 Tg SO₂ (1.57 for raik-02 and 1.5 for raik-03) into the atmosphere, the only difference lies in the vertical distribution. In raik-02 (StratProfile profile from De Leeuw et al., 2021), 69 % of the volcanic SO₂ mass (1.09 Tg) is emitted into the stratosphere, with the primary peak occurring at 12–13 km altitude. In contrast, for raik-03 (equivalent to the VolRes1.5 profile in De Leeuw et al., 2021), only 43 % of the SO₂ mass (0.64 Tg) is emitted into the stratosphere, with the primary peak located around 10 km altitude in the upper troposphere. Additionally, raik-04 is based on the setup from raik-02 with additional emissions stemming from the Ulawun volcano. Table 5 lists the set-ups applied in all sensitivity simulations.

For a detailed study focusing solely on the SO₂ mass burden originating from volcanic eruptions, the results of the RD1SD-base-01 reference simulation have been subtracted from the sensitivity simulations. This approach isolates the SO₂ mass specifically attributable to volcanic activity from other anthropogenic or outgassing volcanic emissions, thereby enabling a more precise analysis of its impact. Figure 4 illustrates the SO₂ mass emitted at different altitudes in both simulations. The altitude of the SO₂ emissions significantly influences their atmospheric distribution and dispersion patterns, impacting their climate effects, lifetime, and oxidation rates (Höpfner et al., 2015).

Figure 4 Shown is the estimated total emitted SO₂ mass for the Raikoke eruption on 21 and 22 June 2019 for two different EMAC set-ups. In the first one, represented by raik-02 (black line), most of the SO₂ mass (69 %) is emitted into the stratosphere (De Leeuw et al., 2021). Conversely, in the second set-up (raik-03,red line), most of the SO₂ mass (57 %) is emitted into the troposphere (De Leeuw et al., 2021).

Figure 5 The daily evolution of the global SO₂ mass (Tg of SO₂) after the 2019 Raikoke and Ulawun volcanic eruptions retrieved from TROPOMI and different EMAC sensitivity simulations. raik-02 and raik-03 denote simulations with SO₂ mass from Raikoke emitted mostly in the stratosphere and troposphere, respectively. raik-04 is based on raik-02 with additionally emitted SO₂ mass into the stratosphere originating from the Ulawun volcano. The blue shading represents the uncertainty estimate for the TROPOMI product.

The combined SO₂ mass in the troposphere and stratosphere represents the total SO₂ mass burden from the Raikoke eruption. This mass burden is then compared with the global SO₂ mass burden derived from the TROPOMI/Sentinel-5P satellite after the Raikoke eruption and until mid-July 2019 (De Leeuw et al., 2021) (see Fig. 5). The raik-02 simulation (black line), which assumes a larger emission of SO₂ into the stratosphere, aligns more closely with TROPOMI's derived SO₂ data than the raik-03 simulation (orange line), in which a larger proportion of SO₂ is released into the troposphere.

Both simulations accurately capture the SO₂ mass burden peak at approximately 1.8 Tg. The peak values of the EMAC simulated Raikoke SO₂ mass are slightly larger than the total emissions presented in Table 5, due to the application of AK_15km to the EMAC results. Moreover, raik-02 (black curve) shows a better long-term agreement with TROPOMI estimates than raik-03 (orange curve), however both consistently remaining within the estimated uncertainty range (as derived by Theys et al., 2017). For TROPOMI, uncertainties of SO₂ in the stratosphere are approximately ± 30 % of the retrieved VCDs (Theys et al., 2017). Conversely, the raik-03 simulation exhibits a more rapid decline of the SO₂ mass in the stratosphere, compared to the TROPOMI data after the Raikoke eruption, suggesting a faster removal of SO₂ from the atmosphere, with an exception for the first two days after the SO₂ mass peak, where raik-03 agrees better with the measurement based estimates than raik-02.

The discrepancies between TROPOMI SO₂ based estimates and EMAC simulation results can be attributed to the vertical injection profiles (which might deviate from reality) or to the rate of sulfur removal from the atmosphere. Cai et al. (2022) suggest that additional injections are required after the initial Raikoke plume to accurately replicate the observed SO₂ mass, underscoring the complexity of modeling volcanic SO₂ emissions and their interactions in the atmosphere. To address these differences, a sensitivity simulation labeled raik-04 (green curve in Fig. 5) was conducted. It is based on the setup of the raik-02 simulation, chosen because it best matches the temporal evolution of TROPOMI derived SO₂ mass, but additionally includes emissions from the Ulawun volcano in the Southern Hemisphere. The Ulawun eruptions on 26 June 2019, at 12:00 UTC and 3 August 2019 at 12:00 UTC were taken into account, with each eruption lasting 6 h. During the first eruption, 0.14 Tg of SO₂ was injected at altitudes between 16 and 19 km in the model. For the second eruption, 0.2 Tg of SO₂ was emitted at altitudes between 11 and 15 km. As shown in Fig. 5, these adjustments improved the temporal evolution of SO₂ mass, slowing the decline in the raik-04 simulation following the Ulawun emission injections on 26 June due to the increased SO₂ mass in the stratosphere.

During the Raikoke eruption and up to 29 June 2019, all simulations consistently show a continuous decrease of the SO₂ mass. TROPOMI data indicates a decline rate of approximately 0.08 Tg (SO₂) d⁻¹, which is slower than the decline simulated by the EMAC model. The raik-03 simulation shows a decline rate of around 0.14 Tg (SO₂) d⁻¹, indicating a lower SO₂ mass compared to TROPOMI. Conversely, the raik-02 simulation exhibits a slower decline rate of 0.1 Tg (SO₂) d⁻¹. However, raik-04 aligns most closely with TROPOMI, with a decline rate of 0.09 Tg (SO₂) d⁻¹. This closer match can be attributed to the additional Ulawun emissions injected into the stratosphere on 26 June 2019, which increases the mass of SO₂ in the atmosphere. Note that for the first two days following the SO₂ mass peak, TROPOMI observations indicate a rapid decline rate of approximately 0.15 Tg (SO₂) d⁻¹. This rate aligns most closely with the raik-03 simulation, where a significant amount of SO₂ is emitted at lower altitudes, resulting in a decline rate of 0.14 Tg (SO₂) d⁻¹. In contrast, the raik-02 and raik-04 simulations, which involve the majority of SO₂ being emitted into the stratosphere, exhibit a slower decrease of about 0.12 Tg (SO₂) d⁻¹. Between 29 June and 15 July 2019, all EMAC simulations show a decrease similar to TROPOMI, at a rate of approximately 0.05 Tg (SO₂) d⁻¹.

Over the period from 22 June to 15 July 2019, raik-04 aligns most closely with TROPOMI observations, by simulating about 3 % lower SO₂ mass than TROPOMI over the entire period. raik-02 indicates a mean relative difference of 10 %, while raik-03 simulates lower values than TROPOMI with a mean relative difference of 25 % due to differences in the decline rate and SO₂ vertical injection profile. The fact that all EMAC simulations fall within the 30 % uncertainty range of the SO₂ total column in TROPOMI provides confidence that EMAC correctly captures the main processes required to represent SO₂ oxidation and deposition.

To summarize, this analysis shows the capability of all sensitivity simulations, to reproduce the TROPOMI measured peak after the Raikoke and Ulawun eruptions. Furthermore, the consistent decay rates between TROPOMI data and sensitivity simulations, particularly raik-04 (which encompasses both, the Raikoke and Ulawun eruptions in the stratosphere), as well as accounting for the deposition of most of the initially emitted SO₂ mass within EMAC, further underline the model's realism in capturing the intricate processes of SO₂ emission, oxidation, and deposition associated with volcanic eruptions. Nevertheless, over extended durations, various factors such as simulated wind patterns, radiative heating effects, and mixing dynamics can introduce deviations between model results and real-world observations. These complexities highlight the ongoing challenges in achieving complete concordance between model simulations and empirical data over prolonged temporal scales.

4.3 Evaluation of SO₂ from anthropogenic and outgassing volcano emissions

EMAC SO₂ emissions from both, anthropogenic and outgassing (non-eruptive) volcanic sources, are derived from prescribed emission inventories, specifically CMIP6 and the AeroCom Project, respectively. Each of these emission inventories is based on distinct assumptions that may not accurately reflect the actual emitted SO₂ masses and injection heights. Consequently, only a relative comparison of SO₂ hotspots (i.e. with large SO₂ emissions) and background regions between the EMAC model results and TROPOMI observations is feasible. Therefore, for the comparison discussed in this section, only the ratios and the relationship between EMAC results and retrieved TROPOMI SO₂ VCDs are investigated.

Table 6 Ratios of SO₂ VCDs between EMAC results and TROPOMI retrievals over different volcano types in 2019. The in EMAC emitted SO₂ in $\mathrm{molec}.{\mathrm{m}}^{-\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$ at different altitudes, between 989 and 577 hPa, are listed.

Table 6 lists the ratios between SO₂ VCDs simulated with the EMAC model compared to those retrieved from TROPOMI, around specific outgassing volcanoes, presented as a yearly mean for 2019. The table indicates that EMAC SO₂ VCDs are generally larger than TROPOMI values over volcanic regions. Note that, when comparing the satellite data to model results folded with the AKs, COBRA data are subject to an uncertainty of approximately 27 % to 32 % on the retrieved SO₂ column, mainly due to instrumental noise (Theys et al., 2022).

The comparison reveals notable discrepancies between the SO₂ VCDs over several volcanic regions, with EMAC-to-TROPOMI ratios of approximately 7 for Etna, 6 for Tajumulco, and 5 for Mt. Fuji. In contrast, a lower ratio of about 2.5 is observed over Nevado Ojos del Salado. Notably, these differences are out of the error margin of the satellite measurements, which typically ranges between 27 % and 32 %.

These differences are attributed to the volcanic emission inventory used within the EMAC model. The larger ratios derived for the first three volcanoes could be attributed to both, the SO₂ emission masses and the emission heights. In the AeroCom inventory (Dentener et al., 2006), the Etna, Tajumulco, and Mt. Fuji volcanoes are considered not only as outgassing, but also explosive volcanoes. The top of Etna is approximately 3300 m, Tajumulco about 4000 m, and Mt. Fuji about 3700 m, with emissions reaching up to around 4500 m (577 hPa). This indicates explosive volcanic activity, with emissions ranging from the “top of the volcano + 500 m” to “top of the volcano + 1500 m”, as detailed in Sect. 2.2. Conversely, the volcano in Chile is categorized solely as an outgassing volcano, thus showing the lower ratio between EMAC results and TROPOMI retrievals. Continuously outgassing volcanoes, within the used model setup, emit from “the height of the volcano × 0.67” up to the height of the volcano.

Figure 6 Geographical map showing the regions selected for this study.

For anthropogenically influenced regions, the EMAC results also shows larger SO₂ VCDs compared to those from TROPOMI, but with smaller ratios in most regions compared to volcanic areas. Figure 6 shows a geographical map depicting the selected areas which have been investigated quantitatively, and Table 7 lists the ratio of EMAC SO₂ VCDs to TROPOMI VCDs in these SO₂ background and hotspot regions. In background regions, a specific area in central Africa (coordinates: 12° N, 15° E to 2° N, 27° E, respectively) and the South Atlantic Ocean (coordinates: 20° S, 20° W to 30° S, 5° W) reveal small discrepancies between the SO₂ VCDs (10 %–20 %). Larger differences are derived in anthropogenically influenced regions: Europe (1.6 factor, 60 % difference), the USA (1.8 factor, 80 % difference), and India (2.5 factor, 150 % difference). In Northeastern China the factor is 3.2 (220 % difference), while in Southeastern China it is 2 (100 % difference).

Table 7 Ratios between SO₂ VCDs from EMAC results and TROPOMI retrievals in different background and SO₂ emission hotspot regions in 2019.

The regional discrepancies of SO₂ VCDs between EMAC and TROPOMI can be attributed to the large SO₂ emissions in these regions, originating from the CMIP6 emission inventory, as well as to their different original spatial resolutions. EMAC, with a coarse resolution of 300 km × 300 km, averages emissions over larger areas, potentially underestimating localized SO₂ peaks from sources like volcanoes or power plants. In contrast, TROPOMI/Sentinel-5P, with a much finer resolution of 3.5 km × 7 km, captures fine-scale variations. This difference might lead to discrepancies, especially in regions with strong emissions, since the emissions in the model become instantaneously diluted by spreading the emitted mass over the model grid-boxes. For instance, in areas like India and China, where large SO₂ emissions take place, the ratios between EMAC and TROPOMI are larger, while background regions such as Africa show minimal differences in 2019, falling within TROPOMI's uncertainty range. In regions like Europe, the USA, and Southeastern China, which have lower SO₂ emissions than China and India for example, there is a better agreement between the two datasets.

Since TROPOMI only provides total VCD values, a detailed analysis of the vertical profile between both datasets is not possible. Therefore, it is difficult to ascertain whether the differences originate near the surface or higher up in the atmosphere. To address this, a comparison of the simulated SO₂ concentrations at the Earth's surface is conducted in the next section.

5 Evaluation of simulated SO₂ with ground-based measurements

To complement the inter-comparison of EMAC results with VCDs derived from TROPOMI data, we next compare the EMAC results with ground-based measurements from observation networks in major SO₂-emitting regions worldwide. For this comparison, it is necessary that both datasets are aligned on the same latitude-longitude grid, as explained in Appendix A4. Specifically, the analysis centers on SO₂ concentrations and sulfur deposition fluxes over the USA, Europe, and at selected observational stations in China and Japan. These regions are chosen due to the availability of extensive and reliable datasets covering a two-decade period, from 2000 to 2019.

It is important to note that for the time series analysis in this section, mean/average values and standard deviation across stations are calculated for each year for both, the EMAC model results and the data from observational networks. The corresponding calculations and formulas can be found in the Appendix B.

The results are presented in a aggregated form, i.e. we map the station data onto the model grid.

5.1 Sulfur concentration and deposition in the USA

For the United States, sulfur species simulated with EMAC near the Earth's surface (i.e. the lowermost grid box) are compared with observation data obtained from the CASTnet network. As detailed in Appendix A3.1, CASTnet provides surface-level observations, including monthly and yearly mean SO₂ concentrations in µg m⁻³, and sulfur wet deposition fluxes in $\mathrm{kg}\left(\mathrm{S}\right){\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$. The SO₂ concentrations and sulfate amounts in precipitation samples are measured, whereas the dry deposition fluxes are simulated based on a multi-layer model.

Figure 7 Geographical distribution of mean SO₂ concentrations for the years between 2000 and 2019 at the Earth's surface as simulated with EMAC (left) and observed at the CASTnet sites (aggregated onto the EMAC grid, right) in µg m⁻³. The red and blue boxes indicate the regions, where the SO₂ emissions from the CMIP6, EDGAR5 and EDGAR8 emission inventories are compared (details see text).

Figure 7 shows the SO₂ concentration measured at the CASTnet sites (right panel, as mentioned above aggregated onto the model grid) and the EMAC simulated concentration (left panel). In both cases, 20 year averages are calculated. It is important to note that the EMAC results are only shown for grid boxes, where observational stations are located, which explains the presence of “empty boxes” in the EMAC model results. The figure indicates that Eastern USA sites exhibit larger SO₂ concentrations compared to the sites in the Western region in both datasets. This disparity is attributed to the higher density of SO₂ emission sources in the Eastern USA compared to the Western part. This is also reported by Hardacre et al. (2021) and Qu et al. (2019).

As shown in Fig. 7, CASTnet (right panel) measures approximately 7 µg m⁻³ at some individual sites in the Eastern USA, with other Eastern sites showing very low SO₂ concentrations of about 0.5 µg m⁻³. On the left panel in Fig. 7, the RD1SD-base-01 simulation driven by the CMIP6 inventory indicates that SO₂ concentrations at some individual Eastern sites are lower than those reported by CASTnet. However, on average, EMAC results show overall consistent SO₂ concentrations between 1.5 µg m⁻³ and a maximum of about 5.5 µg m⁻³ at Eastern sites. In the Western region, both datasets show lower SO₂ concentrations, averaging around 1 µg m⁻³ across all sites. However, the RD1SD-base-01 simulation results in larger SO₂ concentrations reaching up to 3 µg m⁻³ at some Western sites.

Figure 8 Time series of SO₂ concentrations from both, EMAC and CASTnet between 2000 and 2019 in the Western (a), and Eastern USA (b). The middle panels (c and d) show the evolution of wet and dry sulfur deposition fluxes between 2000 and 2019. For the calculation of deposition fluxes, the deposited sulfate and SO₂ were converted to sulfur equivalent. The right panels (e and f) show the comparison between the temporal evolution of CMIP6, EDGAR5 and EDGAR8 anthropogenic sulfur emissions in Western and Eastern USA, as a yearly area integral between 2000 and 2019 (2015 for EDGAR5) of all the emission inventory grid boxes situated in the region marked by the blue and red boxes, respectively, in Fig. 7.

Figure 8 shows the comparison between both, the RD1SD-base-01 simulation results and CASTnet SO₂ concentrations as well as sulfur deposition fluxes across different regions in the USA. Specifically, SO₂ concentrations and sulfur deposition fluxes from Western sites (Fig. 8a and c, respectively), and from Eastern sites (Fig. 8b and d, respectively), are shown. The comparison involves calculating the annual mean of SO₂ concentration in µg m⁻³ and of sulfur deposition fluxes (wet and dry deposition) in $\mathrm{kg}\left(\mathrm{S}\right){\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$, averaged over Eastern, and Western USA sites, respectively. For both, surface SO₂ concentration and sulfur deposition flux, the RD1SD-base-01 simulation driven by the CMIP6 emission inventory effectively captures the decline in both regions of the USA for the period 2000–2019. As shown in Fig. 8, the model tends to simulate larger surface SO₂ concentrations than CASTnet by a factor of 2 in the Western region over the 20 year period (Fig. 8a), while showing approximately a factor of 1.2 larger SO₂ concentrations over the Eastern USA (Fig. 8b). For the western USA, EMAC shows decreasing surface SO₂ concentrations after 2000, which brings the simulated results into better agreement with the observations over time (see Fig. 8a). The large standard deviation derived from the datasets are attributed to the extensive dispersion of sulfur sources across a broad geographical area.

In the Eastern and Western USA, the largest part of sulfur removal occurs via wet deposition. This is effectively simulated by EMAC in agreement with CASTnet observations. For the calculation of deposition flux within the RD1SD-base-01 simulation, the deposited sulfate and SO₂ were converted to a sulfur equivalent. In Fig. 8d, EMAC shows a lower sulfur deposition flux over the Eastern USA, for wet deposition (10 % lower EMAC values compared to CASTnet) and larger values for dry deposition (30 % larger EMAC values compared to CASTnet) over the 20 year period. In the Western USA, EMAC also simulates a 5 % lower sulfur wet deposition flux compared to CASTnet over the entire 20 year period (Fig. 8c). Here, EMAC does not show lower values over all the time range, but indicates larger wet sulfur deposition between 2002 and 2008. For the sulfur dry deposition flux (orange lines), EMAC shows a factor of 2 larger values over the entire 20 year-period compared to CASTnet.

Conrad-Rooney et al. (2023) analysed the wet-deposition of ${\mathrm{SO}}_{\mathrm{4}}^{-\mathrm{2}}$ for urban and suburban site classes separately, whereas we did not distinguish different sites and analysed the total sulfur deposition. Thus, the results are not directly comparable. Nevertheless, Conrad-Rooney et al. (2023, see their Figs. 2D and 5D) estimated for the year 2000 a wet deposition rate (urban stations) of approx. 20 $\mathrm{kg}\left({\mathrm{SO}}_{\mathrm{4}}^{-\mathrm{2}}\right){\mathrm{ha}}^{-\mathrm{1}}$, equivalent to 6.7 $\mathrm{kg}\left(\mathrm{S}\right){\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$. For 2018, they derived a wet deposition rate of 5 (rural) to 7.5 (urban) $\mathrm{kg}\left({\mathrm{SO}}_{\mathrm{4}}^{-\mathrm{2}}\right){\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$, equivalent to 1.7 to 2.5 $\mathrm{kg}\left(\mathrm{S}\right){\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$, respectively. These numbers agree, given the limitations of this comparison, with our results for the eastern USA (see Fig. 8d, blue area). Benish et al. (2022) derived for the USA a decrease of the total sulfur depostion from 5.3 kg (S) ha⁻¹ in 2002 to 1.8 kg (S) ha⁻¹ in 2017 (see their Fig. S14b). These total fluxes are lower than the wet deposition fluxes derived by Conrad-Rooney et al. (2023) and than our results, but the derived trends are consistent.

Since the concentration simulated by the model is directly affected by the prescribed emissions, it is important to assess potential uncertainties of the applied emission inventory (CMIP6). For this, we present the comparison with the EDGAR5 (see Sect. 2.2), and the EDGAR8 emission inventories. The temporal evolution of sulfur emissions in Tg (S) yr⁻¹ from the CMIP6 emission inventory (orange lines in Fig. 8e and f) is compared with those from the EDGAR5 (black lines in Fig. 8e and f) and the EDGAR8 emission inventory (green lines in Fig. 8e and f). This comparison highlights the discrepancies and potential biases between the different emission inventories, which would be reflected in model results, if based on the alternative inventory. For the calculation of sulfur emissions, the emissions from SO₂ were also converted to sulfur equivalent. The anthropogenic emissions (particularly from fossil fuels, ship, road, and aircraft sectors) from both emission inventories are calculated as a yearly area integral over the Western USA region (see Fig. 8e), and over the Eastern USA region (see Fig. 8f). Both regions are shown in Fig. 7, where the red and blue boxes represent the selected Eastern and Western USA regions, respectively.

For the western USA (Fig. 8e), the CMIP6 inventory has 50 % larger sulfur emissions than the EDGAR5 inventory between 2000 and 2015, and 60 % larger sulfur emissions than the EDGAR8 inventory between 2000 and 2019. The picture is different for the Eastern USA (Fig. 8f), where the CMIP6 inventory indicates 10 % less sulfur emissions than the EDGAR5 inventory between 2000 and 2015, and more or less the same sulfur emissions compared to the EDGAR8 inventory over the entire time range. The differences between the emission inventories, particularly the larger anthropogenic sulfur emissions in CMIP6 compared to those in EDGAR5 and EDGAR8 over the Western USA, are a major factor contributing to the larger SO₂ concentrations simulated by EMAC compared to those observed at the CASTnet in Western USA (see Fig. 8a). Consequently, using the EDGAR5 or the EDGAR8 emission inventory over the Western USA would likely result in smaller SO₂ concentrations in the EMAC model, thereby reducing the bias between the CASTnet measurements and the EMAC results in that region.

For the final year of the study, 2019, the used CMIP6 emission inventory used for the RD1SD-base-01 simulation leads overall to a larger SO₂ concentration by a factor of approximately 1.6 compared to the CASTnet measurements.

Figure 9 Geographical distribution of mean SO₂ concentrations for the years between 2000 and 2019 at Earth's surface as simulated with EMAC and observed at the EMEP sites in µg m⁻³. The blue box indicates the region, where the SO₂ emissions from the CMIP6, EDGAR5 and EDGAR8 emission inventories are compared (details see text).

5.2 Sulfur concentration and deposition in Europe

In Europe, 48 observational stations from the EMEP database are analyzed, as detailed in Appendix A3.2. First, the spatial distribution of SO₂ concentration over Europe from both datasets, is shown in Fig. 9. EMAC results are only shown for grid boxes where observational stations are located, which explains the presence of “empty boxes” in the EMAC model results. Here, EMAC (left panel) shows the largest SO₂ concentrations (between 4 and 8 µg m⁻³) over central East Europe, with lower SO₂ concentrations (between 0.3 and 3 µg m⁻³) over sites in the United Kingdom and Western Europe. On the other side, EMEP measures one very large SO₂ concentration (about 7.5 µg m⁻³) over a grid box situated in Serbia, while showing lower SO₂ concentrations at the remaining sites (between 0.2 and 3.3 µg m⁻³).

Figure 10 Temporal evolution of mean annual surface SO₂ concentration and sulfur wet deposition flux from EMAC (blue lines) and EMEP (red lines) between 2000 and 2019 at sites in Europe (top panels). The lower panel shows the comparison between the temporal evolution of CMIP6, EDGAR5 and EDGAR8 anthropogenic sulfur emissions as a yearly area integral between 2000 and 2019 (2015 for EDGAR5) of all the inventory grid boxes situated in the region marked by the blue box in Fig. 9.

Furthermore, Fig. 10 illustrates the decline of SO₂ concentration across Europe, which is well captured by EMAC throughout the 20 year period (top left panel in Fig. 10). Consequently, the temporal reduction in sulfur loss due to wet deposition is also accurately represented by EMAC (top right panel in Fig. 10). However, the used CMIP6 emission inventory tends to produce larger SO₂ concentrations within the RD1SD-base-01 simulation compared to observational data within the European domain (top left panel in Fig. 10). Specifically, SO₂ surface concentrations from the model show a gradual decline between 2000 and 2012, with an annual decrease rate of 0.11 µg m⁻³, whereas ground-based observational data indicate a slower reduction rate of 0.04 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$ during the same period. After 2012, both, model results and observational datasets, exhibit a more pronounced acceleration in the decline of SO₂ concentration, with rates of approximately 0.22 and 0.12 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$, respectively.

Over the entire 20 year period, EMAC driven by the CMIP6 emission inventory consistently shows larger annual mean surface SO₂ concentrations in Europe by a factor of approximately 1.8 compared to the EMEP dataset. Specifically, EMAC indicates larger SO₂ concentrations relative to observational data by a factor of 2 between 2000 and 2012, and with a lower factor of approximately 1.5 between 2012 and 2019. Regarding sulfur wet deposition flux, EMAC also simulates consistently larger values than observed by EMEP by a factor of 1.3 over the entire 20 year period. Notably, a consistency is observed in the decline rates of both datasets, characterized by a yearly mean decrease of about 0.05 µg m⁻³ throughout the duration from 2000 to 2019. This trend, calculated by linear regression, is comparable to the trend of −0.067 µg m⁻³ derived by Aas et al. (2024, see their Table S3 for the perios 2000–2019), which has been derived with a different method and based on a slightly different number of stations. Despite the differences in SO₂ concentration and sulfur deposition flux, it is noteworthy that the model exhibits a good alignment with observational data, as shown in the temporal progression of both, SO₂ concentration and sulfur deposition flux.

Similar to the study of sulfur emissions over the USA (see Sect. 5.1), the prescribed CMIP6 emission inventory in Europe used for the RD1SD-base-01 simulation, shows differences in the temporal evolution of emitted sulfur compared to the EDGAR5 and EDGAR8 emission inventories (see the low panel of Fig. 10). Here, the anthropogenic sulfur emissions (particularly from fossil fuels, ship, road, and aircraft sectors) from all emission inventories are calculated as a yearly area integral between 2000 and 2019 (2015 for EDGAR5) over a region in Europe. The chosen region is marked by a blue box, as shown in Fig. 9. In the lower panel of Fig. 10, an identifiable reduction of anthropogenic sulfur emissions across Europe is evident throughout the temporal evolution of both emission inventories. Specifically, the sulfur emissions in the CMIP6 inventory are about 20 % larger than those from the EDGAR5 and EDGAR8 inventories over the available time frames over Europe. This implies that using the EDGAR5 or EDGAR8 emission inventory would result in 20 % less SO₂ being emitted in the model. Consequently, using the EDGAR5 or EDGAR8 emission inventory in Europe would likely result in lower SO₂ concentrations and therefore lower sulfur wet deposition fluxes in the EMAC model, thereby reducing the bias between the EMEP measurements and the EMAC results.

For the final year, 2019, the used CMIP6 emission inventory used for the RD1SD-base-01 simulation leads to a larger EMAC SO₂ concentration by a factor of approximately 1.45 compared to EMEP measurements.

Figure 11 Geographical distribution of mean SO₂ concentrations for the years between 2000 and 2019 at Earth's surface as simulated with EMAC and observed at the EMEP sites in µg m⁻³. The blue box indicates the region, where the SO₂ emissions from the CMIP6, EDGAR5 and EDGAR8 emission inventories are compared (details see text).

5.3 Sulfur concentration and deposition in China and Japan

For East Asia, a comprehensive investigation of SO₂ and its associated processes, particularly within China, is imperative due to the substantial contribution of Chinese SO₂ emissions, which account for 64 %–71 % of the total emissions across Asia (Kuribayashi et al., 2012). However, as described in Appendix A3.3, it is complicated to find representative monitoring stations providing continuous, long-term datasets of measured SO₂ concentrations and sulfur deposition fluxes. Consequently, only 5 stations in Southeastern China and 9 over Japan from the EANET network have been selected for the present study. This selection was based, as for the previous networks in the USA and Europe, on the availability of long-term measurements (2000 to 2019 for this study) providing both, measured SO₂ concentration in µg m⁻³ and wet deposition sulfur fluxes in kg (S) ha⁻¹ in EANET. The geographical distribution and SO₂ concentrations of selected stations in both, EMAC and EANET datasets at sites situated in China and Japan are shown Fig. 11.

Here, both, EMAC and EANET datasets, indicate larger SO₂ concentrations at sites situated in China compared to those in Japan. As shown on the right panel of Fig. 11, EANET measures large SO₂ concentrations with a maximum of about 12 µg m⁻³ at three sites in China, while showing lower SO₂ concentrations ranging between 1 and 3 µg m⁻³ at the sites in Japan. On the left panel of Fig. 11, the EMAC results show overall more consistent SO₂ concentrations with one site in China showing around 12 µg m⁻³, while SO₂ concentrations at the other sites in China and Japan range between 2 and 8 µg m⁻³.

Figure 12 Temporal evolution of mean annual surface SO₂ concentration and sulfate wet deposition flux from EMAC (blue lines) and EANET (red lines) between 2000 and 2019 at the sites in China and Japan (top panels). The lower panel shows the comparison between the temporal evolution of CMIP6, EDGAR5 and EDGAR8 anthropogenic sulfur emissions as a yearly area integral between 2000 and 2019 (2015 for EDGAR5) of all inventory grid boxes situated in the region marked by the blue box in Fig. 11.

Lu et al. (2010) and Ohara et al. (2007) reported a significant increase of SO₂ emissions in China during the early 2000s, a trend confirmed by CMIP6, EDGAR5 and EDGAR8 emission inventories, as depicted in the lower panel of Fig. 12 (SO₂ emissions were converted to sulfur equivalent). This rise of sulfur emissions is also reflected in the SO₂ concentrations and the wet deposition sulfur fluxes, as illustrated in the same figure (top left and top right panels, respectively). EMAC indicates an overall increase of SO₂ concentration and the deposited sulfur from 2000 till 2014. However, data from EANET reveals a different trend. According to the EANET network, SO₂ concentration increased until 2006 at a rate of 0.7 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$. On the other side, the deposited sulfur also increased until 2006 at a rate of 0.035 $\mathrm{kg}\left(\mathrm{S}\right){\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$. Afterwards, a decline of SO₂ concentration by 0.5 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$ and of the wet deposition sulfur flux by 0.4 $\mathrm{kg}\left(\mathrm{S}\right){\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ until 2014, was measured.

Following 2014, numerous reports from environmental agencies and satellite observations have consistently indicated a significant decrease of China's SO₂ emissions. Studies by Wei et al. (2023) and van der A et al. (2016) have highlighted that SO₂ emissions in China experienced a sharp decline post-2014, nearly vanishing by 2020. This reduction of SO₂ emissions after 2014 is also noticeable in the time evolution of sulfur emissions from the CMIP6 inventory utilized in this study (lower panel in Fig. 12). From 2014 to 2019, a reduction exceeding 50 % has been documented across the China-Japan region, with sulfur emissions decreasing from about 22 to 12 Tg (S) yr⁻¹ in the CMIP6 emission inventory. This decline of sulfur emissions notably influences the SO₂ concentration and the wet deposition sulfur flux over that region. Consequently, a remarkable reduction, from about 5.5 µg m⁻³ of SO₂ concentrations in 2014 to approximately 3 µg m⁻³ by 2019, is detected (see the upper left panel in Fig. 12). Additionally, a good agreement is also evident between EMAC results and EANET measurements regarding the evolution of the deposited sulfur. Here, the wet deposition sulfur flux from both datasets indicate a decrease from 9 kg (S) ha⁻¹ in 2014 to 5.8 kg (S) ha⁻¹ in 2019 (see the upper right panel in Fig. 12). This is an average decline of 0.53 $\mathrm{kg}\left(\mathrm{S}\right){\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ between 2014 and 2019 based on station data in South China and Japan. Xi et al. (2025) reported a China nationwide average decline rate (2013–2023) of −0.244 $\mathrm{kg}\left(\mathrm{S}\right){\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ and furhter showed that the trends in South and Central China are larger (negative) than this average (see their Fig. 2). Thus, our results can be considerd to be consistent. In addition, Yamaga et al. (2021) found “no clear increase or decrease trends in the S deposition amounts throughout the 15 year period” (2003–2017) at 8 stations in Japan. This is also in line with our results (Fig. 12, top right panel) due to the large standard deviations and because the overall decline is largely driven by emission reductions over China.

In contrast to Europe and the USA, a comparative analysis of SO₂ concentration and of the wet deposition sulfur flux between the RD1SD-base-01- simulation results and measurements in China and Japan reveals noticeable differences in their temporal evolution. The RD1SD-base-01 simulation results based on the CMIP6 inventory consistently exhibits lower SO₂ concentrations than EANET between 2000 and 2012, with a factor of approximately 1.5. Subsequently, larger SO₂ concentrations from EMAC of approximately 15 % are detected in the years between 2012 and 2019. This pattern extends to the wet sulfur deposition temporal evolution, where EMAC consistently simulates lower values than EANET from 2000 until 2014, with a factor of approximately 1.6. However, from 2014 onwards, EMAC exhibits a reversal in trend, indicating larger values than EANET measurements by approximately 8 % until 2019.

In the lower panel of Fig. 12, the anthropogenic emissions from, the used CMIP6 emission inventory, the EDGAR5 and EDGAR8 emission inventories are calculated as a yearly area integral over a region covering China and Japan. The chosen region is marked by the blue box shown in Fig. 11. Again, the emissions from SO₂ were converted to sulfur equivalent.

In 2019, the used CMIP6 emission inventory used for the RD1SD-base-01 simulation leads to a SO₂ concentration of 3.6 µg m⁻³ compared to 2 µg m⁻³ measured by the EANET network, giving a bias of 1.6 µg m⁻³ and a ratio of 1.8.

6 Conclusions and Outlook

The results of this study indicate that while the EMAC simulations demonstrate notable strengths in simulating the sulfur cycle, there are also areas that require improvement. Understanding the tropospheric sulfur budget forms the groundwork for the subsequent analyses and examinations. This analysis is conducted using the EMAC RD1SD-base-01 model simulation, for which the CMIP6 SO₂ emission inventory was applied. The model demonstrates a closed sulfur budget, which has been compared with other results from literature. This assessment indicates a consistent representation of the model's sulfur chemistry, such as emissions, transport, chemical kinetics, and deposition. A closed sulfur budget allows for the evaluation of the EMAC model results against other observational data. By utilizing TROPOMI/Sentinel-5P measurements, the global distribution of natural and anthropogenic SO₂ is identified in the EMAC model, showing different magnitudes compared to the SO₂ VCDs measured by the satellites instrument. Notably, the model indicates larger SO₂ VCDs, especially around regions with outgassing volcanic emissions. This discrepancy is attributed to the AeroCom emission inventory (Dentener et al., 2006) used within the EMAC model, which may not accurately reflect current outgassing volcanic activity. Given that the current emission inventory dates back to 2006, a new emission inventory for outgassing volcanic activities in the troposphere (Brühl et al., 2021) should be applied in future EMAC model setups. On the other hand, the enhanced SO₂ signals associated with the eruptions of the Raikoke and Ulawun volcanoes, as observed by the TROPOMI instrument, are successfully reproduced by the EMAC model when additional volcanic SO₂ emissions are explicitly included in the simulation setup. The temporal evolution of the additional global SO₂ mass is well reproduced by EMAC. Thus, future hindcast simulations with EMAC should also include a representation of eruptive volcanic emissions as proposed by Kohl et al. (2024).

Regarding the regions studied here, the biases calculated between the SO₂ VCDs from EMAC results and TROPOMI measurements and those between the EMAC simulated SO₂ concentrations and ground-based measurements over the USA, Europe and China and Japan in 2019, are consistent with each other. This strengthens the earlier hypothesis that the prescribed SO₂ emissions from the CMIP6 inventory used for the RD1SD-base-01 simulation might be overestimated over these regions, as corroborated by a comparison with the EDGAR5 and the EDGAR8 emission inventories. This underscores the importance of further investigating the EMAC model results using various emission inventories to assess the range and sources of uncertainties, leading to a better understanding of the behavior of SO₂ emissions across different regions.

Following our results, additional studies to further reduce the uncertainties of the knowledge about the atmospheric sulfur budget are required. The selection of the emission inventory significantly influences the simulated SO₂ concentrations, which consequently impacts the deposition processes. To enhance the understanding of these impacts, it is recommended to conduct sensitivity simulations using various emission inventories and evaluate the model results against space-, air-, and ground-based measurements. These simulations will help to quantify the uncertainties and variations associated with different inventories in different regions, leading to more accurate simulated SO₂ concentrations and sulfur depositions.

Future research should focus on evaluating the simulated SO₂ VCDs using high-resolution satellite instruments. Unlike the TROPOMI instrument, which provides a daily global coverage, the Geostationary Environment Monitoring Spectrometer (GEMS) launched in 2020 (Kim et al., 2020), the Tropospheric Emissions: Monitoring of Pollution (TEMPO) instrument launched in 2023 (Zoogman et al., 2017), and the Sentinel-4 instrument (Stark et al., 2013) are dedicated to measuring air quality across Asia, North America, and Europe, respectively, every hour. The available high-frequency, near-real-time data provide an excellent basis for validating the model's ability in simulating atmospheric SO₂ and capturing short-term variations and transient events, such as pollution spikes and weather-related changes.

Appendix A: Description of the used observational data

A1 Satellite observations

In this study retrievals from the TROPOMI instrument on board the Sentinel-5P satellite are employed to investigate the SO₂ VCD. Sentinel-5P is the first Copernicus mission specifically designed for atmospheric monitoring, as mentioned by ESA (2017). Here, two dimensional level-2 products from the TROPOMI instrument are used. These products represent the original SO₂ data retrieved from the spectra observed by TROPOMI, including the geographical coordinates and resolution parameters such as scanline and ground pixel. These dimensions constitute the so-called “satellite orbit”.

The retrieval data of TROPOMI is organized vertically into pressure layers from an a-priori profile of a CTM, namely the Tracer Model 5 (TM5 Huijnen et al., 2010). In the case of SO₂, the data is divided into 34 distinct layers, varying approximately from the Earth's surface to 0.1 hPa (i.e. around 60 km).

In this study two distinct products are used:

1. The operational algorithm retrieves first the concentration of SO₂ integrated along the mean light, i.e., the so-called Slant Column Density (SCD).

   VCD cannot be directly measured from the satellite, thus, the conversion of the SCD into VCD becomes essential. This conversion process relies heavily on the air–mass factor “M”:

   $$\begin{array}{}\text{(A1)}& \text{VCD}={\displaystyle \frac{\text{SCD}}{M}},\end{array}$$

   where the air–mass factor is calculated based on the formulation by Palmer et al. (2001), as follows:

   $$\begin{array}{}\text{(A2)}& M=\int m\left(p\right)\cdot s\left(p\right)\mathrm{d}p.\end{array}$$

   Here, m(p) is a weighting function reflecting the sensitivity of the satellite instrument to different altitudes. This function can be determined through pre-calculation or computational methods using a radiative transfer model. For the actual TROPOMI products, m(p) is given by the Linearized Discrete Ordinate Radiative Transfer (LIDORT) model, as introduced by Spurr et al. (2001). The term s(p) represents the vertical shape factor, which describes the normalized vertical profile of the SO₂ mixing ratio as a function of atmospheric pressure (Palmer et al., 2001). This profile can be obtained a-priori from any CCM or CTM. For instance, the CTM Tracer Model 5 (TM5) model is used as an a-priori profile for the TROPOMI/Sentinel-5P retrieval. Therefore, Eq. (A2) could also be written as:

   $$\begin{array}{}\text{(A3)}& M_{\text{TM5}}=\int m\left(p\right)\cdot s_{\text{TM5}}\left(p\right)\mathrm{d}p.\end{array}$$

2. The Covariance-Based Retrieval Algorithm (COBRA) product represents the latest advancement in SO₂ retrieval techniques from TROPOMI onboard the Sentinel-5 Precursor satellite (Theys et al., 2021).

Finally, it is important to note that the TROPOMI level-2 products are provided with the corresponding averaging kernels (AKs) for each case (Theys et al., 2017). These qualify the vertical sensitivity of satellite instruments and are important for ensuring a fair comparison with other types of data, especially atmospheric chemistry model simulation results (Veefkind et al., 2012). A detailed explanation of the AKs is provided in Appendix A2, together with an explanation of how the model data have been prepared for comparison with satellite measurements.

A2 Post-processing of model data for comparison with satellite observations

For a global investigation of atmospheric sulfur chemistry within the EMAC model, a comparison of model results from the SORBIT submodel (Jöckel et al., 2010), with SO₂ products from satellite measurements (TROPOMI on board of Sentinel-5P in this study) is perfomed. TROPOMI SO₂ products are structured based on scanline and ground pixel, with the scanline representing the direction of the satellite's flight and the ground pixel indicating the resolution of the data. The EMAC model operates on a Gaussian lat-lon grid. To facilitate a meaningful comparison, the TROPOMI data must be regridded to match the grid of the EMAC model. This process involves reducing the fine resolution of TROPOMI to align with the coarser resolution of the model. However, before conducting the comparison, the model data needs to be folded with averaging kernels (AKs) from the retrievals to ensure its compatibility with TROPOMI data (Theys et al., 2022). The Averaging Kernel defines the sensitivity of the retrieved column, obtained from satellite-based measurements, to variations in the true profile of the measured trace gases based on a CCM or CTM (Rodgers, 2000). To properly weight the model data, it first needs to be brought onto the same horizontal and vertical grids as the AKs. This involves horizontally mapping the SO₂ mixing ratio from the SORBIT submodel onto the instruments grid resolution using the nearest neighbor method. Subsequently, a vertical linear interpolation is executed to align the 90 pressure levels of the simulated SO₂ mixing ratio with the 34 layers of the a-priori profiles used for the retrievals. Afterwards, the vertically interpolated SO₂ mixing ratio profiles are converted into a partial column for each of those grid-boxes (i.e. DU or $\mathrm{molec}.{\mathrm{cm}}^{-\mathrm{2}}$). The vertically interpolated model data is then ready to be multiplied at each level with the corresponding averaging kernel and vertically integrated to yield the VCD. This step is important, as it translates the model SO₂ VCD into the signal that would be detected by the satellite. Finally, the VCDs of SO₂ from both, the model and the TROPOMI retrieval, are conservatively regridded from the instrument grid to the original EMAC latitude-longitude grid, and can be compared to each other.

Since the retrieved VCDs depend on simulated a-priori vertical profiles (represented as s(p) in Eq. A2), which in turn depend on prescribed, mainly anthropogenic and volcanic SO₂ emissions, the COBRA dataset (see Sect. A1) provides four different VCDs for specific cases:

• The standard case (or “polluted case”) is obtained using profiles of daily forecasts from the global CTM TM5 (Tracer Model 5, version TM5-chem-v3.0; Huijnen et al., 2010).

• The 1 km case is obtained using 1 km thick box profile concentrating between the surface and 1 km (0 to 1 km), and representing a situation of passive degassing volcanoes or anthropogenic near-surface emissions.

• The 7 km case is obtained using 1 km thick box profile centered at 7 km (6.5 to 7.5 km), indicating a case of a moderate volcanic eruption.

• The 15 km case is obtained using 1 km thick box profile centered at 15 km (14.5 to 15.5 km), reflecting an explosive volcanic eruption case.

It is important to note that in order to compare the model VCDs with the four described VCDs cases from TROPOMI, similar assumptions need to be adopted to ensure a valid comparison. By using the averaging kernel, one can ensure that the comparison between the satellite observations and the model output is meaningful, reflecting the same observational biases and sensitivities. In TROPOMI products, to conserve space, only the total column averaging kernel for the TM5 standard “pollution” case is provided as described by Theys et al. (2017):

$$\begin{array}{}\text{(A4)}& \text{AK}\left(p\right)={\displaystyle \frac{m\left(p\right)}{M_{\text{TM5}}}},\end{array}$$

where M_TM5 represents the total air–mass factor of the vertical profile of the TM5 model and is calculated following Eq. (A3). Importantly, m(p) is consistent across all four cases, and the AK is calculated for the four distinct s(p) profiles. Consequently, we can easily recalculate the AK for each situation by scaling the polluted (or standard) averaging kernel by air–mass factor ratios $M_{\text{TM5}}/M_{\text{box}}$, as described by Eskes and Boersma (2003):

$$\begin{array}{}\text{(A5)}& {\text{AK}}_{\text{box}}\left(p\right)=\text{AK}\left(p\right)\cdot {\displaystyle \frac{M_{\text{TM5}}}{M_{\text{box}}}}.\end{array}$$

Here, p represents the pressure level at which the averaging kernel is stored for the TROPOMI product. $M_{\text{TM5}}/M_{\text{box}}$ serves as the scaling factor reported in TROPOMI products as “sulfurdioxide_averaging_kernel_scaling_box_{1,7,15}km”.

In this study the averaging kernels are referred to as AK_polluted, AK_1km, AK_7km and AK_15km, and the resulting VCDs are, respectively expressed as, VCD_AK_polluted, VCD_AK_1km, VCD_AK_7km and VCD_AK_15km. Following, a detailed description of the different SO₂ VCD products used by TROPOMI, is given:

• VCD_AK_polluted is used for the standard “pollution” case. Here the COBRA product (see Sect. A1) is used due to its enhanced sensitivity to detect low SO₂ column densities. Additionally, a quality flag mechanism is employed to filter out potentially erroneous inputs such as cloudy pixels or missing values, which could deteriorate the results. Therefore, only data points with a quality assurance value above 0.5 (qa_value > 0.5) are considered reliable for analysis, as recommended by Theys (2023).

• VCD_AK_1km and VCD_AK_7km are also retrieved from the COBRA product. However, for volcanic activities, the prerequisite of the quality assurance value is no longer applicable. Instead, the only filtering criteria required is the SZA, where just data with SZA < 70° is considered (Theys, 2023). This can lead to higher signal-to-noise ratios in the satellite measurements, improving the quality of the data collected.

• The VCD_AK_15km product is applied specifically for eruptive volcano emissions. It is recommended to use the operational TROPOMI product (see Sect. A1) for such significant eruptions. Here, the filtering flag for volcanoes of SZA < 70° is similarly applied.

A3 Ground-based measurements

For an inter-comparison with model results near the Earth's surface, ground-based measurements are used from three key sulfur-emitting regions: the United States of America (USA), Europe, and East Asia. These regions are selected because of their available and extensive datasets spanning a period of two decades, from 2000 to 2019. Next, a detailed explanation of the data is presented, with each region described separately.

Figure A1 Map of the USA showing positions of Eastern (red points) and Western (blue points) CASTnet measurement sites used in this study.

A3.1 USA

In the USA, data of various trace gases, including SO₂, are obtained from the Clean Air Status and Trends Network (CASTnet). This network (accessible at https://www.epa.gov/castnet, last access: 24 February 2024, Finkelstein et al., 2000), provides surface-level observations including monthly and yearly mean SO₂ concentrations and sulfur deposition fluxes over the USA. Dry deposition fluxes in this data set are derived with the Total Deposition Science Committee (TDep) Measurement Model Fusion method using the measured air concentrations from CASTnet and simualted CMAQ deposition velocities and fluxes. The wet deposition fluxes in this data set are estimted with the TDep Measurement Model Fusion method using the concentrations in precipitation and precipitation amounts measured by the National Atmospheric Deposition Program/National Trends Network (NADP/NTN). In this study the data between 2000 and 2019 are used. Given the large size of the USA's land surface, a categorization of analyzed SO₂ has been undertaken, distinguishing between Eastern and Western sites. A total of 89 sites have been chosen for this study, as they represent data for both, SO₂ concentrations and sulfur deposition fluxes, over the two-decade period. Among these, 29 observation sites positioned West of 100° W longitude represent the Western USA, while the remaining 60 sites East of 100° W are representative for the Eastern regions. The spatial distribution of these site locations is presented in Fig. A1.

In CASTnet, SO₂ and sulfate (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$) concentrations are directly measured on a weekly basis at each of the stations. The concentration of sulfur compounds is multiplied by the volume of precipitation to calculate the deposition fluxes. This calculation provides the amount of sulfur deposited per unit area over a specific time period (in this study this is expressed in $\mathrm{kg}\left(\mathrm{S}\right){\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$). However, measuring sulfur dry deposition fluxes faces some challenges, because it necessitates substantial instrumentation and technical resources (Hardacre et al., 2021). Therefore, deposition velocities are hourly estimated with the Multi-Layer Model (MLM, Meyers et al., 1998; Saylor et al., 2014) and are integrated with measured SO₂ concentrations, land usage, and meteorological data to obtain the SO₂ dry deposition flux. The deposition velocity in the Multi-Layer Model (MLM) is based on the aerodynamic resistance, the quasi-laminar resistance to transport, and the surface uptake resistance (Baumgardner et al., 2002).

A3.2 Europe

The observational data from Europe offers extensive long-term atmospheric SO₂ measurements, obtained from the European Monitoring and Evaluation Program (EMEP) since 1972. This data repository, accessible via the EMEP database (http://ebas.nilu.no/, last access: 27 February 2024; Tørseth et al., 2012), contains observations up to the present day. From EMEP, a total of 48 observational sites, distributed across Europe, are considered for our analyses. These sites not only monitor SO₂ concentrations, but also measure sulfate (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$) amount in precipitation samples (Aas et al., 2019) ranging from 2000 to 2019. The concentration of sulfur compounds in the precipitation is multiplied by the volume of precipitation to calculate the deposition fluxes. This calculation provides the amount of sulfur deposited per unit area over a specific time period (in this study this is expressed in $\mathrm{kg}\left(\mathrm{S}\right){\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$). Figure A2 visually illustrates the distribution of these observational sites across Europe.

Figure A2 Map of SO₂ measuring sites from 2000 to 2019 in Europe from the EMEP network.

SO₂ dry deposition data is, unlike to the CASTnet network, not available from the EMEP network. Consequently, the comparative analysis of sulfur deposition in Europe between observed data and the model results must rely solely on sulfate wet deposition from precipitation.

A3.3 East Asia

In the East Asia region fewer observational stations were available compared to Europe and the USA. 14 stations were selected from the Acid Deposition Monitoring Network in East Asia (EANET), comprising 2 urban, 3 rural, and 9 remote locations in China (specifically Southeastern China) and Japan, as depicted in Fig. A3. Notably, EANET stands out as the only network in East Asia equipped to monitor both, acid deposition and air pollution, with a particular emphasis on SO₂ (Ohizumi, 2023).

Figure A3 Map of SO₂ measuring sites situated in China and Japan from the EANET network.

Same as for European measurement stations, SO₂ dry deposition fluxes are neither measured nor simulated from Asian networks. Fortunately, EANET provided access to both, yearly mean SO₂ near-surface concentrations and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations in precipitation. The concentration of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ in the precipitation is multiplied by the volume of precipitation to calculate the deposition fluxes. This calculation provides the amount of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ deposited per unit area over a specific time period in mmol m⁻². However, here the wet deposition is converted to $\mathrm{kg}\left(\mathrm{S}\right){\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$, to be comparable to the caculated wet deposition fluxes in Europe and the USA.

A4 Post-processing of model and ground-based data for inter-comparison

For this comparison, the measurement locations are aggregated to the 2.8° × 2.8° horizontal grid of the RD1SD-base-01 EMAC simulation. Therefore, the nearest neighbor method is applied to assign the station measurement positions to the model grid. Measurements from multiple stations within one grid-box are averaged. Figure A4 presents an example showing the original positions of the ground-based measurement stations in Europe and after being aggregated onto the model grid.

Figure A4 Visual representation of the original sulfur measurement sites in Europe (white circles) between 2000 and 2019 from EMEP network and after being mapped onto EMAC model's grid (red squares).

Appendix B: Methodology for time series analyses

First, annual averages (2000–2019) for each measurement stations are calcualted from the monthly mean values. Next, these data are aggregated (with the nearest neighbor method) onto the model's grid for a direct comparison (see Sect. A4). For EMAC grid boxes containing multiple observational stations, the mean value of all stations within that box is computed, in order to obtain a single representative value per grid box. Afterwards, the weighted mean over the grid boxes within a specific region is determined by summing all the grid box values weighted by the area of the grid boxes, as follows:

$$\begin{array}{}\text{(B1)}& {\mathit{\mu }}_{\mathrm{w}}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N{w}_i{x}_i}{{\sum }_{i=\mathrm{1}}^N{w}_i}},\end{array}$$

where μ_w is the weighted mean, x_i represents the value in grid box i, w_i is the area (weight) of grid box i, and N is the total number of the grid boxes. Consequently, the weighted standard deviation σ_w is then expressed as:

$$\begin{array}{}\text{(B2)}& {\mathit{\sigma }}_{\mathrm{w}}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N{w}_i(x_i-{\mathit{\mu }}_{\mathrm{w}}{)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^N{w}_i}}}.\end{array}$$

For the spatial analysis, the calculated mean values and the corresponding standard deviations over the entire 20 years, are calculated. Trends between the years 2000 and 2019 are calculated by linear regression.