Evaluation of SO₂, SO₄²⁻ and an updated SO₂ dry deposition parameterization in the United Kingdom Earth System Model

In this study we evaluate simulated surface SO₂ and sulfate (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$) concentrations from the United Kingdom Earth System Model (UKESM1) against observations from ground-based measurement networks in the USA and Europe for the period 1987–2014. We find that UKESM1 captures the historical trend for decreasing concentrations of atmospheric SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ in both Europe and the USA over the period 1987–2014. However, in the polluted regions of the eastern USA and Europe, UKESM1 over-predicts surface SO₂ concentrations by a factor of 3 while under-predicting surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations by 25 %–35 %. In the cleaner western USA, the model over-predicts both surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations by factors of 12 and 1.5 respectively. We find that UKESM1’s bias in surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations is variable according to region and season. We also evaluate UKESM1 against total column SO₂ from the Ozone Monitoring Instrument (OMI) using an updated data product. This comparison provides information about the model's global performance, finding that UKESM1 over-predicts total column SO₂ over much of the globe, including the large source regions of India, China, the USA, and Europe as well as over outflow regions. Finally, we assess the impact of a more realistic treatment of the model's SO₂ dry deposition parameterization. This change increases SO₂ dry deposition to the land and ocean surfaces, thus reducing the atmospheric loading of SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$. In comparison with the ground-based and satellite observations, we find that the modified parameterization reduces the model's over-prediction of surface SO₂ concentrations and total column SO₂. Relative to the ground-based observations, the simulated surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations are also reduced, while the simulated SO₂ dry deposition fluxes increase.

1 Introduction

Anthropogenic sulfur dioxide (SO₂) emissions have been the main driver of the historical aerosol effective radiative forcing (ERF) since the mid-20th century (Boucher et al., 2013). SO₂ is emitted into the atmosphere from a number of anthropogenic (and natural) sources, and once in the atmosphere SO₂ can be oxidized to form sulfate (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$) aerosol, which plays a key role in acid deposition, atmospheric aerosol loading, and cloud properties, thereby directly influencing the Earth's radiative balance. For Earth system models (ESMs) to have a good representation of the historical climate and thereby give us confidence in their future projections, it is extremely important that they can capture the sulfur cycle. The United Kingdom Earth System Model (UKESM1), in common with other ESMs, has a cold bias in the mid-20th century which looks to be associated with an excessively negative aerosol ERF (Sellar et al., 2019; Seland et al., 2020). A key component of the analysis and development of UKESM1 focuses on the model's sulfur cycle and its link to historical aerosol forcing. Mulcahy et al. (2020) conducted an in-depth evaluation of the aerosol species in UKESM1 and its physical model component, HadGEM3-GC3.1, including ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, and uncovered some interesting differences in the sulfur budget between these two models, including differences in the SO₂ lifetimes and oxidant loading. We aim to extend their work by conducting a detailed evaluation of SO₂ and by probing deeper into the process-level uncertainty of the sulfur cycle.

Sources of SO₂ include industry, energy, land-based transport, shipping, volcanoes, biomass burning, and marine dimethyl sulfide (DMS) (Feng et al., 2020; Fioletov et al., 2016; Liu et al., 2018; Janssens-Maenhout et al., 2015; Crippa et al., 2018). Total global emissions of SO₂ increased to a peak value of approximately 180 Tg SO_x (as SO₂) yr⁻¹ in the 1970s, but following emission reduction policies to improve air quality and reduce acid deposition that were implemented in the 1980s (Hoesly et al., 2018), total global emissions had decreased to approximately 120 Tg SO_x yr⁻¹ by 2015 (Aas et al., 2019). This trend is captured in global models, but there is substantial temporal variation at the regional scale (Aas et al., 2019). Legislation has driven reductions in SO₂ emissions and subsequently ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ aerosol across Europe (Tørseth et al., 2012) and North America (Sickles and Shadwick, 2015; Holland et al., 1998). In these regions reductions in SO₂ emissions have had important environmental and health benefits as well as climate impacts. Turnock et al. (2015) found that between 1970 and 2010 surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ aerosol decreased by about 70 % in the observations and also in the simulations. For the same period, top of atmosphere (TOA) aerosol radiative forcing over this region increased by >3 W m⁻² in response to these changes in anthropogenic emissions. Similarly, Leibensperger et al. (2012) reported that over the USA aerosol radiative forcing decreased by 1.0 W m⁻² in the period from 1990 to 2010. Emission reduction policies in China have been implemented since 2013, which has reduced anthropogenic SO₂ emissions (Aas et al., 2019; Zheng et al., 2018; Hoesly et al., 2018; Liu et al., 2018; Krotkov et al., 2016) and subsequently driven decreases in aerosol optical depth (AOD) (Zhao et al., 2017). However, SO₂ emissions from India continue to increase (Aas et al., 2019; Liu et al., 2018; Krotkov et al., 2016).

Good representation of the sulfur cycle in models is essential for constraining uncertainties associated with the impacts of aerosols on the Earth system and thus understanding the global climate. The global atmospheric loading of SO₂ is controlled by the emissions (sources) into the atmosphere and the loss processes, which are oxidation to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, dry deposition, and wet deposition. Global-scale SO₂ emissions are represented in ESMs using emission inventories such as HTAP (Janssens-Maenhout et al., 2015), OMI-HTAP (Liu et al., 2018), EDGAR (Crippa et al., 2018), and CMIP6 (Feng et al., 2020), the latter being developed for use by models participating in the CMIP6 project (Eyring et al., 2016). Although uncertainty in SO₂ emissions is relatively low (Hoesly et al., 2018), in bottom-up inventories such as HTAP and EDGAR there may be uncertainty in the emission and activity factors and in the conversion from country scale to grid scale, and the input data may be incomplete or subject to rapidly changing economic and/or policy conditions (Janssens-Maenhout et al., 2015). In satellite-derived data sets there is uncertainty associated with the retrieval methods and the signal-to-noise ratio, which can make smaller sources and background concentrations more difficult to detect (Fioletov et al., 2016). Yang et al. (2019) have also found that injection height is a larger source of uncertainty in model representation of SO₂ emissions than inventory uncertainty, affecting surface concentrations by 70 %–130 % depending on sector and region, compared with 8 %–14 % from inventory uncertainty. The impact of injection height in UKESM1 was demonstrated by Mulcahy et al. (2020), who found that emitting SO₂ higher into the atmosphere rather than into the lowest model level increased the burden from 0.53 to 0.61 Tg and the lifetime from 2.08 to 2.21 d, although ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ was not significantly affected.

Anthropogenic emissions of SO₂ are generally from point sources such as power stations or smelters. Once emitted, SO₂ has a lifetime of approximately 2 d, although this can vary from 15 to 65 h in summer and winter respectively (Lee et al., 2011). The lifetime of SO₂ depends on both wet and dry deposition of the molecule and the oxidation rate to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$. The ∼2 d lifetime is such that much of the loss via oxidation and deposition occurs locally. SO₂ loss near sources and the impact of environmental conditions on loss processes have been investigated in a number of studies.

SO₂ deposition is highly dependent on the surface type, soil pH, solar radiation level, near-surface relative humidity, and, in particular, whether the underlying surface is wet or dry, with deposition increasing significantly for a wet surface. Wys et al. (1978) calculated diurnal averaged deposition of emitted SO₂ onto an agricultural field of ∼35 % within 300 km of the emission source, with daytime deposition significantly higher than at night. The same study found that ∼15 % of emitted SO₂ was dry deposited onto an arid desert surface within 300 km of the source. Studies over Europe indicate similar rates of deposition and sensitivity to surface type. For example, using flight observations off the eastern coast of the UK, Smith and Jeffrey (1975) estimated that ∼50 % of the SO₂ emitted from UK sources was removed from the atmosphere or converted to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ by the time it was observed in air parcels over the North Sea. This amounts to a loss or conversion of 50 % emitted SO₂ within ∼ 200–300 km of the emission source. Smith and Jeffrey (1975) further partitioned this loss into ∼ 30 %–35 % due to dry deposition and ∼ 10 %–15 % oxidation to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, with wet deposition making only a minor contribution to the total loss. Similar rates of SO₂ loss have been observed in a number of other observational studies of dry deposition (e.g. Payrissat and Beilke, 1975; Garland, 1977; Garland and Branson, 1977; Fowler, 1978; Erisman and Baldocchi, 1994). Studies analysing SO₂ dispersion around US power stations found that fractional oxidation rates to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ are sensitive to the amount of solar radiation, with rates ranging from a winter low of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}$ h⁻¹ to a summer high of $\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{2}}$ h⁻¹ (Altshuller, 1979; Meagher et al., 1983). Representing the SO₂ loss processes is challenging for ESMs because 200–300 km is represented by one to two grid cells, meaning that deposition and oxidation are parameterized on the model grid scale and may not capture temporal and spatial variation. In addition, there is uncertainty associated with the oxidation and deposition processes.

In the atmosphere SO₂ can be oxidized in the gas phase by hydroxyl (OH) radicals and in the aqueous phase by reactions in cloud water and rainwater involving hydrogen peroxide (H₂O₂), ozone (O₃), O₂ catalyzed by transition metal ions, and other oxidants to form ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ (see Turnock et al., 2019, and references within). The oxidation chemistry is necessarily simplified in many models due to the computational cost of detailed chemistry schemes, but studies have shown that oxidant levels can impact the lifetime of aerosol precursor species and ultimately global radiative forcing (Mulcahy et al., 2020; Karset et al., 2018). Uncertainty in aerosol radiative forcing also results from different values of cloud water pH, which alters ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ formation by changing the rate of aqueous-phase oxidation of SO₂ by ozone (Turnock et al., 2019). Observations have shown that cloud pH is both temporally and spatially variable (Aleksic et al., 2009; Murray et al., 2013; Schwab et al., 2016; Li et al., 2017), although measurements are very sparse. Typically, this variation is not accounted for in global chemistry climate models, including UKESM1 and its predecessor HadGEM3-GC3.1, both of which use a temporally and spatially constant cloud pH of 5.0. Turnock et al. (2019) found that increasing the cloud pH by 1.0 in HadGEM3-GC3.1 reduced the total SO₂ column by up to 50 % over Europe, North America, and East Asia for 1970–1974 and 2005–2009. The impact on ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ was variable due to the different SO₂ loadings over the different regions and in the different time periods. Overall aerosol radiative forcings varied by up to 4 W m⁻², with larger changes in some regions depending on whether cloud water pH was assumed to have increased or decreased over recent decades.

Loss of SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ to the Earth's surface by deposition can be through dry or wet processes. Dry deposition describes the removal of a gas or particle through direct contact of air with the Earth’s surface and wet deposition describes the incorporation of gases or particles into rain droplets or snow crystals and their subsequent removal through precipitation. Globally, dry deposition removes around 45 % of SO₂ from the atmosphere (Chin et al., 2000). The importance of dry deposition in the global sulfur budget is the reason why we target it for development in UKESM1. Dry deposition of SO₂ in ESMs is generally represented by a resistance in series approach (e.g. Archibald et al., 2020; Wu et al., 2020). Deposition of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ is mainly via wet processes (approximately 90 %, Chin et al., 2000), including nucleation scavenging within the cloud (rain out) and impact scavenging below the cloud (wash out), but dry deposition of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ does occur through gravitational settling. Deposition processes are necessarily parameterized in global models because they occur at sub-grid scales, and this contributes to model uncertainty. Further, observational flux data sets are sparse and frequently temporally and spatially limited, hindering model evaluation of deposition processes at regional to global scales.

Sulfur species are relatively well observed compared with many atmospheric components as their role in air pollution is well established. In the 1970s and 1980s the increasingly detrimental impacts of rising SO₂ emissions on acid deposition, air quality, and human health in Europe and North America led to monitoring networks being set up in these regions (Tørseth et al., 2012; MACTEC-Engineering and Consulting, 2005). Rising pollution in Asia also led to the establishment of the Acid Deposition Monitoring Network in East Asia (EANET) in 2001 (e.g. Wang et al., 2008). However, even with these data sets it is only possible to evaluate model simulations of the recent historical period, and similar data sets are not available for other large source regions such as India, the Middle East, or remote regions. Further, the lack of a range of measurements, including flux observations, hinders detailed process studies at large scales. Since the early 2000s satellite observations of near-surface SO₂ have also become available. Of these, the satellite data sets with the best temporal resolution and spatial coverage for SO₂ are from the Ozone Monitoring Instrument aboard the NASA Earth-observing system Aura spacecraft (Fioletov et al., 2016). Although biases in the SO₂ retrieval from OMI limit its use at high and low latitudes in winter and over areas with low atmospheric SO₂ loading, they do provide valuable information over regions where there are no long-term or even any ground-based observations (Li et al., 2020; Levelt et al., 2018).

This paper is configured as follows: the model (UKESM1), the model simulations, observation data sets, and modifications to UKESM1's SO₂ dry deposition parameterization are described in Sect. 2. In Sect. 3 we evaluate UKESM1 against observations of surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and total column SO₂. In Sect. 4 we assess the impact of the modifications to UKESM1's SO₂ dry deposition parameterization. The discussion and conclusions are presented in Sects. 5 and 6.

2 Methods

2.1 UKESM1

UKESM1 is the latest-generation Earth system (ES) model developed in the UK. UKESM1 has HadGEM3-GC3.1 (Kuhlbrodt et al., 2018; Williams et al., 2018) as its physical–dynamical core. HadGEM3-GC3.1 is comprised of the Global Atmosphere 7.1 (GA7.1) configuration of the Unified Model (UM) (Walters et al., 2019; Mulcahy et al., 2018), the Nucleus for European Modelling of the Ocean (NEMO) model (Storkey et al., 2018), the Los Alamos Sea Ice Model (CICE, Ridley et al., 2018), and the Joint UK Land Environment Simulator (JULES) land surface model (Best et al., 2011). The additional ES process models include the stratospheric–tropospheric (StratTrop) version of the United Kingdom Chemistry and Aerosol (UKCA) model (Archibald et al., 2020), the Model of Ecosystem Dynamics, nutrient Utilisation, Sequestration and Acidification (MEDUSA, Yool et al., 2013), and the terrestrial biogeochemistry component of JULES (Clark et al., 2011). UKESM1 is described in detail, along with its component models and the coupling between them, by Sellar et al. (2019). The aerosol scheme used in UKESM1 (GLOMAP-Mode, Mann et al., 2010), including SO₂ emissions and chemistry, is described in detail by Mulcahy et al. (2020). In UKESM1 the land and atmosphere share a regular latitude–longitude grid with a resolution of 1.25^∘ × 1.875^∘ (approximately 135 km at the mid latitudes). There are 85 vertical levels on a terrain-following hybrid height coordinate with a model lid at 85 km above sea level and 50 of these levels below 18 km. The ocean has a horizontal resolution of 1^∘ and 75 vertical levels. While the atmospheric time step of the model physics is 20 min, due to the inherent computational cost of the chemistry and aerosol components, both of these components are called once per hour.

In UKESM1 the SO₂ emissions including anthropogenic sources are from the CMIP6 inventory (Feng et al., 2020). Large explosive volcanic sources and biomass burning sources are not interactively modelled but are prescribed using the CMIP6 stratospheric aerosol climatology (Sellar et al., 2019) and van Marle et al. (2017) emission inventory respectively. Continuously degassing volcanic sources are also included as present-day, three-dimensional, temporally fixed (i.e. no seasonal variation) fields (Dentener et al., 2006). Emissions from the energy and industrial sectors are all emitted into the first model layer. We summarize how loss of SO₂ from the atmosphere via oxidation and wet and dry deposition is modelled here, but for a detailed description of these processes in UKESM1 the reader is referred to Archibald et al. (2020) and Mulcahy et al. (2020). Gas- and aqueous-phase oxidation of SO₂ to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ is represented by the reactions shown in Table 1 (Pham et al., 1996; Sander et al., 2003; Kreidenweis et al., 2003). Dry deposition of SO₂ is parameterized following the resistance in series approach originally developed by Wesely (1989) (see Sect. 2.2.1). Loss via wet deposition is the SO₂ that is scavenged and subsequently converted to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ in rainwater. It is parameterized as a first-order loss rate, calculated as a function of UKESM1's three-dimensional convective and large-scale precipitation (Archibald et al., 2020; O'Connor et al., 2014). Sulfate aerosol is also removed from the atmosphere by dry and wet deposition (Mulcahy et al., 2020). The aerosol dry deposition and sedimentation are represented by a resistance in series approach similar to that used for gaseous species but which also accounts for aerosol size (Mann et al., 2010). Wet deposition is parameterized in UKESM1 by an in-cloud convective plume-scavenging scheme following the approach described by Kipling et al. (2013) and by nucleation scavenging (Mulcahy et al., 2020).

Table 1 Summary of SO₂ oxidation chemistry in UKESM1.

2.2 SO₂ dry deposition parameterization

The UKESM1 parameterization of SO₂ dry deposition follows that described in Wesely (1989). This scheme uses the widely accepted approach of calculating the flux of a depositing gas as a function of a deposition velocity multiplied by the concentration gradient of the gas between a reference height (z, e.g. the lowest model level) and the receptor surface (Eq. 1). The deposition velocity is calculated by analogy with electrical resistance and is inversely proportional to three resistances to deposition, representing the three stages of gaseous transport to a receptor surface. These are (i) aerodynamic resistance (R_a) to gas transport through the near-surface turbulent layer, (ii) viscous resistance to gas transfer across a quasi-laminar layer surrounding the receptor surface (R_b), and (iii) structural resistance to deposition of the receptor surface itself (R_c). For a detailed description of this approach, see Wesely (1989), Erisman and Baldocchi (1994), and Zhang et al. (2003). The deposition velocity is calculated for each fractional surface type in a given model grid box, as is a resulting loss rate (flux) of SO₂ from the atmosphere. The loss rates to each fractional surface are combined, resulting in the total loss of SO₂ from the model atmosphere due to dry deposition. The deposition velocity is given by Eq. (2). If the surface is covered by vegetation, R_a is generally calculated at a zero-plane displacement height $z=z-d$, where d is usually 0.6–0.8 times the vegetation height in metres. The UKESM1 calculations of R_a and R_b follow standard approaches (see Eqs. 3 and 4). The aerodynamic resistance, R_a, is calculated from the wind profile taking into account atmospheric stability and the surface roughness, where z₀ is the roughness length, ψ is the Businger dimensionless stability function, κ is von Karman's constant, and u_* is the friction velocity. The quasi-laminar sub-layer resistance, R_b, is calculated with Sc the Schmidt number and Pr the Prandtl number.

The surface or canopy resistance to deposition, R_c, is the most difficult of the three resistances to parameterize as it is sensitive to biochemical details of the individual receptor surfaces. R_c is typically a function of the following receptor-specific resistances: (i) canopy stomatal resistance (R_stom) combined with the mesophyll resistance (R_m) of a given plant, (ii) canopy cuticle or external leaf resistance (R_cut), and (iii) soil resistance (R_soil), combined with an in-canopy resistance (R_inc), describing the turbulent transport of a gas through the plant foliage to the ground. The stomatal resistance, leaf cuticle resistance, and soil resistance are assumed to operate in parallel. For surfaces not covered with vegetation (e.g. open water, bare soil, or snow-covered surfaces), R_c is made equal to one of R_water, R_soil, and R_snow. The receptor-specific resistances are combined as shown in Eq. (5) to calculate R_c. In UKESM1, R_stom follows the approach outlined in Wesely (1989), based on the original work of Baldocchi et al. (1987). R_stom is first calculated for water vapour for each vegetation type, and R_stom for other gases is then derived by scaling R_stom for water vapour by the ratio of the diffusion coefficient for the gas in question and that of water vapour. Due to a general lack of knowledge, R_m values are assumed to be zero for all gases. In UKESM1, R_inc and R_soil are combined into a single value (referred to hereafter as R_soil). In UKESM1, R_c for SO₂ for the 13 fractional land cover types is initially set to the standard surface resistance (R_surf) values given in Table A2.

$$\begin{array}{}\text{(1)}& {\displaystyle }F=V_d\times C\text{(2)}& {\displaystyle }{V}_d={\displaystyle \frac{\mathrm{1}}{R_{\mathrm{a}}+R_{\mathrm{b}}+R_{\mathrm{c}}}}\text{(3)}& {\displaystyle }{R}_{\mathrm{a}}=\left(\ln \left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right)\mathrm{-}\mathit{\psi }\right)/\left(\mathit{\kappa }{u}_{*}\right)\text{(4)}& {\displaystyle }{R}_{\mathrm{b}}=(\mathit{Sc}/\mathit{Pr}{)}^{\frac{\mathrm{2}}{\mathrm{3}}}/(\mathit{\kappa }{u}_{*})\text{(5)}& {\displaystyle }{R}_{\mathrm{c}}={\left[{\displaystyle \frac{\mathrm{1}}{R_{\mathrm{stom}}+R_{\mathrm{m}}}}+{\displaystyle \frac{\mathrm{1}}{R_{\mathrm{inc}}+R_{\mathrm{soil}}}}+{\displaystyle \frac{\mathrm{1}}{R_{\mathrm{cut}}}}\right]}^{-\mathrm{1}}\end{array}$$

2.2.1 Modifications to UKESM1's SO₂ dry deposition parameterization

In this study we investigate two changes to the SO₂ dry deposition parameterization in UKESM1. Firstly, we account for a key omission in UKESM1 in that for R_cut and R_soil no account is taken as to whether the receptor surface is wet or dry nor of the near-surface relative humidity. Observational studies suggest that SO₂ dry deposition (through a decrease in R_c) is significantly more efficient over wet surfaces compared with dry surfaces as well as for increasing values of near-surface relative humidity due to the high solubility of SO₂ in water (e.g. Garland and Branson, 1977; Fowler, 1978; Erisman and Baldocchi, 1994; Erisman et al., 1994). We apply the findings from these studies to extend the calculation of R_cut for SO₂ in UKESM1 to be a function both of whether the model vegetation is wet or dry and of the near-surface relative humidity. This change allows a surface to remain wet after rainfall for a period of 3 h, where previously it would have been “dry” immediately after the rainfall event. R_soil for SO₂ is also made a function of near-surface relative humidity. These changes are referred to as R_surf-mod and will impact SO₂ dry deposition over land surfaces. We include a more detailed description of the modifications to UKESM1's SO₂ parameterization in Appendix A. Secondly, we change the surface resistance term for SO₂ dry deposition to water (R_water) from an erroneously high value of 148 to 1 s m⁻¹ to better reflect the high solubility of SO₂ in water. While lower than the value of 20 s m⁻¹ used by Zhang et al. (2003), it reflects the small, observed value of 0.004 s m⁻¹ from Garland (1977). This change is referred to as R_water-mod and will impact SO₂ dry deposition predominantly over the ocean.

In addition to the primary changes to the SO₂ dry deposition parameterization (R_water-mod and R_surf-mod), we also include two secondary modifications. These are (1) an update in the calculation of the stability parameter ($z/L$) to better describe dry deposition under very stable atmospheric conditions and (2) a bug fix in the DMS chemistry. The stability parameter ($z/L$) describes the flux profile relationship and is important for calculating R_a in Eq. (3). Note that the Monin–Obukhov length (L) is derived locally in the UKCA code using local values of air density, temperature, and friction velocity, where the friction velocity is computed in the UM turbulence scheme and so is consistent across subroutines. Here we update the calculation of the stability parameter from that given by Dyer (1974) to that described by Holtslag and Bruin (1988). We also reduce the reference height for dry deposition (z) from 50 to 10 m. The reference height is the height below which there is no turbulence under very stable conditions and is also important for calculating R_a. Following Ganzeveld and Lelieveld (1995), the reference height should be half the average height of the lowest model layer, which in UKESM1 is 20 m. The changes to $z/L$ and z act to reduce the rate at which the deposition velocity decreases under very stable conditions, although we note that there is also an impact on the calculation of aerodynamic resistance under unstable conditions. The DMSO bug fix corrects the equation for DMS oxidation by OH (see Reaction R1) in UKCA's StratTrop mechanism, where the products incorrectly contain more sulfur atoms than the reactants. We substitute Reaction R1 with Reactions R2 and R3. This reduces the SO₂ yield to a maximum of 0.84, which may be further reduced as DMSO deposits to the Earth's surface. However, the changes in simulated SO₂ are actually only of the order of 1 % because anthropogenic sources are not affected by this change. Although the secondary changes incorporate important updates into the model, their impact on the atmospheric SO₂ loading in UKESM1 is small in comparison to that driven by R_water-mod and R_surf-mod, and we do not discuss it here.

$$\begin{array}{}\text{(R1)}& {\displaystyle }{\mathrm{C}}_{\mathrm{2}}{\mathrm{H}}_{\mathrm{6}}\mathrm{S}+\mathrm{OH}\to {\mathrm{SO}}_{\mathrm{2}}+{\mathrm{CH}}_{\mathrm{3}}{\mathrm{SO}}_{\mathrm{2}}\mathrm{OH}\text{(R2)}& {\displaystyle }{\mathrm{C}}_{\mathrm{2}}{\mathrm{H}}_{\mathrm{6}}\mathrm{S}+\mathrm{OH}\to \mathrm{0.6}{\mathrm{SO}}_{\mathrm{2}}+\mathrm{0.4}{\mathrm{C}}_{\mathrm{2}}{\mathrm{H}}_{\mathrm{6}}{\mathrm{SO}}_{\mathrm{2}}+{\mathrm{CH}}_{\mathrm{3}}{\mathrm{O}}_{\mathrm{2}}\text{(R3)}& {\displaystyle }{\mathrm{C}}_{\mathrm{2}}{\mathrm{H}}_{\mathrm{6}}{\mathrm{SO}}_{\mathrm{2}}+\mathrm{OH}\to \mathrm{0.6}{\mathrm{SO}}_{\mathrm{2}}+\mathrm{0.4}{\mathrm{CH}}_{\mathrm{3}}{\mathrm{SO}}_{\mathrm{2}}\mathrm{OH}.\end{array}$$

2.3 Model simulations

For this evaluation we initially use four simulations from the 19-member ensemble of historical simulations that were conducted for UKESM1's contribution to CMIP6 (Sellar et al., 2019; Tang et al., 2019). The historical simulations cover the period from 1850 to the end of 2014, thus modelling the evolution of climate and composition since the pre-industrial era. These simulations are forced by transient external forcings of solar variability, land use, well-mixed greenhouse gases, and other trace gas emissions and aerosols. The volcanic forcing due to the stratospheric injection of SO₂ from volcanic eruptions is prescribed as a zonal mean climatology of the stratospheric aerosol optical properties over the historical period. All forcings and how they are implemented in UKESM1 are described fully in Sellar et al. (2019). Each historical ensemble member was initialized from a different date in the pre-industrial control simulation (Yool et al., 2020). We use monthly mean output for surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations and SO₂ dry deposition flux. We use a second four-member ensemble of historical simulations to evaluate the impact of the changes to the SO₂ dry deposition parameterization described in Sect. 2.2.1; hereafter, this ensemble is referred to as UKESM1-SO2. The UKESM1-SO2 historical simulations are set up and run as for the UKESM1 historical simulations. To calculate the detailed SO₂ budget, we utilize the atmosphere-only configuration of UKESM1 and UKESM1-SO2. This is referred to as the Atmospheric Model Intercomparison Project (AMIP) configuration and allows us to generate the diagnostics required for the budget analysis that were not output in the historical simulations at a reduced computational cost. The UKESM1 AMIP configuration is driven by observed sea surface temperature (SST) and sea ice. It does not include the additional dynamic ocean and land surface components (Eyring et al., 2016). Instead, the required vegetation (vegetation fractions, leaf area index, canopy height) and surface ocean biology fields (DMS and chlorophyll) are taken from a single UKESM1 historical member and are prescribed as ancillary data, thereby maintaining traceability to the fully coupled model. For the SO₂ budget calculations AMIP simulations were run from 1979 to the end of 1983.

Table 2 Summary of model configurations used in this study.

2.4 Ground-based observations

We compare the modelled surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations to observations from the Clean Air Status and Trends Network (CASTNet, https://www.epa.gov/castnet, last access: 14 December 2021, Finkelstein et al., 2000) and the European Monitoring and Evaluation Program (EMEP, http://ebas.nilu.no/, last access: 14 December 2021; Tørseth et al., 2012). CASTNet provides surface observations of mean seasonal SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations which are available from 1987 to the present at 97 sites situated in the USA. In this study we used observations from the CASTNet sites designated as “western reference” or “eastern reference”. The reference sites have been reporting measurements since at least 1990 and are used for determining long-term trends (e.g. Clarke et al., 1997; Holland et al., 1998; MACTEC-Engineering and Consulting, 2005; Baumgardner et al., 2002). There are 16 western reference sites and 33 eastern reference sites which are located in the continental USA to the west and east of 100^∘ W respectively (MACTEC-Engineering and Consulting, 2005). The eastern region is significantly more polluted than the western region due to the larger number of SO₂ sources there. We therefore keep the western and eastern data sets separate to assess how UKESM1 performs in the two regions. Hereafter, we refer to the eastern and western USA regions as USA–E and USA–W respectively. For this evaluation we used the mean seasonal surface concentrations for SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, which are measured with filter pack samplers at weekly sampling intervals. Details of the quality control procedures and of how the mean seasonal concentrations are calculated are given in Baumgardner et al. (2002).

We also evaluate simulated SO₂ dry deposition flux from UKESM1 against observations from CASTNet, using the same eastern and western reference sites that were used to evaluate surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations. The CASTNet deposition fluxes are derived using modelled deposition velocities rather than directly measured fluxes, which are difficult to obtain due to the requirement for extensive instrumentation and technical resources. Direct measurements of SO₂ dry deposition flux are therefore temporally and spatially limited and not suitable for evaluating long-term trends. To derive the SO₂ dry deposition fluxes, measurements of SO₂ concentration are combined with routine meteorological measurements, information on the land use type, and LAI at the measurement site. These data are then combined with modelled deposition velocities from the Multi Layer Model (MLM, Meyers et al., 1998; Saylor et al., 2014). The methodology used to derive the SO₂ dry deposition fluxes for CASTNet is described in Clarke et al. (1997) and Baumgardner et al. (2002). While the modelled SO₂ dry deposition fluxes can be under-predicted by approximately 30 % (Clarke et al., 1997), it is considered to be the best available approach to regional-scale assessment of dry deposition (Finkelstein et al., 2000; Baumgardner et al., 2002; Sickles and Shadwick, 2007a, b). This approach has been used to determine SO₂ dry deposition fluxes for CASTNet since 1987 (e.g. Clarke et al., 1997; Baumgardner et al., 2002) and to assess global- and regional-scale models (Vet et al., 2014; Tan et al., 2018; Tang et al., 2018).

Surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations have been monitored at EMEP sites for the period 1972–present (Tørseth et al., 2012). In this study we have used observations of surface concentrations of SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ from 48 and 42 sites respectively. We have selected sites where there are at least 10 years of continuous measurements and with a few exceptions have used sites where SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ were co-located. We use monthly mean observations for both species. No SO₂ dry deposition data were available from EMEP. The locations of the CASTNet and EMEP sites used in this study are shown in Fig. 1.

Figure 1 Map of the locations of the CASTNet and EMEP measurement sites used in this study.

2.5 Data processing

For this evaluation we calculated seasonal averages for the modelled surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations and for the EMEP observational data. The seasonal periods were defined as December–January–February (DJF), March–April–May (MAM), June–July–August (JJA), and September–October–November (SON). The CASTNet data for all variables were available as seasonal averages for these periods. Model grid cell output was co-located with the CASTNet and EMEP measurement sites. In some cases, this resulted in model data from a particular grid cell being compared with more than one measurement site. For the time series analysis, regional means and standard deviations were calculated across the sites in the USA–E, USA–W, and Europe regions. Although there is spatial variation in the surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations across Europe, for example concentrations are relatively low in Scandinavia but are much higher in south-eastern Europe, it is less easy to classify “clean” and “polluted” regions at the global model scale. Therefore we classify Europe as a single region. For the spatial analysis and calculation of time series statistics, we calculate mean values over the whole time series, i.e. 1987–2014 and for two time slices at the start (1990–1995) and end of the time series (2009–2014). For the time slices we only used sites that had at least 3 out of 5 years of data available. We investigate the two time slices to assess the model's performance during the different pollution levels at the start and end of the time series. We determine the rate of change (trend) in the surface concentrations by calculating the linear regression for 1987–2014, 1990–1995, and 2009–2014.

2.6 Satellite observations

Total column SO₂ (TCSO₂) measurements came from the Ozone Monitoring Instrument (OMI) and were obtained from the Goddard Earth Sciences Data and Information Services Centre (https://aura.gesdisc.eosdis.nasa.gov/data/Aura_OMI_Level2/OMSO2.003, last access: 14 December 2021) (Li et al., 2020). OMI is situated on board NASA's polar-orbiting Aura satellite launched in 2004 with a local overpass time of approximately 13:45 local solar time (LST). OMI has a nadir footprint of 13 km × 24 km and a spectral viewing range of 270 to 500 nm (Levelt et al., 2018). The TCSO₂ product is quality controlled for cloud radiation fraction > 0.0 and < 0.5, solar zenith angle < 65^∘, the South Atlantic Anomaly flag = 0, ice cover flag = 0, the air mass factor (AMF) > 0.3, and TCSO₂ > −1.0 Dobson unit (DU). Background TCSO₂ average values tend to be positive near-zero quantities (i.e. just above 0.0), where some soundings are slightly negative. If only positive TCSO₂ were incorporated into the background averages, this would positively skew the true value.

For a robust comparison between model simulations and satellite data, both data sets typically require spatiotemporal co-location to reduce sampling (representation) errors. To achieve this, high temporal resolution (e.g. 3-hourly or 6-hourly) model outputs of three-dimensional tracer and pressure fields are required over the analysis period to capture e.g. diurnal variability (Pope et al., 2016; Monks et al., 2017). However, this is difficult when using standard climate model simulations, including those used in this study, which typically output monthly means due to their long-term climate focus. In this comparison we performed tests to show that using the monthly mean model output was suitable given the relatively uniform diurnal cycle of SO₂ emissions. For these tests we made initial comparisons of model output and satellite data using 6-hourly and monthly mean output from the same UKESM1 simulation. We found that the temporal sampling of the model was not overly critical for SO₂, i.e. that modelled SO₂ has a sufficiently long lifetime to dampen the influence of diurnal sampling of the model. Further details on the TCSO₂ product, how it was processed to obtain TCSO₂ values, and the assessment of the temporal resolution are given in Pope and Chipperfield (2021).

3 Evaluation of trends and biases in modelled SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations

3.1 Time series analysis of surface concentrations of SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$

UKESM1 simulations of surfaces of SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations are compared with observations from the CASTNet and EMEP networks for the period 1987–2014 in Fig. 2. The statistics summarizing the model bias and trends over this period as well as for the 1990–1995 and 2009–2014 time slices are shown in Tables 3 and 4. We find that UKESM1 captures the historical reduction in surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations. This is in agreement with Aas et al. (2019), who reported that an ensemble of global aerosol models was generally able to capture the recent historical declines in these two species over the USA and Europe. UKESM1 over-predicts surface SO₂ concentrations in all three regions, but the direction of the model's bias in surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations is spatially variable.

Figure 2 Time series of observed and modelled mean annual surface SO₂ (a–c) and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ (d–f) concentrations for USA–W (a, d; N sites = 16), USA–E (b, e; N sites = 33), and Europe (c, f; N sites = 48 for SO₂ and 42 for ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$). Each point in the time series represents the mean across the measurement sites in the region. Note that the vertical scale for SO₂ (a–c) is a factor of 6 larger than that for ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ (d–f).

Figure 2a–c show that in both Europe and USA–E the model over-predicts surface SO₂ concentrations, particularly at the start of the time series, which then decrease too rapidly. Over Europe, the observed surface SO₂ concentrations decrease at a rate of 0.72 µg m⁻³ yr⁻¹ in the period 1990–1995 which slows to 0.38 µg m⁻³ yr⁻¹ by 2009–2014. However, the modelled surface SO₂ concentrations decrease by 2.52 µg m⁻³ yr⁻¹ for 1990–1995, slowing to only 1.31 µg m⁻³ yr⁻¹ by 2009–2014. Over USA–E the observed surface SO₂ concentrations decrease at similar rates at the start and end of the time series. UKESM1 is better able to capture the trend at the start of the time series but simulates an overly rapid reduction in surface SO₂ concentrations after 2005 (see Fig. 2b). Over USA–E UKESM1 simulates the sharp drop in surface SO₂ concentrations that occurred in 1995 following the implementation of Phase 1 of the USA's Clean Air Act Amendments (McHale et al., 2021). However, the model then simulates relatively high surface SO₂ concentrations for the period 1996–1999 rather than the sustained lower surface SO₂ concentrations that are observed after 1995. Over USA–W the observed surface SO₂ concentrations remain steady from 1987 to 2014 due to there being fewer sources in this region. Figure 2a shows that UKESM1 simulates the steady surface SO₂ concentrations at the start of the time series, albeit with a positive bias. However, after 1995 UKESM1 simulates decreasing surface SO₂ concentrations in USA–W, which brings the modelled values into better agreement with the observations but introduces an artificial trend into the modelled time series.

UKESM1 over-predicts the annual mean surface SO₂ concentrations in the polluted regions of Europe and USA–E by a factor of 3.2 to 3.4 over the period 1987–2014, although the absolute bias is higher USA–E (see Table 3). While the absolute magnitude of the bias in mean annual surface SO₂ concentration is less in USA–W compared with the polluted regions, proportionally it is much larger, with the model simulating surface SO₂ concentrations by more than 10 times the observed values. The absolute magnitude of the bias in mean annual surface SO₂ concentration decreased from 1990–1995 to 2009–2014 in all three regions (see Table 3), reflecting the model's more rapid decrease in surface SO₂ concentrations relative to the observations. However, the normalized mean bias (NMB) values were slightly higher in 2009–2014 compared with 1990–1995.

Table 3 Statistics for mean annual surface SO₂ concentrations at USA–W, USA–E, and Europe. The mean and trend values are in µg m⁻³ and µg m⁻³ yr⁻¹ respectively.

Table 4 Statistics for mean annual surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations at USA–W, USA–E, and Europe. The mean and trend values are in µg m⁻³ and µg m⁻³ yr⁻¹ respectively.

Figures 2e–f show that UKESM1 captures both the magnitude and trends in surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations better than the surface SO₂ concentrations. The model simulated surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations decreasing at a rate of 0.13 µg m⁻³ yr⁻¹ (USA–E) and 0.09 µg m⁻³ yr⁻¹ (Europe) compared with the observed trend of ≈0.10 µg m⁻³ yr⁻¹ (see Table 4). UKESM1 does under-predict mean annual surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations in the polluted regions of USA–E and Europe, but the model bias is relatively small compared with the large over-prediction of mean annual surface SO₂ concentration (see Tables 3 and 4). We also find that there is a large range associated with the modelled and observed data and that the mean surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations lie within these ranges. The model bias remained relatively constant over the period from 1987 to 2014, ranging from −0.96 to −0.80 µg m⁻³ for the periods 1990–1995 and 2009–2014 over USA–E and from −0.91 to −0.69 µg m⁻³ for the same periods over Europe (see Table 4). The picture is different in USA–W, where, in contrast to USA–E and Europe, UKESM1 over-predicts mean annual surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentration by an average of 150 % for 1987–2014. Both the absolute model bias and the NMB are worse in 1990–1995 than in 2009–2014, which may be attributed to the model simulating a much faster decrease in mean annual surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations compared with the observed trend for the later period (see Table 4).

3.2 Spatial evaluation of surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations

Figures 3 and 4 show the spatial distribution of modelled and observed mean annual surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations for the periods 1990–1995 and 2009–2014 for the USA and Europe. We find that UKESM1 captures the spatial distribution of surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations over each region, simulating higher concentrations in USA–E, central Europe, and eastern Europe, where there are numerous large sources, and lower concentrations in USA–W and northern Europe, which have far fewer sources (see Fig. B1). These figures also show the localized versus dispersed nature of the surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations, with high SO₂ concentrations located within two to three grid boxes (200–400 km) of the emission sources (see Fig. B1), while ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ is distributed more widely. Figure 3 shows that modelled and observed surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations across the USA are lower in 2009–2014 compared with 1990–1995, demonstrating the widespread impact of emission reduction policies. However, the disparity between higher concentrations in USA–E and lower concentrations in USA–W is still apparent for both species in the later period.

Figure 3 Mean annual surface SO₂ concentration (a, c) and surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentration (b, d) for 1990–1995 (a, b) and 2009–2014 (c, d) for modelled output. Observations from 49 CASTNet measurement sites are plotted as black-edged circles on the same colour scale.

Figure 4 Mean annual surface SO₂ concentration (a, c) and surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentration (b, d) for 1990–1995 (a, b) and 2009–2014 (c, d) for modelled output. Observations from EMEP measurement sites (N sites = 48 for SO₂ and 42 for ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$) are plotted as black-edged circles on the same colour scale.

In Europe the highest mean annual surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations were observed in central and eastern Europe and the south-east (SE) of England (see Fig. 4). Lower concentrations of both species were observed in northern and western regions, e.g. Scandinavia and the western coast of Ireland. Figure 4c shows that mean annual surface SO₂ concentrations were generally lower in 2009–2014 compared with 1990–1995, especially in central and eastern Europe, due to the impact of air quality legislation. However, for 2009–2014 modelled and observed levels of SO₂ remain high in the south-eastern region (see Fig. 4c). UKESM1 reproduces the spatial distribution of mean annual SO₂ concentrations across Europe but has large positive biases over most of the region. The largest model biases were in eastern and south-eastern Europe during the period 1990–1995, where UKESM1 simulates mean annual surface SO₂ concentrations of up to 100 µg m⁻³ compared with observed values of 10–30 µg m⁻³. Figure 4b shows that UKESM1 captures the spatial distribution between low mean annual surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations in northern and western regions of Europe and high mean annual surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations in central Europe. However, the model under-predicts mean annual surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations at many locations across Europe, with the largest biases occurring across Denmark and the regions surrounding the Baltic Sea.

3.3 Spatial distribution of model bias in surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations

Figures 5 and 6 show that the direction of the model biases, whether positive or negative, is generally consistent across a region. UKESM1 over-predicts mean annual surface SO₂ concentration for all sites in both the USA–E and USA–W regions and, with the exception of some Scandinavian sites, across Europe. Figures 5 and 6 also show that mean annual surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations are generally under-predicted across the USA–E and European sites while being over-predicted at the USA–W sites.

Figure 5 Geographic distribution of mean bias (UKESM1 – obs) in mean annual surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentration at 49 CASTNet measurement sites. The mean annual surface concentrations are calculated over the period 1987 to 2014. Absolute mean bias (MB) is shown in (a) and (b) and normalized mean bias (NMB) is shown in (c) and (d). Note that different scales are used for the SO₂ model bias (a) and normalized mean bias (c).

Figure 6 Geographic distribution of mean bias (UKESM1 – obs) in mean annual surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentration at EMEP measurement sites (N sites = 48 for SO₂ and 42 for ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$). The mean annual surface concentrations are calculated over the period 1987 to 2014. Absolute mean bias is shown in (a) and (b) and NMB is shown in (c) and (d). Note that different scales are used for the SO₂ model bias (a) and normalized mean bias (c).

The model's over-prediction of mean annual surface SO₂ concentration is largest close to the sources. For example, UKESM1 over-predicts surface SO₂ concentrations by up to 50 µg m⁻³ in the central USA–E area but only by around 10 µg m⁻³ at the surrounding sites, whilst in Europe the largest biases of around 10–20 µg m⁻³ tend to be in central and eastern areas. The model's tendency to over-predict SO₂ concentrations to a greater extent close to the sources is also shown in the plots of NMB (see Figs. 5c and 6c) and is likely why there is a larger range in the modelled mean annual surface SO₂ concentrations averaged across the USA–E and European measurement sites compared with the observational values (see Fig. 2b and c). Figures 5b and 6b show that UKESM1 generally under-predicts surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentration across the USA–E and European sites, with model biases of −1 to −3 µg m⁻³. We find that in USA–E the largest negative biases in surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations are not necessarily co-located with the largest positive biases in SO₂, instead occurring at sites several hundred kilometres from the large point sources. Similarly, in Europe the largest biases in surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations occur further north than the large biases in surface SO₂ concentration, and at certain sites located near large sources, UKESM1 over-predicts surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations. The plots of NMB show that there is less spatial variation in the model bias for surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations, reflecting UKESM1's ability to capture the more distributed nature of atmospheric ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$. In USA–W UKESM1 over-predicts both SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations at almost all of the sites (see Fig. 5). For SO₂ this is a consequence of the sparsely distributed measurement sites being located in rural regions remote from any of the sources in USA–W (Clarke et al., 1997; MACTEC-Engineering and Consulting, 2005; Baumgardner et al., 2002) (see also Fig. B1). This results in some very large NMB values in USA–W (see Fig. 5a).

3.4 Seasonal cycles

Figure 7 shows modelled and observed surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations averaged seasonally for the period 1987–2014. The comparison between the model and observations is summarized for DJF and JJA in Table 5. The higher wintertime SO₂ concentrations are driven by greater emissions from coal-fired power plants and domestic heating and less oxidation. Conversely, there are fewer emissions and higher oxidant concentrations in summertime. These cycles drive correspondingly low ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations in winter and high ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations in summer. Overall, we find that the model bias in surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations depends on the season as well as the region and pollution levels.

Figure 7 Modelled and observed seasonal mean surface SO₂ concentration (a–c) and surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentration (d–f) for the period 1987–2014. USA–W (a, d; N sites = 16); USA–E (b, e, N sites = 33), and EMEP (c, e; N sites (SO₂) = 48 and N sites (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$) = 42). The blue-shaded region and the black error bars represent the standard deviation across the sites in the observational network.

Table 5 Statistics for seasonal mean surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations (µg m⁻³) at USA–W, USA–E, and Europe. The mean seasonal values are averaged over the period 1987–2014.

The results presented in Sect. 3.1 and 3.2 show that UKESM1 consistently over-predicts mean annual surface SO₂ concentrations in the USA and Europe. However, Fig. 7a–c show that in the more polluted regions (USA–E and Europe), the magnitude of the bias is seasonal although still with a large positive bias. UKESM1 is able to capture the seasonal cycle in surface SO₂ concentrations over Europe, but the absolute model bias is larger in DJF compared with JJA (see Table 5). In USA–E UKESM1 does not capture the seasonal cycle in surface SO₂ concentrations due to the relatively large model bias in JJA, where the modelled SO₂ concentrations are over 5 times the observed values. In USA–W the modelled and observed surface SO₂ concentrations are slightly higher in DJF compared with JJA, but the model bias is so large in this region that it is difficult to determine whether there is any seasonality in this bias.

UKESM1 clearly captures the seasonal cycle in surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentration over USA–E, simulating the highest values in summer and the lowest values in winter. The model under-predicts surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentration by a factor of 0.7 to 0.8 reasonably consistently throughout the seasonal cycle. In the cleaner USA–W region, UKESM1 is able to capture the seasonal cycle in surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations, with the exception of DJF, where the model over-predicts surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations by a factor of 2.5. In Europe the observed seasonal cycle in surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentration has only a small amplitude, with mean values of 2.81 and 2.85 µg m⁻³ in DJF and JJA respectively.

3.5 Evaluation of total column SO₂ in UKESM1 against satellite observations

Figure 8 shows total column SO₂ (TCSO₂) from UKESM1 and OMI and the difference between them for DJF and JJA. Note that the quality control for solar zenith angle results in no data availability above 65^∘ N or below 65^∘ S in the winter months, and due to OMI's weaker sensitivity to retrieving SO₂ in remote regions, we focus on comparing TCSO₂ over source and outflow regions (Li et al., 2020). Figure 8a–d show that UKESM1 and OMI broadly agree on the location of the main Northern Hemisphere (NH) source regions, including China, India, Europe, and the USA. The model and satellite data both show seasonal cycles in TCSO₂ over the large source regions, with higher values being modelled and observed during the winter months. However, Fig. 8e and f show that the UKESM1 TCSO₂ values were generally larger than the OMI TCSO₂ values in these source regions by 0.6–1.0 DU. Over the background regions UKESM1 over-predicts TCSO₂ values by 0.2–0.5 DU. UKESM1 also has larger volcanic sources and associated outflow, which can be seen over central America, Sicily, Hawaii, and Papua New Guinea, for example. This is likely due to the climatology that UKESM1 uses for continuously degassing volcanoes. In agreement with the ground-based observations, the satellite data show an east–west divide in the USA, with greater TCSO₂ over USA–E compared with USA–W.

Figure 8 Total column SO₂ for UKESM1 (a, b), OMI (c, d), and UKESM1 – OMI (e, f). DJF is shown in the left column and JJA is shown in the right column. Median total column SO₂ in Dobson units is calculated for the period 2005–2014.

Figure 9 Median total column SO₂ calculated for 2005–2014 for USA (a), Europe (b), and South to North-east Asia (c). Total column SO₂ is in Dobson units.

Table 6 Statistics comparing model and OMI total column SO₂ over three regions: USA (60–130^∘ W, 25–50^∘ N), Europe (15^∘ W–40^∘ E, 35–65^∘ N), and South to North-east Asia (75–125^∘ W, 20–45^∘ N) for the period 2005–2014. The metrics are mean bias (MB, DU), root mean square error (RMSE, DU), percentage mean bias (MB %), and correlation (R).

Note that the median value is reported for each metric.

Figure 9 shows modelled and observed TCSO₂ over the period from 2005 to 2014 for three source regions, the USA (60–30^∘ W, 25–50^∘ N), Europe (15–40^∘ E, 35–65^∘ N), and South to North-east (SNE) Asia (75–125^∘ E, 20–45^∘ N). Overall, we find that the observed TCSO₂ is reasonably stable over this period in all three regions and that there are clear seasonal cycles showing peak TCSO₂ during the NH winter in the USA and Europe and slightly earlier ones in SNE Asia. However, Fig. 9a and b show that modelled TCSO₂ decreases over the USA and Europe from 2005 to 2010, with UKESM1 over-predicting TCSO₂ by up to 0.1 DU at the start of the time series. After 2011, UKESM1 is in much better agreement with the observed TCSO₂ over both regions. Figure 9c shows that UKESM1 consistently over-predicts TCSO₂ over SNE Asia by 1.5–2.0 DU during the period from 2005 to 2014. UKESM1 does simulate a seasonal cycle in TCSO₂ in all three regions. In Europe, UKESM1 is able to predict the peak wintertime TCSO₂ values, although between 2005 and 2010 the model has a positive bias of up to 0.1 DU. However, in the USA and to a lesser extent in SNE Asia (Fig. 9a and c respectively), UKESM1 mistimes the peak TCSO₂ values. In the USA the model simulates the highest values in the summer rather than in winter, and in SNE Asia the modelled peak TCSO₂ values appear shifted by several months earlier relative to the observations.

The observations of TCSO₂ and surface SO₂ concentrations over Europe both show the impact of emission control policies on keeping atmospheric SO₂ levels low (Figs. 2c and 9b). However, for the USA, we investigate the surface SO₂ concentrations separately over the USA–E and USA–W regions, whereas the TCSO₂ is averaged over the continental USA as a whole. As a result, TCSO₂ does not appear to decrease over the period 2005–2014 in the same way surface SO₂ concentrations are reduced in USA–E (Figs. 2b and 9a). Using a mean TCSO₂ for the continental USA as a whole also results in lower values relative to Europe. This is in contrast to the ground-based observations, where surface SO₂ concentrations over Europe are intermediate between USA–E and USA–W for the period 2005–2014. The relatively high TCSO₂ over Europe is also due to the inclusion of a number of large eastern European sources which are not well represented in the ground-based observations. We find that the TCSO₂ and surface SO₂ concentration observations both agree on the seasonal cycle, showing higher values in winter compared with the summer.

UKESM1 over-predicts surface SO₂ concentration to a greater extent than TCSO₂ for the USA and Europe. The model over-predicts surface SO₂ concentration by factors of 2.2–11.6 for USA–E and USA–W and 2.4 for Europe compared with values of 1.2 (USA) and 1.1 (Europe) for modelled TCSO₂ (see Tables 3 and 6). However, the R values for surface SO₂ concentration (>0.8) are much better than those for TCSO₂ (0.26–0.53), particularly over the USA. The low R value for the USA reflects the poor seasonal agreement in TCSO₂ in this region. The comparison against both observational data sets shows that the modelled atmospheric SO₂ is too high, both at the surface and through the column. In addition, UKESM1 simulates larger trends in TCSO₂ and surface SO₂, particularly prior to 2010, than are seen in the observations. Both observational data sets also show that UKESM1's over-prediction of atmospheric SO₂ in the large source regions is generally greater in the winter months compared with the summer months (see Fig. 8e and f). Exceptions occur in USA–E, where UKESM1 fails to capture the seasonal cycle in atmospheric SO₂, over-predicting surface SO₂ concentration and TCSO₂ to a greater extent in JJA compared with DJF. Notably, in the southern USA–E region and the Iberian Peninsula, UKESM1 actually under-predicts TCSO₂ by up to 0.1 DU in DJF, which does not occur in the comparison with surface SO₂ concentrations, even if the model and observations are compared at individual CASTNet and EMEP sites.

4 Impact of changes to the SO₂ dry deposition parameterization in UKESM1

4.1 Global-scale impacts

Figure 10 shows SO₂ dry deposition velocity simulated by UKESM1 and how this is affected by changes to the SO₂ dry deposition parameterization. In UKESM1 the mean annual deposition velocities range from approximately $\mathrm{0.5}\times {\mathrm{10}}^{-\mathrm{3}}$ to $\mathrm{4.0}\times {\mathrm{10}}^{-\mathrm{3}}$ m s⁻¹ (see Fig. 10a). Figure 10b, e, and h show that SO₂ dry deposition velocity increases over almost all land and ocean regions in UKESM1-SO2 by approximately 0.02–0.08 m s⁻¹. This represents an increase of a factor of 2–4 relative to UKESM1. The largest increase in SO₂ dry deposition velocity to the ocean is over the Southern Ocean in winter, where it increases by more than a factor of 4 (see Fig. 10i). Over land surfaces the largest increases occur over South America, north-western America and Canada, and north-eastern Europe/western Russia (see Fig. 10f and i). The increased SO₂ dry deposition velocities in UKESM1-SO2 relative to UKESM1 indicate that the changes to the SO₂ dry deposition parameterization are behaving as expected. The reduction in R_c for SO₂ dry deposition to water increases the dry deposition velocity to oceans. Similarly, when land surfaces are allowed to remain wet for a longer period after rainfall events, R_cut and R_soil are reduced for a longer period of time in UKESM1-SO2 relative to UKESM1 and SO₂ dry deposition velocity to the canopy (leaf or soil) increases. In the summer months SO₂ dry deposition velocities are larger over land surfaces compared with in the winter months, while over oceans values are larger in the winter compared with the summer. Although increased rainfall in winter drives wetter surface conditions, the leaf canopy is reduced or absent, and at high latitudes surfaces are likely to have snow cover, which has a relatively high R_c (see Table A1) compared with vegetated surfaces. The higher SO₂ dry deposition velocities over the ocean during the winter months are likely due to higher wind speeds.

Figure 10 SO₂ dry deposition velocity for UKESM1 (a, d, g), absolute difference in SO₂ dry deposition velocity for UKESM1-SO2 – UKESM1 (b, e, h), and percentage difference in SO₂ dry deposition velocity for UKESM1-SO2 – UKESM1 (c, f, i). Annual means are shown in (a)–(c), the DJF mean is shown in (d)–(f), and the JJA mean is shown in (g)–(i). Each plot shows mean data for the period 1987–2014. Deposition velocity is calculated from the simulated surface SO₂ and dry deposition flux values and represents the “bulk” deposition velocity for each model grid cell.

Figure 11 SO₂ dry deposition (a–c), surface SO₂ concentration (d–f), and surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentration (g–i) for UKESM1 (a, d, g), absolute difference for UKESM1-SO2 – UKESM1 (b, e, h), and percentage difference for UKESM1-SO2 – UKESM1 (c, f, i). In the plots showing absolute difference, areas where SO₂ dry deposition is reduced in UKESM1-SO2 compared with UKESM1 are shown in blue (b), and areas where surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations are greater in UKESM1-SO2 compared with UKESM1 are shown in red (e, h). The blue shading in (b) and the red shading in (e) and (h) is not to scale as the absolute differences are very small. Each plot shows mean data for the period 1987–2014.

The increased SO₂ dry deposition velocities in UKESM1-SO2 drive an increase in the SO₂ dry deposition flux of nearly 45 % (from 29.49 to 42.56 T g yr⁻¹) and subsequently reduce the SO₂ lifetime by approximately 25 % (see Table C1). Overall the global SO₂ burden is reduced from 0.54 to 0.41 T g yr⁻¹. The full budget breakdown for SO₂ is given in Appendix C (Table C1). The spatial distribution of these changes is illustrated in Fig. 11, which shows that SO₂ dry deposition increases over most regions, with a corresponding decrease in surface SO₂ concentrations. The lower SO₂ burden then drives lower oxidation fluxes (Table C1), reducing the surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations (Fig. 11g–i).

The largest absolute increases in SO₂ dry deposition occur over the main source regions of USA–E, eastern China, and central and eastern Europe (Fig. 11b). However, Fig. 11c shows that the largest relative changes in dry deposition are over the ocean. Although the absolute changes over the ocean are small ($<\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ kg m⁻² yr⁻¹), the large global surface area of the ocean means that this increase is important for the global sulfur cycle. In UKESM1-SO2 R_surf-mod increases the length of time land surfaces remain wet after rainfall, thus increasing SO₂ dry deposition over most land surfaces (as seen in Fig. 11b) due to the high solubility of SO₂ in water. Note that R_surf-mod also makes SO₂ dry deposition a function of the soil moisture content, and this aspect of the change drives decreases in SO₂ dry deposition of 20 %–40 % over desert regions, such as the Sahara and high northern latitudes (see Fig. 11b and c). Conversely R_surf-mod does not impact ocean surfaces, and the increases in SO₂ dry deposition over oceans are driven by R_water-mod.

The largest absolute reductions in mean annual surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations are over the source regions, corresponding to the locations of the largest increases in SO₂ dry deposition. This is expected because dry deposition of SO₂ to the surface is directly proportional to the surface concentration (Eq. 1). Figure 11e shows that mean annual surface SO₂ concentration was reduced by up to 20 µg m⁻³ in the eastern USA, eastern China, and central and eastern Europe, which corresponds to a percentage decrease of 30 %–50 %. Note that this is similar to the percentage decrease in mean annual surface SO₂ concentration over remote and ocean regions, although the absolute fluxes are much larger over the source regions. With the exception of some areas in the Sahara and the Middle East, we do not see increases in mean annual surface SO₂ concentration where dry deposition fluxes decrease (albeit by very low amounts), such as the Arctic. We suggest that this is because these areas are remote and contain no SO₂ sources, and by reducing SO₂ in the source regions, we reduce overall atmospheric SO₂ loading, and therefore less is transported to remote areas. Figure 11h and i show that mean annual surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations are also reduced over the main source regions, although the reductions over the USA are relatively small compared with the other large source regions (0.5 µg m⁻³ compared with up to 3 µg m⁻³ in central and eastern Europe and China). Proportionally, the reduction in mean annual surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentration is smaller than that for mean annual SO₂ concentration, with decreases generally less than 5 % over most source regions.

4.2 Evaluation of UKESM1-SO2 against ground-based observations

In Fig. 12 and Table 7 we evaluate UKESM1-SO2 against the ground-based observations of mean annual surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations for the USA and European regions over the period from 1987 to 2014. UKESM1-SO2 is also evaluated against SO₂ dry deposition flux from CASTNet over the same period. Figure 12a and b show that SO₂ dry deposition flux is increased in UKESM1-SO2 relative to UKESM1, with this increase being more pronounced over USA–E compared with USA–W (see also Table 7). The increase in SO₂ dry deposition in UKESM1-SO2 does enhance the model's over-prediction of this parameter relative to the CASTNet data, with NMB increasing from 1.15 in UKESM1 to 3.0 in UKESM1-SO2. Mean annual SO₂ dry deposition fluxes are very low over USA–W due to the much lower concentrations of SO₂ in this region. UKESM1 does over-predict mean annual SO₂ dry deposition flux in this region too, but the absolute bias changes very little in UKESM1-SO2 (see Table 7). Figure 12c–e show that model bias in mean annual surface SO₂ concentration is reduced in UKESM1-SO2 compared with UKESM1 in all three regions over the period 1987–2014. The largest absolute reduction is in USA–E, where mean annual surface SO₂ concentration decreases from 20.05 to 14.10 µg m⁻³; however, the largest reduction in NMB is in USA–W, where it decreases from 12.54 to 9.52 (see Table 7). The model's over-prediction of mean annual surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentration is reduced over USA–W in UKESM1-SO2 compared with UKESM1, with NMB decreasing from 0.72 to 0.41. However, the model's under-prediction of mean annual surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentration over USA–E and Europe increases, with NMB = −0.31 in UKESM1 and NMB = −0.47 in UKESM1-SO2 (see Table 7). We find that the changes to the surface concentration and dry deposition flux occur almost uniformly over the seasonal cycle and so do not change the patterns in seasonal bias that are described in Sect. 3.4 for surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations.

Figure 12 Time series of observed and modelled mean annual SO₂ dry deposition flux (a, b), surface SO₂ concentration (c–e), and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentration (f–h) for USA–W (a, c, f; N sites = 16), USA–E (b, d, g; N sites = 33), and Europe (e, h; N SO₂ sites = 48 and N ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ sites = 42). No SO₂ dry deposition flux observations are available for Europe.

Table 7 Statistics for mean annual surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations and SO₂ dry deposition flux for UKESM1 and UKESM1-SO2. The mean seasonal values are averaged over the period 1987–2014. The units for SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations are µg m⁻³ and the units for SO₂ dry deposition flux are kg m⁻² yr⁻¹ for SO₂.

4.3 Evaluation of UKESM1-SO2 against TCSO₂ observations

The impacts of the modifications to the SO₂ dry deposition parameterization on TCSO₂ are shown in Fig. 13 and Table 6. Figure 13a and b show that TCSO₂ over the source regions is lower in UKESM1-SO2 relative to UKESM1 by 0.1–0.5 DU in DJF and JJA. This results in UKESM1-SO2 having a smaller positive bias in TCSO₂ of +0.3 to +0.5 DU compared with that of +0.6 to +1.0 DU for UKESM1 (Fig. 13c and d; Sect. 3.5 and Fig. 8). This represents a decrease in the global TCSO₂ model–OMI bias of 20 %–30 %. Over the outflow regions (e.g. off the US eastern seaboard), TCSO₂ has reduced by 30 %–50 % and, over the source regions, this varies by 30 %–50 % for South to North-east Asia, 20 %–30 % for Europe, and 10 %–30 % for the USA (see Table 6). However, Fig. 13c and d also show that the inter-model differences are smaller than the existing model–satellite difference; i.e. the hatched regions are sporadic with limited coverage.

5 Discussion

The evaluation of UKESM1 against ground-based observations of SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations from the USA and Europe, as well as SO₂ dry deposition fluxes from the USA, shows that the model is able to represent recent historical changes in these variables. UKESM1 is also able to capture the spatial patterns in surface SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations and SO₂ dry deposition, simulating larger values close to the sources and lower values away from the sources. However, UKESM1 generally over-predicts surface SO₂ concentrations and dry deposition fluxes while under-predicting surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations for the period 1987–2014. Further, we find that UKESM1 over-predicts the rate at which the surface SO₂ concentrations decrease over this period.

We also make use of the updated TCSO₂ product from OMI to evaluate UKESM1, finding that the model captures spatial patterns in TCSO₂ at the global scale. Importantly, this evaluation allows us to identify model bias in regions without long-term ground-based networks, showing that UKESM1 over-predicts TCSO₂ over all source and outflow regions. We find that although the ground-based and satellite observations are subject to different uncertainties, UKESM1's relative over-prediction of both surface SO₂ concentration and TCSO₂ is similar in the USA and Europe. This suggests that our finding of positive bias in modelled atmospheric SO₂ is robust. We have also demonstrated that a more realistic treatment of SO₂ dry deposition in UKESM1 reduces the model's atmospheric loading of SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$. However, we find that UKESM1's under-prediction of surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations and over-prediction of SO₂ dry deposition fluxes increases when the changes are included in the model, suggesting that there are further uncertainties in UKESM1's representation of the complex sulfur cycle processes. Additionally, the spatial and temporal differences in the model bias suggest that the drivers of model bias are regionally and seasonally dependent.

Broadly, a model's over-prediction of atmospheric SO₂ can be driven by too little removal of SO₂ (via deposition or oxidation) or overly high emissions. UKESM1 uses SO₂ emissions from CMIP6 (Eyring et al., 2016). In comparing these emissions with the HTAP-OMI (Liu et al., 2018) and EDGAR (Crippa et al., 2018) data sets, Pope and Chipperfield (2021) showed that the total SO₂ emissions in CMIP6 (115 Tg yr⁻¹) are moderately larger than the HTAP-OMI (100 Tg yr⁻¹) and EDGAR data sets (102 Tg yr⁻¹). In general the CMIP6 emissions are larger than the HTAP-OMI and EDGAR emissions over all major source regions by up to $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{10}}$ kg m⁻² yr⁻¹. Exceptions occur at a number of point locations, which are likely from OMI rather than HTAP (Liu et al., 2018). Additionally, Smith (2021) has reported that SO₂ emissions over the western USA are too high, possibly because the proxy data used to spatially distribute emissions do not take into account lower sulfur coal used in power plants and industry in this region.

Emission injection height is also an important constraint on near-surface SO₂ concentrations, as demonstrated by Yang et al. (2019). This study found that uncertainty in industrial emission height resulted in modelled near‐surface SO₂ concentrations varying between 70 % and 130 % over most land regions, higher than the overall uncertainty of 8 %–14 % attributed to SO₂ emission rates. The SO₂ injection height in UKESM1 was investigated by Mulcahy et al. (2020), who used injection heights prescribed as for HadGEM-GC3.1, where 50 % of energy and industry sector emissions are injected into the atmosphere at a height of 500 m. Mulcahy et al. (2020) showed that the introduction of a vertical profile for SO₂ emissions in UKESM1 had negligible impact on surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations at measurement sites in Europe and the USA, suggesting an important role for the aerosol chemistry in these regions. Emitting the SO₂ at higher altitudes will act to reduce surface SO₂ concentrations and therefore a model's bias against surface observations of both SO₂ concentration and SO₂ dry deposition flux. However, Pope and Chipperfield (2021) showed that UKESM1's bias in TCSO₂ increased when emission injection heights were increased. This also suggests that the CMIP6 emissions are too high and that using a vertical profile for the emissions to some extent shifts the model's bias in SO₂ to higher altitudes. However, we note that using varying emission heights for SO₂ did not affect column densities in the GEOS-5/GOCART model (Buchard et al., 2014). We suggest undertaking model experiments with different emission inventories and injection height profiles to cast light on the role of SO₂ emissions in model bias in the sulfur cycle.

The two main removal pathways for SO₂ are oxidation to sulfate and dry deposition to the Earth's surface. In this study we have evaluated a more realistic treatment of SO₂ dry deposition in UKESM1 that accounts for the high solubility of SO₂ in water. We find that this reduces the dry deposition lifetime and consequently reduces the overall SO₂ burden and lifetime. This reduces positive model bias in surface SO₂ concentrations in the USA and Europe and in TCSO₂ across most of the globe. The changes have the largest impact over source regions because dry deposition flux is directly proportional to the atmospheric concentration of SO₂. There is also a reduction in the SO₂ oxidation lifetime which likely drives the reduced atmospheric loading of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$. However, the true impact of the changes to the dry deposition parameterization are confounded by model uncertainty in other aspects of the complex sulfur cycle as well as the inherent difficulties associated with evaluating a global model against point observations.

Figure 13 Difference in TCSO₂ between UKESM1 and UKESM-SO2 for DJF (a) and JJA (b). Difference in TCSO₂ between UKESM1-SO2 and OMI for DJF (c) and JJA (d). Model and satellite data are averaged over the period 2005–2014. The hatched regions show where the inter-model differences are smaller than the existing model–satellite difference.

Dry deposition is a highly parameterized process and often poorly represented, particularly in global models. Similarly to UKESM1, Vet et al. (2014) showed that SO₂ dry deposition fluxes were over-predicted by the 23-member model ensemble used for the TF-HTAP exercise relative to inferential data sets from measurement networks (including CASTNet). A key uncertainty highlighted in this study was that associated with the inferred dry deposition fluxes; from CASTNet these could be up to 30 % lower than direct observations of SO₂ dry deposition flux (Baumgardner et al., 2002). However, the fluxes simulated by UKESM1 are a factor of 2 to 10 higher than the inferred data from USA–E and USA–W, indicating that the modelled deposition fluxes are almost certainly too large. Additional sources of uncertainty in model simulations of SO₂ dry deposition may include land surface cover, changes in the atmospheric SO₂ : NH₃ ratio, and the ratio between wet and dry deposition of SO₂. Mulcahy et al. (2020) showed that wet deposition, and to a lesser extent dry deposition, of SO₂ was considerably lower in UKESM1 compared with HadGEM-GC3.1. Paulot et al. (2017) also report that poor representation of wet deposition likely contributed to bias in modelled surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations. Overly low wet deposition would also contribute to the model's over-prediction of surface SO₂ concentration. Dry deposition is sensitive to land surface type, which may not be well captured in global models. In this study the UKESM1 configuration uses 13 land cover classes, including 11 plant functional types (Archibald et al., 2020). This is reasonable for a global model, but inevitably detail is lost. Vet et al. (2014) also suggest that SO₂ dry deposition may depend on the atmospheric NH₃ loading. Long-term measurements at a UK site showed that SO₂ dry deposition velocity has increased with time, which was attributed to changing ratios of SO₂ : NH₃ as SO₂ concentrations have decreased more quickly than NH₃ concentrations (Vet et al., 2014; Fowler et al., 2009). Currently nitrate chemistry is not represented in UKESM1, although it is planned for future model versions, and NH₃ has not been evaluated in the model, so it is unknown how these factors may contribute to the model's bias in SO₂ dry deposition flux and SO₂ concentrations.

The role of uncertainty in sulfur cycle chemistry becomes apparent when we consider UKESM1's bias in surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations in combination with the biases in SO₂ concentrations and dry deposition. In Europe and USA–E, UKESM1 under-predicts surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations, despite the large positive biases in SO₂ concentrations through much of the period from 1987 to 2014. Note that there are exceptions close to certain point sources, particularly in Europe. However, in the cleaner USA–W region, surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations are consistently over-predicted. We suggest that in USA–W overly high SO₂ emissions (Smith, 2021), possibly in combination with overly low emission heights and the associated biases in dry deposition, drive the model's over-prediction of surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations in this region. However, in the polluted regions of Europe and USA–E, the under-prediction of surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations, despite large over-predictions of surface SO₂ concentrations, suggests that sulfur cycle chemistry is not correctly represented, and Mulcahy et al. (2020) showed that there are global and regional differences in oxidant concentrations and in the SO₂ lifetime between UKESM1 and the HadGEM-GC3.1 model, with the latter better able to capture surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations. In reducing the SO₂ burden, we further reduce its overall lifetime and oxidation lifetime relative to HadGEM-GC3.1. This highlights the requirement for a more detailed investigation of SO₂ oxidation in UKESM1, particularly in polluted regions.

Uncertainty in UKESM1's sulfur chemistry also appears to be seasonally dependent. In USA–E and Europe UKESM1 over-predicts SO₂ (surface concentration and TCSO₂) and under-predicts surface ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ at all times of the year, but there is seasonal variation in these biases. In Europe, the model bias is largest in DJF, which drives a stronger seasonal cycle relative to the observations. This is in agreement with Turnock et al. (2015), who investigated ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ in an earlier HadGEM3-UKCA configuration. The seasonality in UKESM1's bias over Europe may be due to the model under-predicting in-cloud SO₂ oxidation via O₃ (Table 1). In winter and at higher latitudes this reaction is likely to be the dominant removal pathway due to lower availability of H₂O₂ (Turnock et al., 2019). However, the picture is different in USA–E, where UKESM1 over-predicts SO₂ to a larger extent in JJA compared with DJF, and there is little seasonality in the model's under-prediction ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$. In USA–E it is likely that uncertainty in UKESM1's representation of SO₂ oxidation via both O₃ and H₂O₂ (Table 1) contributes to bias in SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$. The lower average latitude of the USA–E sites compared with the European sites means that the O₃ oxidation pathway is more important in this region, and it has been demonstrated that ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations are sensitive to both pathways in winter (Paulot et al., 2017). From this study it is not clear what is driving the relatively large positive bias in SO₂ in JJA over USA–E, and we stress the need for closer examination of the sulfur chemistry in UKESM1.

Recent studies have also shown that cloud water pH may be an important factor in the aqueous-phase oxidation of SO₂ to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ (Turnock et al., 2019). While a temporally and spatially uniform cloud pH of 5 is currently used in UKESM1, observations of this quantity show that it varies in space and time. Observations at an American site showed that mean cloud pH increased from 4 to 4.8 between 1994 and 2014 (Schwab et al., 2016), cloud pH measured at Mt Tai in North China from 2007 to 2008 was between 3.56 and 7.64, and measurement campaigns between 1985 and 2008 at various European, North American, and East Asian locations reported values between 3.34 and 5.29 (Guo et al., 2012). Turnock et al. (2019) showed that varying this value in the HadGEM3-UKCA model can have a large impact on SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations. Over source regions, including Europe and North America, increasing the cloud water pH by 1.0 reduced the annual mean global SO₂ column burden by approximately 50 %, as more SO₂ was oxidized in cloud droplets, and consequently there were small increases in the annual mean sulfate column burden over these regions. Conversely, outside of polluted regions, increasing the cloud water pH reduced the sulfate column burden by 10 % to 40 % globally. These results indicate that having a more realistic treatment of cloud water pH could reduce UKESM1's biases in TCSO₂ and potentially in SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations remote from source regions. However, it is unlikely to be a dominant removal pathway at the surface, and any impact on surface SO₂ concentrations, especially close to sources, would likely be minimal.

In an Earth system model such as UKESM1, there are inevitably some compromises in the complexity of the chemistry and aerosol schemes, as these are computationally expensive. While the sulfur chemistry represented in the UKCA-StratTrop model used in UKESM1 accounts for important SO₂ and DMS oxidation reactions as well as simulating oxidants (rather than the offline oxidant scheme used in HadGEM-GC3.1), it cannot be complete. In the recently developed CRI-Strat scheme the sulfur chemistry reactions are as for UKCA-StratTrop, but there is a more comprehensive treatment of non-methane volatile compounds (NMVOCs) (Archer-Nicholls et al., 2020), resulting in higher surface ozone concentrations, particularly over polluted areas in summer, compared with UKCA-StratTrop. As demonstrated by Mulcahy et al. (2020), the increased oxidants in UKESM1 relative to HadGEM-GC3.1 likely contribute to reducing the SO₂ lifetime from 4.29 to 3.86 d. CRI-Strat is compatible with UKESM1, and the higher oxidant loading may reduce SO₂ oxidation lifetime further with a concurrent increase in ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$. Model bias in remote ocean regions may also result from the necessarily simplified DMS oxidation chemistry in UKESM1. A more detailed representation of DMS chemistry over the Southern Ocean was investigated by Revell et al. (2019), who found that surface SO₂ concentrations increased over the Southern Ocean, possibly due to including reactions between DMS and halogen species, while ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations decreased, likely as a result of there being more DMS oxidation reactions.

Model resolution is also likely to be an important source of model bias in this study. In evaluating UKESM1 against the CASTNet and EMEP data sets we are comparing a simulated value generated from a model grid box at a scale of ≈ 200–300 km with a point observation. This may be a particular problem for the surface SO₂ concentrations and SO₂ dry deposition fluxes because in reality a large fraction (20 %–40 %) of SO₂ emitted from point sources is lost in the first 100 km, which is sub-grid scale relative to the model grid boxes (Smith and Jeffrey, 1975; Wys et al., 1978). In UKESM1 the sub-grid scale loss is not captured because SO₂ is evenly emitted across the grid box and deposition is subsequently calculated. Potentially this drives overall large model biases compared with ground-based observations that are not necessarily close to the point sources. In addition, the model resolution can not capture complex orography, meaning that transport of SO₂ may not be well represented. This could be problematic in mountainous areas of USA–W (VanCuren and Gustin, 2015). High-resolution model studies would be beneficial to address the both importance of orography and to investigate the SO₂ losses close to sources. We also suggest that evaluating UKESM1 against chemical re-analysis fields of SO₂ could reduce some of the bias that occurs with using point observations, but there are uncertainties associated with this approach too (Ukhov et al., 2020).

6 Conclusions

We evaluate UKESM1 against ground-based and satellite observations of selected sulfur species over the recent historical period. We find that UKESM1 is able to capture the temporal and spatial patterns in surface SO₂, surface sulfate concentration, SO₂ dry deposition, and TCSO₂. However, compared with observations, we find that the model is biased, depending on the variable, region, and species. We address one possible source of bias by introducing a more realistic treatment of SO₂ dry deposition, a key loss process for this species. This change reduces model bias in surface SO₂ concentrations and TCSO₂. However, it is apparent that other biases exist within the complex sulfur cycle, and we highlight some key areas for further investigation and development.

Our evaluation suggests that uncertainty in UKESM1's sulfur chemistry is also an important driver of the biases seen in this study, particularly over polluted regions. Two priorities for further investigation into the oxidation of SO₂ in UKESM1 are (i) an evaluation of the CRI-Strat scheme and (ii) a more realistic treatment of cloud water pH. The model's necessarily limited DMS chemistry may also contribute to bias in atmospheric sulfur loading over remote ocean areas. The impact of the nitrate scheme currently in development for UKESM1 will also be investigated in relation to the sulfur cycle. Another aspect of UKESM1's sulfur cycle that would benefit from more detailed analysis is the ratio between wet and dry deposition and how this compares with observations. We also suggest that high-resolution studies to investigate SO₂ deposition close to sources would be beneficial for a better understanding of these processes. Finally, we suggest that the SO₂ emissions may be too high through a possible combination of overly high emissions in the CMIP6 inventory and injection of the emissions into the surface layer only.

The sulfur cycle is a key area of analysis and development for UKESM1 given its importance as a driver of historical aerosol forcing. UKESM1 is relatively unique amongst models in CMIP6 in that it has a fully interactive atmospheric chemistry scheme coupled to a two-moment (mass and number) aerosol scheme. Given the complexity of the model's chemistry-aerosol treatment within the ES framework, the model's performance here is encouraging and provides confidence in UKESM1, particularly in regard to capturing the historical trend. However, there is always space for improvement, and to the more realistic treatment of SO₂ dry deposition will therefore be incorporated into the planned release of UKESM1.1. This latest model version will be documented in a forthcoming publication which will also address the impact of the SO₂ dry deposition changes on aerosol loading and climate.

Appendix A

Here we describe in detail the changes to UKESM1's parameterization of dry deposition of SO₂ to the surface. These modifications account for the impact of surface wetness due to rainfall or humidity on SO₂ dry deposition to vegetated or soil surfaces. We first derive a parameter, z_wet, which designates whether a model grid box at time step N is wet or dry. We assume that on entering time step N the grid box in question has a dry surface (z_wet = 0). If, during time step N, precipitation is greater than a threshold value (set to 0.5 mm d⁻¹), the grid box then becomes classed as wet (z_wet = 1). Once precipitation stops, the grid box is assumed to dry out over a specified period. Assuming no precipitation falls during this drying period, at the end of the period the grid box will be classed as dry (z_wet = 0). If, during the drying period, a new precipitation event occurs at the grid box of intensity greater than 0.5 mm d⁻¹, z_wet is reset to 1. If the precipitation event is less than 0.5 mm d⁻¹ but greater than 0.0 mm d⁻¹, z_wet is not reset to 1, and neither is it decreased in value (i.e. no drying is assumed to occur over that time step). A grid box is classed as wet whenever z_wet > 0. We tested a range of time periods from 3 h to 1 d and found only minor sensitivity to this parameter. For UKESM1-SO2 we used a drying period of 3 h. In future work we will investigate making this parameter a function of surface evaporation and downwelling solar radiation.

If a grid box is classed as wet, then R_soil and R_cut (for all vegetation types) are set equal to 1 s m⁻¹. Through Eq. (5), R_c will then tend towards a value of 1 s m⁻¹, equating to minimal resistance to SO₂ deposition. In these situations, the amount of SO₂ deposited will primarily be limited by the efficiency of gas transport to the receptor surface, i.e. by R_a and R_b in Eq. (2). If the grid box is classed as dry, R_cut is calculated following Eq. (9) in Erisman et al. (1994), whereby R_cut is a decreasing function of near-surface relative humidity. If a grid box surface temperature lies between −1 and −5 ^∘C, R_cut is reset to 200 s m⁻¹ and below −5 ^∘C to 500 s m⁻¹, irrespective of the near-surface relative humidity. For dry grid boxes, R_soil uses a value of 213.5 s m⁻¹ for all surface temperatures (see Table A2). For dry surfaces, R_cut also approaches a value of 1 s m⁻¹ as near-surface relative humidity approaches a value of 1. As for a wet surface, under these conditions there is minimal resistance to SO₂ deposition, and R_c will tend towards 1 s m⁻¹. Thus SO₂ dry deposition will primarily be limited by R_a and R_b.

Table A1 Summary of the representation of R_cut and R_soil in UKESM1 and UKESM1-SO2. Units for resistance values are s m⁻¹.

^a See values in Table A2. ^b Following Erisman et al. (1994). ^c Irrespective of near-surface relative humidity. n/a – not applicable.

Table A2 Standard surface resistance (R_surf) values for SO₂ for land use types in UKESM1. These values were calculated based on the data published in Zhang et al. (2003) for SO₂.

Appendix B

Figure B1 Global emissions of SO₂ used in UKESM1. Mean seasonal emissions are shown for the 1990–1995 time slice for DJF (a) and JJA (b) and the 2009–2014 time slice for DJF (c) and JJA (d).

Appendix C

Table C1 Global SO₂ budget for UKESM1 and UKESM1-SO2. Units for production and loss fluxes are in Tg [S] yr⁻¹, burdens are in Tg, and lifetimes are in days. The values are calculated from a 2-year AMIP simulation covering the period 1981–1983 inclusive.